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

[Andrey Frolov]

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.