2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
5 * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
40 * Implements Interactive Molecular Dynamics
42 * Re-implementation of basic IMD functions to work with VMD,
43 * see imdsocket.h for references to the IMD API.
45 * \author Martin Hoefling, Carsten Kutzner <ckutzne@gwdg.de>
51 #include "gromacs/utility/enumerationhelpers.h"
59 #include "gromacs/commandline/filenm.h"
60 #include "gromacs/domdec/domdec_struct.h"
61 #include "gromacs/domdec/ga2la.h"
62 #include "gromacs/fileio/confio.h"
63 #include "gromacs/fileio/gmxfio.h"
64 #include "gromacs/fileio/xvgr.h"
65 #include "gromacs/gmxlib/network.h"
66 #include "gromacs/imd/imdsocket.h"
67 #include "gromacs/math/units.h"
68 #include "gromacs/math/vec.h"
69 #include "gromacs/mdlib/broadcaststructs.h"
70 #include "gromacs/mdlib/groupcoord.h"
71 #include "gromacs/mdlib/sighandler.h"
72 #include "gromacs/mdlib/stat.h"
73 #include "gromacs/mdrunutility/handlerestart.h"
74 #include "gromacs/mdrunutility/multisim.h"
75 #include "gromacs/mdtypes/commrec.h"
76 #include "gromacs/mdtypes/enerdata.h"
77 #include "gromacs/mdtypes/imdmodule.h"
78 #include "gromacs/mdtypes/inputrec.h"
79 #include "gromacs/mdtypes/md_enums.h"
80 #include "gromacs/mdtypes/mdrunoptions.h"
81 #include "gromacs/mdtypes/state.h"
82 #include "gromacs/pbcutil/pbc.h"
83 #include "gromacs/timing/wallcycle.h"
84 #include "gromacs/topology/mtop_util.h"
85 #include "gromacs/topology/topology.h"
86 #include "gromacs/utility/enumerationhelpers.h"
87 #include "gromacs/utility/fatalerror.h"
88 #include "gromacs/utility/logger.h"
89 #include "gromacs/utility/smalloc.h"
90 #include "gromacs/utility/stringutil.h"
95 /*! \brief How long shall we wait in seconds until we check for a connection again? */
96 constexpr int c_loopWait = 1;
98 /*! \brief How long shall we check for the IMD_GO? */
99 constexpr int c_connectWait = 1;
101 /*! \brief IMD Header Size. */
102 constexpr int c_headerSize = 8;
104 /*! \brief IMD Protocol Version. */
105 constexpr int c_protocolVersion = 2;
107 enum class IMDMessageType : int;
111 * IMD (interactive molecular dynamics) energy record.
113 * As in the original IMD implementation. Energies in kcal/mol.
114 * NOTE: We return the energies in GROMACS / SI units,
115 * so they also show up as SI in VMD.
120 int32_t tstep; /**< time step */
121 float T_abs; /**< absolute temperature */
122 float E_tot; /**< total energy */
123 float E_pot; /**< potential energy */
124 float E_vdw; /**< van der Waals energy */
125 float E_coul; /**< Coulomb interaction energy */
126 float E_bond; /**< bonds energy */
127 float E_angle; /**< angles energy */
128 float E_dihe; /**< dihedrals energy */
129 float E_impr; /**< improper dihedrals energy */
134 * \brief IMD (interactive molecular dynamics) communication structure.
136 * This structure defines the IMD communication message header & protocol version.
140 int32_t type; /**< Type of IMD message, see IMDType_t above */
141 int32_t length; /**< Length */
146 * \brief Implementation type for the IMD session
148 * \todo Make this class (or one extracted from it) model
150 * \todo Use RAII for files and allocations
152 class ImdSession::Impl
156 Impl(const MDLogger& mdlog);
159 /*! \brief Prepare the socket on the MASTER. */
160 void prepareMasterSocket();
161 /*! \brief Disconnect the client. */
162 void disconnectClient();
163 /*! \brief Prints an error message and disconnects the client.
165 * Does not terminate mdrun!
167 void issueFatalError(const char* msg);
168 /*! \brief Check whether we got an incoming connection. */
170 /*! \brief Wrap imd_tryconnect in order to make it blocking.
172 * Used when the simulation should wait for an incoming connection.
175 /*! \brief Make sure that our array holding the forces received via IMD is large enough. */
176 void prepareVmdForces();
177 /*! \brief Reads forces received via IMD. */
178 void readVmdForces();
179 /*! \brief Prepares the MD force arrays. */
180 void prepareMDForces();
181 /*! \brief Copy IMD forces to MD forces.
183 * Do conversion from Cal->Joule and from
184 * Angstrom -> nm and from a pointer array to arrays to 3*N array.
186 void copyToMDForces() const;
187 /*! \brief Return true if any of the forces or indices changed. */
188 bool bForcesChanged() const;
189 /*! \brief Update the old_f_ind and old_forces arrays to contain the current values. */
190 void keepOldValues();
191 /*! \brief Write the applied pull forces to logfile.
193 * Call on master only!
195 void outputForces(double time);
196 /*! \brief Synchronize the nodes. */
197 void syncNodes(const t_commrec* cr, double t);
198 /*! \brief Reads header from the client and decides what to do. */
200 /*! \brief Open IMD output file and write header information.
202 * Call on master only.
204 void openOutputFile(const char* fn, int nat_total, const gmx_output_env_t* oenv, StartingBehavior startingBehavior);
205 /*! \brief Creates the molecule start-end position array of molecules in the IMD group. */
206 void prepareMoleculesInImdGroup(const gmx_mtop_t& top_global);
207 /*! \brief Removes shifts of molecules diffused outside of the box. */
208 void removeMolecularShifts(const matrix box) const;
209 /*! \brief Initialize arrays used to assemble the positions from the other nodes. */
210 void prepareForPositionAssembly(const t_commrec* cr, gmx::ArrayRef<const gmx::RVec> coords);
211 /*! \brief Interact with any connected VMD session */
212 bool run(int64_t step, bool bNS, const matrix box, gmx::ArrayRef<const gmx::RVec> coords, double t);
214 // TODO rename all the data members to have underscore suffixes
216 //! True if tpr and mdrun input combine to permit IMD sessions
217 bool sessionPossible = false;
218 //! Output file for IMD data, mainly forces.
219 FILE* outf = nullptr;
221 //! Number of atoms that can be pulled via IMD.
223 //! Part of the atoms that are local.
225 //! Global indices of the IMD atoms.
227 //! Local indices of the IMD atoms.
228 int* ind_loc = nullptr;
229 //! Allocation size for ind_loc.
231 //! Positions for all IMD atoms assembled on the master node.
233 //! Shifts for all IMD atoms, to make molecule(s) whole.
234 ivec* xa_shifts = nullptr;
235 //! Extra shifts since last DD step.
236 ivec* xa_eshifts = nullptr;
237 //! Old positions for all IMD atoms on master.
238 rvec* xa_old = nullptr;
239 //! Position of each local atom in the collective array.
240 int* xa_ind = nullptr;
242 //! Global IMD frequency, known to all ranks.
244 //! New frequency from IMD client, master only.
246 //! Default IMD frequency when disconnected.
247 int defaultNstImd = -1;
249 //! Port to use for network socket.
251 //! The IMD socket on the master node.
252 IMDSocket* socket = nullptr;
253 //! The IMD socket on the client.
254 IMDSocket* clientsocket = nullptr;
255 //! Length we got with last header.
258 //! Shall we block and wait for connection?
259 bool bWConnect = false;
260 //! Set if MD is terminated.
261 bool bTerminated = false;
262 //! Set if MD can be terminated.
263 bool bTerminatable = false;
264 //! Set if connection is present.
265 bool bConnected = false;
266 //! Set if we received new forces.
267 bool bNewForces = false;
268 //! Set if pulling from VMD is allowed.
269 bool bForceActivated = false;
271 //! Pointer to energies we send back.
272 IMDEnergyBlock* energies = nullptr;
274 //! Number of VMD forces.
275 int32_t vmd_nforces = 0;
276 //! VMD forces indices.
277 int32_t* vmd_f_ind = nullptr;
278 //! The VMD forces flat in memory.
279 float* vmd_forces = nullptr;
280 //! Number of actual MD forces; this gets communicated to the clients.
283 int* f_ind = nullptr;
284 //! The IMD pulling forces.
287 //! Buffer for force sending.
288 char* forcesendbuf = nullptr;
289 //! Buffer for coordinate sending.
290 char* coordsendbuf = nullptr;
291 //! Send buffer for energies.
292 char* energysendbuf = nullptr;
293 //! Buffer to make molecules whole before sending.
294 rvec* sendxbuf = nullptr;
296 //! Molecules block in IMD group.
299 /* The next block is used on the master node only to reduce the output
300 * without sacrificing information. If any of these values changes,
301 * we need to write output */
302 //! Old value for nforces.
304 //! Old values for force indices.
305 int* old_f_ind = nullptr;
306 //! Old values for IMD pulling forces.
307 rvec* old_forces = nullptr;
310 const MDLogger& mdlog;
311 //! Commmunication object
312 const t_commrec* cr = nullptr;
313 //! Wallcycle counting manager.
314 gmx_wallcycle* wcycle = nullptr;
315 //! Energy output handler
316 gmx_enerdata_t* enerd = nullptr;
320 * \brief Implement interactive molecular dynamics.
322 * \todo Some aspects of this module provides forces (when the user
323 * pulls on things in VMD), so in future it should have a class that
324 * models IForceProvider and is contributed to the collection of such
327 class InteractiveMolecularDynamics final : public IMDModule
330 IMdpOptionProvider* mdpOptionProvider() override { return nullptr; }
331 IMDOutputProvider* outputProvider() override { return nullptr; }
332 void initForceProviders(ForceProviders* /* forceProviders */) override {}
333 void subscribeToSimulationSetupNotifications(MDModulesNotifiers* /* notifiers */) override {}
334 void subscribeToPreProcessingNotifications(MDModulesNotifiers* /* notifiers */) override {}
337 std::unique_ptr<IMDModule> createInteractiveMolecularDynamicsModule()
339 return std::make_unique<InteractiveMolecularDynamics>();
342 /*! \brief Enum for types of IMD messages.
344 * We use the same records as the NAMD/VMD IMD implementation.
346 enum class IMDMessageType : int
348 Disconnect, /**< client disconnect */
349 Energies, /**< energy data */
350 FCoords, /**< atomic coordinates */
351 Go, /**< start command for the simulation */
352 Handshake, /**< handshake to determine little/big endianness */
353 Kill, /**< terminates the simulation */
354 Mdcomm, /**< force data */
355 Pause, /**< pauses the simulation */
356 TRate, /**< sets the IMD transmission and processing rate */
357 IOerror, /**< I/O error */
358 Count /**< number of entries */
362 /*! \brief Names of the IMDType for error messages.
364 static const char* enumValueToString(IMDMessageType enumValue)
366 constexpr gmx::EnumerationArray<IMDMessageType, const char*> imdMessageTypeNames = {
367 "IMD_DISCONNECT", "IMD_ENERGIES", "IMD_FCOORDS", "IMD_GO", "IMD_HANDSHAKE",
368 "IMD_KILL", "IMD_MDCOMM", "IMD_PAUSE", "IMD_TRATE", "IMD_IOERROR"
370 return imdMessageTypeNames[enumValue];
374 /*! \brief Fills the header with message and the length argument. */
375 static void fill_header(IMDHeader* header, IMDMessageType type, int32_t length)
377 /* We (ab-)use htonl network function for the correct endianness */
378 header->type = imd_htonl(static_cast<int>(type));
379 header->length = imd_htonl(length);
383 /*! \brief Swaps the endianess of the header. */
384 static void swap_header(IMDHeader* header)
386 /* and vice versa... */
387 header->type = imd_ntohl(static_cast<int>(header->type));
388 header->length = imd_ntohl(header->length);
392 /*! \brief Reads multiple bytes from socket. */
393 static int32_t imd_read_multiple(IMDSocket* socket, char* datptr, int32_t toread)
395 int32_t leftcount, countread;
399 /* Try to read while we haven't reached toread */
400 while (leftcount != 0)
402 if ((countread = imdsock_read(socket, datptr, leftcount)) < 0)
404 /* interrupted function call, try again... */
409 /* this is an unexpected error, return what we got */
412 return toread - leftcount;
415 /* nothing read, finished */
417 else if (countread == 0)
421 leftcount -= countread;
425 /* return nr of bytes read */
426 return toread - leftcount;
430 /*! \brief Writes multiple bytes to socket in analogy to imd_read_multiple. */
431 static int32_t imd_write_multiple(IMDSocket* socket, const char* datptr, int32_t towrite)
433 int32_t leftcount, countwritten;
437 while (leftcount != 0)
439 if ((countwritten = imdsock_write(socket, datptr, leftcount)) <= 0)
447 return towrite - leftcount;
450 leftcount -= countwritten;
451 datptr += countwritten;
454 return towrite - leftcount;
458 /*! \brief Handshake with IMD client. */
459 static int imd_handshake(IMDSocket* socket)
464 fill_header(&header, IMDMessageType::Handshake, 1);
465 header.length = c_protocolVersion; /* client wants unswapped version */
467 return static_cast<int>(imd_write_multiple(socket, reinterpret_cast<char*>(&header), c_headerSize)
472 /*! \brief Send energies using the energy block and the send buffer. */
473 static int imd_send_energies(IMDSocket* socket, const IMDEnergyBlock* energies, char* buffer)
478 recsize = c_headerSize + sizeof(IMDEnergyBlock);
479 fill_header(reinterpret_cast<IMDHeader*>(buffer), IMDMessageType::Energies, 1);
480 memcpy(buffer + c_headerSize, energies, sizeof(IMDEnergyBlock));
482 return static_cast<int>(imd_write_multiple(socket, buffer, recsize) != recsize);
486 /*! \brief Receive IMD header from socket, sets the length and returns the IMD message. */
487 static IMDMessageType imd_recv_header(IMDSocket* socket, int32_t* length)
492 if (imd_read_multiple(socket, reinterpret_cast<char*>(&header), c_headerSize) != c_headerSize)
494 return IMDMessageType::IOerror;
496 swap_header(&header);
497 *length = header.length;
499 return static_cast<IMDMessageType>(header.type);
503 /*! \brief Receive force indices and forces.
505 * The number of forces was previously communicated via the header.
507 static bool imd_recv_mdcomm(IMDSocket* socket, int32_t nforces, int32_t* forcendx, float* forces)
509 /* reading indices */
510 int retsize = sizeof(int32_t) * nforces;
511 int retbytes = imd_read_multiple(socket, reinterpret_cast<char*>(forcendx), retsize);
512 if (retbytes != retsize)
517 /* reading forces as float array */
518 retsize = 3 * sizeof(float) * nforces;
519 retbytes = imd_read_multiple(socket, reinterpret_cast<char*>(forces), retsize);
520 return (retbytes == retsize);
523 /* GROMACS specific functions for the IMD implementation */
524 void write_IMDgroup_to_file(bool bIMD,
526 const t_state* state,
527 const gmx_mtop_t& sys,
529 const t_filenm fnm[])
536 IMDatoms = gmx_mtop_global_atoms(sys);
537 write_sto_conf_indexed(opt2fn("-imd", nfile, fnm),
540 state->x.rvec_array(),
541 state->v.rvec_array(),
550 void ImdSession::dd_make_local_IMD_atoms(const gmx_domdec_t* dd)
552 if (!impl_->sessionPossible)
557 dd_make_local_group_indices(
558 dd->ga2la, impl_->nat, impl_->ind, &impl_->nat_loc, &impl_->ind_loc, &impl_->nalloc_loc, impl_->xa_ind);
562 /*! \brief Send positions from rvec.
564 * We need a separate send buffer and conversion to Angstrom.
566 static int imd_send_rvecs(IMDSocket* socket, int nat, rvec* x, char* buffer)
571 int tuplesize = 3 * sizeof(float);
574 /* Required size for the send buffer */
575 size = c_headerSize + 3 * sizeof(float) * nat;
578 fill_header(reinterpret_cast<IMDHeader*>(buffer), IMDMessageType::FCoords, static_cast<int32_t>(nat));
579 for (i = 0; i < nat; i++)
581 sendx[0] = static_cast<float>(x[i][0]) * gmx::c_nm2A;
582 sendx[1] = static_cast<float>(x[i][1]) * gmx::c_nm2A;
583 sendx[2] = static_cast<float>(x[i][2]) * gmx::c_nm2A;
584 memcpy(buffer + c_headerSize + i * tuplesize, sendx, tuplesize);
587 return static_cast<int>(imd_write_multiple(socket, buffer, size) != size);
591 void ImdSession::Impl::prepareMasterSocket()
593 if (imdsock_winsockinit() == -1)
595 gmx_fatal(FARGS, "%s Failed to initialize winsock.\n", IMDstr);
598 /* The rest is identical, first create and bind a socket and set to listen then. */
599 GMX_LOG(mdlog.warning).appendTextFormatted("%s Setting up incoming socket.", IMDstr);
600 socket = imdsock_create();
603 gmx_fatal(FARGS, "%s Failed to create socket.", IMDstr);
607 int ret = imdsock_bind(socket, port);
610 gmx_fatal(FARGS, "%s binding socket to port %d failed with error %d.\n", IMDstr, port, ret);
613 if (imd_sock_listen(socket))
615 gmx_fatal(FARGS, "%s socket listen failed with error %d.\n", IMDstr, ret);
618 if (imdsock_getport(socket, &port))
620 gmx_fatal(FARGS, "%s Could not determine port number.\n", IMDstr);
623 GMX_LOG(mdlog.warning).appendTextFormatted("%s Listening for IMD connection on port %d.", IMDstr, port);
627 void ImdSession::Impl::disconnectClient()
629 /* Write out any buffered pulling data */
632 /* we first try to shut down the clientsocket */
633 imdsock_shutdown(clientsocket);
634 if (!imdsock_destroy(clientsocket))
636 GMX_LOG(mdlog.warning).appendTextFormatted("%s Failed to destroy socket.", IMDstr);
639 /* then we reset the IMD step to its default, and reset the connection boolean */
640 nstimd_new = defaultNstImd;
641 clientsocket = nullptr;
646 void ImdSession::Impl::issueFatalError(const char* msg)
648 GMX_LOG(mdlog.warning).appendTextFormatted("%s %s", IMDstr, msg);
650 GMX_LOG(mdlog.warning).appendTextFormatted("%s disconnected.", IMDstr);
654 bool ImdSession::Impl::tryConnect()
656 if (imdsock_tryread(socket, 0, 0) > 0)
658 /* yes, we got something, accept on clientsocket */
659 clientsocket = imdsock_accept(socket);
662 GMX_LOG(mdlog.warning)
663 .appendTextFormatted("%s Accepting the connection on the socket failed.", IMDstr);
667 /* handshake with client */
668 if (imd_handshake(clientsocket))
670 issueFatalError("Connection failed.");
674 GMX_LOG(mdlog.warning)
675 .appendTextFormatted("%s Connection established, checking if I got IMD_GO orders.", IMDstr);
677 /* Check if we get the proper "GO" command from client. */
678 if (imdsock_tryread(clientsocket, c_connectWait, 0) != 1
679 || imd_recv_header(clientsocket, &(length)) != IMDMessageType::Go)
681 issueFatalError("No IMD_GO order received. IMD connection failed.");
694 void ImdSession::Impl::blockConnect()
696 /* do not wait for connection, when e.g. ctrl+c is pressed and we will terminate anyways. */
697 if (!(gmx_get_stop_condition() == StopCondition::None))
702 GMX_LOG(mdlog.warning)
703 .appendTextFormatted("%s Will wait until I have a connection and IMD_GO orders.", IMDstr);
705 /* while we have no clientsocket... 2nd part: we should still react on ctrl+c */
706 while ((!clientsocket) && (gmx_get_stop_condition() == StopCondition::None))
709 imd_sleep(c_loopWait);
714 void ImdSession::Impl::prepareVmdForces()
716 srenew((vmd_f_ind), vmd_nforces);
717 srenew((vmd_forces), 3 * vmd_nforces);
721 void ImdSession::Impl::readVmdForces()
723 /* the length of the previously received header tells us the nr of forces we will receive */
724 vmd_nforces = length;
725 /* prepare the arrays */
727 /* Now we read the forces... */
728 if (!(imd_recv_mdcomm(clientsocket, vmd_nforces, vmd_f_ind, vmd_forces)))
730 issueFatalError("Error while reading forces from remote. Disconnecting");
735 void ImdSession::Impl::prepareMDForces()
737 srenew((f_ind), nforces);
738 srenew((f), nforces);
742 void ImdSession::Impl::copyToMDForces() const
745 real conversion = gmx::c_cal2Joule * gmx::c_nm2A;
748 for (i = 0; i < nforces; i++)
750 /* Copy the indices, a copy is important because we may update the incoming forces
751 * whenever we receive new forces while the MD forces are only communicated upon
752 * IMD communication.*/
753 f_ind[i] = vmd_f_ind[i];
755 /* Convert to rvecs and do a proper unit conversion */
756 f[i][0] = vmd_forces[3 * i] * conversion;
757 f[i][1] = vmd_forces[3 * i + 1] * conversion;
758 f[i][2] = vmd_forces[3 * i + 2] * conversion;
763 /*! \brief Returns true if any component of the two rvecs differs. */
764 static inline bool rvecs_differ(const rvec v1, const rvec v2)
766 for (int i = 0; i < DIM; i++)
777 bool ImdSession::Impl::bForcesChanged() const
779 /* First, check whether the number of pulled atoms changed */
780 if (nforces != old_nforces)
785 /* Second, check whether any of the involved atoms changed */
786 for (int i = 0; i < nforces; i++)
788 if (f_ind[i] != old_f_ind[i])
794 /* Third, check whether all forces are the same */
795 for (int i = 0; i < nforces; i++)
797 if (rvecs_differ(f[i], old_forces[i]))
803 /* All old and new forces are identical! */
808 void ImdSession::Impl::keepOldValues()
810 old_nforces = nforces;
812 for (int i = 0; i < nforces; i++)
814 old_f_ind[i] = f_ind[i];
815 copy_rvec(f[i], old_forces[i]);
820 void 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]);
846 void ImdSession::Impl::syncNodes(const t_commrec* cr, double t)
848 /* Notify the other nodes whether we are still connected. */
851 block_bc(cr->mpi_comm_mygroup, bConnected);
854 /* ...if not connected, the job is done here. */
860 /* Let the other nodes know whether we got a new IMD synchronization frequency. */
863 block_bc(cr->mpi_comm_mygroup, nstimd_new);
866 /* Now we all set the (new) nstimd communication time step */
869 /* We're done if we don't allow pulling at all */
870 if (!(bForceActivated))
875 /* OK, let's check if we have received forces which we need to communicate
876 * to the other nodes */
882 new_nforces = vmd_nforces;
884 /* make the "new_forces" negative, when we did not receive new ones */
887 new_nforces = vmd_nforces * -1;
891 /* make new_forces known to the clients */
894 block_bc(cr->mpi_comm_mygroup, new_nforces);
897 /* When new_natoms < 0 then we know that these are still the same forces
898 * so we don't communicate them, otherwise... */
904 /* set local VMD and nforces */
905 vmd_nforces = new_nforces;
906 nforces = vmd_nforces;
908 /* now everybody knows the number of forces in f_ind, so we can prepare
909 * the target arrays for indices and forces */
912 /* we first update the MD forces on the master by converting the VMD forces */
916 /* We also write out forces on every update, so that we know which
917 * forces are applied for every step */
924 /* In parallel mode we communicate the to-be-applied forces to the other nodes */
927 nblock_bc(cr->mpi_comm_mygroup, nforces, f_ind);
928 nblock_bc(cr->mpi_comm_mygroup, nforces, f);
931 /* done communicating the forces, reset bNewForces */
936 void ImdSession::Impl::readCommand()
938 bool IMDpaused = false;
940 while (clientsocket && (imdsock_tryread(clientsocket, 0, 0) > 0 || IMDpaused))
942 IMDMessageType itype = imd_recv_header(clientsocket, &(length));
943 /* let's see what we got: */
946 /* IMD asks us to terminate the simulation, check if the user allowed this */
947 case IMDMessageType::Kill:
950 GMX_LOG(mdlog.warning)
951 .appendTextFormatted(
952 " %s Terminating connection and running simulation (if "
953 "supported by integrator).",
957 gmx_set_stop_condition(StopCondition::Next);
961 GMX_LOG(mdlog.warning)
962 .appendTextFormatted(
963 " %s Set -imdterm command line switch to allow mdrun "
964 "termination from within IMD.",
970 /* the client doen't want to talk to us anymore */
971 case IMDMessageType::Disconnect:
972 GMX_LOG(mdlog.warning).appendTextFormatted(" %s Disconnecting client.", IMDstr);
976 /* we got new forces, read them and set bNewForces flag */
977 case IMDMessageType::Mdcomm:
982 /* the client asks us to (un)pause the simulation. So we toggle the IMDpaused state */
983 case IMDMessageType::Pause:
986 GMX_LOG(mdlog.warning).appendTextFormatted(" %s Un-pause command received.", IMDstr);
991 GMX_LOG(mdlog.warning).appendTextFormatted(" %s Pause command received.", IMDstr);
997 /* the client sets a new transfer rate, if we get 0, we reset the rate
998 * to the default. VMD filters 0 however */
999 case IMDMessageType::TRate:
1000 nstimd_new = (length > 0) ? length : defaultNstImd;
1001 GMX_LOG(mdlog.warning)
1002 .appendTextFormatted(" %s Update frequency will be set to %d.", IMDstr, nstimd_new);
1005 /* Catch all rule for the remaining IMD types which we don't expect */
1007 GMX_LOG(mdlog.warning)
1008 .appendTextFormatted(
1009 " %s Received unexpected %s.", IMDstr, enumValueToString(itype));
1010 issueFatalError("Terminating connection");
1017 void ImdSession::Impl::openOutputFile(const char* fn,
1019 const gmx_output_env_t* oenv,
1020 const gmx::StartingBehavior startingBehavior)
1022 /* Open log file of applied IMD forces if requested */
1026 "%s For a log of the IMD pull forces explicitly specify '-if' on the command "
1028 "%s (Not possible with energy minimization.)\n",
1034 /* If we append to an existing file, all the header information is already there */
1035 if (startingBehavior == StartingBehavior::RestartWithAppending)
1037 outf = gmx_fio_fopen(fn, "a+");
1041 outf = gmx_fio_fopen(fn, "w+");
1042 if (nat == nat_total)
1045 "# Note that you can select an IMD index group in the .mdp file if a subset of "
1046 "the atoms suffices.\n");
1049 xvgr_header(outf, "IMD Pull Forces", "Time (ps)", "# of Forces / Atom IDs / Forces (kJ/mol)", exvggtNONE, oenv);
1051 fprintf(outf, "# Can display and manipulate %d (of a total of %d) atoms via IMD.\n", nat, nat_total);
1052 fprintf(outf, "# column 1 : time (ps)\n");
1054 "# column 2 : total number of atoms feeling an IMD pulling force at that "
1057 "# cols. 3.-6 : global atom number of pulled atom, x-force, y-force, z-force "
1060 "# then follow : atom-ID, f[x], f[y], f[z] for more atoms in case the force on "
1061 "multiple atoms is changed simultaneously.\n");
1063 "# Note that the force on any atom is always equal to the last value for that "
1064 "atom-ID found in the data.\n");
1068 /* To reduce the output file size we remember the old values and output only
1069 * when something changed */
1070 snew(old_f_ind, nat); /* One can never pull on more atoms */
1071 snew(old_forces, nat);
1075 ImdSession::Impl::Impl(const MDLogger& mdlog) : mdlog(mdlog)
1080 ImdSession::Impl::~Impl()
1084 gmx_fio_fclose(outf);
1090 void ImdSession::Impl::prepareMoleculesInImdGroup(const gmx_mtop_t& top_global)
1092 /* check whether index is sorted */
1093 for (int i = 0; i < nat - 1; i++)
1095 if (ind[i] > ind[i + 1])
1097 gmx_fatal(FARGS, "%s IMD index is not sorted. This is currently not supported.\n", IMDstr);
1101 RangePartitioning gmols = gmx_mtop_molecules(top_global);
1104 snew(lmols.index, gmols.numBlocks() + 1);
1107 for (int i = 0; i < gmols.numBlocks(); i++)
1109 auto mol = gmols.block(i);
1111 for (int ii = 0; ii < nat; ii++)
1113 if (mol.isInRange(ind[ii]))
1120 lmols.index[lmols.nr + 1] = lmols.index[lmols.nr] + count;
1124 srenew(lmols.index, lmols.nr + 1);
1125 lmols.nalloc_index = lmols.nr + 1;
1130 /*! \brief Copied and modified from groupcoord.c shift_positions_group(). */
1131 static void shift_positions(const matrix box,
1132 rvec x[], /* The positions [0..nr] */
1133 const ivec is, /* The shift [0..nr] */
1134 int nr) /* The number of positions */
1138 /* Loop over the group's atoms */
1141 for (i = 0; i < nr; i++)
1147 x[i][XX] = x[i][XX] - tx * box[XX][XX] - ty * box[YY][XX] - tz * box[ZZ][XX];
1148 x[i][YY] = x[i][YY] - ty * box[YY][YY] - tz * box[ZZ][YY];
1149 x[i][ZZ] = x[i][ZZ] - tz * box[ZZ][ZZ];
1154 for (i = 0; i < nr; i++)
1160 x[i][XX] = x[i][XX] - tx * box[XX][XX];
1161 x[i][YY] = x[i][YY] - ty * box[YY][YY];
1162 x[i][ZZ] = x[i][ZZ] - tz * box[ZZ][ZZ];
1168 void ImdSession::Impl::removeMolecularShifts(const matrix box) const
1170 /* for each molecule also present in IMD group */
1171 for (int i = 0; i < mols.nr; i++)
1173 /* first we determine the minimum and maximum shifts for each molecule */
1175 ivec largest, smallest, shift;
1176 clear_ivec(largest);
1177 clear_ivec(smallest);
1180 copy_ivec(xa_shifts[mols.index[i]], largest);
1181 copy_ivec(xa_shifts[mols.index[i]], smallest);
1183 for (int ii = mols.index[i] + 1; ii < mols.index[i + 1]; ii++)
1185 if (xa_shifts[ii][XX] > largest[XX])
1187 largest[XX] = xa_shifts[ii][XX];
1189 if (xa_shifts[ii][XX] < smallest[XX])
1191 smallest[XX] = xa_shifts[ii][XX];
1194 if (xa_shifts[ii][YY] > largest[YY])
1196 largest[YY] = xa_shifts[ii][YY];
1198 if (xa_shifts[ii][YY] < smallest[YY])
1200 smallest[YY] = xa_shifts[ii][YY];
1203 if (xa_shifts[ii][ZZ] > largest[ZZ])
1205 largest[ZZ] = xa_shifts[ii][ZZ];
1207 if (xa_shifts[ii][ZZ] < smallest[ZZ])
1209 smallest[ZZ] = xa_shifts[ii][ZZ];
1213 /* check if we what we can subtract/add to the positions
1214 * to put them back into the central box */
1215 if (smallest[XX] > 0)
1217 shift[XX] = smallest[XX];
1219 if (smallest[YY] > 0)
1221 shift[YY] = smallest[YY];
1223 if (smallest[ZZ] > 0)
1225 shift[ZZ] = smallest[ZZ];
1228 if (largest[XX] < 0)
1230 shift[XX] = largest[XX];
1232 if (largest[YY] < 0)
1234 shift[YY] = largest[YY];
1236 if (largest[ZZ] < 0)
1238 shift[ZZ] = largest[ZZ];
1241 /* is there a shift at all? */
1242 if ((shift[XX]) || (shift[YY]) || (shift[ZZ]))
1244 int molsize = mols.index[i + 1] - mols.index[i];
1245 /* shift the positions */
1246 shift_positions(box, &(xa[mols.index[i]]), shift, molsize);
1252 void ImdSession::Impl::prepareForPositionAssembly(const t_commrec* cr, gmx::ArrayRef<const gmx::RVec> coords)
1256 snew(xa_shifts, nat);
1257 snew(xa_eshifts, nat);
1260 /* Save the original (whole) set of positions such that later the
1261 * molecule can always be made whole again */
1264 for (int i = 0; i < nat; i++)
1267 copy_rvec(coords[ii], xa_old[i]);
1271 if (!DOMAINDECOMP(cr))
1276 /* xa_ind[i] needs to be set to i for serial runs */
1277 for (int i = 0; i < nat; i++)
1283 /* Communicate initial coordinates xa_old to all processes */
1284 if (cr && havePPDomainDecomposition(cr))
1286 gmx_bcast(nat * sizeof(xa_old[0]), xa_old, cr->mpi_comm_mygroup);
1291 /*! \brief Check for non-working integrator / parallel options. */
1292 static void imd_check_integrator_parallel(const t_inputrec* ir, const t_commrec* cr)
1296 if (((ir->eI) == IntegrationAlgorithm::Steep) || ((ir->eI) == IntegrationAlgorithm::CG)
1297 || ((ir->eI) == IntegrationAlgorithm::LBFGS) || ((ir->eI) == IntegrationAlgorithm::NM))
1300 "%s Energy minimization via steep, CG, lbfgs and nm in parallel is currently "
1301 "not supported by IMD.\n",
1307 std::unique_ptr<ImdSession> makeImdSession(const t_inputrec* ir,
1308 const t_commrec* cr,
1309 gmx_wallcycle* wcycle,
1310 gmx_enerdata_t* enerd,
1311 const gmx_multisim_t* ms,
1312 const gmx_mtop_t& top_global,
1313 const MDLogger& mdlog,
1314 gmx::ArrayRef<const gmx::RVec> coords,
1316 const t_filenm fnm[],
1317 const gmx_output_env_t* oenv,
1318 const ImdOptions& options,
1319 const gmx::StartingBehavior startingBehavior)
1321 std::unique_ptr<ImdSession> session(new ImdSession(mdlog));
1322 auto* impl = session->impl_.get();
1324 /* We will allow IMD sessions only if supported by the binary and
1325 explicitly enabled in the .tpr file */
1326 if (!GMX_IMD || !ir->bIMD)
1331 // TODO As IMD is intended for interactivity, and the .tpr file
1332 // opted in for that, it is acceptable to write more terminal
1333 // output than in a typical simulation. However, all the GMX_LOG
1334 // statements below should go to both the log file and to the
1335 // terminal. This is probably be implemented by adding a logging
1336 // stream named like ImdInfo, to separate it from warning and to
1337 // send it to both destinations.
1339 if (EI_DYNAMICS(ir->eI))
1341 impl->defaultNstImd = ir->nstcalcenergy;
1343 else if (EI_ENERGY_MINIMIZATION(ir->eI))
1345 impl->defaultNstImd = 1;
1349 GMX_LOG(mdlog.warning)
1350 .appendTextFormatted(
1351 "%s Integrator '%s' is not supported for Interactive Molecular Dynamics, "
1352 "running normally instead",
1354 enumValueToString(ir->eI));
1359 GMX_LOG(mdlog.warning)
1360 .appendTextFormatted(
1361 "%s Cannot use IMD for multiple simulations or replica exchange, running "
1367 bool createSession = false;
1368 /* It seems we have a .tpr file that defines an IMD group and thus allows IMD connections.
1369 * Check whether we can actually provide the IMD functionality for this setting: */
1372 /* Check whether IMD was enabled by one of the command line switches: */
1373 if (options.wait || options.terminatable || options.pull)
1375 GMX_LOG(mdlog.warning)
1376 .appendTextFormatted(
1377 "%s Enabled. This simulation will accept incoming IMD connections.", IMDstr);
1378 createSession = true;
1382 GMX_LOG(mdlog.warning)
1383 .appendTextFormatted(
1384 "%s None of the -imd switches was used.\n"
1385 "%s This run will not accept incoming IMD connections",
1389 } /* end master only */
1391 /* Let the other nodes know whether we want IMD */
1394 block_bc(cr->mpi_comm_mygroup, createSession);
1397 /*... if not we are done.*/
1404 /* check if we're using a sane integrator / parallel combination */
1405 imd_check_integrator_parallel(ir, cr);
1409 *****************************************************
1410 * From here on we assume that IMD is turned on *
1411 *****************************************************
1414 int nat_total = top_global.natoms;
1416 /* Initialize IMD session. If we read in a pre-IMD .tpr file, ir->imd->nat
1417 * will be zero. For those cases we transfer _all_ atomic positions */
1418 impl->sessionPossible = true;
1419 impl->nat = ir->imd->nat > 0 ? ir->imd->nat : nat_total;
1420 if (options.port >= 1)
1422 impl->port = options.port;
1425 impl->wcycle = wcycle;
1426 impl->enerd = enerd;
1428 /* We might need to open an output file for IMD forces data */
1431 impl->openOutputFile(opt2fn("-if", nfile, fnm), nat_total, oenv, startingBehavior);
1434 /* Make sure that we operate with a valid atom index array for the IMD atoms */
1435 if (ir->imd->nat > 0)
1437 /* Point to the user-supplied array of atom numbers */
1438 impl->ind = ir->imd->ind;
1442 /* Make a dummy (ind[i] = i) array of all atoms */
1443 snew(impl->ind, nat_total);
1444 for (int i = 0; i < nat_total; i++)
1450 /* read environment on master and prepare socket for incoming connections */
1453 /* we allocate memory for our IMD energy structure */
1454 int32_t recsize = c_headerSize + sizeof(IMDEnergyBlock);
1455 snew(impl->energysendbuf, recsize);
1457 /* Shall we wait for a connection? */
1460 impl->bWConnect = true;
1461 GMX_LOG(mdlog.warning)
1462 .appendTextFormatted(
1463 "%s Pausing simulation while no IMD connection present (-imdwait).", IMDstr);
1466 /* Will the IMD clients be able to terminate the simulation? */
1467 if (options.terminatable)
1469 impl->bTerminatable = true;
1470 GMX_LOG(mdlog.warning)
1471 .appendTextFormatted(
1472 "%s Allow termination of the simulation from IMD client (-imdterm).", IMDstr);
1475 /* Is pulling from IMD client allowed? */
1478 impl->bForceActivated = true;
1479 GMX_LOG(mdlog.warning)
1480 .appendTextFormatted("%s Pulling from IMD remote is enabled (-imdpull).", IMDstr);
1483 /* Initialize send buffers with constant size */
1484 snew(impl->sendxbuf, impl->nat);
1485 snew(impl->energies, 1);
1486 int32_t bufxsize = c_headerSize + 3 * sizeof(float) * impl->nat;
1487 snew(impl->coordsendbuf, bufxsize);
1490 /* do we allow interactive pulling? If so let the other nodes know. */
1493 block_bc(cr->mpi_comm_mygroup, impl->bForceActivated);
1496 /* setup the listening socket on master process */
1499 GMX_LOG(mdlog.warning).appendTextFormatted("%s Setting port for connection requests to %d.", IMDstr, impl->port);
1500 impl->prepareMasterSocket();
1501 /* Wait until we have a connection if specified before */
1502 if (impl->bWConnect)
1504 impl->blockConnect();
1508 GMX_LOG(mdlog.warning).appendTextFormatted("%s -imdwait not set, starting simulation.", IMDstr);
1511 /* Let the other nodes know whether we are connected */
1512 impl->syncNodes(cr, 0);
1514 /* Initialize arrays used to assemble the positions from the other nodes */
1515 impl->prepareForPositionAssembly(cr, coords);
1517 /* Initialize molecule blocks to make them whole later...*/
1520 impl->prepareMoleculesInImdGroup(top_global);
1527 bool ImdSession::Impl::run(int64_t step, bool bNS, const matrix box, gmx::ArrayRef<const gmx::RVec> coords, double t)
1530 if (!sessionPossible)
1535 wallcycle_start(wcycle, WallCycleCounter::Imd);
1537 /* read command from client and check if new incoming connection */
1540 /* If not already connected, check for new connections */
1553 /* Let's see if we have new IMD messages for us */
1560 /* is this an IMD communication step? */
1561 bool imdstep = do_per_step(step, nstimd);
1563 /* OK so this is an IMD step ... */
1566 /* First we sync all nodes to let everybody know whether we are connected to VMD */
1570 /* If a client is connected, we collect the positions
1571 * and put molecules back into the box before transfer */
1572 if ((imdstep && bConnected) || bNS) /* independent of imdstep, we communicate positions at each NS step */
1574 /* Transfer the IMD positions to the master node. Every node contributes
1575 * its local positions x and stores them in the assembled xa array. */
1576 communicate_group_positions(
1577 cr, xa, xa_shifts, xa_eshifts, true, as_rvec_array(coords.data()), nat, nat_loc, ind_loc, xa_ind, xa_old, box);
1579 /* If connected and master -> remove shifts */
1580 if ((imdstep && bConnected) && MASTER(cr))
1582 removeMolecularShifts(box);
1586 wallcycle_stop(wcycle, WallCycleCounter::Imd);
1591 bool ImdSession::run(int64_t step, bool bNS, const matrix box, gmx::ArrayRef<const gmx::RVec> coords, double t)
1593 return impl_->run(step, bNS, box, coords, t);
1596 void ImdSession::fillEnergyRecord(int64_t step, bool bHaveNewEnergies)
1598 if (!impl_->sessionPossible || !impl_->clientsocket)
1603 IMDEnergyBlock* ene = impl_->energies;
1607 /* In MPI-parallel simulations the energies are not accessible a at every time step.
1608 * We update them if we have new values, otherwise, the energy values from the
1609 * last global communication step are still on display in the viewer. */
1610 if (bHaveNewEnergies)
1612 ene->T_abs = static_cast<float>(impl_->enerd->term[F_TEMP]);
1613 ene->E_pot = static_cast<float>(impl_->enerd->term[F_EPOT]);
1614 ene->E_tot = static_cast<float>(impl_->enerd->term[F_ETOT]);
1615 ene->E_bond = static_cast<float>(impl_->enerd->term[F_BONDS]);
1616 ene->E_angle = static_cast<float>(impl_->enerd->term[F_ANGLES]);
1617 ene->E_dihe = static_cast<float>(impl_->enerd->term[F_PDIHS]);
1618 ene->E_impr = static_cast<float>(impl_->enerd->term[F_IDIHS]);
1619 ene->E_vdw = static_cast<float>(impl_->enerd->term[F_LJ]);
1620 ene->E_coul = static_cast<float>(impl_->enerd->term[F_COUL_SR]);
1625 void ImdSession::sendPositionsAndEnergies()
1627 if (!impl_->sessionPossible || !impl_->clientsocket)
1632 if (imd_send_energies(impl_->clientsocket, impl_->energies, impl_->energysendbuf))
1634 impl_->issueFatalError("Error sending updated energies. Disconnecting client.");
1637 if (imd_send_rvecs(impl_->clientsocket, impl_->nat, impl_->xa, impl_->coordsendbuf))
1639 impl_->issueFatalError("Error sending updated positions. Disconnecting client.");
1644 void ImdSession::updateEnergyRecordAndSendPositionsAndEnergies(bool bIMDstep, int64_t step, bool bHaveNewEnergies)
1646 if (!impl_->sessionPossible)
1651 wallcycle_start(impl_->wcycle, WallCycleCounter::Imd);
1653 /* Update time step for IMD and prepare IMD energy record if we have new energies. */
1654 fillEnergyRecord(step, bHaveNewEnergies);
1658 /* Send positions and energies to VMD client via IMD */
1659 sendPositionsAndEnergies();
1662 wallcycle_stop(impl_->wcycle, WallCycleCounter::Imd);
1665 void ImdSession::applyForces(gmx::ArrayRef<gmx::RVec> force)
1667 if (!impl_->sessionPossible || !impl_->bForceActivated)
1672 wallcycle_start(impl_->wcycle, WallCycleCounter::Imd);
1674 for (int i = 0; i < impl_->nforces; i++)
1676 /* j are the indices in the "System group".*/
1677 int j = impl_->ind[impl_->f_ind[i]];
1679 /* check if this is a local atom and find out locndx */
1681 const t_commrec* cr = impl_->cr;
1682 if (PAR(cr) && (locndx = cr->dd->ga2la->findHome(j)))
1687 rvec_inc(force[j], impl_->f[i]);
1690 wallcycle_stop(impl_->wcycle, WallCycleCounter::Imd);
1693 ImdSession::ImdSession(const MDLogger& mdlog) : impl_(new Impl(mdlog)) {}
1695 ImdSession::~ImdSession() = default;