I would like to calculate the kinetic energy and/or temperature locally in a lipid bilayer. I have already run these simulations and I have .trr files with velocities. There were three temperature groups in the simulation, each coupled with a separate thermostat: membrane, water+ions, protein. I would like to be able to get molecule-level detail for the lipids. If I could read off the velocity values from the .trr, somehow, I could at least calculate kinetic energies for the molecules. What are my options here?
I figured this out on my own and I’d like to write up what I did in case it’s helpful to others. Here’s what I did:
1- Write a script to edit the index file to create a group for each molecule of interest. This could be done by hand with gmx make_ndx but it’s way easier to write a script to do it, if you have more than a handful of molecules you want to calculate temperatures of.
2- Run gmx traj with the -ot flag, and with the index file you created fed as -n index_custom.ndx. This allows you to choose a group. At this point you would choose a molecule of interest. Again it’s hard to do this by hand, so writing a script (in this case a csh script) with a loop that has something like the following helps:
3- Scale the temperatures by 3*N/N_constraints to get the real temperatures. The number of constraints for temperature groups are printed out when gmx grompp is run. In my case, one temperature group was the membrane. It has 200 lipids. In my case, each lipid molecule has the same number of degrees of freedom, and the same number of constraints, as they are identical. Therefore, I can just use the ratio for the entire membrane to scale temperatures of the individual molecules.