[gmx-developers] generalized reaction-field implementation
Jan Henning Peters
jan.henning.peters at fu-berlin.de
Tue Jul 12 15:42:57 CEST 2016
Dear GROMACS developers,
I've been trying to implement the generalized reaction field
electrostatic potential using user tables, but noticed that the value
for the parameter kappa I arrived at strongly differs from that used by
GROMACS (according to md.log) for the same parameters. To understand the
difference, I had a closer look at the implementation (in GROMACS 5.0.7,
5.1.2 and the current git master), and I am a bit confused about it -
maybe someone here can enlighten me.
In the manual (sec 4.1.4), kappa is defined by
kappa^2 = (2*I*F^2)/(epsilon_0*epsilon_rf*R*T)
where F is Faraday's constant, R is the ideal gas constant and T is the
absolute temperature (this is the same as in Tironi et al.). In the
source code (function calc_rffac in src/gromacs/mdlib), this is
implemented as
*kappa = std::sqrt(2*I/(EPSILON0*eps_rf*BOLTZ*Temp));
now EPSILON0 is not the SI value of the dielectric constant, but
defined (in src/gromacs/math/units.h) as
#define EPSILON0
(EPSILON0_SI*NANO*KILO)/(E_CHARGE*E_CHARGE*AVOGADRO)
while
#define BOLTZ (RGAS/KILO)
substituting this into the source code, we get
*kappa =
std::sqrt(2*I*E_CHARGE*E_CHARGE*AVOGADRO*KILO/(EPSILON0_SI*NANO*KILO*eps_rf*RGAS*Temp));
and since the Faraday constant is FARADAY = (E_CHARGE*AVOGADRO),
this results in
*kappa =
std::sqrt(2*I*E_CHARGE*FARADAY/(EPSILON0_SI*NANO*eps_rf*RGAS*Temp));
I have not substituted RGAS, as this seems indeed to be the ideal gas
constant from the original formula (also, since
RGAS=(BOLTZMANN*AVOGADRO), this would actually divide the other FARADAY
into E_CHARGE).
Comparing this with the original formula, the resulting value for kappa
appears to be off by a factor of sqrt(AVOGADRO/NANO) or about 7.76e15.
Suprisingly, this difference seems to result only in a small numerical
change of the resulting potential, but unless I missed something, it
still seems to be not correct.
Regards,
Jan
More information about the gromacs.org_gmx-developers
mailing list