Fixing handling of perturbation mass changes.

Fixes redmine #1232

in force.c, sum_dhdl

 * moved F_DKDL to match the order in efpt_names.  Not required, but
   harmonizes the code (lack of clarity probably helped cause the
   problems before), has no code effect.

 * no longer treating the F_DKDL term separately from the other
   derivative components.  Will be added to F_DVDL if the mass-lambda
   term is not separately specified.  Results in a bit of a misnomer
   (F_DVDL becomes the derivative of the entire hamiltonian), but
   makes it much easier to collapse all molecular perturbation terms
   into a single component for output, where it is no longer really
   F_DVDL.  I think that's better than always printing out a F_DVDL
   and a F_DKDL for everything where F_DKDL will probably usually
   be zero.

in md_support.c, compute_globals

 * Synchronize the behaviors of the dhdls by writing first to the linear component
   corresponding to the mass, and then later transferring it to F_DKDL

in group.h, struct gmx_ekindata_t
 * add pointer to per-thread accumulation variable for dekindl

in tgroup.c, sum_ekin

 * For velocity verlet integrators, computes the dekindl correctly as
   the derivatives of the current ekin.  Shouldn't really affect the results
   in any significant way, since the average contribution will be the same
   regardless, but this is more consistent.

in tgroup.c, init_ekindata

 * reduce use of numeric constants in allocating memory

 * initialize new ekindata_t member

in update.c, calc_ke_part_normal

 * zero the accumulator for dekindl before using it,
   fixing bug introduced in 7b6508e8

in update.c, in calc_ke_part_normal and calc_ke_part_visc

 * sign error in mass change; if mass B is greater than mass A, then the
   change in free energy is positive, not negative.

Change-Id: I9deaf546bca66d400e0eb2c4015abeeda302dd1d

