Refactored pull data structures

Moved the pull work variables into a new pull_t struct.
Only the pull input parameters remain in (renamed) pull_params_t.
For the moment struct pull_t pull_work is still present in t_inputrec.
This should later be moved to a more appropriate place.

Change-Id: I78f6f9736d52514f0327c86e43d4fb2ac88d60cd

diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp
index 0b6001e6676d202d83eb81e759ade83500a44c32..c61e2893bc7a38cb3852bd7423fe71c48b9c017b 100644 (file)
--- a/src/gromacs/domdec/domdec.cpp
+++ b/src/gromacs/domdec/domdec.cpp
@@ -9858,7 +9858,7 @@ void dd_partition_system(FILE                *fplog,
     if (ir->bPull)
     {
         /* Update the local pull groups */
-        dd_make_local_pull_groups(dd, ir->pull, mdatoms);
+        dd_make_local_pull_groups(dd, ir->pull_work, mdatoms);
     }
 
     if (ir->bRot)

diff --git a/src/gromacs/fileio/tpxio.c b/src/gromacs/fileio/tpxio.c
index 84efd2070b36912dc0b2ee8c2854030700cffd76..cf77d95e6b27cf4f7624769bd01a7bd3aa31985e 100644 (file)
--- a/src/gromacs/fileio/tpxio.c
+++ b/src/gromacs/fileio/tpxio.c
@@ -632,7 +632,7 @@ static void do_fepvals(t_fileio *fio, t_lambda *fepvals, gmx_bool bRead, int fil
     }
 }
 
-static void do_pull(t_fileio *fio, t_pull *pull, gmx_bool bRead,
+static void do_pull(t_fileio *fio, pull_params_t *pull, gmx_bool bRead,
                     int file_version, int ePullOld)
 {
     int  eGeomOld;

diff --git a/src/gromacs/gmxlib/txtdump.c b/src/gromacs/gmxlib/txtdump.c
index 0732b9dc891d02134ff1442bfbebce86ba3f5117..4d04ff87592bdd4544f93ecfb3031fae0eaa90ca 100644 (file)
--- a/src/gromacs/gmxlib/txtdump.c
+++ b/src/gromacs/gmxlib/txtdump.c
@@ -760,7 +760,7 @@ static void pr_fepvals(FILE *fp, int indent, t_lambda *fep, gmx_bool bMDPformat)
     PS("dhdl-derivatives", DHDLDERIVATIVESTYPE(fep->dhdl_derivatives));
 };
 
-static void pr_pull(FILE *fp, int indent, t_pull *pull)
+static void pr_pull(FILE *fp, int indent, pull_params_t *pull)
 {
     int g;
 

diff --git a/src/gromacs/gmxlib/typedefs.c b/src/gromacs/gmxlib/typedefs.c
index ca19c6c39ca02cef6b8aba0695d3750e036763c1..12c413d79e2ab4c9fad8ed474dc70bfbd1161bf9 100644 (file)
--- a/src/gromacs/gmxlib/typedefs.c
+++ b/src/gromacs/gmxlib/typedefs.c
@@ -80,21 +80,21 @@ static void done_pull_group(t_pull_group *pgrp)
     if (pgrp->nat > 0)
     {
         sfree(pgrp->ind);
-        sfree(pgrp->ind_loc);
         sfree(pgrp->weight);
-        sfree(pgrp->weight_loc);
     }
 }
 
-static void done_pull(t_pull *pull)
+static void done_pull_params(pull_params_t *pull)
 {
     int i;
 
     for (i = 0; i < pull->ngroup+1; i++)
     {
         done_pull_group(pull->group);
-        done_pull_group(pull->dyna);
     }
+
+    sfree(pull->group);
+    sfree(pull->coord);
 }
 
 void done_inputrec(t_inputrec *ir)
@@ -145,7 +145,7 @@ void done_inputrec(t_inputrec *ir)
 
     if (ir->pull)
     {
-        done_pull(ir->pull);
+        done_pull_params(ir->pull);
         sfree(ir->pull);
     }
 }

