[gmx-developers] generalized reaction-field implementation

[Jan Henning Peters]

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