[gmx-developers] generalized reaction-field implementation

[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?


-- 
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