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

[Justin Lemkul]

On 12/4/17 1:33 AM, Sailesh Bataju wrote:
> 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.

Using direction-periodic requires an NVT ensemble, so you would be 
changing this for some windows and not others. That's not a good idea. 
The simplest solution is to use a larger box. Based on the error 
message, your box is already very small so increasing it slightly 
shouldn't be a big deal, but you would have to repeat the simulations.

-Justin

-- 
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Virginia Tech Department of Biochemistry

303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.biochem.vt.edu/people/faculty/JustinLemkul.html

==================================================