GROMACS version: 2020.1-Ubuntu-2020.1-1
GROMACS modification: No
Here post your question:
So I am doing a simulation of protein solvated in gas - CO2 RCSB PDB - 5WSK: Structure of Ribulose-1,5-bisphosphate carboxylase/oxygenase from wheat. I took topology and PDB of this gas from Virtual Sites. I took the PDB model from PDB and then removed water and all non-standard residues.
Then I ran these gmx commands using the files provided. (I used OPLS/AA force field)
gmx pdb2gmx -f rubisco_mfix.pdb -o 1.gro -p topol.top -water none
gmx editconf -f 1.gro -o box.gro -c -d 3.0 -bt dodecahedron
gmx insert-molecules -f box.gro -ci co2.pdb -nmol 40 -try 20 -o boxco2.gro
Then I included topol.itp (which is itp of co2 with ff included) in topology (topol.top) and I added CO2 molecules to the system.
I didn’t neutralize this system, because I think that you don’t need to neutralize gas solvated systems (If I have to, how can I do it?) (With all gmx grommp I was getting a Warning about this, so I used -maxwarn 1)
Then I ran these gmx commands (mdp files included) :
gmx grompp -f minimize.mdp -c water_ions.gro -p topol.top -o energymin.tpr
gmx mdrun -v -s energymin.tpr -deffnm em
gmx grompp -f nvt.mdp -c em.gro -r em.gro -p topol.top -o nvt.tpr
gmx mdrun -deffnm nvt
All of these runs were successful. Then I attempted NPT run and it returned this error :
gmx grompp -f npt.mdp -c nvt.gro -r nvt.gro -t nvt.cpt -p topol.top -o npt.tpr
gmx mdrun -deffnm npt
Step 0, time 0 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 1.241085, max 56.916912 (between atoms 70025 and 70026)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
69980 69981 90.0 0.2132 11.0419 0.2132
69985 69986 90.0 0.2132 6.3848 0.2132
69990 69991 93.1 0.2132 10.8002 0.2132
69995 69996 99.1 0.2132 9.3370 0.2132
70000 70001 100.8 0.2132 8.7460 0.2132
70010 70011 113.2 0.2132 7.8905 0.2132
70015 70016 90.0 0.2132 11.0835 0.2132
70020 70021 90.0 0.2132 8.9943 0.2132
70025 70026 90.0 0.2132 12.3463 0.2132
70030 70031 90.0 0.2132 9.1664 0.2132
70035 70036 90.0 0.2132 5.8325 0.2132
70040 70041 90.0 0.2132 3.0451 0.2132
70045 70046 90.0 0.2132 7.3950 0.2132
70050 70051 90.0 0.2132 4.0567 0.2132
70055 70056 90.0 0.2132 11.7571 0.2132
70060 70061 90.0 0.2132 7.8807 0.2132
70070 70071 90.0 0.2132 9.5708 0.2132
70075 70076 90.0 0.2132 6.0452 0.2132
70080 70081 90.0 0.2132 7.4494 0.2132
70085 70086 90.0 0.2132 2.7008 0.2132
70090 70091 90.0 0.2132 5.1878 0.2132
70095 70096 94.6 0.2132 10.1661 0.2132
70100 70101 90.0 0.2132 5.2383 0.2132
70105 70106 114.0 0.2132 7.9965 0.2132
70110 70111 90.0 0.2132 2.7225 0.2132
70115 70116 100.1 0.2132 10.9836 0.2132
70120 70121 90.0 0.2132 4.2884 0.2132
70125 70126 90.0 0.2132 5.1484 0.2132
70130 70131 90.0 0.2132 7.7973 0.2132
70135 70136 90.0 0.2132 9.1673 0.2132
70140 70141 90.0 0.2132 7.0320 0.2132
70145 70146 90.0 0.2132 6.3205 0.2132
70150 70151 90.0 0.2132 7.4730 0.2132
70155 70156 90.0 0.2132 11.0906 0.2132
70160 70161 90.0 0.2132 2.3694 0.2132
70165 70166 90.0 0.2132 10.8323 0.2132
70170 70171 90.0 0.2132 9.2923 0.2132
70175 70176 90.0 0.2132 8.7614 0.2132
70175 70176 90.0 0.2132 8.7614 0.2132
70175 70176 90.0 0.2132 8.7614 0.2132
70175 70176 90.0 0.2132 8.7614 0.2132
Atoms given in this message are atoms of my CO2 Virtual Sites (checked in .gro file)Then it produced other LINCS errors in further steps.
My question is - how do I solve this problem? I am pretty sure that problem lies in my CO2 virtual sites model.
; minimize.mdp
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization 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 ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = PME
rvdw = 1.0
rlist = 1.0
rcoulomb = 1.0
fourierspacing = 0.12
pme-order = 4
ewald-rtol = 1e-5
;nvt.mdp
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.001 ; 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 = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; 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 Non-Protein ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 300 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
;npt.mdp
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 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 ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; 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 Non-Protein ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 300 300 ; 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
; Velocity generation
gen_vel = no ; Velocity generation is off
;topol.itp
#include "oplsaa.ff/forcefield.itp"
; Here we introduce a special atom type to be used for each mass center
[ atomtypes ]
; name bond_type mass charge ptype sigma epsilon
MCO MCO 0.000 0.000 A 0.000 0.000
[ moleculetype ]
; name nrexcl
CO2 2
[ atoms ]
; nr type resnr residue atom cgnr charge mass typeB chargeB massB
1 opls_272 1 CO2 O1 1 -0.350 0.0000
2 opls_271 1 CO2 C 1 0.700 0.0000
3 opls_272 1 CO2 O2 1 -0.350 0.0000
4 MCO 1 CO2 M1 1 0.000 22.0049
5 MCO 1 CO2 M2 1 0.000 22.0049
[ constraints ]
; There are no bonds in this system
; Instead, we fix the distance between the mass centers such that
; the virtual sites can be reconstructed
4 5 1 0.213173
[ virtual_sites2 ]
; the M--O distance is 0.125 - 0.1065865 = 0.0184135
; the M--M distance is 0.213173
; therefore, the fraction of the distance along the M--M length is (0.213173+0.0184135)/0.213173 = 1.0851116
; thus placing the virtual O sites beyond the M--M distance
; site ai aj funct a
1 4 5 1 1.0851116 ; relative to mass center 4, extends beyond mass center 5
2 4 5 1 0.5000 ; right in the middle
3 5 4 1 1.0851116 ; as in the case of site 1