AMBER Archive (2006)

Subject: RE: AMBER: Problems with equilibrating the starting density

From: Ross Walker (ross_at_rosswalker.co.uk)
Date: Fri Nov 03 2006 - 13:41:55 CST

Dear Markus,

There is so much wrong with your attempted runs that you should probably
consider working through the tutorials on http://amber.scripps.edu/ first of
all to make sure you have an understanding of how Amber works and MD in
general.

I will attempt to outline some of the problems but without understanding
exactly what you are doing you can end up wasting a significant amount of
time running simulations that are not going to be meaningful.

1) Unless you have a good understanding of how to script things in Unix I
would avoid doing it. Setup the input files yourself as standalone files and
manually run the command lines yourself. Do not try to hack somebody elses
script that you don't necessarily understand. You will only end up confusing
yourself more and wasting time. Manually creating things yourself will also
give you the necessary time to think about the options you are selecting and
whether they are appropriate. I suggest that you try the tutorials to see
how you can run Amber jobs and then you can create your own scripts where
you understand how they work.

2) You don't mention if you did any minimisation or not before starting your
MD simulation. However, if you didn't and you are taking a leap generated
solvated structure you will almost certainly have bad contacts that need to
be minimised to stop your system blowing up when you start MD. See the
tutorials for more details.

3) There are numerous 'issues' with your input files. E.g.

> Titel
> &cntrl
> imin=0, irest=0, ntx=1,
> ntt=1, temp0=300.0, tautp=1.0, taup=0.5,

ntt=3 is the recommended method these days for temperature equilibration.
Although ntt=1 will likely not cause you problems.

taup=0.5 is dangerously short and will likely lead to instabilities in the
simulation. I believe the recommended 'lower' bound on this number is 1.0
and sometimes larger values up to 5.0 may be needed for stability.

> ntb=2, ntp=2,

When starting from 0K you should never use constant pressure. This is
because the calculation of pressure at low temperatures is very inaccurate
and can lead to the barostat doing massive overcorrections which lead to
instability. Instead you should typically heat to around 300K or so with
constant volume and only then switch to constant pressure. Ideally (although
other's opinions may differ here) you probably want to run heating at
constant volume, switch to constant pressure until the pressure and density
have equilibrated. Then switch back to constant volume for the remainder of
the equilibration and then when you are looking at doing production MD turn
off the thermostat all together and run with constant volume.

ntp=2 -> This is anisotropic position scaling and is really only designed to
be used with non isotropic systems such as membranes. For proteins in
solution you should use ntp=1.

> cut=20,

This is a huge cut off for a periodic boundary simulation using PME. If you
are using periodic boundaries with PME (the default) then you only need a
cutoff of around 8 angstroms. More than this and you pretty much just waste
computer time. Such a large cut off can also cause you problems with small
boxes as atoms in the central box may see their own images in the direct
space sum. This will cause sander to quit with an error message.

Short answer, if you are doing gas phase simulations set this to 99999.0
(i.e. no cut off). If you are doing implicit solvent GB set this to 9999.0
if you can afford it or else something larger than about 16.0. If you are
doing implicit solvent set it to 8.0 unless for some reason you want to turn
off PME (I can't see why you'd want to).

> ntc=1, ntf=1,

If you are using tipXp water you should be using shake since these water
models were parameterised for use with shake. If you are using flexible
water then the above is correct. I assume this should probably in your case
be ntc=2 and ntf=2. You also do not specify a time step (dt). Fortunately
the default is 0.001 which is okay for no shake but you should probably turn
on shake and specify dt=0.002 (2fs).

> iwrap=0,
> nstlim=35000,
> ntwe=35001, ntwx=100, ntpr=100,

Specifying ntwe=35001 is kind of redundant. If you don't want the energy
file set it to 0 rather than 'greater than nstlim'. In this way if you
change nstlim you don't also need to remember to change ntwe.

> ntwr=100,

This is extremely short for writing a restart file and will likely adversely
effect performance, especially in parallel. A value of around 10,000 or so
is probably better unless you expect your job to crash frequently due to
issues with your machines.

> ntwprt=0,
> &end
> &wt
> type='END'
> &end

Unless you have weight restraints you probably shouldn't specify an &wt
namelist. While there is no problem with having it even though it is not
needed it does increase the chances for you to make typos in the input file
that cause sander to quit when it is started.

> /d2/amber8/amber8/exe/sander \
> -O -i ddrug011.in -o ddrug011.out -p drug011.top \
> -c drug011.rst -x ddrug011.trj -inf ddrug011.inf \
> -r ddrug011.rst -ref drug011.rst

Unless you have restraints turned on (ntr=1) you don't need to specify a
reference file (-ref).

> cat > dkdrug011.in << schluss
> Titel
> &cntrl
> imin=0, irest=0, ntx=1,
> ntt=1, tempi=300.0, temp0=300.0,

For the second step of MD, unless for some reason you want to throw away the
velocities from the previous MD stage and instead reset them to a random
gaussian distribution at 300K then you should set irest=1 and ntx=5.
Similarly if irest=1 you don't need to specify tempi since it is just
ignored and the initial temperature will be whatever the final temperature
was at the end of the previous run since the coordinates and velocities are
the same.

> ntb=0, ntp=0,

In you second stage you suddenly turn off periodic boundaries - why???

> /d2/amber8/amber8/exe/sander \
> -O -i dkdrug011.in -o dkdrug011.out -p drug011.top \
> -c ddrug011.rst -x dkdrug011.trj -inf dkdrug011.inf \
> -r dkdrug011.rst -ref ddrug011.rst

again no ref is needed as you don't have restraints.

> cat > kdrug011.in << schluss
> Titel
> &cntrl
> imin=0, irest=1, ntx=5,
> ntt=0, tempi=300.0,
> ntb=0, ntp=0,

Again this is a gas phase simulation - a block of water flying through a
vacuum. The water will very quickly evaporate since 300K is above the
boiling point of water in a vacuum. Hence the results from this simulation
are likely to be meaningless.

Onto the actual error you see:

> File name = fort.9 formatted, sequential access record = 1
> In source file _ew_box.f, at line number 1745
> PGFIO-F-217/list-directed read/unit=9/attempt to read past
> end of file.
> File name = fort.9 formatted, sequential access record = 1
> In source file _ew_box.f, at line number 1745

This is from the code trying to read past the end of the inpcrd file. I
suspect that your inpcrd file is not actually for a solvated system or it
doesn't have any box information in it. The reason it doesn't work for the
first simulation is that this one is a periodic boundary simulation and so
it expects box info. Whereas the others are gas phase simulations so don't
look for the box information that isn't there anyway.

Thus you have a large number of problems that you need to address. Do you
really want to be running gas phase simulations? (REALLY, REALLY???). Or do
you want to run explicit solvent PBC simulations with PME? In this case you
need to create yourself In either case you need to fix your inpcrd file to
have box information. I suggest you recreate it in leap in line with how it
is done in the tutorials. You also need to fix all the issues with your
input files or you will be very lucky if you get reasonable simulation
results out.

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