next up previous contents
Next: Reference methods for ground Up: The electron force field, Previous: Electron force field makes   Contents

Wave packet molecular dynamics

In wave packet molecular dynamics, we represent nuclei as classical particles, and electrons as spherical Gaussian wave packets whose positions $ \mathbf{x}$ and extents $ s$ vary over time:

$\displaystyle \Psi \propto \prod_{j} \,\exp\left[-\left(\frac{1}{s^{2}}-\frac{2...
...(\mathbf{r}-\mathbf{x})^{2}\right]\cdot\exp[i \mathbf{p_{x}} \cdot \mathbf{x}].$ (3.6)

In a harmonic potential, Gaussian wave packets stay Gaussian over time, and it is meaningful to talk about the evolution of the coordinates $ x$, $ s$, and momenta $ p_{x}$ and $ p_{s}$. Substituting the wave packet into the time-dependent Schrodinger equation gives the Hamilton equations of motion (Appendix A)
    $\displaystyle \mathbf{\dot{p}_{R}} = -\nabla_{R}V, \ \mathbf{\dot{p}_{x}} = -\nabla_{x}V, \ \dot{p}_{s} = -\partial V/\partial s$  
    $\displaystyle \mathbf{p_{R}} = m_{nuc} \mathbf{R}, \ \mathbf{p_{x}} = m_{elec} \mathbf{x}, \ \mathbf{p_{s}} = (3 m_{elec}/4) s$ (3.7)

assuming a locally harmonic interaction potential. These equations can be viewed as a generalization of Ehrenfest's theorem [22], which states that the average position of a wave packet obeys classical dynamics, with the addition that the size of the wave packet obeys classical dynamics as well. The $ 3/4$ factor in front of the mass multiplying the radial coordinate is related to the dimensionality of the Gaussian packet, and becomes $ 2/4$ for a 2D Gaussian, $ 1/4$ for a 1D Gaussian, and so on.

Heller showed these equations could be applied profitably to anharmonic reaction potentials to describe processes such as collinear $ \mathrm{He + H_{2}}$ scattering [5]. Klakow later applied the same procedure to the anharmonic Coulomb potentials in hydrogen plasma and lithium metal, and found he was able to reproduce equilibrium pair distribution functions well [16]. He also extracted the conductivity of a hydrogen plasma by using a Green-Kubo expression [23,24] relating the fluctuation of current -- obtained from nucleus and electron positions -- to the dissipation of current given an applied potential. Indeed, since electrons are just another particle in these molecular dynamics simulations, we can use relations traditionally applied to nuclear positions and velocities to compute a wide range of electrical properties (Table 3.1).

Table 3.1: Quantities that can be calculated from wave packet molecular dynamics.
Nuclear temperature $ T_{nuc} = \frac{1}{3 k_{B}} \left<\sum_{\textrm{nuc d.o.f.}} \frac{1}{2} m_{i}\dot{x_{i}}^{2}\right>$
Electron temperature $ T_{elec} = \frac{1}{4 k_{B}} \left<\sum_{\textrm{elec d.o.f.}} \frac{1}{2} m_{i}\dot{x_{i}}^{2}\right>$
Pressure $ P = \frac{2\braket{K_{total}}}{V} - \frac{1}{3 V} \left< \sum_{\textrm{all d.o.f.}} x_{i} \frac{\partial E}{\partial x_{i}} \right>$
Current $ \mathbf{J} = \sum_{i\in\textit{nuc}} Z_{i} \dot{\mathbf{R}}_{i} + \sum_{i\in\textrm{elec}} (-1)\cdot \dot{\mathbf{R}}_{i}$
Conductivity $ \sigma = \frac{1}{3 V k_{b} T} \int \int \mathbf{J}(t) \cdot \mathbf{J}(t + \tau) dt d\tau$
Excitation response function $ \upsilon(\tau) = \int \sum_{i\in\textit{elec}} \dot{\mathbf{R}_{i}}(t) \dot{\mathbf{R}_{i}}(t+\tau) dt$
Excitation spectrum $ A(\omega) = \left\vert\int \upsilon(t) \cos(\omega t) dt\right\vert^{2}$
Electron density $ \rho(\mathbf{x}) = \sum_{i \in \textit{elec}} \left\vert\phi(\mathbf{x} - \mathbf{R}_{i}; r_{i})\right\vert^{2}$
Spin polarization $ \zeta(\mathbf{x}) = (\rho^{\uparrow}(\mathbf{x})-\rho^{\downarrow}(\mathbf{x}))/\rho(\mathbf{x})$

In our simulations, we usually set $ m_{electron} = m_{H}$, so that we may use a longer time step $ t = 0.1-0.5 fs$. There are additional reasons that can be given for this choice -- for example, anharmonic potentials tend to damp out radial oscillations in the wavefunction (Appendix A); the Landau theory of Fermi liquids uses heavy quasiparticles that obey fermion statistics [25] -- but it would be best to rerun our simulations with a smaller electron mass to make sure our choice does not have too large an adverse effect.

We are making an assumption of mean field dynamics [26,27] by propagating electron dynamics via the time-dependent Schrodinger equation. In our approach, the electron exists as a superposition of stationary states $ (\Psi_{i}, E_{i})$, and the nuclei move in a potential that is a linear combination of those states:

$\displaystyle E = \sum_{i}\vert c_{i}\vert^{2} E_{i}.$ (3.8)

There are conditions under which this is not a good assumption. For example, if our system has only a few well-separated states, and is prepared into a stationary state through a monochromatic light pulse, that state will tend to stay on one adiabatic path for a long period of time, and switch to other states via conical intersections. In these cases, a surface hopping stochastic dynamics scheme would be more appropriate, as mean field dynamics would split the difference between adjacent paths. A stochastic scheme is also appropriate for cases where the spacing between states changes rapidly during the course of the simulation, for example in the case of an electron scattering off a surface, where the electron goes from a continuum of states in free space to an insulator-like state on the surface.

Most cases we are interested in, however, contain many excited electrons and many closely-spaced states (Figure 3.3), where the entire dynamics of the system take place in a phase space packed with conical intersections with high Massey parameter [28] (see Appendix B). Here, the system is constantly jumping between adiabatic states, and the mean field trajectory is a good approximation of the nonadiabatic motions of the electrons.

Figure 3.3: (a) Excited condensed system evolves through many curve crossings, and can be approximately described by a mean field trajectory. (b) However, the mean field trajectory may incorrectly bisect well-separated adiabatic states.

Mean field dynamics has seen recent interest lately in combination with time-dependent density functional theory. In time-dependent DFT [29], we solve time-dependent Schrodinger equation

$\displaystyle i\hbar\frac{d}{dt}\psi_{j}(\mathbf{r},t) = \left(-\frac{\hbar^{2}}{2 m_{e}}\nabla^{2}+V_{s}[\rho](r,t)\right)\psi_{j}(\mathbf{r}, t)$ (3.9)

with a Kohn-Sham potential that is a functional of the time-varying electron density $ \rho = \sum_{j=1}^{N}\vert\psi_{j}(\mathbf{r}, t)\vert^{2}$:

$\displaystyle V_{s}[\rho] = V_{nuc}(\mathbf{r},t) + \int\frac{\rho(\mathbf{r'}, t)}{\vert\mathbf{r}-\mathbf{r'}\vert} d\mathbf{r'} + V_{xc}[\rho](\mathbf{r}, t).$ (3.10)

An adiabatic approximation is often applied, where it is assumed that the electron exchange-correlation does not depend on the past history of the electron density:

$\displaystyle V_{xc}[\rho](\mathbf{r}, t) = V_{xc}[\rho](\mathbf{r}).$ (3.11)

Even the simplest functional, the adibatic local-density approximation (ALDA), gives good excitation spectra [30], Rydberg states [31], and dispersion coefficients [32]. It has been applied to calculate the mean-field dynamics of sodium dimer [33], lithium cyanide ion [33], and ethylene [34] in response to femtosecond laser pulses, and the chemiadsorption of hydrogen on aluminum (111) surfaces [35].

Electron force field wave packet MD simulations are in some sense complementary to TDDFT simulations. Unlike ALDA, eFF has nonlocality in space, due to pairwise interactions of the electron force field, both in exchange and correlation; and nonlocality in time, due to the inertia of the electrons. eFF may serve as a useful method to investigate the significance of the space/time locality assumptions made in TDDFT methods.

next up previous contents
Next: Reference methods for ground Up: The electron force field, Previous: Electron force field makes   Contents
Julius 2008-04-29