[gmx-developers] high-precision TIP3P parameters?
David van der Spoel
spoel at xray.bmc.uu.se
Sat May 27 14:00:37 CEST 2017
On 27/05/17 12:56, Justin Lemkul wrote:
>
>
> On 5/27/17 4:47 AM, Mark Abraham wrote:
>>
>>
>> On Sat, May 27, 2017 at 10:17 AM Igor Leontyev <ileontyev at ucdavis.edu
>> <mailto:ileontyev at ucdavis.edu>> wrote:
>>
>> Hi,
>> It is unlikely that inaccuracy in 4th digit of the parameters could
>> introduce significant error in simulation results.
>>
>> However, for the energy comparison and debugging which is often
>> made by
>> curious people (including myself) when switching between Amber and
>> Gromacs codes the inaccuracy is really bothering. Accepting suggested
>> correction will match Amber and Gromacs code energies and make the
>> comparison perfect (possibly for the cost of one time issue with
>> regression test).
>>
>>
>> True, but also makes GROMACS different from past versions of itself,
>> which also would worry people doing such a comparison when e.g.
>> contemplating upgrading their version. We're all for doing things
>> correctly, but that should mean everybody actually implementing the
>> published/accepted value, not agreeing with each other on how to be
>> different.
>>
>> It sounds like we're happy enough to change the values in GROMACS to
>> match the paper - can people please raise the issue/question with
>> others who distribute parameter sets?
>>
And by being compatible to Amber we would still differ from the
Jorgensen1983 paper.
We already have different TIP3P params for Amber, Charm and OPLS:
amber99sb-ildn.ff:OW 8 16.00 3.15061e-01 6.36386e-01
charmm27.ff:OWT3 8 15.999400 3.15058e-01 6.36386e-01
oplsaa.ff:opls_111 8 15.99940 3.15061e-01 6.36386e-01
It would probably make sense to update each of the ffnonbonded.itp files
to match the actual force field. If we would like to be really
compatible we would need to correct for the physical constants as well.
>
> From the CHARMM perspective, we just have a standard format that we
> adhere to in terms of the precision of values. So epsilon is -0.1521
> instead of the more precise -0.15207, but I don't see that changing on
> our end. That introduces the possibility that GROMACS and AMBER match
> exactly but GROMACS and CHARMM do not, but that already suggests CHARMM
> and AMBER do not match, a comparison I have never done (though GROMACS
> and CHARMM do match quite nicely).
>
> Another issue is that, when last I checked, all programs (AMBER, CHARMM,
> NAMD, GROMACS) use a different value of 1/(4*pi*eps0) for calculating
> electrostatics so there is always some error there, irrespective of the
> actual topology. Again, these are quantities that will not substantially
> change dynamics (if at all), but it seems as though we're chasing
> minutiae when there are other things that aren't ever going to be standard.
>
> -Justin
>
>> Mark
>>
>> Igor
>>
>> >
>> > Hi,
>> >
>> > These difference are irrelevant for any practical purposes,
>> since the
>> > model (and sampling) is not by far that accurate. But it would
>> still be
>> > nice to have all decimals that we have in files correct.
>> > Is this somehow related to the differences in the kcal/mol to
>> kJ/mol
>> > conversion in different packages and/or at different times?
>> >
>> > Cheers,
>> >
>> > Berk
>> >
>> > On 24/05/17 21:18 , Igor Leontyev wrote:
>> >> Hi all. The issue was discussed in past but nothing was corrected
>> >> since then.
>> >>
>> >>
>>
>> https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-developers/2011-April/005178.html
>>
>> >>
>> >>
>> >> To make gromacs TIP3P energies matching Amber code energies:
>> >>> Suggested correction in amberXX.ff/ffnonbonded.itp:
>> >>> ;OW 8 16.00 0.0000 A 3.15061e-01
>> 6.36386e-01
>> >>> OW 8 16.00 0.0000 A 3.15076e-01 6.35968e-01
>> >>
>> >> Igor
>> >>
>> >>> So, going back to the original TIP3P paper, the parameters
>> are listed
>> >> as A and C parameters (in the formula a/r^12 ? b/r^6).
>> >>>
>> >>> The parameters listed there are A = 582000 A^12 kcal/mol and
>> C = 595
>> >>> A^6 kcal/mol.
>> >>>
>> >>> Converting to sigma/epsilon form, this gives:
>> >>> Sigma = 0.31506561105 nm
>> >>> Epsilon = 0.6362717354 kJ/mol
>> >>>
>> >>> This compares with the values in the most of the distributed ITP
>> >>> files of
>> >>>
>> >>> Sigma = 0.315061 nm (a difference of 4.6 x 10^-6 nm, or a
>> relative
>> >>> difference of of 1.4 x 10^-5)
>> >>> Epsilon = 0.636386 kJ/mol (a difference of 1.1 x 10^4
>> kJ/mol, or a
>> >>> relative difference of 1.8 x 10^-4)
>> >>> (note the GROMACS one is also
>> 0.152100 kcal/mol)
>> >>> I don?t know how much of a difference this would make (likely a
>> >>> difference of the scale of).
>> >>>
>> >>> Now, that?s in the original paper. What values are ACTUALLY
>> >>> specified by other codes?
>> >>>
>> >>> For AMBER, (as produced by tleap), it outputs.
>> >>>
>> >>> A = 581935.564 A^12 kcal/mol
>> >>> C = 594.825035 A^6 kcal/mol
>> >>>
>> >>> Which leads to:
>> >>>
>> >>> Sigma = 0.315075 (a difference of 3.0 x 10^-5 nm, or a relative
>> >>> difference of 3.0 x 10^-5)
>> >>> Epsilon = 0.635968 kJ/mol (a difference of 9.4 x 10^4 kJ/mol
>> or a
>> >>> relative difference of 4.7 x 10^-4)
>> >>> (which is also 0.152000 kcal/mol)
>> >>>
>> >>> So ? I don?t know. Nobody seems to be getting it _exactly_
>> right,
>> >>> and it?s also not so clear how much of a difference it
>> makes. One
>> >>> can estimate that enthalpies would have a relative error of
>> around
>> >>> 10^-4, and relative errors in the density of 3*dsigma = 10^-4
>> to 10^-5.
>> >>>
>> >>> So I don?t know that there is anything actionable to take. Just
>> >>> thought it was interesting to bring out. Something to ponder
>> for now.
>> >>
>>
>> >
>> > I?d say ?probably? irrelevant for practical purposes. It
>> wouldn?t take
>> too much simulation time to determine the density sensitively
>> enough to tell
>> the difference. But, pretty much irrelevant.
>> >
>> > ? Is this somehow related to the differences in the
>> kcal/mol to kJ/mol
>> > ? conversion in different packages and/or at different times?
>> >
>> > I think that it?s due to differences in rounding from the
>> original A and
>> C values vs. rounding from conversion of them to epsilon and sigma
>> values.
>> The epsilons in AMBER and GROMACS are clearly equal to using
>> either 0.152100
>> or 0.15200 kcal/mol as the value (the parameters in the original
>> paper
>> correspond to 0.151207).
>> >
>> > I think there might have been a kJ/kcal issue at some point,
>> but those
>> were fixed a decade or so ago (I seem to recall that popping up
>> when I was a
>> graduate student, but not since then).
>> >
>> > ? These difference are irrelevant for any practical
>> purposes, since the
>> > ? model (and sampling) is not by far that accurate. But it
>> would still be
>> > ? nice to have all decimals that we have in files correct.
>> >
>> > At this stage, it probably involves getting some general set
>> everyone
>> together and deciding in what the ?right? values are. I could
>> potentially do
>> that . . . but as Berk points out, not that big a deal. Unless
>> one is
>> trying to match energy snapshots, and getting frustrated.
>> >
>> > For what it?s worth, the Desmond viparr utility uses for TIP3P:
>> {"sigma":
>> 3.15061, "epsilon": 0.1521}. Both match GROMACS .itp libraries,
>> but not
>> tleap or the original paper.
>> >
>> > Best,
>> > ~~~~~~~~~~~~~~~~
>> > Michael Shirts
>> > Associate Professor
>> > michael.shirts at colorado.edu <mailto:michael.shirts at colorado.edu>
>> > http://www.colorado.edu/lab/shirtsgroup/
>> > Phone: (303) 735-7860
>> > Office: JSCBB C123
>> > Department of Chemical and Biological Engineering
>> > University of Colorado Boulder
>> >
>> > On 5/25/17, 1:51 PM,
>> "gromacs.org_gmx-developers-bounces at maillist.sys.kth.se
>> <mailto:gromacs.org_gmx-developers-bounces at maillist.sys.kth.se> on
>> behalf of
>> Berk Hess" <gromacs.org_gmx-developers-bounces at maillist.sys.kth.se
>> <mailto:gromacs.org_gmx-developers-bounces at maillist.sys.kth.se> on
>> behalf of
>> hess at kth.se <mailto:hess at kth.se>> wrote:
>> >
>> > Hi,
>> >
>> > These difference are irrelevant for any practical purposes,
>> since the
>> > model (and sampling) is not by far that accurate. But it
>> would still be
>> > nice to have all decimals that we have in files correct.
>> > Is this somehow related to the differences in the kcal/mol to
>> kJ/mol
>> > conversion in different packages and/or at different times?
>> >
>> > Cheers,
>> >
>> > Berk
>> >
>> > On 24/05/17 21:18 , Igor Leontyev wrote:
>> > > Hi all. The issue was discussed in past but nothing was
>> corrected
>> > > since then.
>> > >
>> > >
>>
>> https://mailman-1.sys.kth.se/pipermail/gromacs.org_gmx-developers/2011-April/005178.html
>>
>> > >
>> > >
>> > > To make gromacs TIP3P energies matching Amber code energies:
>> > > > Suggested correction in amberXX.ff/ffnonbonded.itp:
>> > > > ;OW 8 16.00 0.0000 A 3.15061e-01
>> 6.36386e-01
>> > > > OW 8 16.00 0.0000 A 3.15076e-01
>> 6.35968e-01
>> > >
>> > > Igor
>> > >
>> > >> So, going back to the original TIP3P paper, the
>> parameters are listed
>> > > as A and C parameters (in the formula a/r^12 ? b/r^6).
>> > >>
>> > >> The parameters listed there are A = 582000 A^12 kcal/mol
>> and C = 595
>> > >> A^6 kcal/mol.
>> > >>
>> > >> Converting to sigma/epsilon form, this gives:
>> > >> Sigma = 0.31506561105 nm
>> > >> Epsilon = 0.6362717354 kJ/mol
>> > >>
>> > >> This compares with the values in the most of the
>> distributed ITP
>> > >> files of
>> > >>
>> > >> Sigma = 0.315061 nm (a difference of 4.6 x 10^-6 nm, or
>> a relative
>> > >> difference of of 1.4 x 10^-5)
>> > >> Epsilon = 0.636386 kJ/mol (a difference of 1.1 x 10^4
>> kJ/mol, or a
>> > >> relative difference of 1.8 x 10^-4)
>> > >> (note the GROMACS one is also
>> 0.152100
>> kcal/mol)
>> > >> I don?t know how much of a difference this would make
>> (likely a
>> > >> difference of the scale of).
>> > >>
>> > >> Now, that?s in the original paper. What values are
>> ACTUALLY
>> > >> specified by other codes?
>> > >>
>> > >> For AMBER, (as produced by tleap), it outputs.
>> > >>
>> > >> A = 581935.564 A^12 kcal/mol
>> > >> C = 594.825035 A^6 kcal/mol
>> > >>
>> > >> Which leads to:
>> > >>
>> > >> Sigma = 0.315075 (a difference of 3.0 x 10^-5 nm, or a
>> relative
>> > >> difference of 3.0 x 10^-5)
>> > >> Epsilon = 0.635968 kJ/mol (a difference of 9.4 x 10^4
>> kJ/mol or a
>> > >> relative difference of 4.7 x 10^-4)
>> > >> (which is also 0.152000 kcal/mol)
>> > >>
>> > >> So ? I don?t know. Nobody seems to be getting it
>> _exactly_ right,
>> > >> and it?s also not so clear how much of a difference it
>> makes. One
>> > >> can estimate that enthalpies would have a relative error
>> of around
>> > >> 10^-4, and relative errors in the density of 3*dsigma =
>> 10^-4 to
>> 10^-5.
>> > >>
>> > >> So I don?t know that there is anything actionable to
>> take. Just
>> > >> thought it was interesting to bring out. Something to
>> ponder for now.
>> --
>> 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
>> <mailto:gmx-developers-request at gromacs.org>.
>>
>>
>>
>
--
David van der Spoel, Ph.D., Professor of Biology
Head of Department, Cell & Molecular Biology, Uppsala University.
Box 596, SE-75124 Uppsala, Sweden. Phone: +46184714205.
http://www.icm.uu.se
More information about the gromacs.org_gmx-developers
mailing list