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" <[EMAIL PROTECTED]> > Subject: Re: [gmx-users] confused about rcoulomb<=rlist > To: [email protected] > Message-ID: > <[EMAIL PROTECTED]> > 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" <[EMAIL PROTECTED]> > > Subject: Re: [gmx-users] confused about rcoulomb<=rlist > > To: "Discussion list for GROMACS users" <[email protected]>, > > "Discussion list for GROMACS development" > <[EMAIL PROTECTED]> > > Cc: Michael Shirts <[EMAIL PROTECTED]>, John Chodera > > <[EMAIL PROTECTED]> > > Message-ID: > > <[EMAIL PROTECTED]> > > 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 [email protected] > > > 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 [EMAIL PROTECTED] > > > Can't post? Read http://www.gromacs.org/mailing_lists/users.php > > > > > > > > > ------------------------------ > > > > _______________________________________________ > > gmx-users mailing list > > [email protected] > > http://www.gromacs.org/mailman/listinfo/gmx-users > > > > > > End of gmx-users Digest, Vol 27, Issue 40 > > ***************************************** > > > > > ------------------------------ _______________________________________________ gmx-users mailing list [email protected] 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 [EMAIL PROTECTED] Can't post? Read http://www.gromacs.org/mailing_lists/users.php
