Add value_ref and value to pull_coord_t

This is only code refactoring.

Change-Id: I64b4cb05f69325a1dd1e17f269894cb538663817

diff --git a/src/gromacs/gmxpreprocess/readpull.c b/src/gromacs/gmxpreprocess/readpull.c
index c92eada29253aaa2f526ac9112c528437c47d810..68752443c20234202e2acb25ded4900df37c7f97 100644 (file)
--- a/src/gromacs/gmxpreprocess/readpull.c
+++ b/src/gromacs/gmxpreprocess/readpull.c
@@ -427,13 +427,9 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l
 {
     t_mdatoms    *md;
     t_pull       *pull;
-    t_pull_coord *pcrd;
-    t_pull_group *pgrp0, *pgrp1;
     t_pbc         pbc;
-    int           c, m;
+    int           c;
     double        t_start;
-    real          init = 0;
-    double        dev, value;
 
     init_pull(NULL, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0);
     md = init_mdatoms(NULL, mtop, ir->efep);
@@ -453,6 +449,11 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l
     fprintf(stderr, "Pull group  natoms  pbc atom  distance at start  reference at t=0\n");
     for (c = 0; c < pull->ncoord; c++)
     {
+        t_pull_coord *pcrd;
+        t_pull_group *pgrp0, *pgrp1;
+        double        value;
+        real          init = 0;
+
         pcrd  = &pull->coord[c];
 
         pgrp0 = &pull->group[pcrd->group[0]];
@@ -468,9 +469,7 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l
             pcrd->init = 0;
         }
 
-        get_pull_coord_distance(pull, c, &pbc, t_start, &dev);
-        /* Calculate the value of the coordinate at t_start */
-        value = pcrd->init + t_start*pcrd->rate + dev;
+        get_pull_coord_value(pull, c, &pbc, &value);
         fprintf(stderr, " %10.3f nm", value);
 
         if (pcrd->bStart)

diff --git a/src/gromacs/legacyheaders/types/inputrec.h b/src/gromacs/legacyheaders/types/inputrec.h
index 01b5451b7b85bba4297cb0ff0a9bbe482fcff938..95cfc402cd00bb5cc3177bea6bfa1dc5c0218e38 100644 (file)
--- a/src/gromacs/legacyheaders/types/inputrec.h
+++ b/src/gromacs/legacyheaders/types/inputrec.h
@@ -139,6 +139,8 @@ typedef struct {
     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 */

diff --git a/src/gromacs/pulling/pull.c b/src/gromacs/pulling/pull.c
index d7cf619295323e85a23bfdcf6d937fd8add1bf97..a252388ef5ad0b9c4238e1050494e1bdc00fb5f2 100644 (file)
--- a/src/gromacs/pulling/pull.c
+++ b/src/gromacs/pulling/pull.c
@@ -74,46 +74,20 @@ 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,
-                                gmx_bool bPrintRefValue, double t,
+                                gmx_bool bPrintRefValue,
                                 gmx_bool bPrintComponents)
 {
-    double distance, distance2;
-    int    m;
-
-    if (pcrd->eGeom == epullgDIST)
-    {
-        /* Geometry=distance: print the length of the distance vector */
-        distance2 = 0;
-        for (m = 0; m < DIM; m++)
-        {
-            if (pcrd->dim[m])
-            {
-                distance2 += pcrd->dr[m]*pcrd->dr[m];
-            }
-        }
-        distance = sqrt(distance2);
-    }
-    else
-    {
-        /* Geometry=direction-like: print distance along the pull vector */
-        distance = 0;
-        for (m = 0; m < DIM; m++)
-        {
-            if (pcrd->dim[m])
-            {
-                distance += pcrd->dr[m]*pcrd->vec[m];
-            }
-        }
-    }
-    fprintf(out, "\t%g", distance);
+    fprintf(out, "\t%g", pcrd->value);
 
     if (bPrintRefValue)
     {
-        fprintf(out, "\t%g", pcrd->init + pcrd->rate*t);
+        fprintf(out, "\t%g", pcrd->value_ref);
     }
 
     if (bPrintComponents)
     {
+        int m;
+
         for (m = 0; m < DIM; m++)
         {
             if (pcrd->dim[m])
@@ -151,7 +125,7 @@ static void pull_print_x(FILE *out, t_pull *pull, double t)
             pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[1]]);
         }
         pull_print_coord_dr(out, pcrd,
-                            pull->bPrintRefValue, t, pull->bPrintComp);
+                            pull->bPrintRefValue, pull->bPrintComp);
     }
     fprintf(out, "\n");
 }
