Merge "Fixing handling of perturbation mass changes." into release-4-6

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;