[gmx-users] Parrinello-Rahman barostat: temperature and potential energy depend on tau_p
Andrey Frolov
andrey.i.frolov at mail.ru
Thu Dec 26 13:18:16 CET 2013
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.
More information about the gromacs.org_gmx-users
mailing list