@@ -469,7 +443,7 @@ static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc)
  */
 static void low_get_pull_coord_dr(const t_pull *pull,
                                   const t_pull_coord *pcrd,
-                                  const t_pbc *pbc, double t,
+                                  const t_pbc *pbc,
                                   dvec xg, dvec xref, double max_dist2,
                                   dvec dr)
 {
@@ -496,7 +470,7 @@ static void low_get_pull_coord_dr(const t_pull *pull,
     {
         for (m = 0; m < DIM; m++)
         {
-            dref[m] = (pcrd->init + pcrd->rate*t)*pcrd->vec[m];
+            dref[m] = pcrd->value_ref*pcrd->vec[m];
         }
         /* Add the reference position, so we use the correct periodic image */
         dvec_inc(xrefr, dref);
@@ -524,9 +498,9 @@ 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(const t_pull *pull,
-                              int coord_ind,
-                              const t_pbc *pbc, double t)
+static void get_pull_coord_dr(t_pull      *pull,
+                              int          coord_ind,
+                              const t_pbc *pbc)
 {
     double        md2;
     t_pull_coord *pcrd;
@@ -573,51 +547,95 @@ static void get_pull_coord_dr(const t_pull *pull,
         }
     }
 
-    low_get_pull_coord_dr(pull, pcrd, pbc, t,
+    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,
                           md2,
                           pcrd->dr);
 }
 
-void get_pull_coord_distance(const t_pull *pull,
-                             int coord_ind,
-                             const t_pbc *pbc, double t,
-                             double *dev)
+static void update_pull_coord_reference_value(t_pull_coord *pcrd, double t)
 {
-    static gmx_bool     bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
-                                            but is fairly benign */
-    const t_pull_coord *pcrd;
-    int                 m;
-    double              ref, drs, inpr;
+    /* With zero rate the reference value is set initially and doesn't change */
+    if (pcrd->rate != 0)
+    {
+        pcrd->value_ref = pcrd->init + pcrd->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)
+{
+    t_pull_coord *pcrd;
+    int           m;
 
     pcrd = &pull->coord[coord_ind];
 
-    get_pull_coord_dr(pull, coord_ind, pbc, t);
+    get_pull_coord_dr(pull, coord_ind, pbc);
 
-    ref = pcrd->init + pcrd->rate*t;
+    switch (pcrd->eGeom)
+    {
+        case epullgDIST:
+            /* Pull along the vector between the com's */
+            pcrd->value = dnorm(pcrd->dr);
+            break;
+        case epullgDIR:
+        case epullgDIRPBC:
+        case epullgDIRRELATIVE:
+        case epullgCYL:
+            /* Pull along vec */
+            pcrd->value = 0;
+            for (m = 0; m < DIM; m++)
+            {
+                pcrd->value += pcrd->vec[m]*pcrd->dr[m];
+            }
+            break;
+    }
+}
+
+/* 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 gmx_bool  bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
+                                         but is fairly benign */
+    t_pull_coord    *pcrd;
+    double           dev = 0;
+
+    pcrd = &pull->coord[coord_ind];
+
+    get_pull_coord_distance(pull, coord_ind, pbc);
+
+    update_pull_coord_reference_value(pcrd, t);
 
     switch (pcrd->eGeom)
     {
         case epullgDIST:
             /* Pull along the vector between the com's */
-            if (ref < 0 && !bWarned)
+            if (pcrd->value_ref < 0 && !bWarned)
             {
-                fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, ref);
+                fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, pcrd->value_ref);
                 bWarned = TRUE;
             }
-            drs = dnorm(pcrd->dr);
-            if (drs == 0)
+            if (pcrd->value == 0)
             {
                 /* With no vector we can not determine the direction for the force,
                  * so we set the force to zero.
                  */
-                *dev = 0;
+                dev = 0;
             }
             else
             {
                 /* Determine the deviation */
-                *dev = drs - ref;
+                dev = pcrd->value - pcrd->value_ref;
             }
             break;
         case epullgDIR:
@@ -625,14 +643,21 @@ void get_pull_coord_distance(const t_pull *pull,
         case epullgDIRRELATIVE:
         case epullgCYL:
             /* Pull along vec */
-            inpr = 0;
-            for (m = 0; m < DIM; m++)
-            {
-                inpr += pcrd->vec[m]*pcrd->dr[m];
-            }
-            *dev = inpr - ref;
+            dev = pcrd->value - pcrd->value_ref;
             break;
     }
+
+    return dev;
+}
+
+void get_pull_coord_value(t_pull      *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)
@@ -661,7 +686,6 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
     dvec          unc_ij; /* xp[i] com of i this step, before constr.   -> unc_ij  */
     dvec         *rnew;   /* current 'new' positions of the groups */
     double       *dr_tot; /* the total update of the coords */
-    double        ref;
     dvec          vec;
     double        inpr;
     double        lambda, rm, invdt = 0;
@@ -696,8 +720,12 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
             continue;
         }
 