diff --git a/include/types/group.h b/include/types/group.h
index 299d295f367d6c9b3c98d9888b82d39b1c335ea6..c6ab45f1e0ce2876b2286c5ec8ad64d33105452d 100644 (file)
--- a/include/types/group.h
+++ b/include/types/group.h
@@ -74,8 +74,9 @@ typedef struct {
     gmx_bool         bNEMD;
     int              ngtc;            /* The number of T-coupling groups      */
     t_grp_tcstat    *tcstat;          /* T-coupling data            */
-    tensor         **ekin_work_alloc; /* Allocated locations of ekin_work   */
+    tensor         **ekin_work_alloc; /* Allocated locations for *_work members */
     tensor         **ekin_work;       /* Work arrays for tcstat per thread    */
+    real           **dekindl_work;    /* Work location for dekindl per thread */
     int              ngacc;           /* The number of acceleration groups    */
     t_grp_acc       *grpstat;         /* Acceleration data                     */
     tensor           ekin;            /* overall kinetic energy               */

diff --git a/src/mdlib/force.c b/src/mdlib/force.c
index 37ba7c3ec2fc18acad5f72901ed7ffe48c46dce5..c331ecb552bcaaef896a00046c73434c416b5bc7 100644 (file)
--- a/src/mdlib/force.c
+++ b/src/mdlib/force.c
@@ -813,6 +813,9 @@ void sum_dhdl(gmx_enerdata_t *enerd, real *lambda, t_lambda *fepvals)
             /* could this be done more readably/compactly? */
             switch (i)
             {
+                case (efptMASS):
+                    index = F_DKDL;
+                    break;
                 case (efptCOUL):
                     index = F_DVDL_COUL;
                     break;
@@ -825,9 +828,6 @@ void sum_dhdl(gmx_enerdata_t *enerd, real *lambda, t_lambda *fepvals)
                 case (efptRESTRAINT):
                     index = F_DVDL_RESTRAINT;
                     break;
-                case (efptMASS):
-                    index = F_DKDL;
-                    break;
                 default:
                     index = F_DVDL;
                     break;
@@ -868,21 +868,13 @@ void sum_dhdl(gmx_enerdata_t *enerd, real *lambda, t_lambda *fepvals)
                                                  so we don't need to add anything to the
                                                  enerd->enerpart_lambda[0] */
 
-        /* we don't need to worry about dvdl contributions to the current lambda, because
-           it's automatically zero */
-
-        /* first kinetic energy term */
-        dlam = (fepvals->all_lambda[efptMASS][i] - lambda[efptMASS]);
-
-        enerd->enerpart_lambda[i+1] += enerd->term[F_DKDL]*dlam;
+        /* we don't need to worry about dvdl_lin contributions to dE at
+           current lambda, because the contributions to the current
+           lambda are automatically zeroed */
 
         for (j = 0; j < efptNR; j++)
         {
-            if (j == efptMASS)
-            {
-                continue;
-            }                            /* no other mass term to worry about */
-
+            /* Note that this loop is over all dhdl components, not just the separated ones */
             dlam = (fepvals->all_lambda[j][i]-lambda[j]);
             enerd->enerpart_lambda[i+1] += dlam*enerd->dvdl_lin[j];
             if (debug)

diff --git a/src/mdlib/md_support.c b/src/mdlib/md_support.c
index 4d04534b87a8fe05e0dc97a8e8deb31b25bfa8c6..4ece30354c7130030a5ae601864948279abbd875 100644 (file)
--- a/src/mdlib/md_support.c
+++ b/src/mdlib/md_support.c
@@ -293,7 +293,7 @@ void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inpu
     gmx_bool bEner, bPres, bTemp, bVV;
     gmx_bool bRerunMD, bStopCM, bGStat, bIterate,
              bFirstIterate, bReadEkin, bEkinAveVel, bScaleEkin, bConstrain;
-    real     ekin, temp, prescorr, enercorr, dvdlcorr;
+    real     ekin, temp, prescorr, enercorr, dvdlcorr, dvdl_ekin;
 
     /* translate CGLO flags to gmx_booleans */
     bRerunMD = flags & CGLO_RERUNMD;
@@ -450,8 +450,9 @@ void compute_globals(FILE *fplog, gmx_global_stat_t gstat, t_commrec *cr, t_inpu
            bSaveEkinOld: If TRUE (in the case of iteration = bIterate is TRUE), we don't reset the ekinscale_nhc.
            If FALSE, we go ahead and erase over it.
          */
-        enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, &(enerd->term[F_DKDL]),
+        enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, &dvdl_ekin,
                                        bEkinAveVel, bIterate, bScaleEkin);
+        enerd->dvdl_lin[efptMASS] = (double) dvdl_ekin;
 
         enerd->term[F_EKIN] = trace(ekind->ekin);
     }

diff --git a/src/mdlib/tgroup.c b/src/mdlib/tgroup.c
index 84c3d80a812485776f82df48c53ff0d32e7cae90..e644f02d3da4ef4f5ca2b71ce5b221337a5d1294 100644 (file)
--- a/src/mdlib/tgroup.c
+++ b/src/mdlib/tgroup.c
@@ -130,15 +130,24 @@ void init_ekindata(FILE *log, gmx_mtop_t *mtop, t_grpopts *opts,
 
     snew(ekind->ekin_work_alloc, nthread);
     snew(ekind->ekin_work, nthread);
+    snew(ekind->dekindl_work, nthread);
 #pragma omp parallel for num_threads(nthread) schedule(static)
     for (thread = 0; thread < nthread; thread++)
     {
-        /* Allocate 2 elements extra on both sides,
-         * so in single precision we have 2*3*3*4=72 bytes buffer
-         * on both sides to avoid cache pollution.
+#define EKIN_WORK_BUFFER_SIZE 2
+        /* Allocate 2 extra elements on both sides, so in single
+         * precision we have
+         * EKIN_WORK_BUFFER_SIZE*DIM*DIM*sizeof(real) = 72/144 bytes
+         * buffer on both sides to avoid cache pollution.
          */
-        snew(ekind->ekin_work_alloc[thread], ekind->ngtc+4);
-        ekind->ekin_work[thread] = ekind->ekin_work_alloc[thread] + 2;
+        snew(ekind->ekin_work_alloc[thread], ekind->ngtc+2*EKIN_WORK_BUFFER_SIZE);
+        ekind->ekin_work[thread] = ekind->ekin_work_alloc[thread] + EKIN_WORK_BUFFER_SIZE;
+        /* Nasty hack so we can have the per-thread accumulation
+         * variable for dekindl in the same thread-local cache lines
+         * as the per-thread accumulation tensors for ekin[fh],
+         * because they are accumulated in the same loop. */
+        ekind->dekindl_work[thread] = &(ekind->ekin_work[thread][ekind->ngtc][0][0]);
+#undef EKIN_WORK_BUFFER_SIZE
     }
 
     ekind->ngacc = opts->ngacc;
@@ -267,7 +276,6 @@ real sum_ekin(t_grpopts *opts, gmx_ekindata_t *ekind, real *dekindlambda,
                 }
             }
             else
-
             {
                 /* Calculate the full step Ekin as the average of the half steps */
                 for (j = 0; (j < DIM); j++)
@@ -308,7 +316,14 @@ real sum_ekin(t_grpopts *opts, gmx_ekindata_t *ekind, real *dekindlambda,
     }
     if (dekindlambda)
     {
-        *dekindlambda = 0.5*(ekind->dekindl + ekind->dekindl_old);
+        if (bEkinAveVel)
+        {
+            *dekindlambda = ekind->dekindl;
+        }
+        else
+        {
+            *dekindlambda = 0.5*(ekind->dekindl + ekind->dekindl_old);
+        }
     }
     return T;
 }

diff --git a/src/mdlib/update.c b/src/mdlib/update.c
index 9cf204480bef101d03f651fedb8a6320e56d0bc1..4c32f9661855df6c92560148cb0b1a93216fe588 100644 (file)
--- a/src/mdlib/update.c
+++ b/src/mdlib/update.c
@@ -1010,12 +1010,13 @@ static void calc_ke_part_normal(rvec v[], t_grpopts *opts, t_mdatoms *md,
         end_t   = md->start + ((thread+1)*md->homenr)/nthread;
 
         ekin_sum    = ekind->ekin_work[thread];
-        dekindl_sum = &ekind->ekin_work[thread][opts->ngtc][0][0];
+        dekindl_sum = ekind->dekindl_work[thread];
 
         for (gt = 0; gt < opts->ngtc; gt++)
         {
             clear_mat(ekin_sum[gt]);
         }
+        *dekindl_sum = 0.0;
 
         ga = 0;
         gt = 0;
@@ -1045,7 +1046,7 @@ static void calc_ke_part_normal(rvec v[], t_grpopts *opts, t_mdatoms *md,
             }
             if (md->nMassPerturbed && md->bPerturbed[n])
             {
-                *dekindl_sum -=
+                *dekindl_sum +=
                     0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
             }
         }
@@ -1068,7 +1069,7 @@ static void calc_ke_part_normal(rvec v[], t_grpopts *opts, t_mdatoms *md,
             }
         }
 
-        ekind->dekindl += ekind->ekin_work[thread][opts->ngtc][0][0];
+        ekind->dekindl += *ekind->dekindl_work[thread];
     }
 
     inc_nrnb(nrnb, eNR_EKIN, md->homenr);
@@ -1133,7 +1134,7 @@ static void calc_ke_part_visc(matrix box, rvec x[], rvec v[],
         }
         if (md->nPerturbed && md->bPerturbed[n])
         {
-            dekindl -= 0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
+            dekindl += 0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
         }
     }
     ekind->dekindl = dekindl;