[gmx-developers] Pressure calculation in Gromacs

Berk Hess hessb at mpip-mainz.mpg.de
Mon Jan 5 12:07:15 CET 2009 

Hi,

There are mainy subtleties to this issue.

One important thing that you did not mention is the time at which you 
take the velocities.
The velocities in Gromacs 4.0 are for a leap-frog integrator, which 
means that
there are at half steps. To obtain Ekin at the full step the kinetic 
energy at minus
and plus a half step are averaged.

On explanation could be that your code is using Ekin determined with the
full step velocities.

Berk

Dmitri Rozmanov wrote:
> Hi,
>
> First of all, thank you very much for such a prompt response. I really 
> appreciate your help and explanation.
>
> Yes, of course I compare exactly the same model. I carefully checked 
> all the parameters of TIP4P I am using for comparison. I use the 
> parameters Gromacs uses.
>
> Then, you can check my pressure calculation yourself
>
> V = 75.7206 nm3
>
> P = 16.6 * 2 * Ekin_trn / V
>
> Then, for translational kinetic molecular energy =
>     2906.0216154     -11.2780088       3.588625
>      -11.2780088    2989.8108419      33.6057171
>        3.588625      33.6057171    2813.6215204
> =  8709.45 kJ/mol
>
> this gives the Kin Stress =
>     -1274.5712879286        4.94649665553416       -1.57395883220682
>         4.94649665553416    -1311.32096032219      -14.7393542769693
>        -1.57395883220682      -14.7393542769693    -1234.04491762744
> --> P = -trace(Stress) / 3. = 1273.31238862608 bar
>
> This is for the case when all the interaction are suppressed.
>
> This, is what our code gives and this is what I expect it should be.
> I cannot get this number in Gromacs.
>
> And yes, our code uses molecular virials, is a sense that we calculate 
> site virial and subtract the rigid molecular corrections.
>
> How it is done in the main SPME paper:Essmann et al., 1995 / page 8579
>
> It is very similar to what Gromacs does, it uses site based virial and 
> then subtracts the correction due to the constraints.
>
> At the moment I fail to see how Gromacs pressure is correct, and why I 
> should not compare site based pressure with the molecular based 
> pressure. Or, to put it in a differently, which pressure is more like 
> the REAL one?
>
> Gromacs gets it like this:
>
> Full kinetic energy (kJ/mol) =
>     5727.552282     -41.4771432     -70.7586576
>      -41.4771432    5916.4805532      28.4006926
>      -70.7586576      28.4006926    5754.1196287
>
> Full kinetic stress (bar) =
>     -2512.08513042746       18.1917352396384       31.0345087838987
>        18.1917352396384    -2594.94843353347      -12.4564480709355
>        31.0345087838987      -12.4564480709355    -2523.73747916459
>
> Constraint virial (kJ/mol) =
>     1457.7536918      18.1317179      -1.9607236
>       17.9054384    1425.9788285      -3.6860869
>       -1.9082659      -3.1906363    1462.1762396
>
> gives constraint stress (bar) =
>     639.365856948194      7.95251037146193     -0.85996676379849
>       7.85326494526879    625.429508978373     -1.61670529312393
>      -0.83695899334818     -1.39940232951192    641.305571510364
>
> With the final stress =
>     -1872.68       26.0407       30.1104
>        26.1022    -1969.55      -14.1334
>        30.3522      -14.4231    -1883.07
>
> so that the Gromacs pressure (bar) = 1908.4
>
> And the difference is, as I said, more than 600 bar.
>
> Where am I wrong?
>
> Thank you again. Best wishes,
>
>     --- Dmitri Rozmanov.
>
> David van der Spoel wrote:
>> Dmitri Rozmanov wrote:
>>> Hello All,
>>>
>>> I am comparing the results of Gromacs with our in-house MD 
>>> simulator. It seems that I have general agreement with energies, and 
>>> forces, but the pressure is somehow different.
>>>
>>> I am doing one point calculation of a system of 2400 TIP4P water 
>>> molecules in a rectangular elongated box.
>>>
>>> I was actually able to separate most of the pressure contributions 
>>> Gromacs calculates and
>>> have some troubles with understanding how the shake/links/settles 
>>> virials are obtained and used.
>>>
>>> All force virials are the same, the trouble is with the constraints 
>>> virial.
>>>
>>> Our code is for rigid dynamics and it explicitly uses translations 
>>> and rotations of the molecules, waters in this case. To debug the 
>>> pressure I removed the charges from TIP4P sites and zeroed LJ params 
>>> as well. So, there are no nonbonded interactions in the system at 
>>> all, and all the molecules are rigid bodies.
>>>
>>> I can assume, that due to the equi-partitioning the kinetic energy 
>>> is split between the rotations and translations equally and only the 
>>> translational part of it will contribute to the pressure. The 
>>> rotational kinetic virial should be canceled in the pressure due to 
>>> constraints in the system. Right?
>>>
>>> So, we should have someting like this: P = 16.6 * 2 * Ekin_trn / 3 / V
>>>
>>> In my topology I have
>>>     EKinetic Trn (kJ/mol) = 8709.4539759232
>>>     EKinetic Rot (kJ/mol) = 8688.6908125891
>>>     EKinetic Tot (kJ/mol) = 17398.1447885123
>>>
>>> Wich gives me P = 1273.3 bar
>>>
>>> This is the number I get in Gromacs if I remove the rotations from 
>>> the starting topology:
>>>     KE = 8709.45 kJ/mol --> P = 1273.3 bar.
>>>
>>> However, I both the translations and rotations are retained I get this:
>>>
>>> KE = 17398.1 kJ/mol --> P = 1883.9 bar = 2543.6 - 659.7
>>>
>>> Where, the first is the pressure due to the total KE, and the second 
>>> is the pressure due to the constraints. The kinetic part is 
>>> correctly about twice the translational only contribution, but
>>> the constraint virial does not cancel the rotational virial!
>>>
>>> The difference is about 600 bar... Well...
>>>
>>> If all the interactions are there, it becomes even worse:
>>>
>>> My calculation P = -112.3 bar, Gromacs gives -2238.98 bar
>>>
>>> The pressure is negative due to the sightly different water model 
>>> used to obtain the topology. This is normal.
>>>
>>> Now I have 2100 bar difference.
>>>
>>> I changed the tip4p.itp in a way that Gromacs uses shake/lincs 
>>> instead of settles, and the pressures are almost the same for all 
>>> three of them.
>>>
>>> So, am I doing something wrong? Which pressure is more correct? If 
>>> Gromacs pressure is not correct, is it a big deal?
>>>
>>
>> gromacs pressure is correct. but you should compare identical models.
>> check the description in the appendix, and note that gromacs uses an 
>> atomic virial rather than a molecular virial which your code might do.
>> note that you can also inspect the virial from the energy file.
>>
>>> Thank you very much. Any suggestions will be appreciated.
>>>
>>> With best wishes and Happy New Year!
>>>
>>>     --- Dmitri.
>>>
>>>
>>>
>>>
>>>
>>>
>>> _______________________________________________
>>> gmx-developers mailing list
>>> gmx-developers at gromacs.org
>>> http://www.gromacs.org/mailman/listinfo/gmx-developers
>>> Please don't post (un)subscribe requests to the list. Use the www 
>>> interface or send it to gmx-developers-request at gromacs.org.
>>
>>
>
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the www 
> interface or send it to gmx-developers-request at gromacs.org.

More information about the gromacs.org_gmx-developers mailing list