GROMACS version:2025.1
i run this umbrella sampling simulaiotn on a small molecule from a water layer to an octanol water layer to determine partionon coeficient. i am able to produce graphs that look decent but when i visualize the molecules it looks as though it is being pulled through the periodic box to the next layer instead of being pulled from the water to the octanol/water layer. is there anyway i can modify this mdp file below so it doesn’t do that? the tutorial im following cna be found here: Free energy techniques – Martini Force Field Initiative . this mdp file is not the original files that come with the downloaded files in the tutorial i was following.
below is what i am using for my pull.mdp file
; RUN CONTROL PARAMETERS
integrator = sd
dt = 0.01
nsteps = 10000
nstcomm = 100
comm-grps =
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 1000
nstenergy = 100
nstxout-compressed = 1000
compressed-x-precision = 100
compressed-x-grps =
;energygrps = B W OCO
; NEIGHBOURLIST and MARTINI
; To achieve faster simulations in combination with the Verlet-neighborlist
; scheme, Martini can be simulated with a straight cutoff. In order to
; do so, the cutoff distance is reduced 1.1 nm.
; Neighborlist length should be optimized depending on your hardware setup:
; updating ever 20 steps should be fine for classic systems, while updating
; every 30-40 steps might be better for GPU based systems.
; The Verlet neighborlist scheme will automatically choose a proper neighborlist
; length, based on a energy drift tolerance.
;
; Coulomb interactions can alternatively be treated using a reaction-field,
; giving slightly better properties.
; Please realize that electrostVatic 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.
;
; With the polarizable W model, the relative electrostatic screening
; (epsilon_r) should have a value of 2.5, representative of a low-dielectric
; apolar solvent. The polarizable W itself will perform the explicit screening
; in aqueous environment. In this case, the use of PME is more realistic.
;cutoff-scheme = group
cutoff-scheme = Verlet
nstlist = 20
;nstlist = 5
ns_type = grid
pbc = xyz
;verlet-buffer-tolerance = 0.005
;verlet-buffer-tolerance = 2e-03
coulombtype = reaction-field
coulomb-modifier = Potential-shift-verlet
rcoulomb = 1.1
epsilon_r = 15 ; 2.5 (with polarizable W)
vdw_type = cutoff
vdw-modifier = Potential-shift-verlet
rvdw = 1.1
; MARTINI and TEMPERATURE/PRESSURE
; normal temperature and pressure coupling schemes can be used.
; It is recommended to couple individual groups in your system separately.
; Good temperature control can be achieved with the velocity rescale (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
Tcoupl = berendsen
tc-grps = System
;tau_t = 1.0 1.0
tau_t = 2.0
ref_t = 300
;nsttcouple = 1
;Pcoupl = parrinello-rahman
Pcoupl = berendsen
Pcoupltype = semiisotropic
tau_p = 12.0 ;parrinello-rahman is more stable with larger tau-p, DdJ, 20130422
compressibility = 0 3e-4
ref_p = 0 1.0
;nstpcouple = 1
gen_vel = no
gen_temp = 320
gen_seed = 473529
; MARTINI and CONSTRAINTS
; for ring systems and stiff bonds constraints are defined
; which are best handled using Lincs.
constraints = hbonds
constraint_algorithm = Lincs; COM PULLING
; Pull type: no, umbrella, constraint or constant_force
pull = yes
; Pull geometry: distance, direction, cylinder or position
pull-coord1-type = umbrella
pull-coord1-geometry = distance
; Select components for the pull vector. default: Y Y Y
pull-coord1-dim = N N Y
; Switch from r1 to r0 in case of dynamic reaction force
pull_constr_tol = 1e-06
pull-coord1-start = no
pull_nstxout = 100
pull_nstfout = 100
; Number of pull groups
pull_ngroups = 2
; Group name, weight (default all 1), vector, init, rate (nm/ps), kJ/(mol*nm^2)
pull-group1-name = W
pull-group2-name = TT1
pull-coord1-groups = 1 2
;pull-coord1-vec = 0 0 1
pull-coord1-init = POS1
pull-coord1-rate = 0
pull-coord1-k = 1000
pull-coord1-kB = 1000
pull-pbc-ref-prev-step-com = yes
pull-group1-pbcatom = 1805
pull-group2-pbcatom = 0
; Free energy control stuff
free-energy = no ;yes
;init-lambda = 0
;delta-lambda = 0
;foreign_lambda =
;sc-alpha = 1.4
;sc-power = 1
;sc-sigma = 0.45
;nstdhdl = 100
;separate-dhdl-file = yes
;dhdl-derivatives = yes
;dh_hist_size = 0
;dh_hist_spacing = 0.1
;couple-moltype = TT1
;couple-lambda0 = vdw-q
;couple-lambda1 = none
;couple-intramol = no
