# Instability of protein

**URL:** <https://gromacs.bioexcel.eu/t/instability-of-protein/120>\
**Category:** User discussions\
**Created:** [May 12, 2020, 6:11pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120 "2020-05-12T18:11:44Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [May 12, 2020, 6:11pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/1 "2020-05-12T18:11:44Z")

</div>

Hi,

I have some troubles with my simulations. I get errors indicating an unstable system.  
_Fatal error:_  
_1 particles communicated to PME node 1 are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension x._  
_This usually means that your system is not well equilibrated._  
_Fatal error:_  
_A charge group moved too far between two domain decomposition steps_  
_This usually means that your system is not well equilibrated_

I have tried to increase equilibration steps, and tried to lower the timesteps, but I keep getting told the system is not well equilibrated.  
Additionally, one of my conformations cannot be placed within its simulation box no matter which pbc i try to use or reference files.

How should I handle the instability of my protein structures, and how can I figure out what is wrong with my structures?

Best regards

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [May 12, 2020, 6:41pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/2 "2020-05-12T18:41:31Z")

</div>

Please provide a full description of what your system is and what steps you took to prepare it, including the output force and energy from energy minimization. Full .mdp files would be helpful.

> [@Irene](#):
>
> Additionally, one of my conformations cannot be placed within its simulation box no matter which pbc i try to use or reference files.

What does this mean?

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [May 12, 2020, 6:58pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/3 "2020-05-12T18:58:02Z")

</div>

I cannot upload my files as I am a new user here.

I am simulating a dimer with a ligand each. I have three different conformations (aa, ab, bb). It’s a system of about 600 residues and a FAD domain each. Before using gromacs, the system has gone through pdb2pqr.

1. pdb2gmx, 2) editconf, 3) solvate, and grompp, 4) genion, and grompp, 5) mdrun minim 6) NVT, 7) NPT, 8) mdrun.

After minim:  
Energy Average Err.Est. RMSD Tot-Drift  
Bond 2153.36 680 4724.93 -3998.66 (kJ/mol)  
Angle 5411.75 260 1187.67 -1571.05 (kJ/mol)  
Proper Dih. 22477.7 360 778.196 -2362.12 (kJ/mol)  
Improper Dih. 351.045 18 50.9008 -120.905 (kJ/mol)  
LJ-14 9069.39 320 900.356 -2009.49 (kJ/mol)  
Coulomb-14 111431 310 742.14 -1995.44 (kJ/mol)  
LJ (SR) 626309 160000 3.16643e+06 -928802 (kJ/mol)  
Coulomb (SR) -3.96935e+06 60000 133707 -402631 (kJ/mol)  
Coul. recip. 21966.2 4400 15334.5 -27734.4 (kJ/mol)  
Potential -3.17018e+06 220000 3.21219e+06 -1.37122e+06 (kJ/mol)  
Pressure -6744.13 110 238.863 -716.279 (bar)

minim.mdp  
_; minim.mdp - used as input into grompp to generate em.tpr_  
_integrator = steep ; Algorithm (steep = steepest descent minimization)_  
_emtol = 1000.0 ; Stop minimization when the maximum force \< 1000.0 kJ/mol/nm_  
_emstep = 0.001 ; Energy step size_  
_nsteps = 50000 ; Maximum number of (minimization) steps to perform_

_; Parameters describing how to find the neighbors of each atom and how to calculate the interactions_  
_nstlist = 1 ; Frequency to update the neighbor list and long range forces_  
_cutoff-scheme = Verlet_  
_ns\_type = grid ; Method to determine neighbor list (simple, grid)_  
_coulombtype = PME ; Treatment of long range electrostatic interactions_  
_rcoulomb = 1.05 ; Short-range electrostatic cut-off_  
_rvdw = 1.05 ; Short-range Van der Waals cut-off_  
_pbc = xyz ; Periodic Boundary Conditions (yes/no)_

