Interactive Molecular Dynamics (IMD)
authorMartin Hoefling <martin.hoefling@gmx.de>
Thu, 7 Jun 2012 11:07:21 +0000 (13:07 +0200)
committerGerrit Code Review <gerrit@gerrit.gromacs.org>
Wed, 26 Mar 2014 16:25:30 +0000 (17:25 +0100)
IMD allows to interact with and to monitor a running molecular dynamics
simulation. The protocol goes back to 2001 ("A system for interactive
molecular dynamics simulation", JE Stone, J Gullingsrud, K Schulten,
P Grayson, in: ACM symposium on interactive 3D graphics, Ed. JF Hughes
and CH Sequin, pp. 191--194, ACM SIGGRAPH). The user can watch the
running simulation (e.g. using VMD) and optionally interact with
it by pulling atoms or residues with a mouse or a force-feedback
device.
Communitcation between GROMACS and VMD is achieved via TCP sockets
and thus enables controlling an mdrun locally or one running on a
remote cluster. Every N steps, mdrun receives the applied forces from
the VMD client and sends the new positions to VMD.
Other features:
- correct PBC treatment, molecules of a (parallel) simulation are made
  whole (with respect to the configuration found in the .tpr file)
- in the .mdp file, one can define an IMD group (including the protein
  but not the water for example is useful). Only the coordinates of
  atoms belonging to this group are then transferred between mdrun and
  VMD. This can be used to reduce the performance impact to an almost
  negligible level.
- adds only two single-line function calls in the main MD loop
- and mdrun test fixture checks whether grompp and mdrun understand
  the IMD options

Change-Id: I235e07e204f2fb77f05c2f06a14b37efca5e70ea

28 files changed:
CMakeLists.txt
manual/special.tex
src/gromacs/CMakeLists.txt
src/gromacs/fileio/mdoutf.h
src/gromacs/fileio/tpxio.c
src/gromacs/gmxlib/mvdata.c
src/gromacs/gmxlib/txtdump.c
src/gromacs/gmxpreprocess/grompp.c
src/gromacs/gmxpreprocess/readir.c
src/gromacs/imd/CMakeLists.txt [new file with mode: 0644]
src/gromacs/imd/imd.c [new file with mode: 0644]
src/gromacs/imd/imd.h [new file with mode: 0644]
src/gromacs/imd/imdsocket.c [new file with mode: 0644]
src/gromacs/imd/imdsocket.h [new file with mode: 0644]
src/gromacs/legacyheaders/mdrun.h
src/gromacs/legacyheaders/types/inputrec.h
src/gromacs/mdlib/domdec.c
src/gromacs/mdlib/minimize.c
src/gromacs/mdlib/sim_util.c
src/gromacs/mdlib/tpi.c
src/gromacs/timing/wallcycle.c
src/gromacs/timing/wallcycle.h
src/programs/mdrun/md.c
src/programs/mdrun/mdrun.cpp
src/programs/mdrun/runner.c
src/programs/mdrun/tests/CMakeLists.txt
src/programs/mdrun/tests/interactiveMD.cpp [new file with mode: 0644]
src/programs/mdrun/tests/spc2-IMD.mdp [new file with mode: 0644]

index e2f4da2c1881f8990fc8d792b7a585da33c803e9..6180d173f7094e86ce6727a415024e91fc45811f 100644 (file)
@@ -326,6 +326,18 @@ if(GMX_SOFTWARE_INVSQRT)
   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()
+
 
 
 ########################################################################
index 4adb8b5f387a02150404af160fcbd30f08b0abdb..83739ed0ba59516c2b7d3e1e4cbea9df594ef63c 100644 (file)
@@ -2064,6 +2064,65 @@ found. Note that these plug-ins are in a binary format, and that format
 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
index d72b95a26fb798534df44e392269de70412d204b..6f5d146b1b1c882d2e8c3c16a483ff4c2b1618c8 100644 (file)
@@ -62,6 +62,7 @@ add_subdirectory(swap)
 add_subdirectory(essentialdynamics)
 add_subdirectory(pulling)
 add_subdirectory(simd)
+add_subdirectory(imd)
 if (NOT GMX_BUILD_MDRUN_ONLY)
     add_subdirectory(legacyheaders)
     add_subdirectory(gmxana)
