[gmx-users] non isotropic kinetic energy

[aherz]

Hey,

I thought a little more about this. Eventually, the water molecules are
rigid bodies.
The rigidity constraint is implemented via constraints (settle or so)
but could also (in theory) be implemented
using some euler angles or what ever scheme you prefer. It should not
matter how you implement
the rigidity, the result of a rigid body in some potential should be
identical (if implemented correctly).

The equipartition theorem must hold for a rigid body (it has 3
translational and 3 (for water) rotational degrees
of freedom which must contribute 1/2*k_b*T per molecule, irrespective
how you implement rigidity.
Also I wouldn't see how you could explain the energy difference if the
molecules were actually implemented
as rigid bodies with e.g. euler angles. Why should the euler angles care
about the water orientation at the interface boundary?

I ran some more analysis assuming the water molecules are rigid bodies.
This shows that the translational energy of the molecules is actually
correct (so 1/2k T) but the average rotational energy is incorrect (with
settle the z part is too high and x/y ok,with lincs x/y is too low and z
too high).

Summarizing, I fail to see why choosing one method over the other for
performing rigid body dynamics should have an influence over
equipartition of energy in the system. I'm still not convinced the
ensemble is correct.

Alex


Berk Hess schrieb:
> Equipartitioning still holds.
> The issue here is that constraints remove degrees of freedom.
> If you remove a degree of freedom, there is not potential or kinetic
> energy
> anymore in that generalized coordinate.
> This does not introduce any problems with ensembles or something like
> that.
>
> I did not have time yet to implement pull code with pbc.
> I will do it this week.
>
> Berk
>
> > Date: Mon, 14 Sep 2009 16:08:55 +0200
> > From: alexander.herz at mytum.de
> > To: gmx-users at gromacs.org
> > Subject: Re: [gmx-users] non isotropic kinetic energy
> >
> > Hm,
> >
> > apparently there are constraint algorithms which don't change the
> ensemble
> > (e.g.
> http://www.informaworld.com/smpp/content~db=all~content=a750934359).
> > The equipartition theorem requires that the energy depends on
> > no other powers of v but v^2
> > (http://www.whfreeman.com/modphysics/PDF/8-2c.pdf)
> > which should be the case for the md systems including constraints
> > (lj pot and estatic only depend on the particles' positions). So as far
> > as I can see the
> > equipartition theorem should hold in td equilibrum for these systems. So
> > the constraints
> > might modify the ensemble which could lead to all kinds of trouble?
> >
> > The water molecules at the interface are oriented but taking the above
> > paragraph into
> > account this shouldn't cause any trouble.
> >
> > By the way,
> > did you find a chance to look at the pulling problem we discussed before
> > last week?
> >
> > Thx,
> > Alex
> >
> >
> > Berk Hess schrieb:
> > > I have been looking at very similar things recently.
> > >
> > > Your result could actually be correct.
> > > If the water orders at the interface you should have a difference
> > > in Ekin components, since the constraints, which remove kinetic energy
> > > are not ordered isotropically in space.
> > > You should be able to quantify this by looking at the order of the
> water
> > > normal.
> > >
> > > Berk
> > >
> > > > Date: Thu, 10 Sep 2009 18:11:16 +0200
> > > > From: alexander.herz at mytum.de
> > > > To: gmx-users at gromacs.org; felix.sedlmeier at ph.tum.de
> > > > Subject: Re: [gmx-users] non isotropic kinetic energy
> > > > CC:
> > > >
> > > > Hi,
> > > >
> > > > we took a couple of days to reexamine this and try out a couple of
> > > things:
> > > >
> > > > -On small interfacial systems (ca. 2k spc/e water molecules
> sandwiched
> > > > by vacuum)
> > > > the difference between ekin_x/y and ekin_z is 20 to 60kJ/mol
> depending
> > > > on the specific system/settings
> > > > (e.g. see attached file of averaged ekin (x/y/z) over time).
> > > > -On larger interfacial systems the difference grows (up to several
> > > > 100kJ/mol for ca. 32k spc/e water).
> > > >
> > > > The following settings make no difference (so the error persists
> even if
> > > > we change the setting):
> > > > -thermostat (we tryed nose-hoover, berendsen and no thermostat)
> > > > -estatics (cut-off instead of pme)
> > > > -type of constraint algo for water (so lincs, shake or settle, we
> > > > modified the spc/e water molec definition to change this)
> > > > -changing system size so that the box length is equal in all 3
> dim (as
> > > > said, magnifies problem)
> > > > -rotating system (error rotates with system)
> > > > -use smaller time step of 1fs (no difference)
> > > >
> > > > We ran a simulation with flexible water (so bonds instead of
> > > > constraints) , here the error vanishes!
> > > > Also, we ran similar sims with amber, where the error does not
> appear at
> > > > all (for constraint water).
> > > >
> > > > We found out about this problem while doing surface tension/contact
> > > > angle calculations.
> > > > The surf tens calced via virial differs from the surf tens
> calculated
> > > > via pressure tensor (unless one corrects
> > > > for the anisotropic kinetic energy). A simple way to see the
> deficit is
> > > > trying to calc the average pressure tensor from the average virial
> > > > tensor (using the simple def from the gromacs manual
> > > P_ij=2/V(Ekin-Vir_ij):
> > > > a) assuming equidistributed kinetic energy (ekin/3 per dim)
> > > > b) using the anisotropic ekin which can be extracted from a trr
> with a
> > > > tool (source attached) if velocities have been written to the trr.
> > > >
> > > > All this appears to hint at a problem with the constraints for the
> > > > waters and might also hint at an incorrect ensemble as a
> consequence.
> > > >
> > > > Alex
> > > >
> > > > sample input file (several variants of this file have been used):
> > > >
> > > > title = pure water
> > > > cpp = /usr/bin/cpp
> > > > integrator = md
> > > > nsteps = 2500000
> > > > nstlist = 50
> > > > nstxout = 1000
> > > > nstvout = 1000
> > > > nstlog = 1000
> > > > dt = 0.002
> > > > nstcomm = 200
> > > > comm-mode = linear
> > > > constraints = h-bonds
> > > > nstenergy = 100
> > > > ns_type = grid
> > > > coulombtype = pme
> > > > fourierspacing = 0.12
> > > > pme_order = 4
> > > > rlist = 0.9
> > > > rvdw = 0.9
> > > > rcoulomb = 0.9
> > > > tcoupl = Berendsen
> > > > tc_grps = SOL
> > > > energygrps = SOL
> > > > tau_t = 0.4
> > > > ref_t = 300
> > > >
> > > >
> > > >
> > > > David van der Spoel schrieb:
> > > > > aherz wrote:
> > > > >> Hi,
> > > > >>
> > > > >> we're running simulations of a water slab sandwiched by vacuum
> > > (ca. 2k
> > > > >> waters, 3x3x12 nm^3 box) in NVT
> > > > >> and look at the pressure tensor. We are surprised to find
> that the
> > > > >> kinetic energy is equivalent (within statistical uncertainty) for
> > > x and
> > > > >> y dim(parallel to the interface) but not in z dim. This
> appears to
> > > > >> violate the equipartition theorem?
> > > > >>
> > > > >> This is gromacs 3.0 and 4.0 using berendsen thermostat.
> > > > >
> > > > > First, try it with a real thermostat, second please report some
> > > numbers.
> > > > >
> > > > >>
> > > > >> Alex
> > > > >> _______________________________________________
> > > > >> gmx-users mailing list gmx-users at gromacs.org
> > > > >> http://lists.gromacs.org/mailman/listinfo/gmx-users
> > > > >> Please search the archive at http://www.gromacs.org/search before
> > > > >> posting!
> > > > >> 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
> > > > >
> > > > >
> > > >
> > >
> > >
> ------------------------------------------------------------------------
> > > See all the ways you can stay connected to friends and family
> > > <http://www.microsoft.com/windows/windowslive/default.aspx>
> > >
> ------------------------------------------------------------------------
> > >
> > > _______________________________________________
> > > gmx-users mailing list gmx-users at gromacs.org
> > > http://lists.gromacs.org/mailman/listinfo/gmx-users
> > > Please search the archive at http://www.gromacs.org/search before
> posting!
> > > 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://lists.gromacs.org/mailman/listinfo/gmx-users
> > Please search the archive at http://www.gromacs.org/search before
> posting!
> > 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
>
> ------------------------------------------------------------------------
> Express yourself instantly with MSN Messenger! MSN Messenger
> <http://clk.atdmt.com/AVE/go/onm00200471ave/direct/01/>
> ------------------------------------------------------------------------
>
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://lists.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> 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