The van der Waals radius must equal the Coulomb radius with Verlet
authorCarsten Kutzner <ckutzne@gwdg.de>
Fri, 14 Mar 2014 10:27:36 +0000 (11:27 +0100)
committerCarsten Kutzner <ckutzne@gwdg.de>
Fri, 14 Mar 2014 10:49:21 +0000 (11:49 +0100)
This patch fixes an issue that can occur when using g_tune_pme
with a Verlet pair-list .tpr input file. If the .tpr itself has
large cutoffs (e.g. 1.2 nm) and one asks g_tune_pme to scale
_down_ the Coulomb radius, the van der Waals radius is not
scaled down with the Coulomb radius (only upscaling worked).
One ends up with an unusable .tpr file because rVdW != rCoul.
This patch ensures that van der Waals and Coulomb radii
are always equal with Verlet pair-lists. Fixes #1460. Thanks
to Joao Rodrigues for reporting the issue!

Change-Id: I5ef30e71a35cd83838040057e16e52e09ea82e9a

src/tools/gmx_tune_pme.c

index 1787953a84cfdc925c25e68aed52109a1aca99e5..44dd73657081e09790b7b037c23bbd97ed3ffdce 100644 (file)
@@ -997,7 +997,7 @@ static void make_benchmark_tprs(
             ir->nkx = ir->nky = ir->nkz = 0;
             calc_grid(NULL, state.box, fourierspacing*fac, &ir->nkx, &ir->nky, &ir->nkz);
 
-            /* Adjust other radii since various conditions neet to be fulfilled */
+            /* Adjust other radii since various conditions need to be fulfilled */
             if (eelPME == ir->coulombtype)
             {
                 /* plain PME, rcoulomb must be equal to rlist */
@@ -1011,8 +1011,16 @@ static void make_benchmark_tprs(
 
             if (bScaleRvdw && evdwCUT == ir->vdwtype)
             {
-                /* For vdw cutoff, rvdw >= rlist */
-                ir->rvdw = max(info->rvdw[0], ir->rlist);
+                if ( ecutsVERLET == ir->cutoff_scheme)
+                {
+                    /* With Verlet, the van der Waals radius must always equal the Coulomb radius */
+                    ir->rvdw = ir->rcoulomb;
+                }
+                else
+                {
+                    /* For vdw cutoff, rvdw >= rlist */
+                    ir->rvdw = max(info->rvdw[0], ir->rlist);
+                }
             }
 
             ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));