GROMACS version: 2023.1
GROMACS modification: No
I am applying an umbrella biasing potential to a ligand that is at the center of a DOPC membrane (consisting of 72 lipids) in order to pull it out and use the Weighted Histogram Analysis Method (WHAM) to extract the PMF.
The membrane system is surrounded by TIP3P water molecules, and it was made through CHARMM-GUI website using the charmm36 force field.
During the time of my initial NVT equilibration and then at my NPT equilibration i observe the hydrophobic mismatch that causes lipid tails to move apart giving unrealistic results. I thought of applying flat-bottomed position restraints after using this script in order to measure the length of the lipids at z axis and also the length of the ligand. I get DOPC lipids length = 4.905 and ligand length = 0.8 which i divide by 2 and add some length. I create .itp files like that: ; position restraint file for DOPC
[ position_restraints ]
; id funct g r k
; Atom number Function (2 = flatbottomed) Direction (8 = cylindrical) Radius Force constant
1 2 8 -3.1 500
2 2 8 -3.1 500
3 2 8 -3.1 500
4 2 8 -3.1 500
5 2 8 -3.1 500
6 2 8 -3.1 500
7 2 8 -3.1 500
I then add the essential lines at topol.top file and i define the restraints at my .mdp files for equilibration. This is my mdp file for nvt equilibration:
NVT equilibration
define= -DPOSRES_LIG_flat_bottomed -DPOSRES_LIPID_flat_bottomed
; Run control
integrator = md
dt = 0.001 ; Decreased it to avoid instability errors
nsteps = 50000
nstcomm = 1
nstcalcenergy = 1
comm-mode = Linear
comm-grps = DOPC POT CLA TIP3 ATE
; Output control
nstxout = 1000
nstvout = 1000
nstenergy = 1000
nstlog = 1000
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs
constraints = all-bonds
lincs_iter = 2 ; Increased it to run my eq with flat bottomed restraints
lincs_order = 6 ; Increased it to run my eq with flat bottomed restraints
; Neighborsearching
ns_type = grid
nstlist = 5
rlist = 1.2
rcoulomb = 1.2
rvdw = 1.2
; Electrostatics
coulombtype = PME
pme_order = 4
fourierspacing = 0.16
; van der Waals
vdwtype = cutoff
; Temperature coupling
tcoupl = V-rescale
tc-grps = system
tau_t = 0.1 ; reduced to avoid errors with big variance and instability
ref_t = 323
; Pressure coupling
pcoupl = no
; Periodic boundary conditions
pbc = xyz
; Dispersion correction
DispCorr = EnerPres
; Velocity generation
gen_vel = yes
gen_temp = 323
gen_seed = -1
But when i try to run my simulation, i get messages like these:
starting mdrun ‘Title’
50000 steps, 50.0 ps.
WARNING: Listed nonbonded interaction between particles 7840 and 7847
at distance 2.284 which is larger than the table limit 2.268 nm.
This is likely either a 1,4 interaction, or a listed interaction inside
a smaller molecule you are decoupling during a free energy calculation.
Since interactions at distances beyond the table cannot be computed,
they are skipped until they are inside the table limit again. You will
only see this message once, even if it occurs for several interactions.
IMPORTANT: This should not happen in a stable simulation, so there is
probably something wrong with your system. Only change the table-extension
distance in the mdp file if you are really sure that is the reason.
Is there any way to fix it or avoid the hydrophobic mismatch? I have also thought of applying position restraints at the phosphorous atoms or the lipids headgroups.
Thank you for your time!