Structural Realignment

User discussions
1

GROMACS version: 2020.2

Hello All,

I am looking to realign my protein after my simulation before calculating RMSD. From my understanding, I would need a tpr file along with the centered xtc file and may use a python script to carry out structural realignment. Does anyone have a script that can do this in addition to more information on this subject?

Thank you!

2 Likes
2

Assuming I understood your question: gmx rms command involves two stages: first select a group for least-squares alignment (fitting) to the reference structure, then calculates rmsd.

3

Thank you. Yes, the selected group is aligned to the reference structure first (tpr), and then rmsd is calculated. My question is more so about a correction to the output xtc file to structurally align it with the reference structure before moving to gmx rms.

In the analysis section of the Lysozyme in water tutorial found here: (http://www.mdtutorials.com/gmx/lysozyme/05_EM.html) we recenter the protein in the output xtc file in the form of a noPBC.xtc file and then we calculate the rmsd in the two steps that you have described.

But based on my current understanding, re-centered and structurally aligned are two different things. As in, if the protein were re-centered then it could still be slightly out of alignment with the reference structure which would therefore influence the rmsd. I am looking to first align the reference structure and the output structure so that they are perfectly superimposed before least-squares fitting and before calculating rmsd.

Any thoughts/advice is greatly appreciated, thank you in advance!

4

This is what gmx trjconv -fit is for.

1 Like
5

I see, is there a certain selection that I have to make: none, rot+trans, rotxy+transxy, translation, transxy, progressive? Or simply include -fit in the command?

6

I tried the follow steps to implement -fit:

  1. gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center
  2. gmx trjconv -s md_0_1.tpr -f md_0_1_noPBC.xtc -o md_0_1_align.xtc -fit rot+trans
  3. gmx rms -s md_0_1.tpr -f md_0_1_align.xtc -o rmsd.xvg -tu ns

But, this gave the error:
Molecule in topology has atom numbers below and above natoms (8748).
You are probably trying to use a trajectory which does not match the first
8748 atoms of the run input file. How can I use gmx convert-tpr to make an appropriate input tpr file?

7

The contents of the .tpr have to match the trajectory. So if you only saved a subset of coordinates in the .xtc file, you need to have that same subset in the .tpr file.

8

I see, thank you. This segment is saved from one part of my simulaiton in md_0_1_align.xtc which came from md_0_1_noPBC.xtc. To do so can I use the command:
gmx convert-tpr -s md_0_1.tpr -f md_0_1_align.xtc -o md_0_1_align.tpr ?

9

No, because convert-tpr does not take -f or any coordinate/trajectory format. All you need is -s and -o and you can name the output whatever you like (doesn’t have anything to do with the coordinate manipulations in the trajectory, really).

10

Right, thank you for the clarification. So then:
gmx convert-tpr -s md_0_1.tpr -o md_0_1_align.tpr

Would give me a matching run input file?

11

Try it and see :)

12

It seems that this leads to the same error, but I did select 0) system for convert-tpr

13

I do see that I could use -nsteps to change the number of nsteps. Maybe to match the number covered in the xtc file?

14

You need to select the group that corresponds to what is in the .xtc file. Your .tpr already has the full system. By selecting “System” again, you are accomplishing nothing.

This is only relevant for mdrun and is not related to your problem.

15

Okay, I see. Thank you for the correction. So say I am interested in the rmsd of the backbone in this particular case. I carrried out the following steps to write the rmsd.xvg file:

  1. gmx trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_noPBC.xtc -pbc mol -center (selecting 1 - Protein and 0 - system to re-center)
  2. gmx trjconv -s md_0_1.tpr -f md_0_1_noPBC.xtc -o md_0_1_align.xtc -fit rot+trans (selecting 4 - Backbone for LSQ and 4 - Backbone for output)
  3. gmx convert-tpr -s md_0_1.tpr -o align.tpr (selecting 4 - Backbone)
  4. gmx rms -s align.tpr -f md_0_1_align.xtc -o rmsd.xvg -tu ns (selecting 4 - Backbone for LSQ and 4 - Backbone for output)

Which seems to have produced a reasonable plot. Did I make any error in my procedure?

16

There’s nothing “wrong” but I see no reason to only select the backbone for output in your trajectory. If you adopt this convention, you have to go back and run trjconv every time you want to analyze something else about the protein. Stripping water is standard procedure and I recommend just saving the whole protein and fitting to whatever group you like later.

17

Thank you so much. This has helped tremendously!