[gmx-developers] Pressure calculation in Gromacs

[Dmitri Rozmanov]

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.
> 
>