[gmx-developers] generalized reaction-field implementation

[Jan Henning Peters]

Am Dienstag, den 12.07.2016, 21:22 +0200 schrieb David van der Spoel:
> 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
> > (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
> >
> Isn't kappa printed to the log file if you use it in the mdp file?
> How does it compare to your expectation?

It is printed out, and it was the difference between this printed out
value and those I calculated that prompted me to look at the source code
in the first place.

On a second look, I think I have it figured out now. The ionic strength
here seems to be given in 1/nm^3 instead of mol/l. This would introduce
another factor of "AVOGADRO" and the formula in the code works out
again.

Sorry for the confusion.

Jan