[gmx-users] Umbrella sampling simulation stopped in middle.

[Sailesh Bataju]

Hi,

I've started to pull a dimer of isobutane molecules as a production md
after npt equilibriation. Here is the md_pull.mdp file for production md.

; Umbrella pulling, production md
; preprocessing parameters
title  = Umbrella Sampling
; run parameters
integrator = md ; leap-frog integrator for integrating Newton's equations
of motion
nsteps = 1000000 ; 10^6 steps to integrate or maximize i.e. total =
5*10^5*1x10^-15 = 0.5 ns
dt = 0.001 ; time step for integration i.e. 1fs=0.001ps
nstcomm = 100 ; frequency for center of mass motion removal
; output control parameters
nstxout = 1000 ; save coordinates every 1000*1x10^-15 = 1ps
nstvout = 1000 ; save velocities every 1ps
nstenergy = 1000 ; save energies every 1ps
nstfout = 1000 ; number of steps that elapse between writing forces to
output trajectory
nstlog = 1000 ; update log file every 1ps
nstxtcout = 1000 ; gives 1000 frames for umbrella sampling generates .xtc
file
energygrps = Alkane SOL ; groups of the molecules present in the system
that writes to the energy file
; neighbour searching parameters
cutoff-scheme = verlet
nstlist = 10 ; frequency to update neighbour list
ns_type = grid ; search neighboring grid cells
rlist = 1.0 ; Cut-off distance for the short-range neighbor list
pbc = xyz ; periodic boundary conditions in all directions
; electrostatics, vdw and ewald parameters
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme-order = 4 ; Interpolation order for PME. 4 equals cubic interpolation
fourierspacing = 0.12 ; grid spacing for FFT
rcoulomb = 1.0 ; distance for coulomb cut-off
epsilon-r = 1.0 ; relative dielectric constant
rvdw = 1.0 ; distance for the LJ or Buckingham cut-off (short-range van der
waal's cut-off)
optimize_fft = yes
; temperature coupling in on
tcoupl = Nose-Hoover ; modified Berendsen thermostat
tc-grps = Alkane SOL ; groups to couple separately to temperature bath
tau-t = 0.2 0.2 ; time constant for coupling
ref-t = 298 298 ; according to experiment performed before
ld-seed = 1225 ; random seed for temperature
; pressure coupling is on
pcoupl = Parrinello-Rahman ; The box is scaled every timestep,Exponential
relaxation pressure coupling with time constant tau-p
pcoupltype = isotropic
tau-p = 1 ; time constant for coupling in ps
compressibility = 4.6e-5 ; For water at 1 atm at 300K is 4.6e-6bar^(-1)
ref-p = 1 ; reference pressure
refcoord-scaling = com
DispCorr = EnerPres ; apply long range dispersion corrections for Energy
and Pressure
; generate velocity at 298K
gen-vel = no ; generate velocity initially
gen-temp = 298 ; generate temperature for Maxwell distrAlkanetion
gen-seed = 12345 ; given seed
; bonds parameters
continuation = yes ; continue after npt
constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
constraint-algorithm = LINCS ; LINear Constraint Solver, holonomic
constraints
; pull code
pull  = yes
pull_coord1_type = umbrella ; Center of mass pulling using an umbrella
potential between the reference group and one or more groups.
pull-ngroups = 2  ; there are two groups that are subject to a biasing
potential
pull-group1-name = IBU_IBU_1 ; The name of the pull group, is looked up in
the index file
pull-group2-name = IBU_IBU_2
pull-coord1-geometry = distance ; pull along the vector connecting the two
groups
pull-ncoords = 1 ; The number of pull coordinates
pull-coord1-groups = 1 2 ; groups 1 and 2 (designated by name above as CA
and CB, respectively) define the reaction coordinate
pull-coord1-k = 840 ; The force constant. For umbrella pulling this is the
harmonic force constant in [kJ mol-1 nm-2]
; 2kcal mol^-1 angostrom^-2 = 840 kJ/mol/nm^2 according to paper
pull-coord1-rate = 0.001 ; 0.001 nm per ps = 1 nm per ns
pull-coord1-dim = Y Y Y ; pulling in all x,y and z direction
pull-coord1-start = yes ; use whatever the distances in the coordinate file
are.

It generated 1000 configuration files without any error. With a separation
of windows of 0.05 nm, i've generated about 23 frames for automated
simulation using python script. All frames went good but last 4 frames
produce error like this:-

Program gmx mdrun, VERSION 5.1.4
Source code file:
/home/akash/Downloads/gromacs-5.1.4/src/gromacs/pulling/pull.cpp,
line: 520

Fatal error:
Distance between pull groups 1 and 2 (1.683443 nm) is larger than 0.49
times the box size (1.717592).
You might want to consider using "pull-geometry = direction-periodic"
instead.

For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors

I think this error has to warn us during production md so that future
simulation runs well. Should I change pull-geometry = direction-periodic
(from distance) to continue simulation? If I do so for the remaining 4
frames, does it produce different results (relative to pull-geometry =
distance)? I've no idea what is the best to be done for the last 4 frames.
Considering the time I spent for the simulation, please introduce me the
best alternative way to go ahead.

Thank You
Sailesh Bataju

-- 
Self-reliant is the great potential for success.