[gmx-users] Fwd: FW: PME settings in free energy calculations
Mark Abraham
mark.j.abraham at gmail.com
Tue Feb 9 11:06:13 CET 2016
Hi,
Runs 2 and 3 are using different rvdw and rvdw-switch, which doesn't seem
like the comparison you should want to make. However that should not affect
the Coulomb calculation at all - and certainly not to this degree. What do
you see when it matches? I notice that somehow you are triggering the use
of the generic nonbonded kernel - is that deliberate by you?
The group scheme plus rvdw != rcoulomb plus a switch could well be a
combination that got broken at some point. I tried to verify such
equivalence before we released 5.0, but I could well have missed some. If
so, the rerun with 4.5.7 would be instructive. (Michael's advice on
alchemistry.org is largely based on 4.5-era GROMACS).
Mark
On Tue, Feb 9, 2016 at 9:51 AM Dries Van Rompaey <dries.vanrompaey at gmail.com>
wrote:
> Hi Mark,
>
> Thanks for the clarification. I did use the rerun option as described in
> the link you provided.
> I uploaded my log files here:
> https://www.dropbox.com/sh/tcmqhn5fszeis27/AABcYOI6_FrjG7qoYJycf_35a?dl=0
>
> Thanks for your time,
> Kind regards
>
> Dries
>
> On 9 February 2016 at 09:20, Mark Abraham <mark.j.abraham at gmail.com>
> wrote:
>
> > Hi,
> >
> > On Fri, Feb 5, 2016 at 1:32 PM Dries Van Rompaey <
> > dries.vanrompaey at gmail.com>
> > wrote:
> >
> > > Hi,
> > >
> > > The corresponding mdp output section is:
> > >
> > > ; OPTIONS FOR ELECTROSTATICS AND VDW
> > > ; Method for doing electrostatics
> > > coulombtype = PME
> > > coulomb-modifier = Potential-shift-Verlet
> > > rcoulomb-switch = 0
> > > rcoulomb = 1.0
> > > ; Relative dielectric constant for the medium and the reaction field
> > > epsilon-r = 1
> > > epsilon-rf = 0
> > > ; Method for doing Van der Waals
> > > vdwtype = cutoff
> > > vdw-modifier = Potential-switch
> > > ; cut-off lengths
> > > rvdw-switch = 0.9
> > > rvdw = 1.0
> > > ; Apply long range dispersion corrections for Energy and Pressure
> > > DispCorr = EnerPres
> > > ; Extension of the potential lookup tables beyond the cut-off
> > > table-extension = 1
> > > ; Separate tables between energy group pairs
> > > energygrp-table =
> > > ; Spacing for the PME/PPPM FFT grid
> > > fourierspacing = 0.08
> > > ; FFT grid size, when a value is 0 fourierspacing will be used
> > > fourier-nx = 0
> > > fourier-ny = 0
> > > fourier-nz = 0
> > > ; EWALD/PME/PPPM parameters
> > > pme_order = 6
> > > ewald_rtol = 1e-06
> > > ewald-rtol-lj = 0.001
> > > lj-pme-comb-rule = Geometric
> > > ewald-geometry = 3d
> > > epsilon_surface = 0
> >
> >
> > That all seems reasonable.
> >
> >
> > > I evaluated the energy with:
> > >
> > > Verlet - settings above (with standard verlet-buffer-tolerance)
> > > @ s0 legend "Coulomb-14"
> > > @ s1 legend "Coulomb (SR)"
> > > @ s2 legend "Coul. recip."
> > > 0.000000 -2221.295898 -546205.250000 9271.518555
> > > Total - 539155.0273
> > >
> > > Verlet - settings above - with verlet-buffer-tolerance -1
> > > 0.000000 -2221.295898 -546205.250000 9271.516602
> > > Total -539155.029296
> > >
> >
> > That kind of agreement is expected, but perhaps might agree exactly if
> you
> > were following
> > http://www.gromacs.org/Documentation/How-tos/Single-Point_Energy
> >
> > Group - settings above
> > > @ s0 legend "Coulomb-14"
> > > @ s1 legend "Coulomb (SR)"
> > > @ s2 legend "Coul. recip."
> > > 0.000000 -2221.296387 -487993.593750 -48931.437500
> > > Total -539146.3276
> > >
> >
> > The difference in those components is expected - in PME you have to avoid
> > counting contributions in the long-range part from pairs that are
> excluded
> > e.g. from the bonded topology. The two schemes handle this in different
> > places. Their sum should match (within the limits of floating-point
> > arithmetic), and a relative error under 1e-5 for summing a few million
> > numbers is acceptable for judging the equivalence of the implementations.
> >
> > Verlet - 3 nm cutoffs, 2.9 vdwswitch, coulomb-modifier none
> > > @ s0 legend "Coulomb-14"
> > > @ s1 legend "Coulomb (SR)"
> > > 0.000000 -2221.295898 -774817.875000
> > > Total -777039.1709
> > >
> > >
> > > Group - 3 nm cutoffs, 2.9 vdwswitch, coulomb-modifier none
> > > @ s0 legend "Coulomb-14"
> > > @ s1 legend "Coulomb (SR)"
> > > 0.000000 -2221.296387 -553398.125000
> > > Total -555619.4214
> > >
> >
> > That ought to agree much better, but It's hard to speculate on where the
> > problem lies without seeing log files. Please upload some to a file
> sharing
> > service and post the links here (list doesn't take attachments)
> >
> > Mark
> >
> >
> > > Thanks,
> > >
> > > Dries
> > >
> > > On 4 February 2016 at 14:48, Michael Shirts <mrshirts at gmail.com>
> wrote:
> > >
> > > > Also, can you be more precise about the inconsistency in the outputs
> > > > you are seeing (with numerical output from GROMACS and analysis
> code).
> > > >
> > > > On Thu, Feb 4, 2016 at 6:17 AM, Michael Shirts <mrshirts at gmail.com>
> > > wrote:
> > > > > Hi, Dries-
> > > > >
> > > > > Can you print out what the corresponding section of mdout.mdp look
> > > like?
> > > > >
> > > > > Have you tried using the cutoffs above with group cutoffs, and
> seen
> > > > > what the difference is?
> > > > >
> > > > > On Thu, Feb 4, 2016 at 4:59 AM, Dries Van Rompaey
> > > > > <dries.vanrompaey at gmail.com> wrote:
> > > > >> Hi,
> > > > >>
> > > > >> I'm using PME with the default Potential-Shift-Verlet modifier. If
> > I'm
> > > > >> interpreting the manual and comments on the mailinglist correctly,
> > the
> > > > >> settings for rcoulomb-switch shouldn't affect this and
> PME-switching
> > > > >> shouldn't be necessary (although there seems to have been some
> > > > discussion
> > > > >> about this, https://gerrit.gromacs.org/#/c/2220/). Has a
> consensus
> > > been
> > > > >> reached elsewhere?
> > > > >>
> > > > >> Current nonbonded settings in my .mdp files are:
> > > > >> ; Electrostatics
> > > > >> coulombtype = PME
> > > > >> rcoulomb = 1.0
> > > > >>
> > > > >> ; van der Waals
> > > > >> vdwtype = cutoff
> > > > >> vdw-modifier = Potential-switch
> > > > >> rvdw-switch = 0.9
> > > > >> rvdw = 1.0
> > > > >>
> > > > >> ; Apply long range dispersion corrections for Energy and Pressure
> > > > >> DispCorr = EnerPres
> > > > >>
> > > > >> ; Spacing for the PME/PPPM FFT grid
> > > > >> fourierspacing = 0.08
> > > > >>
> > > > >> ; EWALD/PME/PPPM parameters
> > > > >> pme_order = 6
> > > > >> ewald_rtol = 1e-06
> > > > >> epsilon_surface = 0
> > > > >>
> > > > >> Thanks in advance,
> > > > >>
> > > > >> Dries
> > > > >>
> > > > >> On 4 February 2016 at 00:56, Michael Shirts <mrshirts at gmail.com>
> > > wrote:
> > > > >>
> > > > >>> Hi, Dries-
> > > > >>>
> > > > >>> Questions like this are probably best answered on the gmx-users
> > list.
> > > > >>> I can't say too much for the Verlet scheme -- I know that it was
> > > > >>> relatively recently adapted for free energies, and there may be
> > some
> > > > >>> combinations of settings that could give unanticipated results.
> > > > >>>
> > > > >>> Pretty much all of our experience about nonbonded calculations
> and
> > > > >>> free energies is collected in the following paper:
> > > > >>> http://pubs.acs.org/doi/abs/10.1021/ct4005068. We used the group
> > > > >>> cutoff scheme since that was the only one that was supported in
> > > > >>> GROMACS at the time.
> > > > >>>
> > > > >>> Since you haven't sent any files, it's hard to tell what is
> > actually
> > > > >>> going on. The one thing that has a tendency to happen with the
> more
> > > > >>> recent update schemes is if you set a potential modifier that is
> a
> > > > >>> switch, but don't set the distance the switch starts at, then it
> is
> > > > >>> automatically set to zero. Check the mdout.mdp to see if this is
> > > > >>> happening.
> > > > >>>
> > > > >>> > From: Van Rompaey Dries <Dries.VanRompaey at uantwerpen.be>
> > > > >>> Date: Wednesday, February 3, 2016 at 12:34 PM
> > > > >>> To: Michael Shirts <michael.shirts at colorado.edu>
> > > > >>> Subject: PME settings in free energy calculations
> > > > >>>
> > > > >>> Dear prof. Shirts,
> > > > >>>
> > > > >>> I'm currently trying to figure out the PME settings for a
> relative
> > > > >>> free binding energy simulation I'm working on. I took the
> > parameters
> > > > >>> from Matteo Aldeghi's
> > > > >>> paper(
> > > > >>>
> > > >
> > >
> >
> http://pubs.rsc.org/en/content/articlelanding/2015/sc/c5sc02678d#!divAbstract
> > > > >>> )
> > > > >>> as a starting point (adapted for the Verlet scheme by scaling the
> > > > >>> fourierspacing to 0.08 and setting coulomb = rvdw = 1.0). I then
> > > tried
> > > > >>> to verify these settings by comparing a single point energy
> > > > >>> calculation with these settings and one with very long coulomb
> > > cutoffs
> > > > >>> as recommended on alchemistry.org.
> > > > >>>
> > > > >>> Unfortunately, I can't seem to get this quite right. I'm getting
> > > > >>> differences in the hundreds of kj/mol, leading me to suspect I'm
> > > doing
> > > > >>> something wrong. I'm calculating the energy values by extracting
> > > > >>> CoulombSR and Coul-recip from the energy.xvg files. I've tried
> > > > >>> calculating the energy with coulomb cutoffs at 3 nm and 10 nm,
> but
> > > > >>> agreement with the PME results remains rather poor. David Mobley
> > > > >>> mentioned you performed extensive research on this topic, and I'm
> > > > >>> hoping you could point me in the right direction.
> > > > >>>
> > > > >>> Thanks in advance,
> > > > >>>
> > > > >>> Dries
> > > > >>>
> > > > >>>
> > > > >>>
> > > > >>> ~~~~~~~~~~~~~~~~
> > > > >>> Michael Shirts
> > > > >>> Associate Professor
> > > > >>> michael.shirts at colorado.edu
> > > > >>> Phone: (303) 735-7860
> > > > >>> Office: JSCBB D317
> > > > >>> Department of Chemical and Biological Engineering
> > > > >>> University of Colorado Boulder
> > > > >>> --
> > > > >>> Gromacs Users mailing list
> > > > >>>
> > > > >>> * Please search the archive at
> > > > >>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_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-users
> > > or
> > > > >>> send a mail to gmx-users-request at gromacs.org.
> > > > >>>
> > > > >> --
> > > > >> Gromacs Users mailing list
> > > > >>
> > > > >> * Please search the archive at
> > > > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_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-users
> > or
> > > > send a mail to gmx-users-request at gromacs.org.
> > > > --
> > > > Gromacs Users mailing list
> > > >
> > > > * Please search the archive at
> > > > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_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-users
> or
> > > > send a mail to gmx-users-request at gromacs.org.
> > > >
> > > --
> > > Gromacs Users mailing list
> > >
> > > * Please search the archive at
> > > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_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-users or
> > > send a mail to gmx-users-request at gromacs.org.
> > >
> > --
> > Gromacs Users mailing list
> >
> > * Please search the archive at
> > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_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-users or
> > send a mail to gmx-users-request at gromacs.org.
> >
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_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-users or
> send a mail to gmx-users-request at gromacs.org.
>
More information about the gromacs.org_gmx-users
mailing list