Hello everyone, my system is a composite, I need to calculate the temperature of a certain group of atoms. I searched a little and have found this post. I also read the doc pages and tried both methods mentioned in that post, i.e. gmx traj -ot and gmx mdrun -rerun with .tpr and .xtc files of that group of atoms.
But the problem I found is the same as that asked in that post, the temperature obtained is consistently lower. I set the temperature of that group to be 300 K for temperature coupling with tc-grps options. But I got about 250 K from gmx traj -ot and about 50 K from gmx mdrun -rerun.
The results of gmx mdrun -rerun for the whole system seem to be correct, but when is applied to a group of the atoms, it turns to be lower and the results seems to be wrong I think.
Could you please tell me what is the correct way to obtain the temperature of a group of the atoms? Any help would be much appreciated!
I have tried to calculate the number of constraints, but sorry I do not know how to do it correctly. I checked with the whole system and two parts of the system, I found if, for all of them, the ratio of [natoms / (natoms - nconstraints)] is 1.2, that would match the temperature I used for temperature coupling. But the information in the log file The number of constraints is XXX gives another value of constraint number. And they do not match with each other. I am sorry maybe this is a naive question, but could you please tell me how should I get the number of nconstraints?
Could you please also teach me why the gmx mdrun -rerun protocol can not get the correct temperature? I am asking because I used it for calculating Potential(A+B)-Potential(A)-Potential(B) to get the bonding energy. If this method can not get the temperature correctly, I am wondering is that a good idea to use it to get the other kind of properties of a group of the system.
For rigid water (SETTLE), each water molecule has 3 constraints. Then it depends on how you constrained other species. If all bonds were constrained, then add the number of bonds to the number of constraints. If only bonds to H were constrained, add the number of H atoms in the non-water part of the system, and you’ll have your result.
The validity of mdrun -rerun here is dependent upon whether you provided a trajectory with velocities or not. An .xtc file does not have velocity information, so the results you get in terms of temperature are flawed. The -rerun functionality should only be used for evaluating energies, anyway. Even a trajectory with velocities can’t give you an accurate temperature because it only contains instantaneous values of velocity for frames that were saved, ignoring all the rest.
Hello,
Thank you very much for the detailed answer!
I do not have water in my system. And I only applied constraint to hydrogen bonds. Following your answer, I found that the number of H atoms matches the number of constraints indicated in the log file. Thanks very much!
For the corrections of the results of gmx traj -ot, I found that if I use the degree of freedom for the correction, the value will be very close to the value obtained from gmx energy for the whole system, as well as the temperature which is set for every parts for the temperature coupling in the .mdp file. I mean I use [3*natoms / (3*natoms - nconstraints)] to calculate the ratio for the correction. I feel it also sounds reasonable to use the degree of freedom. Could you please tell me do you think I am doing this correctly?
I will remember that mdrun -rerun could be only used for evaluating the energies. Thank you very much for the help!
Yes, that sounds more accurate. I dug up the old formula in my last post from a previous forum post and it sounded correct offhand but the DOF approach is correct.
Thank you very much for the help!
I just have one more small question. Since Gromacs knows the number of constraints applied to the system, is that possible for the developers to make an improvement for gmx traj -ot to calculate the temperature of the specific group directly? If it can finish the job directly, that would be a little better.
Submit a feature request through the GitLab site. I imagine this could be incorporated if the user supplies a .tpr file, which is already an optional input. I imagine it’s been (and will be) a low priority since the correction is straightforward. Maybe even something as simple as adding the scaling factor to the help description would suffice.
Sorry for the late reply. Yeah, it is true that this is a low priority issue. And maybe a little more instructions merged into the manual would be sufficient. But it would be better if the manual can contain more content and be easier to understand.
And thank you for all the answers! That helps me a lot.
You should have seen the documentation before I got my hands on it many years ago.
In any case, if you have specific recommendations, please file a GitLab issue. Simply asking for “more content” and being “easier to understand” is an impossible request to satisfy. If there are specific issues that are not properly documented, that is useful feedback.