Pull simulation error: Umbrella sampling simulation of peptide and lipid system (MARITNI CG)

User discussions
1

GROMACS version:2024
GROMACS modification: No
Hi all,

I am trying to simulate a system of lipids and a peptide (MARTINI CG) with an umbrella sampling method. The pull simulation was performed, and the pull.xvg was plotted. The graph looks like below, and the energy steadily increases even after the peptide is outside the membrane.

image
image689×509 20.6 KB

Please let me know why this issue is arising. I have attached the pull simulation script below.

title = Umbrella pulling simulation
define = -DPOSRES -DPOSRES_FC=50 -DBILAYER_LIPIDHEAD_FC=10
; Run parameters
integrator = md
dt = 0.020
tinit = 0
nsteps = 50000000 ; 1000000 ps or 1us
nstcomm = 10
; Output parameters
nstxout = 50000 ; every 100 ps
nstvout = 50000
nstfout = 50000
nstxtcout = 50000 ; every 100 ps
nstenergy = 50000

; 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.1
rcoulomb = 1.1
rvdw = 1.1

; 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 = PROTEIN MEMBRANE SOLUTE
tau_t = 1.0 1.0 1.0
ref_t = 310 310 310

; Pressure coupling is on
; Pressure coupling:
Pcoupl = berendsen
Pcoupltype = semiisotropic
tau_p = 5.0
compressibility = 4e-5 4e-5
ref_p = 1.0 1.0
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 = 2 ; two groups defining one reaction coordinate
pull_group1_name = PROTEIN
pull_group2_name = MEMBRANE
pull_coord1_type = umbrella ; harmonic potential
pull_coord1_geometry = direction; simple distance increase
pull_coord1_vec = 0 0 1
pull_coord1_dim = N N Y ; pull along z
pull_coord1_groups = 1 2 ; groups 1 and 2 define the reaction coordinate
pull_coord1_start = yes ; define initial COM distance > 0
pull_coord1_rate = 0.000015 ;
pull_coord1_k = 1000 ; kJ mol^-1 nm^-2
pull_coord1_init = 0.2

pull_group2_pbcatom =1571
pull-pbc-ref-prev-step-com = yes

Thanks,
Rakesh

2

It shows that the pull force is increasing over time, i.e. the distance between the pull reference position and the group you are pulling. If you look at the trajectory you will probably see that your membrane COM is still relatively unchanged compared to the protein COM. If you continue the simulation you might get a sudden snap as the protein is pulled away from the membrane. That will probably crash your simulation.

You could see if it helps increasing the force constant, e.g., pull_coord1_k = 25000. But it might also lead to crashes.

Another alternative, to avoid sudden motions, would be to try if it helps pulling using a constraint: pull_coord1_type = constraint. Then you can remove the pull_coord1_k line. But pulling using constraints can also lead to instabilities, especially if you pull too fast, since the groups are forced to keep exactly the reference distance at every step.

3

You can add:

pull-print-com = yes
pull-print-ref-value = yes

and have a look at the pullx.xvg output to see how the actual COM coordinates compare to the reference value.