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

zeineb SI CHAIB zeineb-14 at hotmail.com
Fri May 25 13:30:41 CEST 2018

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


; 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


; 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


; 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


; 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

More information about the gromacs.org_gmx-users mailing list