[gmx-users] Equilibrating Water around the Molecule

[Justin A. Lemkul]

On 5/24/12 1:58 PM, Lara Bunte wrote:
> Hello
>
> I want to equilibrate the water around my molecule. For doing this I created a pr.mdp file that makes most sense for me. Sadly I am a beginner with less knowledge, so could you please take a look on this pr.mdp file and say me if this is good or if not, what would you suggest to me? My System is Isoalloxazin in a water box. I use a CHARMM27 force field and a tip3p water model. Here is the file:
>
> title           = isoalloxazin NPT equilibration
> define      = -DPOSRES  ; position restrain the protein
>
> ; Run parameters
> integrator      = md            ; leap-frog integrator
> nsteps          = 50000         ; 2 * 50000 = 100 ps
> dt              = 0.002         ; 2 fs
>
> ; Output control
> nstxout         = 5000          ; save coordinates every 4 ps
> nstvout         = 5000          ; save velocities every 4 ps
> nstenergy       = 5000          ; save energies every 4 ps
> nstlog          = 5000          ; update log file every 4 ps
>
> ; Bond parameters
> continuation    = yes           ; Restarting after NVT
> constraint_algorithm = lincs    ; holonomic constraints
> constraints     = all-bonds     ; all bonds (even heavy atom-H bonds) constrained
> lincs_iter      = 1             ; accuracy of LINCS
> lincs_order     = 4             ; also related to accuracy
>
> ; Neighborsearching
> ns_type         = grid          ; search neighboring grid cells
> nstlist         = 1             ; edited to 1
> rlist           = 1.0           ; short-range neighborlist cutoff (in nm)
> rcoulomb        = 1.0           ; short-range electrostatic cutoff (in nm)
> rvdw            = 1.0           ; short-range van der Waals cutoff (in nm)
>

These cutoff settings are incorrect for CHARMM force fields, and setting 
nstlist=1 is only ever necessary for EM.  For running real MD, it is both 
inefficient (and will slow down your calculation) and unnecessary.

> ; Electrostatics
> coulombtype     = PME           ; Particle Mesh Ewald for long-range electrostatics
> pme_order       = 4             ; cubic interpolation
> fourierspacing  = 0.12          ; grid spacing for FFT
>
> ; Temperature coupling is on
> tcoupl           = nose-hoover

Nose-Hoover is a poor choice for initial equilibration.

> tc-grps         = Isoalloxazin Sol

I suspect "System" would be better here.  A small molecule like isoalloxazine 
does not justify its own thermostat.

> tau_t           = 0.1   0.1     ; time constant, in ps
> ref_t           = 300   300     ; reference temperature, one for each group, in K
>
> ; Pressure coupling is on
> pcoupl          = Parrinello-Rahman     ; Pressure coupling on in NPT

Parrinello-Rahman is not a stable algorithm for initial equilibration.

> pcoupltype      = isotropic     ; uniform scaling of box vectors
> tau_p           = 2.0           ; time constant, in ps
> ref_p           = 1.0           ; reference pressure, in bar
> compressibility = 4.5e-5        ; isothermal compressibility of water, bar^-1
>
> ; Velocity generation
> gen_vel     = no        ; assign velocities from Maxwell distribution

If this is the first MD run, gen_vel needs to be set to "yes."

-Justin

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

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

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