<div dir="ltr">Hi,<div><br></div><div>Great. Is there something we can document/report better?</div><div><br></div><div>Mark</div></div><br><div class="gmail_quote"><div dir="ltr">On Wed, Jul 13, 2016 at 12:08 PM Jan Henning Peters <<a href="mailto:jan.henning.peters@fu-berlin.de">jan.henning.peters@fu-berlin.de</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Am Dienstag, den 12.07.2016, 21:22 +0200 schrieb David van der Spoel:<br>
> On 12/07/16 15:42, Jan Henning Peters wrote:<br>
> > Dear GROMACS developers,<br>
> ><br>
> > I've been trying to implement the generalized reaction field<br>
> > electrostatic potential using user tables, but noticed that the value<br>
> > for the parameter kappa I arrived at strongly differs from that used by<br>
> > GROMACS (according to md.log) for the same parameters. To understand the<br>
> > difference, I had a closer look at the implementation (in GROMACS 5.0.7,<br>
> > 5.1.2 and the current git master), and I am a bit confused about it -<br>
> > maybe someone here can enlighten me.<br>
> ><br>
> > In the manual (sec 4.1.4), kappa is defined by<br>
> ><br>
> > kappa^2 = (2*I*F^2)/(epsilon_0*epsilon_rf*R*T)<br>
> ><br>
> > where F is Faraday's constant, R is the ideal gas constant and T is the<br>
> > absolute temperature (this is the same as in Tironi et al.). In the<br>
> > source code (function calc_rffac in src/gromacs/mdlib), this is<br>
> > implemented as<br>
> ><br>
> > *kappa = std::sqrt(2*I/(EPSILON0*eps_rf*BOLTZ*Temp));<br>
> ><br>
> > now EPSILON0 is not the SI value of the dielectric constant, but<br>
> > defined (in src/gromacs/math/units.h) as<br>
> ><br>
> > #define EPSILON0<br>
> > (EPSILON0_SI*NANO*KILO)/(E_CHARGE*E_CHARGE*AVOGADRO)<br>
> ><br>
> > while<br>
> ><br>
> > #define BOLTZ (RGAS/KILO)<br>
> ><br>
> > substituting this into the source code, we get<br>
> ><br>
> > *kappa =<br>
> > std::sqrt(2*I*E_CHARGE*E_CHARGE*AVOGADRO*KILO/(EPSILON0_SI*NANO*KILO*eps_rf*RGAS*Temp));<br>
> ><br>
> > and since the Faraday constant is FARADAY = (E_CHARGE*AVOGADRO),<br>
> ><br>
> > this results in<br>
> ><br>
> > *kappa =<br>
> > std::sqrt(2*I*E_CHARGE*FARADAY/(EPSILON0_SI*NANO*eps_rf*RGAS*Temp));<br>
> ><br>
> > I have not substituted RGAS, as this seems indeed to be the ideal gas<br>
> > constant from the original formula (also, since<br>
> > RGAS=(BOLTZMANN*AVOGADRO), this would actually divide the other FARADAY<br>
> > into E_CHARGE).<br>
> ><br>
> > Comparing this with the original formula, the resulting value for kappa<br>
> > appears to be off by a factor of sqrt(AVOGADRO/NANO) or about 7.76e15.<br>
> > Suprisingly, this difference seems to result only in a small numerical<br>
> > change of the resulting potential, but unless I missed something, it<br>
> > still seems to be not correct.<br>
> ><br>
> > Regards,<br>
> > Jan<br>
> ><br>
> Isn't kappa printed to the log file if you use it in the mdp file?<br>
> How does it compare to your expectation?<br>
<br>
It is printed out, and it was the difference between this printed out<br>
value and those I calculated that prompted me to look at the source code<br>
in the first place.<br>
<br>
On a second look, I think I have it figured out now. The ionic strength<br>
here seems to be given in 1/nm^3 instead of mol/l. This would introduce<br>
another factor of "AVOGADRO" and the formula in the code works out<br>
again.<br>
<br>
Sorry for the confusion.<br>
<br>
Jan<br>
<br>
--<br>
Gromacs Developers mailing list<br>
<br>
* Please search the archive at <a href="http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List" rel="noreferrer" target="_blank">http://www.gromacs.org/Support/Mailing_Lists/GMX-developers_List</a> before posting!<br>
<br>
* Can't post? Read <a href="http://www.gromacs.org/Support/Mailing_Lists" rel="noreferrer" target="_blank">http://www.gromacs.org/Support/Mailing_Lists</a><br>
<br>
* For (un)subscribe requests visit<br>
<a href="https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers" rel="noreferrer" target="_blank">https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-developers</a> or send a mail to <a href="mailto:gmx-developers-request@gromacs.org" target="_blank">gmx-developers-request@gromacs.org</a>.<br>
</blockquote></div>