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

[Andrey Frolov]

Dear GMX-Users,

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-tp5013532.html
Sent from the GROMACS Users Forum mailing list archive at Nabble.com.