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 ith
particle displacement, ri(t), is
ri(t) =
ri(0) +
ri(1) cos Wt +
ri(2) cos 2Wt ,
where ri(0),
is the distorted equilibrium position of the
ith particle,
and the ri(n)'s
are the time-independent
amplitudes of the different harmonics. The force acting on the
ith 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
Fi(n) +
n2
W2
mi
ri(n) = 0,
n = 0, 1, 2.
The Fourier coefficients, Fi(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
qi(0),
qi(2), and
qi(1)
are now taken to
vary with time and obey the following fictitious equations of motion:
Qi(n) =
Fi(n) +
n2
W2
mi
qi(n) =
mi
d2qi(n)/dt2 ,
n = 0, 1, 2;
so that dynamical simulated annealing can be used to find the
equilibrium values for these quantities
ri(0),
ri(2), and
ri(1)
when all
Qi(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
qi(1) where
ri(1)=a while
qi(0),
qi(2),
and the initial velocities are set to zero.
(2)
Make an MD time step by solving the last equation and update
the quantities, qi(n)
and the
dqi(n)/dt
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 qi(n)
per particle for the
N particles is monitored. (Each object vibrates around its
equilibrium position
ri(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
dqi(n)/dt
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
F0(1)+W2ri(1)=0
for W with
r0(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