[gmx-users] Simulation crashed, fatal error: Bond length not finite and warning: Pressure scaling more than 1%.

[zeineb SI CHAIB]

Dear GMX users,


I'm running a coarse-grained simulation of a homo-dimer in a membrane composed of POPC, POPE, and CHOL (31%, 41%, 28% respectively), using MARTINI force field and GROMACS software.


I followed the usual steps with 1ns minimization, 50 ns NVT equilibration followed by 50ns NPT equilibration.


After running 2.5μs of simulation on a cluster, the system crashed with a fatal error:


Step 108657121  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.


Fatal error:

Bond length not finite.


When I analyzed the pressure and the temperature they seem OK. Pressure average = 1.04 bar and Temperature average 314.835 K


I don't know what I'm missing and can't diagnose the problem.


Any help, please?


Thank you in advance for your help and consideration.




NB: I used the following MDP parameters for the production run (They are the optimal parameters to run with MARTINI FF):


; TIMESTEP IN MARTINI

; Most simulations are numerically stable with dt=40 fs,

; however better energy conservation is achieved using a

; 20-30 fs time step.

; Time steps smaller than 20 fs are not required unless specifically stated in the itp file.


integrator              = md

dt                            = 0.02

nsteps                    = 50000000


nstxout                  = 100

nstvout                  = 100

nstfout                  = 0

nstlog                    = 1000

nstenergy                = 100

nstxout-compressed       = 1000

compressed-x-precision   = 100


continuation           = yes     ; Restarting after NPT


; NEIGHBOR LIST and MARTINI

; To achieve faster simulations in combination with the Verlet-neighbor list

; scheme, Martini can be simulated with a straight cutoff. In order to

; do so, the cutoff distance is reduced 1.1 nm.

; The Verlet neighbor list scheme will automatically choose a proper neighbor list

; length, based on a energy drift tolerance.

;

; Coulomb interactions can alternatively be treated using a reaction-field,

; giving slightly better properties.

; Please realize that electrostatic interactions in the Martini model are

; not considered to be very accurate, to begin with, especially as the

; screening in the system is set to be uniform across the system with

; a screening constant of 15. When using PME, please make sure your

; system properties are still reasonable.


cutoff-scheme            = Verlet

nstlist                  = 20

ns_type                  = grid

pbc                      = xyz

verlet-buffer-tolerance  = 0.005


coulombtype              = reaction-field

rcoulomb                 = 1.1

epsilon_r                = 15 ; 2.5 (with polarizable water)

epsilon_rf               = 0

vdw_type                 = cutoff

vdw-modifier             = Potential-shift-verlet

rvdw                     = 1.1


; MARTINI and TEMPERATURE/PRESSURE

; Good temperature control can be achieved with the V-rescale

; thermostat using a coupling constant of the order of 1 ps. Even better
; temperature control can be achieved by reducing the temperature coupling

; constant to 0.1 ps, although with such tight coupling (approaching

; the time step) one can no longer speak of a weak-coupling scheme.

; We therefore recommend a coupling time constant of at least 0.5 ps.

; The Berendsen thermostat is less suited since it does not give

; a well described thermodynamic ensemble.

;

; Pressure can be controlled with the Parrinello-Rahman barostat,

; with a coupling constant in the range 4-8 ps and typical compressibility

; in the order of 10e-4 - 10e-5 bar-1. Note that, for equilibration purposes,

; the Berendsen barostat probably gives better results, as the Parrinello-

; Rahman is prone to oscillating behaviour. For bilayer systems the pressure

; coupling should be done semiisotropic.


tcoupl              = v-rescale

tc-grps             = Protein POPC_POPE_CHOL W_ION

tau_t               = 1.0  1.0 1.0

ref_t               = 315 315 315  ;used in AA simulation


Pcoupl              = parrinello-rahman

Pcoupltype          = semiisotropic

tau_p               = 12.0 ; parrinello-rahman is more stable with larger tau-p, DdJ, 20130422

compressibility     = 3e-4  3e-4  3e-4

ref_p               = 1.0  1.0  1.0


gen_vel             = no

gen_temp            = 315

gen_seed            = 473529


; MARTINI and CONSTRAINTS

; for ring systems and stiff bonds constraints are defined

; which are best handled using Lincs.


constraints              = none

constraint_algorithm     = Lincs


; COM motion removal

; These options remove motion of the protein/bilayer relative to the solvent/ions

nstcomm        = 100

comm-grps = Protein_POPC_POPE_CHOL W_ION