[gmx-developers] generalized reaction-field implementation

David van der Spoel spoel at xray.bmc.uu.se
Tue Jul 12 21:22:52 CEST 2016

On 12/07/16 15:42, Jan Henning Peters wrote:
> 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
> while
> #define BOLTZ            (RGAS/KILO)
> substituting this into the source code, we get
> *kappa  =
> 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
Isn't kappa printed to the log file if you use it in the mdp file?
How does it compare to your expectation?

David van der Spoel, Ph.D., Professor of Biology
Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205.
spoel at xray.bmc.uu.se    http://folding.bmc.uu.se

More information about the gromacs.org_gmx-developers mailing list