From dfd44cd0513cf12e0cd666ec48b5bdd20f65e497 Mon Sep 17 00:00:00 2001 From: Mark Abraham Date: Sun, 31 Mar 2019 14:29:20 +0200 Subject: [PATCH] Make ImdSession into a Pimpl-ed class with factory function This prepares to make IMD into a proper module. No functionality changes in this commit. Replaced gmx_bool with bool Used fast returns when IMD is inactive, for better readability of code. Refs #2877 Change-Id: Ibbe8c452f6f480e9a357fe1b87da3ab0ae166317 --- src/gromacs/domdec/partition.cpp | 4 +- src/gromacs/domdec/partition.h | 4 +- src/gromacs/gmxpreprocess/grompp.cpp | 2 +- src/gromacs/imd/imd.cpp | 1196 +++++++++++++------------- src/gromacs/imd/imd.h | 223 ++--- src/gromacs/imd/imdsocket.cpp | 4 + src/gromacs/imd/imdsocket.h | 5 + src/gromacs/mdlib/force.h | 6 +- src/gromacs/mdlib/shellfc.cpp | 2 +- src/gromacs/mdlib/shellfc.h | 4 +- src/gromacs/mdlib/sim_util.cpp | 6 +- src/gromacs/mdrun/integrator.h | 4 +- src/gromacs/mdrun/md.cpp | 10 +- src/gromacs/mdrun/minimize.cpp | 34 +- src/gromacs/mdrun/runner.cpp | 11 +- 15 files changed, 777 insertions(+), 738 deletions(-) diff --git a/src/gromacs/domdec/partition.cpp b/src/gromacs/domdec/partition.cpp index ea168fa467..28b3008347 100644 --- a/src/gromacs/domdec/partition.cpp +++ b/src/gromacs/domdec/partition.cpp @@ -3006,7 +3006,7 @@ void dd_partition_system(FILE *fplog, t_state *state_global, const gmx_mtop_t &top_global, const t_inputrec *ir, - t_gmx_IMD *imdSession, + gmx::ImdSession *imdSession, t_state *state_local, PaddedVector *f, gmx::MDAtoms *mdAtoms, @@ -3633,7 +3633,7 @@ void dd_partition_system(FILE *fplog, } /* Update the local atoms to be communicated via the IMD protocol if bIMD is TRUE. */ - dd_make_local_IMD_atoms(dd, imdSession); + imdSession->dd_make_local_IMD_atoms(dd); add_dd_statistics(dd); diff --git a/src/gromacs/domdec/partition.h b/src/gromacs/domdec/partition.h index b7196244f3..cb006f1d12 100644 --- a/src/gromacs/domdec/partition.h +++ b/src/gromacs/domdec/partition.h @@ -58,7 +58,6 @@ struct gmx_vsite_t; struct gmx_wallcycle; struct t_commrec; struct t_forcerec; -struct t_gmx_IMD; struct t_inputrec; struct t_nrnb; class t_state; @@ -66,6 +65,7 @@ class t_state; namespace gmx { class Constraints; +class ImdSession; class MDAtoms; class MDLogger; } // namespace @@ -96,7 +96,7 @@ void dd_partition_system(FILE *fplog, t_state *state_global, const gmx_mtop_t &top_global, const t_inputrec *ir, - t_gmx_IMD *imdSession, + gmx::ImdSession *imdSession, t_state *state_local, PaddedVector *f, gmx::MDAtoms *mdatoms, diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 29a466d0da..9da76a8db2 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -2430,7 +2430,7 @@ int gmx_grompp(int argc, char *argv[]) write_tpx_state(ftp2fn(efTPR, NFILE, fnm), ir, &state, &sys); /* Output IMD group, if bIMD is TRUE */ - write_IMDgroup_to_file(ir->bIMD, ir, &state, &sys, NFILE, fnm); + gmx::write_IMDgroup_to_file(ir->bIMD, ir, &state, &sys, NFILE, fnm); sfree(opts->define); sfree(opts->include); diff --git a/src/gromacs/imd/imd.cpp b/src/gromacs/imd/imd.cpp index 53e3cd4c2b..dd766c1b42 100644 --- a/src/gromacs/imd/imd.cpp +++ b/src/gromacs/imd/imd.cpp @@ -82,6 +82,9 @@ #include "gromacs/utility/smalloc.h" #include "gromacs/utility/stringutil.h" +namespace gmx +{ + /*! \brief How long shall we wait in seconds until we check for a connection again? */ constexpr int c_loopWait = 1; @@ -132,81 +135,187 @@ typedef struct /*! \internal - * \brief IMD (interactive molecular dynamics) main data structure. - * - * Contains private IMD data + * \brief Implementation type for the IMD session * - * \todo Rename this e.g. ImdSession, and make it model + * \todo Make this class (or one extracted from it) model * IForceProvider. + * \todo Use RAII for files and allocations */ -typedef struct t_gmx_IMD +class ImdSession::Impl { - bool sessionPossible; /**< True if tpr and mdrun input combine to permit IMD sessions */ - FILE *outf; /**< Output file for IMD data, mainly forces. */ - - int nat; /**< Number of atoms that can be pulled via IMD. */ - int nat_loc; /**< Part of the atoms that are local. */ - int *ind; /**< Global indices of the IMD atoms. */ - int *ind_loc; /**< Local indices of the IMD atoms. */ - int nalloc_loc; /**< Allocation size for ind_loc. */ - rvec *xa; /**< Positions for all IMD atoms assembled on - the master node. */ - ivec *xa_shifts; /**< Shifts for all IMD atoms, to make - molecule(s) whole. */ - ivec *xa_eshifts; /**< Extra shifts since last DD step. */ - rvec *xa_old; /**< Old positions for all IMD atoms on master. */ - int *xa_ind; /**< Position of each local atom in the - collective array. */ - - int nstimd; /**< Global IMD frequency, known to all nodes. */ - int nstimd_new; /**< New frequency from IMD client, master only. */ - int nstimd_def; /**< Default IMD frequency when disconnected. */ - - int port; /**< Port to use for network socket. */ - IMDSocket *socket; /**< The IMD socket on the master node. */ - IMDSocket *clientsocket; /**< The IMD socket on the client. */ - int length; /**< Length we got with last header. */ - - gmx_bool bWConnect; /**< Shall we block and wait for connection? */ - gmx_bool bTerminated; /**< Set if MD is terminated. */ - gmx_bool bTerminatable; /**< Set if MD can be terminated. */ - gmx_bool bConnected; /**< Set if connection is present. */ - gmx_bool bNewForces; /**< Set if we received new forces. */ - gmx_bool bForceActivated; /**< Set if pulling from VMD is allowed. */ - - IMDEnergyBlock *energies; /**< Pointer to energies we send back. */ - - int32_t vmd_nforces; /**< Number of VMD forces. */ - int32_t *vmd_f_ind; /**< VMD forces indices. */ - float *vmd_forces; /**< The VMD forces flat in memory. */ - int nforces; /**< Number of actual MD forces; - this gets communicated to the clients. */ - int *f_ind; /**< Force indices. */ - rvec *f; /**< The IMD pulling forces. */ - - char *forcesendbuf; /**< Buffer for force sending. */ - char *coordsendbuf; /**< Buffer for coordinate sending. */ - char *energysendbuf; /**< Send buffer for energies. */ - rvec *sendxbuf; /**< Buffer to make molecules whole before - sending. */ - - t_block mols; /**< Molecules block in IMD group. */ - - /* The next block is used on the master node only to reduce the output - * without sacrificing information. If any of these values changes, - * we need to write output */ - int old_nforces; /**< Old value for nforces. */ - int *old_f_ind; /**< Old values for force indices. */ - rvec *old_forces; /**< Old values for IMD pulling forces. */ - - // TODO make this a const reference when this struct has a constructor - const gmx::MDLogger *mdlog; /**< Logger */ - -} t_gmx_IMD_setup; - + public: + //! Constructor + Impl(const MDLogger &mdlog); + ~Impl(); + + /*! \brief Prepare the socket on the MASTER. */ + void prepareMasterSocket(); + /*! \brief Disconnect the client. */ + void disconnectClient(); + /*! \brief Prints an error message and disconnects the client. + * + * Does not terminate mdrun! + */ + void issueFatalError(const char *msg); + /*! \brief Check whether we got an incoming connection. */ + bool tryConnect(); + /*! \brief Wrap imd_tryconnect in order to make it blocking. + * + * Used when the simulation should wait for an incoming connection. + */ + void blockConnect(); + /*! \brief Make sure that our array holding the forces received via IMD is large enough. */ + void prepareVmdForces(); + /*! \brief Reads forces received via IMD. */ + void readVmdForces(); + /*! \brief Prepares the MD force arrays. */ + void prepareMDForces(); + /*! \brief Copy IMD forces to MD forces. + * + * Do conversion from Cal->Joule and from + * Angstrom -> nm and from a pointer array to arrays to 3*N array. + */ + void copyToMDForces(); + /*! \brief Return true if any of the forces or indices changed. */ + bool bForcesChanged() const; + /*! \brief Update the old_f_ind and old_forces arrays to contain the current values. */ + void keepOldValues(); + /*! \brief Write the applied pull forces to logfile. + * + * Call on master only! + */ + void outputForces(double time); + /*! \brief Synchronize the nodes. */ + void syncNodes(const t_commrec *cr, double t); + /*! \brief Reads header from the client and decides what to do. */ + void readCommand(); + /*! \brief Open IMD output file and write header information. + * + * Call on master only. + */ + void openOutputFile(const char *fn, + int nat_total, + const gmx_output_env_t *oenv, + const ContinuationOptions &continuationOptions); + /*! \brief Creates the molecule start-end position array of molecules in the IMD group. */ + void prepareMoleculesInImdGroup(const gmx_mtop_t *top_global); + /*! \brief Removes shifts of molecules diffused outside of the box. */ + void removeMolecularShifts(const matrix box); + /*! \brief Initialize arrays used to assemble the positions from the other nodes. */ + void prepareForPositionAssembly(const t_commrec *cr, const rvec x[]); + /*! \brief Interact with any connected VMD session */ + bool run(int64_t step, + bool bNS, + const matrix box, + const rvec x[], + double t); + + // TODO rename all the data members to have underscore suffixes + + //! True if tpr and mdrun input combine to permit IMD sessions + bool sessionPossible = false; + //! Output file for IMD data, mainly forces. + FILE *outf = nullptr; + + //! Number of atoms that can be pulled via IMD. + int nat = 0; + //! Part of the atoms that are local. + int nat_loc = 0; + //! Global indices of the IMD atoms. + int *ind = nullptr; + //! Local indices of the IMD atoms. + int *ind_loc = nullptr; + //! Allocation size for ind_loc. + int nalloc_loc = 0; + //! Positions for all IMD atoms assembled on the master node. + rvec *xa = nullptr; + //! Shifts for all IMD atoms, to make molecule(s) whole. + ivec *xa_shifts = nullptr; + //! Extra shifts since last DD step. + ivec *xa_eshifts = nullptr; + //! Old positions for all IMD atoms on master. + rvec *xa_old = nullptr; + //! Position of each local atom in the collective array. + int *xa_ind = nullptr; + + //! Global IMD frequency, known to all ranks. + int nstimd = 1; + //! New frequency from IMD client, master only. + int nstimd_new = 1; + //! Default IMD frequency when disconnected. + int defaultNstImd = -1; + + //! Port to use for network socket. + int port = 0; + //! The IMD socket on the master node. + IMDSocket *socket = nullptr; + //! The IMD socket on the client. + IMDSocket *clientsocket = nullptr; + //! Length we got with last header. + int length = 0; + + //! Shall we block and wait for connection? + bool bWConnect = false; + //! Set if MD is terminated. + bool bTerminated = false; + //! Set if MD can be terminated. + bool bTerminatable = false; + //! Set if connection is present. + bool bConnected = false; + //! Set if we received new forces. + bool bNewForces = false; + //! Set if pulling from VMD is allowed. + bool bForceActivated = false; + + //! Pointer to energies we send back. + IMDEnergyBlock *energies = nullptr; + + //! Number of VMD forces. + int32_t vmd_nforces = 0; + //! VMD forces indices. + int32_t *vmd_f_ind = nullptr; + //! The VMD forces flat in memory. + float *vmd_forces = nullptr; + //! Number of actual MD forces; this gets communicated to the clients. + int nforces = 0; + //! Force indices. + int *f_ind = nullptr; + //! The IMD pulling forces. + rvec *f = nullptr; + + //! Buffer for force sending. + char *forcesendbuf = nullptr; + //! Buffer for coordinate sending. + char *coordsendbuf = nullptr; + //! Send buffer for energies. + char *energysendbuf = nullptr; + //! Buffer to make molecules whole before sending. + rvec *sendxbuf = nullptr; + + //! Molecules block in IMD group. + t_block mols; + + /* The next block is used on the master node only to reduce the output + * without sacrificing information. If any of these values changes, + * we need to write output */ + //! Old value for nforces. + int old_nforces = 0; + //! Old values for force indices. + int *old_f_ind = nullptr; + //! Old values for IMD pulling forces. + rvec *old_forces = nullptr; + + //! Logger + const MDLogger &mdlog; + //! Commmunication object + const t_commrec *cr = nullptr; + //! Wallcycle counting manager. + gmx_wallcycle *wcycle = nullptr; + //! Energy output handler + gmx_enerdata_t *enerd = nullptr; +}; -/*! \internal - * \brief Enum for types of IMD messages. +/*! \brief Enum for types of IMD messages. * * We use the same records as the NAMD/VMD IMD implementation. */ @@ -226,8 +335,7 @@ typedef enum IMDType_t } IMDMessageType; -/*! \internal - * \brief Names of the IMDType for error messages. +/*! \brief Names of the IMDType for error messages. */ static const char *eIMDType_names[IMD_NR + 1] = { "IMD_DISCONNECT", @@ -294,7 +402,7 @@ static int32_t imd_read_multiple(IMDSocket *socket, char *datptr, int32_t toread } leftcount -= countread; datptr += countread; - } /* end while */ + } /* end while */ /* return nr of bytes read */ return toread - leftcount; @@ -323,7 +431,7 @@ static int32_t imd_write_multiple(IMDSocket *socket, const char *datptr, int32_t } leftcount -= countwritten; datptr += countwritten; - } /* end while */ + } /* end while */ return towrite - leftcount; } @@ -377,32 +485,24 @@ static IMDMessageType imd_recv_header(IMDSocket *socket, int32_t *length) * * The number of forces was previously communicated via the header. */ -static int imd_recv_mdcomm(IMDSocket *socket, int32_t nforces, int32_t *forcendx, float *forces) +static bool imd_recv_mdcomm(IMDSocket *socket, int32_t nforces, int32_t *forcendx, float *forces) { - int retsize, retbytes; - - /* reading indices */ - retsize = sizeof(int32_t) * nforces; - retbytes = imd_read_multiple(socket, reinterpret_cast(forcendx), retsize); + int retsize = sizeof(int32_t) * nforces; + int retbytes = imd_read_multiple(socket, reinterpret_cast(forcendx), retsize); if (retbytes != retsize) { - return FALSE; + return false; } /* reading forces as float array */ retsize = 3 * sizeof(float) * nforces; retbytes = imd_read_multiple(socket, reinterpret_cast(forces), retsize); - if (retbytes != retsize) - { - return FALSE; - } - - return TRUE; + return (retbytes == retsize); } /* GROMACS specific functions for the IMD implementation */ -void write_IMDgroup_to_file(gmx_bool bIMD, t_inputrec *ir, const t_state *state, +void write_IMDgroup_to_file(bool bIMD, t_inputrec *ir, const t_state *state, const gmx_mtop_t *sys, int nfile, const t_filenm fnm[]) { t_atoms IMDatoms; @@ -417,18 +517,16 @@ void write_IMDgroup_to_file(gmx_bool bIMD, t_inputrec *ir, const t_state *state, } -void dd_make_local_IMD_atoms(const gmx_domdec_t *dd, t_gmx_IMD *IMDsetup) +void +ImdSession::dd_make_local_IMD_atoms(const gmx_domdec_t *dd) { - const gmx_ga2la_t *ga2la; - - if (IMDsetup->sessionPossible) + if (!impl_->sessionPossible) { - ga2la = dd->ga2la; - - dd_make_local_group_indices( - ga2la, IMDsetup->nat, IMDsetup->ind, &IMDsetup->nat_loc, - &IMDsetup->ind_loc, &IMDsetup->nalloc_loc, IMDsetup->xa_ind); + return; } + + dd_make_local_group_indices(dd->ga2la, impl_->nat, impl_->ind, &impl_->nat_loc, + &impl_->ind_loc, &impl_->nalloc_loc, impl_->xa_ind); } @@ -461,151 +559,112 @@ static int imd_send_rvecs(IMDSocket *socket, int nat, rvec *x, char *buffer) } -/*! \brief Initializes the IMD private data. */ -static void prepareSession(t_gmx_IMD_setup *IMDsetup, - const gmx::MDLogger &mdlog, - int imdatoms, int nstimddef, int imdport) +void +ImdSession::Impl::prepareMasterSocket() { - IMDsetup->sessionPossible = true; - IMDsetup->nat = imdatoms; - IMDsetup->bTerminated = FALSE; - IMDsetup->bTerminatable = FALSE; - IMDsetup->bWConnect = FALSE; - IMDsetup->bConnected = FALSE; - IMDsetup->bForceActivated = FALSE; - IMDsetup->bNewForces = FALSE; - IMDsetup->bForceActivated = FALSE; - IMDsetup->nstimd = 1; - IMDsetup->nstimd_new = 1; - IMDsetup->nstimd_def = nstimddef; - if (imdport < 1) - { - IMDsetup->port = 0; - } - else - { - IMDsetup->port = imdport; - } - IMDsetup->mdlog = &mdlog; -} - - -/*! \brief Prepare the socket on the MASTER. */ -static void imd_prepare_master_socket(t_gmx_IMD_setup *IMDsetup) -{ - int ret; - - if (imdsock_winsockinit() == -1) { gmx_fatal(FARGS, "%s Failed to initialize winsock.\n", IMDstr); } /* The rest is identical, first create and bind a socket and set to listen then. */ - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s Setting up incoming socket.", IMDstr); - IMDsetup->socket = imdsock_create(); - if (!IMDsetup->socket) + GMX_LOG(mdlog.warning).appendTextFormatted("%s Setting up incoming socket.", IMDstr); + socket = imdsock_create(); + if (!socket) { gmx_fatal(FARGS, "%s Failed to create socket.", IMDstr); } /* Bind to port */ - ret = imdsock_bind(IMDsetup->socket, IMDsetup->port); + int ret = imdsock_bind(socket, port); if (ret) { - gmx_fatal(FARGS, "%s binding socket to port %d failed with error %d.\n", IMDstr, IMDsetup->port, ret); + gmx_fatal(FARGS, "%s binding socket to port %d failed with error %d.\n", IMDstr, port, ret); } - if (imd_sock_listen(IMDsetup->socket)) + if (imd_sock_listen(socket)) { gmx_fatal(FARGS, "%s socket listen failed with error %d.\n", IMDstr, ret); } - if (imdsock_getport(IMDsetup->socket, &IMDsetup->port)) + if (imdsock_getport(socket, &port)) { gmx_fatal(FARGS, "%s Could not determine port number.\n", IMDstr); } - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s Listening for IMD connection on port %d.", IMDstr, IMDsetup->port); + GMX_LOG(mdlog.warning).appendTextFormatted("%s Listening for IMD connection on port %d.", IMDstr, port); } -/*! \brief Disconnect the client. */ -static void imd_disconnect(t_gmx_IMD_setup *IMDsetup) +void +ImdSession::Impl::disconnectClient() { /* Write out any buffered pulling data */ - fflush(IMDsetup->outf); + fflush(outf); /* we first try to shut down the clientsocket */ - imdsock_shutdown(IMDsetup->clientsocket); - if (!imdsock_destroy(IMDsetup->clientsocket)) + imdsock_shutdown(clientsocket); + if (!imdsock_destroy(clientsocket)) { - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s Failed to destroy socket.", IMDstr); + GMX_LOG(mdlog.warning).appendTextFormatted("%s Failed to destroy socket.", IMDstr); } /* then we reset the IMD step to its default, and reset the connection boolean */ - IMDsetup->nstimd_new = IMDsetup->nstimd_def; - IMDsetup->clientsocket = nullptr; - IMDsetup->bConnected = FALSE; + nstimd_new = defaultNstImd; + clientsocket = nullptr; + bConnected = false; } -/*! \brief Prints an error message and disconnects the client. - * - * Does not terminate mdrun! - */ -static void imd_fatal(t_gmx_IMD_setup *IMDsetup, - const char *msg) +void +ImdSession::Impl::issueFatalError(const char *msg) { - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s %s", IMDstr, msg); - imd_disconnect(IMDsetup); - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s disconnected.", IMDstr); + GMX_LOG(mdlog.warning).appendTextFormatted("%s %s", IMDstr, msg); + disconnectClient(); + GMX_LOG(mdlog.warning).appendTextFormatted("%s disconnected.", IMDstr); } -/*! \brief Check whether we got an incoming connection. */ -static gmx_bool imd_tryconnect(t_gmx_IMD_setup *IMDsetup) +bool +ImdSession::Impl::tryConnect() { - if (imdsock_tryread(IMDsetup->socket, 0, 0) > 0) + if (imdsock_tryread(socket, 0, 0) > 0) { /* yes, we got something, accept on clientsocket */ - IMDsetup->clientsocket = imdsock_accept(IMDsetup->socket); - if (!IMDsetup->clientsocket) + clientsocket = imdsock_accept(socket); + if (!clientsocket) { - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s Accepting the connection on the socket failed.", IMDstr); - return FALSE; + GMX_LOG(mdlog.warning).appendTextFormatted("%s Accepting the connection on the socket failed.", IMDstr); + return false; } /* handshake with client */ - if (imd_handshake(IMDsetup->clientsocket)) + if (imd_handshake(clientsocket)) { - imd_fatal(IMDsetup, "Connection failed."); - return FALSE; + issueFatalError("Connection failed."); + return false; } - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s Connection established, checking if I got IMD_GO orders.", IMDstr); + GMX_LOG(mdlog.warning).appendTextFormatted("%s Connection established, checking if I got IMD_GO orders.", IMDstr); /* Check if we get the proper "GO" command from client. */ - if (imdsock_tryread(IMDsetup->clientsocket, c_connectWait, 0) != 1 || imd_recv_header(IMDsetup->clientsocket, &(IMDsetup->length)) != IMD_GO) + if (imdsock_tryread(clientsocket, c_connectWait, 0) != 1 || imd_recv_header(clientsocket, &(length)) != IMD_GO) { - imd_fatal(IMDsetup, "No IMD_GO order received. IMD connection failed."); + issueFatalError("No IMD_GO order received. IMD connection failed."); } /* IMD connected */ - IMDsetup->bConnected = TRUE; + bConnected = true; - return TRUE; + return true; } - return FALSE; + return false; } -/*! \brief Wrap imd_tryconnect in order to make it blocking. - * - * Used when the simulation should wait for an incoming connection. - */ -static void imd_blockconnect(t_gmx_IMD_setup *IMDsetup) +void +ImdSession::Impl::blockConnect() { /* do not wait for connection, when e.g. ctrl+c is pressed and we will terminate anyways. */ if (!(static_cast(gmx_get_stop_condition()) == gmx_stop_cond_none)) @@ -613,193 +672,177 @@ static void imd_blockconnect(t_gmx_IMD_setup *IMDsetup) return; } - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s Will wait until I have a connection and IMD_GO orders.", IMDstr); + GMX_LOG(mdlog.warning).appendTextFormatted("%s Will wait until I have a connection and IMD_GO orders.", IMDstr); /* while we have no clientsocket... 2nd part: we should still react on ctrl+c */ - while ((!IMDsetup->clientsocket) && (static_cast(gmx_get_stop_condition()) == gmx_stop_cond_none)) + while ((!clientsocket) && (static_cast(gmx_get_stop_condition()) == gmx_stop_cond_none)) { - imd_tryconnect(IMDsetup); + tryConnect(); imd_sleep(c_loopWait); } } -/*! \brief Make sure that our array holding the forces received via IMD is large enough. */ -static void imd_prepare_vmd_Forces(t_gmx_IMD_setup *IMDsetup) +void +ImdSession::Impl::prepareVmdForces() { - srenew((IMDsetup->vmd_f_ind), IMDsetup->vmd_nforces); - srenew((IMDsetup->vmd_forces), 3*IMDsetup->vmd_nforces); + srenew((vmd_f_ind), vmd_nforces); + srenew((vmd_forces), 3*vmd_nforces); } -/*! \brief Reads forces received via IMD. */ -static void imd_read_vmd_Forces(t_gmx_IMD_setup *IMDsetup) +void +ImdSession::Impl::readVmdForces() { /* the length of the previously received header tells us the nr of forces we will receive */ - IMDsetup->vmd_nforces = IMDsetup->length; + vmd_nforces = length; /* prepare the arrays */ - imd_prepare_vmd_Forces(IMDsetup); + prepareVmdForces(); /* Now we read the forces... */ - if (!(imd_recv_mdcomm(IMDsetup->clientsocket, IMDsetup->vmd_nforces, IMDsetup->vmd_f_ind, IMDsetup->vmd_forces))) + if (!(imd_recv_mdcomm(clientsocket, vmd_nforces, vmd_f_ind, vmd_forces))) { - imd_fatal(IMDsetup, "Error while reading forces from remote. Disconnecting"); + issueFatalError("Error while reading forces from remote. Disconnecting"); } } -/*! \brief Prepares the MD force arrays. */ -static void imd_prepare_MD_Forces(t_gmx_IMD_setup *IMDsetup) +void +ImdSession::Impl::prepareMDForces() { - srenew((IMDsetup->f_ind), IMDsetup->nforces); - srenew((IMDsetup->f ), IMDsetup->nforces); + srenew((f_ind), nforces); + srenew((f ), nforces); } -/*! \brief Copy IMD forces to MD forces. - * - * Do conversion from Cal->Joule and from - * Angstrom -> nm and from a pointer array to arrays to 3*N array. - */ -static void imd_copyto_MD_Forces(t_gmx_IMD_setup *IMDsetup) +void +ImdSession::Impl::copyToMDForces() { int i; real conversion = CAL2JOULE * NM2A; - for (i = 0; i < IMDsetup->nforces; i++) + for (i = 0; i < nforces; i++) { /* Copy the indices, a copy is important because we may update the incoming forces * whenever we receive new forces while the MD forces are only communicated upon * IMD communication.*/ - IMDsetup->f_ind[i] = IMDsetup->vmd_f_ind[i]; + f_ind[i] = vmd_f_ind[i]; /* Convert to rvecs and do a proper unit conversion */ - IMDsetup->f[i][0] = IMDsetup->vmd_forces[3*i ] * conversion; - IMDsetup->f[i][1] = IMDsetup->vmd_forces[3*i + 1] * conversion; - IMDsetup->f[i][2] = IMDsetup->vmd_forces[3*i + 2] * conversion; + f[i][0] = vmd_forces[3*i ] * conversion; + f[i][1] = vmd_forces[3*i + 1] * conversion; + f[i][2] = vmd_forces[3*i + 2] * conversion; } } -/*! \brief Return TRUE if any of the forces or indices changed. */ -static gmx_bool bForcesChanged(const t_gmx_IMD_setup *IMDsetup) +bool +ImdSession::Impl::bForcesChanged() const { - int i; - - /* First, check whether the number of pulled atoms changed */ - if (IMDsetup->nforces != IMDsetup->old_nforces) + if (nforces != old_nforces) { - return TRUE; + return true; } /* Second, check whether any of the involved atoms changed */ - for (i = 0; i < IMDsetup->nforces; i++) + for (int i = 0; i < nforces; i++) { - if (IMDsetup->f_ind[i] != IMDsetup->old_f_ind[i]) + if (f_ind[i] != old_f_ind[i]) { - return TRUE; + return true; } } + // TODO replace this with rvecs_differ /* Third, check whether all forces are the same */ - for (i = 0; i < IMDsetup->nforces; i++) + for (int i = 0; i < nforces; i++) { - if (IMDsetup->f[i][XX] != IMDsetup->old_forces[i][XX]) + if (f[i][XX] != old_forces[i][XX]) { - return TRUE; + return true; } - if (IMDsetup->f[i][YY] != IMDsetup->old_forces[i][YY]) + if (f[i][YY] != old_forces[i][YY]) { - return TRUE; + return true; } - if (IMDsetup->f[i][ZZ] != IMDsetup->old_forces[i][ZZ]) + if (f[i][ZZ] != old_forces[i][ZZ]) { - return TRUE; + return true; } } /* All old and new forces are identical! */ - return FALSE; + return false; } -/*! \brief Fill the old_f_ind and old_forces arrays with the new, old values. */ -static void keep_old_values(t_gmx_IMD_setup *IMDsetup) +void +ImdSession::Impl::keepOldValues() { - int i; + old_nforces = nforces; - - IMDsetup->old_nforces = IMDsetup->nforces; - - for (i = 0; i < IMDsetup->nforces; i++) + for (int i = 0; i < nforces; i++) { - IMDsetup->old_f_ind[i] = IMDsetup->f_ind[i]; - copy_rvec(IMDsetup->f[i], IMDsetup->old_forces[i]); + old_f_ind[i] = f_ind[i]; + copy_rvec(f[i], old_forces[i]); } } -/*! \brief Returns TRUE if any component of the two rvecs differs. */ -static inline gmx_bool rvecs_differ(const rvec v1, const rvec v2) +/*! \brief Returns true if any component of the two rvecs differs. */ +static inline bool rvecs_differ(const rvec v1, const rvec v2) { - int i; - - - for (i = 0; i < DIM; i++) + for (int i = 0; i < DIM; i++) { if (v1[i] != v2[i]) { - return TRUE; + return true; } } - return FALSE; + return false; } -/*! \brief Write the applied pull forces to logfile. - * - * Call on master only! - */ -static void output_imd_forces(t_gmx_IMD_setup *IMDsetup, double time) +void +ImdSession::Impl::outputForces(double time) { - if (bForcesChanged(IMDsetup)) + if (!bForcesChanged()) { - /* Write time and total number of applied IMD forces */ - fprintf(IMDsetup->outf, "%14.6e%6d", time, IMDsetup->nforces); + return; + } + + /* Write time and total number of applied IMD forces */ + fprintf(outf, "%14.6e%6d", time, nforces); - /* Write out the global atom indices of the pulled atoms and the forces itself, - * write out a force only if it has changed since the last output */ - for (int i = 0; i < IMDsetup->nforces; i++) + /* Write out the global atom indices of the pulled atoms and the forces itself, + * write out a force only if it has changed since the last output */ + for (int i = 0; i < nforces; i++) + { + if (rvecs_differ(f[i], old_forces[i])) { - if (rvecs_differ(IMDsetup->f[i], IMDsetup->old_forces[i])) - { - fprintf(IMDsetup->outf, "%9d", IMDsetup->ind[IMDsetup->f_ind[i]] + 1); - fprintf(IMDsetup->outf, "%12.4e%12.4e%12.4e", IMDsetup->f[i][0], IMDsetup->f[i][1], IMDsetup->f[i][2]); - } + fprintf(outf, "%9d", ind[f_ind[i]] + 1); + fprintf(outf, "%12.4e%12.4e%12.4e", f[i][0], f[i][1], f[i][2]); } - fprintf(IMDsetup->outf, "\n"); - - keep_old_values(IMDsetup); } + fprintf(outf, "\n"); + + keepOldValues(); } -/*! \brief Synchronize the nodes. */ -static void imd_sync_nodes(t_gmx_IMD_setup *IMDsetup, const t_commrec *cr, double t) +void +ImdSession::Impl::syncNodes(const t_commrec *cr, double t) { - int new_nforces = 0; - - /* Notify the other nodes whether we are still connected. */ if (PAR(cr)) { - block_bc(cr, IMDsetup->bConnected); + block_bc(cr, bConnected); } /* ...if not connected, the job is done here. */ - if (!IMDsetup->bConnected) + if (!bConnected) { return; } @@ -807,30 +850,31 @@ static void imd_sync_nodes(t_gmx_IMD_setup *IMDsetup, const t_commrec *cr, doubl /* Let the other nodes know whether we got a new IMD synchronization frequency. */ if (PAR(cr)) { - block_bc(cr, IMDsetup->nstimd_new); + block_bc(cr, nstimd_new); } /* Now we all set the (new) nstimd communication time step */ - IMDsetup->nstimd = IMDsetup->nstimd_new; + nstimd = nstimd_new; /* We're done if we don't allow pulling at all */ - if (!(IMDsetup->bForceActivated)) + if (!(bForceActivated)) { return; } /* OK, let's check if we have received forces which we need to communicate * to the other nodes */ + int new_nforces = 0; if (MASTER(cr)) { - if (IMDsetup->bNewForces) + if (bNewForces) { - new_nforces = IMDsetup->vmd_nforces; + new_nforces = vmd_nforces; } /* make the "new_forces" negative, when we did not receive new ones */ else { - new_nforces = IMDsetup->vmd_nforces * -1; + new_nforces = vmd_nforces * -1; } } @@ -842,93 +886,93 @@ static void imd_sync_nodes(t_gmx_IMD_setup *IMDsetup, const t_commrec *cr, doubl /* When new_natoms < 0 then we know that these are still the same forces * so we don't communicate them, otherwise... */ - if (new_nforces >= 0) + if (new_nforces < 0) { - /* set local VMD and nforces */ - IMDsetup->vmd_nforces = new_nforces; - IMDsetup->nforces = IMDsetup->vmd_nforces; + return; + } - /* now everybody knows the number of forces in f_ind, so we can prepare - * the target arrays for indices and forces */ - imd_prepare_MD_Forces(IMDsetup); + /* set local VMD and nforces */ + vmd_nforces = new_nforces; + nforces = vmd_nforces; - /* we first update the MD forces on the master by converting the VMD forces */ - if (MASTER(cr)) - { - imd_copyto_MD_Forces(IMDsetup); - /* We also write out forces on every update, so that we know which - * forces are applied for every step */ - if (IMDsetup->outf) - { - output_imd_forces(IMDsetup, t); - } - } + /* now everybody knows the number of forces in f_ind, so we can prepare + * the target arrays for indices and forces */ + prepareMDForces(); - /* In parallel mode we communicate the to-be-applied forces to the other nodes */ - if (PAR(cr)) + /* we first update the MD forces on the master by converting the VMD forces */ + if (MASTER(cr)) + { + copyToMDForces(); + /* We also write out forces on every update, so that we know which + * forces are applied for every step */ + if (outf) { - nblock_bc(cr, IMDsetup->nforces, IMDsetup->f_ind); - nblock_bc(cr, IMDsetup->nforces, IMDsetup->f ); + outputForces(t); } + } - /* done communicating the forces, reset bNewForces */ - IMDsetup->bNewForces = FALSE; + /* In parallel mode we communicate the to-be-applied forces to the other nodes */ + if (PAR(cr)) + { + nblock_bc(cr, nforces, f_ind); + nblock_bc(cr, nforces, f ); } + + /* done communicating the forces, reset bNewForces */ + bNewForces = false; } -/*! \brief Reads header from the client and decides what to do. */ -static void imd_readcommand(t_gmx_IMD_setup *IMDsetup) +void +ImdSession::Impl::readCommand() { - gmx_bool IMDpaused = FALSE; - IMDMessageType itype; - + bool IMDpaused = false; - while (IMDsetup->clientsocket && (imdsock_tryread(IMDsetup->clientsocket, 0, 0) > 0 || IMDpaused)) + while (clientsocket && (imdsock_tryread(clientsocket, 0, 0) > 0 || IMDpaused)) { - itype = imd_recv_header(IMDsetup->clientsocket, &(IMDsetup->length)); + IMDMessageType itype = imd_recv_header(clientsocket, &(length)); /* let's see what we got: */ switch (itype) { /* IMD asks us to terminate the simulation, check if the user allowed this */ case IMD_KILL: - if (IMDsetup->bTerminatable) + if (bTerminatable) { - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted(" %s Terminating connection and running simulation (if supported by integrator).", IMDstr); - IMDsetup->bTerminated = TRUE; - IMDsetup->bWConnect = FALSE; + GMX_LOG(mdlog.warning).appendTextFormatted(" %s Terminating connection and running simulation (if supported by integrator).", IMDstr); + bTerminated = true; + bWConnect = false; gmx_set_stop_condition(gmx_stop_cond_next); } else { - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted(" %s Set -imdterm command line switch to allow mdrun termination from within IMD.", IMDstr); + GMX_LOG(mdlog.warning).appendTextFormatted(" %s Set -imdterm command line switch to allow mdrun termination from within IMD.", IMDstr); } break; /* the client doen't want to talk to us anymore */ case IMD_DISCONNECT: - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted(" %s Disconnecting client.", IMDstr); - imd_disconnect(IMDsetup); + GMX_LOG(mdlog.warning).appendTextFormatted(" %s Disconnecting client.", IMDstr); + disconnectClient(); break; /* we got new forces, read them and set bNewForces flag */ case IMD_MDCOMM: - imd_read_vmd_Forces(IMDsetup); - IMDsetup->bNewForces = TRUE; + readVmdForces(); + bNewForces = true; break; /* the client asks us to (un)pause the simulation. So we toggle the IMDpaused state */ case IMD_PAUSE: if (IMDpaused) { - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted(" %s Un-pause command received.", IMDstr); - IMDpaused = FALSE; + GMX_LOG(mdlog.warning).appendTextFormatted(" %s Un-pause command received.", IMDstr); + IMDpaused = false; } else { - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted(" %s Pause command received.", IMDstr); - IMDpaused = TRUE; + GMX_LOG(mdlog.warning).appendTextFormatted(" %s Pause command received.", IMDstr); + IMDpaused = true; } break; @@ -936,99 +980,86 @@ static void imd_readcommand(t_gmx_IMD_setup *IMDsetup) /* the client sets a new transfer rate, if we get 0, we reset the rate * to the default. VMD filters 0 however */ case IMD_TRATE: - IMDsetup->nstimd_new = (IMDsetup->length > 0) ? IMDsetup->length : IMDsetup->nstimd_def; - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted(" %s Update frequency will be set to %d.", IMDstr, IMDsetup->nstimd_new); + nstimd_new = (length > 0) ? length : defaultNstImd; + GMX_LOG(mdlog.warning).appendTextFormatted(" %s Update frequency will be set to %d.", IMDstr, nstimd_new); break; /* Catch all rule for the remaining IMD types which we don't expect */ default: - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted(" %s Received unexpected %s.", IMDstr, enum_name(static_cast(itype), IMD_NR, eIMDType_names)); - imd_fatal(IMDsetup, "Terminating connection"); + GMX_LOG(mdlog.warning).appendTextFormatted(" %s Received unexpected %s.", IMDstr, enum_name(static_cast(itype), IMD_NR, eIMDType_names)); + issueFatalError("Terminating connection"); break; } /* end switch */ } /* end while */ } -/*! \brief Open IMD output file and write header information. - * - * Call on master only. - */ -static FILE *open_imd_out(const char *fn, - t_gmx_IMD_setup *IMDsetup, - int nat_total, - const gmx_output_env_t *oenv, - const gmx::ContinuationOptions &continuationOptions) +void +ImdSession::Impl::openOutputFile(const char *fn, + int nat_total, + const gmx_output_env_t *oenv, + const ContinuationOptions &continuationOptions) { - FILE *fp; - - /* Open log file of applied IMD forces if requested */ - if (fn && oenv) + if (!fn || !oenv) { - /* If we append to an existing file, all the header information is already there */ - if (continuationOptions.appendFiles) - { - fp = gmx_fio_fopen(fn, "a+"); - } - else - { - fp = gmx_fio_fopen(fn, "w+"); - if (IMDsetup->nat == nat_total) - { - fprintf(fp, "# Note that you can select an IMD index group in the .mdp file if a subset of the atoms suffices.\n"); - } - - xvgr_header(fp, "IMD Pull Forces", "Time (ps)", "# of Forces / Atom IDs / Forces (kJ/mol)", exvggtNONE, oenv); + fprintf(stdout, "%s For a log of the IMD pull forces explicitly specify '-if' on the command line.\n" + "%s (Not possible with energy minimization.)\n", IMDstr, IMDstr); + return; + } - fprintf(fp, "# Can display and manipulate %d (of a total of %d) atoms via IMD.\n", IMDsetup->nat, nat_total); - fprintf(fp, "# column 1 : time (ps)\n"); - fprintf(fp, "# column 2 : total number of atoms feeling an IMD pulling force at that time\n"); - fprintf(fp, "# cols. 3.-6 : global atom number of pulled atom, x-force, y-force, z-force (kJ/mol)\n"); - fprintf(fp, "# then follow : atom-ID, f[x], f[y], f[z] for more atoms in case the force on multiple atoms is changed simultaneously.\n"); - fprintf(fp, "# Note that the force on any atom is always equal to the last value for that atom-ID found in the data.\n"); - fflush(fp); + /* If we append to an existing file, all the header information is already there */ + if (continuationOptions.appendFiles) + { + outf = gmx_fio_fopen(fn, "a+"); + } + else + { + outf = gmx_fio_fopen(fn, "w+"); + if (nat == nat_total) + { + fprintf(outf, "# Note that you can select an IMD index group in the .mdp file if a subset of the atoms suffices.\n"); } - /* To reduce the output file size we remember the old values and output only - * when something changed */ - snew(IMDsetup->old_f_ind, IMDsetup->nat); /* One can never pull on more atoms */ - snew(IMDsetup->old_forces, IMDsetup->nat); + xvgr_header(outf, "IMD Pull Forces", "Time (ps)", "# of Forces / Atom IDs / Forces (kJ/mol)", exvggtNONE, oenv); - return fp; + fprintf(outf, "# Can display and manipulate %d (of a total of %d) atoms via IMD.\n", nat, nat_total); + fprintf(outf, "# column 1 : time (ps)\n"); + fprintf(outf, "# column 2 : total number of atoms feeling an IMD pulling force at that time\n"); + fprintf(outf, "# cols. 3.-6 : global atom number of pulled atom, x-force, y-force, z-force (kJ/mol)\n"); + 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"); + 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"); + fflush(outf); } - fprintf(stdout, "%s For a log of the IMD pull forces explicitly specify '-if' on the command line.\n" - "%s (Not possible with energy minimization.)\n", IMDstr, IMDstr); - - return nullptr; + /* To reduce the output file size we remember the old values and output only + * when something changed */ + snew(old_f_ind, nat); /* One can never pull on more atoms */ + snew(old_forces, nat); } -void IMD_finalize(t_gmx_IMD *IMDsetup) +ImdSession::Impl::Impl(const MDLogger &mdlog) + : mdlog(mdlog) { - if (IMDsetup->sessionPossible && IMDsetup->outf) + init_block(&mols); +} + +ImdSession::Impl::~Impl() +{ + if (outf) { - gmx_fio_fclose(IMDsetup->outf); + gmx_fio_fclose(outf); } + done_block(&mols); } -/*! \brief Creates the molecule start-end position array of molecules in the IMD group. */ -static void init_imd_prepare_mols_in_imdgroup(t_gmx_IMD_setup *IMDsetup, const gmx_mtop_t *top_global) +void +ImdSession::Impl::prepareMoleculesInImdGroup(const gmx_mtop_t *top_global) { - int i, ii; - t_block lmols; - int nat; - int *ind; - - nat = IMDsetup->nat; - ind = IMDsetup->ind; - - lmols.nr = 0; - /* check whether index is sorted */ - for (i = 0; i < nat-1; i++) + for (int i = 0; i < nat-1; i++) { if (ind[i] > ind[i+1]) { @@ -1036,15 +1067,17 @@ static void init_imd_prepare_mols_in_imdgroup(t_gmx_IMD_setup *IMDsetup, const g } } - gmx::RangePartitioning gmols = gmx_mtop_molecules(*top_global); + RangePartitioning gmols = gmx_mtop_molecules(*top_global); + t_block lmols; + lmols.nr = 0; snew(lmols.index, gmols.numBlocks() + 1); lmols.index[0] = 0; - for (i = 0; i < gmols.numBlocks(); i++) + for (int i = 0; i < gmols.numBlocks(); i++) { auto mol = gmols.block(i); int count = 0; - for (ii = 0; ii < nat; ii++) + for (int ii = 0; ii < nat; ii++) { if (mol.inRange(ind[ii])) { @@ -1059,7 +1092,7 @@ static void init_imd_prepare_mols_in_imdgroup(t_gmx_IMD_setup *IMDsetup, const g } srenew(lmols.index, lmols.nr+1); lmols.nalloc_index = lmols.nr+1; - IMDsetup->mols = lmols; + mols = lmols; } @@ -1102,55 +1135,48 @@ static void shift_positions( } -/*! \brief Removes shifts of molecules diffused outside of the box. */ -static void imd_remove_molshifts(t_gmx_IMD_setup *IMDsetup, const matrix box) +void ImdSession::Impl::removeMolecularShifts(const matrix box) { - int i, ii, molsize; - ivec largest, smallest, shift; - t_block mols; - - - mols = IMDsetup->mols; - /* for each molecule also present in IMD group */ - for (i = 0; i < mols.nr; i++) + for (int i = 0; i < mols.nr; i++) { /* first we determine the minimum and maximum shifts for each molecule */ + ivec largest, smallest, shift; clear_ivec(largest); clear_ivec(smallest); clear_ivec(shift); - copy_ivec(IMDsetup->xa_shifts[mols.index[i]], largest); - copy_ivec(IMDsetup->xa_shifts[mols.index[i]], smallest); + copy_ivec(xa_shifts[mols.index[i]], largest); + copy_ivec(xa_shifts[mols.index[i]], smallest); - for (ii = mols.index[i]+1; ii < mols.index[i+1]; ii++) + for (int ii = mols.index[i]+1; ii < mols.index[i+1]; ii++) { - if (IMDsetup->xa_shifts[ii][XX] > largest[XX]) + if (xa_shifts[ii][XX] > largest[XX]) { - largest[XX] = IMDsetup->xa_shifts[ii][XX]; + largest[XX] = xa_shifts[ii][XX]; } - if (IMDsetup->xa_shifts[ii][XX] < smallest[XX]) + if (xa_shifts[ii][XX] < smallest[XX]) { - smallest[XX] = IMDsetup->xa_shifts[ii][XX]; + smallest[XX] = xa_shifts[ii][XX]; } - if (IMDsetup->xa_shifts[ii][YY] > largest[YY]) + if (xa_shifts[ii][YY] > largest[YY]) { - largest[YY] = IMDsetup->xa_shifts[ii][YY]; + largest[YY] = xa_shifts[ii][YY]; } - if (IMDsetup->xa_shifts[ii][YY] < smallest[YY]) + if (xa_shifts[ii][YY] < smallest[YY]) { - smallest[YY] = IMDsetup->xa_shifts[ii][YY]; + smallest[YY] = xa_shifts[ii][YY]; } - if (IMDsetup->xa_shifts[ii][ZZ] > largest[ZZ]) + if (xa_shifts[ii][ZZ] > largest[ZZ]) { - largest[ZZ] = IMDsetup->xa_shifts[ii][ZZ]; + largest[ZZ] = xa_shifts[ii][ZZ]; } - if (IMDsetup->xa_shifts[ii][ZZ] < smallest[ZZ]) + if (xa_shifts[ii][ZZ] < smallest[ZZ]) { - smallest[ZZ] = IMDsetup->xa_shifts[ii][ZZ]; + smallest[ZZ] = xa_shifts[ii][ZZ]; } } @@ -1186,54 +1212,51 @@ static void imd_remove_molshifts(t_gmx_IMD_setup *IMDsetup, const matrix box) /* is there a shift at all? */ if ((shift[XX]) || (shift[YY]) || (shift[ZZ])) { - molsize = mols.index[i+1]-mols.index[i]; + int molsize = mols.index[i+1]-mols.index[i]; /* shift the positions */ - shift_positions(box, &(IMDsetup->xa[mols.index[i]]), shift, molsize); + shift_positions(box, &(xa[mols.index[i]]), shift, molsize); } } } -/*! \brief Initialize arrays used to assemble the positions from the other nodes. */ -static void init_imd_prepare_for_x_assembly(const t_commrec *cr, const rvec x[], t_gmx_IMD_setup *IMDsetup) +void +ImdSession::Impl::prepareForPositionAssembly(const t_commrec *cr, const rvec x[]) { - int i, ii; - - - snew(IMDsetup->xa, IMDsetup->nat); - snew(IMDsetup->xa_ind, IMDsetup->nat); - snew(IMDsetup->xa_shifts, IMDsetup->nat); - snew(IMDsetup->xa_eshifts, IMDsetup->nat); - snew(IMDsetup->xa_old, IMDsetup->nat); + snew(xa, nat); + snew(xa_ind, nat); + snew(xa_shifts, nat); + snew(xa_eshifts, nat); + snew(xa_old, nat); /* Save the original (whole) set of positions such that later the * molecule can always be made whole again */ if (MASTER(cr)) { - for (i = 0; i < IMDsetup->nat; i++) + for (int i = 0; i < nat; i++) { - ii = IMDsetup->ind[i]; - copy_rvec(x[ii], IMDsetup->xa_old[i]); + int ii = ind[i]; + copy_rvec(x[ii], xa_old[i]); } } if (!PAR(cr)) { - IMDsetup->nat_loc = IMDsetup->nat; - IMDsetup->ind_loc = IMDsetup->ind; + nat_loc = nat; + ind_loc = ind; /* xa_ind[i] needs to be set to i for serial runs */ - for (i = 0; i < IMDsetup->nat; i++) + for (int i = 0; i < nat; i++) { - IMDsetup->xa_ind[i] = i; + xa_ind[i] = i; } } /* Communicate initial coordinates xa_old to all processes */ if (PAR(cr)) { - gmx_bcast(IMDsetup->nat * sizeof(IMDsetup->xa_old[0]), IMDsetup->xa_old, cr); + gmx_bcast(nat * sizeof(xa_old[0]), xa_old, cr); } } @@ -1250,29 +1273,28 @@ static void imd_check_integrator_parallel(const t_inputrec *ir, const t_commrec } } -t_gmx_IMD *init_IMD(const t_inputrec *ir, - const t_commrec *cr, - const gmx_multisim_t *ms, - const gmx_mtop_t *top_global, - const gmx::MDLogger &mdlog, - const rvec x[], - int nfile, - const t_filenm fnm[], - const gmx_output_env_t *oenv, - const gmx::MdrunOptions &mdrunOptions) +std::unique_ptr +makeImdSession(const t_inputrec *ir, + const t_commrec *cr, + gmx_wallcycle *wcycle, + gmx_enerdata_t *enerd, + const gmx_multisim_t *ms, + const gmx_mtop_t *top_global, + const MDLogger &mdlog, + const rvec x[], + int nfile, + const t_filenm fnm[], + const gmx_output_env_t *oenv, + const MdrunOptions &mdrunOptions) { - int i; - int nat_total; - int32_t bufxsize; - - t_gmx_IMD_setup *IMDsetup; - snew(IMDsetup, 1); + std::unique_ptr session(new ImdSession(mdlog)); + auto impl = session->impl_.get(); /* We will allow IMD sessions only if supported by the binary and explicitly enabled in the .tpr file */ if (!GMX_IMD || !ir->bIMD) { - return IMDsetup; + return session; } // TODO As IMD is intended for interactivity, and the .tpr file @@ -1283,24 +1305,23 @@ t_gmx_IMD *init_IMD(const t_inputrec *ir, // stream named like ImdInfo, to separate it from warning and to // send it to both destinations. - int nstImd; if (EI_DYNAMICS(ir->eI)) { - nstImd = ir->nstcalcenergy; + impl->defaultNstImd = ir->nstcalcenergy; } else if (EI_ENERGY_MINIMIZATION(ir->eI)) { - nstImd = 1; + impl->defaultNstImd = 1; } else { GMX_LOG(mdlog.warning).appendTextFormatted("%s Integrator '%s' is not supported for Interactive Molecular Dynamics, running normally instead", IMDstr, ei_names[ir->eI]); - return IMDsetup; + return session; } if (isMultiSim(ms)) { GMX_LOG(mdlog.warning).appendTextFormatted("%s Cannot use IMD for multiple simulations or replica exchange, running normally instead", IMDstr); - return IMDsetup; + return session; } const auto &options = mdrunOptions.imdOptions; @@ -1321,7 +1342,7 @@ t_gmx_IMD *init_IMD(const t_inputrec *ir, GMX_LOG(mdlog.warning).appendTextFormatted("%s None of the -imd switches was used.\n" "%s This run will not accept incoming IMD connections", IMDstr, IMDstr); } - } /* end master only */ + } /* end master only */ /* Let the other nodes know whether we want IMD */ if (PAR(cr)) @@ -1332,7 +1353,7 @@ t_gmx_IMD *init_IMD(const t_inputrec *ir, /*... if not we are done.*/ if (!createSession) { - return IMDsetup; + return session; } @@ -1346,32 +1367,39 @@ t_gmx_IMD *init_IMD(const t_inputrec *ir, ***************************************************** */ - nat_total = top_global->natoms; + int nat_total = top_global->natoms; - /* Initialize IMD session. If we read in a pre-IMD .tpr file, imd->nat + /* Initialize IMD session. If we read in a pre-IMD .tpr file, ir->imd->nat * will be zero. For those cases we transfer _all_ atomic positions */ - prepareSession(IMDsetup, mdlog, ir->imd->nat > 0 ? ir->imd->nat : nat_total, - nstImd, options.port); + impl->sessionPossible = true; + impl->nat = ir->imd->nat > 0 ? ir->imd->nat : nat_total; + if (options.port >= 1) + { + impl->port = options.port; + } + impl->cr = cr; + impl->wcycle = wcycle; + impl->enerd = enerd; /* We might need to open an output file for IMD forces data */ if (MASTER(cr)) { - IMDsetup->outf = open_imd_out(opt2fn("-if", nfile, fnm), IMDsetup, nat_total, oenv, mdrunOptions.continuationOptions); + impl->openOutputFile(opt2fn("-if", nfile, fnm), nat_total, oenv, mdrunOptions.continuationOptions); } /* Make sure that we operate with a valid atom index array for the IMD atoms */ if (ir->imd->nat > 0) { /* Point to the user-supplied array of atom numbers */ - IMDsetup->ind = ir->imd->ind; + impl->ind = ir->imd->ind; } else { /* Make a dummy (ind[i] = i) array of all atoms */ - snew(IMDsetup->ind, nat_total); - for (i = 0; i < nat_total; i++) + snew(impl->ind, nat_total); + for (int i = 0; i < nat_total; i++) { - IMDsetup->ind[i] = i; + impl->ind[i] = i; } } @@ -1380,89 +1408,83 @@ t_gmx_IMD *init_IMD(const t_inputrec *ir, { /* we allocate memory for our IMD energy structure */ int32_t recsize = c_headerSize + sizeof(IMDEnergyBlock); - snew(IMDsetup->energysendbuf, recsize); + snew(impl->energysendbuf, recsize); /* Shall we wait for a connection? */ if (options.wait) { - IMDsetup->bWConnect = TRUE; - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s Pausing simulation while no IMD connection present (-imdwait).", IMDstr); + impl->bWConnect = true; + GMX_LOG(mdlog.warning).appendTextFormatted("%s Pausing simulation while no IMD connection present (-imdwait).", IMDstr); } /* Will the IMD clients be able to terminate the simulation? */ if (options.terminatable) { - IMDsetup->bTerminatable = TRUE; - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s Allow termination of the simulation from IMD client (-imdterm).", IMDstr); + impl->bTerminatable = true; + GMX_LOG(mdlog.warning).appendTextFormatted("%s Allow termination of the simulation from IMD client (-imdterm).", IMDstr); } /* Is pulling from IMD client allowed? */ if (options.pull) { - IMDsetup->bForceActivated = TRUE; - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s Pulling from IMD remote is enabled (-imdpull).", IMDstr); + impl->bForceActivated = true; + GMX_LOG(mdlog.warning).appendTextFormatted("%s Pulling from IMD remote is enabled (-imdpull).", IMDstr); } /* Initialize send buffers with constant size */ - snew(IMDsetup->sendxbuf, IMDsetup->nat); - snew(IMDsetup->energies, 1); - bufxsize = c_headerSize + 3 * sizeof(float) * IMDsetup->nat; - snew(IMDsetup->coordsendbuf, bufxsize); + snew(impl->sendxbuf, impl->nat); + snew(impl->energies, 1); + int32_t bufxsize = c_headerSize + 3 * sizeof(float) * impl->nat; + snew(impl->coordsendbuf, bufxsize); } /* do we allow interactive pulling? If so let the other nodes know. */ if (PAR(cr)) { - block_bc(cr, IMDsetup->bForceActivated); + block_bc(cr, impl->bForceActivated); } /* setup the listening socket on master process */ if (MASTER(cr)) { - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s Setting port for connection requests to %d.", IMDstr, IMDsetup->port); - imd_prepare_master_socket(IMDsetup); + GMX_LOG(mdlog.warning).appendTextFormatted("%s Setting port for connection requests to %d.", IMDstr, impl->port); + impl->prepareMasterSocket(); /* Wait until we have a connection if specified before */ - if (IMDsetup->bWConnect) + if (impl->bWConnect) { - imd_blockconnect(IMDsetup); + impl->blockConnect(); } else { - GMX_LOG(IMDsetup->mdlog->warning).appendTextFormatted("%s -imdwait not set, starting simulation.", IMDstr); + GMX_LOG(mdlog.warning).appendTextFormatted("%s -imdwait not set, starting simulation.", IMDstr); } } /* Let the other nodes know whether we are connected */ - imd_sync_nodes(IMDsetup, cr, 0); + impl->syncNodes(cr, 0); /* Initialize arrays used to assemble the positions from the other nodes */ - init_imd_prepare_for_x_assembly(cr, x, IMDsetup); + impl->prepareForPositionAssembly(cr, x); /* Initialize molecule blocks to make them whole later...*/ if (MASTER(cr)) { - init_imd_prepare_mols_in_imdgroup(IMDsetup, top_global); + impl->prepareMoleculesInImdGroup(top_global); } - return IMDsetup; + return session; } -gmx_bool do_IMD(t_gmx_IMD *IMDsetup, - int64_t step, - const t_commrec *cr, - gmx_bool bNS, - const matrix box, - const rvec x[], - double t, - gmx_wallcycle *wcycle) +bool ImdSession::Impl::run(int64_t step, + bool bNS, + const matrix box, + const rvec x[], + double t) { - gmx_bool imdstep = FALSE; - - /* IMD at all? */ - if (!IMDsetup->sessionPossible) + if (!sessionPossible) { - return FALSE; + return false; } wallcycle_start(wcycle, ewcIMD); @@ -1471,50 +1493,50 @@ gmx_bool do_IMD(t_gmx_IMD *IMDsetup, if (MASTER(cr)) { /* If not already connected, check for new connections */ - if (!IMDsetup->clientsocket) + if (!clientsocket) { - if (IMDsetup->bWConnect) + if (bWConnect) { - imd_blockconnect(IMDsetup); + blockConnect(); } else { - imd_tryconnect(IMDsetup); + tryConnect(); } } /* Let's see if we have new IMD messages for us */ - if (IMDsetup->clientsocket) + if (clientsocket) { - imd_readcommand(IMDsetup); + readCommand(); } } /* is this an IMD communication step? */ - imdstep = do_per_step(step, IMDsetup->nstimd); + bool imdstep = do_per_step(step, nstimd); /* OK so this is an IMD step ... */ if (imdstep) { /* First we sync all nodes to let everybody know whether we are connected to VMD */ - imd_sync_nodes(IMDsetup, cr, t); + syncNodes(cr, t); } /* If a client is connected, we collect the positions * and put molecules back into the box before transfer */ - if ((imdstep && IMDsetup->bConnected) + if ((imdstep && bConnected) || bNS) /* independent of imdstep, we communicate positions at each NS step */ { /* Transfer the IMD positions to the master node. Every node contributes * its local positions x and stores them in the assembled xa array. */ - communicate_group_positions(cr, IMDsetup->xa, IMDsetup->xa_shifts, IMDsetup->xa_eshifts, - TRUE, x, IMDsetup->nat, IMDsetup->nat_loc, - IMDsetup->ind_loc, IMDsetup->xa_ind, IMDsetup->xa_old, box); + communicate_group_positions(cr, xa, xa_shifts, xa_eshifts, + true, x, nat, nat_loc, + ind_loc, xa_ind, xa_old, box); /* If connected and master -> remove shifts */ - if ((imdstep && IMDsetup->bConnected) && MASTER(cr)) + if ((imdstep && bConnected) && MASTER(cr)) { - imd_remove_molshifts(IMDsetup, box); + removeMolecularShifts(box); } } @@ -1523,107 +1545,121 @@ gmx_bool do_IMD(t_gmx_IMD *IMDsetup, return imdstep; } +bool ImdSession::run(int64_t step, + bool bNS, + const matrix box, + const rvec x[], + double t) +{ + return impl_->run(step, bNS, box, x, t); +} -void IMD_fill_energy_record(t_gmx_IMD *IMDsetup, const gmx_enerdata_t *enerd, - int64_t step, gmx_bool bHaveNewEnergies) +void ImdSession::fillEnergyRecord(int64_t step, + bool bHaveNewEnergies) { - if (IMDsetup->sessionPossible) + if (!impl_->sessionPossible || !impl_->clientsocket) { - if (IMDsetup->clientsocket) - { - IMDEnergyBlock *ene = IMDsetup->energies; + return; + } - ene->tstep = step; + IMDEnergyBlock *ene = impl_->energies; - /* In MPI-parallel simulations the energies are not accessible a at every time step. - * We update them if we have new values, otherwise, the energy values from the - * last global communication step are still on display in the viewer. */ - if (bHaveNewEnergies) - { - ene->T_abs = static_cast(enerd->term[F_TEMP ]); - ene->E_pot = static_cast(enerd->term[F_EPOT ]); - ene->E_tot = static_cast(enerd->term[F_ETOT ]); - ene->E_bond = static_cast(enerd->term[F_BONDS ]); - ene->E_angle = static_cast(enerd->term[F_ANGLES ]); - ene->E_dihe = static_cast(enerd->term[F_PDIHS ]); - ene->E_impr = static_cast(enerd->term[F_IDIHS ]); - ene->E_vdw = static_cast(enerd->term[F_LJ ]); - ene->E_coul = static_cast(enerd->term[F_COUL_SR]); - } - } + ene->tstep = step; + + /* In MPI-parallel simulations the energies are not accessible a at every time step. + * We update them if we have new values, otherwise, the energy values from the + * last global communication step are still on display in the viewer. */ + if (bHaveNewEnergies) + { + ene->T_abs = static_cast(impl_->enerd->term[F_TEMP ]); + ene->E_pot = static_cast(impl_->enerd->term[F_EPOT ]); + ene->E_tot = static_cast(impl_->enerd->term[F_ETOT ]); + ene->E_bond = static_cast(impl_->enerd->term[F_BONDS ]); + ene->E_angle = static_cast(impl_->enerd->term[F_ANGLES ]); + ene->E_dihe = static_cast(impl_->enerd->term[F_PDIHS ]); + ene->E_impr = static_cast(impl_->enerd->term[F_IDIHS ]); + ene->E_vdw = static_cast(impl_->enerd->term[F_LJ ]); + ene->E_coul = static_cast(impl_->enerd->term[F_COUL_SR]); } } -void IMD_send_positions(t_gmx_IMD *IMDsetup) +void +ImdSession::sendPositionsAndEnergies() { - if (IMDsetup->sessionPossible && IMDsetup->clientsocket) + if (!impl_->sessionPossible || !impl_->clientsocket) { + return; + } - if (imd_send_energies(IMDsetup->clientsocket, IMDsetup->energies, IMDsetup->energysendbuf)) - { - imd_fatal(IMDsetup, "Error sending updated energies. Disconnecting client."); - } + if (imd_send_energies(impl_->clientsocket, impl_->energies, impl_->energysendbuf)) + { + impl_->issueFatalError("Error sending updated energies. Disconnecting client."); + } - if (imd_send_rvecs(IMDsetup->clientsocket, IMDsetup->nat, IMDsetup->xa, IMDsetup->coordsendbuf)) - { - imd_fatal(IMDsetup, "Error sending updated positions. Disconnecting client."); - } + if (imd_send_rvecs(impl_->clientsocket, impl_->nat, impl_->xa, impl_->coordsendbuf)) + { + impl_->issueFatalError("Error sending updated positions. Disconnecting client."); } } -void IMD_prep_energies_send_positions(t_gmx_IMD *IMDsetup, - gmx_bool bIMDstep, - const gmx_enerdata_t *enerd, - int64_t step, gmx_bool bHaveNewEnergies, - gmx_wallcycle *wcycle) +void +ImdSession::updateEnergyRecordAndSendPositionsAndEnergies(bool bIMDstep, + int64_t step, + bool bHaveNewEnergies) { - if (IMDsetup->sessionPossible) + if (!impl_->sessionPossible) { - wallcycle_start(wcycle, ewcIMD); + return; + } - /* Update time step for IMD and prepare IMD energy record if we have new energies. */ - IMD_fill_energy_record(IMDsetup, enerd, step, bHaveNewEnergies); + wallcycle_start(impl_->wcycle, ewcIMD); - if (bIMDstep) - { - /* Send positions and energies to VMD client via IMD */ - IMD_send_positions(IMDsetup); - } + /* Update time step for IMD and prepare IMD energy record if we have new energies. */ + fillEnergyRecord(step, bHaveNewEnergies); - wallcycle_stop(wcycle, ewcIMD); + if (bIMDstep) + { + /* Send positions and energies to VMD client via IMD */ + sendPositionsAndEnergies(); } + + wallcycle_stop(impl_->wcycle, ewcIMD); } -void IMD_apply_forces(t_gmx_IMD *IMDsetup, const t_commrec *cr, rvec *f, - gmx_wallcycle *wcycle) +void ImdSession::applyForces(rvec *f) { - if (IMDsetup->sessionPossible) + if (!impl_->sessionPossible || !impl_->bForceActivated) { - wallcycle_start(wcycle, ewcIMD); - - /* Are forces allowed at all? If not we're done */ - if (!IMDsetup->bForceActivated) - { - return; - } + return; + } - for (int i = 0; i < IMDsetup->nforces; i++) - { - /* j are the indices in the "System group".*/ - int j = IMDsetup->ind[IMDsetup->f_ind[i]]; + wallcycle_start(impl_->wcycle, ewcIMD); - /* check if this is a local atom and find out locndx */ - const int *locndx; - if (PAR(cr) && (locndx = cr->dd->ga2la->findHome(j))) - { - j = *locndx; - } + for (int i = 0; i < impl_->nforces; i++) + { + /* j are the indices in the "System group".*/ + int j = impl_->ind[impl_->f_ind[i]]; - rvec_inc(f[j], IMDsetup->f[i]); + /* check if this is a local atom and find out locndx */ + const int *locndx; + const t_commrec *cr = impl_->cr; + if (PAR(cr) && (locndx = cr->dd->ga2la->findHome(j))) + { + j = *locndx; } - wallcycle_start(wcycle, ewcIMD); + rvec_inc(f[j], impl_->f[i]); } + + wallcycle_stop(impl_->wcycle, ewcIMD); } + +ImdSession::ImdSession(const MDLogger &mdlog) + : impl_(new Impl(mdlog)) +{} + +ImdSession::~ImdSession() = default; + +} // namespace gmx diff --git a/src/gromacs/imd/imd.h b/src/gromacs/imd/imd.h index 65524ec0dc..bd91747f4f 100644 --- a/src/gromacs/imd/imd.h +++ b/src/gromacs/imd/imd.h @@ -50,9 +50,9 @@ /*! \libinternal \file * - * \brief - * This file contains datatypes and function declarations necessary for mdrun - * to interface with VMD via the Interactive Molecular Dynamics protocol. + * \brief This file declares the class that coordinates with VMD via + * the Interactive Molecular Dynamics protocol, along with supporting + * free functions. * * \author Martin Hoefling, Carsten Kutzner * @@ -63,8 +63,11 @@ #ifndef GMX_IMD_IMD_H #define GMX_IMD_IMD_H +#include + #include "gromacs/math/vectypes.h" #include "gromacs/utility/basedefinitions.h" +#include "gromacs/utility/classhelpers.h" struct gmx_domdec_t; struct gmx_enerdata_t; @@ -74,16 +77,15 @@ struct gmx_output_env_t; struct gmx_wallcycle; struct t_commrec; struct t_filenm; -struct t_gmx_IMD; struct t_IMD; struct t_inputrec; class t_state; namespace gmx { +class ImdSession; class MDLogger; struct MdrunOptions; -} static const char IMDstr[] = "IMD:"; /**< Tag output from the IMD module with this string. */ @@ -100,31 +102,18 @@ static const char IMDstr[] = "IMD:"; /**< Tag output from the IMD module with t * \param nfile Number of files. * \param fnm Filename struct. */ -void write_IMDgroup_to_file(gmx_bool bIMD, t_inputrec *ir, const t_state *state, +void write_IMDgroup_to_file(bool bIMD, t_inputrec *ir, const t_state *state, const gmx_mtop_t *sys, int nfile, const t_filenm fnm[]); -/*! \brief Make a selection of the home atoms for the IMD group. - * - * Should be called at every domain decomposition. - * Each node checks which of the atoms from "ind" are local and puts its local - * atom numbers into the "ind_local" array. Furthermore, in "xa_ind" it is - * stored at which position each local atom belongs in the assembled/collective - * array, so that on the master node all positions can be merged into the - * assembled array correctly. - * - * \param dd Structure containing domain decomposition data. - * \param imdSession The IMD session - */ -void dd_make_local_IMD_atoms(const gmx_domdec_t *dd, t_gmx_IMD *imdSession); - - -/*! \brief Returns an initialized IMD session. +/*! \brief Makes and returns an initialized IMD session, which may be inactive. * * This function is called before the main MD loop over time steps. * * \param ir The inputrec structure containing the MD input parameters * \param cr Information structure for MPI communication. + * \param wcycle Count wallcycles of IMD routines for diagnostic output. + * \param enerd Contains the GROMACS energies for the different interaction types. * \param ms Handler for multi-simulations. * \param top_global The topology of the whole system. * \param mdlog Logger @@ -133,93 +122,109 @@ void dd_make_local_IMD_atoms(const gmx_domdec_t *dd, t_gmx_IMD *imdSession); * \param fnm Struct containing file names etc. * \param oenv Output options. * \param mdrunOptions Options for mdrun. - * - * \returns A pointer to an initialized IMD session. - */ -t_gmx_IMD *init_IMD(const t_inputrec *ir, - const t_commrec *cr, - const gmx_multisim_t *ms, - const gmx_mtop_t *top_global, - const gmx::MDLogger &mdlog, - const rvec x[], - int nfile, const t_filenm fnm[], const gmx_output_env_t *oenv, - const gmx::MdrunOptions &mdrunOptions); - - -/*! \brief IMD required in this time step? - * Also checks for new IMD connection and syncs the nodes. - * - * \param IMDsetup The IMD session. - * \param step The time step. - * \param cr Information structure for MPI communication. - * \param bNS Is this a neighbor searching step? - * \param box The simulation box. - * \param x The local atomic positions on this node. - * \param t The time. - * \param wcycle Count wallcycles of IMD routines for diagnostic output. - * - * \returns Whether or not we have to do IMD communication at this step. - */ -gmx_bool do_IMD(t_gmx_IMD *IMDsetup, int64_t step, const t_commrec *cr, - gmx_bool bNS, - const matrix box, const rvec x[], double t, - gmx_wallcycle *wcycle); - - -/*! \brief Add external forces from a running interactive molecular dynamics session. - * - * \param IMDsetup The IMD session. - * \param cr Information structure for MPI communication. - * \param f The forces. - * \param wcycle Count wallcycles of IMD routines for diagnostic output. - */ -void IMD_apply_forces(t_gmx_IMD *IMDsetup, - const t_commrec *cr, rvec *f, - gmx_wallcycle *wcycle); - - -/*! \brief Copy energies and convert to float from enerdata to the IMD energy record. - * - * We do no conversion, so units in client are SI! - * - * \param IMDsetup The IMD session. - * \param enerd Contains the GROMACS energies for the different interaction types. - * \param step The time step. - * \param bHaveNewEnergies Only copy energies if we have done global summing of them before. - * - */ -void IMD_fill_energy_record(t_gmx_IMD *IMDsetup, const gmx_enerdata_t *enerd, - int64_t step, gmx_bool bHaveNewEnergies); - - -/*! \brief Send positions and energies to the client. - * - * \param IMDsetup The IMD session. */ -void IMD_send_positions(t_gmx_IMD *IMDsetup); - - -/*! \brief Calls IMD_prepare_energies() and then IMD_send_positions(). - * - * \param IMDsetup The IMD session. - * \param bIMDstep If true, transfer the positions. Otherwise just update the time step and potentially the energy record. - * \param enerd Contains the GROMACS energies for the different interaction types. - * \param step The time step. - * \param bHaveNewEnergies Only update the energy record if we have done global summing of the energies. - * \param wcycle Count wallcycles of IMD routines for diagnostic output. - * - */ -void IMD_prep_energies_send_positions(t_gmx_IMD *IMDsetup, gmx_bool bIMDstep, - const gmx_enerdata_t *enerd, - int64_t step, gmx_bool bHaveNewEnergies, - gmx_wallcycle *wcycle); - -/*! \brief Finalize IMD and do some cleaning up. - * - * Currently, IMD finalize closes the force output file. - * - * \param IMDsetup The IMD session. - */ -void IMD_finalize(t_gmx_IMD *IMDsetup); +std::unique_ptr +makeImdSession(const t_inputrec *ir, + const t_commrec *cr, + gmx_wallcycle *wcycle, + gmx_enerdata_t *enerd, + const gmx_multisim_t *ms, + const gmx_mtop_t *top_global, + const MDLogger &mdlog, + const rvec x[], + int nfile, const t_filenm fnm[], const gmx_output_env_t *oenv, + const MdrunOptions &mdrunOptions); + +class ImdSession +{ + private: + //! Private constructor, to force the use of makeImdSession() + ImdSession(const MDLogger &mdlog); + public: + ~ImdSession(); + + /*! \brief Make a selection of the home atoms for the IMD group. + * + * Should be called at every domain decomposition. + * Each node checks which of the atoms from "ind" are local and puts its local + * atom numbers into the "ind_local" array. Furthermore, in "xa_ind" it is + * stored at which position each local atom belongs in the assembled/collective + * array, so that on the master node all positions can be merged into the + * assembled array correctly. + * + * \param dd Structure containing domain decomposition data. + */ + void dd_make_local_IMD_atoms(const gmx_domdec_t *dd); + + /*! \brief Prepare to do IMD if required in this time step + * Also checks for new IMD connection and syncs the nodes. + * + * \param step The time step. + * \param bNS Is this a neighbor searching step? + * \param box The simulation box. + * \param x The local atomic positions on this node. + * \param t The time. + * + * \returns Whether or not we have to do IMD communication at this step. + */ + bool run(int64_t step, + bool bNS, + const matrix box, + const rvec x[], + double t); + + /*! \brief Add external forces from a running interactive molecular dynamics session. + * + * \param f The forces. + */ + void applyForces(rvec *f); + + /*! \brief Copy energies and convert to float from enerdata to the IMD energy record. + * + * We do no conversion, so units in client are SI! + * + * \param step The time step. + * \param bHaveNewEnergies Only copy energies if we have done global summing of them before. + */ + void fillEnergyRecord(int64_t step, + bool bHaveNewEnergies); + + /*! \brief Send positions and energies to the client. */ + void sendPositionsAndEnergies(); + + /*! \brief Updates the energy record sent to VMD if needed, and sends positions and energies. + * + * \param bIMDstep If true, transfer the positions. Otherwise just update the time step and potentially the energy record. + * \param step The time step. + * \param bHaveNewEnergies Update the energy record if we have done global summing of the energies. + */ + void updateEnergyRecordAndSendPositionsAndEnergies(bool bIMDstep, + int64_t step, + bool bHaveNewEnergies); + + private: + //! Implementation type. + class Impl; + //! Implementation object. + PrivateImplPointer impl_; + + public: + // Befriend the factory function. + friend std::unique_ptr + makeImdSession(const t_inputrec *ir, + const t_commrec *cr, + gmx_wallcycle *wcycle, + gmx_enerdata_t *enerd, + const gmx_multisim_t *ms, + const gmx_mtop_t *top_global, + const MDLogger &mdlog, + const rvec x[], + int nfile, + const t_filenm fnm[], + const gmx_output_env_t *oenv, + const MdrunOptions &mdrunOptions); +}; + +} // namespace gmx #endif diff --git a/src/gromacs/imd/imdsocket.cpp b/src/gromacs/imd/imdsocket.cpp index c5bb1d19eb..67e6abde81 100644 --- a/src/gromacs/imd/imdsocket.cpp +++ b/src/gromacs/imd/imdsocket.cpp @@ -84,6 +84,8 @@ typedef int socklen_t; #endif +namespace gmx +{ /*! \internal * @@ -437,3 +439,5 @@ int imdsock_tryread(IMDSocket *sock, int timeoutsec, int timeoutusec) return ret; } + +} // namespace gmx diff --git a/src/gromacs/imd/imdsocket.h b/src/gromacs/imd/imdsocket.h index e67ad87639..3cf846dbd7 100644 --- a/src/gromacs/imd/imdsocket.h +++ b/src/gromacs/imd/imdsocket.h @@ -51,6 +51,9 @@ #ifndef GMX_IMD_IMDSOCKET_H #define GMX_IMD_IMDSOCKET_H +namespace gmx +{ + struct IMDSocket; @@ -181,4 +184,6 @@ int imdsock_destroy(IMDSocket *sock); */ int imdsock_tryread(IMDSocket *sock, int timeoutsec, int timeoutusec); +} // namespace gmx + #endif diff --git a/src/gromacs/mdlib/force.h b/src/gromacs/mdlib/force.h index 67d9d47404..b6193b4bbf 100644 --- a/src/gromacs/mdlib/force.h +++ b/src/gromacs/mdlib/force.h @@ -56,7 +56,6 @@ struct t_blocka; struct t_commrec; struct t_fcdata; struct t_forcerec; -struct t_gmx_IMD; struct t_graph; struct t_idef; struct t_inputrec; @@ -67,8 +66,9 @@ struct t_nrnb; namespace gmx { class Awh; -class PpForceWorkload; class ForceWithVirial; +class ImdSession; +class PpForceWorkload; class MDLogger; } @@ -96,7 +96,7 @@ void do_force(FILE *log, const t_inputrec *inputrec, gmx::Awh *awh, gmx_enfrot *enforcedRotation, - t_gmx_IMD *imdSession, + gmx::ImdSession *imdSession, int64_t step, t_nrnb *nrnb, gmx_wallcycle *wcycle, diff --git a/src/gromacs/mdlib/shellfc.cpp b/src/gromacs/mdlib/shellfc.cpp index 067eed78c0..095a564977 100644 --- a/src/gromacs/mdlib/shellfc.cpp +++ b/src/gromacs/mdlib/shellfc.cpp @@ -971,7 +971,7 @@ void relax_shell_flexcon(FILE *fplog, gmx_enfrot *enforcedRotation, int64_t mdstep, const t_inputrec *inputrec, - t_gmx_IMD *imdSession, + gmx::ImdSession *imdSession, gmx_bool bDoNS, int force_flags, gmx_localtop_t *top, diff --git a/src/gromacs/mdlib/shellfc.h b/src/gromacs/mdlib/shellfc.h index c9f868048c..1842314e40 100644 --- a/src/gromacs/mdlib/shellfc.h +++ b/src/gromacs/mdlib/shellfc.h @@ -51,7 +51,6 @@ struct gmx_shellfc_t; struct gmx_mtop_t; struct t_forcerec; struct t_fcdata; -struct t_gmx_IMD; struct t_graph; struct t_inputrec; class t_state; @@ -59,6 +58,7 @@ class t_state; namespace gmx { class Constraints; +class ImdSession; class PpForceWorkload; } @@ -82,7 +82,7 @@ void relax_shell_flexcon(FILE *log, gmx_enfrot *enforcedRotation, int64_t mdstep, const t_inputrec *inputrec, - t_gmx_IMD *imdSession, + gmx::ImdSession *imdSession, gmx_bool bDoNS, int force_flags, gmx_localtop_t *top, diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index ba3874c2ed..3627015193 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -519,7 +519,7 @@ computeSpecialForces(FILE *fplog, const t_inputrec *inputrec, gmx::Awh *awh, gmx_enfrot *enforcedRotation, - t_gmx_IMD *imdSession, + gmx::ImdSession *imdSession, int64_t step, double t, gmx_wallcycle_t wcycle, @@ -587,7 +587,7 @@ computeSpecialForces(FILE *fplog, /* Add forces from interactive molecular dynamics (IMD), if any */ if (inputrec->bIMD && computeForces) { - IMD_apply_forces(imdSession, cr, f, wcycle); + imdSession->applyForces(f); } } @@ -805,7 +805,7 @@ void do_force(FILE *fplog, const t_inputrec *inputrec, gmx::Awh *awh, gmx_enfrot *enforcedRotation, - t_gmx_IMD *imdSession, + gmx::ImdSession *imdSession, int64_t step, t_nrnb *nrnb, gmx_wallcycle_t wcycle, diff --git a/src/gromacs/mdrun/integrator.h b/src/gromacs/mdrun/integrator.h index 529fcbbd35..dc8ccd98ff 100644 --- a/src/gromacs/mdrun/integrator.h +++ b/src/gromacs/mdrun/integrator.h @@ -65,7 +65,6 @@ struct t_commrec; struct t_fcdata; struct t_forcerec; struct t_filenm; -struct t_gmx_IMD; struct t_inputrec; struct t_nrnb; class t_state; @@ -77,6 +76,7 @@ class BoxDeformation; class Constraints; class PpForceWorkload; class IMDOutputProvider; +class ImdSession; class MDLogger; class MDAtoms; class StopHandlerBuilder; @@ -141,7 +141,7 @@ struct Integrator //! Contains user input mdp options. t_inputrec *inputrec; //! The Interactive Molecular Dynamics session. - t_gmx_IMD *imdSession; + ImdSession *imdSession; //! Full system topology. gmx_mtop_t *top_global; //! Helper struct for force calculations. diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 2a7f4e6f2c..307e3182f0 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -192,8 +192,7 @@ void gmx::Integrator::do_md() gmx_bool bPMETune = FALSE; gmx_bool bPMETunePrinting = FALSE; - /* Interactive MD */ - gmx_bool bIMDstep = FALSE; + bool bInteractiveMDstep = false; /* Domain decomposition could incorrectly miss a bonded interaction, but checking for that requires a global @@ -1050,7 +1049,7 @@ void gmx::Integrator::do_md() mdrunOptions.writeConfout, bSumEkinhOld); /* Check if IMD step and do IMD communication, if bIMD is TRUE. */ - bIMDstep = do_IMD(imdSession, step, cr, bNS, state->box, state->x.rvec_array(), t, wcycle); + bInteractiveMDstep = imdSession->run(step, bNS, state->box, state->x.rvec_array(), t); /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */ if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))) @@ -1467,7 +1466,7 @@ void gmx::Integrator::do_md() nrnb, fr->pmedata, pme_loadbal, wcycle, walltime_accounting); /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */ - IMD_prep_energies_send_positions(imdSession, bIMDstep, enerd, step, bCalcEner, wcycle); + imdSession->updateEnergyRecordAndSendPositionsAndEnergies(bInteractiveMDstep, step, bCalcEner); } /* End of main MD loop */ @@ -1514,9 +1513,6 @@ void gmx::Integrator::do_md() finish_swapcoords(ir->swap); } - /* IMD cleanup, if bIMD is TRUE. */ - IMD_finalize(imdSession); - walltime_accounting_set_nsteps_done(walltime_accounting, step_rel); global_stat_destroy(gstat); diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp index 380342400b..4447ac014a 100644 --- a/src/gromacs/mdrun/minimize.cpp +++ b/src/gromacs/mdrun/minimize.cpp @@ -353,7 +353,7 @@ static void init_em(FILE *fplog, const char *title, const t_commrec *cr, t_inputrec *ir, - t_gmx_IMD *imdSession, + gmx::ImdSession *imdSession, t_state *state_global, gmx_mtop_t *top_global, em_state_t *ems, gmx_localtop_t *top, t_nrnb *nrnb, @@ -708,7 +708,7 @@ static void em_dd_partition_system(FILE *fplog, const gmx::MDLogger &mdlog, int step, const t_commrec *cr, gmx_mtop_t *top_global, t_inputrec *ir, - t_gmx_IMD *imdSession, + gmx::ImdSession *imdSession, em_state_t *ems, gmx_localtop_t *top, gmx::MDAtoms *mdAtoms, t_forcerec *fr, gmx_vsite_t *vsite, gmx::Constraints *constr, @@ -770,7 +770,7 @@ class EnergyEvaluator //! User input options. t_inputrec *inputrec; //! The Interactive Molecular Dynamics session. - t_gmx_IMD *imdSession; + gmx::ImdSession *imdSession; //! Manages flop accounting. t_nrnb *nrnb; //! Manages wall cycle accounting. @@ -1570,7 +1570,7 @@ Integrator::do_cg() do_log = do_per_step(step, inputrec->nstlog); do_ene = do_per_step(step, inputrec->nstenergy); - IMD_fill_energy_record(imdSession, enerd, step, TRUE); + imdSession->fillEnergyRecord(step, TRUE); if (do_log) { @@ -1582,9 +1582,9 @@ Integrator::do_cg() } /* Send energies and positions to the IMD client if bIMD is TRUE. */ - if (MASTER(cr) && do_IMD(imdSession, step, cr, TRUE, state_global->box, state_global->x.rvec_array(), 0, wcycle)) + if (MASTER(cr) && imdSession->run(step, TRUE, state_global->box, state_global->x.rvec_array(), 0)) { - IMD_send_positions(imdSession); + imdSession->sendPositionsAndEnergies(); } /* Stop when the maximum force lies below tolerance. @@ -1594,8 +1594,6 @@ Integrator::do_cg() } /* End of the loop */ - IMD_finalize(imdSession); - if (converged) { step--; /* we never took that last step in this case */ @@ -2308,7 +2306,7 @@ Integrator::do_lbfgs() do_log = do_per_step(step, inputrec->nstlog); do_ene = do_per_step(step, inputrec->nstenergy); - IMD_fill_energy_record(imdSession, enerd, step, TRUE); + imdSession->fillEnergyRecord(step, TRUE); if (do_log) { @@ -2320,9 +2318,9 @@ Integrator::do_lbfgs() } /* Send x and E to IMD client, if bIMD is TRUE. */ - if (do_IMD(imdSession, step, cr, TRUE, state_global->box, state_global->x.rvec_array(), 0, wcycle) && MASTER(cr)) + if (imdSession->run(step, TRUE, state_global->box, state_global->x.rvec_array(), 0) && MASTER(cr)) { - IMD_send_positions(imdSession); + imdSession->sendPositionsAndEnergies(); } // Reset stepsize in we are doing more iterations @@ -2335,8 +2333,6 @@ Integrator::do_lbfgs() } /* End of the loop */ - IMD_finalize(imdSession); - if (converged) { step--; /* we never took that last step in this case */ @@ -2530,7 +2526,7 @@ Integrator::do_steep() mdatoms->tmass, enerd, nullptr, nullptr, nullptr, nullBox, nullptr, nullptr, vir, pres, nullptr, mu_tot, constr); - IMD_fill_energy_record(imdSession, enerd, count, TRUE); + imdSession->fillEnergyRecord(count, TRUE); const bool do_dr = do_per_step(steps_accepted, inputrec->nstdisreout); const bool do_or = do_per_step(steps_accepted, inputrec->nstorireout); @@ -2604,19 +2600,17 @@ Integrator::do_steep() } /* Send IMD energies and positions, if bIMD is TRUE. */ - if (do_IMD(imdSession, count, cr, TRUE, state_global->box, - MASTER(cr) ? state_global->x.rvec_array() : nullptr, - 0, wcycle) && + if (imdSession->run(count, TRUE, state_global->box, + MASTER(cr) ? state_global->x.rvec_array() : nullptr, + 0) && MASTER(cr)) { - IMD_send_positions(imdSession); + imdSession->sendPositionsAndEnergies(); } count++; } /* End of the loop */ - IMD_finalize(imdSession); - /* Print some data... */ if (MASTER(cr)) { diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp index 2fc235969f..de8ecc19ab 100644 --- a/src/gromacs/mdrun/runner.cpp +++ b/src/gromacs/mdrun/runner.cpp @@ -1500,10 +1500,9 @@ int Mdrunner::mdrunner() init_enerdata(mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].nr, inputrec->fepvals->n_lambda, enerd); /* Set up interactive MD (IMD) */ - t_gmx_IMD *imdSession = - init_IMD(inputrec, cr, ms, &mtop, mdlog, - MASTER(cr) ? globalState->x.rvec_array() : nullptr, - filenames.size(), filenames.data(), oenv, mdrunOptions); + auto imdSession = makeImdSession(inputrec, cr, wcycle, enerd, ms, &mtop, mdlog, + MASTER(cr) ? globalState->x.rvec_array() : nullptr, + filenames.size(), filenames.data(), oenv, mdrunOptions); if (DOMAINDECOMP(cr)) { @@ -1533,7 +1532,7 @@ int Mdrunner::mdrunner() enforcedRotation ? enforcedRotation->getLegacyEnfrot() : nullptr, deform.get(), mdModules_->outputProvider(), - inputrec, imdSession, &mtop, + inputrec, imdSession.get(), &mtop, fcd, globalState.get(), &observablesHistory, @@ -1576,7 +1575,7 @@ int Mdrunner::mdrunner() // clean up cycle counter wallcycle_destroy(wcycle); - // Free PME data +// Free PME data if (pmedata) { gmx_pme_destroy(pmedata); -- 2.22.0