High RMSD value and Fluctuations

User discussions
1

GROMACS version: 2019.3
GROMACS modification: No

Dear all,

I know that there are lots of topics about high RMSD value, sudden jumps or fluctuations. I’ve been scrolling through the topics and trying to apply several trjconv commands into my system but could not even get a smooth RMSD value and I am not sure if my system has crashed. I have carried out protein+ligand MD simulation and I had smooth RMSD values so far except the last system I set up for my simulation. I couldn’t fix my RMSD result even though I have tried several commands that I have listed out of order.

gmx trjconv -f md_0_1.xtc -s md_0_1.tpr -center -pbc mol

trjconv -s md_0_1.tpr -f md_0_1.xtc -o md_0_1_fit-rot_trans.xtc -ur compact -fit rot+trans

gmx trjconv -s prod_md.tpr -f prod_md.xtc -o prod_md_nopbc_cluster.xtc -pbc cluster

Sadly, these commands did not work, I think I am missing some points since I’m new to GROMACS. Sorry for bringing this topic up again but I got quite confused.

Here is the RMSD result of my system:

Screenshot_3
Screenshot_3695×796 7.54 KB

2

Would anyone suggest any solution of this problem? Any help would be appreciated.

3

What happens when you watch the trajectory? It looks like there is a legitimate, steady increase but with a few dips. What is happening in the frames with the sudden dips?

4

Hi Mr. Jalemkul

First trajectory image belongs to the original .xtc file:

Screenshot_13
Screenshot_131241×946 138 KB

And the second one belongs to the .xtc file that I have modified with above-mentioned commands.

Screenshot_14
Screenshot_141241×946 111 KB

It seems that my ligand is not within the protein and the commands did not work.

5

Both frames appear to be simple PBC issues. What group are you choosing for the trjconv selections related to centering and fitting? You should be choosing the protein, or perhaps just the backbone. If you do, the ligand will wrap to the protein.

6

You can take a look at some code we developed to fix these problems by ensuring assemblies are whole. (GitHub - BartBruininks/mdvwhole: Density based object completion over PBC.). The selection you would need would be a selection which includes both the protein and the ligand. If you try it and you have any issues feel free to get in contact with me.

Cheers,

Bart

7

Hi, @jalemkul
I have same type of issues with my all simulations.
Ligand/peptide shows a lot of movement. Rather vibrating at the docked site, it dances outside the protein periphery.

Is it normal in MD simulations? or there is something wrong in MDP files? or Ligand didn’t bind strongly/close enough in the docked complex?

After

gmx mdrun -deffnm md_0_10

I have recentered and performed rotational and translational fitting

gmx trjconv -s md_0_10.tpr -f md_0_10.xtc -o md_0_10_center.xtc -center -pbc mol -ur compact

gmx trjconv -s md_0_10.tpr -f md_0_10_center.xtc -o md_0_10_fit.xtc -fit rot+trans

Any suggestions/guidance will be appreciated.

Please have a look in the MD movie
https://drive.google.com/file/d/1f75WcRbMkLDoWq9xR7sbQvGsvF1wzHlv/view?usp=share_link

Ligand and the docked sites are in Red and Cyan color respectively.

Please find my MDP files attached
production.mdp (2.4 KB)
eq1-nvt.mdp (2.8 KB)
eq2-npt.mdp (2.7 KB)
ions.mdp (1.0 KB)
minim.mdp (1.1 KB)

8

Your ligand clearly dissociates and it isn’t an artifact. What is the ligand and how did you parametrize it? It looks rather small and very flexible, so it will be extremely sensitive to the quality of the topology you generate for it.

9

@jalemkul Thanks for kind consideration.

Ligand is Gama Amino Butyric Acid (https://pubchem.ncbi.nlm.nih.gov/compound/119).

I used SwissParam Server (https://www.swissparam.ch/) to generate its PDB from mol2 file and then used gmx editconf -f LIG.pdb -o LIG.gro to generate it’s gro file.

out of the scope question:
PDB must be renumbered before MD simulation? Also, chain must not have any missing residues

Thanks

10

Can you post the topology that it generated? That is a simple molecule so the topology should be straightforward.

No.

Correct.

11

Thank you @jalemkul

please find the latest Topol file attached. Also find the other LIG associated files in the drive (https://drive.google.com/file/d/1sfuIDwD_YpjZsN8QZg2aaS1rkbfsrXZy/view?usp=share_link).

topol.top (1.7 MB)

12

The Google Drive link is inaccessible. The information in topol.top is not terribly useful because what I need is gaba_lig.itp.

13

Please find the link for the gaba_lig.itp. Can’t upload here as this portal doesn’t allow us to upload certain files.

14

That topology makes no sense at all. The charges are totally wrong, with most of the C and H atoms bearing zero charge. The charged amino and carboxy groups are over-polarized with respect to typical CHARMM groups so altogether it probably renders the molecule wildly imbalanced and prone to dissociation.

You can generate a topology for this species very easily by analogy from existing CHARMM parameters (by analogy to LYS and GLU, from CHARMM36):

[ atoms ]
; nr type resnr resid atom cgnr charge mass 
   1 NH3  1  LIG N       1 -0.3000  14.0067
   2 CT2  1  LIG CA      2  0.2100  12.0110
   3 CT2  1  LIG CB      3 -0.1800  12.0110
   4 CT2  1  LIG CG      4 -0.2800  12.0110
   5 CC   1  LIG CD      5  0.6200  12.0110
   6 OC   1  LIG OE1     6 -0.7600  15.9994
   7 OC   1  LIG OE2     7 -0.7600  15.9994
   8 HA2  1  LIG HA2     8  0.0500   1.0079
   9 HA2  1  LIG HA3     9  0.0500   1.0079
  10 HA2  1  LIG HB2    10  0.0900   1.0079
  11 HA2  1  LIG HB3    11  0.0900   1.0079
  12 HA2  1  LIG HG2    12  0.0900   1.0079
  13 HA2  1  LIG HG3    13  0.0900   1.0079
  14 HC   1  LIG H1     14  0.3300   1.0079
  15 HC   1  LIG H2     15  0.3300   1.0079
  16 HC   1  LIG H3     16  0.3300   1.0079
15

@jalemkul Thanks for your kind response

I have generated ligand topology as per your protocol using ‘https://cgenff.umaryland.edu/userAccount/userLogin.php’ and cgenff_charmm2gmx_py3_nx2.py script.

it still disassociates after some time. Total MD run is 60 ns.
Please find the associated files in the drive.
https://drive.google.com/drive/folders/1iqwPVY1OQ_3Ygfr1p-1K_enPioIoVPTq?usp=share_link

16

Why are you protonating the carboxylate group? Only under extremely acidic conditions would this be the case.

17

Thanks for your kind consideration. I did it as a part of molecular docking protocol where I have to add hydrogen. Not sure though if I can avoid protonation of terminal carboxylate while adding polar hydrogen.
Please guide me over it.
Thanks

18

You need to model real conditions, and GABA is a zwitterion at physiological pH. So you need to protonate it properly, or the model lacks value.

19

Thank you so much :) Can I use pbd2pqr module for the protonation?

20

Sure, there are any number of tools that can accomplish this task. Use whichever one you’re most comfortable with.