[gmx-users] cutt offs

[baptista at itqb.unl.pt]

Dear Ilya, David and Eric,

I don't see any thermodynamic inconsistency in using different cutoffs for
different molecules in the system. From a theoretical viewpoint, choosing
a cutoff is just part of specifying the Hamiltonian of the system, which
is what really _defines_ your system. Thus, if you try out different
cutoffs you are actually simulating different systems, whose thermodynamic
properties would differ from those of the real system in different ways.
Which one should be the best? It may be tempting to say that the best
thing would be to run the simulations using the largest possible cutoff,
since it will correspond to the least "approximated" energy, but this is
not quite so. Each force field provides an _effective_ energy function
that is parameterized with a fixed set of cutoffs (or at least it
should...), and thus these should be the ones used in any subsequent
simulation, unless one has very strong reasons to go beyond them (eg, when
you suspect that the force field cannot handle properly your particular
system, as in Ilya's example of a very large cavity). In any case, the use
of different cutoffs for different molecules (eg, protein and solvent)
will not send you to "thermodynamic neverland". You stay in the ensemble
implied by your choice of thermal and pressure baths, regardless of your
choice of cutoffs. That is, changing the Hamiltonian changes the nature of
the system, but not the thermodynamics.

I don't think that David's example of the two halves of a box having
different cutoffs is a particular illustration of the issue raised by
Ilya. Here you don't have different cutoffs for different _molecules_, but
rather for different _regions_. Thus, whatever manipulation you do to
your Hamiltonian, you will never be able to write it down in a way that is
independent of the reference frame. It is not possible to assign the
different cutoffs to different "pseudo-molecules", even if the latter are
artificial constructs made up from collecting assorted Hamiltonian terms.
The only formal way of making physical sense of this situation is to
assume that there is some kind of _external_ field that imposes the
different cutoffs in each half. Obviously, the nature of this "field" is
not electrostatic, nor gravitational, nor of any other usual type. This
field is just a formal construct that allows us to think of the situation
in physical terms. The consequence of introducing this field is that the
equilibrium will be attained between both halves when the _total_ chemical
potential is the same, as implied by thermodynamics. For example, if we
imagine for a moment that the field was electrostatic or gravitational,
the two halves would have the same _electrochemical_ or _gravitochemical_
potential. In that sense, I agree with David that the _strictly_chemical_
potentials need not be the same in both halves, but the _total_ chemical
potentials must certainly be identical. (Note that the splitting of the
total chemical potential can be highly artificial, as in the
electrochemical case; see Guggenheim or Denbigh on this.) Therefore, I
don't see any problem in defining the ensemble being modeled: eg, if you
are using thermal and pressure baths, the ensemble being modeled would be
the isothermal-isobaric _with_ the applied field (traditionally, external
fields are not included in the designation of the ensemble, but rather
added as an additional external "parameter"). Of course, the situation is
very strange in terms of the actual simulation, with atoms experiencing a
sudden change of interactions when they move from one half to the other.
In particular, I would expect some instabilities, which could even result
in a total crash of the numerical solution of the equations of motion.
However, even in this case (which is _not_ the one discussed by Ilya), I
don't see any theoretical problem. Again: you don't screw up
thermodynamics just by messing up the Hamiltonian; you just get a
different system.

As for the range of different interactions, I think that Ilya's point is
right, in the sense that we _could_ use a lower cutoff for dipole-dipole
interactions than for charge-dipole ones and still get comparable errors,
because the former falls off faster than the latter. Actually, the now
standard procedure of using group-based electrostatic cutoffs came from
the fact that most charge groups can be arranged to be neutral, and
therefore most of the interactions in the system become dipole-dipole,
some charge-dipole, and only a few charge-charge. By arranging things in
this way the truncation error is somewhat minimized.

And if you really want to make a simulation with such a big cavity, I
don't think you have much choice: since there is no danger of
thermodynamic mess up, be bold and make a hack to increase the cutoff of
your LJ particle! Then check if the small-cavity case still gives the same
properties as before. If so, I guess you could proceed (at your own
risk!). Of course, you should do all this with your fingers crossed! ;-)

Best,
Antonio

-- 
Antonio M. Baptista
Instituto de Tecnologia Quimica e Biologica, Universidade Nova de Lisboa
Av. da Republica, EAN, ITQB II, Piso 6, Apartado 127
2781-901 Oeiras, Portugal
phone: +351-214469619(NEW!)   email: baptista at itqb.unl.pt
fax:   +351-214411277         WWW:   http://www.itqb.unl.pt/~baptista
--------------------------------------------------------------------------


