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?
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.
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.
pdb2gmx, 2) editconf, 3) solvate, and grompp, 4) genion, and grompp, 5) mdrun minim 6) NVT, 7) NPT, 8) mdrun.
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.
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.
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
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.
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.
