AMBER Archive (2006)

Subject: RE: AMBER: MD and MMPBSA

From: Ross Walker (ross_at_rosswalker.co.uk)
Date: Wed Mar 08 2006 - 23:26:34 CST

Hi Alexandra,

To answer some of your questions:

> 1-I think leap is not numbering the water residues correctly.
> The complex has
> 1024 residues and leap added 13992 water residues. In the pdb file
> created from
> the crd and top files, after the residue 9999, the residues numbering
> restart in
> 0. This leads to some residues having repeated numbers,
> although the atom
> numbers are correct. Ambmask recognize only the residues

This is a limitation of the pdb format which assigns only 4 characters for
the residue number which limits the maximum residue number to 9999. Pdb's
were not designed to have the structure of a protein surrounded by explicit
water. This 'should' only be a cosmetic issue in the pdb file, the prmtop
file will have full numbering and so you shouldn't encounter any problems
here. Although I have never used mm_pbsa on a system with more than 9999
total residues so can't say for sure.

If your calculations are all based on the prmtop file then you will be fine
since this does not have the 9999 residue limitation.

> 2- I have some doubts about the restraints used in minimizations and
> MD. I want
> to perform the minimization (explicit solvent) in three
> steps: first to
> hydrogens and solvent, second to lateral chains and the last to all.
> What force
> constant value one´s shall use? I am confused with that
> because some tutorials
> use 500, while others 5000. Is it better to come back to restraint the
> hydrogens together with solvent in the MD equilibration phase? I am
> asking that
> because I saw a tutorial, available in the internet, in which
> that was done.

Okay, I should really update the DNA tutorial to make this clearer.
Typically you do not want to use a restraint as large as 500 (5000 is way
too much) since such a large restraint gives you a big force constant for
the restraint and thus high frequency oscillations during MD. With a
restraint force constant of 500 these oscillations may be of such a high
frequency that even with a 1fs time step the simulation is not stable. Thus
for MD simulations a restraint value of 5 to 50 is more appropriate. For
minimisation you do not have to worry about the issue of a time step and so
large restraint values can be used without causing instabilities. A value of
500 will keep the restrained system rigidly fixed. Since the real aim of the
restraints in this specific situation is to allow the system to relax away
from the crystal structure in a 'controlled' fashion so a small restraint
force (say 10) would probably suffice.

> 3- I also have a doubt about the initial temperature (temp0)

Initial temperature is tempi - target temperature is temp0. I know this
seems kind of backward but hey, I gave up wondering why such things are the
case long ago... ;-)

> value that one´s
> shall use in the equilibration phase. I saw some tutorials in which
> temp0=0 but
> in others, temp0=300.0 What shall be used?

This really depends on how good / stable your starting structure is. For a
small system in gas phase or implicit solvent you can probably 'get away'
with starting it at 300K. However, for most protein simulations and explicit
solvent simulations the initial structures will not be as stable. This is
largely due to the addition of hydrogen atoms at standard bond lengths,
crystal packing frustration that is present in the crystal structure but not
indicitive of the system in solution and steric clashes due to the limited
resolution of crystal structures. In this situation 'hot starting' the
system can lead to some atoms having large velocities. This gives a very
sudden deviation from the equilibrium geometry and can lead to sudden
unfolding of the protein or the system 'blowing up'. Hence for large systems
in explicit solvent it is advisable to set the initial temperature to 0 K
and then slowly heat the system using NMR weight restraints to control the
temperature if required. In this way the system can relax away from a high
energy starting configuration in a controlled fashion. This approach does of
course require more initial startup time as you have to heat the system
slowly - so if you are brave set temp0=300.0 and don't worry about it. If
you are more cautious set it to a smaller value and increment it slowly.

> My md.in file for equilibration phase is:
>
> &cntrl
> imin = 0, irest=0, ntx=1,
> nstlim=1000, dt=0.002, scee = 1.2, nrespa=2,
> ntt=3, gamma_ln=1.0, tempi=0.0, temp0=300.0,
> ntb = 1, ntp=0,
> ntc=2, ntf=2,
> ntr = 1,
> ntwe=100, ntwx=100, ntpr=200,
> restraint_wt = 50.0
> restraintmask='!:WA= & !@H=',
> cut = 10.0
> /

This looks good to me, although you are normally okay with a cut of 8.0
angstroms with periodic boundary pme simulations (ntb=1) and you may want to
set the value of ntpr to be the same as ntwx - this can just make later
analysis slightly easier. Note 2ps is pretty short for an equilibration, you
probably want to make this quite a bit longer. Plot the potential energy as
a function of time and you should be able to see if the system has
equilibrated.

Note, also you can use NMR weight restraints to obtain a more controlled
heating profile if necessary:

  nmropt=1,
 /
 &wt type='TEMP0', istep1=0,istep2=25000,
                           value1=0.0, value2=300.0, /
 &wt type='END' /

This will linearly ramp the value of temp0 from 0.0 to 300.0k over 25000
steps.

I hope this helps.
All the best
Ross

/\
\/
|\oss Walker

| HPC Consultant and Staff Scientist |
| San Diego Supercomputer Center |
| Tel: +1 858 822 0854 | EMail:- ross_at_rosswalker.co.uk |
| http://www.rosswalker.co.uk | PGP Key available on request |

Note: Electronic Mail is not secure, has no guarantee of delivery, may not
be read every day, and should not be used for urgent or sensitive issues.

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu