Questions and problems?

BLOWUP

This section talks about possible cases where SHAKE failures occur and/or the box seems to blow up or bonds stretch way too far.

NOTE: See related material in linmin and equilibration.

Sometimes right when I start my periodic boundary condition (PBC) simulation, I get SHAKE problems or the box expands...

Pre-4.1 EDIT did not do a good job of trimming waters at the faces of the periodic box, so waters could overlap to the point of virtual superimposition at the box boundary. This leads to strong enough forces to blow things apart. As of 4.1 EDIT, .4A 'breathing' room has been added to the box, so this problem should become relatively rare. One can achieve the same effect with earlier versions of AMBER by manually editing the box coordinates at the end of the prmtop file, adding .4A or so to X, Y, and Z. Note that whatever the degree of excess (or deficient) 'packing' at the box boundary, correct density is likely to be obtained only by running constant pressure dynamics.

Minimization (even as little as 50 steps) can resolve many steric overlap situations to the point where dynamics will start. Hint: turn off SHAKE in the minimization only! to proper density.

When I do a dynamics run with constant pressure PBC's everything runs O.K. for about 5 psec, then I get the following message:
        COORDINATE RESETTING CANNOT BE ACCOMPLISHED,
	     DEVIATION IS TOO LARGE
	     NITER, NIT, LL, I AND J ARE :    0    0 5533 5799 5800

	 %MINMD-F-ERSTOP, fatal error
This seems to happen only with constant pressure PBC's, and needless to say it is annoying. Can you tell me what the error message means, and what to do to get around the problem?

Coordinate resetting is done by SHAKE to keep constant bond lengths, etc. Constant pressure is achieved by changing the box size. All the coordinates get scaled and then shake has to repair the bond distances. If the scaling is excessive, then shake can't cope. One possibility is that scaling blows up when running with IFTRES=0 (no cutoff applied to solute nonbonded interactions) and a counterion (i.e. solute atom) crosses the box boundary -- i.e. jumps into the other side of the box, where it happens to suddenly see another solute atom. The resulting jump in virial causes major expansion in the box, stretching bonds too far for shake to recover. Assuming you have IFTRES=0, you could watch your trajectory (e.g. using Terry Lybrand's MD Display program on SGI after stripping the waters with carnal) and see if a solute atom is translated across the box boundary right before the blowup. Having a bigger box is the only solution I know of for this problem, and it just puts off the blowup, which is inevitable with IFTRES=0 and solute moieties are mobile, if you think about it.

From the manual... 

IFTRES..

            = 0  ALL intramolecular solute - solute nonbonded interactions
                 are calculated regardless of whether the interatomic distance
                 is greater than the nonbonded cutoff.  Solute-solute imaging
                 is turned off.

             NOTE: For simulations of highly charged solutes in a water
             bath, it can be useful to calculate ALL solute - solute
             nonbonded interactions in order to reduce electrostatic
             problems.  This is especially important for highly charged
             systems like nucleic acids.  Note that this option is
             intended for small solutes, and will generate many more
             nonbonded pairs than the normal method if the solute is
             large.   Counterions added in EDIT are considered part of
             the solute.
Other than that, you may just have not equilibrated carefully enough, e.g. by warming gradually. Look at the pressure fluctuations leading up to the problem. The pressure algorithm is not perfect, and pressure fluctuates from -300..300 atm even in a well-equilibrated system. 
Subject: Re: gibbs div-by-0
To: ross@cgl.ucsf.EDU (Bill Ross)
Date: Tue, 13 Oct 92 14:15:11 EDT

Bill Ross writes:
> 
> I'm starting w/ a restrt from a sander run w/ 1 fs steps & setting
> DT=2fs; things start to blow up immediately & after a while I get 
> floating divide-by-0 in nnbond. Is the timestep change a likely
> cause? Also, does this mean anything given dt=.002?
> 
>      DTE   =    .00100     DTA   =    .00100
> 
> Bill
> 
>  &cntrl
>   nrun=9999, nstlim=600,
>   ntx=7, ntt=1, temp0=300.0, tautp=0.2, npscal=1,
>   dt=0.002, nsnb=25, idiel=1, ntb=2, ntp=1, iftres=0,
>   isldyn=-2, almda=1.0, nstmeq=200, nstmul=400, idifrg=0,idwide=1,
>   ntpr=1, ntwx=200, cut=8.0, intprt=0, init=4,
>  &end
>
You have obviously turned off SHAKE in going from md to Gibbs. (See how low the bond energy is on the first step and how much bigger it is on the second). That, coupled with changing the timestep, is probably the culprit. You should not try to run timesteps of 0.002 without any SHAKE. 
> I've encountered a strange phenomenon using Gibbs4 and I 
> don't have an idea what's happening.
> ...the simulations run fine and are stable, at least that's
> what the energies show, but after ~ 90% of the simulation the
> energies explode within 2-3 hundreds of steps and the whole
> systems flies apart. This happens both in vacuum and with 
> periodic boundaries: very funny to watch your water vaporize ...
You're seeing instabilities in the MD integration. These can result in wild oscillations in the system, and such phenomena as the waters apparently "vaporizing"...

The fact that the instability arises near an endpoint of the simulation suggests two possibilities:

  1. Instabilities arising from charges on atoms which are disappearing (i.e. whose vdw parameters are approaching zero). This is a well-known potential problem, and is usually avoided by using an "electrostatic decoupling" approach within Amber. This problem can sometimes be exacerbated by using a very small "mixed van der Walls" parameter, as specified by the IDSX0 input parameter. Sometimes, setting the mixed "vanised" vdW distance to a larger value will avoid the instability.
  2. Instabilities arising from shrinking bonds. If you attempt to shrink the bonds to dummy atoms down to very small values (< ~0.2-0.4 A), you can encounter instabilities near the endpoint of the simulation. Such instabilities can be avoided by either increasing, somewhat, the shortest bond length, or else by decreasing the timestep.
Decreasing the timestep will often avoid allow the simulation to be run without instability. But this can mean reducing the timestep to 0.5 fsec from the 1.0-2.0 fsec timesteps usually used. That will make the calculation more expensive. One approach is to run most of the simulation with the larger timestep, then restart the calculation with a small timestep as you approach the "lambda" region where instabilities are a problem. (dap)

SHAKE blowup in a box of n-octanol

> while trying to run a plain dynamics of 189 n-octanol molecules,
> I always get the error message
> COORDINATE RESETTING CANNOT BE ACCOMPLISHED
What is your time step? Are you evaluating all solute-solute interactions? I wonder if all-atom representation would be appropriate. Is each octanol a single residue, or multiple? If you have a long molecule viewed as a single residue, and all solute-solute interactions are not evaluated or the cutoff allows some 1st-atom/1st-atom distances to fail the residues for pairlist inclusion when last atoms are actually close, sudden high forces could result when the 1st atoms come within cutoff and the (even overlapping) 'tail' atoms suddenly 'see' each other and voila! a hydrogen is suddenly propelled into infinity beyond the ability of shake to recover it. It is tempting to try IFTRES=0 (all solute pairs evaluated), except in an all-solute periodic system, the solute doesn't see the periodic images, so you'd get the same sort of blowup when an octanol is translated across the box. Dividing the molecule up into multiple residues or using a very long cutoff is the solution, (the former saves compute time).
I am using SHAKE on H's, and my other bonds are growing very long, so my molecule looks like a spider monster from space.

Earlier, you were probably running with SHAKE on all bonds, and then switched to SHAKE on H only. At that point, most likely you changed NTC from 3->2 but forgot to update NTF similarly so that forces on other bonds would be evaluated.