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

[Dr. Vitaly Chaban]

I would suggest to repost this in the bug management system.

http://bugzilla.gromacs.org/


Dr. Vitaly V. Chaban


On Tue, Jan 7, 2014 at 12:37 PM, Andrey Frolov <andrey.i.frolov at mail.ru> wrote:
> Dear GMX-Users,
>
> Still no reply. I repost the questions. Any suggestions?
>
>
> 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.
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.