[gmx-users] Re: confused about rcoulomb<=rlist

[Janne Hirvi]

Hello!

I was delighted to see this conversation about the use of rcoulomb and rlist. I
made my own three test simulations with bulk water (1372 SPC/E molecules & 2fs
& 200ps) using double precision compilation with parameters

nstlist       = 5
vdwtype       = Switch
rvdw_switch   = 1.0
rvdw          = 1.1
rlist         = 1.2
coulombtype   = Pme
rcoulomb      = a)1.0        b)1.2        c)1.2
ewald_rtol    = a)1.0E-6     b)1.0E-6     c)3.632E-9

and found out that energy conservation is actually very good (drift~4kJ/mol) in
cases a) and c), but poor in case b) (~209kJ/mol) as you predicted. However I
was surprised that with larger rcoulomb one should use that much smaller
ewald_rtol to see energy conservation (I had made earlier some tests with
little bit larger value)?

Actually earlier I have used rcoulomb=1.2 and ewald_rtol=1.0E-5 with pure wan
der Waals cut-off (rlist=rvdw=1.2) and of course the energy conservation has
been bad. I would need some comments from the experts like you about these
simulations:

1. Do you see a "huge" problem when I havent used perfect PME parameters with
pure wan der Waals cut-off, while the vdw cut-off already destroyes the energy
conservation? 

2. How about the energy conservation in general with temperature coupling. Wan
der Waals cut-off and not perfect PME parameters will create clear drift in NVE
simulations but is it acceptable to remove it with temperature coupling?


Thanks,

Janne

