[gmx-users] Parrinello-Rahman barostat: temperature and potential energy depend on tau_p

[Justin Lemkul]

On 1/7/14, 6:37 AM, Andrey Frolov wrote:
> Dear GMX-Users,
>
> Still no reply. I repost the questions. Any suggestions?
>
>

It would be useful to know whether the problem is in the barostat itself or in 
the combination of integrator and barostat.  If you have time to test various 
combinations (at least the standard 'md' integrator), that would be useful and 
would speed up finding the source of the problem.  In either case, posting a 
test case and this description of the problem to redmine.gromacs.org is the best 
way to achieve a resolution in a future version.

-Justin

> Andrey Frolov wrote
>>
>> I am simulating pure supercritical CO2 fluid in NPT ensemble. Reference
>> temperature is 308.15K. I use “sd” (Langevin dynamics) to control
>> temperature with tau_t = 2ps. After 100 ps of equilibration  in NPT with
>> Berendsen thermostat, I run production simulation of 10 ns with
>> Parrinello-Rahman barostat. I checked several values of tau_p (from 2 to
>> 100 ps) with the compressibility set to 4.5e-5 bar^-1.
>>
>> 1) Temperature as given in the mdrun log file and edr files is higher than
>> the reference temperature. Moreover, the temperature depends on tau_p: the
>> larger tau_p the lower the resulting T:
>>   tau_p                               Quantity                      Average
>> Err.Est.    RMSD  Tot-Drift
>>         2.0                            Temperature                 437.688
>> 6.8    82.8611    5.83615  (K)
>>         5.0                            Temperature                 321.009
>> 0.49    10.1346    1.51488  (K)
>>        10.0                            Temperature                 314.071
>> 0.18    7.06755   0.588431  (K)
>>        15.0                            Temperature                 312.508
>> 0.2    6.48075   0.725419  (K)
>>        20.0                            Temperature                 311.831
>> 0.12    6.48166 -0.0137431  (K)
>>        25.0                            Temperature                 311.783
>> 0.093    6.48474  -0.713862  (K)
>>        30.0                            Temperature                 311.289
>> 0.099    6.41436   0.764874  (K)
>>        50.0                            Temperature                 311.392
>> 0.15     6.5612   0.585978  (K)
>>       100.0                            Temperature                 311.278
>> 0.11    6.51524   0.488872  (K)
>>
>> What could cause such behavior? Could it be that the kinetic energy
>> associated with additional degrees of freedom (motion of the the
>> simulation box edge length) added by Parrinello-Rahman barostat is NOT
>> subtracted from the kinetic energy of the system when the temperature is
>> calculated?
>>
>> 2) Potential energy also depends on tau_p:
>>   tau_p                               Quantity                  Average
>> Err.Est.    RMSD  Tot-Drift
>>         2.0                       Potential                   374.131
>> 120     1535.4     89.344  (kJ/mol)
>>         5.0                       Potential                  -1757.13
>> 11    225.176    24.8182  (kJ/mol)
>>        10.0                       Potential                  -1879.24
>> 4.6     186.86    5.23366  (kJ/mol)
>>        15.0                       Potential                  -1912.12
>> 5.3    183.542    25.0223  (kJ/mol)
>>        20.0                       Potential                  -1927.11
>> 4.1    198.066   -9.76647  (kJ/mol)
>>        25.0                       Potential                  -1922.43
>> 7    187.913   -22.8972  (kJ/mol)
>>        30.0                       Potential                  -1939.59
>> 5.9    175.813    32.3813  (kJ/mol)
>>        50.0                       Potential                  -1935.01
>> 8.4      182.9    11.7818  (kJ/mol)
>>       100.0                       Potential                  -1930.76
>> 8.3      199.6    40.4494  (kJ/mol)
>>
>> I would appreciate if someone can explain the dependence of potential
>> energy on “tau_p”.
>> Could it be explained by the implementation problem of Parinello-Rahman
>> with leap-frog algorithm discussed in the manual:  “... using the
>> leap-frog algorithm, the pressure at time t is not available until after
>> the time step has completed, and so the pressure from the previous step
>> must be used, which makes the [Parrinello-Rahman] algorithm not directly
>> reversible, and may not be appropriate for high precision thermodynamic
>> calculations.”?
>>
>> Densities and pressures are identical among these simulations. As I can
>> see from the density VS time plots the period of density oscillations
>> corresponds roughly to 4*tau_p.
>>
>> Also, it seems that the potential energy goes to plateau when tau_p > 20
>> ps. Does it show that for my system I should not use tau_p less than 20 ps
>> at given compressibility value and that with tau_p > 20 ps the true NPT
>> ensemble is sampled?
>>
>> Thank you very much.
>> Andrey
>
>
>
> -----
> ----
> Andrey I. Frolov, PhD
> Junior researcher
> Institute of Solution Chemistry
> Russian Academy of Sciences
> --
> View this message in context: http://gromacs.5086.x6.nabble.com/Parrinello-Rahman-barostat-temperature-and-potential-energy-depend-on-tau-p-tp5013532p5013626.html
> Sent from the GROMACS Users Forum mailing list archive at Nabble.com.
>

-- 
==================================================

Justin A. Lemkul, Ph.D.
Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441

==================================================