Date: Wed, 12 Jun 2002 12:35:53 -0700 (PDT)
From: Thomas Cheatham
Subject: Re: vlimit exceeded


> I tried to run MD by constraint some atoms (NTR=1). I got an error message
> "vlimit exceeded" at beginning of dynamics. An input and output are as below
...
> Constrain G8
> 5000.0
> ATOM 236 250
> END
> END
...
> vlimit exceeded for step 2; vmax = 271.116268
>
> COORDINATE RESETTING CANNOT BE ACCOMPLISHED,
> DEVIATION IS TOO LARGE
> NITER, NIT, LL, I AND J ARE : 0 1 65 181 182

To follow up on Carlos Simmerling's answer with a little more information,
at issue is the time step and high frequency motion. The NTR=1 option
turns on a harmonic potential to restrain coordinates to a reference
frame. The higher the force constant, the tighter the restraint, but also
the higher the frequency of the effective vibration during dynamics
(thermal motion). With force constants greater than ~20-250 kcal/mol
(depending on the system) for positional restraints, it is not possible to
integrate these high frequency motions with a 2 fs time step. In
principle you could simply run with a smaller time step, but this is
impractical. However, large restraint force constants are not necessary
to maintain the structure since the harmonic potential is fairly steep so
large motion away from the reference coordinates is not really possible
without a very large (energetic) push. A value of 1.0-5.0 kcal/mol
for the restraint constant is typically fine for dynamics and will
maintain the structure.

[Note that this issue with dynamics/vibration is irrelevant for
minimization where you can use large force constants like 5000.0]

As a rule of thumb, generally you want the integration time step in MD
simulation, dt, to be dt ~< T / 10 where T is the shortest period (or
fastest vibration). So for example, H- stretches are in the > 3000 cm-1
range, whereas bends (N-H, C-H, etc) are in the 1000-1600 cm-1 range.

v ~ 3000 cm-1 x (3x10**10 cm/sec) = 90x10**12/sec so
T = 1/v = 0.011 ps or 11 fs

--> dt ~1 fs for bond stretches (unless SHAKE is applied)

v ~ 1600 cm-1 x (3x10**10 cm/sec) = 48x10**12/sec so
T = 1/v = 0.021 ps or 21 fs

--> dt ~ 2 fs for bends

Note that the force constants for bond vibration in parm94.dat are
~200-600 kcal/mol and for bending are ~25-75 kcal/mol. With positional
restraints due to the (effectively) larger instantaneously violations,
smaller force constants are typically recommended.

I hope this helps.