[gmx-developers] generalized reaction-field implementation

[Mark Abraham]

Hi,

Great. Is there something we can document/report better?

Mark

On Wed, Jul 13, 2016 at 12:08 PM Jan Henning Peters <
jan.henning.peters at fu-berlin.de> wrote:

> 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
>
> --
> Gromacs Developers mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers
> or send a mail to gmx-developers-request at gromacs.org.
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-developers/attachments/20160713/cabd6904/attachment-0003.html>