> Message: 1
> Date: Wed, 12 Jul 2006 12:25:05 -0400
> From: "Michael Shirts" <mrshirts at gmail.com>
> Subject: Re: [gmx-users] confused about rcoulomb<=rlist
> To: gmx-users at gromacs.org
> Message-ID:
> 	<662e7a930607120925i1eb6128cr3c8642722b0d033a at mail.gmail.com>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
> 
> Hi, all-
> 
> I admit I've been wondering about this as well -- I was involved in
> some of the previous discussions.  It seems like the simplest thing is
> to have the real space cutoff be at rcoulomb, and if one wants the
> real space cutoff to go out to rlist, then one sets rcoulomb = rlist.
> But one is not forced to.    Am I missing something?
> 
> I think part of the previous discussion was how to get the best energy
> conservation -- Berk argued (correctly, I think -- at least it appears
> in agreement with simulations I ran) that setting the real space
> cutoff = rlist gives the best energy conservation.  But this is
> independent of the question of user control -- I think the user should
> be able to choose real space cutoff < rlist, for the reasons David
> Mobley lists above.
> 
> Best,
> Michael Shirts
> Research Fellow
> Chemistry Department
> Columbia University
> 
> > Message: 8
> > Date: Wed, 12 Jul 2006 09:08:01 -0700
> > From: "David Mobley" <dmobley at gmail.com>
> > Subject: Re: [gmx-users] confused about rcoulomb<=rlist
> > To: "Discussion list for GROMACS users" <gmx-users at gromacs.org>,
> >         "Discussion list for GROMACS development"
> <gmx-developers at gromacs.org>
> > Cc: Michael Shirts <mrshirts at gmail.com>, John Chodera
> >         <jchodera at gmail.com>
> > Message-ID:
> >         <bc2c99750607120908k684349ebqe904ff21ba3d256c at mail.gmail.com>
> > Content-Type: text/plain; charset=ISO-8859-1; format=flowed
> >
> > Berk,
> >
> > > I have added the check.
> > >
> > > The problem was that Gromacs did not truncate the potential
> > > at rcoulomb with PME, only rlist was used. However rcoulomb
> > > was used in the determination of beta.
> > > To avoid that people thought that rcoulomb had an effect
> > > on the cut-off, I have added the check.
> >
> > OK, so, in other words, rcoulomb now only affects beta with PME; rlist
> > is used for the electrostatic cutoff.
> >
> > Is there a straightforward way this could be changed in future
> > versions so that the cutoff used for electrostatics could be different
> > from rlist? In particular, I am concerned that this means that one
> > can't vary the vdw cutoff independently of the electrostatics cutoff
> > beyond some limit. That is, rlist needs to be >rvdw with shift/switch
> > vdw potentials:
> >
> > WARNING 1 [file equil_constp0.7.mdp, line unknown]:
> >   For energy conservation with switch/shift potentials, rlist should be
> 0.1
> >   to 0.3 nm larger than rcoulomb/rvdw.
> >
> > Suppose I want to run with a 1.4 nm rvdw -- then I now have to use
> > rlist=1.5 nm and rcoulomb=1.5 nm? That doesn't seem good. One would
> > like to be able to change rvdw independently of the electrostatics --
> > for example, to explore the effects of long range van der Waals
> > cutoffs. In fact, Michael Shirts and I and others had been trying to
> > do exactly that in GROMACS 3.1.4 and 3.3 (adjust rvdw and rlist
> > independently of rcoulomb in binding free energy calculations, in
> > order to see how strongly the choice of vdW cutoff affects computed
> > binding free energies). If changing rlist also changes the real space
> > cutoff for electrostatics, things are substantially more complicated
> > and it seems unclear how to even address this issue.
> >
> > Maybe a solution would be to NOT use rlist for beta and for the PME
> > potential, and instead to use rcoulomb or some extra variable... Then
> > one could change the neighbor list (rlist) and rvdw without having to
> > tweak the electrostatic interactions...
> >
> > > What is the optimal cut-off scheme is a different issue.
> > > Indeed one would always want the force to go smoothly
> > > at the cut-off, or before the cut-off in case one has charge groups
> > > or nstlist>1.
> > > However for PME one can not have 'exact' electrostatics while
> > > the particle-particle force is zero at the cut-off since the reciprocal
> > > space requires an error function contribution in real space.
> > > Therefore the real space interaction has infinite range, but decays
> > > very rapidly.
> > > In PME the real space interaction is erfc(beta r)/r, which in Gromacs
> > > is applied to all atom pairs in the neighborlist. The cut-off error
> > > made here is very small, since ewald_rtol=1e-5 (I also often use 1e-6).
> > > The alternative would be to somehow make the direct space interaction
> > > go exactly to zero before the cut-off. This would lead to larger
> > > errors, as one then misses more of the electrostatic interactions.
> > > The only advantage would be that the integration could be more
> > > reversible (assuming that all other algorithms are well reversible).
> > >
> > > Given the accuracies of all other algorithms I would say it does not
> > > make sense to remove some electrostatic interactions at the cut-off
> > > when they are already 1e-5 or 1-e6 smaller than the Coulomb
> > > interaction at that distance.
> >
> > So, if I understand properly, you're saying that the notion of an
> > electrostatics "cutoff" doesn't really make sense with PME, so my
> > logic about having the neighborlist extend past the "cutoff" doesn't
> > really apply. Rather, one just wants to ensure that the "cutoff" is
> > sufficiently large that electrostatic interactions are fairly small
> > and thus any inaccuracies due to neighbor list issues will be
> > basically negligible...?
> >
> > Thanks,
> > David
> >
> >
> > > Berk.
> > >
> > >
> > > _______________________________________________
> > > gmx-users mailing list    gmx-users at gromacs.org
> > > http://www.gromacs.org/mailman/listinfo/gmx-users
> > > Please don't post (un)subscribe requests to the list. Use the
> > > www interface or send it to gmx-users-request at gromacs.org.
> > > Can't post? Read http://www.gromacs.org/mailing_lists/users.php
> > >
> >
> >
> > ------------------------------
> >
> > _______________________________________________
> > gmx-users mailing list
> > gmx-users at gromacs.org
> > http://www.gromacs.org/mailman/listinfo/gmx-users
> >
> >
> > End of gmx-users Digest, Vol 27, Issue 40
> > *****************************************
> >
> 
> 
> ------------------------------