[gmx-developers] Pressure calculation in Gromacs
Dmitri Rozmanov
rozmanov at gmail.com
Fri Jan 2 17:43:30 CET 2009
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.
>
>
More information about the gromacs.org_gmx-developers
mailing list