set(PKG_CFLAGS "${PKG_CFLAGS} -DGMX_SOFTWARE_INVSQRT")
endif()
+if(WIN32 AND NOT CYGWIN)
+ set(GMX_WSOCKLIB_PATH CACHE PATH "Path to winsock (wsock32.lib) library.")
+ mark_as_advanced(GMX_WSOCKLIB_PATH)
+ find_library(WSOCK32_LIBRARY NAMES wsock32 PATHS ${GMX_WSOCKLIB_PATH})
+ if(WSOCK32_LIBRARY)
+ list(APPEND GMX_EXTRA_LIBRARIES ${WSOCK32_LIBRARY})
+ add_definitions(-DGMX_HAVE_WINSOCK)
+ else()
+ message(STATUS "No winsock found. Cannot use interactive molecular dynamics (IMD).")
+ endif(WSOCK32_LIBRARY)
+endif()
+
########################################################################
must match the architecture of the machine attempting to use them.
+\section{\normindex{Interactive Molecular Dynamics}}
+{\gromacs} supports the interactive molecular dynamics (IMD) protocol as implemented
+by \href{http://www.ks.uiuc.edu/Research/vmd}{VMD} to control a running simulation
+in NAMD. IMD allows to monitor a running {\gromacs} simulation from a VMD client.
+In addition, the user can interact with the simulation by pulling on atoms, residues
+or fragments with a mouse or a force-feedback device. Additional information about
+the {\gromacs} implementation and an exemplary {\gromacs} IMD system can be found
+\href{http://www.mpibpc.mpg.de/grubmueller/interactivemd}{on this homepage}.
+
+\subsection{Simulation input preparation}
+The {\gromacs} implementation allows transmission and interaction with a part of the
+running simulation only, e.g.\ in cases where no water molecules should be transmitted
+or pulled. The group is specified via the {\tt .mdp} option {\tt IMD-group}. When
+{\tt IMD-group} is empty, the IMD protocol is disabled and cannot be enabled via the
+switches in {\tt mdrun}. To interact with the entire system, {\tt IMD-group} can
+be set to {\tt System}. When using {\tt grompp}, a {\tt .gro} file
+to be used as VMD input is written out ({\tt -imd} switch of {\tt grompp}).
+
+\subsection{Starting the simulation}
+Communication between VMD and {\gromacs} is achieved via TCP sockets and thus enables
+controlling an {\tt mdrun} running locally or on a remote cluster. The port for the
+connection can be specified with the {\tt -imdport} switch of {\tt mdrun}, 8888 is
+the default. If a port number of 0 or smaller is provided, {\gromacs} automatically
+assigns a free port to use with IMD.
+
+Every $N$ steps, the {\tt mdrun} client receives the applied forces from VMD and sends the new
+positions to the client. VMD permits increasing or decreasing the communication frequency
+interactively.
+By default, the simulation starts and runs even if no IMD client is connected. This
+behavior is changed by the {\tt -imdwait} switch of {\tt mdrun}. After startup and
+whenever the client has disconnected, the integration stops until reconnection of
+the client.
+When the {\tt -imdterm} switch is used, the simulation can be terminated by pressing
+the stop button in VMD. This is disabled by default.
+Finally, to allow interacting with the simulation (i.e. pulling from VMD) the {\tt -imdpull}
+switch has to be used.
+Therefore, a simulation can only be monitored but not influenced from the VMD client
+when none of {\tt -imdwait}, {\tt -imdterm} or {\tt -imdpull} are set. However, since
+the IMD protocol requires no authentication, it is not advisable to run simulations on
+a host directly reachable from an insecure environment. Secure shell forwarding of TCP
+can be used to connect to running simulations not directly reachable from the interacting host.
+Note that the IMD command line switches of {\tt mdrun} are hidden by default and show
+up in the help text only with {\tt gmx mdrun -h -hidden}.
+
+\subsection{Connecting from VMD}
+In VMD, first the structure corresponding to the IMD group has to be loaded ({\it File
+$\rightarrow$ New Molecule}). Then the IMD connection window has to be used
+({\it Extensions $\rightarrow$ Simulation $\rightarrow$ IMD Connect (NAMD)}). In the IMD
+connection window, hostname and port have to be specified and followed by pressing
+{\it Connect}. {\it Detach Sim} allows disconnecting without terminating the simulation, while
+{\it Stop Sim} ends the simulation on the next neighbor searching step (if allowed by
+{\tt -imdterm}).
+
+The timestep transfer rate allows adjusting the communication frequency between simulation
+and IMD client. Setting the keep rate loads every $N^\mathrm{th}$ frame into VMD instead
+of discarding them when a new one is received. The displayed energies are in SI units
+in contrast to energies displayed from NAMD simulations.
+
+
% LocalWords: PMF pmf kJ mol Jarzynski BT bilayer rup mdp AFM fepmf fecalc rb
% LocalWords: posre grompp fs Verlet dihedrals hydrogens hydroxyl sulfhydryl
% LocalWords: vsitehydro planarity chirality pdb gmx virtualize virtualized xz
add_subdirectory(essentialdynamics)
add_subdirectory(pulling)
add_subdirectory(simd)
+add_subdirectory(imd)
if (NOT GMX_BUILD_MDRUN_ONLY)
add_subdirectory(legacyheaders)
add_subdirectory(gmxana)
#define MDOF_F (1<<2)
#define MDOF_X_COMPRESSED (1<<3)
#define MDOF_CPT (1<<4)
+#define MDOF_IMD (1<<5)
+
#endif /* GMX_FILEIO_MDOUTF_H */
* below that does the right thing according to the value of
* file_version. */
enum tpxv {
- tpxv_ComputationalElectrophysiology = 96, /**< support for ion/water position swaps (computational electrophysiology) */
- tpxv_Use64BitRandomSeed, /**< change ld_seed from int to gmx_int64_t */
- tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials /**< potentials for supporting coarse-grained force fields */
+ tpxv_ComputationalElectrophysiology = 96, /**< support for ion/water position swaps (computational electrophysiology) */
+ tpxv_Use64BitRandomSeed, /**< change ld_seed from int to gmx_int64_t */
+ tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, /**< potentials for supporting coarse-grained force fields */
+ tpxv_InteractiveMolecularDynamics /**< interactive molecular dynamics (IMD) */
};
/*! \brief Version number of the file format written to run input
*
* When developing a feature branch that needs to change the run input
* file format, change tpx_tag instead. */
-static const int tpx_version = tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials;
+static const int tpx_version = tpxv_InteractiveMolecularDynamics;
+
/* This number should only be increased when you edit the TOPOLOGY section
* or the HEADER of the tpx format.
}
}
+static void do_imd(t_fileio *fio, t_IMD *imd, gmx_bool bRead)
+{
+ gmx_fio_do_int(fio, imd->nat);
+ if (bRead)
+ {
+ snew(imd->ind, imd->nat);
+ }
+ gmx_fio_ndo_int(fio, imd->ind, imd->nat);
+}
+
static void do_fepvals(t_fileio *fio, t_lambda *fepvals, gmx_bool bRead, int file_version)
{
/* i is defined in the ndo_double macro; use g to iterate. */
ir->bRot = FALSE;
}
+ /* Interactive molecular dynamics */
+ if (file_version >= tpxv_InteractiveMolecularDynamics)
+ {
+ gmx_fio_do_int(fio, ir->bIMD);
+ if (TRUE == ir->bIMD)
+ {
+ if (bRead)
+ {
+ snew(ir->imd, 1);
+ }
+ do_imd(fio, ir->imd, bRead);
+ }
+ }
+ else
+ {
+ /* We don't support IMD sessions for old .tpr files */
+ ir->bIMD = FALSE;
+ }
+
/* grpopts stuff */
gmx_fio_do_int(fio, ir->opts.ngtc);
if (file_version >= 69)
nblock_bc(cr, adress->n_energy_grps, adress->group_explicit);
}
}
+
+static void bc_imd(const t_commrec *cr, t_IMD *imd)
+{
+ int g;
+
+ block_bc(cr, *imd);
+ snew_bc(cr, imd->ind, imd->nat);
+ nblock_bc(cr, imd->nat, imd->ind);
+}
+
static void bc_fepvals(const t_commrec *cr, t_lambda *fep)
{
gmx_bool bAlloc = TRUE;
snew_bc(cr, inputrec->rot, 1);
bc_rot(cr, inputrec->rot);
}
+ if (inputrec->bIMD)
+ {
+ snew_bc(cr, inputrec->imd, 1);
+ bc_imd(cr, inputrec->imd);
+ }
for (i = 0; (i < DIM); i++)
{
bc_cosines(cr, &(inputrec->ex[i]));
}
+static void pr_imd(FILE *fp, int indent, t_IMD *imd)
+{
+ PI("IMD_atoms", imd->nat);
+ pr_ivec_block(fp, indent, "atom", imd->ind, imd->nat, TRUE);
+}
+
+
void pr_inputrec(FILE *fp, int indent, const char *title, t_inputrec *ir,
gmx_bool bMDPformat)
{
pr_rot(fp, indent, ir->rot);
}
+ PS("interactiveMD", EBOOL(ir->bIMD));
+ if (ir->bIMD)
+ {
+ pr_imd(fp, indent, ir->imd);
+ }
+
PS("disre", EDISRETYPE(ir->eDisre));
PS("disre-weighting", EDISREWEIGHTING(ir->eDisreWeighting));
PS("disre-mixed", EBOOL(ir->bDisreMixed));
#include "mtop_util.h"
#include "genborn.h"
#include "calc_verletbuf.h"
-
#include "tomorse.h"
+#include "gromacs/imd/imd.h"
+
static int rm_interactions(int ifunc, int nrmols, t_molinfo mols[])
{
warninp_t wi;
char warn_buf[STRLEN];
unsigned int useed;
+ t_atoms IMDatoms; /* Atoms to be operated on interactively (IMD) */
t_filenm fnm[] = {
{ efMDP, NULL, NULL, ffREAD },
{ efTPX, "-o", NULL, ffWRITE },
{ efTRN, "-t", NULL, ffOPTRD },
{ efEDR, "-e", NULL, ffOPTRD },
- { efTRN, "-ref", "rotref", ffOPTRW }
+ /* This group is needed by the VMD viewer as the start configuration for IMD sessions: */
+ { efGRO, "-imd", "imdgroup", ffOPTWR },
+ { efTRN, "-ref", "rotref", ffOPTRW }
};
#define NFILE asize(fnm)
done_warning(wi, FARGS);
write_tpx_state(ftp2fn(efTPX, NFILE, fnm), ir, &state, sys);
+ /* Output IMD group, if bIMD is TRUE */
+ write_IMDgroup_to_file(ir->bIMD, ir, &state, sys, NFILE, fnm);
+
done_atomtype(atype);
done_mtop(sys, TRUE);
done_inputrec_strings();
acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
- wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN];
+ wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
+ imd_grp[STRLEN];
char fep_lambda[efptNR][STRLEN];
char lambda_weights[STRLEN];
char **pull_grp;
is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
}
+ /* Interactive MD */
+ ir->bIMD = FALSE;
+ CCTYPE("Group to display and/or manipulate in interactive MD session");
+ STYPE ("IMD-group", is->imd_grp, NULL);
+ if (is->imd_grp[0] != '\0')
+ {
+ snew(ir->imd, 1);
+ ir->bIMD = TRUE;
+ }
+
/* Refinement */
CCTYPE("NMR refinement stuff");
CTYPE ("Distance restraints type: No, Simple or Ensemble");
}
+void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
+{
+ int ig = -1, i;
+
+
+ ig = search_string(IMDgname, grps->nr, gnames);
+ IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
+
+ if (IMDgroup->nat > 0)
+ {
+ fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
+ IMDgname, IMDgroup->nat);
+ snew(IMDgroup->ind, IMDgroup->nat);
+ for (i = 0; i < IMDgroup->nat; i++)
+ {
+ IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
+ }
+ }
+}
+
+
void do_index(const char* mdparin, const char *ndx,
gmx_mtop_t *mtop,
gmx_bool bVerbose,
make_swap_groups(ir->swap, swapgrp, splitgrp0, splitgrp1, solgrp, grps, gnames);
}
+ /* Make indices for IMD session */
+ if (ir->bIMD)
+ {
+ make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
+ }
+
nacc = str_nelem(is->acc, MAXPTR, ptr1);
nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
if (nacg*DIM != nacc)
--- /dev/null
+#
+# This file is part of the GROMACS molecular simulation package.
+#
+# Copyright (c) 2014, by the GROMACS development team, led by
+# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+# and including many others, as listed in the AUTHORS file in the
+# top-level source directory and at http://www.gromacs.org.
+#
+# GROMACS is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public License
+# as published by the Free Software Foundation; either version 2.1
+# of the License, or (at your option) any later version.
+#
+# GROMACS is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with GROMACS; if not, see
+# http://www.gnu.org/licenses, or write to the Free Software Foundation,
+# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+# If you want to redistribute modifications to GROMACS, please
+# consider that scientific software is very special. Version
+# control is crucial - bugs must be traceable. We will be happy to
+# consider code for inclusion in the official distribution, but
+# derived work must not be called official GROMACS. Details are found
+# in the README & COPYING files - if they are missing, get the
+# official version at http://www.gromacs.org.
+#
+# To help us fund GROMACS development, we humbly ask that you cite
+# the research papers on the package. Check out http://www.gromacs.org.
+
+file(GLOB IMD_SOURCES *.cpp *.c)
+set(LIBGROMACS_SOURCES ${LIBGROMACS_SOURCES} ${IMD_SOURCES} PARENT_SCOPE)
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \libinternal \file
+ *
+ * \brief
+ * Implements functions of imd.h.
+ *
+ * Re-implementation of basic IMD functions from NAMD/VMD from scratch,
+ * see imdsocket.h for references to the IMD API.
+ *
+ * \author Martin Hoefling, Carsten Kutzner <ckutzne@gwdg.de>
+ *
+ * \inlibraryapi
+ * \ingroup module_imd
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include "imd.h"
+#include "imdsocket.h"
+#include "smalloc.h"
+#include "network.h"
+#include "mdrun.h"
+#include "sighandler.h"
+#include "gmx_ga2la.h"
+#include "xvgr.h"
+#include "gromacs/mdlib/groupcoord.h"
+#include "gromacs/fileio/confio.h"
+#include "mtop_util.h"
+#include "names.h"
+#include "gromacs/timing/wallcycle.h"
+
+#include <string.h>
+
+#ifdef GMX_NATIVE_WINDOWS
+#include <Windows.h>
+#else
+#include <unistd.h>
+#endif
+
+/*! \brief How long shall we wait in seconds until we check for a connection again? */
+#define IMDLOOPWAIT 1
+
+/*! \brief How long shall we check for the IMD_GO? */
+#define IMDCONNECTWAIT 2
+
+/*! \brief IMD Header Size. */
+#define HEADERSIZE 8
+/*! \brief IMD Protocol Version. */
+#define IMDVERSION 2
+
+/*! \brief Broadcast d to all nodes */
+#define block_bc(cr, d) gmx_bcast( sizeof(d), &(d), (cr))
+
+/*! \brief Broadcast nr elements of d to all nodes */
+#define nblock_bc(cr, nr, d) gmx_bcast((nr)*sizeof((d)[0]), (d), (cr))
+
+
+/*! \libinternal
+ * \brief
+ * IMD (interactive molecular dynamics) energy record.
+ *
+ * As in the original IMD implementation. Energies in kcal/mol.
+ * NOTE: We return the energies in GROMACS / SI units,
+ * so they also show up as SI in VMD.
+ *
+ */
+typedef struct
+{
+ gmx_int32_t tstep; /**< time step */
+ float T_abs; /**< absolute temperature */
+ float E_tot; /**< total energy */
+ float E_pot; /**< potential energy */
+ float E_vdw; /**< van der Waals energy */
+ float E_coul; /**< Coulomb interaction energy */
+ float E_bond; /**< bonds energy */
+ float E_angle; /**< angles energy */
+ float E_dihe; /**< dihedrals energy */
+ float E_impr; /**< improper dihedrals energy */
+} IMDEnergyBlock;
+
+
+/*! \libinternal
+ * \brief IMD (interactive molecular dynamics) communication structure.
+ *
+ * This structure defines the IMD communication message header & protocol version.
+ */
+typedef struct
+{
+ gmx_int32_t type; /**< Type of IMD message, see IMDType_t above */
+ gmx_int32_t length; /**< Length */
+} IMDHeader;
+
+
+/*! \libinternal
+ * \brief IMD (interactive molecular dynamics) main data structure.
+ *
+ * Contains private IMD data
+ */
+typedef struct gmx_IMD
+{
+ FILE *outf; /**< Output file for IMD data, mainly forces. */
+
+ int nat; /**< Number of atoms that can be pulled via IMD. */
+ int nat_loc; /**< Part of the atoms that are local. */
+ atom_id *ind; /**< Global indices of the IMD atoms. */
+ atom_id *ind_loc; /**< Local indices of the IMD atoms. */
+ int nalloc_loc; /**< Allocation size for ind_loc. */
+ rvec *xa; /**< Positions for all IMD atoms assembled on
+ the master node. */
+ ivec *xa_shifts; /**< Shifts for all IMD atoms, to make
+ molecule(s) whole. */
+ ivec *xa_eshifts; /**< Extra shifts since last DD step. */
+ rvec *xa_old; /**< Old positions for all IMD atoms on master. */
+ int *xa_ind; /**< Position of each local atom in the
+ collective array. */
+
+ int nstimd; /**< Global IMD frequency, known to all nodes. */
+ int nstimd_new; /**< New frequency from IMD client, master only. */
+ int nstimd_def; /**< Default IMD frequency when disconnected. */
+
+ int port; /**< Port to use for network socket. */
+ IMDSocket *socket; /**< The IMD socket on the master node. */
+ IMDSocket *clientsocket; /**< The IMD socket on the client. */
+ int length; /**< Length we got with last header. */
+
+ gmx_bool bWConnect; /**< Shall we block and wait for connection? */
+ gmx_bool bTerminated; /**< Set if MD is terminated. */
+ gmx_bool bTerminatable; /**< Set if MD can be terminated. */
+ gmx_bool bConnected; /**< Set if connection is present. */
+ gmx_bool bNewForces; /**< Set if we received new forces. */
+ gmx_bool bForceActivated; /**< Set if pulling from VMD is allowed. */
+
+ IMDEnergyBlock *energies; /**< Pointer to energies we send back. */
+
+ gmx_int32_t vmd_nforces; /**< Number of VMD forces. */
+ gmx_int32_t *vmd_f_ind; /**< VMD forces indices. */
+ float *vmd_forces; /**< The VMD forces flat in memory. */
+ int nforces; /**< Number of actual MD forces;
+ this gets communicated to the clients. */
+ atom_id *f_ind; /**< Force indices. */
+ rvec *f; /**< The IMD pulling forces. */
+
+ char *forcesendbuf; /**< Buffer for force sending. */
+ char *coordsendbuf; /**< Buffer for coordinate sending. */
+ char *energysendbuf; /**< Send buffer for energies. */
+ rvec *sendxbuf; /**< Buffer to make molecules whole before
+ sending. */
+
+ t_block mols; /**< Molecules block in IMD group. */
+
+ /* The next block is used on the master node only to reduce the output
+ * without sacrificing information. If any of these values changes,
+ * we need to write output */
+ int old_nforces; /**< Old value for nforces. */
+ atom_id *old_f_ind; /**< Old values for force indices. */
+ rvec *old_forces; /**< Old values for IMD pulling forces. */
+
+} t_gmx_IMD_setup;
+
+
+/*! \libinternal
+ * \brief Enum for types of IMD messages.
+ *
+ * We use the same records as the NAMD/VMD IMD implementation.
+ */
+typedef enum IMDType_t
+{
+ IMD_DISCONNECT, /**< client disconnect */
+ IMD_ENERGIES, /**< energy data */
+ IMD_FCOORDS, /**< atomic coordinates */
+ IMD_GO, /**< start command for the simulation */
+ IMD_HANDSHAKE, /**< handshake to determine little/big endianness */
+ IMD_KILL, /**< terminates the simulation */
+ IMD_MDCOMM, /**< force data */
+ IMD_PAUSE, /**< pauses the simulation */
+ IMD_TRATE, /**< sets the IMD transmission and processing rate */
+ IMD_IOERROR, /**< I/O error */
+ IMD_NR /**< number of entries */
+} IMDMessageType;
+
+
+/*! \libinternal
+ * \brief Names of the IMDType for error messages.
+ */
+const char *eIMDType_names[IMD_NR + 1] = {
+ "IMD_DISCONNECT",
+ "IMD_ENERGIES",
+ "IMD_FCOORDS",
+ "IMD_GO",
+ "IMD_HANDSHAKE",
+ "IMD_KILL",
+ "IMD_MDCOMM",
+ "IMD_PAUSE",
+ "IMD_TRATE",
+ "IMD_IOERROR",
+ NULL
+};
+
+
+#ifdef GMX_IMD
+
+/*! \brief Fills the header with message and the length argument. */
+static void fill_header(IMDHeader *header, IMDMessageType type, gmx_int32_t length)
+{
+ /* We (ab-)use htonl network function for the correct endianness */
+ header->type = htonl((gmx_int32_t) type);
+ header->length = htonl(length);
+}
+
+
+/*! \brief Swaps the endianess of the header. */
+static void swap_header(IMDHeader *header)
+{
+ /* and vice versa... */
+ header->type = ntohl(header->type);
+ header->length = ntohl(header->length);
+}
+
+
+/*! \brief Reads multiple bytes from socket. */
+static gmx_int32_t imd_read_multiple(IMDSocket *socket, char *datptr, gmx_int32_t toread)
+{
+ gmx_int32_t leftcount, countread;
+
+
+ leftcount = toread;
+ /* Try to read while we haven't reached toread */
+ while (leftcount != 0)
+ {
+ if ((countread = imdsock_read(socket, datptr, leftcount)) < 0)
+ {
+ /* interrupted function call, try again... */
+ if (errno == EINTR)
+ {
+ countread = 0;
+ }
+ /* this is an unexpected error, return what we got */
+ else
+ {
+ return toread - leftcount;
+ }
+
+ /* nothing read, finished */
+ }
+ else if (countread == 0)
+ {
+ break;
+ }
+ leftcount -= countread;
+ datptr += countread;
+ } /* end while */
+
+ /* return nr of bytes read */
+ return toread - leftcount;
+}
+
+
+/*! \brief Writes multiple bytes to socket in analogy to imd_read_multiple. */
+static gmx_int32_t imd_write_multiple(IMDSocket *socket, const char *datptr, gmx_int32_t towrite)
+{
+ gmx_int32_t leftcount, countwritten;
+
+
+ leftcount = towrite;
+ while (leftcount != 0)
+ {
+ if ((countwritten = imdsock_write(socket, datptr, leftcount)) <= 0)
+ {
+ if (errno == EINTR)
+ {
+ countwritten = 0;
+ }
+ else
+ {
+ return towrite - leftcount;
+ }
+ }
+ leftcount -= countwritten;
+ datptr += countwritten;
+ } /* end while */
+
+ return towrite - leftcount;
+}
+
+
+/*! \brief Handshake with IMD client. */
+static int imd_handshake(IMDSocket *socket)
+{
+ IMDHeader header;
+
+
+ fill_header(&header, IMD_HANDSHAKE, 1);
+ header.length = IMDVERSION; /* client wants unswapped version */
+
+ return (imd_write_multiple(socket, (char *) &header, HEADERSIZE) != HEADERSIZE);
+}
+
+
+/*! \brief Send energies using the energy block and the send buffer. */
+static int imd_send_energies(IMDSocket *socket, const IMDEnergyBlock *energies, char *buffer)
+{
+ gmx_int32_t recsize;
+
+
+ recsize = HEADERSIZE + sizeof(IMDEnergyBlock);
+ fill_header((IMDHeader *) buffer, IMD_ENERGIES, 1);
+ memcpy(buffer + HEADERSIZE, energies, sizeof(IMDEnergyBlock));
+
+ return (imd_write_multiple(socket, buffer, recsize) != recsize);
+}
+
+
+/*! \brief Receive IMD header from socket, sets the length and returns the IMD message. */
+static IMDMessageType imd_recv_header(IMDSocket *socket, gmx_int32_t *length)
+{
+ IMDHeader header;
+
+
+ if (imd_read_multiple(socket, (char *) &header, HEADERSIZE) != HEADERSIZE)
+ {
+ return IMD_IOERROR;
+ }
+ swap_header(&header);
+ *length = header.length;
+
+ return (IMDMessageType) header.type;
+}
+
+
+/*! \brief Receive force indices and forces.
+ *
+ * The number of forces was previously communicated via the header.
+ */
+static int imd_recv_mdcomm(IMDSocket *socket, gmx_int32_t nforces, gmx_int32_t *forcendx, float *forces)
+{
+ int retsize, retbytes;
+
+
+ /* reading indices */
+ retsize = sizeof(gmx_int32_t) * nforces;
+ retbytes = imd_read_multiple(socket, (char *) forcendx, retsize);
+ if (retbytes != retsize)
+ {
+ return FALSE;
+ }
+
+ /* reading forces as float array */
+ retsize = 3 * sizeof(float) * nforces;
+ retbytes = imd_read_multiple(socket, (char *) forces, retsize);
+ if (retbytes != retsize)
+ {
+ return FALSE;
+ }
+
+ return TRUE;
+}
+
+#endif
+
+/* GROMACS specific functions for the IMD implementation */
+
+extern void write_IMDgroup_to_file(gmx_bool bIMD, t_inputrec *ir, t_state *state,
+ gmx_mtop_t *sys, int nfile, const t_filenm fnm[])
+{
+ t_atoms IMDatoms;
+
+
+ if (bIMD)
+ {
+ IMDatoms = gmx_mtop_global_atoms(sys);
+ write_sto_conf_indexed(opt2fn("-imd", nfile, fnm), "IMDgroup", &IMDatoms,
+ state->x, state->v, ir->ePBC, state->box, ir->imd->nat, ir->imd->ind);
+ }
+}
+
+
+extern void dd_make_local_IMD_atoms(gmx_bool bIMD, gmx_domdec_t *dd, t_IMD *imd)
+{
+ gmx_ga2la_t ga2la;
+ t_gmx_IMD_setup *IMDsetup;
+
+
+ if (bIMD)
+ {
+ IMDsetup = imd->setup;
+ ga2la = dd->ga2la;
+
+ dd_make_local_group_indices(
+ ga2la, IMDsetup->nat, IMDsetup->ind, &IMDsetup->nat_loc,
+ &IMDsetup->ind_loc, &IMDsetup->nalloc_loc, IMDsetup->xa_ind);
+ }
+}
+
+
+#ifdef GMX_IMD
+/*! \brief Send positions from rvec.
+ *
+ * We need a separate send buffer and conversion to Angstrom.
+ */
+static int imd_send_rvecs(IMDSocket *socket, int nat, rvec *x, char *buffer)
+{
+ gmx_int32_t size;
+ int i;
+ float sendx[3];
+ int tuplesize = 3 * sizeof(float);
+
+
+ /* Required size for the send buffer */
+ size = HEADERSIZE + 3 * sizeof(float) * nat;
+
+ /* Prepare header */
+ fill_header((IMDHeader *) buffer, IMD_FCOORDS, (gmx_int32_t) nat);
+ for (i = 0; i < nat; i++)
+ {
+ sendx[0] = (float) x[i][0] * NM2A;
+ sendx[1] = (float) x[i][1] * NM2A;
+ sendx[2] = (float) x[i][2] * NM2A;
+ memcpy(buffer + HEADERSIZE + i * tuplesize, sendx, tuplesize);
+ }
+
+ return (imd_write_multiple(socket, buffer, size) != size);
+}
+
+
+/*! \brief Initializes the IMD private data. */
+static t_gmx_IMD_setup* imd_create(int imdatoms, int nstimddef, int imdport)
+{
+ t_gmx_IMD_setup *IMDsetup = NULL;
+
+
+ snew(IMDsetup, 1);
+ IMDsetup->nat = imdatoms;
+ IMDsetup->bTerminated = FALSE;
+ IMDsetup->bTerminatable = FALSE;
+ IMDsetup->bWConnect = FALSE;
+ IMDsetup->bConnected = FALSE;
+ IMDsetup->bForceActivated = FALSE;
+ IMDsetup->bNewForces = FALSE;
+ IMDsetup->bForceActivated = FALSE;
+ IMDsetup->nstimd = 1;
+ IMDsetup->nstimd_new = 1;
+ IMDsetup->nstimd_def = nstimddef;
+ if (imdport < 1)
+ {
+ IMDsetup->port = 0;
+ fprintf(stderr, "%s You chose a port number < 1. Will automatically assign a free port.\n", IMDstr);
+ }
+ else
+ {
+ IMDsetup->port = imdport;
+ }
+
+ return IMDsetup;
+}
+
+
+/*! \brief Prepare the socket on the MASTER. */
+static void imd_prepare_master_socket(t_gmx_IMD_setup *IMDsetup)
+{
+ int ret;
+
+
+#ifdef GMX_NATIVE_WINDOWS
+ /* Winsock requires separate initialization */
+ fprintf(stderr, "%s Initializing winsock.\n", IMDstr);
+#ifdef GMX_HAVE_WINSOCK
+ if (imdsock_winsockinit())
+ {
+ gmx_fatal(FARGS, "%s Failed to initialize winsock.\n", IMDstr);
+ }
+#endif
+#endif
+
+ /* The rest is identical, first create and bind a socket and set to listen then. */
+ fprintf(stderr, "%s Setting up incoming socket.\n", IMDstr);
+ IMDsetup->socket = imdsock_create();
+ if (!IMDsetup->socket)
+ {
+ gmx_fatal(FARGS, "%s Failed to create socket.", IMDstr);
+ }
+
+ /* Bind to port */
+ ret = imdsock_bind(IMDsetup->socket, IMDsetup->port);
+ if (ret)
+ {
+ gmx_fatal(FARGS, "%s binding socket to port %d failed with error %d.\n", IMDstr, IMDsetup->port, ret);
+ }
+
+ if (imd_sock_listen(IMDsetup->socket))
+ {
+ gmx_fatal(FARGS, "%s socket listen failed with error %d.\n", IMDstr, ret);
+ }
+
+ if (imdsock_getport(IMDsetup->socket, &IMDsetup->port))
+ {
+ gmx_fatal(FARGS, "%s Could not determine port number.\n", IMDstr, ret);
+ }
+
+ fprintf(stderr, "%s Listening for IMD connection on port %d.\n", IMDstr, IMDsetup->port);
+}
+
+
+/*! \brief Disconnect the client. */
+static void imd_disconnect(t_gmx_IMD_setup *IMDsetup)
+{
+ /* Write out any buffered pulling data */
+ fflush(IMDsetup->outf);
+
+ /* we first try to shut down the clientsocket */
+ imdsock_shutdown(IMDsetup->clientsocket);
+ if (!imdsock_destroy(IMDsetup->clientsocket))
+ {
+ fprintf(stderr, "%s Failed to destroy socket.\n", IMDstr);
+ }
+
+ /* then we reset the IMD step to its default, and reset the connection boolean */
+ IMDsetup->nstimd_new = IMDsetup->nstimd_def;
+ IMDsetup->clientsocket = NULL;
+ IMDsetup->bConnected = FALSE;
+}
+
+
+/*! \brief Prints an error message and disconnects the client.
+ *
+ * Does not terminate mdrun!
+ */
+static void imd_fatal(t_gmx_IMD_setup *IMDsetup, const char *msg)
+{
+ fprintf(stderr, "%s %s", IMDstr, msg);
+ imd_disconnect(IMDsetup);
+ fprintf(stderr, "%s disconnected.\n", IMDstr);
+}
+
+
+/*! \brief Check whether we got an incoming connection. */
+static gmx_bool imd_tryconnect(t_gmx_IMD_setup *IMDsetup)
+{
+ if (imdsock_tryread(IMDsetup->socket, 0, 0) > 0)
+ {
+ /* yes, we got something, accept on clientsocket */
+ IMDsetup->clientsocket = imdsock_accept(IMDsetup->socket);
+ if (!IMDsetup->clientsocket)
+ {
+ fprintf(stderr, "%s Accepting the connection on the socket failed.\n", IMDstr);
+ return FALSE;
+ }
+
+ /* handshake with client */
+ if (imd_handshake(IMDsetup->clientsocket))
+ {
+ imd_fatal(IMDsetup, "Connection failed.\n");
+ return FALSE;
+ }
+
+ fprintf(stderr, "%s Connection established, checking if I got IMD_GO orders.\n", IMDstr);
+
+ /* Check if we get the proper "GO" command from client. */
+ if (imdsock_tryread(IMDsetup->clientsocket, IMDCONNECTWAIT, 0) != 1 || imd_recv_header(IMDsetup->clientsocket, &(IMDsetup->length)) != IMD_GO)
+ {
+ imd_fatal(IMDsetup, "No IMD_GO order received. IMD connection failed.\n");
+ }
+
+ /* IMD connected */
+ IMDsetup->bConnected = TRUE;
+
+ return TRUE;
+ }
+
+ return FALSE;
+}
+
+
+/*! \brief Wrap imd_tryconnect in order to make it blocking.
+ *
+ * Used when the simulation should wait for an incoming connection.
+ */
+static void imd_blockconnect(t_gmx_IMD_setup *IMDsetup)
+{
+ /* do not wait for connection, when e.g. ctrl+c is pressed and we will terminate anyways. */
+ if (!(int) gmx_get_stop_condition() == gmx_stop_cond_none)
+ {
+ return;
+ }
+
+ fprintf(stderr, "%s Will wait until I have a connection and IMD_GO orders.\n", IMDstr);
+
+ /* while we have no clientsocket... 2nd part: we should still react on ctrl+c */
+ while ((!IMDsetup->clientsocket) && ((int) gmx_get_stop_condition() == gmx_stop_cond_none))
+ {
+ imd_tryconnect(IMDsetup);
+#ifdef GMX_NATIVE_WINDOWS
+ /* for whatever reason, it is called Sleep on windows */
+ Sleep(IMDLOOPWAIT);
+#else
+ sleep(IMDLOOPWAIT);
+#endif
+ }
+}
+
+
+/*! \brief Make sure that our array holding the forces received via IMD is large enough. */
+static void imd_prepare_vmd_Forces(t_gmx_IMD_setup *IMDsetup)
+{
+ srenew((IMDsetup->vmd_f_ind), IMDsetup->vmd_nforces);
+ srenew((IMDsetup->vmd_forces), 3*IMDsetup->vmd_nforces);
+}
+
+
+/*! \brief Reads forces received via IMD. */
+static void imd_read_vmd_Forces(t_gmx_IMD_setup *IMDsetup)
+{
+ /* the length of the previously received header tells us the nr of forces we will receive */
+ IMDsetup->vmd_nforces = IMDsetup->length;
+ /* prepare the arrays */
+ imd_prepare_vmd_Forces(IMDsetup);
+ /* Now we read the forces... */
+ if (!(imd_recv_mdcomm(IMDsetup->clientsocket, IMDsetup->vmd_nforces, IMDsetup->vmd_f_ind, IMDsetup->vmd_forces)))
+ {
+ imd_fatal(IMDsetup, "Error while reading forces from remote. Disconnecting\n");
+ }
+}
+
+
+/*! \brief Prepares the MD force arrays. */
+static void imd_prepare_MD_Forces(t_gmx_IMD_setup *IMDsetup)
+{
+ srenew((IMDsetup->f_ind), IMDsetup->nforces);
+ srenew((IMDsetup->f ), IMDsetup->nforces);
+}
+
+
+/*! \brief Copy IMD forces to MD forces.
+ *
+ * Do conversion from Cal->Joule and from
+ * Angstrom -> nm and from a pointer array to arrays to 3*N array.
+ */
+static void imd_copyto_MD_Forces(t_gmx_IMD_setup *IMDsetup)
+{
+ int i;
+ real conversion = CAL2JOULE * NM2A;
+
+
+ for (i = 0; i < IMDsetup->nforces; i++)
+ {
+ /* Copy the indices, a copy is important because we may update the incoming forces
+ * whenever we receive new forces while the MD forces are only communicated upon
+ * IMD communication.*/
+ IMDsetup->f_ind[i] = IMDsetup->vmd_f_ind[i];
+
+ /* Convert to rvecs and do a proper unit conversion */
+ IMDsetup->f[i][0] = IMDsetup->vmd_forces[3*i ] * conversion;
+ IMDsetup->f[i][1] = IMDsetup->vmd_forces[3*i + 1] * conversion;
+ IMDsetup->f[i][2] = IMDsetup->vmd_forces[3*i + 2] * conversion;
+ }
+}
+
+
+/*! \brief Return TRUE if any of the forces or indices changed. */
+static gmx_bool bForcesChanged(t_gmx_IMD_setup *IMDsetup)
+{
+ int i;
+
+
+ /* First, check whether the number of pulled atoms changed */
+ if (IMDsetup->nforces != IMDsetup->old_nforces)
+ {
+ return TRUE;
+ }
+
+ /* Second, check whether any of the involved atoms changed */
+ for (i = 0; i < IMDsetup->nforces; i++)
+ {
+ if (IMDsetup->f_ind[i] != IMDsetup->old_f_ind[i])
+ {
+ return TRUE;
+ }
+ }
+
+ /* Third, check whether all forces are the same */
+ for (i = 0; i < IMDsetup->nforces; i++)
+ {
+ if (IMDsetup->f[i][XX] != IMDsetup->old_forces[i][XX])
+ {
+ return TRUE;
+ }
+ if (IMDsetup->f[i][YY] != IMDsetup->old_forces[i][YY])
+ {
+ return TRUE;
+ }
+ if (IMDsetup->f[i][ZZ] != IMDsetup->old_forces[i][ZZ])
+ {
+ return TRUE;
+ }
+ }
+
+ /* All old and new forces are identical! */
+ return FALSE;
+}
+
+
+/*! \brief Fill the old_f_ind and old_forces arrays with the new, old values. */
+static void keep_old_values(t_gmx_IMD_setup *IMDsetup)
+{
+ int i;
+
+
+ IMDsetup->old_nforces = IMDsetup->nforces;
+
+ for (i = 0; i < IMDsetup->nforces; i++)
+ {
+ IMDsetup->old_f_ind[i] = IMDsetup->f_ind[i];
+ copy_rvec(IMDsetup->f[i], IMDsetup->old_forces[i]);
+ }
+}
+
+
+/*! \brief Returns TRUE if any component of the two rvecs differs. */
+static inline gmx_bool rvecs_differ(const rvec v1, const rvec v2)
+{
+ int i;
+
+
+ for (i = 0; i < DIM; i++)
+ {
+ if (v1[i] != v2[i])
+ {
+ return TRUE;
+ }
+ }
+
+ return FALSE;
+}
+
+
+/*! \brief Write the applied pull forces to logfile.
+ *
+ * Call on master only!
+ */
+static void output_imd_forces(t_inputrec *ir, double time)
+{
+ t_gmx_IMD_setup *IMDsetup;
+ int i;
+
+
+ IMDsetup = ir->imd->setup;
+
+ if (bForcesChanged(IMDsetup))
+ {
+ /* Write time and total number of applied IMD forces */
+ fprintf(IMDsetup->outf, "%14.6e%6d", time, IMDsetup->nforces);
+
+ /* Write out the global atom indices of the pulled atoms and the forces itself,
+ * write out a force only if it has changed since the last output */
+ for (i = 0; i < IMDsetup->nforces; i++)
+ {
+ if (rvecs_differ(IMDsetup->f[i], IMDsetup->old_forces[i]))
+ {
+ fprintf(IMDsetup->outf, "%9d", IMDsetup->ind[IMDsetup->f_ind[i]] + 1);
+ fprintf(IMDsetup->outf, "%12.4e%12.4e%12.4e", IMDsetup->f[i][0], IMDsetup->f[i][1], IMDsetup->f[i][2]);
+ }
+ }
+ fprintf(IMDsetup->outf, "\n");
+
+ keep_old_values(IMDsetup);
+ }
+}
+
+
+/*! \brief Synchronize the nodes. */
+static void imd_sync_nodes(t_inputrec *ir, t_commrec *cr, double t)
+{
+ int new_nforces = 0;
+ t_gmx_IMD_setup *IMDsetup;
+ int start, end, i;
+
+
+ IMDsetup = ir->imd->setup;
+
+ /* Notify the other nodes whether we are still connected. */
+ if (PAR(cr))
+ {
+ block_bc(cr, IMDsetup->bConnected);
+ }
+
+ /* ...if not connected, the job is done here. */
+ if (!IMDsetup->bConnected)
+ {
+ return;
+ }
+
+ /* Let the other nodes know whether we got a new IMD synchronization frequency. */
+ if (PAR(cr))
+ {
+ block_bc(cr, IMDsetup->nstimd_new);
+ }
+
+ /* Now we all set the (new) nstimd communication time step */
+ IMDsetup->nstimd = IMDsetup->nstimd_new;
+
+ /* We're done if we don't allow pulling at all */
+ if (!(IMDsetup->bForceActivated))
+ {
+ return;
+ }
+
+ /* OK, let's check if we have received forces which we need to communicate
+ * to the other nodes */
+ if (MASTER(cr))
+ {
+ if (IMDsetup->bNewForces)
+ {
+ new_nforces = IMDsetup->vmd_nforces;
+ }
+ /* make the "new_forces" negative, when we did not receive new ones */
+ else
+ {
+ new_nforces = IMDsetup->vmd_nforces * -1;
+ }
+ }
+
+ /* make new_forces known to the clients */
+ if (PAR(cr))
+ {
+ block_bc(cr, new_nforces);
+ }
+
+ /* When new_natoms < 0 then we know that these are still the same forces
+ * so we don't communicate them, otherwise... */
+ if (new_nforces >= 0)
+ {
+ /* set local VMD and nforces */
+ IMDsetup->vmd_nforces = new_nforces;
+ IMDsetup->nforces = IMDsetup->vmd_nforces;
+
+ /* now everybody knows the number of forces in f_ind, so we can prepare
+ * the target arrays for indices and forces */
+ imd_prepare_MD_Forces(IMDsetup);
+
+ /* we first update the MD forces on the master by converting the VMD forces */
+ if (MASTER(cr))
+ {
+ imd_copyto_MD_Forces(IMDsetup);
+ /* We also write out forces on every update, so that we know which
+ * forces are applied for every step */
+ if (IMDsetup->outf)
+ {
+ output_imd_forces(ir, t);
+ }
+ }
+
+ /* In parallel mode we communicate the to-be-applied forces to the other nodes */
+ if (PAR(cr))
+ {
+ nblock_bc(cr, IMDsetup->nforces, IMDsetup->f_ind);
+ nblock_bc(cr, IMDsetup->nforces, IMDsetup->f );
+ }
+
+ /* done communicating the forces, reset bNewForces */
+ IMDsetup->bNewForces = FALSE;
+ }
+}
+
+
+/*! \brief Reads header from the client and decides what to do. */
+static void imd_readcommand(t_gmx_IMD_setup *IMDsetup)
+{
+ gmx_bool IMDpaused = FALSE;
+ IMDMessageType itype;
+
+
+ while (IMDsetup->clientsocket && (imdsock_tryread(IMDsetup->clientsocket, 0, 0) > 0 || IMDpaused))
+ {
+ itype = imd_recv_header(IMDsetup->clientsocket, &(IMDsetup->length));
+ /* let's see what we got: */
+ switch (itype)
+ {
+ /* IMD asks us to terminate the simulation, check if the user allowed this */
+ case IMD_KILL:
+ if (IMDsetup->bTerminatable)
+ {
+ fprintf(stderr, " %s Terminating connection and running simulation (if supported by integrator).\n", IMDstr);
+ IMDsetup->bTerminated = TRUE;
+ IMDsetup->bWConnect = FALSE;
+ gmx_set_stop_condition(gmx_stop_cond_next);
+ }
+ else
+ {
+ fprintf(stderr, " %s Set -imdterm command line switch to allow mdrun termination from within IMD.\n", IMDstr);
+ }
+
+ break;
+
+ /* the client doen't want to talk to us anymore */
+ case IMD_DISCONNECT:
+ fprintf(stderr, " %s Disconnecting client.\n", IMDstr);
+ imd_disconnect(IMDsetup);
+ break;
+
+ /* we got new forces, read them and set bNewForces flag */
+ case IMD_MDCOMM:
+ imd_read_vmd_Forces(IMDsetup);
+ IMDsetup->bNewForces = TRUE;
+ break;
+
+ /* the client asks us to (un)pause the simulation. So we toggle the IMDpaused state */
+ case IMD_PAUSE:
+ if (IMDpaused)
+ {
+ fprintf(stderr, " %s Un-pause command received.\n", IMDstr);
+ IMDpaused = FALSE;
+ }
+ else
+ {
+ fprintf(stderr, " %s Pause command received.\n", IMDstr);
+ IMDpaused = TRUE;
+ }
+
+ break;
+
+ /* the client sets a new transfer rate, if we get 0, we reset the rate
+ * to the default. VMD filters 0 however */
+ case IMD_TRATE:
+ IMDsetup->nstimd_new = (IMDsetup->length > 0) ? IMDsetup->length : IMDsetup->nstimd_def;
+ fprintf(stderr, " %s Update frequency will be set to %d.\n", IMDstr, IMDsetup->nstimd_new);
+ break;
+
+ /* Catch all rule for the remaining IMD types which we don't expect */
+ default:
+ fprintf(stderr, " %s Received unexpected %s.\n", IMDstr, ENUM_NAME((int)itype, IMD_NR, eIMDType_names));
+ imd_fatal(IMDsetup, "Terminating connection\n");
+ break;
+ } /* end switch */
+ } /* end while */
+}
+
+
+/*! \brief Open IMD output file and write header information.
+ *
+ * Call on master only.
+ */
+static FILE *open_imd_out(
+ const char *fn,
+ t_gmx_IMD_setup *IMDsetup,
+ int nat_total,
+ output_env_t oenv,
+ unsigned long Flags)
+{
+ FILE *fp;
+
+
+ /* Open log file of applied IMD forces if requested */
+ if (fn && oenv)
+ {
+ /* If we append to an existing file, all the header information is already there */
+ if (Flags & MD_APPENDFILES)
+ {
+ fp = gmx_fio_fopen(fn, "a+");
+ }
+ else
+ {
+ fp = gmx_fio_fopen(fn, "w+");
+ if (IMDsetup->nat == nat_total)
+ {
+ fprintf(fp, "# Note that you can select an IMD index group in the .mdp file if a subset of the atoms suffices.\n");
+ }
+
+ xvgr_header(fp, "IMD Pull Forces", "Time (ps)", "# of Forces / Atom IDs / Forces (kJ/mol)", exvggtNONE, oenv);
+
+ fprintf(fp, "# Can display and manipulate %d (of a total of %d) atoms via IMD.\n", IMDsetup->nat, nat_total);
+ fprintf(fp, "# column 1 : time (ps)\n");
+ fprintf(fp, "# column 2 : total number of atoms feeling an IMD pulling force at that time\n");
+ fprintf(fp, "# cols. 3.-6 : global atom number of pulled atom, x-force, y-force, z-force (kJ/mol)\n");
+ fprintf(fp, "# then follow : atom-ID, f[x], f[y], f[z] for more atoms in case the force on multiple atoms is changed simultaneously.\n");
+ fprintf(fp, "# Note that the force on any atom is always equal to the last value for that atom-ID found in the data.\n");
+ fflush(fp);
+ }
+
+ /* To reduce the output file size we remember the old values and output only
+ * when something changed */
+ snew(IMDsetup->old_f_ind, IMDsetup->nat); /* One can never pull on more atoms */
+ snew(IMDsetup->old_forces, IMDsetup->nat);
+
+ return fp;
+ }
+
+ fprintf(stdout, "%s For a log of the IMD pull forces explicitly specify '-if' on the command line.\n"
+ "%s (Not possible with energy minimization.)\n", IMDstr, IMDstr);
+
+ return NULL;
+}
+#endif
+
+
+extern void IMD_finalize(gmx_bool bIMD, t_IMD *imd)
+{
+ if (bIMD)
+ {
+ if (imd->setup->outf)
+ {
+ gmx_fio_fclose(imd->setup->outf);
+ }
+ }
+}
+
+
+#ifdef GMX_IMD
+/*! \brief Creates the molecule start-end position array of molecules in the IMD group. */
+static void init_imd_prepare_mols_in_imdgroup(t_gmx_IMD_setup *IMDsetup, gmx_mtop_t *top_global)
+{
+ int i, ii;
+ int gstart, gend, count;
+ t_block gmols, lmols;
+ int nat;
+ atom_id *ind;
+
+ gmols = top_global->mols;
+ nat = IMDsetup->nat;
+ ind = IMDsetup->ind;
+
+ lmols.nr = 0;
+
+ /* check whether index is sorted */
+ for (i = 0; i < nat-1; i++)
+ {
+ if (ind[i] > ind[i+1])
+ {
+ gmx_fatal(FARGS, "%s IMD index is not sorted. This is currently not supported.\n", IMDstr);
+ }
+ }
+
+ snew(lmols.index, gmols.nr+1);
+ lmols.index[0] = 0;
+
+ for (i = 0; i < gmols.nr; i++)
+ {
+ gstart = gmols.index[i];
+ gend = gmols.index[i+1];
+ count = 0;
+ for (ii = 0; ii < nat; ii++)
+ {
+ if ((ind[ii] >= gstart) && (ind[ii] < gend))
+ {
+ count += 1;
+ }
+ }
+ if (count > 0)
+ {
+ lmols.index[lmols.nr+1] = lmols.index[lmols.nr]+count;
+ lmols.nr += 1;
+ }
+ }
+ srenew(lmols.index, lmols.nr+1);
+ lmols.nalloc_index = lmols.nr+1;
+ IMDsetup->mols = lmols;
+}
+
+
+/*! \brief Copied and modified from groupcoord.c shift_positions_group(). */
+static void shift_positions(
+ matrix box,
+ rvec x[], /* The positions [0..nr] */
+ ivec is, /* The shift [0..nr] */
+ int nr) /* The number of positions */
+{
+ int i, tx, ty, tz;
+
+ /* Loop over the group's atoms */
+ if (TRICLINIC(box))
+ {
+ for (i = 0; i < nr; i++)
+ {
+ tx = is[XX];
+ ty = is[YY];
+ tz = is[ZZ];
+
+ x[i][XX] = x[i][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
+ x[i][YY] = x[i][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
+ x[i][ZZ] = x[i][ZZ]-tz*box[ZZ][ZZ];
+ }
+ }
+ else
+ {
+ for (i = 0; i < nr; i++)
+ {
+ tx = is[XX];
+ ty = is[YY];
+ tz = is[ZZ];
+
+ x[i][XX] = x[i][XX]-tx*box[XX][XX];
+ x[i][YY] = x[i][YY]-ty*box[YY][YY];
+ x[i][ZZ] = x[i][ZZ]-tz*box[ZZ][ZZ];
+ }
+ }
+}
+
+
+/*! \brief Removes shifts of molecules diffused outside of the box. */
+static void imd_remove_molshifts(t_gmx_IMD_setup *IMDsetup, matrix box)
+{
+ int i, ii, molsize;
+ ivec largest, smallest, shift;
+ t_block mols;
+
+
+ mols = IMDsetup->mols;
+
+ /* for each molecule also present in IMD group */
+ for (i = 0; i < mols.nr; i++)
+ {
+ /* first we determine the minimum and maximum shifts for each molecule */
+
+ clear_ivec(largest);
+ clear_ivec(smallest);
+ clear_ivec(shift);
+
+ copy_ivec(IMDsetup->xa_shifts[mols.index[i]], largest);
+ copy_ivec(IMDsetup->xa_shifts[mols.index[i]], smallest);
+
+ for (ii = mols.index[i]+1; ii < mols.index[i+1]; ii++)
+ {
+ if (IMDsetup->xa_shifts[ii][XX] > largest[XX])
+ {
+ largest[XX] = IMDsetup->xa_shifts[ii][XX];
+ }
+ if (IMDsetup->xa_shifts[ii][XX] < smallest[XX])
+ {
+ smallest[XX] = IMDsetup->xa_shifts[ii][XX];
+ }
+
+ if (IMDsetup->xa_shifts[ii][YY] > largest[YY])
+ {
+ largest[YY] = IMDsetup->xa_shifts[ii][YY];
+ }
+ if (IMDsetup->xa_shifts[ii][YY] < smallest[YY])
+ {
+ smallest[YY] = IMDsetup->xa_shifts[ii][YY];
+ }
+
+ if (IMDsetup->xa_shifts[ii][ZZ] > largest[ZZ])
+ {
+ largest[ZZ] = IMDsetup->xa_shifts[ii][ZZ];
+ }
+ if (IMDsetup->xa_shifts[ii][ZZ] < smallest[ZZ])
+ {
+ smallest[ZZ] = IMDsetup->xa_shifts[ii][ZZ];
+ }
+
+ }
+
+ /* check if we what we can subtract/add to the positions
+ * to put them back into the central box */
+ if (smallest[XX] > 0)
+ {
+ shift[XX] = smallest[XX];
+ }
+ if (smallest[YY] > 0)
+ {
+ shift[YY] = smallest[YY];
+ }
+ if (smallest[ZZ] > 0)
+ {
+ shift[ZZ] = smallest[ZZ];
+ }
+
+ if (largest[XX] < 0)
+ {
+ shift[XX] = largest[XX];
+ }
+ if (largest[YY] < 0)
+ {
+ shift[YY] = largest[YY];
+ }
+ if (largest[ZZ] < 0)
+ {
+ shift[ZZ] = largest[ZZ];
+ }
+
+ /* is there a shift at all? */
+ if ((shift[XX]) || (shift[YY]) || (shift[ZZ]))
+ {
+ molsize = mols.index[i+1]-mols.index[i];
+ /* shift the positions */
+ shift_positions(box, &(IMDsetup->xa[mols.index[i]]), shift, molsize);
+ }
+
+ }
+}
+
+
+/*! \brief Initialize arrays used to assemble the positions from the other nodes. */
+static void init_imd_prepare_for_x_assembly(t_commrec *cr, rvec x[], t_gmx_IMD_setup *IMDsetup)
+{
+ int i, ii;
+
+
+ snew(IMDsetup->xa, IMDsetup->nat);
+ snew(IMDsetup->xa_ind, IMDsetup->nat);
+ snew(IMDsetup->xa_shifts, IMDsetup->nat);
+ snew(IMDsetup->xa_eshifts, IMDsetup->nat);
+ snew(IMDsetup->xa_old, IMDsetup->nat);
+
+ /* Save the original (whole) set of positions such that later the
+ * molecule can always be made whole again */
+ if (MASTER(cr))
+ {
+ for (i = 0; i < IMDsetup->nat; i++)
+ {
+ ii = IMDsetup->ind[i];
+ copy_rvec(x[ii], IMDsetup->xa_old[i]);
+ }
+ }
+
+ if (!PAR(cr))
+ {
+ IMDsetup->nat_loc = IMDsetup->nat;
+ IMDsetup->ind_loc = IMDsetup->ind;
+
+ /* xa_ind[i] needs to be set to i for serial runs */
+ for (i = 0; i < IMDsetup->nat; i++)
+ {
+ IMDsetup->xa_ind[i] = i;
+ }
+ }
+
+ /* Communicate initial coordinates xa_old to all processes */
+#ifdef GMX_MPI
+ if (PAR(cr))
+ {
+ gmx_bcast(IMDsetup->nat * sizeof(IMDsetup->xa_old[0]), IMDsetup->xa_old, cr);
+ }
+#endif
+}
+#endif
+
+
+/*! \brief Check for non-working integrator / parallel options. */
+static void imd_check_integrator_parallel(t_inputrec *ir, t_commrec *cr)
+{
+ if (PAR(cr))
+ {
+ if (((ir->eI) == eiSteep) || ((ir->eI) == eiCG) || ((ir->eI) == eiLBFGS) || ((ir->eI) == eiNM))
+ {
+ gmx_fatal(FARGS, "%s Energy minimization via steep, CG, lbfgs and nm in parallel is currently not supported by IMD.\n", IMDstr);
+ return;
+ }
+ }
+}
+
+
+extern void init_IMD(
+ t_inputrec *ir,
+ t_commrec *cr,
+ gmx_mtop_t *top_global,
+ FILE *fplog,
+ int defnstimd,
+ rvec x[],
+ int nfile,
+ const t_filenm fnm[],
+ output_env_t oenv,
+ int imdport,
+ unsigned long Flags)
+{
+ int i;
+ int nat_total;
+ t_gmx_IMD_setup *IMDsetup;
+ gmx_int32_t bufxsize;
+ gmx_bool bIMD = FALSE;
+
+
+ /* We will allow IMD sessions only if explicitly enabled in the .tpr file */
+ if (FALSE == ir->bIMD)
+ {
+ return;
+ }
+
+ /* It seems we have a .tpr file that defines an IMD group and thus allows IMD sessions.
+ * Check whether we can actually provide the IMD functionality for this setting: */
+ if (MASTER(cr))
+ {
+ /* Check whether IMD was enabled by one of the command line switches: */
+ if ((Flags & MD_IMDWAIT) || (Flags & MD_IMDTERM) || (Flags & MD_IMDPULL))
+ {
+ /* Multiple simulations or replica exchange */
+ if (MULTISIM(cr))
+ {
+ fprintf(stderr, "%s Cannot use IMD for multiple simulations or replica exchange.\n", IMDstr);
+ }
+ /* OK, IMD seems to be allowed and turned on... */
+ else
+ {
+ fprintf(stderr, "%s Enabled. This simulation will accept incoming IMD connections.\n", IMDstr);
+ bIMD = TRUE;
+ }
+ }
+ else
+ {
+ fprintf(stderr, "%s None of the -imd switches was used.\n"
+ "%s This run will not accept incoming IMD connections\n", IMDstr, IMDstr);
+ }
+ } /* end master only */
+
+ /* Disable IMD if not all the needed functionality is there! */
+#if defined(GMX_NATIVE_WINDOWS) && !defined(GMX_HAVE_WINSOCK)
+ bIMD = FALSE;
+ fprintf(stderr, "Disabling IMD because the winsock library was not found at compile time.\n");
+#endif
+
+ /* Let the other nodes know whether we want IMD */
+ if (PAR(cr))
+ {
+ block_bc(cr, bIMD);
+ }
+ /* ... and update our local inputrec accordingly. */
+ ir->bIMD = bIMD;
+
+ /*... if not we are done.*/
+ if (!ir->bIMD)
+ {
+ return;
+ }
+
+
+ /* check if we're using a sane integrator / parallel combination */
+ imd_check_integrator_parallel(ir, cr);
+
+
+ /*******************************************************/
+ /** From here on we assume that IMD is turned on **/
+ /*******************************************************/
+
+#ifdef GMX_IMD
+ nat_total = top_global->natoms;
+
+ /* Initialize IMD setup structure. If we read in a pre-IMD .tpr file, imd->nat
+ * will be zero. For those cases we transfer _all_ atomic positions */
+ ir->imd->setup = imd_create(ir->imd->nat > 0 ? ir->imd->nat : nat_total,
+ defnstimd, imdport);
+ IMDsetup = ir->imd->setup;
+
+ /* We might need to open an output file for IMD forces data */
+ if (MASTER(cr))
+ {
+ IMDsetup->outf = open_imd_out(opt2fn("-if", nfile, fnm), ir->imd->setup, nat_total, oenv, Flags);
+ }
+
+ /* Make sure that we operate with a valid atom index array for the IMD atoms */
+ if (ir->imd->nat > 0)
+ {
+ /* Point to the user-supplied array of atom numbers */
+ IMDsetup->ind = ir->imd->ind;
+ }
+ else
+ {
+ /* Make a dummy (ind[i] = i) array of all atoms */
+ snew(IMDsetup->ind, nat_total);
+ for (i = 0; i < nat_total; i++)
+ {
+ IMDsetup->ind[i] = i;
+ }
+ }
+
+ /* read environment on master and prepare socket for incoming connections */
+ if (MASTER(cr))
+ {
+ /* we allocate memory for our IMD energy structure */
+ gmx_int32_t recsize = HEADERSIZE + sizeof(IMDEnergyBlock);
+ snew(IMDsetup->energysendbuf, recsize);
+
+ /* Shall we wait for a connection? */
+ if (Flags & MD_IMDWAIT)
+ {
+ IMDsetup->bWConnect = TRUE;
+ fprintf(stderr, "%s Pausing simulation while no IMD connection present (-imdwait).\n", IMDstr);
+ }
+
+ /* Will the IMD clients be able to terminate the simulation? */
+ if (Flags & MD_IMDTERM)
+ {
+ IMDsetup->bTerminatable = TRUE;
+ fprintf(stderr, "%s Allow termination of the simulation from IMD client (-imdterm).\n", IMDstr);
+ }
+
+ /* Is pulling from IMD client allowed? */
+ if (Flags & MD_IMDPULL)
+ {
+ IMDsetup->bForceActivated = TRUE;
+ fprintf(stderr, "%s Pulling from IMD remote is enabled (-imdpull).\n", IMDstr);
+ }
+
+ /* Initialize send buffers with constant size */
+ snew(IMDsetup->sendxbuf, IMDsetup->nat);
+ snew(IMDsetup->energies, 1);
+ bufxsize = HEADERSIZE + 3 * sizeof(float) * IMDsetup->nat;
+ snew(IMDsetup->coordsendbuf, bufxsize);
+ }
+
+ /* do we allow interactive pulling? If so let the other nodes know. */
+ if (PAR(cr))
+ {
+ block_bc(cr, IMDsetup->bForceActivated);
+ }
+
+ /* setup the listening socket on master process */
+ if (MASTER(cr))
+ {
+ fprintf(fplog, "%s Setting port for connection requests to %d.\n", IMDstr, IMDsetup->port);
+ fprintf(stderr, "%s Turning on IMD - port for incoming requests is %d.\n", IMDstr, IMDsetup->port);
+ imd_prepare_master_socket(IMDsetup);
+ /* Wait until we have a connection if specified before */
+ if (IMDsetup->bWConnect)
+ {
+ imd_blockconnect(IMDsetup);
+ }
+ else
+ {
+ fprintf(stderr, "%s -imdwait not set, starting simulation.\n", IMDstr);
+ }
+ }
+ /* Let the other nodes know whether we are connected */
+ imd_sync_nodes(ir, cr, 0);
+
+ /* Initialize arrays used to assemble the positions from the other nodes */
+ init_imd_prepare_for_x_assembly(cr, x, IMDsetup);
+
+ /* Initialize molecule blocks to make them whole later...*/
+ if (MASTER(cr))
+ {
+ init_imd_prepare_mols_in_imdgroup(IMDsetup, top_global);
+ }
+#else
+ gmx_incons("init_IMD: this GROMACS version was not compiled with IMD support!");
+#endif
+}
+
+
+extern gmx_bool do_IMD(
+ gmx_bool bIMD,
+ gmx_int64_t step,
+ t_commrec *cr,
+ gmx_bool bNS,
+ matrix box,
+ rvec x[],
+ t_inputrec *ir,
+ double t,
+ gmx_wallcycle_t wcycle)
+{
+ gmx_bool imdstep = FALSE;
+ t_gmx_IMD_setup *IMDsetup;
+
+
+ /* IMD at all? */
+ if (!bIMD)
+ {
+ return FALSE;
+ }
+
+#ifdef GMX_IMD
+ wallcycle_start(wcycle, ewcIMD);
+
+ IMDsetup = ir->imd->setup;
+
+ /* read command from client and check if new incoming connection */
+ if (MASTER(cr))
+ {
+ /* If not already connected, check for new connections */
+ if (!IMDsetup->clientsocket)
+ {
+ if (IMDsetup->bWConnect)
+ {
+ imd_blockconnect(IMDsetup);
+ }
+ else
+ {
+ imd_tryconnect(IMDsetup);
+ }
+ }
+
+ /* Let's see if we have new IMD messages for us */
+ if (IMDsetup->clientsocket)
+ {
+ imd_readcommand(IMDsetup);
+ }
+ }
+
+ /* is this an IMD communication step? */
+ imdstep = do_per_step(step, IMDsetup->nstimd);
+
+ /* OK so this is an IMD step ... */
+ if (imdstep)
+ {
+ /* First we sync all nodes to let everybody know whether we are connected to VMD */
+ imd_sync_nodes(ir, cr, t);
+ }
+
+ /* If a client is connected, we collect the positions
+ * and put molecules back into the box before transfer */
+ if ((imdstep && IMDsetup->bConnected)
+ || bNS) /* independent of imdstep, we communicate positions at each NS step */
+ {
+ /* Transfer the IMD positions to the master node. Every node contributes
+ * its local positions x and stores them in the assembled xa array. */
+ communicate_group_positions(cr, IMDsetup->xa, IMDsetup->xa_shifts, IMDsetup->xa_eshifts,
+ TRUE, x, IMDsetup->nat, IMDsetup->nat_loc,
+ IMDsetup->ind_loc, IMDsetup->xa_ind, IMDsetup->xa_old, box);
+
+ /* If connected and master -> remove shifts */
+ if ((imdstep && IMDsetup->bConnected) && MASTER(cr))
+ {
+ imd_remove_molshifts(IMDsetup, box);
+ }
+ }
+
+ wallcycle_stop(wcycle, ewcIMD);
+#else
+ gmx_incons("do_IMD called without IMD support!");
+#endif
+
+ return imdstep;
+}
+
+
+extern void IMD_fill_energy_record(gmx_bool bIMD, t_IMD *imd, gmx_enerdata_t *enerd,
+ gmx_int64_t step, gmx_bool bHaveNewEnergies)
+{
+ IMDEnergyBlock *ene;
+ t_gmx_IMD IMDsetup;
+
+
+ if (bIMD)
+ {
+#ifdef GMX_IMD
+ IMDsetup = imd->setup;
+
+ if (IMDsetup->clientsocket)
+ {
+ ene = IMDsetup->energies;
+
+ ene->tstep = step;
+
+ /* In MPI-parallel simulations the energies are not accessible a at every time step.
+ * We update them if we have new values, otherwise, the energy values from the
+ * last global communication step are still on display in the viewer. */
+ if (bHaveNewEnergies)
+ {
+ ene->T_abs = (float) enerd->term[F_TEMP ];
+ ene->E_pot = (float) enerd->term[F_EPOT ];
+ ene->E_tot = (float) enerd->term[F_ETOT ];
+ ene->E_bond = (float) enerd->term[F_BONDS ];
+ ene->E_angle = (float) enerd->term[F_ANGLES ];
+ ene->E_dihe = (float) enerd->term[F_PDIHS ];
+ ene->E_impr = (float) enerd->term[F_IDIHS ];
+ ene->E_vdw = (float) enerd->term[F_LJ ];
+ ene->E_coul = (float) (enerd->term[F_COUL_SR] + enerd->term[F_COUL_LR]);
+ }
+ }
+#else
+ gmx_incons("IMD_fill_energy_record called without IMD support.");
+#endif
+ }
+}
+
+
+extern void IMD_send_positions(t_IMD *imd)
+{
+#ifdef GMX_IMD
+ t_gmx_IMD IMDsetup;
+
+
+ IMDsetup = imd->setup;
+
+ if (IMDsetup->clientsocket)
+ {
+
+ if (imd_send_energies(IMDsetup->clientsocket, IMDsetup->energies, IMDsetup->energysendbuf))
+ {
+ imd_fatal(IMDsetup, "Error sending updated energies. Disconnecting client.\n");
+ }
+
+ if (imd_send_rvecs(IMDsetup->clientsocket, IMDsetup->nat, IMDsetup->xa, IMDsetup->coordsendbuf))
+ {
+ imd_fatal(IMDsetup, "Error sending updated positions. Disconnecting client.\n");
+ }
+ }
+#else
+ gmx_incons("IMD_send_positions called without IMD support.");
+#endif
+}
+
+
+extern void IMD_prep_energies_send_positions(gmx_bool bIMD, gmx_bool bIMDstep,
+ t_IMD *imd, gmx_enerdata_t *enerd,
+ gmx_int64_t step, gmx_bool bHaveNewEnergies,
+ gmx_wallcycle_t wcycle)
+{
+ if (bIMD)
+ {
+#ifdef GMX_IMD
+ wallcycle_start(wcycle, ewcIMD);
+
+ /* Update time step for IMD and prepare IMD energy record if we have new energies. */
+ IMD_fill_energy_record(TRUE, imd, enerd, step, bHaveNewEnergies);
+
+ if (bIMDstep)
+ {
+ /* Send positions and energies to VMD client via IMD */
+ IMD_send_positions(imd);
+ }
+
+ wallcycle_stop(wcycle, ewcIMD);
+#else
+ gmx_incons("IMD_prep_energies_send_positions called without IMD support.");
+#endif
+ }
+}
+
+
+extern int IMD_get_step(t_gmx_IMD IMDsetup)
+{
+ return IMDsetup->nstimd;
+}
+
+
+extern void IMD_apply_forces(gmx_bool bIMD, t_IMD *imd, t_commrec *cr, rvec *f,
+ gmx_wallcycle_t wcycle)
+{
+ int i, j;
+ int locndx;
+ t_gmx_IMD_setup *IMDsetup;
+
+
+ if (bIMD)
+ {
+#ifdef GMX_IMD
+ wallcycle_start(wcycle, ewcIMD);
+
+ IMDsetup = imd->setup;
+
+ /* Are forces allowed at all? If not we're done */
+ if (!IMDsetup->bForceActivated)
+ {
+ return;
+ }
+
+ for (i = 0; i < IMDsetup->nforces; i++)
+ {
+ /* j are the indices in the "System group".*/
+ j = IMDsetup->ind[IMDsetup->f_ind[i]];
+
+ /* check if this is a local atom and find out locndx */
+ if (PAR(cr) && ga2la_get_home(cr->dd->ga2la, j, &locndx))
+ {
+ j = locndx;
+ }
+
+ rvec_inc(f[j], IMDsetup->f[i]);
+ }
+
+ wallcycle_start(wcycle, ewcIMD);
+#else
+ gmx_incons("IMD_apply_forces called without IMD support.");
+#endif
+ }
+}
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \defgroup module_imd Interactive molecular dynamics (IMD)
+ * \ingroup group_mdrun
+ *
+ * \brief
+ * Allows mdrun to interface with VMD via the interactive molecular dynamics
+ * (IMD) protocol.
+ *
+ * \author Martin Hoefling, Carsten Kutzner <ckutzne@gwdg.de>
+ *
+ */
+
+/*! \libinternal \file
+ *
+ * \brief
+ * This file contains datatypes and function declarations necessary for mdrun
+ * to interface with VMD via the interactive molecular dynamics protocol.
+ *
+ * \author Martin Hoefling, Carsten Kutzner <ckutzne@gwdg.de>
+ *
+ * \inlibraryapi
+ * \ingroup module_imd
+ */
+
+#ifndef GMX_IMD_IMD_H
+#define GMX_IMD_IMD_H
+
+#include "typedefs.h"
+#include "../fileio/filenm.h"
+
+#ifdef GMX_NATIVE_WINDOWS
+#include <Windows.h>
+#define NOFLAGS 0
+#endif
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+static const char IMDstr[] = "IMD:"; /**< Tag output from the IMD module with this string. */
+
+
+/*! \brief Writes out the group of atoms selected for interactive manipulation.
+ *
+ * Called by grompp.
+ * The resulting file has to be read in by VMD if one wants it to connect to mdrun.
+ *
+ * \param bIMD Only springs into action if bIMD is TRUE. Otherwise returns directly.
+ * \param ir Structure containing MD input parameters, among those
+ * the IMD data structure.
+ * \param state The current state of the MD system.
+ * \param sys The global, complete system topology.
+ * \param nfile Number of files.
+ * \param fnm Filename struct.
+ */
+extern void write_IMDgroup_to_file(gmx_bool bIMD, t_inputrec *ir, t_state *state,
+ gmx_mtop_t *sys, int nfile, const t_filenm fnm[]);
+
+
+/*! \brief Make a selection of the home atoms for the IMD group.
+ *
+ * Should be called at every domain decomposition.
+ * Each node checks which of the atoms from "ind" are local and puts its local
+ * atom numbers into the "ind_local" array. Furthermore, in "xa_ind" it is
+ * stored at which position each local atom belongs in the assembled/collective
+ * array, so that on the master node all positions can be merged into the
+ * assembled array correctly.
+ *
+ * \param bIMD Only springs into action if bIMD is TRUE. Otherwise returns directly.
+ * \param dd Structure containing domain decomposition data.
+ * \param imd The IMD group of atoms.
+ */
+extern void dd_make_local_IMD_atoms(gmx_bool bIMD, gmx_domdec_t *dd, t_IMD *imd);
+
+
+/*! \brief Initializes (or disables) IMD.
+ *
+ * This function is called before the main MD loop over time steps,
+ * and it must be called prior to any call to dd_partition_system if in parallel.
+ *
+ * \param ir The inputrec structure containing the MD input parameters
+ * including a pointer to the IMD data structure.
+ * \param cr Information structure for MPI communication.
+ * \param top_global The topology of the whole system.
+ * \param fplog General output file, normally md.log.
+ * \param defnstimd Default IMD update (=communication) frequency.
+ * \param x The starting positions of the atoms.
+ * \param nfile Number of files.
+ * \param fnm Struct containing file names etc.
+ * \param oenv Output options.
+ * \param imdport Port to use for IMD connections.
+ * \param Flags Flags passed over from main, used to determine
+ * whether or not we are appending.
+ */
+extern void init_IMD(t_inputrec *ir, t_commrec *cr, gmx_mtop_t *top_global,
+ FILE *fplog, int defnstimd, rvec x[],
+ int nfile, const t_filenm fnm[], output_env_t oenv,
+ int imdport, unsigned long Flags);
+
+
+/*! \brief IMD required in this time step?
+ * Also checks for new IMD connection and syncs the nodes.
+ *
+ * \param bIMD Only springs into action if bIMD is TRUE. Otherwise returns directly.
+ * \param step The time step.
+ * \param cr Information structure for MPI communication.
+ * \param bNS Is this a neighbor searching step?
+ * \param box The simulation box.
+ * \param x The local atomic positions on this node.
+ * \param ir The inputrec structure containing the MD input parameters
+ * including a pointer to the IMD data structure.
+ * \param t The time.
+ * \param wcycle Count wallcycles of IMD routines for diagnostic output.
+ *
+ * \returns Whether or not we have to do IMD communication at this step.
+ */
+extern gmx_bool do_IMD(gmx_bool bIMD, gmx_int64_t step, t_commrec *cr, gmx_bool bNS,
+ matrix box, rvec x[], t_inputrec *ir, double t,
+ gmx_wallcycle_t wcycle);
+
+
+/*! \brief Get the IMD update frequency.
+ *
+ * \param IMDsetup Opaque pointer to IMD private data.
+ *
+ * \returns The current IMD update/communication frequency
+ */
+extern int IMD_get_step(t_gmx_IMD IMDsetup);
+
+
+/*! \brief Add external forces from a running interactive molecular dynamics session.
+ *
+ * \param bIMD Returns directly if bIMD is FALSE.
+ * \param imd The IMD data structure.
+ * \param cr Information structure for MPI communication.
+ * \param f The forces.
+ * \param wcycle Count wallcycles of IMD routines for diagnostic output.
+ */
+extern void IMD_apply_forces(gmx_bool bIMD, t_IMD *imd, t_commrec *cr, rvec *f,
+ gmx_wallcycle_t wcycle);
+
+
+/*! \brief Copy energies and convert to float from enerdata to the IMD energy record.
+ *
+ * We do no conversion, so units in client are SI!
+ *
+ * \param bIMD Only springs into action if bIMD is TRUE. Otherwise returns directly.
+ * \param imd The IMD data structure.
+ * \param enerd Contains the GROMACS energies for the different interaction types.
+ * \param step The time step.
+ * \param bHaveNewEnergies Only copy energies if we have done global summing of them before.
+ *
+ */
+extern void IMD_fill_energy_record(gmx_bool bIMD, t_IMD *imd, gmx_enerdata_t *enerd,
+ gmx_int64_t step, gmx_bool bHaveNewEnergies);
+
+
+/*! \brief Send positions and energies to the client.
+ *
+ * \param imd The IMD data structure.
+ */
+extern void IMD_send_positions(t_IMD *imd);
+
+
+/*! \brief Calls IMD_prepare_energies() and then IMD_send_positions().
+ *
+ * \param bIMD Returns directly if bIMD is FALSE.
+ * \param bIMDstep If true, transfer the positions. Otherwise just update the time step and potentially the energy record.
+ * \param imd The IMD data structure.
+ * \param enerd Contains the GROMACS energies for the different interaction types.
+ * \param step The time step.
+ * \param bHaveNewEnergies Only update the energy record if we have done global summing of the energies.
+ * \param wcycle Count wallcycles of IMD routines for diagnostic output.
+ *
+ */
+extern void IMD_prep_energies_send_positions(gmx_bool bIMD, gmx_bool bIMDstep,
+ t_IMD *imd, gmx_enerdata_t *enerd,
+ gmx_int64_t step, gmx_bool bHaveNewEnergies,
+ gmx_wallcycle_t wcycle);
+
+/*! \brief Finalize IMD and do some cleaning up.
+ *
+ * Currently, IMD finalize closes the force output file.
+ *
+ * \param bIMD Returns directly if bIMD is FALSE.
+ * \param imd The IMD data structure.
+ */
+extern void IMD_finalize(gmx_bool bIMD, t_IMD *imd);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \libinternal \file
+ *
+ * \brief
+ * Implements functions of imdsocket.h.
+ *
+ * This file re-implements vmdsock.c functions from original IMD API from scratch,
+ * see imdsocket.h for references to the IMD API.
+ *
+ * \author Martin Hoefling, Carsten Kutzner <ckutzne@gwdg.de>
+ *
+ * \inlibraryapi
+ * \ingroup module_imd
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <string.h>
+#include "smalloc.h"
+#include "gmx_fatal.h"
+#include "imdsocket.h"
+#include "imd.h"
+
+
+#ifdef GMX_NATIVE_WINDOWS
+#ifdef GMX_HAVE_WINSOCK
+/*! \brief Define socklen type on Windows. */
+typedef int socklen_t;
+
+/*! \brief Define a function to initialize winsock. */
+extern int imdsock_winsockinit()
+{
+ int ret = -1;
+
+
+ WSADATA wsd;
+
+ /* We use winsock 1.1 compatibility for now. Though I guess no one will try on Windows 95. */
+ ret = WSAStartup(MAKEWORD(1, 1), &wsd);
+ return ret;
+}
+#endif
+#else
+/* On UNIX, we can use nice errors from errno.h */
+#include <errno.h>
+#include <unistd.h>
+#endif
+
+
+/*! \brief Simple error handling. */
+#ifdef GMX_NATIVE_WINDOWS
+#define ERR_ARGS __FILE__, __LINE__, NULL
+#else
+#define ERR_ARGS __FILE__, __LINE__, strerror(errno)
+#endif
+
+
+/*! \brief Currently only 1 client connection is supported. */
+#define MAXIMDCONNECTIONS 1
+
+
+/*! \brief Print a nice error message on UNIX systems, using errno.h. */
+static void print_IMD_error(char *file, int line, char *msg)
+{
+ fprintf(stderr, "%s Error in file %s on line %d.\n", IMDstr, file, line);
+
+ if (NULL != msg)
+ {
+ fprintf(stderr, "%s\n", msg);
+ }
+}
+
+
+extern IMDSocket* imdsock_create()
+{
+ IMDSocket *sock = NULL;
+
+
+#ifdef GMX_IMD
+ snew(sock, 1);
+ /* Try to create socket: */
+ if ((sock->sockfd = socket(PF_INET, SOCK_STREAM, 0)) == -1)
+ {
+ print_IMD_error(ERR_ARGS);
+ sfree(sock);
+
+ return NULL;
+ }
+ else
+#endif
+ {
+ return sock;
+ }
+}
+
+
+extern int imdsock_bind(IMDSocket *sock, int port)
+{
+ int ret = -1;
+
+
+#ifdef GMX_IMD
+ memset(&(sock->address), 0, sizeof(sock->address));
+ sock->address.sin_family = PF_INET;
+ sock->address.sin_port = htons(port);
+
+ /* Try to bind to address and port ...*/
+ ret = bind(sock->sockfd, (struct sockaddr *) &sock->address, sizeof(sock->address));
+#endif
+
+ if (ret)
+ {
+ print_IMD_error(ERR_ARGS);
+ }
+
+ return ret;
+}
+
+
+extern int imd_sock_listen(IMDSocket *sock)
+{
+ int ret = -1;
+
+
+#ifdef GMX_IMD
+ /* Try to set to listening state */
+ ret = listen(sock->sockfd, MAXIMDCONNECTIONS);
+#endif
+
+ if (ret)
+ {
+ print_IMD_error(ERR_ARGS);
+ }
+
+ return ret;
+}
+
+
+extern IMDSocket* imdsock_accept(IMDSocket *sock)
+{
+ int ret = -1;
+
+#ifdef GMX_IMD
+ socklen_t length;
+
+
+ length = sizeof(sock->address);
+ ret = accept(sock->sockfd, (struct sockaddr *) &sock->address, &length);
+
+ /* successful, redirect to distinct clientsocket */
+ if (ret >= 0)
+ {
+ IMDSocket *newsock;
+
+ snew(newsock, 1);
+ newsock->address = sock->address;
+ newsock->sockfd = ret;
+
+ return newsock;
+ }
+ else
+#endif
+ {
+ print_IMD_error(ERR_ARGS);
+
+ return NULL;
+ }
+}
+
+
+extern int imdsock_getport(IMDSocket *sock, int *port)
+{
+ int ret = -1;
+#ifdef GMX_IMD
+ struct sockaddr_in sin;
+ socklen_t len;
+
+
+ len = sizeof(sin);
+ ret = getsockname(sock->sockfd, (struct sockaddr *) &(sock->address), &len);
+ if (ret)
+ {
+ fprintf(stderr, "%s getsockname failed with error %d.\n", IMDstr, ret);
+ print_IMD_error(ERR_ARGS);
+ }
+ else
+ {
+ *port = ntohs(sock->address.sin_port);
+ }
+#else
+ gmx_incons("imdsock_getport called without IMD support.");
+#endif
+
+ return ret;
+}
+
+
+extern int imdsock_write(IMDSocket *sock, const char *buffer, int length)
+{
+ /* No read and write on windows, we have to use send and recv instead... */
+#ifdef GMX_NATIVE_WINDOWS
+#ifdef GMX_HAVE_WINSOCK
+ return send(sock->sockfd, (const char *) buffer, length, NOFLAGS);
+#else
+ return -1;
+#endif
+#else
+ return write(sock->sockfd, buffer, length);
+#endif
+}
+
+
+extern int imdsock_read(IMDSocket *sock, char *buffer, int length)
+{
+ /* See above... */
+#ifdef GMX_NATIVE_WINDOWS
+#ifdef GMX_HAVE_WINSOCK
+ return recv(sock->sockfd, (char *) buffer, length, NOFLAGS);
+#else
+ return -1;
+#endif
+#else
+ return read(sock->sockfd, buffer, length);
+#endif
+}
+
+
+extern void imdsock_shutdown(IMDSocket *sock)
+{
+ int ret = -1;
+
+
+ /* is the socket already NULL? */
+ if (sock == NULL)
+ {
+ return;
+ }
+
+#ifdef GMX_IMD
+ /* If not, try to properly shut down. */
+ ret = shutdown(sock->sockfd, 1);
+#endif
+
+ if (ret == -1)
+ {
+ fprintf(stderr, "%s Failed to shutdown socket. Did the client already disconnect?\n", IMDstr);
+ print_IMD_error(ERR_ARGS);
+ }
+}
+
+
+extern int imdsock_destroy(IMDSocket *sock)
+{
+ int ret = -1;
+
+
+ if (sock == NULL)
+ {
+ return 1;
+ }
+
+#ifdef GMX_NATIVE_WINDOWS
+ /* On Windows, this function is called closesocket */
+#ifdef GMX_HAVE_WINSOCK
+ ret = closesocket(sock->sockfd);
+#endif
+#else
+ ret = close(sock->sockfd);
+#endif
+
+ if (ret == -1)
+ {
+ sfree(sock);
+ print_IMD_error(ERR_ARGS);
+
+ return 0;
+ }
+
+ return 1;
+}
+
+
+extern int imdsock_tryread(IMDSocket *sock, int timeoutsec, int timeoutusec)
+{
+ int ret = -1;
+
+
+#ifdef GMX_IMD
+ fd_set readfds;
+ /* Create new time structure with sec and usec. */
+ struct timeval *tval;
+
+
+ snew(tval, 1);
+
+ /* clear the set */
+ FD_ZERO(&readfds);
+ /* add the socket to the read set */
+ FD_SET(sock->sockfd, &readfds);
+
+ /* set the timeout */
+ tval->tv_sec = timeoutsec;
+ tval->tv_usec = timeoutusec;
+ do
+ {
+ /* check the set for read readiness. */
+ ret = select(sock->sockfd + 1, &readfds, NULL, NULL, tval);
+ /* redo on system interrupt */
+ }
+ while (ret < 0 && errno == EINTR);
+
+ sfree(tval);
+#endif
+
+ if (ret < 0)
+ {
+ print_IMD_error(ERR_ARGS);
+ }
+
+ return ret;
+}
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \libinternal \file
+ *
+ * \brief
+ * Implements the parts of the vmdsock.h interface needed for IMD communication.
+ *
+ * \author Martin Hoefling, Carsten Kutzner <ckutzne@gwdg.de>
+ *
+ * For more information, see https://www-s.ks.uiuc.edu/Research/vmd/imd/ for general
+ * IMD information and https://www-s.ks.uiuc.edu/Research/vmd/imd/code/imdapi.tar.gz
+ * for the IMD reference implementation and API.
+ *
+ * \inlibraryapi
+ * \ingroup module_imd
+ */
+
+#ifndef GMX_IMD_IMDSOCKET_H
+#define GMX_IMD_IMDSOCKET_H
+
+/* Check if we can/should use winsock or standard UNIX sockets. */
+#ifdef GMX_NATIVE_WINDOWS
+ #ifdef GMX_HAVE_WINSOCK
+ #include <Winsock.h>
+ #define GMX_IMD
+ #endif
+#else
+#include <sys/socket.h>
+#include <netinet/in.h>
+#define GMX_IMD
+#endif
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+
+/*! \internal
+ *
+ * \brief
+ * IMD (interactive molecular dynamics) socket structure
+ *
+ */
+typedef struct
+{
+#ifdef GMX_IMD
+ struct sockaddr_in address; /**< address of socket */
+#endif
+ int sockfd; /**< socket file descriptor */
+} IMDSocket;
+
+
+
+#if defined(GMX_NATIVE_WINDOWS) && defined(GMX_HAVE_WINSOCK)
+/*! \internal
+ *
+ * \brief Function to initialize winsock
+ *
+ * \returns 0 if successful.
+ */
+extern int imdsock_winsockinit();
+#endif
+
+
+/*! \brief Create an IMD master socket.
+ *
+ * \returns The IMD socket if successful. Otherwise prints an error message and returns NULL.
+ */
+extern IMDSocket *imdsock_create();
+
+
+/*! \brief Bind the IMD socket to address and port.
+ *
+ * Prints out an error message if unsuccessful.
+ * If port == 0, bind() assigns a free port automatically.
+ *
+ *
+ * \param sock The IMD socket.
+ * \param port The port to bind to.
+ *
+ * \returns 0 if successful.
+ */
+extern int imdsock_bind(IMDSocket *sock, int port);
+
+
+/*! \brief Set socket to listening state.
+ *
+ * Prints out an error message if unsuccessful.
+ *
+ * \param sock The IMD socket.
+ *
+ * \returns 0 if successful.
+ *
+ */
+extern int imd_sock_listen(IMDSocket *sock);
+
+
+/*! \brief Accept incoming connection and redirect to client socket.
+ *
+ * Prints out an error message if unsuccessful.
+ *
+ * \param sock The IMD socket.
+ *
+ * \returns IMD socket if successful, NULL otherwise.
+ */
+extern IMDSocket *imdsock_accept(IMDSocket *sock);
+
+
+/*! \brief Get the port number used for IMD connection.
+ *
+ * Prints out an error message if unsuccessful.
+ *
+ * \param sock The IMD socket.
+ * \param port The assigned port number.
+ *
+ * \returns 0 if successful, an error code otherwise.
+ */
+extern int imdsock_getport(IMDSocket *sock, int *port);
+
+
+/*! \brief Write to socket.
+ *
+ *
+ * \param sock The IMD socket.
+ * \param buffer The data to write.
+ * \param length Number of bytes to write.
+ *
+ * \returns The number of bytes written, or -1.
+ */
+extern int imdsock_write(IMDSocket *sock, const char *buffer, int length);
+
+
+/*! \brief Read from socket.
+ *
+ * \param sock The IMD socket.
+ * \param buffer Buffer to put the read data.
+ * \param length Number of bytes to read.
+ *
+ * \returns The number of bytes read, or -1 for errors.
+ */
+extern int imdsock_read(IMDSocket *sock, char *buffer, int length);
+
+
+/*! \brief Shutdown the socket.
+ *
+ * \param sock The IMD socket.
+ *
+ */
+extern void imdsock_shutdown(IMDSocket *sock);
+
+
+/*! \brief Close the socket and free the sock struct memory.
+ *
+ * Writes an error message if unsuccessful.
+ *
+ * \param sock The IMD socket.
+ *
+ * \returns 1 on success, or 0 if unsuccessful.
+ */
+extern int imdsock_destroy(IMDSocket *sock);
+
+
+/*! \brief Try to read from the socket.
+ *
+ * Time out after waiting the interval specified.
+ * Print an error message if unsuccessful.
+ *
+ * \param sock The IMD socket.
+ * \param timeoutsec Time out seconds
+ * \param timeoutusec Time out microseconds
+ *
+ */
+extern int imdsock_tryread(IMDSocket *sock, int timeoutsec, int timeoutusec);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif
#define MD_RESETCOUNTERSHALFWAY (1<<19)
#define MD_TUNEPME (1<<20)
#define MD_TESTVERLET (1<<22)
+#define MD_IMDWAIT (1<<23)
+#define MD_IMDTERM (1<<24)
+#define MD_IMDPULL (1<<25)
/* The options for the domain decomposition MPI task ordering */
enum {
gmx_membed_t membed,
real cpt_period, real max_hours,
const char *deviceOptions,
+ int imdport,
unsigned long Flags,
gmx_walltime_accounting_t walltime_accounting);
gmx_int64_t nsteps_cmdline, int nstepout, int resetstep,
int nmultisim, int repl_ex_nst, int repl_ex_nex,
int repl_ex_seed, real pforce, real cpt_period, real max_hours,
- const char *deviceOptions, unsigned long Flags);
+ const char *deviceOptions, int imdport, unsigned long Flags);
/* Driver routine, that calls the different methods */
#ifdef __cplusplus
gmx_enfrot_t enfrot; /* Stores non-inputrec enforced rotation data */
} t_rot;
+/* Abstract type for IMD only defined in IMD.c */
+typedef struct gmx_IMD *t_gmx_IMD;
+
+typedef struct {
+ int nat; /* Number of interactive atoms */
+ atom_id *ind; /* The global indices of the interactive atoms */
+ t_gmx_IMD setup; /* Stores non-inputrec IMD data */
+} t_IMD;
/* Abstract types for position swapping only defined in swapcoords.c */
typedef struct swap *gmx_swapcoords_t;
t_rot *rot; /* The data for enforced rotation potentials */
int eSwapCoords; /* Do ion/water position exchanges (CompEL)? */
t_swapcoords *swap;
+ gmx_bool bIMD; /* Allow interactive MD sessions for this .tpr? */
+ t_IMD *imd; /* Interactive molecular dynamics */
real cos_accel; /* Acceleration for viscosity calculation */
tensor deform; /* Triclinic deformation velocities (nm/ps) */
int userint1; /* User determined parameters */
#include "gromacs/utility/qsort_threadsafe.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/pulling/pull_rotation.h"
+#include "gromacs/imd/imd.h"
#define DDRANK(dd, rank) (rank)
#define DDMASTERRANK(dd) (dd->masterrank)
dd_make_local_swap_groups(dd, ir->swap);
}
+ /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */
+ dd_make_local_IMD_atoms(ir->bIMD, dd, ir->imd);
+
add_dd_statistics(dd);
/* Make sure we only count the cycles for this DD partitioning */
#include "gromacs/linearalgebra/sparsematrix.h"
#include "gromacs/timing/wallcycle.h"
#include "gromacs/timing/walltime_accounting.h"
+#include "gromacs/imd/imd.h"
typedef struct {
t_state s;
t_graph **graph, t_mdatoms *mdatoms, gmx_global_stat_t *gstat,
gmx_vsite_t *vsite, gmx_constr_t constr,
int nfile, const t_filenm fnm[],
- gmx_mdoutf_t *outf, t_mdebin **mdebin)
+ gmx_mdoutf_t *outf, t_mdebin **mdebin,
+ int imdport, unsigned long gmx_unused Flags)
{
int i;
real dvdl_constr;
init_nrnb(nrnb);
+ /* Interactive molecular dynamics */
+ init_IMD(ir, cr, top_global, fplog, 1, state_global->x,
+ nfile, fnm, NULL, imdport, Flags);
+
if (DOMAINDECOMP(cr))
{
*top = dd_init_local_top(top_global);
em_state_t *state,
t_state *state_global, rvec *f_global)
{
- int mdof_flags;
+ int mdof_flags;
+ gmx_bool bIMDout = FALSE;
+
- if ((bX || bF || confout != NULL) && !DOMAINDECOMP(cr))
+ /* Shall we do IMD output? */
+ if (ir->bIMD)
+ {
+ bIMDout = do_per_step(step, IMD_get_step(ir->imd->setup));
+ }
+
+ if ((bX || bF || bIMDout || confout != NULL) && !DOMAINDECOMP(cr))
{
copy_em_coords(state, state_global);
f_global = state->f;
{
mdof_flags |= MDOF_F;
}
+
+ /* If we want IMD output, set appropriate MDOF flag */
+ if (ir->bIMD)
+ {
+ mdof_flags |= MDOF_IMD;
+ }
+
mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags,
top_global, step, (double)step,
&state->s, state_global, state->f, f_global);
gmx_membed_t gmx_unused membed,
real gmx_unused cpt_period, real gmx_unused max_hours,
const char gmx_unused *deviceOptions,
+ int imdport,
unsigned long gmx_unused Flags,
gmx_walltime_accounting_t walltime_accounting)
{
init_em(fplog, CG, cr, inputrec,
state_global, top_global, s_min, &top, &f, &f_global,
nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
- nfile, fnm, &outf, &mdebin);
+ nfile, fnm, &outf, &mdebin, imdport, Flags);
/* Print to log file */
print_em_start(fplog, cr, walltime_accounting, wcycle, CG);
do_log = do_per_step(step, inputrec->nstlog);
do_ene = do_per_step(step, inputrec->nstenergy);
+
+ /* Prepare IMD energy record, if bIMD is TRUE. */
+ IMD_fill_energy_record(inputrec->bIMD, inputrec->imd, enerd, step, TRUE);
+
if (do_log)
{
print_ebin_header(fplog, step, step, s_min->s.lambda[efptFEP]);
TRUE, mdebin, fcd, &(top_global->groups), &(inputrec->opts));
}
+ /* Send energies and positions to the IMD client if bIMD is TRUE. */
+ if (do_IMD(inputrec->bIMD, step, cr, TRUE, state_global->box, state_global->x, inputrec, 0, wcycle) && MASTER(cr))
+ {
+ IMD_send_positions(inputrec->imd);
+ }
+
/* Stop when the maximum force lies below tolerance.
* If we have reached machine precision, converged is already set to true.
*/
} /* End of the loop */
+ /* IMD cleanup, if bIMD is TRUE. */
+ IMD_finalize(inputrec->bIMD, inputrec->imd);
+
if (converged)
{
step--; /* we never took that last step in this case */
gmx_membed_t gmx_unused membed,
real gmx_unused cpt_period, real gmx_unused max_hours,
const char gmx_unused *deviceOptions,
+ int imdport,
unsigned long gmx_unused Flags,
gmx_walltime_accounting_t walltime_accounting)
{
init_em(fplog, LBFGS, cr, inputrec,
state, top_global, &ems, &top, &f, &f_global,
nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
- nfile, fnm, &outf, &mdebin);
+ nfile, fnm, &outf, &mdebin, imdport, Flags);
/* Do_lbfgs is not completely updated like do_steep and do_cg,
* so we free some memory again.
*/
mdof_flags |= MDOF_F;
}
+ if (inputrec->bIMD)
+ {
+ mdof_flags |= MDOF_IMD;
+ }
+
mdoutf_write_to_trajectory_files(fplog, cr, outf, mdof_flags,
top_global, step, (real)step, state, state, f, f);
TRUE, mdebin, fcd, &(top_global->groups), &(inputrec->opts));
}
+ /* Send x and E to IMD client, if bIMD is TRUE. */
+ if (do_IMD(inputrec->bIMD, step, cr, TRUE, state->box, state->x, inputrec, 0, wcycle) && MASTER(cr))
+ {
+ IMD_send_positions(inputrec->imd);
+ }
+
/* Stop when the maximum force lies below tolerance.
* If we have reached machine precision, converged is already set to true.
*/
} /* End of the loop */
+ /* IMD cleanup, if bIMD is TRUE. */
+ IMD_finalize(inputrec->bIMD, inputrec->imd);
+
if (converged)
{
step--; /* we never took that last step in this case */
gmx_membed_t gmx_unused membed,
real gmx_unused cpt_period, real gmx_unused max_hours,
const char gmx_unused *deviceOptions,
+ int imdport,
unsigned long gmx_unused Flags,
gmx_walltime_accounting_t walltime_accounting)
{
init_em(fplog, SD, cr, inputrec,
state_global, top_global, s_try, &top, &f, &f_global,
nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
- nfile, fnm, &outf, &mdebin);
+ nfile, fnm, &outf, &mdebin, imdport, Flags);
/* Print to log file */
print_em_start(fplog, cr, walltime_accounting, wcycle, SD);
upd_mdebin(mdebin, FALSE, FALSE, (double)count,
mdatoms->tmass, enerd, &s_try->s, inputrec->fepvals, inputrec->expandedvals,
s_try->s.box, NULL, NULL, vir, pres, NULL, mu_tot, constr);
+
+ /* Prepare IMD energy record, if bIMD is TRUE. */
+ IMD_fill_energy_record(inputrec->bIMD, inputrec->imd, enerd, count, TRUE);
+
print_ebin(mdoutf_get_fp_ene(outf), TRUE,
do_per_step(steps_accepted, inputrec->nstdisreout),
do_per_step(steps_accepted, inputrec->nstorireout),
bAbort = TRUE;
}
+ /* Send IMD energies and positions, if bIMD is TRUE. */
+ if (do_IMD(inputrec->bIMD, count, cr, TRUE, state_global->box, state_global->x, inputrec, 0, wcycle) && MASTER(cr))
+ {
+ IMD_send_positions(inputrec->imd);
+ }
+
count++;
} /* End of the loop */
+ /* IMD cleanup, if bIMD is TRUE. */
+ IMD_finalize(inputrec->bIMD, inputrec->imd);
+
/* Print some shit... */
if (MASTER(cr))
{
gmx_membed_t gmx_unused membed,
real gmx_unused cpt_period, real gmx_unused max_hours,
const char gmx_unused *deviceOptions,
+ int imdport,
unsigned long gmx_unused Flags,
gmx_walltime_accounting_t walltime_accounting)
{
state_global, top_global, state_work, &top,
&f, &f_global,
nrnb, mu_tot, fr, &enerd, &graph, mdatoms, &gstat, vsite, constr,
- nfile, fnm, &outf, NULL);
+ nfile, fnm, &outf, NULL, imdport, Flags);
natoms = top_global->natoms;
snew(fneg, natoms);
#include "gromacs/essentialdynamics/edsam.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/pulling/pull_rotation.h"
-
+#include "gromacs/imd/imd.h"
#include "adress.h"
#include "qmmm.h"
wallcycle_stop(wcycle, ewcROTadd);
}
+ /* Add forces from interactive molecular dynamics (IMD), if bIMD == TRUE. */
+ IMD_apply_forces(inputrec->bIMD, inputrec->imd, cr, f, wcycle);
+
if (PAR(cr) && !(cr->duty & DUTY_PME))
{
/* In case of node-splitting, the PP nodes receive the long-range
wallcycle_stop(wcycle, ewcROTadd);
}
+ /* Add forces from interactive molecular dynamics (IMD), if bIMD == TRUE. */
+ IMD_apply_forces(inputrec->bIMD, inputrec->imd, cr, f, wcycle);
+
if (PAR(cr) && !(cr->duty & DUTY_PME))
{
/* In case of node-splitting, the PP nodes receive the long-range
gmx_membed_t gmx_unused membed,
real gmx_unused cpt_period, real gmx_unused max_hours,
const char gmx_unused *deviceOptions,
+ int gmx_unused imdport,
unsigned long gmx_unused Flags,
gmx_walltime_accounting_t walltime_accounting)
{
"PME redist. X/F", "PME spread/gather", "PME 3D-FFT", "PME 3D-FFT Comm.", "PME solve LJ", "PME solve Elec",
"PME wait for PP", "Wait + Recv. PME F", "Wait GPU nonlocal", "Wait GPU local", "NB X/F buffer ops.",
"Vsite spread", "Write traj.", "Update", "Constraints", "Comm. energies",
- "Enforced rotation", "Add rot. forces", "Coordinate swapping", "Test"
+ "Enforced rotation", "Add rot. forces", "Coordinate swapping", "IMD", "Test"
};
static const char *wcsn[ewcsNR] =
ewcMOVEX, ewcGB, ewcFORCE, ewcMOVEF, ewcPMEMESH,
ewcPME_REDISTXF, ewcPME_SPREADGATHER, ewcPME_FFT, ewcPME_FFTCOMM, ewcLJPME, ewcPME_SOLVE,
ewcPMEWAITCOMM, ewcPP_PMEWAITRECVF, ewcWAIT_GPU_NB_NL, ewcWAIT_GPU_NB_L, ewcNB_XF_BUF_OPS,
- ewcVSITESPREAD, ewcTRAJ, ewcUPDATE, ewcCONSTR, ewcMoveE, ewcROT, ewcROTadd, ewcSWAP,
+ ewcVSITESPREAD, ewcTRAJ, ewcUPDATE, ewcCONSTR, ewcMoveE, ewcROT, ewcROTadd, ewcSWAP, ewcIMD,
ewcTEST, ewcNR
};
+
enum {
ewcsDD_REDIST, ewcsDD_GRID, ewcsDD_SETUPCOMM,
ewcsDD_MAKETOP, ewcsDD_MAKECONSTR, ewcsDD_TOPOTHER,
#include "gromacs/timing/walltime_accounting.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/swap/swapcoords.h"
+#include "gromacs/imd/imd.h"
#ifdef GMX_FAHCORE
#include "corewrap.h"
int repl_ex_nst, int repl_ex_nex, int repl_ex_seed, gmx_membed_t membed,
real cpt_period, real max_hours,
const char gmx_unused *deviceOptions,
+ int imdport,
unsigned long Flags,
gmx_walltime_accounting_t walltime_accounting)
{
double cycles_pmes;
gmx_bool bPMETuneTry = FALSE, bPMETuneRunning = FALSE;
+ /* Interactive MD */
+ gmx_bool bIMDstep = FALSE;
+
#ifdef GMX_FAHCORE
/* Temporary addition for FAHCORE checkpointing */
int chkpt_ret;
setup_bonded_threading(fr, &top->idef);
}
+ /* Set up interactive MD (IMD) */
+ init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, state_global->x,
+ nfile, fnm, oenv, imdport, Flags);
+
if (DOMAINDECOMP(cr))
{
/* Distribute the charge groups over the nodes from the master node */
wcycle, &nchkpt,
bCPT, bRerunMD, bLastStep, (Flags & MD_CONFOUT),
bSumEkinhOld);
+ /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
+ bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x, ir, t, wcycle);
/* kludge -- virial is lost with restart for NPT control. Must restart */
if (bStartingFromCpt && bVV)
gs.set[eglsRESETCOUNTERS] = 0;
}
+ /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
+ IMD_prep_energies_send_positions(ir->bIMD && MASTER(cr), bIMDstep, ir->imd, enerd, step, bCalcEner, wcycle);
+
}
/* End of main MD loop */
debug_gmx();
print_replica_exchange_statistics(fplog, repl_ex);
}
+ /* IMD cleanup, if bIMD is TRUE. */
+ IMD_finalize(ir->bIMD, ir->imd);
+
walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
return 0;
"is sufficient, this signal should not be sent to mpirun or",
"the [TT]mdrun[tt] process that is the parent of the others.",
"[PAR]",
+ "Interactive molecular dynamics (IMD) can be activated by using at least one",
+ "of the three IMD switches: The [TT]-imdterm[tt] switch allows to terminate the",
+ "simulation from the molecular viewer (e.g. VMD). With [TT]-imdwait[tt],",
+ "[TT]mdrun[tt] pauses whenever no IMD client is connected. Pulling from the",
+ "IMD remote can be turned on by [TT]-imdpull[tt].",
+ "The port [TT]mdrun[tt] listens to can be altered by [TT]-imdport[tt].The",
+ "file pointed to by [TT]-if[tt] contains atom indices and forces if IMD",
+ "pulling is used."
+ "[PAR]",
"When [TT]mdrun[tt] is started with MPI, it does not run niced by default."
};
t_commrec *cr;
{ efDAT, "-membed", "membed", ffOPTRD },
{ efTOP, "-mp", "membed", ffOPTRD },
{ efNDX, "-mn", "membed", ffOPTRD },
+ { efXVG, "-if", "imdforces", ffOPTWR },
{ efXVG, "-swap", "swapions", ffOPTWR }
};
#define NFILE asize(fnm)
gmx_bool bRerunVSite = FALSE;
gmx_bool bConfout = TRUE;
gmx_bool bReproducible = FALSE;
+ gmx_bool bIMDwait = FALSE;
+ gmx_bool bIMDterm = FALSE;
+ gmx_bool bIMDpull = FALSE;
int npme = -1;
int nstlist = 0;
int repl_ex_nex = 0;
int nstepout = 100;
int resetstep = -1;
- gmx_int64_t nsteps = -2; /* the value -2 means that the mdp option will be used */
+ gmx_int64_t nsteps = -2; /* the value -2 means that the mdp option will be used */
+ int imdport = 8888; /* can be almost anything, 8888 is easy to remember */
rvec realddxyz = {0, 0, 0};
const char *ddno_opt[ddnoNR+1] =
"Number of random exchanges to carry out each exchange interval (N^3 is one suggestion). -nex zero or not specified gives neighbor replica exchange." },
{ "-reseed", FALSE, etINT, {&repl_ex_seed},
"Seed for replica exchange, -1 is generate a seed" },
+ { "-imdport", FALSE, etINT, {&imdport},
+ "HIDDENIMD listening port" },
+ { "-imdwait", FALSE, etBOOL, {&bIMDwait},
+ "HIDDENPause the simulation while no IMD client is connected" },
+ { "-imdterm", FALSE, etBOOL, {&bIMDterm},
+ "HIDDENAllow termination of the simulation from IMD client" },
+ { "-imdpull", FALSE, etBOOL, {&bIMDpull},
+ "HIDDENAllow pulling in the simulation from IMD client" },
{ "-rerunvsite", FALSE, etBOOL, {&bRerunVSite},
"HIDDENRecalculate virtual site coordinates with [TT]-rerun[tt]" },
{ "-confout", FALSE, etBOOL, {&bConfout},
Flags = Flags | (bKeepAndNumCPT ? MD_KEEPANDNUMCPT : 0);
Flags = Flags | (sim_part > 1 ? MD_STARTFROMCPT : 0);
Flags = Flags | (bResetCountersHalfWay ? MD_RESETCOUNTERSHALFWAY : 0);
-
+ Flags = Flags | (bIMDwait ? MD_IMDWAIT : 0);
+ Flags = Flags | (bIMDterm ? MD_IMDTERM : 0);
+ Flags = Flags | (bIMDpull ? MD_IMDPULL : 0);
/* We postpone opening the log file if we are appending, so we can
first truncate the old log file and append to the correct position
nbpu_opt[0], nstlist,
nsteps, nstepout, resetstep,
nmultisim, repl_ex_nst, repl_ex_nex, repl_ex_seed,
- pforce, cpt_period, max_hours, deviceOptions, Flags);
+ pforce, cpt_period, max_hours, deviceOptions, imdport, Flags);
/* Log file has to be closed in mdrunner if we are appending to it
(fplog not set here) */
real cpt_period;
real max_hours;
const char *deviceOptions;
+ int imdport;
unsigned long Flags;
};
mc.nbpu_opt, mc.nstlist_cmdline,
mc.nsteps_cmdline, mc.nstepout, mc.resetstep,
mc.nmultisim, mc.repl_ex_nst, mc.repl_ex_nex, mc.repl_ex_seed, mc.pforce,
- mc.cpt_period, mc.max_hours, mc.deviceOptions, mc.Flags);
+ mc.cpt_period, mc.max_hours, mc.deviceOptions, mc.imdport, mc.Flags);
}
/* called by mdrunner() to start a specific number of threads (including
gmx_int64_t nsteps_cmdline, int nstepout, int resetstep,
int gmx_unused nmultisim, int repl_ex_nst, int repl_ex_nex,
int repl_ex_seed, real pforce, real cpt_period, real max_hours,
- const char *deviceOptions, unsigned long Flags)
+ const char *deviceOptions, int imdport, unsigned long Flags)
{
gmx_bool bForceUseGPU, bTryUseGPU;
double nodetime = 0, realtime;
membed,
cpt_period, max_hours,
deviceOptions,
+ imdport,
Flags,
walltime_accounting);
# files with code for test fixtures
moduletest.cpp
swapcoords.cpp
+ interactiveMD.cpp
# pseudo-library for code for mdrun
$<TARGET_OBJECTS:mdrun_objlib>
)
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ * \brief
+ * Tests utilities for interactive molecular dynamics (IMD) setups.
+ *
+ * \author Carsten Kutzner <ckutzne@gwdg.de>
+ * \ingroup module_mdrun
+ */
+#include "moduletest.h"
+
+namespace gmx
+{
+namespace test
+{
+
+class ImdTestFixture : public MdrunTestFixture
+{
+ protected:
+ ImdTestFixture();
+ ~ImdTestFixture();
+};
+
+
+ImdTestFixture::ImdTestFixture()
+{
+}
+
+ImdTestFixture::~ImdTestFixture()
+{
+}
+
+
+//! Test fixture for mdrun with IMD settings
+typedef gmx::test::ImdTestFixture ImdTest;
+
+/* If GROMACS was compiled with IMD support, this test checks
+ * - whether the IMD-group parameter from the .mdp file is understood,
+ * - whether mdrun understands the IMD-related command line parameters -imdpull, -imdwait, -imdterm,
+ * - whether mdrun finishes without error when IMD is enabled.
+ */
+TEST_F(ImdTest, ImdCanRun)
+{
+ std::string name = "spc2";
+ useTopGroAndNdxFromDatabase(name.c_str());
+ mdpInputFileName = fileManager_.getInputFilePath((name + "-IMD.mdp").c_str());
+
+ EXPECT_EQ(0, callGrompp());
+
+ ::gmx::test::CommandLine imdCaller;
+ imdCaller.append("mdrun");
+
+ imdCaller.addOption("-imdport", 0); // automatically assign a free port
+ imdCaller.append("-imdpull");
+ imdCaller.append("-noimdwait"); // cannot use -imdwait: then mdrun would not return control ...
+ imdCaller.append("-noimdterm");
+
+ // Do an mdrun with IMD enabled
+ ASSERT_EQ(0, callMdrun(imdCaller));
+}
+
+
+
+} // namespace test
+} // namespace gmx
--- /dev/null
+dt = 0.004
+nsteps = 2
+tcoupl = Berendsen
+tc-grps = System
+tau-t = 0.5
+ref-t = 300
+constraints = all-bonds
+cutoff-scheme = Verlet
+IMD-group = SecondWaterMolecule