AMBER Archive (2003)

Subject: RE: failure of minimization

From: Yong Duan (yduan_at_udel.edu)
Date: Mon Feb 24 2003 - 19:42:14 CST

Dear Ioana,

You probably need longer steepest descent minimization before switching
to conjugate gradient.

I've never done simulation with SPC/E. So, my comment may not be that
helpful. But I would try. ntc=2 to minimize for 100-500 steps (SD).
Then, minimize again with NTC=1 (with 100-500 steepest descent). I have
the impression that dxm is not effective in the conjugated gradient
minimization. In general, minimization does not give you "good"
structure, even for systems like yours. But it can allow the system to
relax a bit conservatively to remove the close contacts which can lead
to a lot of instability in the system.

Good luck!

yong

-----Original Message-----
From: Ioana Cozmuta [mailto:ioana_at_nas.nasa.gov]
Sent: Monday, February 24, 2003 2:34 PM
To: Yong Duan
Cc: amber_at_heimdal.compchem.ucsf.edu
Subject: RE: failure of minimization

Hi Yong,

Thank you for your very elaborate e-mail message, I appreciate it!
However, although I intended to write my e-mail to the list very
explicit, I forgot to mention that I do not use the standard TIP3/P
model but the SPC/E model for water. So far I only worked with flexible
water models so I am afraid that the reason why my minimization is
crashing is related to that. I did try to use various values for dmx in
my input file and it does not change anything. Also it seems that when I
copied the input file the ntc = 2 seems to be left out. I did try to use
ntf=1 and ntf=2 but it does not change anything. So unless I have a good
minimized structure I can not run any dynamics, of course. I would
appreciate any help on this matter, I feel that I am left out of
options.

Thank you,
Ioana

On Fri, 21 Feb 2003, Yong Duan wrote:

>
> 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 [mailto:ioana_at_nas.nasa.gov]
> Sent: Friday, February 21, 2003 8:11 PM
> To: amber_at_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
>
>