Fix tabulated wall potential.
authorJakub Krajniak <jkrajniak@gmail.com>
Thu, 25 Feb 2016 15:23:13 +0000 (16:23 +0100)
committerGerrit Code Review <gerrit@gerrit.gromacs.org>
Mon, 14 Mar 2016 14:31:39 +0000 (15:31 +0100)
The dispersion and repulsion of the tabulated wall potential were a factor
too small by a factor 6 and 12 respectively.

Fixes #1912.

Change-Id: I56c4b81e7f40f5d5dc9e0f2a184cf4a1bb01e3e5

src/gromacs/mdlib/wall.c

index 6de4e9e4478e65991a0d5b327c0bde341234699e..882925ce03020f44c537a35cabbe25cce43ab31d 100644 (file)
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2008, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,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.
@@ -177,9 +177,9 @@ real do_walls(t_inputrec *ir, t_forcerec *fr, matrix box, t_mdatoms *md,
                 /* The wall energy groups are always at the end of the list */
                 ggid = gid[i]*ngid + ngid - nwall + w;
                 at   = type[i];
-                /* nbfp now includes the 6.0/12.0 derivative prefactors */
-                Cd = nbfp[ntw[w]+2*at]/6.0;
-                Cr = nbfp[ntw[w]+2*at+1]/12.0;
+                /* nbfp now includes the 6/12 derivative prefactors */
+                Cd = nbfp[ntw[w]+2*at]/6;
+                Cr = nbfp[ntw[w]+2*at+1]/12;
                 if (!((Cd == 0 && Cr == 0) || (egp_flags[ggid] & EGP_EXCL)))
                 {
                     if (w == 0)
@@ -231,8 +231,8 @@ real do_walls(t_inputrec *ir, t_forcerec *fr, matrix box, t_mdatoms *md,
                                 Fp    = Ft + Geps + Heps2;
                                 VV    = Yt + Fp*eps;
                                 FF    = Fp + Geps + 2.0*Heps2;
-                                Vd    = Cd*VV;
-                                Fd    = Cd*FF;
+                                Vd    = 6*Cd*VV;
+                                Fd    = 6*Cd*FF;
                                 /* Repulsion */
                                 nnn   = nnn + 4;
                                 Yt    = VFtab[nnn];
@@ -242,8 +242,8 @@ real do_walls(t_inputrec *ir, t_forcerec *fr, matrix box, t_mdatoms *md,
                                 Fp    = Ft + Geps + Heps2;
                                 VV    = Yt + Fp*eps;
                                 FF    = Fp + Geps + 2.0*Heps2;
-                                Vr    = Cr*VV;
-                                Fr    = Cr*FF;
+                                Vr    = 12*Cr*VV;
+                                Fr    = 12*Cr*FF;
                                 V     = Vd + Vr;
                                 F     = -lamfac*(Fd + Fr)*tabscale;
                             }