Merge branch release-5-1 into release-2016

Change-Id: Ia1a7fad67f0ff11175ea24c46f813a445cd49ed6

diff --git a/src/gromacs/domdec/domdec_topology.cpp b/src/gromacs/domdec/domdec_topology.cpp
index 548e373e8564e9237308a1ace4970eccb2fe027f..46b05ddd71ddf01e40ece4a02723cd8fa42f0143 100644 (file)
--- a/src/gromacs/domdec/domdec_topology.cpp
+++ b/src/gromacs/domdec/domdec_topology.cpp
@@ -2219,8 +2219,6 @@ void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
     if (dd->reverse_top->bExclRequired)
     {
         dd->nbonded_local += nexcl;
-
-        forcerec_set_excl_load(fr, ltop);
     }
 
     ltop->atomtypes  = mtop->atomtypes;

diff --git a/src/gromacs/ewald/long-range-correction.cpp b/src/gromacs/ewald/long-range-correction.cpp
index 7f54c1f9e5a63c1bb886d6aaf1e3bae8fff55910..ae5ff5ed23c2116d688420cff6852c0a1bbaab58 100644 (file)
--- a/src/gromacs/ewald/long-range-correction.cpp
+++ b/src/gromacs/ewald/long-range-correction.cpp
@@ -3,7 +3,7 @@
  *
  * 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
+ * Copyright (c) 2013,2014,2015,2016, 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.
@@ -48,6 +48,7 @@
 #include "gromacs/mdtypes/forcerec.h"
 #include "gromacs/mdtypes/md_enums.h"
 #include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
 
 /* There's nothing special to do here if just masses are perturbed,
  * but if either charge or type is perturbed then the implementation
@@ -59,8 +60,10 @@
  * perturbations. The parameter vectors for LJ-PME are likewise
  * undefined when LJ-PME is not active. This works because
  * bHaveChargeOrTypePerturbed handles the control flow. */