index 889c03999565173e3929dfcf4cd80af27695b59f..783d3da63ab2fd4018ada8abe3e809b173113939 100644 (file)
@@ -92,5 +92,7 @@ void mdoutf_write_to_trajectory_files(FILE *fplog, t_commrec *cr,
 #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 */
index cb221d483aed7b545160c08034d68e7e18802d6f..db3ae238b648a971286132d25c2944b32b732623 100644 (file)
@@ -84,9 +84,10 @@ static const char *tpx_tag = TPX_TAG_RELEASE;
  * 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
@@ -100,7 +101,8 @@ enum tpxv {
  *
  * 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.
@@ -407,6 +409,16 @@ static void do_simtempvals(t_fileio *fio, t_simtemp *simtemp, int n_lambda, gmx_
     }
 }
 
+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. */
@@ -1557,6 +1569,25 @@ static void do_inputrec(t_fileio *fio, t_inputrec *ir, gmx_bool bRead,
         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)
index 299d2729e6e47217b22989143be355b0954116ad..7fad00743d224516026e9fc689fadadeb2a44195 100644 (file)
@@ -555,6 +555,16 @@ static void bc_adress(const t_commrec *cr, t_adress *adress)
         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;
@@ -707,6 +717,11 @@ static void bc_inputrec(const t_commrec *cr, t_inputrec *inputrec)
         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]));
index 90ff1083508a444c91bfd4863e998bf9c5f800e1..9528bc41e9ca146a23508137127befd3889fe0f9 100644 (file)
@@ -844,6 +844,13 @@ static void pr_swap(FILE *fp, int indent, t_swapcoords *swap)
 }
 
 
+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)
 {
@@ -1000,6 +1007,12 @@ void pr_inputrec(FILE *fp, int indent, const char *title, t_inputrec *ir,
             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));
index b337cd6b612d6aeb19e6fe5979d0fea7a458a6f0..4deecf51371dac8f67f7225f6a00c1d70f8f073d 100644 (file)
@@ -85,8 +85,9 @@
 #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[])
 {
@@ -1489,6 +1490,7 @@ int gmx_grompp(int argc, char *argv[])
     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  },
@@ -1502,7 +1504,9 @@ int gmx_grompp(int argc, char *argv[])
         { 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)
 
@@ -2045,6 +2049,9 @@ int gmx_grompp(int argc, char *argv[])
     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();
