Make ImdSession into a Pimpl-ed class with factory function
authorMark Abraham <mark.j.abraham@gmail.com>
Sun, 31 Mar 2019 12:29:20 +0000 (14:29 +0200)
committerMark Abraham <mark.j.abraham@gmail.com>
Wed, 10 Apr 2019 05:39:05 +0000 (07:39 +0200)
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

15 files changed:
src/gromacs/domdec/partition.cpp
src/gromacs/domdec/partition.h
src/gromacs/gmxpreprocess/grompp.cpp
src/gromacs/imd/imd.cpp
src/gromacs/imd/imd.h
src/gromacs/imd/imdsocket.cpp
src/gromacs/imd/imdsocket.h
src/gromacs/mdlib/force.h
src/gromacs/mdlib/shellfc.cpp
src/gromacs/mdlib/shellfc.h
src/gromacs/mdlib/sim_util.cpp
src/gromacs/mdrun/integrator.h
src/gromacs/mdrun/md.cpp
src/gromacs/mdrun/minimize.cpp
src/gromacs/mdrun/runner.cpp

index ea168fa467667f7e3287cdf5c845712e54e02ecd..28b3008347154979530f0b80963f26d70ffd6a43 100644 (file)
@@ -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<gmx::RVec> *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);
 
index b7196244f3a89816a9a7d5bcf22172348e78311c..cb006f1d12cb72bb1c6dc63269c1d0d5f9b55ba6 100644 (file)
@@ -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<gmx::RVec> *f,
                          gmx::MDAtoms            *mdatoms,
index 29a466d0dae8b8e0a8b45d69760e36c8ba6e0ba6..9da76a8db29d4ac8b5c6b76a1d11bde0c4c8020e 100644 (file)
@@ -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);
index 53e3cd4c2bc7f5595639babe3755b65637f39bc9..dd766c1b4206e5441cc0f3630235c1813d6af33a 100644 (file)
@@ -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<char *>(forcendx), retsize);
+    int retsize  = sizeof(int32_t) * nforces;
+    int retbytes = imd_read_multiple(socket, reinterpret_cast<char *>(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<char *>(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<int>(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<int>(gmx_get_stop_condition()) == gmx_stop_cond_none))
+    while ((!clientsocket) && (static_cast<int>(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<int>(itype), IMD_NR, eIMDType_names));
-                imd_fatal(IMDsetup, "Terminating connection");
+                GMX_LOG(mdlog.warning).appendTextFormatted(" %s Received unexpected %s.", IMDstr, enum_name(static_cast<int>(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<ImdSession>
+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<ImdSession> 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<float>(enerd->term[F_TEMP   ]);
-                ene->E_pot   = static_cast<float>(enerd->term[F_EPOT   ]);
-                ene->E_tot   = static_cast<float>(enerd->term[F_ETOT   ]);
-                ene->E_bond  = static_cast<float>(enerd->term[F_BONDS  ]);
-                ene->E_angle = static_cast<float>(enerd->term[F_ANGLES ]);
-                ene->E_dihe  = static_cast<float>(enerd->term[F_PDIHS  ]);
-                ene->E_impr  = static_cast<float>(enerd->term[F_IDIHS  ]);
-                ene->E_vdw   = static_cast<float>(enerd->term[F_LJ     ]);
-                ene->E_coul  = static_cast<float>(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<float>(impl_->enerd->term[F_TEMP   ]);
+        ene->E_pot   = static_cast<float>(impl_->enerd->term[F_EPOT   ]);
+        ene->E_tot   = static_cast<float>(impl_->enerd->term[F_ETOT   ]);
+        ene->E_bond  = static_cast<float>(impl_->enerd->term[F_BONDS  ]);
+        ene->E_angle = static_cast<float>(impl_->enerd->term[F_ANGLES ]);
+        ene->E_dihe  = static_cast<float>(impl_->enerd->term[F_PDIHS  ]);
+        ene->E_impr  = static_cast<float>(impl_->enerd->term[F_IDIHS  ]);
+        ene->E_vdw   = static_cast<float>(impl_->enerd->term[F_LJ     ]);
+        ene->E_coul  = static_cast<float>(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
index 65524ec0dc4861e0f94fd5609b0a4bb0fd3aa6f4..bd91747f4f284363b57db20cb6c35504a01a56eb 100644 (file)
@@ -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 <ckutzne@gwdg.de>
  *
 #ifndef GMX_IMD_IMD_H
 #define GMX_IMD_IMD_H
 
+#include <memory>
+
 #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<ImdSession>
+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> impl_;
+
+    public:
+        // Befriend the factory function.
+        friend std::unique_ptr<ImdSession>
+        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
index c5bb1d19eb6556b83782b607674f916357e81ab1..67e6abde81f6786b4b53e6b4b0bb50c0cba8517d 100644 (file)
@@ -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
index e67ad8763957fb29f745c5502673c01a288afa1d..3cf846dbd70888581b356b87ba916c4de2222c37 100644 (file)
@@ -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
index 67d9d474047d7dc9308eb1707f65c5ab511fe725..b6193b4bbff75486fde28f97d2d3927502f53c48 100644 (file)
@@ -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,
index 067eed78c00c211c7b60e794bb7627b0a41ce54f..095a56497766b639eb8f8b78e22f4adc6d3a6051 100644 (file)
@@ -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,
index c9f868048c9125b7378e3731ff2f0ea7b47c307b..1842314e40ffd4838997c58bf2e2b4e0c7ab670e 100644 (file)
@@ -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,
index ba3874c2edb4792a5245f97423fc030183e8936b..3627015193e2ee7194318e7d318af090c93860d0 100644 (file)
@@ -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,
index 529fcbbd357f5d7f4da9d009979f9120759e8e8a..dc8ccd98ff681fc0ab1d899701372e6beff21a8c 100644 (file)
@@ -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.
index 2a7f4e6f2c4d3fbe0e4bdd7e9291f56b3389a12d..307e3182f040b87d889a69147d835697566069cb 100644 (file)
@@ -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);
index 380342400bdf09d8aeb807609bd86312fffe6a39..4447ac014a6dc85841cd6059331ee5245b4cb887 100644 (file)
@@ -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))
     {
index 2fc235969f000cafaa0695a530416778fedbea11..de8ecc19ab1b69af5987d5b367fae30527967c80 100644 (file)
@@ -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);