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
Thanks for the reply, Hess. I have a stacked bilayer, hence the larger box size.
I am calculating the free energy profile of a molecule through the lipid bilayer.
Like you suggested, I did change the pull_coord1_geometry = cylinder ; and pull-cylinder-r = 1.5 ;
However, I am getting the following error, "Fatal error: A cylinder pull group is not supported when using absolute reference!"
Currently I am specifying a single atom of my molecule being pulled using “pull-group1-name”
The simulation is currently running without any errors and completed 800000 steps.
The following is the modified script. (Below, Italics)
I tried changing the pull_coord1_geometry to cylinder, but it also didnt work.
So, I kept it as ‘direction’ and removed the position restraints on all the atoms and kept only the following, define = -DHEAVY_H
I also removed “pull_coord1_origin = 2.5 2.5 10.0” and introduced “pull_pbc_ref_prev_step_com = yes” and “pull_group1_pbcatom = 12148” (nearest atom to MOL at t=0, below MOL)
The “pullx.xvg” file is increasing linearly as a function of steps, and “pullf.xvg” remains almost constant as a function of steps.
Note: Here ‘MOL’ is the molecule (39 atoms) that I am trying to pull along the bilayer.
Requesting your suggestions, if any corrections are required to this.
; Pull code pull = yes pull_ncoords = 1 ; only one reaction coordinate pull_ngroups = 2 ; two groups defining one reaction coordinate pull_group1_name = Bilayer pull_group2_name = MOL pull_coord1_type = umbrella ; harmonic potential pull_coord1_geometry = direction ; simple distance increase pull_coord1_dim = N N Y pull_coord1_groups = 1 2 pull_pbc_ref_prev_step_com = yes pull_group1_pbcatom = 12148 pull_coord1_vec = 0 0 1 pull_coord1_start = yes ; define initial COM distance > 0 pull_coord1_rate = 0.0001 ; 0.002 nm per ps = 0.25 nm per ns pull_coord1_k = 5000 ; kJ mol^-1 nm^-2
I have one more important query/observation.
My forcefields were taken from the literature. When ’ define = -DHEAVY_H’ was selected, the atomic mass of hydrogen was increased (3 times), and that of oxygen and other heavy atoms were reduced (by a number 3). Under such conditions my pull simulations ran properly. I was able to pull my molecule across the bilayer. Can you please tell me the concept of atomic mass reduction in GROMACS. When is it used, and what are its implications?
Thanks again
I wouldn’t think that mass repartitioning should affect this issue directly, as you seemed to be pulling on whole molecules. But if the forces are high it could be that a hydrogen atom becomes unstable. Your force constant of 5000 is not that high though.
How would this repartitioning affect the PMF?
Since the purpose of pull simulation (for PMF calculation) is to sample the molecule across a bilayer, the Hydrogen mass repartitioning (HMR) should not affect my PMF values. Because, from the sampling using pull simulations I further do umbrella simulations, where I am not using any HMR.
Please correct me if I am wrong.
Hi Hess, I completed the simulation, followed by umbrella sampling and pmf calculation. I have some follow up questions regarding the output. Can you please give your insights. The link to the question is pasted above.
And what about the two output files, Hess? “profile.avg” and “pmfintegrated.xvg”. How is the “pmfintegrated.xvg” file obtained?
Does “profile.xvg” gives the free energy profile of a molecule across the bilayer?
profile.xvg is the free energy profile. If your reaction coordinate is the position of a molecule with respect to the bilayer, then the answer to your question is yes. The pmfintegrated.xvg file is undocumented. Is the the PMF obtained by integrating the force, not using WHAM then, I assume.