From: Erik Lindahl <erik@kth.se>
Date: Wed, 25 Jun 2014 17:41:12 +0000 (+0200)
Subject: Fix LJPME Lorentz-Berthelot combination rules
X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=79ed2b8e54c0e205ea95fef94fbed417b357aa57;p=alexxy%2Fgromacs.git

Fix LJPME Lorentz-Berthelot combination rules

Christian found a bug in the way we evaluate
arithmetic combination rules for LJ-PME reciprocal
space. With this change, both group and verlet kernels
produce virtually identical total energies for both
types of reciprocal-space combination rules. It does
not affect the (default) geometric LJPME case.

Change-Id: I4bb0de1df9588dce3ca94e39a370a0b994b25b53
---

diff --git a/src/gromacs/mdlib/forcerec.c b/src/gromacs/mdlib/forcerec.c
index 3b2d69ca04..2ea1e549c8 100644
--- a/src/gromacs/mdlib/forcerec.c
+++ b/src/gromacs/mdlib/forcerec.c
@@ -218,8 +218,8 @@ static real *mk_nbfp_combination_rule(const gmx_ffparams_t *idef, int comb_rule)
                 sigmaj = pow(c12j / c6j, 1.0/6.0);
                 epsi   = c6i * c6i / c12i;
                 epsj   = c6j * c6j / c12j;
-                c6     = epsi * epsj * pow(0.5*(sigmai+sigmaj), 6);
-                c12    = epsi * epsj * pow(0.5*(sigmai+sigmaj), 12);
+                c6     = sqrt(epsi * epsj) * pow(0.5*(sigmai+sigmaj), 6);
+                c12    = sqrt(epsi * epsj) * pow(0.5*(sigmai+sigmaj), 12);
             }
             C6(nbfp, atnr, i, j)   = c6*6.0;
             C12(nbfp, atnr, i, j)  = c12*12.0;