[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.
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
More information about the gromacs.org_gmx-users