[gmx-users] Fwd: FW: PME settings in free energy calculations
Dries Van Rompaey
dries.vanrompaey at gmail.com
Tue Feb 9 13:30:26 CET 2016
Hi,
I'm not deliberately triggering the generic kernel. I unfortunately have no
idea what's causing that. I uploaded logfiles for a setup with matching
cutoffs as setup3alt and setup2alt (I got quite a lot of warnings for this
setup, which is why I initially changed the cutoffs slightly).
I'm a bit baffled what could cause the large difference between the verlet
setup with long cutoffs and the verlet setup with shorter cutoffs and PME..
After all, that's the one that should be fairly close in order to prove my
electrostatics treatment is appropriate. This difference seems to be a lot
more pronounced than the difference to the group scheme.
Thanks,
Dries
On 9 February 2016 at 11:06, Mark Abraham <mark.j.abraham at gmail.com> wrote:
> 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.
> >
> --
> 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