GROMACS version: 2023
GROMACS modification: No
I am trying to calculate the PMF of a molecule through a lipid bilayer.
I have an equilibrated bilayer.
I inserted my molecule of interest (single molecule) using ‘insert-molecule’ commands. Then energy minimization was performed (with the molecule in the bilayer), followed by growing the molecule.
Then I restrained the positions of the molecules in the bilayer, and only single atom of the molecule. Ran npt.
After that I again ran npt simulation (without position restraint) for 30ns.
The issue comes after this. When I take this structure (after 30ns NPT) and perform pull simulation (mdp parameters listed below), I am getting LINCS errors.
I even changed the configuration of the molecule in the bilayer and repeated. I am getting same LINCS error.
Please let me know what can be done to avoid LINCS error (error message also pasted below)
Pull MDP:
title = Umbrella pulling simulation
define = -DPOSRES_CHL1 -DPOSRES_CP20 -DPOSRES_CP22 -DPOSRES_CP24 -DPOSRES_CERNS -DPOSRES_CP26 -DPOSRES_CP28 -DPOSRES_CP30 -DPOSRES_CEOS -DPOSRES_FAH20 -DPOSRES_FAH22 -DPOSRES_FAH24 -DPOSRES_FAH26 -DPOSRES_FAH28 -DPOSRES_FAH30
; Run parameters
integrator = md
dt = 0.002
tinit = 0
nsteps = 22000000 ; 44000 ps
; mode for center of mass motion removal
comm-mode = Linear
nstcomm = 1
; Output parameters
nstxout = 0
nstvout = 0
nstfout = 0
nstxtcout = 500 ; every 1 ps
nstenergy = 500
nstlog = 500
; Bond parameters
constraint_algorithm = lincs
constraints = all-bonds
continuation = yes ; continuing from NPT
; Single-range cutoff scheme
cutoff-scheme = Verlet
nstlist = 20
ns_type = grid
rlist = 1.2
rcoulomb = 1.2
rvdw = 1.2
; PME electrostatics parameters
coulombtype = PME
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
ewald_rtol = 1e-5
optimize_fft = yes
; Berendsen temperature coupling is on in two groups
Tcoupl = V-rescale
tc_grps = Bilayer PA
tau_t = 1.0 1.0
ref_t = 303.15 303.15
; Pressure coupling is on
pcoupl = Berendsen ; Pressure coupling on in NPT
pcoupltype = semiisotropic ; uniform scaling of x-y box vectors, independent z
tau_p = 1.0 ; time constant, in ps
ref_p = 1.0 1.0 ; reference pressure, x-y, z (in bar)
compressibility = 4.5e-5 4.5e-5 ; isothermal compressibility, bar^-1
refcoord_scaling = com
; Generate velocities is off
gen_vel = no
; Periodic boundary conditions are on in all directions
pbc = xyz
; Long-range dispersion correction
DispCorr = EnerPres
; Pull code
pull = yes
pull_ncoords = 1 ; only one reaction coordinate
pull_ngroups = 1 ; two groups defining one reaction coordinate
pull_group1_name = PA_C1
pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = direction ; simple distance increase
pull_coord1_dim = N N Y
pull_coord1_origin = 2.5 2.5 10.42
pull_coord1_groups = 0 1
pull_coord1_vec = 0 0 1
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.00025 ; 0.002 nm per ps = 0.25 nm per ns
pull_coord1_k = 5000 ; kJ mol^-1 nm^-2
Error:
Step Time
0 0.00000
Step 0 Warning: pressure scaling more than 1%, mu: 0.621468 0.621468 0.341818
Energies (kJ/mol)
Angle U-B Proper Dih. Improper Dih. LJ-14
2.00986e+01 1.88588e+05 1.19815e+05 7.95625e+02 1.56401e+04
Coulomb-14 LJ (SR) Disper. corr. Coulomb (SR) Coul. recip.
-3.23907e+04 -1.48157e+05 -1.13909e+04 -1.27353e+05 2.49864e+03
Position Rest. COM Pull En. Potential Kinetic En. Total Energy
2.83918e-07 2.95902e-13 8.06531e+03 4.15371e+05 4.23437e+05
Conserved En. Temperature Pres. DC (bar) Pressure (bar) Constr. rmsd
1.07388e+07 6.29547e+02 -2.75762e+02 -1.57249e+05 2.65514e-04
Program: gmx mdrun, version 2023.3-spack
Source file: src/gromacs/mdlib/constr.cpp (line 240)
Fatal error:
Too many LINCS warnings (74491)
If you know what you are doing you can adjust the lincs warning threshold in
your mdp file
or set the environment variable GMX_MAXCONSTRWARN to -1,
but normally it is better to fix the problem