Calculating the system temperature from the velocity information from .gro file

User discussions
, ,
1

GROMACS version:2023.4
GROMACS modification: Yes/No
Here post your question

Hello. I am a beginner of GROMACS.
I have a question about how to calculate the temperature from the velocity information in the .gro file. This is the format of .gro file that I currently have which looks like

Title of the system t= 0.00000 step= 0
7581
1ALA BB 1 3.866 1.935 3.832 -0.1483 0.1788 0.2353
1ALA SC1 2 3.821 1.769 4.041 -0.5410 -0.3384 -0.2435
2ASP BB 3 3.974 2.251 3.932 -0.0148 -0.0533 -0.3008
2ASP SC1 4 3.939 2.276 3.559 0.2182 -0.2081 0.1115
3ASN BB 5 4.248 2.410 3.856 -0.2157 -0.1441 0.0451
3ASN SC1 6 4.379 2.277 3.565 -0.1033 0.0424 0.1119
4TYR BB 7 4.072 2.726 3.864 -0.0584 -0.2129 0.0969
……….

The last three columns indicate the velocity (x,y,z). By using these velocity values, I want to calculate the temperature of the system. I am trying to use the three equations described in the below link [Molecular Dynamics — GROMACS 2021.1 documentation 23].
If you look at this link, I am trying to use the three equations (7), (8) and (9).
I think the velocity given in the .gro file, the unit should be nm/ps so when I calculate the velocity, I converted the Boltzmann’s constant to match the unit which I think is (1.380649 * 10*-29 nm^2 * kg * ps^(-2) * k^(-1)). The mass depends on the type of the protein, water and ions but I’m pretty sure those should be correct. 1.67377e-2718 kg for water, 20.991.67377e-27 kg for ions… Etc.
For the Nc, I assigned 0 and 3 for Ncom.
By using those given values, I successfully calculated the temperature of the whole system.
(I calculated the temperature for each of the time steps I assigned ). The thing is, it is very weird because the values I got from the calculations are ⅓ of the correct answers. I ran MD simulation by using gmx energy to see if what I did was correct but the values from the MD simulations are 3 times bigger than values from my calculations consistently. So I assume I am basically on the right track but maybe I am missing/miscalculating something?

For example these are the values I got from my calculations (first column is timestep and second colum is the temp in K)

0.0 93.37701473340636
20.0 93.89294647489554
40.0 96.20400246375763
60.0 97.36930355169025
80.0 95.73484956466677
100.0 90.39040962452691
120.0 105.2233903793852
……….

And these are the values I got from MD simulation

0.000000 302.080750
2.000000 300.867462
4.000000 298.133179
6.000000 301.865265
8.000000 301.926117
10.000000 302.514252
12.000000 296.901550
14.000000 301.065033
……….

I would be grateful if you could give me some advice on why I got different values but differences are very consistent. If you need more info, please feel free to comment! Thanks

2

I don’t really understand what you mean with the masses you specified. I’d suggest calculating the kinetic energy (in J) using the mass and velocity of each atom. Then you just use Boltzmann’s constant 1.380649e-23 J/K. That should reduce the risk of mistakes.

3

J K−1 is equal to kg m2 s−2 K−1. The reason why I did the unit conversion is because the velocity I am trying to use is not in m/s.
So in my .gro file, there are 6 types of atoms
3 types of protein, 2 types of ions and 1 is water
the mass of the three different proteins are (36, 54, 72 amu)
2 ions are 40, 22.99amu
and water for 18amu
I converted every mass into kg so that I can use boltzmann’s constant

4

I thought the easiest unit conversion would be the velocity to m/s, 1 nm/ps = 1000 m/s.

But the main source for the errors would be the number of degrees of freedom, I guess. How many constraints do you have? Do you have more than one temperature coupling group?

5

Ndf = 3N−Nc−Ncom.

for the N, it should be the number of atoms so N=7581 for my case and I did
Nc=0 and Ncom=3.
The fun thing is, if I do not multiply 3 with N, I get the right answer. But I cant do that cause Ndf = 3N−Nc−Ncom. is the correct formula