GROMACS version: 5.1.5
GROMACS modification: No
Respected community members, I am a beginner in Gromacs. I have performed a 50 ns protein in water simulation. Now I would like to calculate the Total energy of my protein with respect to the time. I have followed what was mentioned in the Analysis part of the Protein-Ligand Tutorials.
First, I downloaded the mdp file that was linked to that page. I editted two main things in that file. Firstly, I changed energygrps = Protein, and secondly I changed the tc-grps = Protein Non-Protein (Since in that mdp file it was Protein_JZ4 Water_and_ions). I saved my editted mdp file as e.mdp. Then I used the grompp to convert e.mdp into e.tpr file using the following command:
gmx grompp -f e.mdp -c npt.gro -t npt.cpt -p topol.top -o e.tpr
Next, I used mdrun with the -rerun option to recalculate the energy from the existing simulation trajectory by executing the following command:
gmx mdrun -deffnm e -rerun md.xtc -nb cpu
After executing the command I got the following warning message:
WARNING: Some frames do not contain velocities.
Ekin, temperature and pressure are incorrect,
the virial will be incorrect when constraints are present.
After this an e.edr file was successfully generated and then I used the following command
gmx energy -f ie.edr -o energy.xvg
After this I selected the Total Energy (13) option from the interactive menu and got the Total Energy Graph. However, I feel due to the above Warning message my result is not accurate (since the velocities are missing). What can I do now to avoid the above warning and get the accurate result. Please guide me. Thanks in advance.
This does not cause only the protein’s energy to be written. energygrps are only for decomposition of short-range nonbonded energies.
In this case, the total energy reported is the same as the potential energy of the system, since there is no kinetic energy.
There are a couple of points to make here. The “total energy of the protein” is a physically meaningless quantity. If you want to calculate it, the only way to do so is to strip only the protein coordinates out of the trajectory and make a matching .tpr file with convert-tpr and use mdrun -rerun. But again I emphasize, this quantity means nothing. Force fields are not parametrized to reproduce absolute potential energies of configurations, and kinetic energy is a function of the temperature of the system.
Thanks Dr. Lemkul for replying to my question. Sir can you please share the command by which I can strip the protein coordinates alone and make a new .tpr file. Many thanks for you help sir.
Dr. Lemkul as you mentioned, I have converted my md.tpr file by using the following command:
gmx convert-tpr -s md_0_50.tpr -o protein.tpr
In the interactive menu I selected option 1 for Protein.
Then i used the command: gmx mdrun -deffnm protein -rerun md_0_50.xtc -nb cpu
Upon running this command I encountered the following error:
Fatal error: Number of atoms in trajectory (49891) does not match the run input file (3363).
Now I understand that the number of atoms in my md_0_50.xtc file do not match the number of molecules in the protein.tpr file since the latter file contains only protein. Thus, I tried to extract only the protein trajectory from md_0_50.xtc file by using the command:
gmx trjconv -f md_0_50.xtc -s md_0_50.tpr -o protein.xtc
Then, from the interactive menu I selected option 1 for protein and got the protein.xtc file. One question that I would like to ask here is my selection for -s in the above command is md_0_50.tpr. Is this correct? At this point I would also like to mention that my original md_0_50.xtc was generated but using md.mdp that did not specify any energygrps. Is this OK?
Then again I executed the command gmx mdrun -deffnm protein -rerun protein.xtc -nb cpu
This time there was successful execution but I got the same WARNING again stating
WARNING: Some frames do not contain velocities.
Ekin, temperature and pressure are incorrect,
the virial will be incorrect when constraints are present.
Dr. Lemkul I suspect after doing all this the velocities are still missing. Now how can I proceed from here onwards? Kindly give your valuable suggestion. Thanks in advance.
A new update… I figured out that my working directory does not have an md_0_50.trr file that’s why the velocities are missing. Actually this file was never generated in my original simulation. So what I have done is I have run my simulation again from scratch with an updated md.mdp file with the following changes so that the md_0_50.trr file gets generated.
nstxout = 5000 ; save coordinates every 10.0 ps
nstvout = 5000 ; save velocities every 10.0 ps
Now, after the simulation got over I have used the following command to generate a new .xtc file from my .trr file. I have used the following command to do this:
gmx trjconv -s md_0_50.tpr -f md_0_50.trr -o protein.xtc
In the interactive menu I selected option 1 for Protein. Now my xtc file contains only protein coordinates and velocities.
Since I will use -rerun so I generated a new .tpr file that contains the protein coordinates alone. I generated this .trp file by using the command gmx convert-tpr -s md_0_50.tpr -o protein.tpr. From the interactive menu I selected the option 1 for protein.
So now I have both the protein.xtc file (containing coordinates and velocity data of proteins) and protein.tpr file.
Next I executed -rerun by applying the command: gmx mdrun -deffnm protein protein.xtc -nb cpu
The resulting command gave me the output file protein.edr
Then by using the command gmx energy -f protein.edr -o Total_energy.xvg I got the value and graph.
Since I am a beginner in gromacs I would be extremely grateful to any expert in Gromacs who can go through my procedure and confirm that whatever I have done is correct. Please help me. Thanks in advance.
The procedure is correct, but again I caution that the value you obtain has no physical meaning.
Many thanks Dr. lemkul for your confirmation. I am aware that this quantity may not be that much significant. I was just interested in finding this out as I was comparing my results with the previously published data. Once again thanks for your kind support.
There is another update in this post. Recently i have run another protein in water simulation for 100 ns. My md.mdp file had
nstxout = 5000 ; save coordinates every 10.0 ps
nstvout = 5000 ; save velocities every 10.0 ps
So that the .trr file gets generated.
And then inorder to calculate the total energy i made a new tpr file containing only protein coordinates (convert-tpr -s md_0_100.tpr -o protein.tpr. Then i also made another xtc file from the trr file (gmx trjconv -s md_0_100.tpr -f md_0_100.trr -o protein.xtc. (selected protein when prompted).
Then when i gave the command
gmx mdrun -deffnm protein -rerun protein.xtc -nb cpu, i got the same warning as follows
\\
WARNING: Some frames do not contain velocities
Ekin, temperature and pressure are incorrect,
the virial will be incorrect
I am not able to understand why is this happening? I have .trr file also which has the velocities then why it says “some frames do not contain velocities”. Please give your valuable suggestion. Thanks in advance
The .trr file is the only one with velocities. An .xtc file only ever has coordinates.
Then in that case how can i avoid this error? Please suggest some solution sir.
If you need velocities, use the file that contains them (the .trr file). Conversion to .xtc format makes no sense if you need velocities, because you lose them there.
Sir thank you for your response. I already have the protein.tpr file that contains the protein coordinates. Then i used the command
gmx trjconv -s md_0_100.tpr -f md_0_100.trr -o protein.trr and selected (protein) #1 when prompted. As a result of this i have the protein.trr file. Now, when I executed the command
gmx mdrun -deffnm protein -rerun protein.trr -nb cpu I got the following error:
Back Off! I just backed up protein.trr to ./#protein.trr.1#
starting md rerun ‘Protein in water’, reading coordinates from input trajectory ‘protein.trr’
Last frame -1 time 0.000
Program gmx mdrun, VERSION 5.1.4
Source code file: /home/guest/gromacs/gromacs-5.1.4/src/programs/mdrun/md.cpp, line: 685
Fatal error:
Number of atoms in trajectory (-1) does not match the run input file (3372).
Please give some suggestion. Thanks in advance.
Dr. Lemkul please give your valuable suggestion so that I can overcome this problem.
No idea. What does gmx check tell you about the md_0_100.trr and protein.trr files?
Dr. Lemkul I used the command: gmx check -f md_0_100.trr -f2 protein.trr. Lot of output was printed on the terminal and I am pasting the output printed in the beginning and at the very end:
gmx check -f md_0_100.trr -f2 protein.trr
Comparing trajectory files md_0_100.trr and protein.trr
trn version: GMX_trn_file (single precision)
Reading frame 0 time 0.000
natoms (49891 - 3363)
Reading frame 1 time 10.000
natoms (49891 - 3363)
Reading frame 2 time 20.000
natoms (49891 - 3363)
Reading frame 3 time 30.000
natoms (49891 - 3363)
Reading frame 4 time 40.000
natoms (49891 - 3363)
Reading frame 5 time 50.000
natoms (49891 - 3363)
Reading frame 6 time 60.000
natoms (49891 - 3363)
Reading frame 7 time 70.000
natoms (49891 - 3363)
Reading frame 8 time 80.000
natoms (49891 - 3363)
Reading frame 9 time 90.000
natoms (49891 - 3363)
Reading frame 10 time 100.000
natoms (49891 - 3363)
Reading frame 11 time 110.000
natoms (49891 - 3363)
Reading frame 12 time 120.000
natoms (49891 - 3363)
Reading frame 13 time 130.000
natoms (49891 - 3363)
Reading frame 14 time 140.000
natoms (49891 - 3363)
Reading frame 15 time 150.000
natoms (49891 - 3363)
Reading frame 16 time 160.000
natoms (49891 - 3363)
Reading frame 17 time 170.000
natoms (49891 - 3363)
Reading frame 18 time 180.000
natoms (49891 - 3363)
Reading frame 19 time 190.000
natoms (49891 - 3363)
Reading frame 20 time 200.000
natoms (49891 - 3363)
natoms (49891 - 3363)
Reading frame 10000 time 100000.000
natoms (49891 - 3363)
Last frame 10000 time 100000.000
Last frame 10000 time 100000.000
Both files read correctly.
The protein.trr file is obviously smaller as it only contains information related to protein since I generated this by using the command gmx trjconv -s md_0_100.tpr -f md_0_100.trr -o protein.trr and selected (protein) #1. Kindly look into this and give your valuable suggestions sir. Many thanks for your support.
If gmx check reports normal contents of protein.trr then there is some problem in mdrun but you should not be trying to do this kind of calculation with a wildly outdated version of the software (5.1.4) that isn’t even supported any more. If you can reproduce the issue with the latest version, file a bug report.