[gmx-developers] Pressure calculation in Gromacs
Dmitri Rozmanov
dmitri.rozmanov at ucalgary.ca
Tue Jan 6 00:51:24 CET 2009
Thank you Berk for the tips...
I will further look into it and will let you know if I find anything new.
Thank you again.
--- Dmitri.
Berk Hess wrote:
> 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.
>
>
--
-----------------------------------------------
Dmitri Rozmanov
PhD Student
Department of Chemistry
University of Calgary
2500 University Dr. NW,
Calgary, Alberta T2N 1N4, Canada
tel.: (403) 220-7561
e-mail: dmitri.rozmanov at ucalgary.ca
More information about the gromacs.org_gmx-developers
mailing list