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