**Dynamical Simulated Annealing**

To construct the stationary localized mode eigenvector for the
nonlinear 3-D diatomic lattice with long range interactions, we build
on techniques which have been used successfully to identify 1-D
anharmonic modes [1], one of which is the rotating wave
approximation. A Fourier extension of that idea for vibration in a
stationary periodic mode with fundamental frequency, *W*,
which
includes the static and second harmonic for the **i**^{th}
particle displacement, **r**_{i}(t), is
**r**_{i}(t) =
**r**_{i}^{(0)} +
**r**_{i}^{(1)} cos *Wt* +
**r**_{i}^{(2)} cos 2*Wt* ,
where **r**_{i}^{(0)},
is the distorted equilibrium position of the
**i**^{th} particle,
and the **r**_{i}^{(n)}'s
are the time-independent
amplitudes of the different harmonics. The force acting on the
**i**^{th} particle can also be represented
by a similar Fourier series.
Substituting the coordinate and force Fourier series into the
classical equations of motion and equating the terms with different
harmonics gives a time-independent system of nonlinear equations
**F**_{i}^{(n)} +
*n*^{2}
*W*^{2}
*m*_{i}
**r**_{i}^{(n)} = 0,
*n* = 0, 1, 2.
The Fourier coefficients, **F**_{i}^{(n)},
are determined in the usual way [1].
The procedure for finding the anharmonic localized mode eigenvector
relies on first generating a fictitious dynamics for the Fourier
amplitudes and then applying a version of dynamical simulated
annealing [2]. The coefficients
**q**_{i}^{(0)},
**q**_{i}^{(2)}, and
**q**_{i}^{(1)}
are now taken to
vary with time and obey the following fictitious equations of motion:

**Q**_{i}^{(n)} =
**F**_{i}^{(n)} +
*n*^{2}
*W*^{2}
*m*_{i}
**q**_{i}^{(n)} =
*m*_{i}
d^{2}**q**_{i}^{(n)}/dt^{2} ,
*n* = 0, 1, 2;
so that dynamical simulated annealing can be used to find the
equilibrium values for these quantities
**r**_{i}^{(0)},
**r**_{i}^{(2)}, and
**r**_{i}^{(1)}
when all
**Q**_{i}^{(n)}=0.
Our iterative procedure to obtain the eigenvector is
as follows:

(1)
Guess an initial shape for the gap mode eigenvector. In our case we
use the mass-defect eigenvector associated with a harmonic gap mode of
frequency *W*.
The rescaled mass-defect gap mode eigenvector
gives initial amplitudes
**q**_{i}^{(1)} where
**r**_{i}^{(1)}=*a* while
**q**_{i}^{(0)},
**q**_{i}^{(2)},
and the initial velocities are set to zero.
(2)
Make an MD time step by solving the last equation and update
the quantities, **q**_{i}^{(n)}
and the
d**q**_{i}^{(n)}/d*t*
for fixed central particle amplitude *a*.
The classical equations of motion are integrated using the
"leap-frog" algorithm [3] with a time step of 1.35 fs.
(3)
Apply a form of dynamical simulated annealing [2]. During a MD
simulation run, the oscillatory values of the kinetic energies of each
of these three objects **q**_{i}^{(n)}
per particle for the
*N* particles is monitored. (Each object vibrates around its
equilibrium position
**r**_{i}^{(n)}).
When the total kinetic
energy of the system of objects passes through its maximal value in
the MD fictitious time-evolution, all object's velocities
d**q**_{i}^{(n)}/d*t*
are set to zero.
In this way we
incrementally move the system closer to its equilibrium configuration.
(4)
After 40 MD steps as described in (2) and checking for
condition (3) the intrinsic gap mode frequency *W* is
updated by solving the single equation
*F*_{0}^{(1)}+*W*^{2}*r*_{i}^{(1)}=0
for *W* with
*r*_{0}^{(1)}=*a*.
(5)
Verify the correctness and stability of the resultant IGM eigenvector
with a regular MD simulation. Repeat (1)
through (4) until the lifetime of the mode remains unchanged
when the resultant eigenvector is used as an initial condition for an
MD simulation. As long as the frequency of the mode remains in the
gap we find the procedure described here converges to a localized
eigenvector. No change in the results is observed when this procedure
is repeated with a 1/2 time step.
A system of up to 9,000 nonlinear
equations is solved successfully using the developed method. This
numerical technique is general enough to apply to any nonlinear system
which allows a classical molecular dynamics treatment, including the
rigid-ion and the shell models.

[1] S.A.Kiselev, S.R.Bickham, and A.J.Sievers, Phys. Rev. B
**48**, 13508 (1993).

[2] R.Car and M.Parrinello, Phys. Rev. Lett. **55**, 2471 (1985).

[3] R.W.Hockney and J.W.Eastwood,
*Computer Simulation Using Particles*
(McGraw-Hill Inc., New York, 1981).

Last modified: April 3, 1997