npt  
\*title = OPLS Lysozyme NPT equilibration \*  
_define = ; position restrain the protein_  
_; Run parameters_  
_integrator = md ; leap-frog integrator_  
_nsteps = 10000 ; 2 \* 100000 = 200 ps_  
_dt = 0.002 ; 2 fs_  
_; Output control_  
_nstxout = 500 ; save coordinates every 1.0 ps_  
_nstvout = 500 ; save velocities every 1.0 ps_  
_nstenergy = 500 ; save energies every 1.0 ps_  
_nstlog = 500 ; update log file every 1.0 ps_  
_; Bond parameters_  
\*continuation = yes ; Restarting after NVT \*  
\*constraint\_algorithm = lincs ; holonomic constraints \*  
_constraints = h-bonds ; all bonds (even heavy atom-H bonds) constrained_  
_lincs\_iter = 1 ; accuracy of LINCS_  
_lincs\_order = 4 ; also related to accuracy_  
_; Neighborsearching_  
_cutoff-scheme = Verlet_  
_ns\_type = grid ; search neighboring grid cells_  
_nstlist = 20 ; 20 fs, largely irrelevant with Verlet scheme_  
_rcoulomb = 1.05 ; short-range electrostatic cutoff (in nm)_  
_rvdw = 1.05 ; short-range van der Waals cutoff (in nm)_  
_; Electrostatics_  
_coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics_  
_pme\_order = 4 ; cubic interpolation_  
_fourierspacing = 0.16 ; grid spacing for FFT_  
_; Temperature coupling is on_  
_tcoupl = V-rescale ; modified Berendsen thermostat_  
_tc-grps = Protein\_FAD Water\_Ion ; two coupling groups - more accurate_  
_tau\_t = 0.1 0.1 ; time constant, in ps_  
_ref\_t = 310.15 310.15 ; reference temperature, one for each group, in K_  
_; Pressure coupling is on_  
_pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT_  
_pcoupltype = isotropic ; uniform scaling of box vectors_  
_tau\_p = 2.0 ; time constant, in ps_  
_ref\_p = 1.0 ; reference pressure, in bar_  
_compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1_  
_refcoord\_scaling = com_  
_; Periodic boundary conditions_  
_pbc = xyz ; 3-D PBC_  
_; Dispersion correction_  
_DispCorr = EnerPres ; account for cut-off vdW scheme_  
_; Velocity generation_  
_gen\_vel = no ; Velocity generation is off_

For the pbc part, then the protein is split in different parts of the simulation box due to moving out of the box. It doesn’t help using nojump, mol, res, center, etc. Nothing can get the protein back together as one.

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [May 12, 2020, 7:14pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/4 "2020-05-12T19:14:50Z")

</div>

> [@Irene](#):
>
> I cannot upload my files as I am a new user here.
> 
> I am simulating a dimer with a ligand each. I have three different conformations (aa, ab, bb). It’s a system of about 600 residues and a FAD domain each. Before using gromacs, the system has gone through pdb2pqr.
> 
> 1. pdb2gmx, 2) editconf, 3) solvate, and grompp, 4) genion, and grompp, 5) mdrun minim 6) NVT, 7) NPT, 8) mdrun.

Please post actual commands. We’re all familiar with general workflows but specifics of what you have done are important; the `editconf` command is particularly important given your last point (see below).

