AMBER Archive (2003)

Subject: RE: failure of minimization

From: Ioana Cozmuta (ioana_at_nas.nasa.gov)
Date: Mon Feb 24 2003 - 13:33:58 CST

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
>
>