2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
39 * Implements functions of imd.h.
41 * Re-implementation of basic IMD functions from NAMD/VMD from scratch,
42 * see imdsocket.h for references to the IMD API.
44 * \author Martin Hoefling, Carsten Kutzner <ckutzne@gwdg.de>
57 #if GMX_NATIVE_WINDOWS
63 #include "gromacs/commandline/filenm.h"
64 #include "gromacs/domdec/domdec_struct.h"
65 #include "gromacs/domdec/ga2la.h"
66 #include "gromacs/fileio/confio.h"
67 #include "gromacs/fileio/gmxfio.h"
68 #include "gromacs/fileio/xvgr.h"
69 #include "gromacs/gmxlib/network.h"
70 #include "gromacs/imd/imdsocket.h"
71 #include "gromacs/math/units.h"
72 #include "gromacs/math/vec.h"
73 #include "gromacs/mdlib/broadcaststructs.h"
74 #include "gromacs/mdlib/groupcoord.h"
75 #include "gromacs/mdlib/mdrun.h"
76 #include "gromacs/mdlib/sighandler.h"
77 #include "gromacs/mdlib/sim_util.h"
78 #include "gromacs/mdtypes/inputrec.h"
79 #include "gromacs/mdtypes/md_enums.h"
80 #include "gromacs/pbcutil/pbc.h"
81 #include "gromacs/timing/wallcycle.h"
82 #include "gromacs/topology/mtop_util.h"
83 #include "gromacs/topology/topology.h"
84 #include "gromacs/utility/fatalerror.h"
85 #include "gromacs/utility/smalloc.h"
87 /*! \brief How long shall we wait in seconds until we check for a connection again? */
90 /*! \brief How long shall we check for the IMD_GO? */
91 #define IMDCONNECTWAIT 2
93 /*! \brief IMD Header Size. */
95 /*! \brief IMD Protocol Version. */
101 * IMD (interactive molecular dynamics) energy record.
103 * As in the original IMD implementation. Energies in kcal/mol.
104 * NOTE: We return the energies in GROMACS / SI units,
105 * so they also show up as SI in VMD.
110 gmx_int32_t tstep; /**< time step */
111 float T_abs; /**< absolute temperature */
112 float E_tot; /**< total energy */
113 float E_pot; /**< potential energy */
114 float E_vdw; /**< van der Waals energy */
115 float E_coul; /**< Coulomb interaction energy */
116 float E_bond; /**< bonds energy */
117 float E_angle; /**< angles energy */
118 float E_dihe; /**< dihedrals energy */
119 float E_impr; /**< improper dihedrals energy */
124 * \brief IMD (interactive molecular dynamics) communication structure.
126 * This structure defines the IMD communication message header & protocol version.
130 gmx_int32_t type; /**< Type of IMD message, see IMDType_t above */
131 gmx_int32_t length; /**< Length */
136 * \brief IMD (interactive molecular dynamics) main data structure.
138 * Contains private IMD data
140 typedef struct t_gmx_IMD
142 FILE *outf; /**< Output file for IMD data, mainly forces. */
144 int nat; /**< Number of atoms that can be pulled via IMD. */
145 int nat_loc; /**< Part of the atoms that are local. */
146 int *ind; /**< Global indices of the IMD atoms. */
147 int *ind_loc; /**< Local indices of the IMD atoms. */
148 int nalloc_loc; /**< Allocation size for ind_loc. */
149 rvec *xa; /**< Positions for all IMD atoms assembled on
151 ivec *xa_shifts; /**< Shifts for all IMD atoms, to make
152 molecule(s) whole. */
153 ivec *xa_eshifts; /**< Extra shifts since last DD step. */
154 rvec *xa_old; /**< Old positions for all IMD atoms on master. */
155 int *xa_ind; /**< Position of each local atom in the
158 int nstimd; /**< Global IMD frequency, known to all nodes. */
159 int nstimd_new; /**< New frequency from IMD client, master only. */
160 int nstimd_def; /**< Default IMD frequency when disconnected. */
162 int port; /**< Port to use for network socket. */
163 IMDSocket *socket; /**< The IMD socket on the master node. */
164 IMDSocket *clientsocket; /**< The IMD socket on the client. */
165 int length; /**< Length we got with last header. */
167 gmx_bool bWConnect; /**< Shall we block and wait for connection? */
168 gmx_bool bTerminated; /**< Set if MD is terminated. */
169 gmx_bool bTerminatable; /**< Set if MD can be terminated. */
170 gmx_bool bConnected; /**< Set if connection is present. */
171 gmx_bool bNewForces; /**< Set if we received new forces. */
172 gmx_bool bForceActivated; /**< Set if pulling from VMD is allowed. */
174 IMDEnergyBlock *energies; /**< Pointer to energies we send back. */
176 gmx_int32_t vmd_nforces; /**< Number of VMD forces. */
177 gmx_int32_t *vmd_f_ind; /**< VMD forces indices. */
178 float *vmd_forces; /**< The VMD forces flat in memory. */
179 int nforces; /**< Number of actual MD forces;
180 this gets communicated to the clients. */
181 int *f_ind; /**< Force indices. */
182 rvec *f; /**< The IMD pulling forces. */
184 char *forcesendbuf; /**< Buffer for force sending. */
185 char *coordsendbuf; /**< Buffer for coordinate sending. */
186 char *energysendbuf; /**< Send buffer for energies. */
187 rvec *sendxbuf; /**< Buffer to make molecules whole before
190 t_block mols; /**< Molecules block in IMD group. */
192 /* The next block is used on the master node only to reduce the output
193 * without sacrificing information. If any of these values changes,
194 * we need to write output */
195 int old_nforces; /**< Old value for nforces. */
196 int *old_f_ind; /**< Old values for force indices. */
197 rvec *old_forces; /**< Old values for IMD pulling forces. */
203 * \brief Enum for types of IMD messages.
205 * We use the same records as the NAMD/VMD IMD implementation.
207 typedef enum IMDType_t
209 IMD_DISCONNECT, /**< client disconnect */
210 IMD_ENERGIES, /**< energy data */
211 IMD_FCOORDS, /**< atomic coordinates */
212 IMD_GO, /**< start command for the simulation */
213 IMD_HANDSHAKE, /**< handshake to determine little/big endianness */
214 IMD_KILL, /**< terminates the simulation */
215 IMD_MDCOMM, /**< force data */
216 IMD_PAUSE, /**< pauses the simulation */
217 IMD_TRATE, /**< sets the IMD transmission and processing rate */
218 IMD_IOERROR, /**< I/O error */
219 IMD_NR /**< number of entries */
224 * \brief Names of the IMDType for error messages.
226 const char *eIMDType_names[IMD_NR + 1] = {
243 /*! \brief Byte swap in case we are little-endian */
244 static gmx_int32_t gmx_htonl(gmx_int32_t src)
248 if (*(char *)&num == 1)
254 gmx_int32_t dest = 0;
256 dest |= (src & 0xFF000000) >> 24;
257 dest |= (src & 0x00FF0000) >> 8;
258 dest |= (src & 0x0000FF00) << 8;
259 dest |= (src & 0x000000FF) << 24;
265 /*! \brief Byte-unswap 32 bit word in case we are little-endian */
266 static gmx_int32_t gmx_ntohl(gmx_int32_t src)
268 return gmx_htonl(src);
271 /*! \brief Fills the header with message and the length argument. */
272 static void fill_header(IMDHeader *header, IMDMessageType type, gmx_int32_t length)
274 /* We (ab-)use htonl network function for the correct endianness */
275 header->type = gmx_htonl((gmx_int32_t) type);
276 header->length = gmx_htonl(length);
280 /*! \brief Swaps the endianess of the header. */
281 static void swap_header(IMDHeader *header)
283 /* and vice versa... */
284 header->type = gmx_ntohl(header->type);
285 header->length = gmx_ntohl(header->length);
289 /*! \brief Reads multiple bytes from socket. */
290 static gmx_int32_t imd_read_multiple(IMDSocket *socket, char *datptr, gmx_int32_t toread)
292 gmx_int32_t leftcount, countread;
296 /* Try to read while we haven't reached toread */
297 while (leftcount != 0)
299 if ((countread = imdsock_read(socket, datptr, leftcount)) < 0)
301 /* interrupted function call, try again... */
306 /* this is an unexpected error, return what we got */
309 return toread - leftcount;
312 /* nothing read, finished */
314 else if (countread == 0)
318 leftcount -= countread;
322 /* return nr of bytes read */
323 return toread - leftcount;
327 /*! \brief Writes multiple bytes to socket in analogy to imd_read_multiple. */
328 static gmx_int32_t imd_write_multiple(IMDSocket *socket, const char *datptr, gmx_int32_t towrite)
330 gmx_int32_t leftcount, countwritten;
334 while (leftcount != 0)
336 if ((countwritten = imdsock_write(socket, datptr, leftcount)) <= 0)
344 return towrite - leftcount;
347 leftcount -= countwritten;
348 datptr += countwritten;
351 return towrite - leftcount;
355 /*! \brief Handshake with IMD client. */
356 static int imd_handshake(IMDSocket *socket)
361 fill_header(&header, IMD_HANDSHAKE, 1);
362 header.length = IMDVERSION; /* client wants unswapped version */
364 return (imd_write_multiple(socket, (char *) &header, HEADERSIZE) != HEADERSIZE);
368 /*! \brief Send energies using the energy block and the send buffer. */
369 static int imd_send_energies(IMDSocket *socket, const IMDEnergyBlock *energies, char *buffer)
374 recsize = HEADERSIZE + sizeof(IMDEnergyBlock);
375 fill_header((IMDHeader *) buffer, IMD_ENERGIES, 1);
376 memcpy(buffer + HEADERSIZE, energies, sizeof(IMDEnergyBlock));
378 return (imd_write_multiple(socket, buffer, recsize) != recsize);
382 /*! \brief Receive IMD header from socket, sets the length and returns the IMD message. */
383 static IMDMessageType imd_recv_header(IMDSocket *socket, gmx_int32_t *length)
388 if (imd_read_multiple(socket, (char *) &header, HEADERSIZE) != HEADERSIZE)
392 swap_header(&header);
393 *length = header.length;
395 return (IMDMessageType) header.type;
399 /*! \brief Receive force indices and forces.
401 * The number of forces was previously communicated via the header.
403 static int imd_recv_mdcomm(IMDSocket *socket, gmx_int32_t nforces, gmx_int32_t *forcendx, float *forces)
405 int retsize, retbytes;
408 /* reading indices */
409 retsize = sizeof(gmx_int32_t) * nforces;
410 retbytes = imd_read_multiple(socket, (char *) forcendx, retsize);
411 if (retbytes != retsize)
416 /* reading forces as float array */
417 retsize = 3 * sizeof(float) * nforces;
418 retbytes = imd_read_multiple(socket, (char *) forces, retsize);
419 if (retbytes != retsize)
429 /* GROMACS specific functions for the IMD implementation */
430 void write_IMDgroup_to_file(gmx_bool bIMD, t_inputrec *ir, t_state *state,
431 gmx_mtop_t *sys, int nfile, const t_filenm fnm[])
438 IMDatoms = gmx_mtop_global_atoms(sys);
439 write_sto_conf_indexed(opt2fn("-imd", nfile, fnm), "IMDgroup", &IMDatoms,
440 as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), ir->ePBC, state->box, ir->imd->nat, ir->imd->ind);
445 void dd_make_local_IMD_atoms(gmx_bool bIMD, gmx_domdec_t *dd, t_IMD *imd)
448 t_gmx_IMD_setup *IMDsetup;
452 IMDsetup = imd->setup;
455 dd_make_local_group_indices(
456 ga2la, IMDsetup->nat, IMDsetup->ind, &IMDsetup->nat_loc,
457 &IMDsetup->ind_loc, &IMDsetup->nalloc_loc, IMDsetup->xa_ind);
463 /*! \brief Send positions from rvec.
465 * We need a separate send buffer and conversion to Angstrom.
467 static int imd_send_rvecs(IMDSocket *socket, int nat, rvec *x, char *buffer)
472 int tuplesize = 3 * sizeof(float);
475 /* Required size for the send buffer */
476 size = HEADERSIZE + 3 * sizeof(float) * nat;
479 fill_header((IMDHeader *) buffer, IMD_FCOORDS, (gmx_int32_t) nat);
480 for (i = 0; i < nat; i++)
482 sendx[0] = (float) x[i][0] * NM2A;
483 sendx[1] = (float) x[i][1] * NM2A;
484 sendx[2] = (float) x[i][2] * NM2A;
485 memcpy(buffer + HEADERSIZE + i * tuplesize, sendx, tuplesize);
488 return (imd_write_multiple(socket, buffer, size) != size);
492 /*! \brief Initializes the IMD private data. */
493 static t_gmx_IMD_setup* imd_create(int imdatoms, int nstimddef, int imdport)
495 t_gmx_IMD_setup *IMDsetup = NULL;
499 IMDsetup->nat = imdatoms;
500 IMDsetup->bTerminated = FALSE;
501 IMDsetup->bTerminatable = FALSE;
502 IMDsetup->bWConnect = FALSE;
503 IMDsetup->bConnected = FALSE;
504 IMDsetup->bForceActivated = FALSE;
505 IMDsetup->bNewForces = FALSE;
506 IMDsetup->bForceActivated = FALSE;
507 IMDsetup->nstimd = 1;
508 IMDsetup->nstimd_new = 1;
509 IMDsetup->nstimd_def = nstimddef;
513 fprintf(stderr, "%s You chose a port number < 1. Will automatically assign a free port.\n", IMDstr);
517 IMDsetup->port = imdport;
524 /*! \brief Prepare the socket on the MASTER. */
525 static void imd_prepare_master_socket(t_gmx_IMD_setup *IMDsetup)
530 #if GMX_NATIVE_WINDOWS
531 /* Winsock requires separate initialization */
532 fprintf(stderr, "%s Initializing winsock.\n", IMDstr);
533 #ifdef GMX_HAVE_WINSOCK
534 if (imdsock_winsockinit())
536 gmx_fatal(FARGS, "%s Failed to initialize winsock.\n", IMDstr);
541 /* The rest is identical, first create and bind a socket and set to listen then. */
542 fprintf(stderr, "%s Setting up incoming socket.\n", IMDstr);
543 IMDsetup->socket = imdsock_create();
544 if (!IMDsetup->socket)
546 gmx_fatal(FARGS, "%s Failed to create socket.", IMDstr);
550 ret = imdsock_bind(IMDsetup->socket, IMDsetup->port);
553 gmx_fatal(FARGS, "%s binding socket to port %d failed with error %d.\n", IMDstr, IMDsetup->port, ret);
556 if (imd_sock_listen(IMDsetup->socket))
558 gmx_fatal(FARGS, "%s socket listen failed with error %d.\n", IMDstr, ret);
561 if (imdsock_getport(IMDsetup->socket, &IMDsetup->port))
563 gmx_fatal(FARGS, "%s Could not determine port number.\n", IMDstr, ret);
566 fprintf(stderr, "%s Listening for IMD connection on port %d.\n", IMDstr, IMDsetup->port);
570 /*! \brief Disconnect the client. */
571 static void imd_disconnect(t_gmx_IMD_setup *IMDsetup)
573 /* Write out any buffered pulling data */
574 fflush(IMDsetup->outf);
576 /* we first try to shut down the clientsocket */
577 imdsock_shutdown(IMDsetup->clientsocket);
578 if (!imdsock_destroy(IMDsetup->clientsocket))
580 fprintf(stderr, "%s Failed to destroy socket.\n", IMDstr);
583 /* then we reset the IMD step to its default, and reset the connection boolean */
584 IMDsetup->nstimd_new = IMDsetup->nstimd_def;
585 IMDsetup->clientsocket = NULL;
586 IMDsetup->bConnected = FALSE;
590 /*! \brief Prints an error message and disconnects the client.
592 * Does not terminate mdrun!
594 static void imd_fatal(t_gmx_IMD_setup *IMDsetup, const char *msg)
596 fprintf(stderr, "%s %s", IMDstr, msg);
597 imd_disconnect(IMDsetup);
598 fprintf(stderr, "%s disconnected.\n", IMDstr);
602 /*! \brief Check whether we got an incoming connection. */
603 static gmx_bool imd_tryconnect(t_gmx_IMD_setup *IMDsetup)
605 if (imdsock_tryread(IMDsetup->socket, 0, 0) > 0)
607 /* yes, we got something, accept on clientsocket */
608 IMDsetup->clientsocket = imdsock_accept(IMDsetup->socket);
609 if (!IMDsetup->clientsocket)
611 fprintf(stderr, "%s Accepting the connection on the socket failed.\n", IMDstr);
615 /* handshake with client */
616 if (imd_handshake(IMDsetup->clientsocket))
618 imd_fatal(IMDsetup, "Connection failed.\n");
622 fprintf(stderr, "%s Connection established, checking if I got IMD_GO orders.\n", IMDstr);
624 /* Check if we get the proper "GO" command from client. */
625 if (imdsock_tryread(IMDsetup->clientsocket, IMDCONNECTWAIT, 0) != 1 || imd_recv_header(IMDsetup->clientsocket, &(IMDsetup->length)) != IMD_GO)
627 imd_fatal(IMDsetup, "No IMD_GO order received. IMD connection failed.\n");
631 IMDsetup->bConnected = TRUE;
640 /*! \brief Wrap imd_tryconnect in order to make it blocking.
642 * Used when the simulation should wait for an incoming connection.
644 static void imd_blockconnect(t_gmx_IMD_setup *IMDsetup)
646 /* do not wait for connection, when e.g. ctrl+c is pressed and we will terminate anyways. */
647 if (!((int) gmx_get_stop_condition() == gmx_stop_cond_none))
652 fprintf(stderr, "%s Will wait until I have a connection and IMD_GO orders.\n", IMDstr);
654 /* while we have no clientsocket... 2nd part: we should still react on ctrl+c */
655 while ((!IMDsetup->clientsocket) && ((int) gmx_get_stop_condition() == gmx_stop_cond_none))
657 imd_tryconnect(IMDsetup);
658 #if GMX_NATIVE_WINDOWS
659 /* for whatever reason, it is called Sleep on windows */
668 /*! \brief Make sure that our array holding the forces received via IMD is large enough. */
669 static void imd_prepare_vmd_Forces(t_gmx_IMD_setup *IMDsetup)
671 srenew((IMDsetup->vmd_f_ind), IMDsetup->vmd_nforces);
672 srenew((IMDsetup->vmd_forces), 3*IMDsetup->vmd_nforces);
676 /*! \brief Reads forces received via IMD. */
677 static void imd_read_vmd_Forces(t_gmx_IMD_setup *IMDsetup)
679 /* the length of the previously received header tells us the nr of forces we will receive */
680 IMDsetup->vmd_nforces = IMDsetup->length;
681 /* prepare the arrays */
682 imd_prepare_vmd_Forces(IMDsetup);
683 /* Now we read the forces... */
684 if (!(imd_recv_mdcomm(IMDsetup->clientsocket, IMDsetup->vmd_nforces, IMDsetup->vmd_f_ind, IMDsetup->vmd_forces)))
686 imd_fatal(IMDsetup, "Error while reading forces from remote. Disconnecting\n");
691 /*! \brief Prepares the MD force arrays. */
692 static void imd_prepare_MD_Forces(t_gmx_IMD_setup *IMDsetup)
694 srenew((IMDsetup->f_ind), IMDsetup->nforces);
695 srenew((IMDsetup->f ), IMDsetup->nforces);
699 /*! \brief Copy IMD forces to MD forces.
701 * Do conversion from Cal->Joule and from
702 * Angstrom -> nm and from a pointer array to arrays to 3*N array.
704 static void imd_copyto_MD_Forces(t_gmx_IMD_setup *IMDsetup)
707 real conversion = CAL2JOULE * NM2A;
710 for (i = 0; i < IMDsetup->nforces; i++)
712 /* Copy the indices, a copy is important because we may update the incoming forces
713 * whenever we receive new forces while the MD forces are only communicated upon
714 * IMD communication.*/
715 IMDsetup->f_ind[i] = IMDsetup->vmd_f_ind[i];
717 /* Convert to rvecs and do a proper unit conversion */
718 IMDsetup->f[i][0] = IMDsetup->vmd_forces[3*i ] * conversion;
719 IMDsetup->f[i][1] = IMDsetup->vmd_forces[3*i + 1] * conversion;
720 IMDsetup->f[i][2] = IMDsetup->vmd_forces[3*i + 2] * conversion;
725 /*! \brief Return TRUE if any of the forces or indices changed. */
726 static gmx_bool bForcesChanged(t_gmx_IMD_setup *IMDsetup)
731 /* First, check whether the number of pulled atoms changed */
732 if (IMDsetup->nforces != IMDsetup->old_nforces)
737 /* Second, check whether any of the involved atoms changed */
738 for (i = 0; i < IMDsetup->nforces; i++)
740 if (IMDsetup->f_ind[i] != IMDsetup->old_f_ind[i])
746 /* Third, check whether all forces are the same */
747 for (i = 0; i < IMDsetup->nforces; i++)
749 if (IMDsetup->f[i][XX] != IMDsetup->old_forces[i][XX])
753 if (IMDsetup->f[i][YY] != IMDsetup->old_forces[i][YY])
757 if (IMDsetup->f[i][ZZ] != IMDsetup->old_forces[i][ZZ])
763 /* All old and new forces are identical! */
768 /*! \brief Fill the old_f_ind and old_forces arrays with the new, old values. */
769 static void keep_old_values(t_gmx_IMD_setup *IMDsetup)
774 IMDsetup->old_nforces = IMDsetup->nforces;
776 for (i = 0; i < IMDsetup->nforces; i++)
778 IMDsetup->old_f_ind[i] = IMDsetup->f_ind[i];
779 copy_rvec(IMDsetup->f[i], IMDsetup->old_forces[i]);
784 /*! \brief Returns TRUE if any component of the two rvecs differs. */
785 static gmx_inline gmx_bool rvecs_differ(const rvec v1, const rvec v2)
790 for (i = 0; i < DIM; i++)
802 /*! \brief Write the applied pull forces to logfile.
804 * Call on master only!
806 static void output_imd_forces(t_inputrec *ir, double time)
808 t_gmx_IMD_setup *IMDsetup;
812 IMDsetup = ir->imd->setup;
814 if (bForcesChanged(IMDsetup))
816 /* Write time and total number of applied IMD forces */
817 fprintf(IMDsetup->outf, "%14.6e%6d", time, IMDsetup->nforces);
819 /* Write out the global atom indices of the pulled atoms and the forces itself,
820 * write out a force only if it has changed since the last output */
821 for (i = 0; i < IMDsetup->nforces; i++)
823 if (rvecs_differ(IMDsetup->f[i], IMDsetup->old_forces[i]))
825 fprintf(IMDsetup->outf, "%9d", IMDsetup->ind[IMDsetup->f_ind[i]] + 1);
826 fprintf(IMDsetup->outf, "%12.4e%12.4e%12.4e", IMDsetup->f[i][0], IMDsetup->f[i][1], IMDsetup->f[i][2]);
829 fprintf(IMDsetup->outf, "\n");
831 keep_old_values(IMDsetup);
836 /*! \brief Synchronize the nodes. */
837 static void imd_sync_nodes(t_inputrec *ir, t_commrec *cr, double t)
840 t_gmx_IMD_setup *IMDsetup;
843 IMDsetup = ir->imd->setup;
845 /* Notify the other nodes whether we are still connected. */
848 block_bc(cr, IMDsetup->bConnected);
851 /* ...if not connected, the job is done here. */
852 if (!IMDsetup->bConnected)
857 /* Let the other nodes know whether we got a new IMD synchronization frequency. */
860 block_bc(cr, IMDsetup->nstimd_new);
863 /* Now we all set the (new) nstimd communication time step */
864 IMDsetup->nstimd = IMDsetup->nstimd_new;
866 /* We're done if we don't allow pulling at all */
867 if (!(IMDsetup->bForceActivated))
872 /* OK, let's check if we have received forces which we need to communicate
873 * to the other nodes */
876 if (IMDsetup->bNewForces)
878 new_nforces = IMDsetup->vmd_nforces;
880 /* make the "new_forces" negative, when we did not receive new ones */
883 new_nforces = IMDsetup->vmd_nforces * -1;
887 /* make new_forces known to the clients */
890 block_bc(cr, new_nforces);
893 /* When new_natoms < 0 then we know that these are still the same forces
894 * so we don't communicate them, otherwise... */
895 if (new_nforces >= 0)
897 /* set local VMD and nforces */
898 IMDsetup->vmd_nforces = new_nforces;
899 IMDsetup->nforces = IMDsetup->vmd_nforces;
901 /* now everybody knows the number of forces in f_ind, so we can prepare
902 * the target arrays for indices and forces */
903 imd_prepare_MD_Forces(IMDsetup);
905 /* we first update the MD forces on the master by converting the VMD forces */
908 imd_copyto_MD_Forces(IMDsetup);
909 /* We also write out forces on every update, so that we know which
910 * forces are applied for every step */
913 output_imd_forces(ir, t);
917 /* In parallel mode we communicate the to-be-applied forces to the other nodes */
920 nblock_bc(cr, IMDsetup->nforces, IMDsetup->f_ind);
921 nblock_bc(cr, IMDsetup->nforces, IMDsetup->f );
924 /* done communicating the forces, reset bNewForces */
925 IMDsetup->bNewForces = FALSE;
930 /*! \brief Reads header from the client and decides what to do. */
931 static void imd_readcommand(t_gmx_IMD_setup *IMDsetup)
933 gmx_bool IMDpaused = FALSE;
934 IMDMessageType itype;
937 while (IMDsetup->clientsocket && (imdsock_tryread(IMDsetup->clientsocket, 0, 0) > 0 || IMDpaused))
939 itype = imd_recv_header(IMDsetup->clientsocket, &(IMDsetup->length));
940 /* let's see what we got: */
943 /* IMD asks us to terminate the simulation, check if the user allowed this */
945 if (IMDsetup->bTerminatable)
947 fprintf(stderr, " %s Terminating connection and running simulation (if supported by integrator).\n", IMDstr);
948 IMDsetup->bTerminated = TRUE;
949 IMDsetup->bWConnect = FALSE;
950 gmx_set_stop_condition(gmx_stop_cond_next);
954 fprintf(stderr, " %s Set -imdterm command line switch to allow mdrun termination from within IMD.\n", IMDstr);
959 /* the client doen't want to talk to us anymore */
961 fprintf(stderr, " %s Disconnecting client.\n", IMDstr);
962 imd_disconnect(IMDsetup);
965 /* we got new forces, read them and set bNewForces flag */
967 imd_read_vmd_Forces(IMDsetup);
968 IMDsetup->bNewForces = TRUE;
971 /* the client asks us to (un)pause the simulation. So we toggle the IMDpaused state */
975 fprintf(stderr, " %s Un-pause command received.\n", IMDstr);
980 fprintf(stderr, " %s Pause command received.\n", IMDstr);
986 /* the client sets a new transfer rate, if we get 0, we reset the rate
987 * to the default. VMD filters 0 however */
989 IMDsetup->nstimd_new = (IMDsetup->length > 0) ? IMDsetup->length : IMDsetup->nstimd_def;
990 fprintf(stderr, " %s Update frequency will be set to %d.\n", IMDstr, IMDsetup->nstimd_new);
993 /* Catch all rule for the remaining IMD types which we don't expect */
995 fprintf(stderr, " %s Received unexpected %s.\n", IMDstr, enum_name((int)itype, IMD_NR, eIMDType_names));
996 imd_fatal(IMDsetup, "Terminating connection\n");
1003 /*! \brief Open IMD output file and write header information.
1005 * Call on master only.
1007 static FILE *open_imd_out(const char *fn,
1008 t_gmx_IMD_setup *IMDsetup,
1010 const gmx_output_env_t *oenv,
1011 unsigned long Flags)
1016 /* Open log file of applied IMD forces if requested */
1019 /* If we append to an existing file, all the header information is already there */
1020 if (Flags & MD_APPENDFILES)
1022 fp = gmx_fio_fopen(fn, "a+");
1026 fp = gmx_fio_fopen(fn, "w+");
1027 if (IMDsetup->nat == nat_total)
1029 fprintf(fp, "# Note that you can select an IMD index group in the .mdp file if a subset of the atoms suffices.\n");
1032 xvgr_header(fp, "IMD Pull Forces", "Time (ps)", "# of Forces / Atom IDs / Forces (kJ/mol)", exvggtNONE, oenv);
1034 fprintf(fp, "# Can display and manipulate %d (of a total of %d) atoms via IMD.\n", IMDsetup->nat, nat_total);
1035 fprintf(fp, "# column 1 : time (ps)\n");
1036 fprintf(fp, "# column 2 : total number of atoms feeling an IMD pulling force at that time\n");
1037 fprintf(fp, "# cols. 3.-6 : global atom number of pulled atom, x-force, y-force, z-force (kJ/mol)\n");
1038 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");
1039 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");
1043 /* To reduce the output file size we remember the old values and output only
1044 * when something changed */
1045 snew(IMDsetup->old_f_ind, IMDsetup->nat); /* One can never pull on more atoms */
1046 snew(IMDsetup->old_forces, IMDsetup->nat);
1051 fprintf(stdout, "%s For a log of the IMD pull forces explicitly specify '-if' on the command line.\n"
1052 "%s (Not possible with energy minimization.)\n", IMDstr, IMDstr);
1059 void IMD_finalize(gmx_bool bIMD, t_IMD *imd)
1063 if (imd->setup->outf)
1065 gmx_fio_fclose(imd->setup->outf);
1072 /*! \brief Creates the molecule start-end position array of molecules in the IMD group. */
1073 static void init_imd_prepare_mols_in_imdgroup(t_gmx_IMD_setup *IMDsetup, gmx_mtop_t *top_global)
1076 int gstart, gend, count;
1077 t_block gmols, lmols;
1081 gmols = top_global->mols;
1082 nat = IMDsetup->nat;
1083 ind = IMDsetup->ind;
1087 /* check whether index is sorted */
1088 for (i = 0; i < nat-1; i++)
1090 if (ind[i] > ind[i+1])
1092 gmx_fatal(FARGS, "%s IMD index is not sorted. This is currently not supported.\n", IMDstr);
1096 snew(lmols.index, gmols.nr+1);
1099 for (i = 0; i < gmols.nr; i++)
1101 gstart = gmols.index[i];
1102 gend = gmols.index[i+1];
1104 for (ii = 0; ii < nat; ii++)
1106 if ((ind[ii] >= gstart) && (ind[ii] < gend))
1113 lmols.index[lmols.nr+1] = lmols.index[lmols.nr]+count;
1117 srenew(lmols.index, lmols.nr+1);
1118 lmols.nalloc_index = lmols.nr+1;
1119 IMDsetup->mols = lmols;
1123 /*! \brief Copied and modified from groupcoord.c shift_positions_group(). */
1124 static void shift_positions(
1126 rvec x[], /* The positions [0..nr] */
1127 ivec is, /* The shift [0..nr] */
1128 int nr) /* The number of positions */
1132 /* Loop over the group's atoms */
1135 for (i = 0; i < nr; i++)
1141 x[i][XX] = x[i][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1142 x[i][YY] = x[i][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1143 x[i][ZZ] = x[i][ZZ]-tz*box[ZZ][ZZ];
1148 for (i = 0; i < nr; i++)
1154 x[i][XX] = x[i][XX]-tx*box[XX][XX];
1155 x[i][YY] = x[i][YY]-ty*box[YY][YY];
1156 x[i][ZZ] = x[i][ZZ]-tz*box[ZZ][ZZ];
1162 /*! \brief Removes shifts of molecules diffused outside of the box. */
1163 static void imd_remove_molshifts(t_gmx_IMD_setup *IMDsetup, matrix box)
1166 ivec largest, smallest, shift;
1170 mols = IMDsetup->mols;
1172 /* for each molecule also present in IMD group */
1173 for (i = 0; i < mols.nr; i++)
1175 /* first we determine the minimum and maximum shifts for each molecule */
1177 clear_ivec(largest);
1178 clear_ivec(smallest);
1181 copy_ivec(IMDsetup->xa_shifts[mols.index[i]], largest);
1182 copy_ivec(IMDsetup->xa_shifts[mols.index[i]], smallest);
1184 for (ii = mols.index[i]+1; ii < mols.index[i+1]; ii++)
1186 if (IMDsetup->xa_shifts[ii][XX] > largest[XX])
1188 largest[XX] = IMDsetup->xa_shifts[ii][XX];
1190 if (IMDsetup->xa_shifts[ii][XX] < smallest[XX])
1192 smallest[XX] = IMDsetup->xa_shifts[ii][XX];
1195 if (IMDsetup->xa_shifts[ii][YY] > largest[YY])
1197 largest[YY] = IMDsetup->xa_shifts[ii][YY];
1199 if (IMDsetup->xa_shifts[ii][YY] < smallest[YY])
1201 smallest[YY] = IMDsetup->xa_shifts[ii][YY];
1204 if (IMDsetup->xa_shifts[ii][ZZ] > largest[ZZ])
1206 largest[ZZ] = IMDsetup->xa_shifts[ii][ZZ];
1208 if (IMDsetup->xa_shifts[ii][ZZ] < smallest[ZZ])
1210 smallest[ZZ] = IMDsetup->xa_shifts[ii][ZZ];
1215 /* check if we what we can subtract/add to the positions
1216 * to put them back into the central box */
1217 if (smallest[XX] > 0)
1219 shift[XX] = smallest[XX];
1221 if (smallest[YY] > 0)
1223 shift[YY] = smallest[YY];
1225 if (smallest[ZZ] > 0)
1227 shift[ZZ] = smallest[ZZ];
1230 if (largest[XX] < 0)
1232 shift[XX] = largest[XX];
1234 if (largest[YY] < 0)
1236 shift[YY] = largest[YY];
1238 if (largest[ZZ] < 0)
1240 shift[ZZ] = largest[ZZ];
1243 /* is there a shift at all? */
1244 if ((shift[XX]) || (shift[YY]) || (shift[ZZ]))
1246 molsize = mols.index[i+1]-mols.index[i];
1247 /* shift the positions */
1248 shift_positions(box, &(IMDsetup->xa[mols.index[i]]), shift, molsize);
1255 /*! \brief Initialize arrays used to assemble the positions from the other nodes. */
1256 static void init_imd_prepare_for_x_assembly(t_commrec *cr, rvec x[], t_gmx_IMD_setup *IMDsetup)
1261 snew(IMDsetup->xa, IMDsetup->nat);
1262 snew(IMDsetup->xa_ind, IMDsetup->nat);
1263 snew(IMDsetup->xa_shifts, IMDsetup->nat);
1264 snew(IMDsetup->xa_eshifts, IMDsetup->nat);
1265 snew(IMDsetup->xa_old, IMDsetup->nat);
1267 /* Save the original (whole) set of positions such that later the
1268 * molecule can always be made whole again */
1271 for (i = 0; i < IMDsetup->nat; i++)
1273 ii = IMDsetup->ind[i];
1274 copy_rvec(x[ii], IMDsetup->xa_old[i]);
1280 IMDsetup->nat_loc = IMDsetup->nat;
1281 IMDsetup->ind_loc = IMDsetup->ind;
1283 /* xa_ind[i] needs to be set to i for serial runs */
1284 for (i = 0; i < IMDsetup->nat; i++)
1286 IMDsetup->xa_ind[i] = i;
1290 /* Communicate initial coordinates xa_old to all processes */
1294 gmx_bcast(IMDsetup->nat * sizeof(IMDsetup->xa_old[0]), IMDsetup->xa_old, cr);
1301 /*! \brief Check for non-working integrator / parallel options. */
1302 static void imd_check_integrator_parallel(t_inputrec *ir, t_commrec *cr)
1306 if (((ir->eI) == eiSteep) || ((ir->eI) == eiCG) || ((ir->eI) == eiLBFGS) || ((ir->eI) == eiNM))
1308 gmx_fatal(FARGS, "%s Energy minimization via steep, CG, lbfgs and nm in parallel is currently not supported by IMD.\n", IMDstr);
1314 void init_IMD(t_inputrec *ir,
1316 gmx_mtop_t *top_global,
1321 const t_filenm fnm[],
1322 const gmx_output_env_t *oenv,
1324 unsigned long Flags)
1328 t_gmx_IMD_setup *IMDsetup;
1329 gmx_int32_t bufxsize;
1330 gmx_bool bIMD = FALSE;
1333 /* We will allow IMD sessions only if explicitly enabled in the .tpr file */
1334 if (FALSE == ir->bIMD)
1339 /* It seems we have a .tpr file that defines an IMD group and thus allows IMD sessions.
1340 * Check whether we can actually provide the IMD functionality for this setting: */
1343 /* Check whether IMD was enabled by one of the command line switches: */
1344 if ((Flags & MD_IMDWAIT) || (Flags & MD_IMDTERM) || (Flags & MD_IMDPULL))
1346 /* Multiple simulations or replica exchange */
1349 fprintf(stderr, "%s Cannot use IMD for multiple simulations or replica exchange.\n", IMDstr);
1351 /* OK, IMD seems to be allowed and turned on... */
1354 fprintf(stderr, "%s Enabled. This simulation will accept incoming IMD connections.\n", IMDstr);
1360 fprintf(stderr, "%s None of the -imd switches was used.\n"
1361 "%s This run will not accept incoming IMD connections\n", IMDstr, IMDstr);
1363 } /* end master only */
1365 /* Disable IMD if not all the needed functionality is there! */
1366 #if GMX_NATIVE_WINDOWS && !defined(GMX_HAVE_WINSOCK)
1368 fprintf(stderr, "Disabling IMD because the winsock library was not found at compile time.\n");
1371 /* Let the other nodes know whether we want IMD */
1376 /* ... and update our local inputrec accordingly. */
1379 /*... if not we are done.*/
1386 /* check if we're using a sane integrator / parallel combination */
1387 imd_check_integrator_parallel(ir, cr);
1391 *****************************************************
1392 * From here on we assume that IMD is turned on *
1393 *****************************************************
1397 nat_total = top_global->natoms;
1399 /* Initialize IMD setup structure. If we read in a pre-IMD .tpr file, imd->nat
1400 * will be zero. For those cases we transfer _all_ atomic positions */
1401 ir->imd->setup = imd_create(ir->imd->nat > 0 ? ir->imd->nat : nat_total,
1402 defnstimd, imdport);
1403 IMDsetup = ir->imd->setup;
1405 /* We might need to open an output file for IMD forces data */
1408 IMDsetup->outf = open_imd_out(opt2fn("-if", nfile, fnm), ir->imd->setup, nat_total, oenv, Flags);
1411 /* Make sure that we operate with a valid atom index array for the IMD atoms */
1412 if (ir->imd->nat > 0)
1414 /* Point to the user-supplied array of atom numbers */
1415 IMDsetup->ind = ir->imd->ind;
1419 /* Make a dummy (ind[i] = i) array of all atoms */
1420 snew(IMDsetup->ind, nat_total);
1421 for (i = 0; i < nat_total; i++)
1423 IMDsetup->ind[i] = i;
1427 /* read environment on master and prepare socket for incoming connections */
1430 /* we allocate memory for our IMD energy structure */
1431 gmx_int32_t recsize = HEADERSIZE + sizeof(IMDEnergyBlock);
1432 snew(IMDsetup->energysendbuf, recsize);
1434 /* Shall we wait for a connection? */
1435 if (Flags & MD_IMDWAIT)
1437 IMDsetup->bWConnect = TRUE;
1438 fprintf(stderr, "%s Pausing simulation while no IMD connection present (-imdwait).\n", IMDstr);
1441 /* Will the IMD clients be able to terminate the simulation? */
1442 if (Flags & MD_IMDTERM)
1444 IMDsetup->bTerminatable = TRUE;
1445 fprintf(stderr, "%s Allow termination of the simulation from IMD client (-imdterm).\n", IMDstr);
1448 /* Is pulling from IMD client allowed? */
1449 if (Flags & MD_IMDPULL)
1451 IMDsetup->bForceActivated = TRUE;
1452 fprintf(stderr, "%s Pulling from IMD remote is enabled (-imdpull).\n", IMDstr);
1455 /* Initialize send buffers with constant size */
1456 snew(IMDsetup->sendxbuf, IMDsetup->nat);
1457 snew(IMDsetup->energies, 1);
1458 bufxsize = HEADERSIZE + 3 * sizeof(float) * IMDsetup->nat;
1459 snew(IMDsetup->coordsendbuf, bufxsize);
1462 /* do we allow interactive pulling? If so let the other nodes know. */
1465 block_bc(cr, IMDsetup->bForceActivated);
1468 /* setup the listening socket on master process */
1471 fprintf(fplog, "%s Setting port for connection requests to %d.\n", IMDstr, IMDsetup->port);
1472 fprintf(stderr, "%s Turning on IMD - port for incoming requests is %d.\n", IMDstr, IMDsetup->port);
1473 imd_prepare_master_socket(IMDsetup);
1474 /* Wait until we have a connection if specified before */
1475 if (IMDsetup->bWConnect)
1477 imd_blockconnect(IMDsetup);
1481 fprintf(stderr, "%s -imdwait not set, starting simulation.\n", IMDstr);
1484 /* Let the other nodes know whether we are connected */
1485 imd_sync_nodes(ir, cr, 0);
1487 /* Initialize arrays used to assemble the positions from the other nodes */
1488 init_imd_prepare_for_x_assembly(cr, x, IMDsetup);
1490 /* Initialize molecule blocks to make them whole later...*/
1493 init_imd_prepare_mols_in_imdgroup(IMDsetup, top_global);
1496 gmx_incons("init_IMD: this GROMACS version was not compiled with IMD support!");
1501 gmx_bool do_IMD(gmx_bool bIMD,
1509 gmx_wallcycle_t wcycle)
1511 gmx_bool imdstep = FALSE;
1512 t_gmx_IMD_setup *IMDsetup;
1522 wallcycle_start(wcycle, ewcIMD);
1524 IMDsetup = ir->imd->setup;
1526 /* read command from client and check if new incoming connection */
1529 /* If not already connected, check for new connections */
1530 if (!IMDsetup->clientsocket)
1532 if (IMDsetup->bWConnect)
1534 imd_blockconnect(IMDsetup);
1538 imd_tryconnect(IMDsetup);
1542 /* Let's see if we have new IMD messages for us */
1543 if (IMDsetup->clientsocket)
1545 imd_readcommand(IMDsetup);
1549 /* is this an IMD communication step? */
1550 imdstep = do_per_step(step, IMDsetup->nstimd);
1552 /* OK so this is an IMD step ... */
1555 /* First we sync all nodes to let everybody know whether we are connected to VMD */
1556 imd_sync_nodes(ir, cr, t);
1559 /* If a client is connected, we collect the positions
1560 * and put molecules back into the box before transfer */
1561 if ((imdstep && IMDsetup->bConnected)
1562 || bNS) /* independent of imdstep, we communicate positions at each NS step */
1564 /* Transfer the IMD positions to the master node. Every node contributes
1565 * its local positions x and stores them in the assembled xa array. */
1566 communicate_group_positions(cr, IMDsetup->xa, IMDsetup->xa_shifts, IMDsetup->xa_eshifts,
1567 TRUE, x, IMDsetup->nat, IMDsetup->nat_loc,
1568 IMDsetup->ind_loc, IMDsetup->xa_ind, IMDsetup->xa_old, box);
1570 /* If connected and master -> remove shifts */
1571 if ((imdstep && IMDsetup->bConnected) && MASTER(cr))
1573 imd_remove_molshifts(IMDsetup, box);
1577 wallcycle_stop(wcycle, ewcIMD);
1579 gmx_incons("do_IMD called without IMD support!");
1586 void IMD_fill_energy_record(gmx_bool bIMD, t_IMD *imd, gmx_enerdata_t *enerd,
1587 gmx_int64_t step, gmx_bool bHaveNewEnergies)
1589 IMDEnergyBlock *ene;
1590 t_gmx_IMD *IMDsetup;
1596 IMDsetup = imd->setup;
1598 if (IMDsetup->clientsocket)
1600 ene = IMDsetup->energies;
1604 /* In MPI-parallel simulations the energies are not accessible a at every time step.
1605 * We update them if we have new values, otherwise, the energy values from the
1606 * last global communication step are still on display in the viewer. */
1607 if (bHaveNewEnergies)
1609 ene->T_abs = (float) enerd->term[F_TEMP ];
1610 ene->E_pot = (float) enerd->term[F_EPOT ];
1611 ene->E_tot = (float) enerd->term[F_ETOT ];
1612 ene->E_bond = (float) enerd->term[F_BONDS ];
1613 ene->E_angle = (float) enerd->term[F_ANGLES ];
1614 ene->E_dihe = (float) enerd->term[F_PDIHS ];
1615 ene->E_impr = (float) enerd->term[F_IDIHS ];
1616 ene->E_vdw = (float) enerd->term[F_LJ ];
1617 ene->E_coul = (float) enerd->term[F_COUL_SR];
1621 gmx_incons("IMD_fill_energy_record called without IMD support.");
1627 void IMD_send_positions(t_IMD *imd)
1630 t_gmx_IMD *IMDsetup;
1633 IMDsetup = imd->setup;
1635 if (IMDsetup->clientsocket)
1638 if (imd_send_energies(IMDsetup->clientsocket, IMDsetup->energies, IMDsetup->energysendbuf))
1640 imd_fatal(IMDsetup, "Error sending updated energies. Disconnecting client.\n");
1643 if (imd_send_rvecs(IMDsetup->clientsocket, IMDsetup->nat, IMDsetup->xa, IMDsetup->coordsendbuf))
1645 imd_fatal(IMDsetup, "Error sending updated positions. Disconnecting client.\n");
1649 gmx_incons("IMD_send_positions called without IMD support.");
1654 void IMD_prep_energies_send_positions(gmx_bool bIMD, gmx_bool bIMDstep,
1655 t_IMD *imd, gmx_enerdata_t *enerd,
1656 gmx_int64_t step, gmx_bool bHaveNewEnergies,
1657 gmx_wallcycle_t wcycle)
1662 wallcycle_start(wcycle, ewcIMD);
1664 /* Update time step for IMD and prepare IMD energy record if we have new energies. */
1665 IMD_fill_energy_record(TRUE, imd, enerd, step, bHaveNewEnergies);
1669 /* Send positions and energies to VMD client via IMD */
1670 IMD_send_positions(imd);
1673 wallcycle_stop(wcycle, ewcIMD);
1675 gmx_incons("IMD_prep_energies_send_positions called without IMD support.");
1680 int IMD_get_step(t_gmx_IMD *IMDsetup)
1682 return IMDsetup->nstimd;
1685 void IMD_apply_forces(gmx_bool bIMD, t_IMD *imd, t_commrec *cr, rvec *f,
1686 gmx_wallcycle_t wcycle)
1690 t_gmx_IMD_setup *IMDsetup;
1696 wallcycle_start(wcycle, ewcIMD);
1698 IMDsetup = imd->setup;
1700 /* Are forces allowed at all? If not we're done */
1701 if (!IMDsetup->bForceActivated)
1706 for (i = 0; i < IMDsetup->nforces; i++)
1708 /* j are the indices in the "System group".*/
1709 j = IMDsetup->ind[IMDsetup->f_ind[i]];
1711 /* check if this is a local atom and find out locndx */
1712 if (PAR(cr) && ga2la_get_home(cr->dd->ga2la, j, &locndx))
1717 rvec_inc(f[j], IMDsetup->f[i]);
1720 wallcycle_start(wcycle, ewcIMD);
1722 gmx_incons("IMD_apply_forces called without IMD support.");