> [@](#):
>
> After minim:  
> Energy Average Err.Est. RMSD Tot-Drift  
> Bond 2153.36 680 4724.93 -3998.66 (kJ/mol)  
> Angle 5411.75 260 1187.67 -1571.05 (kJ/mol)  
> Proper Dih. 22477.7 360 778.196 -2362.12 (kJ/mol)  
> Improper Dih. 351.045 18 50.9008 -120.905 (kJ/mol)  
> LJ-14 9069.39 320 900.356 -2009.49 (kJ/mol)  
> Coulomb-14 111431 310 742.14 -1995.44 (kJ/mol)  
> LJ (SR) 626309 160000 3.16643e+06 -928802 (kJ/mol)  
> Coulomb (SR) -3.96935e+06 60000 133707 -402631 (kJ/mol)  
> Coul. recip. 21966.2 4400 15334.5 -27734.4 (kJ/mol)  
> Potential -3.17018e+06 220000 3.21219e+06 -1.37122e+06 (kJ/mol)  
> Pressure -6744.13 110 238.863 -716.279 (bar)

Averages from energy minimization are useless because the system is inherently being forced to change. Did you achieve a maximum force below the specified tolerance?

> [@](#):
>
> For the pbc part, then the protein is split in different parts of the simulation box due to moving out of the box. It doesn’t help using nojump, mol, res, center, etc. Nothing can get the protein back together as one.

If there is no way to reconstruct your protein, this implies that your box is too small and your buffer set using `editconf -d` was too small. You could have clashes between your periodic images that are causing the instability upon introducing the barostat. Please provide your actual `editconf` command. Images (before and after minimization, after NVT and NPT) would also be useful.

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [May 12, 2020, 7:26pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/5 "2020-05-12T19:26:29Z")

</div>

Here’s my script  
name=“free\_Ec\_FOFR\_1”

```
### PDB2GMX ###
{ echo "1"; echo "1"; } | gmx pdb2gmx -f Ec_FOFR_propka.pdb -o 01_${name}.pdb

### EDITCONF ###
gmx editconf -f 01_${name}.pdb -o 02_${name}_newbox.pdb -c -d 1.7 -bt dodecahedron # I have tried -d 1.5 1.7 2.0 2.2

### SOLVATE ###
gmx solvate -cp 02_${name}_newbox.pdb -cs spc216.gro -o 03_${name}_solv.pdb -p topol.top

### GROMP ###
gmx grompp -f input/ions.mdp -c 03_${name}_solv.pdb -p topol.top -o 04_${name}_ions.tpr -maxwarn 1 # warning that the system is charged and on PME.

### GENION ###
echo "15" | gmx genion -s 04_${name}_ions.tpr -o 04_${name}_ions.pdb -p topol.top -pname NA -nname CL -conc 0.15 -neutral

### GROMPP ###
gmx grompp -f input/minim.mdp -c 04_${name}_ions.pdb -p topol.top -o 05_${name}_em.tpr

### MDRUN ###
gmx mdrun -s 05_${name}_em.tpr -c 05_${name}_em.pdb -mp topol.top -o 05_${name}_em.trr -g 05_${name}_em.log -nt 4 -e 05_${name}_em.edr

### ENERGY ###
echo "10" | gmx energy -f 05_${name}_em.edr -o 05_${name}_potential.xvg

### MAKE NDX ###
{ echo "1 | 20"; echo " 16 | 19"; echo "q"; } | gmx make_ndx -f 05_${name}_em.pdb

### NVT ###
gmx grompp -f input/nvt.mdp -c 05_${name}_em.pdb -r 05_${name}_em.pdb -p topol.top -o 06_${name}_nvt.tpr -n index.ndx
gmx mdrun -s 06_${name}_nvt.tpr -c 06_${name}_nvt.pdb -mp topol.top -o 06_${name}_nvt.trr -e 06_${name}_nvt.edr -g 06_${name}_nvt.log -cpo 06_${name}_nvt.cpt
echo "15" | gmx energy -f 06_${name}_nvt.edr -o 06_${name}_nvt_temperature.xvg

### NPT ###
gmx grompp -f input/npt.mdp -c 06_${name}_nvt.pdb -r 06_${name}_nvt.pdb -t 06_${name}_nvt.cpt -p topol.top -o 07_${name}_npt.tpr -n index.ndx
gmx mdrun -s 07_${name}_npt.tpr -c 07_${name}_npt.pdb -mp topol.top -o 07_${name}_npt.trr -e 07_${name}_npt.edr -g 07_${name}_npt.log -cpo 07_${name}_npt.cpt
echo "17" | gmx energy -f 07_${name}_npt.edr -o 07_${name}_npt_pressure.xvg

### MD ###
gmx grompp -f input/md.mdp -c 07_${name}_npt.pdb -r 07_${name}_npt.pdb -t 07_${name}_npt.cpt -p topol.top -o 08_${name}_md.tpr -n index.ndx
gmx mdrun -deffnm 08_${name}_md -c 08_${name}_md.pdb -mp topol.top

```

The maximum potential force is the start of 4.0423e+11 kJ/mol and the minimum is -2717776 kJ/mol.

I will try to make the box 2.5

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [May 13, 2020, 3:26pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/8 "2020-05-13T15:26:00Z")

</div>

> [@Irene](#):
>
> > gmx editconf -f 01\_{name}.pdb -o 02\_{name}\_newbox.pdb -c -d 1.7 -bt dodecahedron # I have tried -d 1.5 1.7 2.0 2.2

This should set an appropriate periodic distance.

> [@](#):
>
> I will try to make the box 2.5

There is no reason to do that.

A dodecahedron is automatically represented as a compact, triclinic cell when the coordinates are visualized directly. Using `gmx trjconv -pbc mol -ur compact` should account for any visualization artifacts and will reconstruct the dodecahedral shape of the unit cell.

If you are still having problems, please upload images of the system.

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [May 16, 2020, 11:55am UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/9 "2020-05-16T11:55:58Z")

</div>

Sorry for my late feedback, I couldn’t access my files as the server was being maintained.  
I have tried multiple ways;  
`echo 24 24 | gmx trjconv -s $tpr -f $xtc -n $ndx -center -pbc mol -ur compact -b 90000 -dump 90005 -o L_FOFR_frame_90000_1.pdb`

 ![Skærmbillede 2020-05-16 kl. 13.40.11](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/1X/a801e7484ba6e27fc399d90de480a6599857ad8d.png)

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [May 16, 2020, 11:56am UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/10 "2020-05-16T11:56:40Z")

</div>

and this way  
`echo 24 24 | gmx trjconv -s $pdb -f $xtc -n $ndx -center -pbc nojump -ur compact -b 90000 -dump 90005 -o L_FOFR_frame_90000_1.pdb`

 ![Skærmbillede 2020-05-16 kl. 12.30.50](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/1X/337f4cc4cd3a69bbe601514ab3a24b5dfb8c6b47.png)

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [May 16, 2020, 11:57am UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/11 "2020-05-16T11:57:21Z")

</div>

Neither can set my protein back. It was supposed to look like

![Skærmbillede 2020-05-16 kl. 13.54.02](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/1X/422c969e4ce5b5d56d3cdcd569cca6d0dac918a7.png)

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [May 17, 2020, 11:03am UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/12 "2020-05-17T11:03:21Z")

</div>

What is group 24 in your index file? Is the protein a single chain or a dimer?

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [May 17, 2020, 11:18am UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/13 "2020-05-17T11:18:12Z")

</div>

Group 24 is my protein and ligand; the dimer

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [May 17, 2020, 2:32pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/14 "2020-05-17T14:32:02Z")

</div>

Center on one chain of the dimer with `trjconv -center`. The output you are getting is centered based on the center of geometry of the selection being coincident with the center of the box, but that’s inconvenient for visualization.

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [May 17, 2020, 2:34pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/15 "2020-05-17T14:34:08Z")

</div>

so I should make an index of just one of the monomers in the dimer and center it to that? The protein is said to only have one chain after simulation.

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [May 17, 2020, 2:46pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/16 "2020-05-17T14:46:24Z")

</div>

> [@Irene](#):
>
> so I should make an index of just one of the monomers in the dimer and center it to that?

Yes. This is almost always what one needs to do when there are multiple species (even a single protein and ligand) involved.

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [May 17, 2020, 2:48pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/17 "2020-05-17T14:48:05Z")

</div>

how do i make the index? both monomers go from 1-310, but are ‘chainless’?

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [May 17, 2020, 2:59pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/18 "2020-05-17T14:59:02Z")

</div>

There should not be duplicate residue numbers in a `.gro` file so you can use non-redundant residue numbers or specify the chain you want by atom number with `make_ndx`.

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [May 17, 2020, 3:29pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/19 "2020-05-17T15:29:46Z")

</div>

> [@jalemkul](#):
>
> There should not be duplicate residue numbers in a `.gro` file so you can use non-redundant residue numbers or specify the chain you want by atom number with `make_ndx` .

I have generated a .gro file. the monomers are named 1-310 in the file as well. Can I make a selection of atom 1-4822? Then I can get one chain like that.

I have tried various versions of “a 1-4821” but the group ends up empty

---

<div class="post-metadata">

**Author:** ![jalemkul](https://dub1.discourse-cdn.com/flex017/user_avatar/gromacs.bioexcel.eu/jalemkul/32/18_2.png) [@jalemkul](https://gromacs.bioexcel.eu/u/jalemkul)\
**Post date:** [May 17, 2020, 7:20pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/20 "2020-05-17T19:20:12Z")

</div>

If you get an empty group, then your atoms aren’t numbered from 1, but everything you’ve described is odd and shouldn’t be the case with coordinate files produced by GROMACS utilities. Everything should be numbered from 1 and there shouldn’t be any repeated residue numbers. Try `editconf -resnr 1` and work with the resulting file.

---

<div class="post-metadata">

**Author:** ![Irene](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@Irene](https://gromacs.bioexcel.eu/u/Irene)\
**Post date:** [May 17, 2020, 8:02pm UTC](https://gromacs.bioexcel.eu/t/instability-of-protein/120/21 "2020-05-17T20:02:08Z")

</div>

![Skærmbillede 2020-05-17 kl. 22.00.40](https://europe1.discourse-cdn.com/flex017/uploads/bioexcel1/original/1X/030b7f2b00f2232c33d7f6f95a81f06cb38adaf0.png)

This is part of the .gro file.

I will try with editconf. Thanks