[gmx-users] Fwd: FW: PME settings in free energy calculations
Mark Abraham
mark.j.abraham at gmail.com
Tue Feb 9 09:20:29 CET 2016
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.
>
More information about the gromacs.org_gmx-users
mailing list