Instability of protein


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

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.

What does this mean?

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 - 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)

*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.

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).

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?

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.

Here’s my script

### 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

### GROMP ###
gmx grompp -f input/ions.mdp -c 03_${name}_solv.pdb -p -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 -pname NA -nname CL -conc 0.15 -neutral

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

### MDRUN ###
gmx mdrun -s 05_${name}_em.tpr -c 05_${name}_em.pdb -mp -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 -o 06_${name}_nvt.tpr -n index.ndx
gmx mdrun -s 06_${name}_nvt.tpr -c 06_${name}_nvt.pdb -mp -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 -o 07_${name}_npt.tpr -n index.ndx
gmx mdrun -s 07_${name}_npt.tpr -c 07_${name}_npt.pdb -mp -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 -o 08_${name}_md.tpr -n index.ndx
gmx mdrun -deffnm 08_${name}_md -c 08_${name}_md.pdb -mp

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

This should set an appropriate periodic distance.

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.

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

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

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

Skærmbillede 2020-05-16 kl. 13.54.02

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

Group 24 is my protein and ligand; the dimer

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.

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.

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

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

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

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.

This is part of the .gro file.

I will try with editconf. Thanks