-        get_pull_coord_dr(pull, c, pbc, t);
-        /* Store the difference vector at time t for printing */
+        /* Note that get_pull_coord_distance also sets pcrd->dr and pcrd->value.
+         * We don't modify dr and value anymore, so these values are also used
+         * for printing.
+         */
+        get_pull_coord_distance(pull, c, pbc);
+
         if (debug)
         {
             fprintf(debug, "Pull coord %d dr %f %f %f\n",
@@ -745,17 +773,17 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
                 continue;
             }
 
+            update_pull_coord_reference_value(pcrd, t);
+
             pgrp0 = &pull->group[pcrd->group[0]];
             pgrp1 = &pull->group[pcrd->group[1]];
 
             /* Get the current difference vector */
-            low_get_pull_coord_dr(pull, pcrd, pbc, t,
+            low_get_pull_coord_dr(pull, pcrd, pbc,
                                   rnew[pcrd->group[1]],
                                   rnew[pcrd->group[0]],
                                   -1, unc_ij);
 
-            ref = pcrd->init + pcrd->rate*t;
-
             if (debug)
             {
                 fprintf(debug, "Pull coord %d, iteration %d\n", c, niter);
@@ -766,9 +794,9 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
             switch (pcrd->eGeom)
             {
                 case epullgDIST:
-                    if (ref <= 0)
+                    if (pcrd->value_ref <= 0)
                     {
-                        gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, ref);
+                        gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, pcrd->value_ref);
                     }
 
                     {
@@ -776,7 +804,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
 
                         c_a = diprod(r_ij[c], r_ij[c]);
                         c_b = diprod(unc_ij, r_ij[c])*2;
-                        c_c = diprod(unc_ij, unc_ij) - dsqr(ref);
+                        c_c = diprod(unc_ij, unc_ij) - dsqr(pcrd->value_ref);
 
                         if (c_b < 0)
                         {
@@ -814,7 +842,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
                     }
                     /* Select only the component along the vector */
                     dsvmul(a, vec, unc_ij);
-                    lambda = a - ref;
+                    lambda = a - pcrd->value_ref;
                     if (debug)
                     {
                         fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
@@ -834,15 +862,15 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
 
                 g0 = pcrd->group[0];
                 g1 = pcrd->group[1];
-                low_get_pull_coord_dr(pull, pcrd, pbc, t, rnew[g1], rnew[g0], -1, tmp);
-                low_get_pull_coord_dr(pull, pcrd, pbc, t, dr1, dr0, -1, tmp3);
+                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,
                         "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
                         rnew[g0][0], rnew[g0][1], rnew[g0][2],
                         rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
                 fprintf(debug,
                         "Pull ref %8s %8s %8s   %8s %8s %8s d: %8.5f\n",
-                        "", "", "", "", "", "", ref);
+                        "", "", "", "", "", "", pcrd->value_ref);
                 fprintf(debug,
                         "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
                         dr0[0], dr0[1], dr0[2],
@@ -865,9 +893,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
                 continue;
             }
 
-            ref  = pcrd->init + pcrd->rate*t;
-
-            low_get_pull_coord_dr(pull, pcrd, pbc, t,
+            low_get_pull_coord_dr(pull, pcrd, pbc,
                                   rnew[pcrd->group[1]],
                                   rnew[pcrd->group[0]],
                                   -1, unc_ij);
@@ -875,7 +901,8 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
             switch (pcrd->eGeom)
             {
                 case epullgDIST:
-                    bConverged = fabs(dnorm(unc_ij) - ref) < pull->constr_tol;
+                    bConverged =
+                        fabs(dnorm(unc_ij) - pcrd->value_ref) < pull->constr_tol;
                     break;
                 case epullgDIR:
                 case epullgDIRPBC:
@@ -887,7 +914,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) - ref) < pull->constr_tol;
+                        fabs(diprod(unc_ij, vec) - pcrd->value_ref) < pull->constr_tol;
                     break;
             }
 
@@ -897,7 +924,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
                 {
                     fprintf(debug, "NOT CONVERGED YET: Group %d:"
                             "d_ref = %f, current d = %f\n",
-                            g, ref, dnorm(unc_ij));
+                            g, pcrd->value_ref, dnorm(unc_ij));
                 }
 
                 bConverged_all = FALSE;
@@ -1017,7 +1044,7 @@ static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda,
             continue;
         }
 
