AMBER Archive (2009)

Subject: RE: [AMBER] NPT md with PME method for a water box

From: Ross Walker (ross_at_rosswalker.co.uk)
Date: Wed Jan 14 2009 - 12:07:13 CST

Hi Jeffrey,

> I am going to perform md simulation of a water box under NPT
> ensemble. First, I built the water box by packmol and then created the
> .prmtop, inpcrd files in xleap. In order to use the pme method for
> electrostatic interaction, the setBox command was used to add the
> periodic boundary conditions to the unit.
>
> My questions are:
>
> 1. minimization:
> input:
> -------
> energy minimization for starting structures
>
> &cntrl
> imin=1, maxcyc=2000, ncyc=2000, ntpr=20,
> ntx=1,irest=0,
> ntb=1, cut=10.0
> /
> &ewald
> eedmeth=5,
> /
> ----------

For a periodic simulation, where PME is in use you should not set eedmeth=5
- this is using a distance dependent dielectric and thus does not include
the erf correction for PME and is likely the cause of your problems. Just
remove the entire &ewald namelist and leave it at the defaults.

> output:
> -------
> NSTEP ENERGY RMS GMAX NAME NUMBER
> 1020 -1.3947E+09 1.2216E+11 8.6032E+12 H1 677
>
> BOND = 30799.0575 ANGLE = 0.0000 DIHED =
> 0.0000
> VDWAALS = 13173.5998 EEL = ************* HBOND =
> 0.0000
> 1-4 VDW = 0.0000 1-4 EEL = 0.0000 RESTRAINT =
> 0.0000
>
> Where the energy jumped from -1.3947E+09 to -5.1545E+12 for only 20
> steps and the ele interaction was so not printed. Are these behaviors
> reasonable?

The *'s mean the value is too large to fit in the output space. Typically
this means you have atoms sitting on top of each other. This is because they
either have the same or similar coordinates or you box dimensions are
specified incorrectly. This normally also means that VDWAALS is very high
but in the case of TIP3P water the H atoms do not have VDW radii and so you
can have two H's sitting on top of each other or a H and an O which would
generate a very large EEL value but not affect the VDWAALS. However, in your
case the eedmeth=5 confuses matters so it could be that eedmeth=5 is causing
it and/or that you also have atoms sitting on top of each other.
 
> 2. equilibration and production runs:
> Which temperature method should I choose? NTT=1 or NTT=3? It is said in
> the amber mannual that "Unless you are sure you know what you are doing,
> please don't use ntt=1!", but in the benchmark (ambermd.org) input files,
> I notice that the temperature method is set to ntt=1. If ntt=3 is
> adopted, what the best values of gamma_ln at heating, equilibrating,
> production stages?

There is lots of debate on this. You probably want to do something along the
lines of:

1) Minimize with ntb=1

2) Heat slowly over 20ps or so with ntb=1 and ntt=3

3) Equilibrate with ntb=2 and ntt=3 - this will allow the pressure to
stabilize.

4) production with ntb=1 and ntt=0 - this would be NVE. A lot of people also
do ntb=1 and ntt=1 but then set the thermostat coupling constant fairly long
(tautp=10.0 or so) such that it approximates NVE without the issues related
to slow energy loss ove rmultiple ns if various tolerances aren't tightened
up enough for NVE.

There are many many combinations, however, and Ph.D's can be (and have been)
written on this subject.

Note you may find that for initial stability you need to set your box size
slightly larger than the actual box and then allow AMBER to correct it when
you equilibrate with NPT (ntb=2).

Good luck,
Ross

/\
\/
|\oss Walker

| Assistant Research Professor |
| 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.

_______________________________________________
AMBER mailing list
AMBER_at_ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber