[gmx-users] non isotropic kinetic energy

aherz alexander.herz at mytum.de
Mon Sep 14 16:08:55 CEST 2009


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




More information about the gromacs.org_gmx-users mailing list