[gmx-users] Fwd: FW: PME settings in free energy calculations
Dries Van Rompaey
dries.vanrompaey at gmail.com
Tue Feb 9 09:50:36 CET 2016
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.
>
More information about the gromacs.org_gmx-users
mailing list