From: "Yong Duan"
Subject: RE: failure of minimization
Date: Fri, 21 Feb 2003 21:34:04 -0500


Dear Ioana:

If the first step energy is not a number, it could mean there are a few
particles were very close to each other in your initial model. In
principle, the minimizer should be able to handle it. In reality, it
depends on the choice of parameters. For example, one may wish to use
very conservative step size limits in these cases. The default dxm=0.5
is probably too large for your case, even though it is quite
approperiate for typical systems. I personally would use dxm=0.1 to
start the minimization.

Secondly, when you set NTC=1, the solvent model (water) is no-longer the
default TIP3P model since you allow the bond lengths to change. To
understand why this could be a problem, one needs to know that the van
der Waals force of the hydrogen atoms in TIP3P model is zero, i.e., the
hydrogen atoms does not feel the van der Waals force which is mainly
responsible to separate atoms. This is not a bug or anything like that.
It is part of the TIP3P model. The things that prevent hydrogen atoms
from getting on top of other atoms is the bonding force between H and O
and the van der Waals force of the oxygen atom. In typical simulations,
the bond lengths of TIP3P model is fixed. In your case, when the bond
lengths are allowed to change and the fact that a dxm=0.5 was used, it
is likely that one (or more) of the hydrogen atoms went over the small
energy barrier, formed mainly by a combination of bonding harmonic force
(between H and O of water), van der Waals force between O and other
atoms. Now, the strong attractive electrostatic force it feels pulls it
toward that negatively charged particle (likely O of H2O and Cl-). The
relatively small harmonic force (bonds) is not sufficient to pull it
back. Consequently, the hydrogen atoms will go all the way to simply sit
on top of the particle and making an enormously negative electrostatic
energy.

In summary, the combination of NTC=1 and dxm=0.5 was probably the cause
of what you saw. Since NTC=1 is not consistent with the intended TIP3P
model, I would suggest you use the NTC=2 instead. Dxm=0.5 is quite
subjective. It may actually work if you have NTC=2, even though I
personally would prefer dxm=0.1.

Your ealier question regarding resizing the united cell is interesting.
I think Dave Case answered pretty well.
Even within the framework of classical dynamics, one may find it a
little hard to formulate a method to rescale the box size. This is
because the size is adjusted in reaction to the pressure. It is rather
difficult to "calculate" the pressure (based on virial expansion) at 0K
because, within the framework of classical dynamics, pressure arises
because of collision. But at 0K one should not expect collision (within
classical dynamics). A somewhat non-orthodoxical approach is to use the
virial term only (without the correction of the temperature factor). But
then, one has to justify this unphysical approach.

However, if you are really concerned about the equilibration of your
system, a good approach is to run a relatively short MD simulation at
low temperature (e.g., 10-100K). This allows the system to adjust the
box size and the approach is physically correct.

Hope this helps!

yong

-----Original Message-----
From: Ioana Cozmuta
Sent: Friday, February 21, 2003 8:11 PM
To: amber@heimdal.compchem.ucsf.edu
Subject: failure of minimization


Hi,

I have a box of ionic solution (KCl, 44 ions and about 1700 water
molecules, size of 40A). I did build this in Leap with some initial
random coordinates of the ions (something I thoughth was ok as I was
expecting the minimizer would be able to handle).

Here is my input file:

Minimization of the KCl cell
&cntrl
imin = 1, maxcyc = 500, ncyc = 200
ntpr = 1, ntx = 1
scnb = 2.0
scee = 1.2
ntf = 1, ntc = 1, igb = 0
ntr = 1
cut = 12.00
ntc = 1
&end
Group input for restrained atoms
1.0
RES 1 1756
END
END

I get the following message:

***** Processor 1 ***** System must be very
inhomogeneous. ***** Readjusting recip sizes. END In this slab, Atoms
found: 5180 Allocated: 3885

Initially my non-bonded energies are not a number but the minimizer
seems to be able to handle the vdw (it is reduced to -24.6510). However
the electrostatic energy remains NaN.

If I load the same initial structure in Cerius2 (accelrys software) and
I perform a minimization there, indeed I have a maximum force in the
system of about 10^20 kcal/mol/A and the energy component due to stress
in the order of 10^19kcal/mol. However after 500 steps of minimization I
get reasonable values (negative energies and max force ~10^1).

I would appreciate any suggestions on this.

Thank you in advance,
Ioana