On Mon, 29 Dec 2003, Ilya Chorny wrote:

> The reason I was asking was not to debate the subject. I am growing a LJ
> cavity in water. When the cavity gets really big beyond the standard LJ
> cutoff I run into problems as one might imagine. I was looking for a way to
> overcome this problem.
>
> By the way. A dipole-dipole interaction falls of faster than a chrage-charge
> or charge-dipole.
>
> Thanks for your help!!!!
>
> Ilya
>
>
> -----Original Message-----
> From: gmx-users-admin at gromacs.org [mailto:gmx-users-admin at gromacs.org]On
> Behalf Of Eric Jakobsson
> Sent: Monday, December 29, 2003 4:01 PM
> To: gmx-users at gromacs.org
> Subject: RE: [gmx-users] cutt offs
>
>
> Absolutely not.
>
> Probably too big a gap in perception here to profitably continue such a
> debate.
>
> At 01:20 PM 12/29/2003 -0800, you wrote:
> >Your example is extreme. What about the case of a single charge in
> >water. The water-water Coulomb cut off should be less than the charge
> >water Coulomb cutoff.
> >
> >Ilya
> >
> >
> >-----Original Message-----
> >From: gmx-users-admin at gromacs.org [mailto:gmx-users-admin at gromacs.org]
> >On Behalf Of David van der Spoel
> >Sent: Monday, December 29, 2003 12:40 PM
> >To: gmx-users at gromacs.org
> >Subject: RE: [gmx-users] cutt offs
> >
> >On Mon, 29 Dec 2003, Ilya Chorny wrote:
> >
> > >Its not obvious to me why the thermodynamics would get screwed up. Its
> > >intuitively the same as having different cut offs for LJ vs Coul.
> >That is an approximation, and since LJ potential goes to zero faster it
> >can be justified (although that does not mean it is correct).
> >Treating different parts of the system with different
> >cut-off is worse. Let's say we divided a 6x6x6 nm system in two pieces,
> >that with X coordinate < 3 nm and that with X > 3 nm. Now if we use a
> >different cut-off for the first half of the system than for the other
> >half, we create artifacts on the border between the two, and what's
> >more,
> >I suspect that for such a system in "equilibrium" the chemical potential
> >
> >is not the same in both halves. It will hence be impossible to define
> >what
> >kind of ensemble is modeled. Treating cut-offs differently for the
> >protein
> >and water is the same kind of division, comparable in nature to
> >simulations in implicit solvent.
> >
> > >
> > >>You would not want to do that anyway--it would really put you in
> > >>thermodynamic neverland.
> > >I agree, but IIRC this is what Amber does quite frequently, no cut-off
> > >within the protein, and short cut-off for protein-water and water-water
> >
> > >(e.g. the 1 us folding sim of Duan and Kollman).
> >
> >And to come back to that simulation, they also determined solvation free
> >
> >energies from that, but from just the protein coordinates, dismiising
> >all the
> >solvent.
> >
> >Just my (possibly biased) opinion...
> >
> >--
> >David.
> >________________________________________________________________________
> >David van der Spoel, PhD, Assist. Prof., Molecular Biophysics group,
> >Dept. of Cell and Molecular Biology, Uppsala University.
> >Husargatan 3, Box 596,          75124 Uppsala, Sweden
> >phone:  46 18 471 4205          fax: 46 18 511 755
> >spoel at xray.bmc.uu.se    spoel at gromacs.org   http://xray.bmc.uu.se/~spoel
> >++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> >
> >_______________________________________________
> >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.
> >
> >
> >_______________________________________________
> >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.
>
> ---------------------------------
> Eric Jakobsson, Ph.D.
> Professor, Department of Molecular and Integrative Physiology, and of
> Biochemistry
> Senior Research Scientist, National Center for Supercomputing Applications
> Professor, Beckman Institute for Advanced Science and Technology
> 4021 Beckman Institute, mc251
> 405 N. Mathews Avenue
> University of Illinois, Urbana, IL 61801
> ph. 217-244-2896       fax 217-244-2909
>
> (Currently on leave to NIH to serve as Director of Center for
> Bioinformatics and Computational Biology at the National Institute of
> General Medical Sciences and Chair of the NIH Biomedical Information
> Science and Technology Initiative Consortium.  Usual schedule is four days
> a week at NIH and three days a week at Illlinois.)
>
>
>
> _______________________________________________
> 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.
>
>
> _______________________________________________
> 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.
>