2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016,2017,2018,2019, 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 Interactive Molecular Dynamics
41 * Re-implementation of basic IMD functions to work with VMD,
42 * see imdsocket.h for references to the IMD API.
44 * \author Martin Hoefling, Carsten Kutzner <ckutzne@gwdg.de>
57 #include "gromacs/commandline/filenm.h"
58 #include "gromacs/domdec/domdec_struct.h"
59 #include "gromacs/domdec/ga2la.h"
60 #include "gromacs/fileio/confio.h"
61 #include "gromacs/fileio/gmxfio.h"
62 #include "gromacs/fileio/xvgr.h"
63 #include "gromacs/gmxlib/network.h"
64 #include "gromacs/imd/imdsocket.h"
65 #include "gromacs/math/units.h"
66 #include "gromacs/math/vec.h"
67 #include "gromacs/mdlib/broadcaststructs.h"
68 #include "gromacs/mdlib/groupcoord.h"
69 #include "gromacs/mdlib/sighandler.h"
70 #include "gromacs/mdlib/stat.h"
71 #include "gromacs/mdtypes/enerdata.h"
72 #include "gromacs/mdtypes/imdmodule.h"
73 #include "gromacs/mdtypes/inputrec.h"
74 #include "gromacs/mdtypes/md_enums.h"
75 #include "gromacs/mdtypes/mdrunoptions.h"
76 #include "gromacs/mdtypes/state.h"
77 #include "gromacs/pbcutil/pbc.h"
78 #include "gromacs/timing/wallcycle.h"
79 #include "gromacs/topology/mtop_util.h"
80 #include "gromacs/topology/topology.h"
81 #include "gromacs/utility/fatalerror.h"
82 #include "gromacs/utility/logger.h"
83 #include "gromacs/utility/smalloc.h"
84 #include "gromacs/utility/stringutil.h"
89 /*! \brief How long shall we wait in seconds until we check for a connection again? */
90 constexpr int c_loopWait = 1;
92 /*! \brief How long shall we check for the IMD_GO? */
93 constexpr int c_connectWait = 1;
95 /*! \brief IMD Header Size. */
96 constexpr int c_headerSize = 8;
98 /*! \brief IMD Protocol Version. */
99 constexpr int c_protocolVersion = 2;
103 * IMD (interactive molecular dynamics) energy record.
105 * As in the original IMD implementation. Energies in kcal/mol.
106 * NOTE: We return the energies in GROMACS / SI units,
107 * so they also show up as SI in VMD.
112 int32_t tstep; /**< time step */
113 float T_abs; /**< absolute temperature */
114 float E_tot; /**< total energy */
115 float E_pot; /**< potential energy */
116 float E_vdw; /**< van der Waals energy */
117 float E_coul; /**< Coulomb interaction energy */
118 float E_bond; /**< bonds energy */
119 float E_angle; /**< angles energy */
120 float E_dihe; /**< dihedrals energy */
121 float E_impr; /**< improper dihedrals energy */
126 * \brief IMD (interactive molecular dynamics) communication structure.
128 * This structure defines the IMD communication message header & protocol version.
132 int32_t type; /**< Type of IMD message, see IMDType_t above */
133 int32_t length; /**< Length */
138 * \brief Implementation type for the IMD session
140 * \todo Make this class (or one extracted from it) model
142 * \todo Use RAII for files and allocations
144 class ImdSession::Impl
148 Impl(const MDLogger &mdlog);
151 /*! \brief Prepare the socket on the MASTER. */
152 void prepareMasterSocket();
153 /*! \brief Disconnect the client. */
154 void disconnectClient();
155 /*! \brief Prints an error message and disconnects the client.
157 * Does not terminate mdrun!
159 void issueFatalError(const char *msg);
160 /*! \brief Check whether we got an incoming connection. */
162 /*! \brief Wrap imd_tryconnect in order to make it blocking.
164 * Used when the simulation should wait for an incoming connection.
167 /*! \brief Make sure that our array holding the forces received via IMD is large enough. */
168 void prepareVmdForces();
169 /*! \brief Reads forces received via IMD. */
170 void readVmdForces();
171 /*! \brief Prepares the MD force arrays. */
172 void prepareMDForces();
173 /*! \brief Copy IMD forces to MD forces.
175 * Do conversion from Cal->Joule and from
176 * Angstrom -> nm and from a pointer array to arrays to 3*N array.
178 void copyToMDForces();
179 /*! \brief Return true if any of the forces or indices changed. */
180 bool bForcesChanged() const;
181 /*! \brief Update the old_f_ind and old_forces arrays to contain the current values. */
182 void keepOldValues();
183 /*! \brief Write the applied pull forces to logfile.
185 * Call on master only!
187 void outputForces(double time);
188 /*! \brief Synchronize the nodes. */
189 void syncNodes(const t_commrec *cr, double t);
190 /*! \brief Reads header from the client and decides what to do. */
192 /*! \brief Open IMD output file and write header information.
194 * Call on master only.
196 void openOutputFile(const char *fn,
198 const gmx_output_env_t *oenv,
199 const ContinuationOptions &continuationOptions);
200 /*! \brief Creates the molecule start-end position array of molecules in the IMD group. */
201 void prepareMoleculesInImdGroup(const gmx_mtop_t *top_global);
202 /*! \brief Removes shifts of molecules diffused outside of the box. */
203 void removeMolecularShifts(const matrix box);
204 /*! \brief Initialize arrays used to assemble the positions from the other nodes. */
205 void prepareForPositionAssembly(const t_commrec *cr, const rvec x[]);
206 /*! \brief Interact with any connected VMD session */
207 bool run(int64_t step,
213 // TODO rename all the data members to have underscore suffixes
215 //! True if tpr and mdrun input combine to permit IMD sessions
216 bool sessionPossible = false;
217 //! Output file for IMD data, mainly forces.
218 FILE *outf = nullptr;
220 //! Number of atoms that can be pulled via IMD.
222 //! Part of the atoms that are local.
224 //! Global indices of the IMD atoms.
226 //! Local indices of the IMD atoms.
227 int *ind_loc = nullptr;
228 //! Allocation size for ind_loc.
230 //! Positions for all IMD atoms assembled on the master node.
232 //! Shifts for all IMD atoms, to make molecule(s) whole.
233 ivec *xa_shifts = nullptr;
234 //! Extra shifts since last DD step.
235 ivec *xa_eshifts = nullptr;
236 //! Old positions for all IMD atoms on master.
237 rvec *xa_old = nullptr;
238 //! Position of each local atom in the collective array.
239 int *xa_ind = nullptr;
241 //! Global IMD frequency, known to all ranks.
243 //! New frequency from IMD client, master only.
245 //! Default IMD frequency when disconnected.
246 int defaultNstImd = -1;
248 //! Port to use for network socket.
250 //! The IMD socket on the master node.
251 IMDSocket *socket = nullptr;
252 //! The IMD socket on the client.
253 IMDSocket *clientsocket = nullptr;
254 //! Length we got with last header.
257 //! Shall we block and wait for connection?
258 bool bWConnect = false;
259 //! Set if MD is terminated.
260 bool bTerminated = false;
261 //! Set if MD can be terminated.
262 bool bTerminatable = false;
263 //! Set if connection is present.
264 bool bConnected = false;
265 //! Set if we received new forces.
266 bool bNewForces = false;
267 //! Set if pulling from VMD is allowed.
268 bool bForceActivated = false;
270 //! Pointer to energies we send back.
271 IMDEnergyBlock *energies = nullptr;
273 //! Number of VMD forces.
274 int32_t vmd_nforces = 0;
275 //! VMD forces indices.
276 int32_t *vmd_f_ind = nullptr;
277 //! The VMD forces flat in memory.
278 float *vmd_forces = nullptr;
279 //! Number of actual MD forces; this gets communicated to the clients.
282 int *f_ind = nullptr;
283 //! The IMD pulling forces.
286 //! Buffer for force sending.
287 char *forcesendbuf = nullptr;
288 //! Buffer for coordinate sending.
289 char *coordsendbuf = nullptr;
290 //! Send buffer for energies.
291 char *energysendbuf = nullptr;
292 //! Buffer to make molecules whole before sending.
293 rvec *sendxbuf = nullptr;
295 //! Molecules block in IMD group.
298 /* The next block is used on the master node only to reduce the output
299 * without sacrificing information. If any of these values changes,
300 * we need to write output */
301 //! Old value for nforces.
303 //! Old values for force indices.
304 int *old_f_ind = nullptr;
305 //! Old values for IMD pulling forces.
306 rvec *old_forces = nullptr;
309 const MDLogger &mdlog;
310 //! Commmunication object
311 const t_commrec *cr = nullptr;
312 //! Wallcycle counting manager.
313 gmx_wallcycle *wcycle = nullptr;
314 //! Energy output handler
315 gmx_enerdata_t *enerd = nullptr;
319 * \brief Implement interactive molecular dynamics.
321 * \todo Some aspects of this module provides forces (when the user
322 * pulls on things in VMD), so in future it should have a class that
323 * models IForceProvider and is contributed to the collection of such
326 class InteractiveMolecularDynamics final : public IMDModule
329 IMdpOptionProvider *mdpOptionProvider() override { return nullptr; }
330 IMDOutputProvider *outputProvider() override { return nullptr; }
331 void initForceProviders(ForceProviders * /* forceProviders */) override {}
334 std::unique_ptr<IMDModule> createInteractiveMolecularDynamicsModule()
336 return std::make_unique<InteractiveMolecularDynamics>();
339 /*! \brief Enum for types of IMD messages.
341 * We use the same records as the NAMD/VMD IMD implementation.
343 typedef enum IMDType_t
345 IMD_DISCONNECT, /**< client disconnect */
346 IMD_ENERGIES, /**< energy data */
347 IMD_FCOORDS, /**< atomic coordinates */
348 IMD_GO, /**< start command for the simulation */
349 IMD_HANDSHAKE, /**< handshake to determine little/big endianness */
350 IMD_KILL, /**< terminates the simulation */
351 IMD_MDCOMM, /**< force data */
352 IMD_PAUSE, /**< pauses the simulation */
353 IMD_TRATE, /**< sets the IMD transmission and processing rate */
354 IMD_IOERROR, /**< I/O error */
355 IMD_NR /**< number of entries */
359 /*! \brief Names of the IMDType for error messages.
361 static const char *eIMDType_names[IMD_NR + 1] = {
377 /*! \brief Fills the header with message and the length argument. */
378 static void fill_header(IMDHeader *header, IMDMessageType type, int32_t length)
380 /* We (ab-)use htonl network function for the correct endianness */
381 header->type = imd_htonl(static_cast<int32_t>(type));
382 header->length = imd_htonl(length);
386 /*! \brief Swaps the endianess of the header. */
387 static void swap_header(IMDHeader *header)
389 /* and vice versa... */
390 header->type = imd_ntohl(header->type);
391 header->length = imd_ntohl(header->length);
395 /*! \brief Reads multiple bytes from socket. */
396 static int32_t imd_read_multiple(IMDSocket *socket, char *datptr, int32_t toread)
398 int32_t leftcount, countread;
402 /* Try to read while we haven't reached toread */
403 while (leftcount != 0)
405 if ((countread = imdsock_read(socket, datptr, leftcount)) < 0)
407 /* interrupted function call, try again... */
412 /* this is an unexpected error, return what we got */
415 return toread - leftcount;
418 /* nothing read, finished */
420 else if (countread == 0)
424 leftcount -= countread;
428 /* return nr of bytes read */
429 return toread - leftcount;
433 /*! \brief Writes multiple bytes to socket in analogy to imd_read_multiple. */
434 static int32_t imd_write_multiple(IMDSocket *socket, const char *datptr, int32_t towrite)
436 int32_t leftcount, countwritten;
440 while (leftcount != 0)
442 if ((countwritten = imdsock_write(socket, datptr, leftcount)) <= 0)
450 return towrite - leftcount;
453 leftcount -= countwritten;
454 datptr += countwritten;
457 return towrite - leftcount;
461 /*! \brief Handshake with IMD client. */
462 static int imd_handshake(IMDSocket *socket)
467 fill_header(&header, IMD_HANDSHAKE, 1);
468 header.length = c_protocolVersion; /* client wants unswapped version */
470 return static_cast<int>(imd_write_multiple(socket, reinterpret_cast<char *>(&header), c_headerSize) != c_headerSize);
474 /*! \brief Send energies using the energy block and the send buffer. */
475 static int imd_send_energies(IMDSocket *socket, const IMDEnergyBlock *energies, char *buffer)
480 recsize = c_headerSize + sizeof(IMDEnergyBlock);
481 fill_header(reinterpret_cast<IMDHeader *>(buffer), IMD_ENERGIES, 1);
482 memcpy(buffer + c_headerSize, energies, sizeof(IMDEnergyBlock));
484 return static_cast<int>(imd_write_multiple(socket, buffer, recsize) != recsize);
488 /*! \brief Receive IMD header from socket, sets the length and returns the IMD message. */
489 static IMDMessageType imd_recv_header(IMDSocket *socket, int32_t *length)
494 if (imd_read_multiple(socket, reinterpret_cast<char *>(&header), c_headerSize) != c_headerSize)
498 swap_header(&header);
499 *length = header.length;
501 return static_cast<IMDMessageType>(header.type);
505 /*! \brief Receive force indices and forces.
507 * The number of forces was previously communicated via the header.
509 static bool imd_recv_mdcomm(IMDSocket *socket, int32_t nforces, int32_t *forcendx, float *forces)
511 /* reading indices */
512 int retsize = sizeof(int32_t) * nforces;
513 int retbytes = imd_read_multiple(socket, reinterpret_cast<char *>(forcendx), retsize);
514 if (retbytes != retsize)
519 /* reading forces as float array */
520 retsize = 3 * sizeof(float) * nforces;
521 retbytes = imd_read_multiple(socket, reinterpret_cast<char *>(forces), retsize);
522 return (retbytes == retsize);
525 /* GROMACS specific functions for the IMD implementation */
526 void write_IMDgroup_to_file(bool bIMD, t_inputrec *ir, const t_state *state,
527 const gmx_mtop_t *sys, int nfile, const t_filenm fnm[])
534 IMDatoms = gmx_mtop_global_atoms(sys);
535 write_sto_conf_indexed(opt2fn("-imd", nfile, fnm), "IMDgroup", &IMDatoms,
536 state->x.rvec_array(), state->v.rvec_array(), ir->ePBC, state->box, ir->imd->nat, ir->imd->ind);
542 ImdSession::dd_make_local_IMD_atoms(const gmx_domdec_t *dd)
544 if (!impl_->sessionPossible)
549 dd_make_local_group_indices(dd->ga2la, impl_->nat, impl_->ind, &impl_->nat_loc,
550 &impl_->ind_loc, &impl_->nalloc_loc, impl_->xa_ind);
554 /*! \brief Send positions from rvec.
556 * We need a separate send buffer and conversion to Angstrom.
558 static int imd_send_rvecs(IMDSocket *socket, int nat, rvec *x, char *buffer)
563 int tuplesize = 3 * sizeof(float);
566 /* Required size for the send buffer */
567 size = c_headerSize + 3 * sizeof(float) * nat;
570 fill_header(reinterpret_cast<IMDHeader *>(buffer), IMD_FCOORDS, static_cast<int32_t>(nat));
571 for (i = 0; i < nat; i++)
573 sendx[0] = static_cast<float>(x[i][0]) * NM2A;
574 sendx[1] = static_cast<float>(x[i][1]) * NM2A;
575 sendx[2] = static_cast<float>(x[i][2]) * NM2A;
576 memcpy(buffer + c_headerSize + i * tuplesize, sendx, tuplesize);
579 return static_cast<int>(imd_write_multiple(socket, buffer, size) != size);
584 ImdSession::Impl::prepareMasterSocket()
586 if (imdsock_winsockinit() == -1)
588 gmx_fatal(FARGS, "%s Failed to initialize winsock.\n", IMDstr);
591 /* The rest is identical, first create and bind a socket and set to listen then. */
592 GMX_LOG(mdlog.warning).appendTextFormatted("%s Setting up incoming socket.", IMDstr);
593 socket = imdsock_create();
596 gmx_fatal(FARGS, "%s Failed to create socket.", IMDstr);
600 int ret = imdsock_bind(socket, port);
603 gmx_fatal(FARGS, "%s binding socket to port %d failed with error %d.\n", IMDstr, port, ret);
606 if (imd_sock_listen(socket))
608 gmx_fatal(FARGS, "%s socket listen failed with error %d.\n", IMDstr, ret);
611 if (imdsock_getport(socket, &port))
613 gmx_fatal(FARGS, "%s Could not determine port number.\n", IMDstr);
616 GMX_LOG(mdlog.warning).appendTextFormatted("%s Listening for IMD connection on port %d.", IMDstr, port);
621 ImdSession::Impl::disconnectClient()
623 /* Write out any buffered pulling data */
626 /* we first try to shut down the clientsocket */
627 imdsock_shutdown(clientsocket);
628 if (!imdsock_destroy(clientsocket))
630 GMX_LOG(mdlog.warning).appendTextFormatted("%s Failed to destroy socket.", IMDstr);
633 /* then we reset the IMD step to its default, and reset the connection boolean */
634 nstimd_new = defaultNstImd;
635 clientsocket = nullptr;
641 ImdSession::Impl::issueFatalError(const char *msg)
643 GMX_LOG(mdlog.warning).appendTextFormatted("%s %s", IMDstr, msg);
645 GMX_LOG(mdlog.warning).appendTextFormatted("%s disconnected.", IMDstr);
650 ImdSession::Impl::tryConnect()
652 if (imdsock_tryread(socket, 0, 0) > 0)
654 /* yes, we got something, accept on clientsocket */
655 clientsocket = imdsock_accept(socket);
658 GMX_LOG(mdlog.warning).appendTextFormatted("%s Accepting the connection on the socket failed.", IMDstr);
662 /* handshake with client */
663 if (imd_handshake(clientsocket))
665 issueFatalError("Connection failed.");
669 GMX_LOG(mdlog.warning).appendTextFormatted("%s Connection established, checking if I got IMD_GO orders.", IMDstr);
671 /* Check if we get the proper "GO" command from client. */
672 if (imdsock_tryread(clientsocket, c_connectWait, 0) != 1 || imd_recv_header(clientsocket, &(length)) != IMD_GO)
674 issueFatalError("No IMD_GO order received. IMD connection failed.");
688 ImdSession::Impl::blockConnect()
690 /* do not wait for connection, when e.g. ctrl+c is pressed and we will terminate anyways. */
691 if (!(static_cast<int>(gmx_get_stop_condition()) == gmx_stop_cond_none))
696 GMX_LOG(mdlog.warning).appendTextFormatted("%s Will wait until I have a connection and IMD_GO orders.", IMDstr);
698 /* while we have no clientsocket... 2nd part: we should still react on ctrl+c */
699 while ((!clientsocket) && (static_cast<int>(gmx_get_stop_condition()) == gmx_stop_cond_none))
702 imd_sleep(c_loopWait);
708 ImdSession::Impl::prepareVmdForces()
710 srenew((vmd_f_ind), vmd_nforces);
711 srenew((vmd_forces), 3*vmd_nforces);
716 ImdSession::Impl::readVmdForces()
718 /* the length of the previously received header tells us the nr of forces we will receive */
719 vmd_nforces = length;
720 /* prepare the arrays */
722 /* Now we read the forces... */
723 if (!(imd_recv_mdcomm(clientsocket, vmd_nforces, vmd_f_ind, vmd_forces)))
725 issueFatalError("Error while reading forces from remote. Disconnecting");
731 ImdSession::Impl::prepareMDForces()
733 srenew((f_ind), nforces);
734 srenew((f ), nforces);
739 ImdSession::Impl::copyToMDForces()
742 real conversion = CAL2JOULE * NM2A;
745 for (i = 0; i < nforces; i++)
747 /* Copy the indices, a copy is important because we may update the incoming forces
748 * whenever we receive new forces while the MD forces are only communicated upon
749 * IMD communication.*/
750 f_ind[i] = vmd_f_ind[i];
752 /* Convert to rvecs and do a proper unit conversion */
753 f[i][0] = vmd_forces[3*i ] * conversion;
754 f[i][1] = vmd_forces[3*i + 1] * conversion;
755 f[i][2] = vmd_forces[3*i + 2] * conversion;
760 /*! \brief Returns true if any component of the two rvecs differs. */
761 static inline bool rvecs_differ(const rvec v1, const rvec v2)
763 for (int i = 0; i < DIM; i++)
775 ImdSession::Impl::bForcesChanged() const
777 /* First, check whether the number of pulled atoms changed */
778 if (nforces != old_nforces)
783 /* Second, check whether any of the involved atoms changed */
784 for (int i = 0; i < nforces; i++)
786 if (f_ind[i] != old_f_ind[i])
792 /* Third, check whether all forces are the same */
793 for (int i = 0; i < nforces; i++)
795 if (rvecs_differ(f[i], old_forces[i]))
801 /* All old and new forces are identical! */
807 ImdSession::Impl::keepOldValues()
809 old_nforces = nforces;
811 for (int i = 0; i < nforces; i++)
813 old_f_ind[i] = f_ind[i];
814 copy_rvec(f[i], old_forces[i]);
820 ImdSession::Impl::outputForces(double time)
822 if (!bForcesChanged())
827 /* Write time and total number of applied IMD forces */
828 fprintf(outf, "%14.6e%6d", time, nforces);
830 /* Write out the global atom indices of the pulled atoms and the forces itself,
831 * write out a force only if it has changed since the last output */
832 for (int i = 0; i < nforces; i++)
834 if (rvecs_differ(f[i], old_forces[i]))
836 fprintf(outf, "%9d", ind[f_ind[i]] + 1);
837 fprintf(outf, "%12.4e%12.4e%12.4e", f[i][0], f[i][1], f[i][2]);
847 ImdSession::Impl::syncNodes(const t_commrec *cr, double t)
849 /* Notify the other nodes whether we are still connected. */
852 block_bc(cr, bConnected);
855 /* ...if not connected, the job is done here. */
861 /* Let the other nodes know whether we got a new IMD synchronization frequency. */
864 block_bc(cr, nstimd_new);
867 /* Now we all set the (new) nstimd communication time step */
870 /* We're done if we don't allow pulling at all */
871 if (!(bForceActivated))
876 /* OK, let's check if we have received forces which we need to communicate
877 * to the other nodes */
883 new_nforces = vmd_nforces;
885 /* make the "new_forces" negative, when we did not receive new ones */
888 new_nforces = vmd_nforces * -1;
892 /* make new_forces known to the clients */
895 block_bc(cr, new_nforces);
898 /* When new_natoms < 0 then we know that these are still the same forces
899 * so we don't communicate them, otherwise... */
905 /* set local VMD and nforces */
906 vmd_nforces = new_nforces;
907 nforces = vmd_nforces;
909 /* now everybody knows the number of forces in f_ind, so we can prepare
910 * the target arrays for indices and forces */
913 /* we first update the MD forces on the master by converting the VMD forces */
917 /* We also write out forces on every update, so that we know which
918 * forces are applied for every step */
925 /* In parallel mode we communicate the to-be-applied forces to the other nodes */
928 nblock_bc(cr, nforces, f_ind);
929 nblock_bc(cr, nforces, f );
932 /* done communicating the forces, reset bNewForces */
938 ImdSession::Impl::readCommand()
940 bool IMDpaused = false;
942 while (clientsocket && (imdsock_tryread(clientsocket, 0, 0) > 0 || IMDpaused))
944 IMDMessageType itype = imd_recv_header(clientsocket, &(length));
945 /* let's see what we got: */
948 /* IMD asks us to terminate the simulation, check if the user allowed this */
952 GMX_LOG(mdlog.warning).appendTextFormatted(" %s Terminating connection and running simulation (if supported by integrator).", IMDstr);
955 gmx_set_stop_condition(gmx_stop_cond_next);
959 GMX_LOG(mdlog.warning).appendTextFormatted(" %s Set -imdterm command line switch to allow mdrun termination from within IMD.", IMDstr);
964 /* the client doen't want to talk to us anymore */
966 GMX_LOG(mdlog.warning).appendTextFormatted(" %s Disconnecting client.", IMDstr);
970 /* we got new forces, read them and set bNewForces flag */
976 /* the client asks us to (un)pause the simulation. So we toggle the IMDpaused state */
980 GMX_LOG(mdlog.warning).appendTextFormatted(" %s Un-pause command received.", IMDstr);
985 GMX_LOG(mdlog.warning).appendTextFormatted(" %s Pause command received.", IMDstr);
991 /* the client sets a new transfer rate, if we get 0, we reset the rate
992 * to the default. VMD filters 0 however */
994 nstimd_new = (length > 0) ? length : defaultNstImd;
995 GMX_LOG(mdlog.warning).appendTextFormatted(" %s Update frequency will be set to %d.", IMDstr, nstimd_new);
998 /* Catch all rule for the remaining IMD types which we don't expect */
1000 GMX_LOG(mdlog.warning).appendTextFormatted(" %s Received unexpected %s.", IMDstr, enum_name(static_cast<int>(itype), IMD_NR, eIMDType_names));
1001 issueFatalError("Terminating connection");
1009 ImdSession::Impl::openOutputFile(const char *fn,
1011 const gmx_output_env_t *oenv,
1012 const ContinuationOptions &continuationOptions)
1014 /* Open log file of applied IMD forces if requested */
1017 fprintf(stdout, "%s For a log of the IMD pull forces explicitly specify '-if' on the command line.\n"
1018 "%s (Not possible with energy minimization.)\n", IMDstr, IMDstr);
1022 /* If we append to an existing file, all the header information is already there */
1023 if (continuationOptions.appendFiles)
1025 outf = gmx_fio_fopen(fn, "a+");
1029 outf = gmx_fio_fopen(fn, "w+");
1030 if (nat == nat_total)
1032 fprintf(outf, "# Note that you can select an IMD index group in the .mdp file if a subset of the atoms suffices.\n");
1035 xvgr_header(outf, "IMD Pull Forces", "Time (ps)", "# of Forces / Atom IDs / Forces (kJ/mol)", exvggtNONE, oenv);
1037 fprintf(outf, "# Can display and manipulate %d (of a total of %d) atoms via IMD.\n", nat, nat_total);
1038 fprintf(outf, "# column 1 : time (ps)\n");
1039 fprintf(outf, "# column 2 : total number of atoms feeling an IMD pulling force at that time\n");
1040 fprintf(outf, "# cols. 3.-6 : global atom number of pulled atom, x-force, y-force, z-force (kJ/mol)\n");
1041 fprintf(outf, "# then follow : atom-ID, f[x], f[y], f[z] for more atoms in case the force on multiple atoms is changed simultaneously.\n");
1042 fprintf(outf, "# Note that the force on any atom is always equal to the last value for that atom-ID found in the data.\n");
1046 /* To reduce the output file size we remember the old values and output only
1047 * when something changed */
1048 snew(old_f_ind, nat); /* One can never pull on more atoms */
1049 snew(old_forces, nat);
1053 ImdSession::Impl::Impl(const MDLogger &mdlog)
1059 ImdSession::Impl::~Impl()
1063 gmx_fio_fclose(outf);
1070 ImdSession::Impl::prepareMoleculesInImdGroup(const gmx_mtop_t *top_global)
1072 /* check whether index is sorted */
1073 for (int i = 0; i < nat-1; i++)
1075 if (ind[i] > ind[i+1])
1077 gmx_fatal(FARGS, "%s IMD index is not sorted. This is currently not supported.\n", IMDstr);
1081 RangePartitioning gmols = gmx_mtop_molecules(*top_global);
1084 snew(lmols.index, gmols.numBlocks() + 1);
1087 for (int i = 0; i < gmols.numBlocks(); i++)
1089 auto mol = gmols.block(i);
1091 for (int ii = 0; ii < nat; ii++)
1093 if (mol.inRange(ind[ii]))
1100 lmols.index[lmols.nr+1] = lmols.index[lmols.nr]+count;
1104 srenew(lmols.index, lmols.nr+1);
1105 lmols.nalloc_index = lmols.nr+1;
1110 /*! \brief Copied and modified from groupcoord.c shift_positions_group(). */
1111 static void shift_positions(
1113 rvec x[], /* The positions [0..nr] */
1114 const ivec is, /* The shift [0..nr] */
1115 int nr) /* The number of positions */
1119 /* Loop over the group's atoms */
1122 for (i = 0; i < nr; i++)
1128 x[i][XX] = x[i][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
1129 x[i][YY] = x[i][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
1130 x[i][ZZ] = x[i][ZZ]-tz*box[ZZ][ZZ];
1135 for (i = 0; i < nr; i++)
1141 x[i][XX] = x[i][XX]-tx*box[XX][XX];
1142 x[i][YY] = x[i][YY]-ty*box[YY][YY];
1143 x[i][ZZ] = x[i][ZZ]-tz*box[ZZ][ZZ];
1149 void ImdSession::Impl::removeMolecularShifts(const matrix box)
1151 /* for each molecule also present in IMD group */
1152 for (int i = 0; i < mols.nr; i++)
1154 /* first we determine the minimum and maximum shifts for each molecule */
1156 ivec largest, smallest, shift;
1157 clear_ivec(largest);
1158 clear_ivec(smallest);
1161 copy_ivec(xa_shifts[mols.index[i]], largest);
1162 copy_ivec(xa_shifts[mols.index[i]], smallest);
1164 for (int ii = mols.index[i]+1; ii < mols.index[i+1]; ii++)
1166 if (xa_shifts[ii][XX] > largest[XX])
1168 largest[XX] = xa_shifts[ii][XX];
1170 if (xa_shifts[ii][XX] < smallest[XX])
1172 smallest[XX] = xa_shifts[ii][XX];
1175 if (xa_shifts[ii][YY] > largest[YY])
1177 largest[YY] = xa_shifts[ii][YY];
1179 if (xa_shifts[ii][YY] < smallest[YY])
1181 smallest[YY] = xa_shifts[ii][YY];
1184 if (xa_shifts[ii][ZZ] > largest[ZZ])
1186 largest[ZZ] = xa_shifts[ii][ZZ];
1188 if (xa_shifts[ii][ZZ] < smallest[ZZ])
1190 smallest[ZZ] = xa_shifts[ii][ZZ];
1195 /* check if we what we can subtract/add to the positions
1196 * to put them back into the central box */
1197 if (smallest[XX] > 0)
1199 shift[XX] = smallest[XX];
1201 if (smallest[YY] > 0)
1203 shift[YY] = smallest[YY];
1205 if (smallest[ZZ] > 0)
1207 shift[ZZ] = smallest[ZZ];
1210 if (largest[XX] < 0)
1212 shift[XX] = largest[XX];
1214 if (largest[YY] < 0)
1216 shift[YY] = largest[YY];
1218 if (largest[ZZ] < 0)
1220 shift[ZZ] = largest[ZZ];
1223 /* is there a shift at all? */
1224 if ((shift[XX]) || (shift[YY]) || (shift[ZZ]))
1226 int molsize = mols.index[i+1]-mols.index[i];
1227 /* shift the positions */
1228 shift_positions(box, &(xa[mols.index[i]]), shift, molsize);
1236 ImdSession::Impl::prepareForPositionAssembly(const t_commrec *cr, const rvec x[])
1240 snew(xa_shifts, nat);
1241 snew(xa_eshifts, nat);
1244 /* Save the original (whole) set of positions such that later the
1245 * molecule can always be made whole again */
1248 for (int i = 0; i < nat; i++)
1251 copy_rvec(x[ii], xa_old[i]);
1260 /* xa_ind[i] needs to be set to i for serial runs */
1261 for (int i = 0; i < nat; i++)
1267 /* Communicate initial coordinates xa_old to all processes */
1270 gmx_bcast(nat * sizeof(xa_old[0]), xa_old, cr);
1275 /*! \brief Check for non-working integrator / parallel options. */
1276 static void imd_check_integrator_parallel(const t_inputrec *ir, const t_commrec *cr)
1280 if (((ir->eI) == eiSteep) || ((ir->eI) == eiCG) || ((ir->eI) == eiLBFGS) || ((ir->eI) == eiNM))
1282 gmx_fatal(FARGS, "%s Energy minimization via steep, CG, lbfgs and nm in parallel is currently not supported by IMD.\n", IMDstr);
1287 std::unique_ptr<ImdSession>
1288 makeImdSession(const t_inputrec *ir,
1289 const t_commrec *cr,
1290 gmx_wallcycle *wcycle,
1291 gmx_enerdata_t *enerd,
1292 const gmx_multisim_t *ms,
1293 const gmx_mtop_t *top_global,
1294 const MDLogger &mdlog,
1297 const t_filenm fnm[],
1298 const gmx_output_env_t *oenv,
1299 const MdrunOptions &mdrunOptions)
1301 std::unique_ptr<ImdSession> session(new ImdSession(mdlog));
1302 auto impl = session->impl_.get();
1304 /* We will allow IMD sessions only if supported by the binary and
1305 explicitly enabled in the .tpr file */
1306 if (!GMX_IMD || !ir->bIMD)
1311 // TODO As IMD is intended for interactivity, and the .tpr file
1312 // opted in for that, it is acceptable to write more terminal
1313 // output than in a typical simulation. However, all the GMX_LOG
1314 // statements below should go to both the log file and to the
1315 // terminal. This is probably be implemented by adding a logging
1316 // stream named like ImdInfo, to separate it from warning and to
1317 // send it to both destinations.
1319 if (EI_DYNAMICS(ir->eI))
1321 impl->defaultNstImd = ir->nstcalcenergy;
1323 else if (EI_ENERGY_MINIMIZATION(ir->eI))
1325 impl->defaultNstImd = 1;
1329 GMX_LOG(mdlog.warning).appendTextFormatted("%s Integrator '%s' is not supported for Interactive Molecular Dynamics, running normally instead", IMDstr, ei_names[ir->eI]);
1334 GMX_LOG(mdlog.warning).appendTextFormatted("%s Cannot use IMD for multiple simulations or replica exchange, running normally instead", IMDstr);
1338 const auto &options = mdrunOptions.imdOptions;
1340 bool createSession = false;
1341 /* It seems we have a .tpr file that defines an IMD group and thus allows IMD connections.
1342 * Check whether we can actually provide the IMD functionality for this setting: */
1345 /* Check whether IMD was enabled by one of the command line switches: */
1346 if (options.wait || options.terminatable || options.pull)
1348 GMX_LOG(mdlog.warning).appendTextFormatted("%s Enabled. This simulation will accept incoming IMD connections.", IMDstr);
1349 createSession = true;
1353 GMX_LOG(mdlog.warning).appendTextFormatted("%s None of the -imd switches was used.\n"
1354 "%s This run will not accept incoming IMD connections", IMDstr, IMDstr);
1356 } /* end master only */
1358 /* Let the other nodes know whether we want IMD */
1361 block_bc(cr, createSession);
1364 /*... if not we are done.*/
1371 /* check if we're using a sane integrator / parallel combination */
1372 imd_check_integrator_parallel(ir, cr);
1376 *****************************************************
1377 * From here on we assume that IMD is turned on *
1378 *****************************************************
1381 int nat_total = top_global->natoms;
1383 /* Initialize IMD session. If we read in a pre-IMD .tpr file, ir->imd->nat
1384 * will be zero. For those cases we transfer _all_ atomic positions */
1385 impl->sessionPossible = true;
1386 impl->nat = ir->imd->nat > 0 ? ir->imd->nat : nat_total;
1387 if (options.port >= 1)
1389 impl->port = options.port;
1392 impl->wcycle = wcycle;
1393 impl->enerd = enerd;
1395 /* We might need to open an output file for IMD forces data */
1398 impl->openOutputFile(opt2fn("-if", nfile, fnm), nat_total, oenv, mdrunOptions.continuationOptions);
1401 /* Make sure that we operate with a valid atom index array for the IMD atoms */
1402 if (ir->imd->nat > 0)
1404 /* Point to the user-supplied array of atom numbers */
1405 impl->ind = ir->imd->ind;
1409 /* Make a dummy (ind[i] = i) array of all atoms */
1410 snew(impl->ind, nat_total);
1411 for (int i = 0; i < nat_total; i++)
1417 /* read environment on master and prepare socket for incoming connections */
1420 /* we allocate memory for our IMD energy structure */
1421 int32_t recsize = c_headerSize + sizeof(IMDEnergyBlock);
1422 snew(impl->energysendbuf, recsize);
1424 /* Shall we wait for a connection? */
1427 impl->bWConnect = true;
1428 GMX_LOG(mdlog.warning).appendTextFormatted("%s Pausing simulation while no IMD connection present (-imdwait).", IMDstr);
1431 /* Will the IMD clients be able to terminate the simulation? */
1432 if (options.terminatable)
1434 impl->bTerminatable = true;
1435 GMX_LOG(mdlog.warning).appendTextFormatted("%s Allow termination of the simulation from IMD client (-imdterm).", IMDstr);
1438 /* Is pulling from IMD client allowed? */
1441 impl->bForceActivated = true;
1442 GMX_LOG(mdlog.warning).appendTextFormatted("%s Pulling from IMD remote is enabled (-imdpull).", IMDstr);
1445 /* Initialize send buffers with constant size */
1446 snew(impl->sendxbuf, impl->nat);
1447 snew(impl->energies, 1);
1448 int32_t bufxsize = c_headerSize + 3 * sizeof(float) * impl->nat;
1449 snew(impl->coordsendbuf, bufxsize);
1452 /* do we allow interactive pulling? If so let the other nodes know. */
1455 block_bc(cr, impl->bForceActivated);
1458 /* setup the listening socket on master process */
1461 GMX_LOG(mdlog.warning).appendTextFormatted("%s Setting port for connection requests to %d.", IMDstr, impl->port);
1462 impl->prepareMasterSocket();
1463 /* Wait until we have a connection if specified before */
1464 if (impl->bWConnect)
1466 impl->blockConnect();
1470 GMX_LOG(mdlog.warning).appendTextFormatted("%s -imdwait not set, starting simulation.", IMDstr);
1473 /* Let the other nodes know whether we are connected */
1474 impl->syncNodes(cr, 0);
1476 /* Initialize arrays used to assemble the positions from the other nodes */
1477 impl->prepareForPositionAssembly(cr, x);
1479 /* Initialize molecule blocks to make them whole later...*/
1482 impl->prepareMoleculesInImdGroup(top_global);
1489 bool ImdSession::Impl::run(int64_t step,
1496 if (!sessionPossible)
1501 wallcycle_start(wcycle, ewcIMD);
1503 /* read command from client and check if new incoming connection */
1506 /* If not already connected, check for new connections */
1519 /* Let's see if we have new IMD messages for us */
1526 /* is this an IMD communication step? */
1527 bool imdstep = do_per_step(step, nstimd);
1529 /* OK so this is an IMD step ... */
1532 /* First we sync all nodes to let everybody know whether we are connected to VMD */
1536 /* If a client is connected, we collect the positions
1537 * and put molecules back into the box before transfer */
1538 if ((imdstep && bConnected)
1539 || bNS) /* independent of imdstep, we communicate positions at each NS step */
1541 /* Transfer the IMD positions to the master node. Every node contributes
1542 * its local positions x and stores them in the assembled xa array. */
1543 communicate_group_positions(cr, xa, xa_shifts, xa_eshifts,
1544 true, x, nat, nat_loc,
1545 ind_loc, xa_ind, xa_old, box);
1547 /* If connected and master -> remove shifts */
1548 if ((imdstep && bConnected) && MASTER(cr))
1550 removeMolecularShifts(box);
1554 wallcycle_stop(wcycle, ewcIMD);
1559 bool ImdSession::run(int64_t step,
1565 return impl_->run(step, bNS, box, x, t);
1568 void ImdSession::fillEnergyRecord(int64_t step,
1569 bool bHaveNewEnergies)
1571 if (!impl_->sessionPossible || !impl_->clientsocket)
1576 IMDEnergyBlock *ene = impl_->energies;
1580 /* In MPI-parallel simulations the energies are not accessible a at every time step.
1581 * We update them if we have new values, otherwise, the energy values from the
1582 * last global communication step are still on display in the viewer. */
1583 if (bHaveNewEnergies)
1585 ene->T_abs = static_cast<float>(impl_->enerd->term[F_TEMP ]);
1586 ene->E_pot = static_cast<float>(impl_->enerd->term[F_EPOT ]);
1587 ene->E_tot = static_cast<float>(impl_->enerd->term[F_ETOT ]);
1588 ene->E_bond = static_cast<float>(impl_->enerd->term[F_BONDS ]);
1589 ene->E_angle = static_cast<float>(impl_->enerd->term[F_ANGLES ]);
1590 ene->E_dihe = static_cast<float>(impl_->enerd->term[F_PDIHS ]);
1591 ene->E_impr = static_cast<float>(impl_->enerd->term[F_IDIHS ]);
1592 ene->E_vdw = static_cast<float>(impl_->enerd->term[F_LJ ]);
1593 ene->E_coul = static_cast<float>(impl_->enerd->term[F_COUL_SR]);
1599 ImdSession::sendPositionsAndEnergies()
1601 if (!impl_->sessionPossible || !impl_->clientsocket)
1606 if (imd_send_energies(impl_->clientsocket, impl_->energies, impl_->energysendbuf))
1608 impl_->issueFatalError("Error sending updated energies. Disconnecting client.");
1611 if (imd_send_rvecs(impl_->clientsocket, impl_->nat, impl_->xa, impl_->coordsendbuf))
1613 impl_->issueFatalError("Error sending updated positions. Disconnecting client.");
1619 ImdSession::updateEnergyRecordAndSendPositionsAndEnergies(bool bIMDstep,
1621 bool bHaveNewEnergies)
1623 if (!impl_->sessionPossible)
1628 wallcycle_start(impl_->wcycle, ewcIMD);
1630 /* Update time step for IMD and prepare IMD energy record if we have new energies. */
1631 fillEnergyRecord(step, bHaveNewEnergies);
1635 /* Send positions and energies to VMD client via IMD */
1636 sendPositionsAndEnergies();
1639 wallcycle_stop(impl_->wcycle, ewcIMD);
1642 void ImdSession::applyForces(rvec *f)
1644 if (!impl_->sessionPossible || !impl_->bForceActivated)
1649 wallcycle_start(impl_->wcycle, ewcIMD);
1651 for (int i = 0; i < impl_->nforces; i++)
1653 /* j are the indices in the "System group".*/
1654 int j = impl_->ind[impl_->f_ind[i]];
1656 /* check if this is a local atom and find out locndx */
1658 const t_commrec *cr = impl_->cr;
1659 if (PAR(cr) && (locndx = cr->dd->ga2la->findHome(j)))
1664 rvec_inc(f[j], impl_->f[i]);
1667 wallcycle_stop(impl_->wcycle, ewcIMD);
1670 ImdSession::ImdSession(const MDLogger &mdlog)
1671 : impl_(new Impl(mdlog))
1674 ImdSession::~ImdSession() = default;