<html>
<head>
<style>
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 10pt;
font-family:Verdana
}
</style>
</head>
<body class='hmmessage'>
I have been looking at very similar things recently.<br><br>Your result could actually be correct.<br>If the water orders at the interface you should have a difference<br>in Ekin components, since the constraints, which remove kinetic energy<br>are not ordered isotropically in space.<br>You should be able to quantify this by looking at the order of the water<br>normal.<br><br>Berk<br><br>> Date: Thu, 10 Sep 2009 18:11:16 +0200<br>> From: alexander.herz@mytum.de<br>> To: gmx-users@gromacs.org; felix.sedlmeier@ph.tum.de<br>> Subject: Re: [gmx-users] non isotropic kinetic energy<br>> CC: <br>> <br>> Hi,<br>> <br>> we took a couple of days to reexamine this and try out a couple of things:<br>> <br>> -On small interfacial systems (ca. 2k spc/e water molecules sandwiched<br>> by vacuum)<br>> the difference between ekin_x/y and ekin_z is 20 to 60kJ/mol depending<br>> on the specific system/settings<br>> (e.g. see attached file of averaged ekin (x/y/z) over time).<br>> -On larger interfacial systems the difference grows (up to several<br>> 100kJ/mol for ca. 32k spc/e water).<br>> <br>> The following settings make no difference (so the error persists even if<br>> we change the setting):<br>> -thermostat (we tryed nose-hoover, berendsen and no thermostat)<br>> -estatics (cut-off instead of pme)<br>> -type of constraint algo for water (so lincs, shake or settle, we<br>> modified the spc/e water molec definition to change this)<br>> -changing system size so that the box length is equal in all 3 dim (as<br>> said, magnifies problem)<br>> -rotating system (error rotates with system)<br>> -use smaller time step of 1fs (no difference)<br>> <br>> We ran a simulation with flexible water (so bonds instead of<br>> constraints) , here the error vanishes!<br>> Also, we ran similar sims with amber, where the error does not appear at<br>> all (for constraint water).<br>> <br>> We found out about this problem while doing surface tension/contact<br>> angle calculations.<br>> The surf tens calced via virial differs from the surf tens calculated<br>> via pressure tensor (unless one corrects<br>> for the anisotropic kinetic energy). A simple way to see the deficit is<br>> trying to calc the average pressure tensor from the average virial<br>> tensor (using the simple def from the gromacs manual P_ij=2/V(Ekin-Vir_ij):<br>> a) assuming equidistributed kinetic energy (ekin/3 per dim)<br>> b) using the anisotropic ekin which can be extracted from a trr with a<br>> tool (source attached) if velocities have been written to the trr.<br>> <br>> All this appears to hint at a problem with the constraints for the<br>> waters and might also hint at an incorrect ensemble as a consequence.<br>> <br>> Alex<br>> <br>> sample input file (several variants of this file have been used):<br>> <br>> title = pure water<br>> cpp = /usr/bin/cpp<br>> integrator = md<br>> nsteps = 2500000<br>> nstlist = 50<br>> nstxout = 1000<br>> nstvout = 1000<br>> nstlog = 1000<br>> dt = 0.002<br>> nstcomm = 200<br>> comm-mode = linear<br>> constraints = h-bonds<br>> nstenergy = 100<br>> ns_type = grid<br>> coulombtype = pme<br>> fourierspacing = 0.12<br>> pme_order = 4<br>> rlist = 0.9<br>> rvdw = 0.9<br>> rcoulomb = 0.9<br>> tcoupl = Berendsen<br>> tc_grps = SOL<br>> energygrps = SOL<br>> tau_t = 0.4<br>> ref_t = 300<br>> <br>> <br>> <br>> David van der Spoel schrieb:<br>> > aherz wrote:<br>> >> Hi,<br>> >><br>> >> we're running simulations of a water slab sandwiched by vacuum (ca. 2k<br>> >> waters, 3x3x12 nm^3 box) in NVT<br>> >> and look at the pressure tensor. We are surprised to find that the<br>> >> kinetic energy is equivalent (within statistical uncertainty) for x and<br>> >> y dim(parallel to the interface) but not in z dim. This appears to<br>> >> violate the equipartition theorem?<br>> >><br>> >> This is gromacs 3.0 and 4.0 using berendsen thermostat.<br>> ><br>> > First, try it with a real thermostat, second please report some numbers.<br>> ><br>> >><br>> >> Alex<br>> >> _______________________________________________<br>> >> gmx-users mailing list gmx-users@gromacs.org<br>> >> http://lists.gromacs.org/mailman/listinfo/gmx-users<br>> >> Please search the archive at http://www.gromacs.org/search before<br>> >> posting!<br>> >> Please don't post (un)subscribe requests to the list. Use the www<br>> >> interface or send it to gmx-users-request@gromacs.org.<br>> >> Can't post? Read http://www.gromacs.org/mailing_lists/users.php<br>> ><br>> ><br>> <br><br /><hr />See all the ways you can stay connected <a href='http://www.microsoft.com/windows/windowslive/default.aspx' target='_new'>to friends and family</a></body>
</html>