[gmx-users] Simulated Annealing Protocol...

[Justin A. Lemkul]

rama david wrote:
> Dear Gromacs Specialists,
>                 I am very novice to Molecular Simulation study.
>  I am using GROMACS 4.5.4 version .
> I completed some GROMACS tutorials , I not found any
> tutorial on Simulated Annealing..
> If Any one know the link please give me it..
> 
> I make my  protocol to work on simulated annealing as follow ...
> (I am not writing in detail sorry for that )
> 1.  pdb2gmx ...
> 2. editconf
> 3. Solvent Addition
> 4.  Ion addition
> 5. Energy minimisation
> 6. simulated annealing
>  mdp for simulated annealing is as follow...
> 
> ; 7.3.3 Run Control
> title                   = simulated run
> 
> integrator              = md                    ; md integrator
> tinit                      = 0                     ; [ps] starting time for run
> dt                         = 0.002                 ; [ps] time step
> for integration
> nsteps                  = 5000000               ; maximum number of
> steps to integrate, 0.002 * 5000000 10 ns
> comm_mode         = Linear                ; remove center of mass translation
> comm_grps           = Protein Non-Protein   ; group(s) for center of
> mass motion removal
> 
> ; 7.3.8 Output Control
> nstxout                 = 100       ; [steps] freq to write
> coordinates to trajectory
> nstvout                 = 100       ; [steps] freq to write velocities
> to trajectory
> nstfout                 = 100       ; [steps] freq to write forces to trajectory
> nstlog                  = 100           ; [steps] freq to write
> energies to log file
> nstenergy               = 100           ; [steps] freq to write
> energies to energy file
> nstxtcout               = 500           ; [steps] freq to write
> coordinates to xtc trajectory
> xtc_precision           = 500          ; [real] precision to write xtc
> trajectory
> xtc_grps                = System        ; group(s) to write to xtc trajectory
> energygrps              = System        ; group(s) to write to energy file
> 
> ; 7.3.9 Neighbor Searching
> nstlist                 = 5             ; [steps] freq to update neighbor list
> ns_type                 = grid          ; method of updating neighbor list
> pbc                     = xyz           ; periodic boundary conditions
> in all directions
> rlist                   = 1.0           ; [nm] cut-off distance for
> the short-range neighbor list
> 
> ; 7.3.10 Electrostatics
> coulombtype             = PME           ; Particle-Mesh Ewald electrostatics
> rcoulomb                = 1.0           ; [nm] distance for Coulomb cut-off
> 
> ; 7.3.11 VdW
> vdwtype                 = cut-off       ; twin-range cut-off with
> rlist where rvdw >= rlist
> rvdw                    = 1.4           ; [nm] distance for LJ cut-off
> DispCorr                = EnerPres      ; apply long range dispersion
> corrections for energy
> 
> ; 7.3.13 Ewald
> fourierspacing          = 0.16          ; [nm] grid spacing for FFT
> grid when using PME
> pme_order               = 4             ; interpolation order for PME, 4 = cubic
> ewald_rtol              = 1e-5          ; relative strength of
> Ewald-shifted potential at rcoulomb
> 
> ; 7.3.14 Temperature Coupling
> tcoupl                  = berendson                   ; temperature
> coupling with Nose-Hoover ensemble
> tc_grps                 = Protein    Non-Protein        ; groups to
> couple seperately to temperature bath
> tau_t                   = 0.1        0.1                ; [ps] time
> constant for coupling
> ref_t                   = 5           5               ; [K] reference
> temperature for coupling
> annealing               = single single
> annealing_npoint        = 2  2
> annealing_time          = 0 20 0 20
> annealing_temp          = 5 333 5 333
> ; 7.3.15 Pressure Coupling
> pcoupl                  =  parrinello-Rahman    ; pressure coupling
> where box vectors are variable
> pcoupltype              = isotropic             ; pressure coupling in
> x-y-z directions
> tau_p                   = 2.0                   ; [ps] time constant
> for coupling
> compressibility         = 4.5e-5                ; [bar^-1] compressibility
> ref_p                   = 1.0                   ; [bar] reference
> pressure for coupling
> 
> ; 7.3.17 Velocity Generation
> gen_vel                 = yes           ; velocity generation turned off
> gen_temp                = 5
> ; 7.3.18 Bonds
> constraints                  = all-bonds     ; convert all bonds to constraints
> constraint_algorithm    = LINCS         ; LINear Constraint Solver
> continuation               = yes           ; apply constraints to the
> start configuration
> lincs_order                 = 4             ; highest order in the
> expansion of the contraint coupling matrix
> lincs_iter                    = 1             ; number of iterations
> to correct for rotational lengthening
> lincs_warnangle         = 30            ; [degrees] maximum angle that
> a bond can rotate before LINCS will complain
> 
> My Queries
> 
> 1. Is My mdp file ok ???..please give me a nice protocol..

It is a bad idea to generate velocities and use P-R pressure coupling at the 
same time.  Such an approach is extremely fragile and your simulation is likely 
to crash.

I assume that since you have specified heating over 20 ps of the 10 ns 
trajectory that you wish to have the temperature maintained from then on.  I 
believe Gromacs will handle this elegantly.  Pay attention to everything grompp 
tells you, in any case.

> 2.. Should I have to do position restrained MD before SA(simulated annealing)
>     (If yes then what the temp. should I have to used for NVT and NPT
> (as in mdp file has lower 5 K and high 300 k  ) )

The answer to this question depends on what you're trying to accomplish with SA. 
  Position restraints are used to disfavor major structural changes in the 
solute as the solvent is equilibrated around it.  If SA is your initial step in 
equilibration, I would think position restraints can't hurt.

My general protocol is NVT at the initial temperature, SA to heat to the desired 
temperature, NPT at the final conditions, and then data collection.  You should 
base your decision upon what you are wishing to observe and by interpreting how 
others have done similar things.

-Justin

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

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

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