-void ewald_LRcorrection(int start, int end,
-                        t_commrec *cr, int thread, t_forcerec *fr,
+void ewald_LRcorrection(int numAtomsLocal,
+                        t_commrec *cr,
+                        int numThreads, int thread,
+                        t_forcerec *fr,
                         real *chargeA, real *chargeB,
                         real *C6A, real *C6B,
                         real *sigmaA, real *sigmaB,
@@ -75,6 +78,22 @@ void ewald_LRcorrection(int start, int end,
                         real lambda_q, real lambda_lj,
                         real *dvdlambda_q, real *dvdlambda_lj)
 {
+    int numAtomsToBeCorrected;
+    if (calc_excl_corr)
+    {
+        /* We need to correct all exclusion pairs (cutoff-scheme = group) */
+        numAtomsToBeCorrected = excl->nr;
+
+        GMX_RELEASE_ASSERT(numAtomsToBeCorrected >= numAtomsLocal, "We might need to do self-corrections");
+    }
+    else
+    {
+        /* We need to correct only self interactions */
+        numAtomsToBeCorrected = numAtomsLocal;
+    }
+    int         start =  (numAtomsToBeCorrected* thread     )/numThreads;
+    int         end   =  (numAtomsToBeCorrected*(thread + 1))/numThreads;
+
     int         i, i1, i2, j, k, m, iv, jv, q;
     int        *AA;
     double      Vexcl_q, dvdl_excl_q, dvdl_excl_lj; /* Necessary for precision */
@@ -305,7 +324,7 @@ void ewald_LRcorrection(int start, int end,
                 }
             }
             /* Dipole correction on force */
-            if (dipole_coeff != 0)
+            if (dipole_coeff != 0 && i < numAtomsLocal)
             {
                 for (j = 0; (j < DIM); j++)
                 {
@@ -452,7 +471,7 @@ void ewald_LRcorrection(int start, int end,
                 }
             }
             /* Dipole correction on force */
-            if (dipole_coeff != 0)
+            if (dipole_coeff != 0 && i < numAtomsLocal)
             {
                 for (j = 0; (j < DIM); j++)
                 {
@@ -532,7 +551,6 @@ void ewald_LRcorrection(int start, int end,
     if (debug)
     {
         fprintf(debug, "Long Range corrections for Ewald interactions:\n");
-        fprintf(debug, "start=%d,natoms=%d\n", start, end-start);
         fprintf(debug, "q2sum = %g, Vself_q=%g c6sum = %g, Vself_lj=%g\n",
                 L1_q*fr->q2sum[0]+lambda_q*fr->q2sum[1], L1_q*Vself_q[0]+lambda_q*Vself_q[1], L1_lj*fr->c6sum[0]+lambda_lj*fr->c6sum[1], L1_lj*Vself_lj[0]+lambda_lj*Vself_lj[1]);
         fprintf(debug, "Electrostatic Long Range correction: Vexcl=%g\n", Vexcl_q);

diff --git a/src/gromacs/ewald/long-range-correction.h b/src/gromacs/ewald/long-range-correction.h
index fe076ad81a12eadbd527689576afc0d1374225d6..37ca4fb4d0422eda75f69b14f09d6eb74a5fc4f9 100644 (file)
--- a/src/gromacs/ewald/long-range-correction.h
+++ b/src/gromacs/ewald/long-range-correction.h
@@ -3,7 +3,7 @@
  *
  * 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
+ * Copyright (c) 2013,2014,2015,2016, 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.
@@ -66,8 +66,10 @@
  * For both cutoff schemes, but only for Coulomb interactions,
  * calculates correction for surface dipole terms. */
 void
-ewald_LRcorrection(int start, int end,
-                   t_commrec *cr, int thread, t_forcerec *fr,
+ewald_LRcorrection(int numAtomsLocal,
+                   t_commrec *cr,
+                   int numThreads, int thread,
+                   t_forcerec *fr,
                    real *chargeA, real *chargeB,
                    real *C6A, real *C6B,
                    real *sigmaA, real *sigmaB,

diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp
index 107e00a8e7b3393ec0eaf42a07f584d29c6ea4eb..90c4b02f5772714ec0be982d49ea8a27a1b365b1 100644 (file)
--- a/src/gromacs/mdlib/force.cpp
+++ b/src/gromacs/mdlib/force.cpp
@@ -449,8 +449,7 @@ void do_force_lowlevel(t_forcerec *fr,      t_inputrec *ir,
                          * exclusion forces) are calculated, so we can store
                          * the forces in the normal, single fr->f_novirsum array.
                          */
-                        ewald_LRcorrection(fr->excl_load[t], fr->excl_load[t+1],
-                                           cr, t, fr,
+                        ewald_LRcorrection(md->homenr, cr, nthreads, t, fr,
                                            md->chargeA, md->chargeB,
                                            md->sqrt_c6A, md->sqrt_c6B,
                                            md->sigmaA, md->sigmaB,

diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp
index 703a1cd471ad50196031a097fa5ad0d9349ed62e..6b5a102f83a2f9adc68960e8b672e1c253384ca1 100644 (file)
--- a/src/gromacs/mdlib/forcerec.cpp
+++ b/src/gromacs/mdlib/forcerec.cpp
@@ -3186,7 +3186,6 @@ void init_forcerec(FILE              *fp,
 
     fr->nthread_ewc = gmx_omp_nthreads_get(emntBonded);
     snew(fr->ewc_t, fr->nthread_ewc);
-    snew(fr->excl_load, fr->nthread_ewc + 1);
 
     /* fr->ic is used both by verlet and group kernels (to some extent) now */
     init_interaction_const(fp, &fr->ic, fr);
@@ -3238,48 +3237,6 @@ void pr_forcerec(FILE *fp, t_forcerec *fr)
     fflush(fp);
 }
 
-void forcerec_set_excl_load(t_forcerec           *fr,
-                            const gmx_localtop_t *top)
-{
-    const int *ind, *a;
-    int        t, i, j, ntot, n, ntarget;
-
-    ind = top->excls.index;
-    a   = top->excls.a;
-
-    ntot = 0;
-    for (i = 0; i < top->excls.nr; i++)
-    {
-        for (j = ind[i]; j < ind[i+1]; j++)
-        {
-            if (a[j] > i)
-            {
-                ntot++;
-            }
-        }
-    }
-
-    fr->excl_load[0] = 0;
-    n                = 0;
-    i                = 0;
-    for (t = 1; t <= fr->nthread_ewc; t++)
-    {
-        ntarget = (ntot*t)/fr->nthread_ewc;
-        while (i < top->excls.nr && n < ntarget)
-        {
-            for (j = ind[i]; j < ind[i+1]; j++)
-            {
-                if (a[j] > i)
-                {
-                    n++;
-                }
-            }
-            i++;
-        }
-        fr->excl_load[t] = i;
-    }
-}
-
 /* Frees GPU memory and destroys the GPU context.
  *
  * Note that this function needs to be called even if GPUs are not used

diff --git a/src/gromacs/mdlib/minimize.cpp b/src/gromacs/mdlib/minimize.cpp
index 3890e24e3734a167e3e6c9efab91a204945dcbc1..591e20a7dc88cbfa338f971b35b97ed014b2a718 100644 (file)
--- a/src/gromacs/mdlib/minimize.cpp
+++ b/src/gromacs/mdlib/minimize.cpp
@@ -401,8 +401,6 @@ void init_em(FILE *fplog, const char *title,
 
         *top      = gmx_mtop_generate_local_top(top_global, ir->efep != efepNO);
 
-        forcerec_set_excl_load(fr, *top);
-
         setup_bonded_threading(fr, &(*top)->idef);
 
         if (ir->ePBC != epbcNONE && !fr->bMolPBC)

diff --git a/src/gromacs/mdlib/update.cpp b/src/gromacs/mdlib/update.cpp
index 0a65e04d0b396f68d1ab5fbb81558fe750e8ce76..cc01f266c007e30ea5691bdc99942026f6c9950c 100644 (file)
--- a/src/gromacs/mdlib/update.cpp
+++ b/src/gromacs/mdlib/update.cpp
@@ -113,7 +113,8 @@ struct gmx_update_t
 };
 
 
-static void do_update_md(int start, int nrend, double dt,
+static void do_update_md(int start, int nrend,
+                         double dt, int nstpcouple,
                          t_grp_tcstat *tcstat,
                          double nh_vxi[],
                          gmx_bool bNEMD, t_grp_acc *gstat, rvec accel[],
@@ -166,7 +167,7 @@ static void do_update_md(int start, int nrend, double dt,
                 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
                 {
                     vnrel = (lg*vrel[d] + dt*(imass*f[n][d] - 0.5*vxi*vrel[d]
-                                              - iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
+                                              - nstpcouple*iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
                     /* do not scale the mean velocities u */
                     vn             = gstat[ga].u[d] + accel[ga][d]*dt + vnrel;
                     v[n][d]        = vn;
@@ -350,7 +351,8 @@ static void do_update_vv_pos(int start, int nrend, double dt,
     }
 } /* do_update_vv_pos */
 
-static void do_update_visc(int start, int nrend, double dt,
+static void do_update_visc(int start, int nrend,
+                           double dt, int nstpcouple,
                            t_grp_tcstat *tcstat,
                            double nh_vxi[],
                            real invmass[],
@@ -398,7 +400,7 @@ static void do_update_visc(int start, int nrend, double dt,
                 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell))
                 {
                     vn  = (lg*vrel[d] + dt*(imass*f[n][d] - 0.5*vxi*vrel[d]
-                                            - iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
+                                            - nstpcouple*iprod(M[d], vrel)))/(1 + 0.5*vxi*dt);
                     if (d == XX)
                     {
                         vn += vc + dt*cosz*cos_accel;
@@ -1669,7 +1671,8 @@ void update_coords(FILE             *fplog,
                 case (eiMD):
                     if (ekind->cosacc.cos_accel == 0)
                     {
-                        do_update_md(start_th, end_th, dt,
+                        do_update_md(start_th, end_th,
+                                     dt, inputrec->nstpcouple,
                                      ekind->tcstat, state->nosehoover_vxi,
                                      ekind->bNEMD, ekind->grpstat, inputrec->opts.acc,
                                      inputrec->opts.nFreeze,
@@ -1680,7 +1683,8 @@ void update_coords(FILE             *fplog,
                     }
                     else
                     {
-                        do_update_visc(start_th, end_th, dt,
+                        do_update_visc(start_th, end_th,
+                                       dt, inputrec->nstpcouple,
                                        ekind->tcstat, state->nosehoover_vxi,
                                        md->invmass, md->ptype,
                                        md->cTC, state->x, upd->xp, state->v, f, M,

diff --git a/src/gromacs/mdtypes/forcerec.h b/src/gromacs/mdtypes/forcerec.h
index a2aa4f8d747cebce5b599234e970f6a16d2411a7..1ba04de5238db7940e6865baa5d6a92a2e1851d4 100644 (file)
--- a/src/gromacs/mdtypes/forcerec.h
+++ b/src/gromacs/mdtypes/forcerec.h
@@ -423,8 +423,6 @@ typedef struct t_forcerec {
     /* Ewald correction thread local virial and energy data */
     int                  nthread_ewc;
     ewald_corr_thread_t *ewc_t;
-    /* Ewald charge correction load distribution over the threads */
-    int                 *excl_load;
 } t_forcerec;
 
 /* Important: Starting with Gromacs-4.6, the values of c6 and c12 in the nbfp array have

diff --git a/src/programs/mdrun/md.cpp b/src/programs/mdrun/md.cpp
index 2a49aba4a37b3d20df5bb7b49f472fbf1f96000c..06c25f321404e012b870cfff79c5044308431b7e 100644 (file)
--- a/src/programs/mdrun/md.cpp
+++ b/src/programs/mdrun/md.cpp
@@ -421,8 +421,6 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     {
         top = gmx_mtop_generate_local_top(top_global, ir->efep != efepNO);
 
-        forcerec_set_excl_load(fr, top);
-
         state    = serial_init_local_state(state_global);
 
         atoms2md(top_global, ir, 0, NULL, top_global->natoms, mdatoms);