-        get_pull_coord_distance(pull, c, pbc, t, &dev);
+        dev = get_pull_coord_deviation(pull, c, pbc, t);
 
         k    = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB;
         dkdl = pcrd->kB - pcrd->k;
@@ -1424,6 +1451,13 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
         {
             pull->group[pcrd->group[g]].bCalcCOM = TRUE;
         }
+
+        /* With non-zero rate the reference value is set at every step */
+        if (pcrd->rate == 0)
+        {
+            /* Initialize the constant reference value */
+            pcrd->value_ref = pcrd->init;
+        }
     }
 
     pull->ePBC = ir->ePBC;

diff --git a/src/gromacs/pulling/pull.h b/src/gromacs/pulling/pull.h
index 72a0f8bfe6a30998ce069507fe52c6b8e0209898..f1f49bd7ace6b6e5bed7898c16f25cf2d504aa41 100644 (file)
--- a/src/gromacs/pulling/pull.h
+++ b/src/gromacs/pulling/pull.h
@@ -59,20 +59,18 @@ extern "C" {
 
 struct t_pbc;
 
-/*! \brief Get the deviation for pull coord coord_ind.
- *
- * This call updates some data in the pull coordinates in \p pull.
+
+/*! \brief Get the value for pull coord coord_ind.
  *
  * \param[in,out] pull      The pull struct.
  * \param[in]     coord_ind Number of the pull coordinate.
  * \param[in]     pbc       Information structure about periodicity.
- * \param[in]     t         Time.
- * \param[out]    dev       The deviation from the reference distance.
+ * \param[out]    value     The value of the pull coordinate.
  */
-void get_pull_coord_distance(const t_pull *pull,
-                             int coord_ind,
-                             const struct t_pbc *pbc, double t,
-                             double *dev);
+void get_pull_coord_value(t_pull             *pull,
+                          int                 coord_ind,
+                          const struct t_pbc *pbc,
+                          double             *value);
 
 
 /*! \brief Set the all the pull forces to zero.

diff --git a/src/gromacs/pulling/pullutil.c b/src/gromacs/pulling/pullutil.c
index 135edf649cc8017e2425e56e997304850f847ec3..abfcda8ae6f39cf86f3d0172e7c4418db7dc14f5 100644 (file)
--- a/src/gromacs/pulling/pullutil.c
+++ b/src/gromacs/pulling/pullutil.c
@@ -158,9 +158,14 @@ 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)
+            {
+                /* With rate=0, value_ref is set initially */
+                pcrd->value_ref = pcrd->init + pcrd->rate*t;
+            }
             for (m = 0; m < DIM; m++)
             {
-                g_x[m] = pgrp->x[m] - pcrd->vec[m]*(pcrd->init + pcrd->rate*t);
+                g_x[m] = pgrp->x[m] - pcrd->vec[m]*pcrd->value_ref;
             }
 
             /* loop over all atoms in the main ref group */
@@ -280,7 +285,7 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md,
             pcrd->cyl_dev  = 0;
             for (m = 0; m < DIM; m++)
             {
-                g_x[m]         = pgrp->x[m] - pcrd->vec[m]*(pcrd->init + pcrd->rate*t);
+                g_x[m]         = pgrp->x[m] - pcrd->vec[m]*pcrd->value_ref;
                 dist           = -pcrd->vec[m]*pull->dbuf_cyl[c*stride+2]*pdyna->mwscale;
                 pdyna->x[m]    = g_x[m] - dist;
                 pcrd->cyl_dev += dist;