From fa3a647a74cf73d22db39355857f3e0e0319256f Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 6 Nov 2009 15:38:50 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3320 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/Makefile | 3 +- src/USER-ATC/README | 3 + src/USER-CG-CMM/README | 4 +- src/USER-IMD/Install.csh | 20 + src/USER-IMD/README | 17 + src/USER-IMD/fix_imd.cpp | 1325 +++++++++++++++++++++++++++++++++ src/USER-IMD/fix_imd.h | 99 +++ src/USER-IMD/style_user_imd.h | 20 + src/style_user_packages.h | 1 + 9 files changed, 1489 insertions(+), 3 deletions(-) create mode 100644 src/USER-IMD/Install.csh create mode 100644 src/USER-IMD/README create mode 100644 src/USER-IMD/fix_imd.cpp create mode 100644 src/USER-IMD/fix_imd.h create mode 100644 src/USER-IMD/style_user_imd.h diff --git a/src/Makefile b/src/Makefile index 0fde23a842..6f37610abb 100755 --- a/src/Makefile +++ b/src/Makefile @@ -16,7 +16,8 @@ OBJ = $(SRC:.cpp=.o) PACKAGE = asphere class2 colloid dipole dpd gpu granular \ kspace manybody meam molecule opt peri poems prd reax xtc -PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-ewaldn user-smd +PACKUSER = user-ackland user-atc user-cd-eam user-cg-cmm user-ewaldn \ + user-imd user-smd PACKALL = $(PACKAGE) $(PACKUSER) diff --git a/src/USER-ATC/README b/src/USER-ATC/README index 373c25cda6..3e84669048 100644 --- a/src/USER-ATC/README +++ b/src/USER-ATC/README @@ -10,3 +10,6 @@ LAMMPS input script. This fix can be employed to either do concurrent coupling of MD with FE-based physics surrogates or on-the-fly post-processing of atomic information to continuum fields. See the doc page for the fix atc command to get started. + +There are example scripts for using this package with LAMMPS in +examples/USER/atc. diff --git a/src/USER-CG-CMM/README b/src/USER-CG-CMM/README index 2bc51cef7e..a7ff592f05 100644 --- a/src/USER-CG-CMM/README +++ b/src/USER-CG-CMM/README @@ -18,8 +18,8 @@ and angle_style cg/cmm. See the documentation files for these commands for details. -There is also an cg-cmm example directory under examples/USER with -sample inputs and outputs. +There are example scripts for using this package with LAMMPS in +examples/USER/cg-cmm. These styles allow coarse grained MD simulations with the parametrization of Shinoda, DeVane, Klein, Mol Sim, 33, 27 (2007) diff --git a/src/USER-IMD/Install.csh b/src/USER-IMD/Install.csh new file mode 100644 index 0000000000..5e0474673e --- /dev/null +++ b/src/USER-IMD/Install.csh @@ -0,0 +1,20 @@ +# Install/unInstall package classes in LAMMPS + +if ($1 == 1) then + + cp -p style_user_imd.h .. + + cp -p fix_imd.cpp .. + + cp -p fix_imd.h .. + +else if ($1 == 0) then + + rm ../style_user_imd.h + touch ../style_user_imd.h + + rm ../fix_imd.cpp + + rm ../fix_imd.h + +endif diff --git a/src/USER-IMD/README b/src/USER-IMD/README new file mode 100644 index 0000000000..65316c0bdb --- /dev/null +++ b/src/USER-IMD/README @@ -0,0 +1,17 @@ +The files in this directory are a user-contributed package for LAMMPS. + +The person who created these files is Axel Kohlmeyer +(axel.kohlmeyer@temple.edu). Contact him directly if you +have questions or for example scripts that use it. + +This package implements a "fix imd" command which can be used in a +LAMMPS input script. IMD stands for interactive molecular dynamics, +and allows realtime visualization and manipulation of MD simulations +through the IMD protocol, initially implemented in VMD and NAMD. + +There are example scripts for using this package with LAMMPS in +examples/USER/imd. + +This software includes code developed by the Theoretical and Computational +Biophysics Group in the Beckman Institute for Advanced Science and +Technology at the University of Illinois at Urbana-Champaign. diff --git a/src/USER-IMD/fix_imd.cpp b/src/USER-IMD/fix_imd.cpp new file mode 100644 index 0000000000..01c84798d2 --- /dev/null +++ b/src/USER-IMD/fix_imd.cpp @@ -0,0 +1,1325 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + The FixIMD class contains code from VMD and NAMD which is copyrighted + by the Board of Trustees of the University of Illinois and is free to + use with LAMMPS according to point 2 of the UIUC license (10% clause): + +" Licensee may, at its own expense, create and freely distribute +complimentary works that interoperate with the Software, directing others to +the TCBG server to license and obtain the Software itself. Licensee may, at +its own expense, modify the Software to make derivative works. Except as +explicitly provided below, this License shall apply to any derivative work +as it does to the original Software distributed by Illinois. Any derivative +work should be clearly marked and renamed to notify users that it is a +modified version and not the original Software distributed by Illinois. +Licensee agrees to reproduce the copyright notice and other proprietary +markings on any derivative work and to include in the documentation of such +work the acknowledgement: + + "This software includes code developed by the Theoretical and Computational + Biophysics Group in the Beckman Institute for Advanced Science and + Technology at the University of Illinois at Urbana-Champaign." + +Licensee may redistribute without restriction works with up to 1/2 of their +non-comment source code derived from at most 1/10 of the non-comment source +code developed by Illinois and contained in the Software, provided that the +above directions for notice and acknowledgement are observed. Any other +distribution of the Software or any derivative work requires a separate +license with Illinois. Licensee may contact Illinois (vmd@ks.uiuc.edu) to +negotiate an appropriate license for such distribution." +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Axel Kohlmeyer (TempleU) + IMD API, hash, and socket code written by: John E. Stone, + Justin Gullingsrud, and James Phillips, (TCBG, Beckman Institute, UIUC) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include + +#if defined(_MSC_VER) +#include +#else +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +#include + +#include "fix_imd.h" +#include "atom.h" +#include "comm.h" +#include "update.h" +#include "respa.h" +#include "domain.h" +#include "error.h" +#include "group.h" +#include "memory.h" + +/********** API definitions of the VMD/NAMD code ************************ + * This code was taken and adapted from VMD-1.8.7/NAMD-2.7 in Sep 2009. * + * If there are any bugs or problems, please contact akohlmey@gmail.com * + ************************************************************************/ + +/*************************************************************************** + *cr + *cr (C) Copyright 1995-2009 The Board of Trustees of the + *cr University of Illinois + *cr All Rights Reserved + *cr + ***************************************************************************/ + +/* part 1: integer hash table */ + +/** hash table top level data structure */ +typedef struct inthash_t { + struct inthash_node_t **bucket; /* array of hash nodes */ + int size; /* size of the array */ + int entries; /* number of entries in table */ + int downshift; /* shift cound, used in hash function */ + int mask; /* used to select bits for hashing */ +} inthash_t; + +/** hash table node data structure */ +typedef struct inthash_node_t { + int data; /* data in hash node */ + int key; /* key for hash lookup */ + struct inthash_node_t *next; /* next node in hash chain */ +} inthash_node_t; + +#define HASH_FAIL -1 +#define HASH_LIMIT 0.5 + +/* initialize new hash table */ +static void inthash_init(inthash_t *tptr, int buckets); +/* lookup entry in hash table */ +static int inthash_lookup(const inthash_t *tptr, int key); +/* generate list of keys for reverse lookups. */ +static int *inthash_keys(inthash_t *tptr); +/* insert an entry into hash table. */ +static int inthash_insert(inthash_t *tptr, int key, int data); +/* delete the hash table */ +static void inthash_destroy(inthash_t *tptr); + +/* part 2: Interactive MD (IMD) API */ + +#include + +#if ( INT_MAX == 2147483647 ) +typedef int int32; +#else +typedef short int32; +#endif + +typedef struct { + int32 type; + int32 length; +} IMDheader; + +#define IMDHEADERSIZE 8 +#define IMDVERSION 2 + +typedef enum IMDType_t { + IMD_DISCONNECT, /**< close IMD connection, leaving sim running */ + IMD_ENERGIES, /**< energy data block */ + IMD_FCOORDS, /**< atom coordinates */ + IMD_GO, /**< start the simulation */ + IMD_HANDSHAKE, /**< endianism and version check message */ + IMD_KILL, /**< kill the simulation job, shutdown IMD */ + IMD_MDCOMM, /**< MDComm style force data */ + IMD_PAUSE, /**< pause the running simulation */ + IMD_TRATE, /**< set IMD update transmission rate */ + IMD_IOERROR /**< indicate an I/O error */ +} IMDType; /**< IMD command message type enumerations */ + +typedef struct { + int32 tstep; /**< integer timestep index */ + float T; /**< Temperature in degrees Kelvin */ + float Etot; /**< Total energy, in Kcal/mol */ + float Epot; /**< Potential energy, in Kcal/mol */ + float Evdw; /**< Van der Waals energy, in Kcal/mol */ + float Eelec; /**< Electrostatic energy, in Kcal/mol */ + float Ebond; /**< Bond energy, Kcal/mol */ + float Eangle; /**< Angle energy, Kcal/mol */ + float Edihe; /**< Dihedral energy, Kcal/mol */ + float Eimpr; /**< Improper energy, Kcal/mol */ +} IMDEnergies; /**< IMD simulation energy report structure */ + +/** Send control messages - these consist of a header with no subsequent data */ +static int imd_handshake(void *); /**< check endianness, version compat */ +/** Receive header and data */ +static IMDType imd_recv_header(void *, int32 *); +/** Receive MDComm-style forces, units are Kcal/mol/angstrom */ +static int imd_recv_mdcomm(void *, int32, int32 *, float *); +/** Receive energies */ +static int imd_recv_energies(void *, IMDEnergies *); +/** Receive atom coordinates. */ +static int imd_recv_fcoords(void *, int32, float *); +/** Prepare IMD data packet header */ +static void imd_fill_header(IMDheader *header, IMDType type, int32 length); +/** Write data to socket */ +static int32 imd_writen(void *s, const char *ptr, int32 n); + +/* part 3: abstracts platform-dependent routines/APIs for using sockets */ + +typedef struct { + struct sockaddr_in addr; /* address of socket provided by bind() */ + int addrlen; /* size of the addr struct */ + int sd; /* socket file descriptor */ +} imdsocket; + +static int imdsock_init(void); +static void *imdsock_create(void); +static int imdsock_bind(void *, int); +static int imdsock_listen(void *); +static void *imdsock_accept(void *); /* return new socket */ +static int imdsock_write(void *, const void *, int); +static int imdsock_read(void *, void *, int); +static int imdsock_selread(void *, int); +static int imdsock_selwrite(void *, int); +static void imdsock_shutdown(void *); +static void imdsock_destroy(void *); + +/*************************************************************** + * End of API definitions of the VMD/NAMD code. * + * The implementation follows at the end of the file. * + ***************************************************************/ + +using namespace LAMMPS_NS; + +/* adapted sort for in-place sorting of map indices. */ +static void id_sort(int *idmap, int left, int right); + +/* struct for packed data communication of coordinates and forces. */ +struct commdata { + int tag; + float x,y,z; +}; + +/*************************************************************** + * create class and parse arguments in LAMMPS script. Syntax: + * fix ID group-ID imd [unwrap (on|off)] [fscale ] + ***************************************************************/ +FixIMD::FixIMD(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg < 4) + error->all("Illegal fix imd command"); + + imd_port = atoi(arg[3]); + if (imd_port < 1024) + error->all("Illegal fix imd parameter. port < 1024."); + + /* default values for optional flags */ + unwrap_flag = 0; + imd_fscale = 1.0; + imd_trate = 1; + + /* parse optional arguments */ + int argsdone = 4; + while (argsdone+1 < narg) { + if (0 == strcmp(arg[argsdone], "unwrap")) { + if (0 == strcmp(arg[argsdone+1], "on")) { + unwrap_flag = 1; + } else { + unwrap_flag = 0; + } + } else if (0 == strcmp(arg[argsdone], "fscale")) { + imd_fscale = atof(arg[argsdone+1]); + } else if (0 == strcmp(arg[argsdone], "trate")) { + imd_trate = atoi(arg[argsdone+1]); + } else { + error->all("Unknown fix imd parameter."); + } + ++argsdone; ++argsdone; + } + + /* sanity check on parameters */ + if (imd_trate < 1) + error->all("Illegal fix imd parameter. trate < 1."); + + if (igroup == group->find("all")) { + num_coords = static_cast (atom->natoms); + } else { + num_coords = static_cast (group->count(igroup)); + if (num_coords <= 0) error->all("Invalid number of group atoms for 'fix imd'"); + } + + MPI_Comm_rank(world,&me); + + /* initialize various imd state variables. */ + clientsock = NULL; + localsock = NULL; + nlevels_respa = 0; + imd_inactive = 0; + imd_terminate = 0; + imd_forces = 0; + force_buf = NULL; + maxbuf = 0; + comm_buf = NULL; + idmap = NULL; + rev_idmap = NULL; + + /* storage required to communicate a single coordinate or force. */ + size_one = sizeof(struct commdata); +} + +/********************************* + * Clean up on deleting the fix. * + *********************************/ +FixIMD::~FixIMD() +{ + inthash_t *hashtable = (inthash_t *)idmap; + memory->sfree(comm_buf); + memory->sfree(force_buf); + inthash_destroy(hashtable); + delete hashtable; + free(rev_idmap); + imdsock_shutdown(clientsock); + imdsock_destroy(localsock); + return; +} + +/* ---------------------------------------------------------------------- */ +int FixIMD::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + mask |= POST_FORCE_RESPA; + return mask; +} + +/* ---------------------------------------------------------------------- */ +void FixIMD::init() +{ + if (me == 0) { + /* set up incoming socket on MPI rank 0. */ + imdsock_init(); + localsock = imdsock_create(); + clientsock = NULL; + if (imdsock_bind(localsock,imd_port)) { + perror("bind to socket failed"); + imdsock_destroy(localsock); + imd_terminate = 1; + } else { + imdsock_listen(localsock); + } + } + MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); + if (imd_terminate) + error->all("LAMMPS Terminated on error in IMD."); + + if (strcmp(update->integrate_style,"respa") == 0) + nlevels_respa = ((Respa *) update->integrate)->nlevels; + + return; +} + +/* ---------------------------------------------------------------------- */ +/* wait for IMD client (e.g. VMD) to respond, initialize communication + * buffers and collect tag/id maps. */ +void FixIMD::setup(int) +{ + /* nme: number of atoms in group on this MPI task + * nmax: max number of atoms in group across all MPI tasks + * nlocal: all local atoms + */ + int i,j; + int nmax,nme,nlocal; + int *mask = atom->mask; + int *tag = atom->tag; + nlocal = atom->nlocal; + nme=0; + for (i=0; i < nlocal; ++i) + if (mask[i] & groupbit) ++nme; + + MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); + maxbuf = nmax*size_one; + comm_buf = memory->smalloc(maxbuf,"imd:comm_buf"); + + /* set up IMD communication. */ + if (me == 0) { + if (screen) + fprintf(screen,"Waiting for IMD connection on port %d.\n",imd_port); + + int retval=0; + do { + retval = imdsock_selread(localsock, 60); + } while (retval <= 0); + clientsock = imdsock_accept(localsock); + + if (!clientsock) { + if (screen) + fprintf(screen, "IMD socket accept error. Dropping connection.\n"); + imd_terminate = 1; + } else { + /* check endianness and IMD protocol version. */ + if (imd_handshake(clientsock)) { + if (screen) + fprintf(screen, "IMD handshake error. Dropping connection.\n"); + imdsock_destroy(clientsock); + imd_terminate = 1; + } else { + int32 length; + if (imdsock_selread(clientsock, 1) != 1 || + imd_recv_header(clientsock, &length) != IMD_GO) { + if (screen) + fprintf(screen, "Incompatible IMD client version? Dropping connection.\n"); + imdsock_destroy(clientsock); + imd_terminate = 1; + } + } + } + } + MPI_Bcast(&imd_terminate, 1, MPI_INT, 0, world); + if (imd_terminate) + error->all("LAMMPS terminated on error in setting up IMD connection."); + + /* initialize and build hashtable. */ + inthash_t *hashtable=new inthash_t; + inthash_init(hashtable, num_coords); + idmap = (void *)hashtable; + + MPI_Status status; + MPI_Request request; + int tmp, ndata; + struct commdata *buf = static_cast(comm_buf); + + if (me == 0) { + int *taglist = new int[num_coords]; + int numtag=0; /* counter to map atom tags to a 0-based consecutive index list */ + + for (i=0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + taglist[numtag] = tag[i]; + ++numtag; + } + } + + /* loop over procs to receive remote data */ + for (i=1; i < comm->nprocs; ++i) { + MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request); + MPI_Send(&tmp, 0, MPI_INT, i, 0, world); + MPI_Wait(&request, &status); + MPI_Get_count(&status, MPI_BYTE, &ndata); + ndata /= size_one; + + for (j=0; j < ndata; ++j) { + taglist[numtag] = buf[j].tag; + ++numtag; + } + } + + /* sort list of tags by value to have consistently the + * same list when running in parallel and build hash table. */ + id_sort(taglist, 0, num_coords-1); + for (i=0; i < num_coords; ++i) { + inthash_insert(hashtable, taglist[i], i); + } + delete[] taglist; + + /* generate reverse index-to-tag map for communicating + * IMD forces back to the proper atoms */ + rev_idmap=inthash_keys(hashtable); + } else { + nme=0; + for (i=0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + buf[nme].tag = tag[i]; + ++nme; + } + } + /* blocking receive to wait until it is our turn to send data. */ + MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, &status); + MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world); + } + + return; +} + +/* ---------------------------------------------------------------------- */ +/* Main IMD protocol handler: + * Send coodinates, energies, and add IMD forces to atoms. */ +void FixIMD::post_force(int vflag) +{ + if (imd_inactive) return; /* IMD client has detached. do nothing. */ + + int *tag = atom->tag; + double **x = atom->x; + int *image = atom->image; + int nlocal = atom->nlocal; + int *mask = atom->mask; + struct commdata *buf; + + /* Check if we need to communicate coordinates to the client. + * Tuning imd_trate allows to keep the overhead for IMD low + * at the expense of a more jumpy display. Rather than using + * end_of_step() we do everything here in one go. + * + * If we don't communicate, only check if we have forces + * stored away and apply them. */ + if (update->ntimestep % imd_trate) { + if (imd_forces > 0) { + double **f = atom->f; + buf = static_cast(force_buf); + + /* XXX. this is in principle O(N**2) = not good. + * however we assume for now that the number of atoms + * that we manipulate via IMD will be small compared + * to the total system size, so we don't hurt too much. */ + for (int j=0; j < imd_forces; ++j) { + for (int i=0; i < nlocal; ++i) { + if (mask[i] & groupbit) { + if (buf[j].tag == tag[i]) { + f[i][0] += buf[j].x; + f[i][1] += buf[j].y; + f[i][2] += buf[j].z; + } + } + } + } + } + return; + } + + /* check and potentially grow local communication buffers. */ + int i, k, nmax, nme=0; + for (i=0; i < nlocal; ++i) + if (mask[i] & groupbit) ++nme; + + MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world); + if (nmax*size_one > maxbuf) { + memory->sfree(comm_buf); + maxbuf = nmax*size_one; + comm_buf = memory->smalloc(maxbuf,"imd:comm_buf"); + } + + MPI_Status status; + MPI_Request request; + int tmp, ndata; + buf = static_cast(comm_buf); + + if (me == 0) { + /* collect data into new array. we bypass the IMD API to save + * us one extra copy of the data. */ + int msglen = 3*sizeof(float)*num_coords+IMDHEADERSIZE; + char *msgdata = new char[msglen]; + imd_fill_header((IMDheader *)msgdata, IMD_FCOORDS, num_coords); + /* array pointer, to the offset where we receive the coordinates. */ + float *recvcoord = (float *) (msgdata+IMDHEADERSIZE); + + /* add local data */ + if (unwrap_flag) { + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + + for (i=0; i> 10 & 1023) - 512) * yprd; + recvcoord[j+2] = x[i][2] + ((image[i] >> 20) - 512) * zprd; + } + } + } + } else { + for (i=0; inprocs; ++i) { + MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request); + MPI_Send(&tmp, 0, MPI_INT, i, 0, world); + MPI_Wait(&request, &status); + MPI_Get_count(&status, MPI_BYTE, &ndata); + ndata /= size_one; + + for (k=0; k 0) || imd_paused) { + /* if something requested to turn off IMD while paused get out */ + if (imd_inactive) break; + + int32 length; + int msg = imd_recv_header(clientsock, &length); + + switch(msg) { + + case IMD_GO: + if (screen) + fprintf(screen, "Ignoring unexpected IMD_GO message.\n"); + break; + + case IMD_IOERROR: + if (screen) + fprintf(screen, "IMD connection lost.\n"); + /* fallthrough */ + + case IMD_DISCONNECT: { + /* disconnect from client. wait for new connection. */ + imd_paused = 0; + imd_forces = 0; + memory->sfree(force_buf); + force_buf = NULL; + imdsock_destroy(clientsock); + clientsock = NULL; + if (screen) { + fprintf(screen, "IMD client detached. LAMMPS run continues.\n"); + fprintf(screen, "Waiting for new IMD connection on port %d.\n",imd_port); + } + int retval=0; + do { + retval = imdsock_selread(localsock, 60); + } while (retval <= 0); + clientsock = imdsock_accept(localsock); + if (!clientsock) { + if (screen) + fprintf(screen, "IMD socket accept error. Dropping connection.\n"); + imd_terminate = 1; + } else { + /* check endianness and IMD protocol version. */ + if (imd_handshake(clientsock)) { + if (screen) + fprintf(screen, "IMD handshake error. Dropping connection.\n"); + imdsock_destroy(clientsock); + imd_terminate = 1; + } else { + int32 length; + if (imdsock_selread(clientsock, 1) != 1 || + imd_recv_header(clientsock, &length) != IMD_GO) { + if (screen) + fprintf(screen, "Incompatible IMD client version? Dropping connection.\n"); + imdsock_destroy(clientsock); + imd_terminate = 1; + } + } + } + if (imd_terminate) imd_inactive = 1; + break; + } + + case IMD_KILL: + /* stop the simulation job and shutdown IMD */ + if (screen) + fprintf(screen, "IMD client requested termination of run.\n"); + imd_inactive = 1; + imd_terminate = 1; + imd_paused = 0; + imdsock_destroy(clientsock); + clientsock = NULL; + break; + + case IMD_PAUSE: + /* pause the running simulation. wait for second IMD_PAUSE to continue. */ + if (imd_paused) { + if (screen) + fprintf(screen, "Continuing run on IMD client request.\n"); + imd_paused = 0; + } else { + if (screen) + fprintf(screen, "Pausing run on IMD client request.\n"); + imd_paused = 1; + } + break; + + case IMD_TRATE: + /* change the IMD transmission data rate */ + if (length > 0) + imd_trate = length; + if (screen) + fprintf(screen, "IMD client requested change of transfer rate. Now it is %d.\n", imd_trate); + break; + + case IMD_ENERGIES: { + IMDEnergies dummy_energies; + imd_recv_energies(clientsock, &dummy_energies); + break; + } + + case IMD_FCOORDS: { + float *dummy_coords = new float[3*length]; + imd_recv_fcoords(clientsock, length, dummy_coords); + delete[] dummy_coords; + break; + } + + case IMD_MDCOMM: { + int32 *imd_tags = new int32[length]; + float *imd_fdat = new float[3*length]; + imd_recv_mdcomm(clientsock, length, imd_tags, imd_fdat); + + if (imd_forces < length) { /* grow holding space for forces, if needed. */ + if (force_buf != NULL) + memory->sfree(force_buf); + force_buf = memory->smalloc(length*size_one, "imd:force_buf"); + } + imd_forces = length; + buf = static_cast(force_buf); + + /* compare data to hash table */ + for (int i=0; i < length; ++i) { + buf[i].tag = rev_idmap[imd_tags[i]]; + buf[i].x = imd_fdat[3*i]; + buf[i].y = imd_fdat[3*i+1]; + buf[i].z = imd_fdat[3*i+2]; + } + delete[] imd_tags; + delete[] imd_fdat; + break; + } + + default: + if (screen) + fprintf(screen, "Unhandled incoming IMD message #%d. length=%d\n", msg, length); + break; + } + } + } else { + /* copy coordinate data into communication buffer */ + nme = 0; + if (unwrap_flag) { + double xprd = domain->xprd; + double yprd = domain->yprd; + double zprd = domain->zprd; + + for (i=0; i> 10 & 1023) - 512) * yprd; + buf[nme].z = x[i][2] + ((image[i] >> 20) - 512) * zprd; + ++nme; + } + } + } else { + for (i=0; iall("LAMMPS terminated on IMD request."); + + if (imd_forces > 0) { + /* check if we need to readjust the forces comm buffer on the receiving nodes. */ + if (me != 0) { + if (old_imd_forces < imd_forces) { /* grow holding space for forces, if needed. */ + if (force_buf != NULL) + memory->sfree(force_buf); + force_buf = memory->smalloc(imd_forces*size_one, "imd:force_buf"); + } + } + MPI_Bcast(force_buf, imd_forces*size_one, MPI_BYTE, 0, world); + } + return; +} + +/* ---------------------------------------------------------------------- */ +void FixIMD::post_force_respa(int vflag, int ilevel, int iloop) +{ + /* only process IMD on the outmost RESPA level. */ + if (ilevel == nlevels_respa-1) post_force(vflag); + return; +} + +/* ---------------------------------------------------------------------- */ +/* local memory usage. approximately. */ +double FixIMD::memory_usage(void) +{ + return static_cast(num_coords+maxbuf+imd_forces)*size_one; +} + +/* End of FixIMD class implementation. */ + +/************************************************************************ + * integer list sort code: + ************************************************************************/ + +/* sort for integer map. initial call id_sort(idmap, 0, natoms - 1); */ +void id_sort(int *idmap, int left, int right) +{ + int pivot, l_hold, r_hold; + + l_hold = left; + r_hold = right; + pivot = idmap[left]; + + while (left < right) { + while ((idmap[right] >= pivot) && (left < right)) + right--; + if (left != right) { + idmap[left] = idmap[right]; + left++; + } + while ((idmap[left] <= pivot) && (left < right)) + left++; + if (left != right) { + idmap[right] = idmap[left]; + right--; + } + } + idmap[left] = pivot; + pivot = left; + left = l_hold; + right = r_hold; + + if (left < pivot) + id_sort(idmap, left, pivot-1); + if (right > pivot) + id_sort(idmap, pivot+1, right); +} + +/***************************************************************************/ + +/* NOTE: the following code is the based on the example implementation + * of the IMD protocol API from VMD and NAMD. The UIUC license allows + * to re-use up to 10% of a project's code to be used in other software */ + +/*************************************************************************** + * DESCRIPTION: + * Socket interface, abstracts machine dependent APIs/routines. + ***************************************************************************/ + +int imdsock_init(void) { +#if defined(_MSC_VER) + int rc = 0; + static int initialized=0; + + if (!initialized) { + WSADATA wsdata; + rc = WSAStartup(MAKEWORD(1,1), &wsdata); + if (rc == 0) + initialized = 1; + } + + return rc; +#else + return 0; +#endif +} + + +void * imdsock_create(void) { + imdsocket * s; + + s = (imdsocket *) malloc(sizeof(imdsocket)); + if (s != NULL) + memset(s, 0, sizeof(imdsocket)); + + if ((s->sd = socket(PF_INET, SOCK_STREAM, 0)) == -1) { + printf("Failed to open socket."); + free(s); + return NULL; + } + + return (void *) s; +} + +int imdsock_bind(void * v, int port) { + imdsocket *s = (imdsocket *) v; + memset(&(s->addr), 0, sizeof(s->addr)); + s->addr.sin_family = PF_INET; + s->addr.sin_port = htons(port); + + return bind(s->sd, (struct sockaddr *) &s->addr, sizeof(s->addr)); +} + +int imdsock_listen(void * v) { + imdsocket *s = (imdsocket *) v; + return listen(s->sd, 5); +} + +void *imdsock_accept(void * v) { + int rc; + imdsocket *new_s = NULL, *s = (imdsocket *) v; +#if defined(ARCH_AIX5) || defined(ARCH_AIX5_64) || defined(ARCH_AIX6_64) + unsigned int len; +#elif defined(SOCKLEN_T) + SOCKLEN_T len; +#elif defined(_POSIX_SOURCE) + socklen_t len; +#else + int len; +#endif + + len = sizeof(s->addr); + rc = accept(s->sd, (struct sockaddr *) &s->addr, &len); + if (rc >= 0) { + new_s = (imdsocket *) malloc(sizeof(imdsocket)); + if (new_s != NULL) { + *new_s = *s; + new_s->sd = rc; + } + } + return (void *)new_s; +} + +int imdsock_write(void * v, const void *buf, int len) { + imdsocket *s = (imdsocket *) v; +#if defined(_MSC_VER) + return send(s->sd, (const char*) buf, len, 0); /* windows lacks the write() call */ +#else + return write(s->sd, buf, len); +#endif +} + +int imdsock_read(void * v, void *buf, int len) { + imdsocket *s = (imdsocket *) v; +#if defined(_MSC_VER) + return recv(s->sd, (char*) buf, len, 0); /* windows lacks the read() call */ +#else + return read(s->sd, buf, len); +#endif + +} + +void imdsock_shutdown(void *v) { + imdsocket * s = (imdsocket *) v; + if (s == NULL) + return; + +#if defined(_MSC_VER) + shutdown(s->sd, SD_SEND); +#else + shutdown(s->sd, 1); /* complete sends and send FIN */ +#endif +} + +void imdsock_destroy(void * v) { + imdsocket * s = (imdsocket *) v; + if (s == NULL) + return; + +#if defined(_MSC_VER) + closesocket(s->sd); +#else + close(s->sd); +#endif + free(s); +} + +int imdsock_selread(void *v, int sec) { + imdsocket *s = (imdsocket *)v; + fd_set rfd; + struct timeval tv; + int rc; + + if (v == NULL) return 0; + + FD_ZERO(&rfd); + FD_SET(s->sd, &rfd); + memset((void *)&tv, 0, sizeof(struct timeval)); + tv.tv_sec = sec; + do { + rc = select(s->sd+1, &rfd, NULL, NULL, &tv); + } while (rc < 0 && errno == EINTR); + return rc; + +} + +int imdsock_selwrite(void *v, int sec) { + imdsocket *s = (imdsocket *)v; + fd_set wfd; + struct timeval tv; + int rc; + + if (v == NULL) return 0; + + FD_ZERO(&wfd); + FD_SET(s->sd, &wfd); + memset((void *)&tv, 0, sizeof(struct timeval)); + tv.tv_sec = sec; + do { + rc = select(s->sd + 1, NULL, &wfd, NULL, &tv); + } while (rc < 0 && errno == EINTR); + return rc; +} + +/* end of socket code. */ +/*************************************************************************/ + +/*************************************************************************/ +/* start of imd API code. */ +/* Only works with aligned 4-byte quantities, will cause a bus error */ +/* on some platforms if used on unaligned data. */ +void swap4_aligned(void *v, long ndata) { + int *data = (int *) v; + long i; + int *N; + for (i=0; i>24)&0xff) | ((*N&0xff)<<24) | + ((*N>>8)&0xff00) | ((*N&0xff00)<<8)); + } +} + + +/** structure used to perform byte swapping operations */ +typedef union { + int32 i; + struct { + unsigned int highest : 8; + unsigned int high : 8; + unsigned int low : 8; + unsigned int lowest : 8; + } b; +} netint; + +static int32 imd_htonl(int32 h) { + netint n; + n.b.highest = h >> 24; + n.b.high = h >> 16; + n.b.low = h >> 8; + n.b.lowest = h; + return n.i; +} + +static int32 imd_ntohl(int32 n) { + netint u; + u.i = n; + return (u.b.highest << 24 | u.b.high << 16 | u.b.low << 8 | u.b.lowest); +} + +static void imd_fill_header(IMDheader *header, IMDType type, int32 length) { + header->type = imd_htonl((int32)type); + header->length = imd_htonl(length); +} + +static void swap_header(IMDheader *header) { + header->type = imd_ntohl(header->type); + header->length= imd_ntohl(header->length); +} + +static int32 imd_readn(void *s, char *ptr, int32 n) { + int32 nleft; + int32 nread; + + nleft = n; + while (nleft > 0) { + if ((nread = imdsock_read(s, ptr, nleft)) < 0) { + if (errno == EINTR) + nread = 0; /* and call read() again */ + else + return -1; + } else if (nread == 0) + break; /* EOF */ + nleft -= nread; + ptr += nread; + } + return n-nleft; +} + +static int32 imd_writen(void *s, const char *ptr, int32 n) { + int32 nleft; + int32 nwritten; + + nleft = n; + while (nleft > 0) { + if ((nwritten = imdsock_write(s, ptr, nleft)) <= 0) { + if (errno == EINTR) + nwritten = 0; + else + return -1; + } + nleft -= nwritten; + ptr += nwritten; + } + return n; +} + +int imd_disconnect(void *s) { + IMDheader header; + imd_fill_header(&header, IMD_DISCONNECT, 0); + return (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE); +} + +int imd_handshake(void *s) { + IMDheader header; + imd_fill_header(&header, IMD_HANDSHAKE, 1); + header.length = IMDVERSION; /* Not byteswapped! */ + return (imd_writen(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE); +} + +/* The IMD receive functions */ + +IMDType imd_recv_header(void *s, int32 *length) { + IMDheader header; + if (imd_readn(s, (char *)&header, IMDHEADERSIZE) != IMDHEADERSIZE) + return IMD_IOERROR; + swap_header(&header); + *length = header.length; + return IMDType(header.type); +} + +int imd_recv_mdcomm(void *s, int32 n, int32 *indices, float *forces) { + if (imd_readn(s, (char *)indices, 4*n) != 4*n) return 1; + if (imd_readn(s, (char *)forces, 12*n) != 12*n) return 1; + return 0; +} + +int imd_recv_energies(void *s, IMDEnergies *energies) { + return (imd_readn(s, (char *)energies, sizeof(IMDEnergies)) + != sizeof(IMDEnergies)); +} + +int imd_recv_fcoords(void *s, int32 n, float *coords) { + return (imd_readn(s, (char *)coords, 12*n) != 12*n); +} + +/************************************************************************ + * integer hash code: + ************************************************************************/ + +/* inthash() - Hash function returns a hash number for a given key. + * tptr: Pointer to a hash table, key: The key to create a hash number for */ +static int inthash(const inthash_t *tptr, int key) { + int hashvalue; + + hashvalue = (((key*1103515249)>>tptr->downshift) & tptr->mask); + if (hashvalue < 0) { + hashvalue = 0; + } + + return hashvalue; +} + +/* + * rebuild_table_int() - Create new hash table when old one fills up. + * + * tptr: Pointer to a hash table + */ +static void rebuild_table_int(inthash_t *tptr) { + inthash_node_t **old_bucket, *old_hash, *tmp; + int old_size, h, i; + + old_bucket=tptr->bucket; + old_size=tptr->size; + + /* create a new table and rehash old buckets */ + inthash_init(tptr, old_size<<1); + for (i=0; inext; + h=inthash(tptr, tmp->key); + tmp->next=tptr->bucket[h]; + tptr->bucket[h]=tmp; + tptr->entries++; + } /* while */ + } /* for */ + + /* free memory used by old table */ + free(old_bucket); + + return; +} + +/* + * inthash_init() - Initialize a new hash table. + * + * tptr: Pointer to the hash table to initialize + * buckets: The number of initial buckets to create + */ +void inthash_init(inthash_t *tptr, int buckets) { + + /* make sure we allocate something */ + if (buckets==0) + buckets=16; + + /* initialize the table */ + tptr->entries=0; + tptr->size=2; + tptr->mask=1; + tptr->downshift=29; + + /* ensure buckets is a power of 2 */ + while (tptr->sizesize<<=1; + tptr->mask=(tptr->mask<<1)+1; + tptr->downshift--; + } /* while */ + + /* allocate memory for table */ + tptr->bucket=(inthash_node_t **) calloc(tptr->size, sizeof(inthash_node_t *)); + + return; +} + +/* + * inthash_lookup() - Lookup an entry in the hash table and return a pointer to + * it or HASH_FAIL if it wasn't found. + * + * tptr: Pointer to the hash table + * key: The key to lookup + */ +int inthash_lookup(const inthash_t *tptr, int key) { + int h; + inthash_node_t *node; + + + /* find the entry in the hash table */ + h=inthash(tptr, key); + for (node=tptr->bucket[h]; node!=NULL; node=node->next) { + if (node->key == key) + break; + } + + /* return the entry if it exists, or HASH_FAIL */ + return(node ? node->data : HASH_FAIL); +} + + +/* + * inthash_keys() - Return a list of keys. + * NOTE: the returned list must be freed with free(3). + */ +int *inthash_keys(inthash_t *tptr) { + + int *keys; + inthash_node_t *node; + + keys = (int *)calloc(tptr->entries, sizeof(int)); + + for (int i=0; i < tptr->size; ++i) { + for (node=tptr->bucket[i]; node != NULL; node=node->next) { + keys[node->data] = node->key; + } + } + + return keys; +} + +/* + * inthash_insert() - Insert an entry into the hash table. If the entry already + * exists return a pointer to it, otherwise return HASH_FAIL. + * + * tptr: A pointer to the hash table + * key: The key to insert into the hash table + * data: A pointer to the data to insert into the hash table + */ +int inthash_insert(inthash_t *tptr, int key, int data) { + int tmp; + inthash_node_t *node; + int h; + + /* check to see if the entry exists */ + if ((tmp=inthash_lookup(tptr, key)) != HASH_FAIL) + return(tmp); + + /* expand the table if needed */ + while (tptr->entries>=HASH_LIMIT*tptr->size) + rebuild_table_int(tptr); + + /* insert the new entry */ + h=inthash(tptr, key); + node=(struct inthash_node_t *) malloc(sizeof(inthash_node_t)); + node->data=data; + node->key=key; + node->next=tptr->bucket[h]; + tptr->bucket[h]=node; + tptr->entries++; + + return HASH_FAIL; +} + +/* + * inthash_destroy() - Delete the entire table, and all remaining entries. + * + */ +void inthash_destroy(inthash_t *tptr) { + inthash_node_t *node, *last; + int i; + + for (i=0; isize; i++) { + node = tptr->bucket[i]; + while (node != NULL) { + last = node; + node = node->next; + free(last); + } + } + + /* free the entire array of buckets */ + if (tptr->bucket != NULL) { + free(tptr->bucket); + memset(tptr, 0, sizeof(inthash_t)); + } +} + + +// Local Variables: +// mode: c++ +// compile-command: "make -j4 openmpi" +// c-basic-offset: 2 +// fill-column: 76 +// indent-tabs-mode: nil +// End: + diff --git a/src/USER-IMD/fix_imd.h b/src/USER-IMD/fix_imd.h new file mode 100644 index 0000000000..2c0d0213e1 --- /dev/null +++ b/src/USER-IMD/fix_imd.h @@ -0,0 +1,99 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + The FixIMD class contains code from VMD and NAMD which is copyrighted + by the Board of Trustees of the University of Illinois and is free to + use with LAMMPS according to point 2 of the UIUC license (10% clause): + +" Licensee may, at its own expense, create and freely distribute +complimentary works that interoperate with the Software, directing others to +the TCBG server to license and obtain the Software itself. Licensee may, at +its own expense, modify the Software to make derivative works. Except as +explicitly provided below, this License shall apply to any derivative work +as it does to the original Software distributed by Illinois. Any derivative +work should be clearly marked and renamed to notify users that it is a +modified version and not the original Software distributed by Illinois. +Licensee agrees to reproduce the copyright notice and other proprietary +markings on any derivative work and to include in the documentation of such +work the acknowledgement: + + "This software includes code developed by the Theoretical and Computational + Biophysics Group in the Beckman Institute for Advanced Science and + Technology at the University of Illinois at Urbana-Champaign." + +Licensee may redistribute without restriction works with up to 1/2 of their +non-comment source code derived from at most 1/10 of the non-comment source +code developed by Illinois and contained in the Software, provided that the +above directions for notice and acknowledgement are observed. Any other +distribution of the Software or any derivative work requires a separate +license with Illinois. Licensee may contact Illinois (vmd@ks.uiuc.edu) to +negotiate an appropriate license for such distribution." +------------------------------------------------------------------------- */ + +#ifndef FIX_IMD_H +#define FIX_IMD_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixIMD : public Fix { + public: + FixIMD(class LAMMPS *, int, char **); + ~FixIMD(); + int setmask(); + void init(); + void setup(int); + void post_force(int); + void post_force_respa(int, int, int); + double memory_usage(); + + private: + int imd_port; + void *localsock; + void *clientsock; + + int num_coords; // total number of atoms controlled by this fix + int size_one; // bytes per atom in communication buffer. + int maxbuf; // size of atom communication buffer. + void *comm_buf; // communication buffer + void *idmap; // hash for mapping atom indices to consistent order. + int *rev_idmap; // list of the hash keys for reverse mapping. + + int imd_forces; // number of forces communicated via IMD. + void *force_buf; // force data buffer + double imd_fscale; // scale factor for forces. in case VMD's units are off. + + int imd_inactive; // true if IMD connection stopped. + int imd_terminate; // true if IMD requests termination of run. + int imd_trate; // IMD transmission rate. + + int unwrap_flag; // true if coordinates need to be unwrapped before sending + + int me; // my MPI rank in this "world". + int nlevels_respa; // flag to determine respa levels. +}; + +} + +#endif + +// Local Variables: +// mode: c++ +// compile-command: "make -j4 openmpi" +// c-basic-offset: 2 +// fill-column: 76 +// indent-tabs-mode: nil +// End: + diff --git a/src/USER-IMD/style_user_imd.h b/src/USER-IMD/style_user_imd.h new file mode 100644 index 0000000000..9b44bbdeb4 --- /dev/null +++ b/src/USER-IMD/style_user_imd.h @@ -0,0 +1,20 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FixInclude +#include "fix_imd.h" +#endif + +#ifdef FixClass +FixStyle(imd,FixIMD) +#endif diff --git a/src/style_user_packages.h b/src/style_user_packages.h index 70c68719e0..87603c7b9b 100644 --- a/src/style_user_packages.h +++ b/src/style_user_packages.h @@ -19,4 +19,5 @@ #include "style_user_cd_eam.h" #include "style_user_cg_cmm.h" #include "style_user_ewaldn.h" +#include "style_user_imd.h" #include "style_user_smd.h"