[gmx-users] Simulation crashed - Large VCM, Pressure scaling more than 1%, Bond length not finite

[zeineb SI CHAIB]

Dear GMX users,

I am running a Coarse-Grained simulation using Martini Force Field. My system consists of a Protein inserted in a lipid bilayer containing POPC, POPE, and CHOL molecules, in the presence of Water and ions.

I followed the conventional steps to prepare my system for the production run:
1- Generating Coarse-Grained structure and topology files
2- Insertion of the protein in the lipid bilayer
3- Adding neutralizing counterions
4- Energy minimization with backbone positions restraints for 1ns.
5- NVT equilibration with backbone positions restraints for 50 ns with v-rescale.
6- NPT equilibration with backbone positions restraints for 50 ns with Parrinello-Rahman barostat.

The simulation runs for ~ 2 microseconds before crashing:
Large VCM (group Protein_POPC_POPE_CHOL): 965789.43750, -1086946.50000, 135001.78125, Temp-cm:  9.21574e+14
Large VCM(group W_ION): -396609.00000, 446514.28125, -55416.26953, Temp-cm:  2.38187e+13

Step 123438101  Warning: Pressure scaling more than 1%. This may mean your system is not yet equilibrated. Use of Parrinello-Rahman pressure coupling during equilibration can lead to simulation instability and is discouraged.

Program: gmx mdrun, version 2018.3
Source file: src/gromacs/mdlib/clincs.cpp (line 2252)
Fatal error: Bond length not finite.

I used the following parameters in the mdp file (They were recommended by Martini People):
integrator               = md
dt                             = 0.02
nsteps                     = 500000000

nstxout                   = 100
nstvout                   = 100
nstfout                    = 0
nstlog                     = 1000
nstenergy              = 100
nstxout-compressed         = 1000
compressed-x-precision   = 100

continuation           = yes

cutoff-scheme            = Verlet
nstlist                           = 20
ns_type                       = grid
pbc                               = xyz
verlet-buffer-tolerance     = 0.005

coulombtype        = reaction-field
rcoulomb              = 1.1
epsilon_r              = 15
epsilon_rf             = 0
vdw_type             = cutoff
vdw-modifier       = Potential-shift-verlet
rvdw                      = 1.1

tcoupl             = v-rescale
tc-grps            = Protein POPC_POPE_CHOL W_ION
tau_t               = 1.0  1.0 1.0
ref_t                = 315 315 315

Pcoupl                    = parrinello-rahman
Pcoupltype            = semiisotropic
tau_p                      = 12.0
compressibility     = 3e-4  3e-4  3e-4
ref_p                       = 1.0  1.0  1.0

gen_vel                 = no
gen_temp            = 315
gen_seed             = 473529

constraints                       = none
constraint_algorithm     = Lincs

nstcomm           = 100
comm-grps       = Protein_POPC_POPE_CHOL W_ION

I had a look on the GROMACS mailing list for similar problems and probable causes:

In one post they said that the Large VCM might be related to an un-equilibrated system but in my case, the system was equilibrated for 50ns (NVT equilibration) and another 50ns for the NPT equilibration so I don't think that this is the case here.

Another suggestion for the VCM thing is that it might be related to the force field and the topology. Since the last 5 residues of my protein have not been solved I thought that it's better to keep the actual C-terminal residue neutral. I would have liked to use the capping residues but they are not parameterized in Martini force field.  Maybe this is why the simulation is crashing? But it runs for more than 2 microseconds before it crashed!
I would like to have your opinion about this point, please.

The error message is really confusing and I can't really understand it : VCM/Pressure scaling more that 1%/Bond length not finite.

I would really appreciate your help on this.

Zeineb