diff --git a/src/gromacs/gmxpreprocess/readir.c b/src/gromacs/gmxpreprocess/readir.c
index 9e39ad6edea032c6d4e8b4e87b9f97a6fb602bbc..3f3e233de560482151e53e59a97482c7062fec34 100644 (file)
--- a/src/gromacs/gmxpreprocess/readir.c
+++ b/src/gromacs/gmxpreprocess/readir.c
@@ -2718,7 +2718,7 @@ static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
 {
     t_grpopts              *opts;
     gmx_groups_t           *groups;
-    t_pull                 *pull;
+    pull_params_t          *pull;
     int                     natoms, ai, aj, i, j, d, g, imin, jmin;
     t_iatom                *ia;
     int                    *nrdf2, *na_vcm, na_tot;

diff --git a/src/gromacs/gmxpreprocess/readir.h b/src/gromacs/gmxpreprocess/readir.h
index 49b6de6f73b9340d5b1dc352a93f1ae446373f9f..ea76d63004759091f6689f6d17237759bceca28f 100644 (file)
--- a/src/gromacs/gmxpreprocess/readir.h
+++ b/src/gromacs/gmxpreprocess/readir.h
@@ -127,16 +127,16 @@ void do_index(const char* mdparin,
 /* Routines In readpull.c */
 
 char **read_pullparams(int *ninp_p, t_inpfile **inp,
-                       t_pull *pull,
+                       pull_params_t *pull,
                        warninp_t wi);
 /* Reads the pull parameters, returns a list of the pull group names */
 
-void make_pull_groups(t_pull *pull,
+void make_pull_groups(pull_params_t *pull,
                       char **pgnames,
                       const t_blocka *grps, char **gnames);
 /* Process the pull group parameters after reading the index groups */
 
-void make_pull_coords(t_pull *pull);
+void make_pull_coords(pull_params_t *pull);
 /* Process the pull coordinates after reading the pull groups */
 
 void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real lambda,

diff --git a/src/gromacs/gmxpreprocess/readpull.c b/src/gromacs/gmxpreprocess/readpull.c
index 68752443c20234202e2acb25ded4900df37c7f97..49793283ccf7f5a9bcc232805e68cae1cd5d74cc 100644 (file)
--- a/src/gromacs/gmxpreprocess/readpull.c
+++ b/src/gromacs/gmxpreprocess/readpull.c
@@ -186,7 +186,7 @@ static void init_pull_coord(t_pull_coord *pcrd,
 }
 
 char **read_pullparams(int *ninp_p, t_inpfile **inp_p,
-                       t_pull *pull,
+                       pull_params_t *pull,
                        warninp_t wi)
 {
     int           ninp, i, nchar, nscan, m, idum;
@@ -306,7 +306,7 @@ char **read_pullparams(int *ninp_p, t_inpfile **inp_p,
     return grpbuf;
 }
 
-void make_pull_groups(t_pull *pull,
+void make_pull_groups(pull_params_t *pull,
                       char **pgnames,
                       const t_blocka *grps, char **gnames)
 {
@@ -374,7 +374,7 @@ void make_pull_groups(t_pull *pull,
     }
 }
 
-void make_pull_coords(t_pull *pull)
+void make_pull_coords(pull_params_t *pull)
 {
     int           c, d;
     t_pull_coord *pcrd;
@@ -425,26 +425,27 @@ void make_pull_coords(t_pull *pull)
 void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real lambda,
                    const output_env_t oenv)
 {
-    t_mdatoms    *md;
-    t_pull       *pull;
-    t_pbc         pbc;
-    int           c;
-    double        t_start;
-
-    init_pull(NULL, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0);
-    md = init_mdatoms(NULL, mtop, ir->efep);
+    pull_params_t *pull;
+    struct pull_t *pull_work;
+    t_mdatoms     *md;
+    t_pbc          pbc;
+    int            c;
+    double         t_start;
+
+    pull      = ir->pull;
+    pull_work = init_pull(NULL, pull, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0);
+    md        = init_mdatoms(NULL, mtop, ir->efep);
     atoms2md(mtop, ir, 0, NULL, mtop->natoms, md);
     if (ir->efep)
     {
         update_mdatoms(md, lambda);
     }
-    pull = ir->pull;
 
     set_pbc(&pbc, ir->ePBC, box);
 
     t_start = ir->init_t + ir->init_step*ir->delta_t;
 
-    pull_calc_coms(NULL, pull, md, &pbc, t_start, x, NULL);
+    pull_calc_coms(NULL, pull_work, md, &pbc, t_start, x, NULL);
 
     fprintf(stderr, "Pull group  natoms  pbc atom  distance at start  reference at t=0\n");
     for (c = 0; c < pull->ncoord; c++)
@@ -469,7 +470,7 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l
             pcrd->init = 0;
         }
 
-        get_pull_coord_value(pull, c, &pbc, &value);
+        get_pull_coord_value(pull_work, c, &pbc, &value);
         fprintf(stderr, " %10.3f nm", value);
 
         if (pcrd->bStart)
@@ -478,4 +479,6 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l
         }
         fprintf(stderr, "     %10.3f nm\n", pcrd->init);
     }
+
+    finish_pull(pull_work);
 }

diff --git a/src/gromacs/legacyheaders/types/inputrec.h b/src/gromacs/legacyheaders/types/inputrec.h
index 95cfc402cd00bb5cc3177bea6bfa1dc5c0218e38..72d85c0119ce6ef5fdd83e12bb2a3ee07ea81bd1 100644 (file)
--- a/src/gromacs/legacyheaders/types/inputrec.h
+++ b/src/gromacs/legacyheaders/types/inputrec.h
@@ -98,31 +98,12 @@ typedef struct {
     gmx_bool    *bTS;
 } t_grpopts;
 
-enum {
-    epgrppbcNONE, epgrppbcREFAT, epgrppbcCOS
-};
-
 typedef struct {
     int         nat;        /* Number of atoms in the pull group */
     atom_id    *ind;        /* The global atoms numbers */
-    int         nat_loc;    /* Number of local pull atoms */
-    int         nalloc_loc; /* Allocation size for ind_loc and weight_loc */
-    atom_id    *ind_loc;    /* Local pull indices */
     int         nweight;    /* The number of weights (0 or nat) */
     real       *weight;     /* Weights (use all 1 when weight==NULL) */
-    real       *weight_loc; /* Weights for the local indices */
-    int         epgrppbc;   /* The type of pbc for this pull group, see enum above */
     atom_id     pbcatom;    /* The reference atom for pbc (global number) */
-
-    /* Variables not present in mdp, but used at run time */
-    gmx_bool    bCalcCOM;   /* Calculate COM? Not if only used as cylinder group */
-    real        mwscale;    /* mass*weight scaling factor 1/sum w m */
-    real        wscale;     /* scaling factor for the weights: sum w m/sum w w m */
-    real        invtm;      /* inverse total mass of the group: 1/wscale sum w m */
-    dvec       *mdw;        /* mass*gradient(weight) for atoms */
-    double     *dv;         /* distance to the other group along vec */
-    dvec        x;          /* center of mass before update */
-    dvec        xp;         /* center of mass after update before constraining */
 } t_pull_group;
 
 typedef struct {
@@ -137,16 +118,6 @@ typedef struct {
     real        rate;       /* Rate of motion (nm/ps) */
     real        k;          /* force constant */
     real        kB;         /* force constant for state B */
-
-    /* Variables not present in mdp, but used at run time */
-    double      value_ref;  /* The reference value, usually init+rate*t */
-    double      value;      /* The current value of the coordinate */
-    dvec        dr;         /* The distance from the reference group */
-    double      vec_len;    /* Length of vec for direction-relative */
-    dvec        ffrad;      /* conversion factor from vec to radial force */
-    double      cyl_dev;    /* The deviation from the reference position */
-    double      f_scal;     /* Scalar force for directional pulling */
-    dvec        f;          /* force due to the pulling/constraining */
 } t_pull_coord;
 
 typedef struct {
@@ -229,27 +200,13 @@ typedef struct {
     gmx_bool       bPrintComp;     /* Print cartesian components for each coord with geometry=distance */
     int            nstxout;        /* Output frequency for pull x */
     int            nstfout;        /* Output frequency for pull f */
-    int            ePBC;           /* the boundary conditions */
-    int            npbcdim;        /* do pbc in dims 0 <= dim < npbcdim */
-    gmx_bool       bRefAt;         /* do we need reference atoms for a group COM ? */
-    int            cosdim;         /* dimension for cosine weighting, -1 if none */
+
     t_pull_group  *group;          /* groups to pull/restrain/etc/ */
     t_pull_coord  *coord;          /* the pull coordinates */
+} pull_params_t;
 
-    /* Variables not present in mdp, but used at run time */
-    gmx_bool       bPotential;   /* Are there coordinates with potential? */
-    gmx_bool       bConstraint;  /* Are there constrained coordinates? */
-    gmx_bool       bCylinder;    /* Is group 0 a cylinder group? */
-    t_pull_group  *dyna;         /* dynamic groups for use with local constraints */
-    gmx_bool       bSetPBCatoms; /* Do we need to set x_pbc for the groups? */
-
-    rvec          *rbuf;         /* COM calculation buffer */
-    dvec          *dbuf;         /* COM calculation buffer */
-    double        *dbuf_cyl;     /* cylinder ref. groups COM calculation buffer */
-
-    FILE          *out_x;        /* output file for pull data */
-    FILE          *out_f;        /* output file for pull data */
-} t_pull;
+/* Abstract type for COM pull caclulations only defined in the pull module */
+struct pull_t;
 
 
 /* Abstract types for enforced rotation only defined in pull_rotation.c       */
@@ -460,7 +417,9 @@ typedef struct {
     real            wall_density[2];         /* Number density for walls                     */
     real            wall_ewald_zfac;         /* Scaling factor for the box for Ewald         */
     gmx_bool        bPull;                   /* Do we do COM pulling?                        */
-    t_pull         *pull;                    /* The data for center of mass pulling          */
+    pull_params_t  *pull;                    /* The data for center of mass pulling          */
+    struct pull_t  *pull_work;               /* The COM pull force calculation data structure; TODO this pointer should live somewhere else */
+
     gmx_bool        bRot;                    /* Calculate enforced rotation potential(s)?    */
     t_rot          *rot;                     /* The data for enforced rotation potentials    */
     int             eSwapCoords;             /* Do ion/water position exchanges (CompEL)?    */

diff --git a/src/gromacs/mdlib/broadcaststructs.cpp b/src/gromacs/mdlib/broadcaststructs.cpp
index c13b7185fe306c6c5fe805cbcb06df36dedd2c75..049fc4cba74efcc2ae11a5d264c403d0b50a024a 100644 (file)
--- a/src/gromacs/mdlib/broadcaststructs.cpp
+++ b/src/gromacs/mdlib/broadcaststructs.cpp
@@ -493,7 +493,7 @@ static void bc_pull_group(const t_commrec *cr, t_pull_group *pgrp)
     }
 }
 
-static void bc_pull(const t_commrec *cr, t_pull *pull)
+static void bc_pull(const t_commrec *cr, pull_params_t *pull)
 {
     int g;
 

diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp
index ef655c062d06dbfcfdb584b9606b2475f1ed8ac1..e4d2f9d3c5528f71cf384cd6a11b3e5a88ddf5bf 100644 (file)
--- a/src/gromacs/mdlib/constr.cpp
+++ b/src/gromacs/mdlib/constr.cpp
@@ -614,7 +614,7 @@ gmx_bool constrain(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
 
     if (econq == econqCoord)
     {
-        if (ir->bPull && ir->pull->bConstraint)
+        if (ir->bPull && pull_have_constraint(ir->pull_work))
         {
             if (EI_DYNAMICS(ir->eI))
             {
@@ -625,7 +625,7 @@ gmx_bool constrain(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
                 t = ir->init_t;
             }
             set_pbc(&pbc, ir->ePBC, box);
-            pull_constraint(ir->pull, md, &pbc, cr, ir->delta_t, t, x, xprime, v, *vir);
+            pull_constraint(ir->pull_work, md, &pbc, cr, ir->delta_t, t, x, xprime, v, *vir);
         }
         if (constr->ed && delta_step > 0)
         {
@@ -1172,7 +1172,7 @@ gmx_constr_t init_constraints(FILE *fplog,
     nset = gmx_mtop_ftype_count(mtop, F_SETTLE);
 
     if (ncon+nset == 0 &&
-        !(ir->bPull && ir->pull->bConstraint) &&
+        !(ir->bPull && pull_have_constraint(ir->pull_work)) &&
         ed == NULL)
     {
         return NULL;

diff --git a/src/gromacs/mdlib/mdebin.c b/src/gromacs/mdlib/mdebin.c
index d46d50f63772830550cdc198011d221c38e6d87d..cef361bfd6f67d066eb3c1c7f08ad37998e320dc 100644 (file)
--- a/src/gromacs/mdlib/mdebin.c
+++ b/src/gromacs/mdlib/mdebin.c
@@ -57,6 +57,7 @@
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/mdebin_bar.h"
 #include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pulling/pull.h"
 #include "gromacs/topology/mtop_util.h"
 #include "gromacs/utility/smalloc.h"
 
@@ -291,7 +292,7 @@ t_mdebin *init_mdebin(ener_file_t       fp_ene,
         }
         else if (i == F_COM_PULL)
         {
-            md->bEner[i] = (ir->bPull && ir->pull->bPotential);
+            md->bEner[i] = (ir->bPull && pull_have_potential(ir->pull_work));
         }
         else if (i == F_ECONSERVED)
         {

diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp
index fd4074dde3c7673e0bbe692ac9b196093a680b94..08bc5071d31d1557ee774831d27cda14ecd7328f 100644 (file)
--- a/src/gromacs/mdlib/sim_util.cpp
+++ b/src/gromacs/mdlib/sim_util.cpp
@@ -332,7 +332,7 @@ static void pull_potential_wrapper(t_commrec *cr,
     set_pbc(&pbc, ir->ePBC, box);
     dvdl                     = 0;
     enerd->term[F_COM_PULL] +=
-        pull_potential(ir->pull, mdatoms, &pbc,
+        pull_potential(ir->pull_work, mdatoms, &pbc,
                        cr, t, lambda[efptRESTRAINT], x, f, vir_force, &dvdl);
     enerd->dvdl_lin[efptRESTRAINT] += dvdl;
     wallcycle_stop(wcycle, ewcPULLPOT);
@@ -1155,9 +1155,9 @@ void do_force_cutsVERLET(FILE *fplog, t_commrec *cr,
         clear_rvec(fr->vir_diag_posres);
     }
 
-    if (inputrec->bPull && inputrec->pull->bConstraint)
+    if (inputrec->bPull && pull_have_constraint(inputrec->pull_work))
     {
-        clear_pull_forces(inputrec->pull);
+        clear_pull_forces(inputrec->pull_work);
     }
 
     /* We calculate the non-bonded forces, when done on the CPU, here.
@@ -1455,7 +1455,7 @@ void do_force_cutsVERLET(FILE *fplog, t_commrec *cr,
         }
     }
 
-    if (inputrec->bPull && inputrec->pull->bPotential)
+    if (inputrec->bPull && pull_have_potential(inputrec->pull_work))
     {
         /* Since the COM pulling is always done mass-weighted, no forces are
          * applied to vsites and this call can be done after vsite spreading.
@@ -1798,9 +1798,9 @@ void do_force_cutsGROUP(FILE *fplog, t_commrec *cr,
 
         clear_rvec(fr->vir_diag_posres);
     }
-    if (inputrec->bPull && inputrec->pull->bConstraint)
+    if (inputrec->bPull && pull_have_constraint(inputrec->pull_work))
     {
-        clear_pull_forces(inputrec->pull);
+        clear_pull_forces(inputrec->pull_work);
     }
 
     /* update QMMMrec, if necessary */
@@ -1919,7 +1919,7 @@ void do_force_cutsGROUP(FILE *fplog, t_commrec *cr,
         }
     }
 
-    if (inputrec->bPull && inputrec->pull->bPotential)
+    if (inputrec->bPull && pull_have_potential(inputrec->pull_work))
     {
         pull_potential_wrapper(cr, inputrec, box, x,
                                f, vir_force, mdatoms, enerd, lambda, t,

diff --git a/src/gromacs/pulling/pull.cpp b/src/gromacs/pulling/pull.cpp
index 4635654dc8236452b6d7ad19fde557a5a4779d4b..3c970da0e43bfa5edd4b00df9e7bf36376bbc3b3 100644 (file)
--- a/src/gromacs/pulling/pull.cpp
+++ b/src/gromacs/pulling/pull.cpp
@@ -59,12 +59,14 @@
 #include "gromacs/legacyheaders/types/commrec.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pulling/pull_internal.h"
 #include "gromacs/topology/mtop_util.h"
 #include "gromacs/utility/futil.h"
 #include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/smalloc.h"
 
-static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp)
+static void pull_print_group_x(FILE *out, ivec dim,
+                               const pull_group_work_t *pgrp)
 {
     int m;
 
@@ -77,7 +79,7 @@ static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp)
     }
 }
 
-static void pull_print_coord_dr(FILE *out, const t_pull_coord *pcrd,
+static void pull_print_coord_dr(FILE *out, const pull_coord_work_t *pcrd,
                                 gmx_bool bPrintRefValue,
                                 gmx_bool bPrintComponents)
 {
@@ -94,7 +96,7 @@ static void pull_print_coord_dr(FILE *out, const t_pull_coord *pcrd,
 
         for (m = 0; m < DIM; m++)
         {
-            if (pcrd->dim[m])
+            if (pcrd->params.dim[m])
             {
                 fprintf(out, "\t%g", pcrd->dr[m]);
             }
@@ -102,39 +104,43 @@ static void pull_print_coord_dr(FILE *out, const t_pull_coord *pcrd,
     }
 }
 
-static void pull_print_x(FILE *out, t_pull *pull, double t)
+static void pull_print_x(FILE *out, struct pull_t *pull, double t)
 {
-    int           c;
-    t_pull_coord *pcrd;
+    int c;
 
     fprintf(out, "%.4f", t);
 
     for (c = 0; c < pull->ncoord; c++)
     {
+        pull_coord_work_t *pcrd;
+
         pcrd = &pull->coord[c];
 
-        if (pull->bPrintCOM1)
+        if (pull->params.bPrintCOM1)
         {
-            if (pcrd->eGeom == epullgCYL)
+            if (pcrd->params.eGeom == epullgCYL)
             {
-                pull_print_group_x(out, pcrd->dim, &pull->dyna[c]);
+                pull_print_group_x(out, pcrd->params.dim, &pull->dyna[c]);
             }
             else
             {
-                pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[0]]);
+                pull_print_group_x(out, pcrd->params.dim,
+                                   &pull->group[pcrd->params.group[0]]);
             }
         }
-        if (pull->bPrintCOM2)
+        if (pull->params.bPrintCOM2)
         {
-            pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[1]]);
+            pull_print_group_x(out, pcrd->params.dim,
+                               &pull->group[pcrd->params.group[1]]);
         }
         pull_print_coord_dr(out, pcrd,
-                            pull->bPrintRefValue, pull->bPrintComp);
+                            pull->params.bPrintRefValue,
+                            pull->params.bPrintComp);
     }
     fprintf(out, "\n");
 }
 
-static void pull_print_f(FILE *out, t_pull *pull, double t)
+static void pull_print_f(FILE *out, struct pull_t *pull, double t)
 {
     int c;
 
@@ -147,20 +153,20 @@ static void pull_print_f(FILE *out, t_pull *pull, double t)
     fprintf(out, "\n");
 }
 
-void pull_print_output(t_pull *pull, gmx_int64_t step, double time)
+void pull_print_output(struct pull_t *pull, gmx_int64_t step, double time)
 {
-    if ((pull->nstxout != 0) && (step % pull->nstxout == 0))
+    if ((pull->params.nstxout != 0) && (step % pull->params.nstxout == 0))
     {
         pull_print_x(pull->out_x, pull, time);
     }
 
-    if ((pull->nstfout != 0) && (step % pull->nstfout == 0))
+    if ((pull->params.nstfout != 0) && (step % pull->params.nstfout == 0))
     {
         pull_print_f(pull->out_f, pull, time);
     }
 }
 
-static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv,
+static FILE *open_pull_out(const char *fn, struct pull_t *pull, const output_env_t oenv,
                            gmx_bool bCoord, unsigned long Flags)
 {
     FILE  *fp;
@@ -199,12 +205,12 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv
                  * the data in print_pull_x above.
                  */
 
-                if (pull->bPrintCOM1)
+                if (pull->params.bPrintCOM1)
                 {
                     /* Legend for reference group position */
                     for (m = 0; m < DIM; m++)
                     {
-                        if (pull->coord[c].dim[m])
+                        if (pull->coord[c].params.dim[m])
                         {
                             sprintf(buf, "%d %s%d%c", c+1, "c", 1, 'X'+m);
                             setname[nsets] = gmx_strdup(buf);
@@ -212,12 +218,12 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv
                         }
                     }
                 }
-                if (pull->bPrintCOM2)
+                if (pull->params.bPrintCOM2)
                 {
                     /* Legend for reference group position */
                     for (m = 0; m < DIM; m++)
                     {
-                        if (pull->coord[c].dim[m])
+                        if (pull->coord[c].params.dim[m])
                         {
                             sprintf(buf, "%d %s%d%c", c+1, "c", 2, 'X'+m);
                             setname[nsets] = gmx_strdup(buf);
@@ -229,18 +235,18 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv
                 sprintf(buf, "%d", c+1);
                 setname[nsets] = gmx_strdup(buf);
                 nsets++;
-                if (pull->bPrintRefValue)
+                if (pull->params.bPrintRefValue)
                 {
                     sprintf(buf, "%c%d", 'r', c+1);
                     setname[nsets] = gmx_strdup(buf);
                     nsets++;
                 }
-                if (pull->bPrintComp)
+                if (pull->params.bPrintComp)
                 {
                     /* The pull coord distance components */
                     for (m = 0; m < DIM; m++)
                     {
-                        if (pull->coord[c].dim[m])
+                        if (pull->coord[c].params.dim[m])
                         {
                             sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
                             setname[nsets] = gmx_strdup(buf);
@@ -272,7 +278,8 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv
 }
 
 /* Apply forces in a mass weighted fashion */
-static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md,
+static void apply_forces_grp(const pull_group_work_t *pgrp,
+                             const t_mdatoms *md,
                              const dvec f_pull, int sign, rvec *f)
 {
     int    i, ii, m;
@@ -280,7 +287,7 @@ static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md,
 
     inv_wm = pgrp->mwscale;
 
-    if (pgrp->nat == 1 && pgrp->nat_loc == 1)
+    if (pgrp->params.nat == 1 && pgrp->nat_loc == 1)
     {
         /* Only one atom and our rank has this atom: we can skip
          * the mass weighting, which means that this code also works
@@ -311,7 +318,8 @@ static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md,
 }
 
 /* Apply forces in a mass weighted fashion to a cylinder group */
-static void apply_forces_cyl_grp(const t_pull_group *pgrp, double dv_corr,
+static void apply_forces_cyl_grp(const pull_group_work_t *pgrp,
+                                 double dv_corr,
                                  const t_mdatoms *md,
                                  const dvec f_pull, double f_scal,
                                  int sign, rvec *f)
@@ -347,10 +355,10 @@ static void apply_forces_cyl_grp(const t_pull_group *pgrp, double dv_corr,
 /* Apply torque forces in a mass weighted fashion to the groups that define
  * the pull vector direction for pull coordinate pcrd.
  */
-static void apply_forces_vec_torque(const t_pull       *pull,
-                                    const t_pull_coord *pcrd,
-                                    const t_mdatoms    *md,
-                                    rvec               *f)
+static void apply_forces_vec_torque(const struct pull_t     *pull,
+                                    const pull_coord_work_t *pcrd,
+                                    const t_mdatoms         *md,
+                                    rvec                    *f)
 {
     double inpr;
     int    m;
@@ -376,21 +384,22 @@ static void apply_forces_vec_torque(const t_pull       *pull,
     }
 
     /* Apply the force to the groups defining the vector using opposite signs */
-    apply_forces_grp(&pull->group[pcrd->group[2]], md, f_perp, -1, f);
-    apply_forces_grp(&pull->group[pcrd->group[3]], md, f_perp,  1, f);
+    apply_forces_grp(&pull->group[pcrd->params.group[2]], md, f_perp, -1, f);
+    apply_forces_grp(&pull->group[pcrd->params.group[3]], md, f_perp,  1, f);
 }
 
 /* Apply forces in a mass weighted fashion */
-static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f)
+static void apply_forces(struct pull_t * pull, t_mdatoms * md, rvec *f)
 {
-    int                 c;
-    const t_pull_coord *pcrd;
+    int c;
 
     for (c = 0; c < pull->ncoord; c++)
     {
+        const pull_coord_work_t *pcrd;
+
         pcrd = &pull->coord[c];
 
-        if (pcrd->eGeom == epullgCYL)
+        if (pcrd->params.eGeom == epullgCYL)
         {
             dvec f_tot;
             int  m;
@@ -403,11 +412,11 @@ static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f)
             {
                 f_tot[m] = pcrd->f[m] + pcrd->f_scal*pcrd->ffrad[m];
             }
-            apply_forces_grp(&pull->group[pcrd->group[1]], md, f_tot, 1, f);
+            apply_forces_grp(&pull->group[pcrd->params.group[1]], md, f_tot, 1, f);
         }
         else
         {
-            if (pcrd->eGeom == epullgDIRRELATIVE)
+            if (pcrd->params.eGeom == epullgDIRRELATIVE)
             {
                 /* We need to apply the torque forces to the pull groups
                  * that define the pull vector.
@@ -415,16 +424,17 @@ static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f)
                 apply_forces_vec_torque(pull, pcrd, md, f);
             }
 
-            if (pull->group[pcrd->group[0]].nat > 0)
+            if (pull->group[pcrd->params.group[0]].params.nat > 0)
             {
-                apply_forces_grp(&pull->group[pcrd->group[0]], md, pcrd->f, -1, f);
+                apply_forces_grp(&pull->group[pcrd->params.group[0]], md, pcrd->f, -1, f);
             }
-            apply_forces_grp(&pull->group[pcrd->group[1]], md, pcrd->f, 1, f);
+            apply_forces_grp(&pull->group[pcrd->params.group[1]], md, pcrd->f, 1, f);
         }
     }
 }
 
-static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc)
+static double max_pull_distance2(const pull_coord_work_t *pcrd,
+                                 const t_pbc             *pbc)
 {
     double max_d2;
     int    m;
@@ -433,7 +443,7 @@ static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc)
 
     for (m = 0; m < pbc->ndim_ePBC; m++)
     {
-        if (pcrd->dim[m] != 0)
+        if (pcrd->params.dim[m] != 0)
         {
             max_d2 = std::min(max_d2, static_cast<double>(norm2(pbc->box[m])));
         }
@@ -445,31 +455,31 @@ static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc)
 /* This function returns the distance based on coordinates xg and xref.
  * Note that the pull coordinate struct pcrd is not modified.
  */
-static void low_get_pull_coord_dr(const t_pull *pull,
-                                  const t_pull_coord *pcrd,
+static void low_get_pull_coord_dr(const struct pull_t *pull,
+                                  const pull_coord_work_t *pcrd,
                                   const t_pbc *pbc,
                                   dvec xg, dvec xref, double max_dist2,
                                   dvec dr)
 {
-    const t_pull_group *pgrp0;
-    int                 m;
-    dvec                xrefr, dref = {0, 0, 0};
-    double              dr2;
+    const pull_group_work_t *pgrp0;
+    int                      m;
+    dvec                     xrefr, dref = {0, 0, 0};
+    double                   dr2;
 
-    pgrp0 = &pull->group[pcrd->group[0]];
+    pgrp0 = &pull->group[pcrd->params.group[0]];
 
     /* Only the first group can be an absolute reference, in that case nat=0 */
-    if (pgrp0->nat == 0)
+    if (pgrp0->params.nat == 0)
     {
         for (m = 0; m < DIM; m++)
         {
-            xref[m] = pcrd->origin[m];
+            xref[m] = pcrd->params.origin[m];
         }
     }
 
     copy_dvec(xref, xrefr);
 
-    if (pcrd->eGeom == epullgDIRPBC)
+    if (pcrd->params.eGeom == epullgDIRPBC)
     {
         for (m = 0; m < DIM; m++)
         {
@@ -483,16 +493,17 @@ static void low_get_pull_coord_dr(const t_pull *pull,
     dr2 = 0;
     for (m = 0; m < DIM; m++)
     {
-        dr[m] *= pcrd->dim[m];
+        dr[m] *= pcrd->params.dim[m];
         dr2   += dr[m]*dr[m];
     }
     if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
     {
         gmx_fatal(FARGS, "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the box size (%f).\nYou might want to consider using \"pull-geometry = direction-periodic\" instead.\n",
-                  pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2));
+                  pcrd->params.group[0], pcrd->params.group[1],
+                  sqrt(dr2), sqrt(max_dist2));
     }
 
-    if (pcrd->eGeom == epullgDIRPBC)
+    if (pcrd->params.eGeom == epullgDIRPBC)
     {
         dvec_inc(dr, dref);
     }
@@ -501,16 +512,16 @@ static void low_get_pull_coord_dr(const t_pull *pull,
 /* This function returns the distance based on the contents of the pull struct.
  * pull->coord[coord_ind].dr, and potentially vector, are updated.
  */
-static void get_pull_coord_dr(t_pull      *pull,
-                              int          coord_ind,
-                              const t_pbc *pbc)
+static void get_pull_coord_dr(struct pull_t *pull,
+                              int            coord_ind,
+                              const t_pbc   *pbc)
 {
-    double        md2;
-    t_pull_coord *pcrd;
+    double             md2;
+    pull_coord_work_t *pcrd;
 
     pcrd = &pull->coord[coord_ind];
 
-    if (pcrd->eGeom == epullgDIRPBC)
+    if (pcrd->params.eGeom == epullgDIRPBC)
     {
         md2 = -1;
     }
@@ -519,21 +530,21 @@ static void get_pull_coord_dr(t_pull      *pull,
         md2 = max_pull_distance2(pcrd, pbc);
     }
 
-    if (pcrd->eGeom == epullgDIRRELATIVE)
+    if (pcrd->params.eGeom == epullgDIRRELATIVE)
     {
         /* We need to determine the pull vector */
-        const t_pull_group *pgrp2, *pgrp3;
-        dvec                vec;
-        int                 m;
+        const pull_group_work_t *pgrp2, *pgrp3;
+        dvec                     vec;
+        int                      m;
 
-        pgrp2 = &pull->group[pcrd->group[2]];
-        pgrp3 = &pull->group[pcrd->group[3]];
+        pgrp2 = &pull->group[pcrd->params.group[2]];
+        pgrp3 = &pull->group[pcrd->params.group[3]];
 
         pbc_dx_d(pbc, pgrp3->x, pgrp2->x, vec);
 
         for (m = 0; m < DIM; m++)
         {
-            vec[m] *= pcrd->dim[m];
+            vec[m] *= pcrd->params.dim[m];
         }
         pcrd->vec_len = dnorm(vec);
         for (m = 0; m < DIM; m++)
@@ -550,36 +561,36 @@ static void get_pull_coord_dr(t_pull      *pull,
     }
 
     low_get_pull_coord_dr(pull, pcrd, pbc,
-                          pull->group[pcrd->group[1]].x,
-                          pcrd->eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x,
+                          pull->group[pcrd->params.group[1]].x,
+                          pcrd->params.eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->params.group[0]].x,
                           md2,
                           pcrd->dr);
 }
 
-static void update_pull_coord_reference_value(t_pull_coord *pcrd, double t)
+static void update_pull_coord_reference_value(pull_coord_work_t *pcrd, double t)
 {
     /* With zero rate the reference value is set initially and doesn't change */
-    if (pcrd->rate != 0)
+    if (pcrd->params.rate != 0)
     {
-        pcrd->value_ref = pcrd->init + pcrd->rate*t;
+        pcrd->value_ref = pcrd->params.init + pcrd->params.rate*t;
     }
 }
 
 /* Calculates pull->coord[coord_ind].value.
  * This function also updates pull->coord[coord_ind].dr.
  */
-static void get_pull_coord_distance(t_pull      *pull,
-                                    int          coord_ind,
-                                    const t_pbc *pbc)
+static void get_pull_coord_distance(struct pull_t *pull,
+                                    int            coord_ind,
+                                    const t_pbc   *pbc)
 {
-    t_pull_coord *pcrd;
-    int           m;
+    pull_coord_work_t *pcrd;
+    int                m;
 
     pcrd = &pull->coord[coord_ind];
 
     get_pull_coord_dr(pull, coord_ind, pbc);
 
-    switch (pcrd->eGeom)
+    switch (pcrd->params.eGeom)
     {
         case epullgDIST:
             /* Pull along the vector between the com's */
@@ -602,15 +613,15 @@ static void get_pull_coord_distance(t_pull      *pull,
 /* Returns the deviation from the reference value.
  * Updates pull->coord[coord_ind].dr, .value and .value_ref.
  */
-static double get_pull_coord_deviation(t_pull      *pull,
-                                       int          coord_ind,
-                                       const t_pbc *pbc,
-                                       double       t)
+static double get_pull_coord_deviation(struct pull_t *pull,
+                                       int            coord_ind,
+                                       const t_pbc   *pbc,
+                                       double         t)
 {
-    static gmx_bool  bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
-                                         but is fairly benign */
-    t_pull_coord    *pcrd;
-    double           dev = 0;
+    static gmx_bool    bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
+                                           but is fairly benign */
+    pull_coord_work_t *pcrd;
+    double             dev = 0;
 
     pcrd = &pull->coord[coord_ind];
 
@@ -618,7 +629,7 @@ static double get_pull_coord_deviation(t_pull      *pull,
 
     update_pull_coord_reference_value(pcrd, t);
 
-    switch (pcrd->eGeom)
+    switch (pcrd->params.eGeom)
     {
         case epullgDIST:
             /* Pull along the vector between the com's */
@@ -652,17 +663,17 @@ static double get_pull_coord_deviation(t_pull      *pull,
     return dev;
 }
 
-void get_pull_coord_value(t_pull      *pull,
-                          int          coord_ind,
-                          const t_pbc *pbc,
-                          double      *value)
+void get_pull_coord_value(struct pull_t *pull,
+                          int            coord_ind,
+                          const t_pbc   *pbc,
+                          double        *value)
 {
     get_pull_coord_distance(pull, coord_ind, pbc);
 
     *value = pull->coord[coord_ind].value;
 }
 
-void clear_pull_forces(t_pull *pull)
+void clear_pull_forces(struct pull_t *pull)
 {
     int i;
 
@@ -678,7 +689,7 @@ void clear_pull_forces(t_pull *pull)
 }
 
 /* Apply constraint using SHAKE */
-static void do_constraint(t_pull *pull, t_pbc *pbc,
+static void do_constraint(struct pull_t *pull, t_pbc *pbc,
                           rvec *x, rvec *v,
                           gmx_bool bMaster, tensor vir,
                           double dt, double t)
@@ -695,8 +706,6 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
     int           niter = 0, g, c, ii, j, m, max_iter = 100;
     double        a;
     dvec          tmp, tmp3;
-    t_pull_group *pgrp0, *pgrp1;
-    t_pull_coord *pcrd;
 
     snew(r_ij,   pull->ncoord);
     snew(dr_tot, pull->ncoord);
@@ -714,9 +723,11 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
     /* Determine the constraint directions from the old positions */
     for (c = 0; c < pull->ncoord; c++)
     {
+        pull_coord_work_t *pcrd;
+
         pcrd = &pull->coord[c];
 
-        if (pcrd->eType != epullCONSTRAINT)
+        if (pcrd->params.eType != epullCONSTRAINT)
         {
             continue;
         }
@@ -733,7 +744,8 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
                     c, pcrd->dr[XX], pcrd->dr[YY], pcrd->dr[ZZ]);
         }
 
-        if (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgDIRPBC)
+        if (pcrd->params.eGeom == epullgDIR ||
+            pcrd->params.eGeom == epullgDIRPBC)
         {
             /* Select the component along vec */
             a = 0;
@@ -765,24 +777,26 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
         /* loop over all constraints */
         for (c = 0; c < pull->ncoord; c++)
         {
-            dvec dr0, dr1;
+            pull_coord_work_t *pcrd;
+            pull_group_work_t *pgrp0, *pgrp1;
+            dvec               dr0, dr1;
 
             pcrd  = &pull->coord[c];
 
-            if (pcrd->eType != epullCONSTRAINT)
+            if (pcrd->params.eType != epullCONSTRAINT)
             {
                 continue;
             }
 
             update_pull_coord_reference_value(pcrd, t);
 
-            pgrp0 = &pull->group[pcrd->group[0]];
-            pgrp1 = &pull->group[pcrd->group[1]];
+            pgrp0 = &pull->group[pcrd->params.group[0]];
+            pgrp1 = &pull->group[pcrd->params.group[1]];
 
             /* Get the current difference vector */
             low_get_pull_coord_dr(pull, pcrd, pbc,
-                                  rnew[pcrd->group[1]],
-                                  rnew[pcrd->group[0]],
+                                  rnew[pcrd->params.group[1]],
+                                  rnew[pcrd->params.group[0]],
                                   -1, unc_ij);
 
             if (debug)
@@ -792,7 +806,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
 
             rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
 
-            switch (pcrd->eGeom)
+            switch (pcrd->params.eGeom)
             {
                 case epullgDIST:
                     if (pcrd->value_ref <= 0)
@@ -866,8 +880,8 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
             {
                 int g0, g1;
 
-                g0 = pcrd->group[0];
-                g1 = pcrd->group[1];
+                g0 = pcrd->params.group[0];
+                g1 = pcrd->params.group[1];
                 low_get_pull_coord_dr(pull, pcrd, pbc, rnew[g1], rnew[g0], -1, tmp);
                 low_get_pull_coord_dr(pull, pcrd, pbc, dr1, dr0, -1, tmp3);
                 fprintf(debug,
@@ -885,30 +899,32 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
             } /* END DEBUG */
 
             /* Update the COMs with dr */
-            dvec_inc(rnew[pcrd->group[1]], dr1);
-            dvec_inc(rnew[pcrd->group[0]], dr0);
+            dvec_inc(rnew[pcrd->params.group[1]], dr1);
+            dvec_inc(rnew[pcrd->params.group[0]], dr0);
         }
 
         /* Check if all constraints are fullfilled now */
         for (c = 0; c < pull->ncoord; c++)
         {
+            pull_coord_work_t *pcrd;
+
             pcrd = &pull->coord[c];
 
-            if (pcrd->eType != epullCONSTRAINT)
+            if (pcrd->params.eType != epullCONSTRAINT)
             {
                 continue;
             }
 
             low_get_pull_coord_dr(pull, pcrd, pbc,
-                                  rnew[pcrd->group[1]],
-                                  rnew[pcrd->group[0]],
+                                  rnew[pcrd->params.group[1]],
+                                  rnew[pcrd->params.group[0]],
                                   -1, unc_ij);
 
-            switch (pcrd->eGeom)
+            switch (pcrd->params.eGeom)
             {
                 case epullgDIST:
                     bConverged =
-                        fabs(dnorm(unc_ij) - pcrd->value_ref) < pull->constr_tol;
+                        fabs(dnorm(unc_ij) - pcrd->value_ref) < pull->params.constr_tol;
                     break;
                 case epullgDIR:
                 case epullgDIRPBC:
@@ -920,7 +936,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
                     inpr = diprod(unc_ij, vec);
                     dsvmul(inpr, vec, unc_ij);
                     bConverged =
-                        fabs(diprod(unc_ij, vec) - pcrd->value_ref) < pull->constr_tol;
+                        fabs(diprod(unc_ij, vec) - pcrd->value_ref) < pull->params.constr_tol;
                     break;
             }
 
@@ -956,8 +972,8 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
     /* update atoms in the groups */
     for (g = 0; g < pull->ngroup; g++)
     {
-        const t_pull_group *pgrp;
-        dvec                dr;
+        const pull_group_work_t *pgrp;
+        dvec                     dr;
 
         pgrp = &pull->group[g];
 
@@ -996,16 +1012,18 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
     /* calculate the constraint forces, used for output and virial only */
     for (c = 0; c < pull->ncoord; c++)
     {
+        pull_coord_work_t *pcrd;
+
         pcrd = &pull->coord[c];
 
-        if (pcrd->eType != epullCONSTRAINT)
+        if (pcrd->params.eType != epullCONSTRAINT)
         {
             continue;
         }
 
-        pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->group[0]].invtm + pull->group[pcrd->group[1]].invtm)*dt*dt);
+        pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->params.group[0]].invtm + pull->group[pcrd->params.group[1]].invtm)*dt*dt);
 
-        if (vir != NULL && pcrd->eGeom != epullgDIRPBC && bMaster)
+        if (vir != NULL && pcrd->params.eGeom != epullgDIRPBC && bMaster)
         {
             double f_invr;
 
@@ -1030,32 +1048,33 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
 }
 
 /* Pulling with a harmonic umbrella potential or constant force */
-static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda,
+static void do_pull_pot(struct pull_t *pull, t_pbc *pbc, double t, real lambda,
                         real *V, tensor vir, real *dVdl)
 {
-    int           c, j, m;
-    double        dev, ndr, invdr = 0;
-    real          k, dkdl;
-    t_pull_coord *pcrd;
+    int    c, j, m;
+    double dev, ndr, invdr = 0;
+    real   k, dkdl;
 
     /* loop over the pull coordinates */
     *V    = 0;
     *dVdl = 0;
     for (c = 0; c < pull->ncoord; c++)
     {
+        pull_coord_work_t *pcrd;
+
         pcrd = &pull->coord[c];
 
-        if (pcrd->eType == epullCONSTRAINT)
+        if (pcrd->params.eType == epullCONSTRAINT)
         {
             continue;
         }
 
         dev = get_pull_coord_deviation(pull, c, pbc, t);
 
-        k    = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB;
-        dkdl = pcrd->kB - pcrd->k;
+        k    = (1.0 - lambda)*pcrd->params.k + lambda*pcrd->params.kB;
+        dkdl = pcrd->params.kB - pcrd->params.k;
 
-        if (pcrd->eGeom == epullgDIST)
+        if (pcrd->params.eGeom == epullgDIST)
         {
             ndr   = dnorm(pcrd->dr);
             if (ndr > 0)
@@ -1081,14 +1100,14 @@ static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda,
             }
         }
 
-        switch (pcrd->eType)
+        switch (pcrd->params.eType)
         {
             case epullUMBRELLA:
             case epullFLATBOTTOM:
                 /* The only difference between an umbrella and a flat-bottom
                  * potential is that a flat-bottom is zero below zero.
                  */
-                if (pcrd->eType == epullFLATBOTTOM && dev < 0)
+                if (pcrd->params.eType == epullFLATBOTTOM && dev < 0)
                 {
                     dev = 0;
                 }
@@ -1106,7 +1125,7 @@ static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda,
                 gmx_incons("Unsupported pull type in do_pull_pot");
         }
 
-        if (pcrd->eGeom == epullgDIST)
+        if (pcrd->params.eGeom == epullgDIST)
         {
             for (m = 0; m < DIM; m++)
             {
@@ -1121,7 +1140,7 @@ static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda,
             }
         }
 
-        if (vir != NULL && pcrd->eGeom != epullgDIRPBC)
+        if (vir != NULL && pcrd->params.eGeom != epullgDIRPBC)
         {
             /* Add the pull contribution to the virial */
             for (j = 0; j < DIM; j++)
@@ -1135,12 +1154,14 @@ static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda,
     }
 }
 
-real pull_potential(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
+real pull_potential(struct pull_t *pull, t_mdatoms *md, t_pbc *pbc,
                     t_commrec *cr, double t, real lambda,
                     rvec *x, rvec *f, tensor vir, real *dvdlambda)
 {
     real V, dVdl;
 
+    assert(pull != NULL);
+
     pull_calc_coms(cr, pull, md, pbc, t, x, NULL);
 
     do_pull_pot(pull, pbc, t, lambda,
@@ -1157,24 +1178,26 @@ real pull_potential(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
     return (MASTER(cr) ? V : 0.0);
 }
 
-void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
+void pull_constraint(struct pull_t *pull, t_mdatoms *md, t_pbc *pbc,
                      t_commrec *cr, double dt, double t,
                      rvec *x, rvec *xp, rvec *v, tensor vir)
 {
+    assert(pull != NULL);
+
     pull_calc_coms(cr, pull, md, pbc, t, x, xp);
 
     do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
 }
 
 static void make_local_pull_group(gmx_ga2la_t ga2la,
-                                  t_pull_group *pg, int start, int end)
+                                  pull_group_work_t *pg, int start, int end)
 {
     int i, ii;
 
     pg->nat_loc = 0;
-    for (i = 0; i < pg->nat; i++)
+    for (i = 0; i < pg->params.nat; i++)
     {
-        ii = pg->ind[i];
+        ii = pg->params.ind[i];
         if (ga2la)
         {
             if (!ga2la_get_home(ga2la, ii, &ii))
@@ -1189,22 +1212,22 @@ static void make_local_pull_group(gmx_ga2la_t ga2la,
             {
                 pg->nalloc_loc = over_alloc_dd(pg->nat_loc+1);
                 srenew(pg->ind_loc, pg->nalloc_loc);
-                if (pg->epgrppbc == epgrppbcCOS || pg->weight)
+                if (pg->epgrppbc == epgrppbcCOS || pg->params.weight != NULL)
                 {
                     srenew(pg->weight_loc, pg->nalloc_loc);
                 }
             }
             pg->ind_loc[pg->nat_loc] = ii;
-            if (pg->weight)
+            if (pg->params.weight != NULL)
             {
-                pg->weight_loc[pg->nat_loc] = pg->weight[i];
+                pg->weight_loc[pg->nat_loc] = pg->params.weight[i];
             }
             pg->nat_loc++;
         }
     }
 }
 
-void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md)
+void dd_make_local_pull_groups(gmx_domdec_t *dd, struct pull_t *pull, t_mdatoms *md)
 {
     gmx_ga2la_t ga2la;
     int         g;
@@ -1229,14 +1252,15 @@ void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md)
 }
 
 static void init_pull_group_index(FILE *fplog, t_commrec *cr,
-                                  int g, t_pull_group *pg,
+                                  int g, pull_group_work_t *pg,
                                   gmx_bool bConstraint, ivec pulldim_con,
-                                  gmx_mtop_t *mtop, t_inputrec *ir, real lambda)
+                                  const gmx_mtop_t *mtop,
+                                  const t_inputrec *ir, real lambda)
 {
     int                   i, ii, d, nfrozen, ndim;
     real                  m, w, mbd;
     double                tmass, wmass, wwmass;
-    gmx_groups_t         *groups;
+    const gmx_groups_t   *groups;
     gmx_mtop_atomlookup_t alook;
     t_atom               *atom;
 
@@ -1247,9 +1271,9 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr,
          * So we store the real masses in the weights.
          * We do not set nweight, so these weights do not end up in the tpx file.
          */
-        if (pg->nweight == 0)
+        if (pg->params.nweight == 0)
         {
-            snew(pg->weight, pg->nat);
+            snew(pg->params.weight, pg->params.nat);
         }
     }
 
@@ -1262,15 +1286,15 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr,
     }
     else
     {
-        pg->nat_loc = pg->nat;
-        pg->ind_loc = pg->ind;
+        pg->nat_loc = pg->params.nat;
+        pg->ind_loc = pg->params.ind;
         if (pg->epgrppbc == epgrppbcCOS)
         {
-            snew(pg->weight_loc, pg->nat);
+            snew(pg->weight_loc, pg->params.nat);
         }
         else
         {
-            pg->weight_loc = pg->weight;
+            pg->weight_loc = pg->params.weight;
         }
     }
 
@@ -1282,9 +1306,9 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr,
     tmass   = 0;
     wmass   = 0;
     wwmass  = 0;
-    for (i = 0; i < pg->nat; i++)
+    for (i = 0; i < pg->params.nat; i++)
     {
-        ii = pg->ind[i];
+        ii = pg->params.ind[i];
         gmx_mtop_atomnr_to_atom(alook, ii, &atom);
         if (bConstraint && ir->opts.nFreeze)
         {
@@ -1305,9 +1329,9 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr,
         {
             m = (1 - lambda)*atom->m + lambda*atom->mB;
         }
-        if (pg->nweight > 0)
+        if (pg->params.nweight > 0)
         {
-            w = pg->weight[i];
+            w = pg->params.weight[i];
         }
         else
         {
@@ -1316,9 +1340,9 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr,
         if (EI_ENERGY_MINIMIZATION(ir->eI))
         {
             /* Move the mass to the weight */
-            w            *= m;
-            m             = 1;
-            pg->weight[i] = w;
+            w                   *= m;
+            m                    = 1;
+            pg->params.weight[i] = w;
         }
         else if (ir->eI == eiBD)
         {
@@ -1337,9 +1361,9 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr,
                     mbd = ir->delta_t/ir->opts.tau_t[groups->grpnr[egcTC][ii]];
                 }
             }
-            w            *= m/mbd;
-            m             = mbd;
-            pg->weight[i] = w;
+            w                   *= m/mbd;
+            m                    = mbd;
+            pg->params.weight[i] = w;
         }
         tmass  += m;
         wmass  += m*w;
@@ -1353,7 +1377,7 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr,
         /* We can have single atom groups with zero mass with potential pulling
          * without cosine weighting.
          */
-        if (pg->nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
+        if (pg->params.nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
         {
             /* With one atom the mass doesn't matter */
             wwmass = 1;
@@ -1361,14 +1385,16 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr,
         else
         {
             gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
-                      pg->weight ? " weighted" : "", g);
+                      pg->params.weight ? " weighted" : "", g);
         }
     }
     if (fplog)
     {
         fprintf(fplog,
-                "Pull group %d: %5d atoms, mass %9.3f", g, pg->nat, tmass);
-        if (pg->weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
+                "Pull group %d: %5d atoms, mass %9.3f",
+                g, pg->params.nat, tmass);
+        if (pg->params.weight ||
+            EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
         {
             fprintf(fplog, ", weighted mass %9.3f", wmass*wmass/wwmass);
         }
@@ -1389,7 +1415,7 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr,
         ndim = 0;
         for (d = 0; d < DIM; d++)
         {
-            ndim += pulldim_con[d]*pg->nat;
+            ndim += pulldim_con[d]*pg->params.nat;
         }
         if (fplog && nfrozen > 0 && nfrozen < ndim)
         {
@@ -1404,35 +1430,92 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr,
     }
 }
 
-void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
-               gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
-               gmx_bool bOutFile, unsigned long Flags)
+struct pull_t *
+init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir,
+          int nfile, const t_filenm fnm[],
+          gmx_mtop_t *mtop, t_commrec *cr, const output_env_t oenv, real lambda,
+          gmx_bool bOutFile, unsigned long Flags)
 {
-    t_pull       *pull;
-    t_pull_group *pgrp;
-    int           c, g, m;
+    struct pull_t *pull;
+    int            g, c, m;
+
+    snew(pull, 1);
 
-    pull = ir->pull;
+    /* Copy the pull parameters */
+    pull->params       = *pull_params;
+    /* Avoid pointer copies */
+    pull->params.group = NULL;
+    pull->params.coord = NULL;
+
+    pull->ncoord       = pull_params->ncoord;
+    pull->ngroup       = pull_params->ngroup;
+    snew(pull->coord, pull->ncoord);
+    snew(pull->group, pull->ngroup);
 
     pull->bPotential  = FALSE;
     pull->bConstraint = FALSE;
     pull->bCylinder   = FALSE;
+
+    for (g = 0; g < pull->ngroup; g++)
+    {
+        pull_group_work_t *pgrp;
+        int                i;
+
+        pgrp = &pull->group[g];
+
+        /* Copy the pull group parameters */
+        pgrp->params = pull_params->group[g];
+
+        /* Avoid pointer copies by allocating and copying arrays */
+        snew(pgrp->params.ind, pgrp->params.nat);
+        for (i = 0; i < pgrp->params.nat; i++)
+        {
+            pgrp->params.ind[i] = pull_params->group[g].ind[i];
+        }
+        if (pgrp->params.nweight > 0)
+        {
+            snew(pgrp->params.ind, pgrp->params.nweight);
+            for (i = 0; i < pgrp->params.nweight; i++)
+            {
+                pgrp->params.weight[i] = pull_params->group[g].weight[i];
+            }
+        }
+    }
+
     for (c = 0; c < pull->ncoord; c++)
     {
-        t_pull_coord *pcrd;
-        int           calc_com_start, calc_com_end, g;
+        pull_coord_work_t *pcrd;
+        int                calc_com_start, calc_com_end, g;
 
         pcrd = &pull->coord[c];
 
-        if (pcrd->eType == epullCONSTRAINT)
+        /* Copy all pull coordinate parameters */
+        pcrd->params = pull_params->coord[c];
+
+        switch (pcrd->params.eGeom)
+        {
+            case epullgDIST:
+            case epullgDIRRELATIVE:
+                /* Direction vector is determined at each step */
+                break;
+            case epullgDIR:
+            case epullgDIRPBC:
+            case epullgCYL:
+                copy_rvec(pull_params->coord[c].vec, pcrd->vec);
+                break;
+            default:
+                gmx_incons("Pull geometry not handled");
+        }
+
+        if (pcrd->params.eType == epullCONSTRAINT)
         {
             /* Check restrictions of the constraint pull code */
-            if (pcrd->eGeom == epullgCYL ||
-                pcrd->eGeom == epullgDIRRELATIVE)
+            if (pcrd->params.eGeom == epullgCYL ||
+                pcrd->params.eGeom == epullgDIRRELATIVE)
             {
                 gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
-                          epull_names[pcrd->eType],
-                          epullg_names[pcrd->eGeom],
+                          epull_names[pcrd->params.eType],
+                          epullg_names[pcrd->params.eGeom],
                           epull_names[epullUMBRELLA]);
             }
 
@@ -1443,26 +1526,26 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
             pull->bPotential = TRUE;
         }
 
-        if (pcrd->eGeom == epullgCYL)
+        if (pcrd->params.eGeom == epullgCYL)
         {
             pull->bCylinder = TRUE;
         }
         /* We only need to calculate the plain COM of a group
          * when it is not only used as a cylinder group.
          */
-        calc_com_start = (pcrd->eGeom == epullgCYL         ? 1 : 0);
-        calc_com_end   = (pcrd->eGeom == epullgDIRRELATIVE ? 4 : 2);
+        calc_com_start = (pcrd->params.eGeom == epullgCYL         ? 1 : 0);
+        calc_com_end   = (pcrd->params.eGeom == epullgDIRRELATIVE ? 4 : 2);
 
         for (g = calc_com_start; g <= calc_com_end; g++)
         {
-            pull->group[pcrd->group[g]].bCalcCOM = TRUE;
+            pull->group[pcrd->params.group[g]].bCalcCOM = TRUE;
         }
 
         /* With non-zero rate the reference value is set at every step */
-        if (pcrd->rate == 0)
+        if (pcrd->params.rate == 0)
         {
             /* Initialize the constant reference value */
-            pcrd->value_ref = pcrd->init;
+            pcrd->value_ref = pcrd->params.init;
         }
     }
 
@@ -1481,8 +1564,8 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
         bAbs = FALSE;
         for (c = 0; c < pull->ncoord; c++)
         {
-            if (pull->group[pull->coord[c].group[0]].nat == 0 ||
-                pull->group[pull->coord[c].group[1]].nat == 0)
+            if (pull->group[pull->coord[c].params.group[0]].params.nat == 0 ||
+                pull->group[pull->coord[c].params.group[1]].params.nat == 0)
             {
                 bAbs = TRUE;
             }
@@ -1507,8 +1590,8 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
         bCos = FALSE;
         for (g = 0; g < pull->ngroup; g++)
         {
-            if (pull->group[g].nat > 1 &&
-                pull->group[g].pbcatom < 0)
+            if (pull->group[g].params.nat > 1 &&
+                pull->group[g].params.pbcatom < 0)
             {
                 /* We are using cosine weighting */
                 fprintf(fplog, "Cosine weighting is used for group %d\n", g);
@@ -1528,9 +1611,11 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
     pull->cosdim   = -1;
     for (g = 0; g < pull->ngroup; g++)
     {
+        pull_group_work_t *pgrp;
+
         pgrp           = &pull->group[g];
         pgrp->epgrppbc = epgrppbcNONE;
-        if (pgrp->nat > 0)
+        if (pgrp->params.nat > 0)
         {
             /* There is an issue when a group is used in multiple coordinates
              * and constraints are applied in different dimensions with atoms
@@ -1552,19 +1637,19 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
             clear_ivec(pulldim_con);
             for (c = 0; c < pull->ncoord; c++)
             {
-                if (pull->coord[c].group[0] == g ||
-                    pull->coord[c].group[1] == g ||
-                    (pull->coord[c].eGeom == epullgDIRRELATIVE &&
-                     (pull->coord[c].group[2] == g ||
-                      pull->coord[c].group[3] == g)))
+                if (pull->coord[c].params.group[0] == g ||
+                    pull->coord[c].params.group[1] == g ||
+                    (pull->coord[c].params.eGeom == epullgDIRRELATIVE &&
+                     (pull->coord[c].params.group[2] == g ||
+                      pull->coord[c].params.group[3] == g)))
                 {
                     for (m = 0; m < DIM; m++)
                     {
-                        if (pull->coord[c].dim[m] == 1)
+                        if (pull->coord[c].params.dim[m] == 1)
                         {
                             pulldim[m] = 1;
 
-                            if (pull->coord[c].eType == epullCONSTRAINT)
+                            if (pull->coord[c].params.eType == epullCONSTRAINT)
                             {
                                 bConstraint    = TRUE;
                                 pulldim_con[m] = 1;
@@ -1581,16 +1666,16 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
             for (m = 0; m < pull->npbcdim; m++)
             {
                 GMX_RELEASE_ASSERT(m <= DIM, "npbcdim must be <= the number of dimensions");
-                if (pulldim[m] == 1 && pgrp->nat > 1)
+                if (pulldim[m] == 1 && pgrp->params.nat > 1)
                 {
-                    if (pgrp->pbcatom >= 0)
+                    if (pgrp->params.pbcatom >= 0)
                     {
                         pgrp->epgrppbc = epgrppbcREFAT;
                         pull->bRefAt   = TRUE;
                     }
                     else
                     {
-                        if (pgrp->weight)
+                        if (pgrp->params.weight != NULL)
                         {
                             gmx_fatal(FARGS, "Pull groups can not have relative weights and cosine weighting at same time");
                         }
@@ -1625,13 +1710,13 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
 
         for (c = 0; c < pull->ncoord; c++)
         {
-            const t_pull_coord *pcrd;
+            const pull_coord_work_t *pcrd;
 
             pcrd = &pull->coord[c];
 
-            if (pcrd->eGeom == epullgCYL)
+            if (pcrd->params.eGeom == epullgCYL)
             {
-                if (pull->group[pcrd->group[0]].nat == 0)
+                if (pull->group[pcrd->params.group[0]].params.nat == 0)
                 {
                     gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n");
                 }
@@ -1647,20 +1732,65 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
     pull->out_f = NULL;
     if (bOutFile)
     {
-        if (pull->nstxout > 0)
+        if (pull->params.nstxout > 0)
         {
             pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv,
                                         TRUE, Flags);
         }
-        if (pull->nstfout > 0)
+        if (pull->params.nstfout > 0)
         {
             pull->out_f = open_pull_out(opt2fn("-pf", nfile, fnm), pull, oenv,
                                         FALSE, Flags);
         }
     }
+
+    return pull;
+}
+
+static void destroy_pull_group(pull_group_work_t *pgrp)
+{
+    if (pgrp->ind_loc != pgrp->params.ind)
+    {
+        sfree(pgrp->ind_loc);
+    }
+    if (pgrp->weight_loc != pgrp->params.weight)
+    {
+        sfree(pgrp->weight_loc);
+    }
+    sfree(pgrp->mdw);
+    sfree(pgrp->dv);
+
+    sfree(pgrp->params.ind);
+    sfree(pgrp->params.weight);
+}
+
+static void destroy_pull(struct pull_t *pull)
+{
+    int i;
+
+    for (i = 0; i < pull->ngroup; i++)
+    {
+        destroy_pull_group(&pull->group[i]);
+    }
+    for (i = 0; i < pull->ncoord; i++)
+    {
+        if (pull->coord[i].params.eGeom == epullgCYL)
+        {
+            destroy_pull_group(&(pull->dyna[i]));
+        }
+    }
+    sfree(pull->group);
+    sfree(pull->dyna);
+    sfree(pull->coord);
+
+    sfree(pull->rbuf);
+    sfree(pull->dbuf);
+    sfree(pull->dbuf_cyl);
+
+    sfree(pull);
 }
 
-void finish_pull(t_pull *pull)
+void finish_pull(struct pull_t *pull)
 {
     if (pull->out_x)
     {
@@ -1670,4 +1800,16 @@ void finish_pull(t_pull *pull)
     {
         gmx_fio_fclose(pull->out_f);
     }
+
+    destroy_pull(pull);
+}
+
+gmx_bool pull_have_potential(const struct pull_t *pull)
+{
+    return pull->bPotential;
+}
+
+gmx_bool pull_have_constraint(const struct pull_t *pull)
+{
+    return pull->bConstraint;
 }

diff --git a/src/gromacs/pulling/pull.h b/src/gromacs/pulling/pull.h
index f1f49bd7ace6b6e5bed7898c16f25cf2d504aa41..c73ca640dc7455ee6fb6536fcf2d9061e5d69a1d 100644 (file)
--- a/src/gromacs/pulling/pull.h
+++ b/src/gromacs/pulling/pull.h
@@ -67,7 +67,7 @@ struct t_pbc;
  * \param[in]     pbc       Information structure about periodicity.
  * \param[out]    value     The value of the pull coordinate.
  */
-void get_pull_coord_value(t_pull             *pull,
+void get_pull_coord_value(struct pull_t      *pull,
                           int                 coord_ind,
                           const struct t_pbc *pbc,
                           double             *value);
@@ -77,7 +77,7 @@ void get_pull_coord_value(t_pull             *pull,
  *
  * \param pull              The pull group.
  */
-void clear_pull_forces(t_pull *pull);
+void clear_pull_forces(struct pull_t *pull);
 
 
 /*! \brief Determine the COM pull forces and add them to f, return the potential
@@ -95,7 +95,7 @@ void clear_pull_forces(t_pull *pull);
  *
  * \returns The pull potential energy.
  */
-real pull_potential(t_pull *pull, t_mdatoms *md, struct t_pbc *pbc,
+real pull_potential(struct pull_t *pull, t_mdatoms *md, struct t_pbc *pbc,
                     t_commrec *cr, double t, real lambda,
                     rvec *x, rvec *f, tensor vir, real *dvdlambda);
 
@@ -114,7 +114,7 @@ real pull_potential(t_pull *pull, t_mdatoms *md, struct t_pbc *pbc,
  * \param[in,out] v      Velocities, which may get a pull correction.
  * \param[in,out] vir    The virial, which, if != NULL, gets a pull correction.
  */
-void pull_constraint(t_pull *pull, t_mdatoms *md, struct t_pbc *pbc,
+void pull_constraint(struct pull_t *pull, t_mdatoms *md, struct t_pbc *pbc,
                      t_commrec *cr, double dt, double t,
                      rvec *x, rvec *xp, rvec *v, tensor vir);
 
@@ -127,56 +127,57 @@ void pull_constraint(t_pull *pull, t_mdatoms *md, struct t_pbc *pbc,
  * \param md             All atoms.
  */
 void dd_make_local_pull_groups(gmx_domdec_t *dd,
-                               t_pull *pull, t_mdatoms *md);
-
-
-/*! \brief Get memory and initialize the fields of pull that still need it, and
- * do runtype specific initialization.
- *
- * \param fplog      General output file, normally md.log.
- * \param ir         The inputrec.
- * \param nfile      Number of files.
- * \param fnm        Standard filename struct.
- * \param mtop       The topology of the whole system.
- * \param cr         Struct for communication info.
- * \param oenv       Output options.
- * \param lambda     FEP lambda.
- * \param bOutFile   Open output files?
- * \param Flags      Flags passed over from main, used to determine
- *                   whether or not we are appending.
+                               struct pull_t *pull, t_mdatoms *md);
+
+
+/*! \brief Allocate, initialize and return a pull work struct.
+ *
+ * \param fplog       General output file, normally md.log.
+ * \param pull_params The pull input parameters containing all pull settings.
+ * \param ir          The inputrec.
+ * \param nfile       Number of files.
+ * \param fnm         Standard filename struct.
+ * \param mtop        The topology of the whole system.
+ * \param cr          Struct for communication info.
+ * \param oenv        Output options.
+ * \param lambda      FEP lambda.
+ * \param bOutFile    Open output files?
+ * \param Flags       Flags passed over from main, used to determine
+ *                    whether or not we are appending.
  */
-void init_pull(FILE              *fplog,
-               t_inputrec        *ir,
-               int                nfile,
-               const t_filenm     fnm[],
-               gmx_mtop_t        *mtop,
-               t_commrec        * cr,
-               const output_env_t oenv,
-               real               lambda,
-               gmx_bool           bOutFile,
-               unsigned long      Flags);
+struct pull_t *init_pull(FILE                *fplog,
+                         const pull_params_t *pull_params,
+                         const t_inputrec    *ir,
+                         int                  nfile,
+                         const t_filenm       fnm[],
+                         gmx_mtop_t          *mtop,
+                         t_commrec          * cr,
+                         const output_env_t   oenv,
+                         real                 lambda,
+                         gmx_bool             bOutFile,
+                         unsigned long        Flags);
 
 
 /*! \brief Close the pull output files.
  *
  * \param pull       The pull group.
  */
-void finish_pull(t_pull *pull);
+void finish_pull(struct pull_t *pull);
 
 
 /*! \brief Print the pull output (x and/or f)
  *
- * \param pull     The pull group.
+ * \param pull     The pull data structure.
  * \param step     Time step number.
  * \param time     Time.
  */
-void pull_print_output(t_pull *pull, gmx_int64_t step, double time);
+void pull_print_output(struct pull_t *pull, gmx_int64_t step, double time);
 
 
 /*! \brief Calculates centers of mass all pull groups.
  *
  * \param[in] cr       Struct for communication info.
- * \param[in] pull     The pull group.
+ * \param[in] pull     The pull data structure.
  * \param[in] md       All atoms.
  * \param[in] pbc      Information struct about periodicity.
  * \param[in] t        Time, only used for cylinder ref.
@@ -185,13 +186,27 @@ void pull_print_output(t_pull *pull, gmx_int64_t step, double time);
  *
  */
 void pull_calc_coms(t_commrec        *cr,
-                    t_pull           *pull,
+                    struct pull_t    *pull,
                     t_mdatoms        *md,
                     struct t_pbc     *pbc,
                     double            t,
                     rvec              x[],
                     rvec             *xp);
 
+
+/*! \brief Returns if we have pull coordinates with potential pulling.
+ *
+ * \param[in] pull     The pull data structure.
+ */
+gmx_bool pull_have_potential(const struct pull_t *pull);
+
+
+/*! \brief Returns if we have pull coordinates with constraint pulling.
+ *
+ * \param[in] pull     The pull data structure.
+ */
+gmx_bool pull_have_constraint(const struct pull_t *pull);
+
 #ifdef __cplusplus
 }
 #endif

diff --git a/src/gromacs/pulling/pull_internal.h b/src/gromacs/pulling/pull_internal.h
new file mode 100644 (file)
index 0000000..ad5c10f
--- /dev/null
+++ b/src/gromacs/pulling/pull_internal.h
@@ -0,0 +1,133 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \libinternal \file
+ *
+ *
+ * \brief
+ * This file contains datatypes and function declarations for internal
+   use in the pull code.
+ *
+ * \author Berk Hess
+ *
+ * \inlibraryapi
+ */
+
+#ifndef GMX_PULLING_PULL_INTERNAL_H
+#define GMX_PULLING_PULL_INTERNAL_H
+
+#include "gromacs/legacyheaders/typedefs.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+enum {
+    epgrppbcNONE, epgrppbcREFAT, epgrppbcCOS
+};
+
+typedef struct
+{
+    t_pull_group  params;
+
+    gmx_bool      bCalcCOM;   /* Calculate COM? Not if only used as cylinder group */
+    int           epgrppbc;   /* The type of pbc for this pull group, see enum above */
+
+    int           nat_loc;    /* Number of local pull atoms */
+    int           nalloc_loc; /* Allocation size for ind_loc and weight_loc */
+    atom_id      *ind_loc;    /* Local pull indices */
+    real         *weight_loc; /* Weights for the local indices */
+
+    real          mwscale;    /* mass*weight scaling factor 1/sum w m */
+    real          wscale;     /* scaling factor for the weights: sum w m/sum w w m */
+    real          invtm;      /* inverse total mass of the group: 1/wscale sum w m */
+    dvec         *mdw;        /* mass*gradient(weight) for atoms */
+    double       *dv;         /* distance to the other group along vec */
+    dvec          x;          /* center of mass before update */
+    dvec          xp;         /* center of mass after update before constraining */
+}
+pull_group_work_t;
+
+typedef struct
+{
+    t_pull_coord  params;    /* Pull coordinate (constant) parameters */
+
+    double        value_ref; /* The reference value, usually init+rate*t */
+    double        value;     /* The current value of the coordinate */
+    dvec          dr;        /* The distance from the reference group */
+    rvec          vec;       /* The pull direction */
+    double        vec_len;   /* Length of vec for direction-relative */
+    dvec          ffrad;     /* conversion factor from vec to radial force */
+    double        cyl_dev;   /* The deviation from the reference position */
+    double        f_scal;    /* Scalar force for directional pulling */
+    dvec          f;         /* force due to the pulling/constraining */
+}
+pull_coord_work_t;
+
+struct pull_t
+{
+    pull_params_t      params;       /* The pull parameters, from inputrec */
+
+    gmx_bool           bPotential;   /* Are there coordinates with potential? */
+    gmx_bool           bConstraint;  /* Are there constrained coordinates? */
+
+    int                ePBC;         /* the boundary conditions */
+    int                npbcdim;      /* do pbc in dims 0 <= dim < npbcdim */
+    gmx_bool           bRefAt;       /* do we need reference atoms for a group COM ? */
+    int                cosdim;       /* dimension for cosine weighting, -1 if none */
+
+    int                ngroup;       /* Number of pull groups */
+    int                ncoord;       /* Number of pull coordinates */
+    pull_group_work_t *group;        /* The pull group param and work data */
+    pull_group_work_t *dyna;         /* Dynamic groups for geom=cylinder */
+    pull_coord_work_t *coord;        /* The pull group param and work data */
+
+    gmx_bool           bCylinder;    /* Is group 0 a cylinder group? */
+    gmx_bool           bSetPBCatoms; /* Do we need to set x_pbc for the groups? */
+
+    rvec              *rbuf;         /* COM calculation buffer */
+    dvec              *dbuf;         /* COM calculation buffer */
+    double            *dbuf_cyl;     /* cylinder ref. groups COM calculation buffer */
+    FILE              *out_x;        /* Output file for pull data */
+    FILE              *out_f;        /* Output file for pull data */
+};
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif

diff --git a/src/gromacs/pulling/pullutil.cpp b/src/gromacs/pulling/pullutil.cpp
index 852acc9521ed43646e46570cd21832c20e205870..7345aad8570be6e214a9df46ebe952e8d655431b 100644 (file)
--- a/src/gromacs/pulling/pullutil.cpp
+++ b/src/gromacs/pulling/pullutil.cpp
@@ -49,11 +49,12 @@
 #include "gromacs/math/vec.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/pulling/pull.h"
+#include "gromacs/pulling/pull_internal.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/futil.h"
 #include "gromacs/utility/smalloc.h"
 
-static void pull_set_pbcatom(t_commrec *cr, t_pull_group *pgrp,
+static void pull_set_pbcatom(t_commrec *cr, pull_group_work_t *pgrp,
                              rvec *x,
                              rvec x_pbc)
 {
@@ -61,7 +62,7 @@ static void pull_set_pbcatom(t_commrec *cr, t_pull_group *pgrp,
 
     if (cr != NULL && DOMAINDECOMP(cr))
     {
-        if (ga2la_get_home(cr->dd->ga2la, pgrp->pbcatom, &a))
+        if (ga2la_get_home(cr->dd->ga2la, pgrp->params.pbcatom, &a))
         {
             copy_rvec(x[a], x_pbc);
         }
@@ -72,11 +73,11 @@ static void pull_set_pbcatom(t_commrec *cr, t_pull_group *pgrp,
     }
     else
     {
-        copy_rvec(x[pgrp->pbcatom], x_pbc);
+        copy_rvec(x[pgrp->params.pbcatom], x_pbc);
     }
 }
 
-static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull,
+static void pull_set_pbcatoms(t_commrec *cr, struct pull_t *pull,
                               rvec *x,
                               rvec *x_pbc)
 {
@@ -85,7 +86,7 @@ static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull,
     n = 0;
     for (g = 0; g < pull->ngroup; g++)
     {
-        if (!pull->group[g].bCalcCOM || pull->group[g].pbcatom == -1)
+        if (!pull->group[g].bCalcCOM || pull->group[g].params.pbcatom == -1)
         {
             clear_rvec(x_pbc[g]);
         }
@@ -103,7 +104,7 @@ static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull,
     }
 }
 
-static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md,
+static void make_cyl_refgrps(t_commrec *cr, struct pull_t *pull, t_mdatoms *md,
                              t_pbc *pbc, double t, rvec *x)
 {
     /* The size and stride per coord for the reduction buffer */
@@ -111,8 +112,6 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md,
     int           c, i, ii, m, start, end;
     rvec          g_x, dx, dir;
     double        inv_cyl_r2;
-    t_pull_coord *pcrd;
-    t_pull_group *pref, *pgrp, *pdyna;
     gmx_ga2la_t   ga2la = NULL;
 
     if (pull->dbuf_cyl == NULL)
@@ -128,13 +127,14 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md,
     start = 0;
     end   = md->homenr;
 
-    inv_cyl_r2 = 1/dsqr(pull->cylinder_r);
+    inv_cyl_r2 = 1/dsqr(pull->params.cylinder_r);
 
     /* loop over all groups to make a reference group for each*/
     for (c = 0; c < pull->ncoord; c++)
     {
-        double sum_a, wmass, wwmass;
-        dvec   radf_fac0, radf_fac1;
+        pull_coord_work_t *pcrd;
+        double             sum_a, wmass, wwmass;
+        dvec               radf_fac0, radf_fac1;
 
         pcrd   = &pull->coord[c];
 
@@ -144,11 +144,13 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md,
         clear_dvec(radf_fac0);
         clear_dvec(radf_fac1);
 
-        if (pcrd->eGeom == epullgCYL)
+        if (pcrd->params.eGeom == epullgCYL)
         {
+            pull_group_work_t *pref, *pgrp, *pdyna;
+
             /* pref will be the same group for all pull coordinates */
-            pref  = &pull->group[pcrd->group[0]];
-            pgrp  = &pull->group[pcrd->group[1]];
+            pref  = &pull->group[pcrd->params.group[0]];
+            pgrp  = &pull->group[pcrd->params.group[1]];
             pdyna = &pull->dyna[c];
             copy_rvec(pcrd->vec, dir);
             pdyna->nat_loc = 0;
@@ -158,10 +160,10 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md,
              * we reduced the other group COM over the ranks. This resolves
              * any PBC issues and we don't need to use a PBC-atom here.
              */
-            if (pcrd->rate != 0)
+            if (pcrd->params.rate != 0)
             {
                 /* With rate=0, value_ref is set initially */
-                pcrd->value_ref = pcrd->init + pcrd->rate*t;
+                pcrd->value_ref = pcrd->params.init + pcrd->params.rate*t;
             }
             for (m = 0; m < DIM; m++)
             {
@@ -169,12 +171,12 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md,
             }
 
             /* loop over all atoms in the main ref group */
-            for (i = 0; i < pref->nat; i++)
+            for (i = 0; i < pref->params.nat; i++)
             {
-                ii = pref->ind[i];
+                ii = pref->params.ind[i];
                 if (ga2la)
                 {
-                    if (!ga2la_get_home(ga2la, pref->ind[i], &ii))
+                    if (!ga2la_get_home(ga2la, pref->params.ind[i], &ii))
                     {
                         ii = -1;
                     }
@@ -260,14 +262,17 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md,
 
     for (c = 0; c < pull->ncoord; c++)
     {
+        pull_coord_work_t *pcrd;
+
         pcrd  = &pull->coord[c];
 
-        if (pcrd->eGeom == epullgCYL)
+        if (pcrd->params.eGeom == epullgCYL)
         {
-            double wmass, wwmass, dist;
+            pull_group_work_t *pdyna, *pgrp;
+            double             wmass, wwmass, dist;
 
             pdyna = &pull->dyna[c];
-            pgrp  = &pull->group[pcrd->group[1]];
+            pgrp  = &pull->group[pcrd->params.group[1]];
 
             wmass          = pull->dbuf_cyl[c*stride+0];
             wwmass         = pull->dbuf_cyl[c*stride+1];
@@ -327,7 +332,7 @@ static double atan2_0_2pi(double y, double x)
 
 /* calculates center of mass of selection index from all coordinates x */
 void pull_calc_coms(t_commrec *cr,
-                    t_pull *pull, t_mdatoms *md, t_pbc *pbc, double t,
+                    struct pull_t *pull, t_mdatoms *md, t_pbc *pbc, double t,
                     rvec x[], rvec *xp)
 {
     int  g;
@@ -377,7 +382,7 @@ void pull_calc_coms(t_commrec *cr,
 
     for (g = 0; g < pull->ngroup; g++)
     {
-        t_pull_group *pgrp;
+        pull_group_work_t *pgrp;
 
         pgrp = &pull->group[g];
 
@@ -467,7 +472,7 @@ void pull_calc_coms(t_commrec *cr,
                  * Note that with constraint pulling the mass does matter, but
                  * in that case a check group mass != 0 has been done before.
                  */
-                if (pgrp->nat == 1 && pgrp->nat_loc == 1 && wmass == 0)
+                if (pgrp->params.nat == 1 && pgrp->nat_loc == 1 && wmass == 0)
                 {
                     int m;
 
@@ -554,10 +559,10 @@ void pull_calc_coms(t_commrec *cr,
 
     for (g = 0; g < pull->ngroup; g++)
     {
-        t_pull_group *pgrp;
+        pull_group_work_t *pgrp;
 
         pgrp = &pull->group[g];
-        if (pgrp->nat > 0 && pgrp->bCalcCOM)
+        if (pgrp->params.nat > 0 && pgrp->bCalcCOM)
         {
             if (pgrp->epgrppbc != epgrppbcCOS)
             {

diff --git a/src/gromacs/tools/compare.c b/src/gromacs/tools/compare.c
index 7be156ab53ed9f33facf9bf66715ffbd33d9908f..fd73cd8c9799a2363e5014504eef53f4f5292607 100644 (file)
--- a/src/gromacs/tools/compare.c
+++ b/src/gromacs/tools/compare.c
@@ -888,7 +888,7 @@ static void cmp_inputrec(FILE *fp, t_inputrec *ir1, t_inputrec *ir2, real ftol,
     cmp_cosines(fp, "et", ir1->et, ir2->et, ftol, abstol);
 }
 
-static void comp_pull_AB(FILE *fp, t_pull *pull, real ftol, real abstol)
+static void comp_pull_AB(FILE *fp, pull_params_t *pull, real ftol, real abstol)
 {
     int i;
 

diff --git a/src/programs/mdrun/md.cpp b/src/programs/mdrun/md.cpp
index c377038774f94a683c65df23c39995bd31dd395f..a839ee8adbbce5c2e0513c7ba12243d98d6eab6c 100644 (file)
--- a/src/programs/mdrun/md.cpp
+++ b/src/programs/mdrun/md.cpp
@@ -1603,7 +1603,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             }
             if (ir->bPull)
             {
-                pull_print_output(ir->pull, step, t);
+                pull_print_output(ir->pull_work, step, t);
             }
 
             if (do_per_step(step, ir->nstlog))

diff --git a/src/programs/mdrun/runner.cpp b/src/programs/mdrun/runner.cpp
index 613611896c8266d8846b4d31175c519460c03909..d88c5cd656b7742142896f0247cb1082e4455064 100644 (file)
--- a/src/programs/mdrun/runner.cpp
+++ b/src/programs/mdrun/runner.cpp
@@ -1250,8 +1250,10 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
         if (inputrec->bPull)
         {
             /* Initialize pull code */
-            init_pull(fplog, inputrec, nfile, fnm, mtop, cr, oenv, inputrec->fepvals->init_lambda,
-                      EI_DYNAMICS(inputrec->eI) && MASTER(cr), Flags);
+            inputrec->pull_work =
+                init_pull(fplog, inputrec->pull, inputrec, nfile, fnm,
+                          mtop, cr, oenv, inputrec->fepvals->init_lambda,
+                          EI_DYNAMICS(inputrec->eI) && MASTER(cr), Flags);
         }
 
         if (inputrec->bRot)
@@ -1298,7 +1300,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
 
         if (inputrec->bPull)
         {
-            finish_pull(inputrec->pull);
+            finish_pull(inputrec->pull_work);
         }
 
         if (inputrec->bRot)