index 72da1996f061af4954a43daaa82be7163fde6be1..14cbadd1bc9afddf902ab571f07b6124a86087e6 100644 (file)
@@ -80,7 +80,8 @@ typedef struct t_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;
@@ -2053,6 +2054,16 @@ void get_ir(const char *mdparin, const char *mdparout,
         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");
@@ -3056,6 +3067,27 @@ static void make_swap_groups(
 }
 
 
+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,
@@ -3385,6 +3417,12 @@ void do_index(const char* mdparin, const char *ndx,
         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)
diff --git a/src/gromacs/imd/CMakeLists.txt b/src/gromacs/imd/CMakeLists.txt
new file mode 100644 (file)
index 0000000..6714e7b
--- /dev/null
@@ -0,0 +1,36 @@
+#
+# 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)
diff --git a/src/gromacs/imd/imd.c b/src/gromacs/imd/imd.c
new file mode 100644 (file)
index 0000000..dd7b3d7
--- /dev/null
@@ -0,0 +1,1700 @@
+/*
+ * 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
+    }
+}
diff --git a/src/gromacs/imd/imd.h b/src/gromacs/imd/imd.h
new file mode 100644 (file)
index 0000000..ed3e708
--- /dev/null
@@ -0,0 +1,230 @@
+/*
+ * 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
diff --git a/src/gromacs/imd/imdsocket.c b/src/gromacs/imd/imdsocket.c
new file mode 100644 (file)
index 0000000..13c602d
--- /dev/null
@@ -0,0 +1,358 @@
+/*
+ * 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;
+}
diff --git a/src/gromacs/imd/imdsocket.h b/src/gromacs/imd/imdsocket.h
new file mode 100644 (file)
index 0000000..5903db3
--- /dev/null
@@ -0,0 +1,216 @@
+/*
+ * 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
index c1a16210d5113c58f94f71fdd97d34317d24af60..e8cb420bd9034fc29675312127115cc0ead20528 100644 (file)
@@ -73,6 +73,9 @@ extern "C" {
 #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 {
@@ -96,6 +99,7 @@ typedef double gmx_integrator_t (FILE *log, t_commrec *cr,
                                  gmx_membed_t membed,
                                  real cpt_period, real max_hours,
                                  const char *deviceOptions,
+                                 int imdport,
                                  unsigned long Flags,
                                  gmx_walltime_accounting_t walltime_accounting);
 
@@ -156,7 +160,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
              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
index 54a2afd7c5a7c7254051ce7f3caaf856e3151fc3..3f463dcd53443547f885c9b8df131a04a3dcb85e 100644 (file)
@@ -271,6 +271,14 @@ typedef struct {
     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;
@@ -443,6 +451,8 @@ typedef struct {
     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                   */
index dc9132a70e3ea11064d33de8462b39de43a20331..3ec4f8120b970de18d9fc3d406f5f4937b55481f 100644 (file)
@@ -79,6 +79,7 @@
 #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)
@@ -9771,6 +9772,9 @@ void dd_partition_system(FILE                *fplog,
         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 */
index 26122e26d4c6f0ce4ed295b0ab80d246e1b43747..6af2bca9d91cecf41e3e665be800600b5f4a919b 100644 (file)
@@ -78,6 +78,7 @@
 #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;
@@ -308,7 +309,8 @@ void init_em(FILE *fplog, const char *title,
              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;
@@ -325,6 +327,10 @@ void init_em(FILE *fplog, const char *title,
 
     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);
@@ -479,9 +485,17 @@ static void write_em_traj(FILE *fplog, t_commrec *cr,
                           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;
@@ -496,6 +510,13 @@ static void write_em_traj(FILE *fplog, t_commrec *cr,
     {
         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);
@@ -939,6 +960,7 @@ double do_cg(FILE *fplog, t_commrec *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)
 {
@@ -978,7 +1000,7 @@ double do_cg(FILE *fplog, t_commrec *cr,
     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);
@@ -1445,6 +1467,10 @@ double do_cg(FILE *fplog, t_commrec *cr,
 
             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]);
@@ -1454,6 +1480,12 @@ double do_cg(FILE *fplog, t_commrec *cr,
                        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.
          */
@@ -1461,6 +1493,9 @@ double do_cg(FILE *fplog, t_commrec *cr,
 
     } /* 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 */
@@ -1553,6 +1588,7 @@ double do_lbfgs(FILE *fplog, t_commrec *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)
 {
@@ -1636,7 +1672,7 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
     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.
      */
@@ -1782,6 +1818,11 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
             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);
 
@@ -2243,6 +2284,12 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
                        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.
          */
@@ -2251,6 +2298,9 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
 
     } /* 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 */
@@ -2335,6 +2385,7 @@ double do_steep(FILE *fplog, t_commrec *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)
 {
@@ -2366,7 +2417,7 @@ double do_steep(FILE *fplog, t_commrec *cr,
     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);
@@ -2443,6 +2494,10 @@ double do_steep(FILE *fplog, t_commrec *cr,
                 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),
@@ -2512,9 +2567,18 @@ double do_steep(FILE *fplog, t_commrec *cr,
             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))
     {
@@ -2562,6 +2626,7 @@ double do_nm(FILE *fplog, t_commrec *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)
 {
@@ -2604,7 +2669,7 @@ double do_nm(FILE *fplog, t_commrec *cr,
             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);
index e2fbf7c59460d760864ec10df7801eaa5cb92217..25586cc8cf085f15bf061a1d73adb605d73ba7de 100644 (file)
@@ -89,7 +89,7 @@
 #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"
 
@@ -1538,6 +1538,9 @@ void do_force_cutsVERLET(FILE *fplog, t_commrec *cr,
         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
@@ -2009,6 +2012,9 @@ void do_force_cutsGROUP(FILE *fplog, t_commrec *cr,
         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
index 3af90cbe10269b7030349b41bbe1d175c85d3e8d..68131e1e6db387123ef565aec7b88da46e16b723 100644 (file)
@@ -125,6 +125,7 @@ double do_tpi(FILE *fplog, t_commrec *cr,
               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)
 {
index 78c0f4d2f52e1e58662738f242d64184a820f115..5b6c93161ac714a0701a53d2185fed99c7c815b6 100644 (file)
@@ -100,7 +100,7 @@ static const char *wcn[ewcNR] =
     "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] =
index f560cae5b647292d43ba19db3e40004798f916f6..45930f441267016930b9c3f66918824de94bf6ab 100644 (file)
@@ -51,10 +51,11 @@ enum {
     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,
index 2b9760ed659696c93f78c4dfb04b0e112de1df27..f6b41af0d0d88c5e5ce3000f4a0fe0dd669849e7 100644 (file)
@@ -94,6 +94,7 @@
 #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"
@@ -146,6 +147,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
              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)
 {
@@ -219,6 +221,9 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     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;
@@ -389,6 +394,10 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         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 */
@@ -1290,6 +1299,8 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                                  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)
@@ -1915,6 +1926,9 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             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();
@@ -1970,6 +1984,9 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         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;
index 4767c1a513de87e5336f1c02d609be430122395b..ca96aadcfb6bb6bf4e91d48e5ec0af0bc454013f 100644 (file)
@@ -382,6 +382,15 @@ int gmx_mdrun(int argc, char *argv[])
         "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;
@@ -419,6 +428,7 @@ int gmx_mdrun(int argc, char *argv[])
         { efDAT, "-membed", "membed",   ffOPTRD },
         { efTOP, "-mp",     "membed",   ffOPTRD },
         { efNDX, "-mn",     "membed",   ffOPTRD },
+        { efXVG, "-if",     "imdforces", ffOPTWR },
         { efXVG, "-swap",   "swapions", ffOPTWR }
     };
 #define NFILE asize(fnm)
@@ -434,6 +444,9 @@ int gmx_mdrun(int argc, char *argv[])
     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;
@@ -444,7 +457,8 @@ int gmx_mdrun(int argc, char *argv[])
     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] =
@@ -560,6 +574,14 @@ int gmx_mdrun(int argc, char *argv[])
           "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},
@@ -726,7 +748,9 @@ int gmx_mdrun(int argc, char *argv[])
     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
@@ -759,7 +783,7 @@ int gmx_mdrun(int argc, char *argv[])
                   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) */
index f3d6ed3dda654263d791ee8584063fe3ec1b3027..f40e4c8422a476781a8c911149682e91ddc079a6 100644 (file)
@@ -149,6 +149,7 @@ struct mdrunner_arglist
     real            cpt_period;
     real            max_hours;
     const char     *deviceOptions;
+    int             imdport;
     unsigned long   Flags;
 };
 
@@ -184,7 +185,7 @@ static void mdrunner_start_fn(void *arg)
              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
@@ -1081,7 +1082,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
              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;
@@ -1759,6 +1760,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
                                       membed,
                                       cpt_period, max_hours,
                                       deviceOptions,
+                                      imdport,
                                       Flags,
                                       walltime_accounting);
 
index 4da6ce4ac7ce1c88340c165a1b33fe861652c9a2..25ebd7e6068cca09883f67a6da54fd7574254500 100644 (file)
@@ -46,6 +46,7 @@ gmx_build_unit_test(
     # files with code for test fixtures
     moduletest.cpp
     swapcoords.cpp
+    interactiveMD.cpp
     # pseudo-library for code for mdrun
     $<TARGET_OBJECTS:mdrun_objlib>
     )
diff --git a/src/programs/mdrun/tests/interactiveMD.cpp b/src/programs/mdrun/tests/interactiveMD.cpp
new file mode 100644 (file)
index 0000000..bb610a7
--- /dev/null
@@ -0,0 +1,98 @@
+/*
+ * 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
diff --git a/src/programs/mdrun/tests/spc2-IMD.mdp b/src/programs/mdrun/tests/spc2-IMD.mdp
new file mode 100644 (file)
index 0000000..72824a6
--- /dev/null
@@ -0,0 +1,9 @@
+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