Direct molecular dynamics Hamiltonian equations

For a complete treatment of a laser-driven molecule, one must solve the many-body, multidimensional time-dependent Schrodinger equation (TDSE). This represents a tremendous task and direct wavepacket simulations of nuclear and electronic motions under an intense laser pulse is presently restricted to a few bodies (at most three or four) and/or to a model of low dimensionality [27]. For a more general treatment, an approximate separation of variables between electrons (fast subsystem) and nuclei (slow subsystem) is customarily made, in the spirit of the BO approximation. To lay out the ideas underlying this approximation as adapted to field-driven molecular dynamics, we will consider from now on a molecule consisting of Nn nuclei (labeled a, p,...) and Ne electrons (labeled /, j,...), with position vectors Ro, and r respectively, defined in the center of mass (rotating) body-fixed coordinate system, in a classical field E(f) of the form Eof t) cos cot). The full semiclassical length gauge Hamiltonian is written, for a system of electrons and nuclei, as [4]... 

At its most basic level, molecular dynamics is about mapping out complicated point sets using trajectories of a system of ordinary differential equations (or, in Chaps. 6-8, a stochastic-differential equation system). The sets are typically defined as the collection of probable states for a certain system. In the case of Hamiltonian dynamics, they are directly associated to a region of the energy landscape. The trajectories are the means by which we efficiently explore the energy surface. In this chapter we address the design of numerical methods to calculate trajectories. 

Here, it is usual to make the Bom-Oppenheimer approximation that allows a classical treatment of the nuclei to be separated from a quantum mechanical description of the electrons. In this case, the wave function becomes just that of the electrons, and the nuclear-nuclear interaction is added to the energy as a sum over point particles. Consequently, the Hamiltonian operator H includes the kinetic energy of the electrons, the electron-electron interactions, and the electron-nuclei interactions. The wave function determined by solving this eigenproblem consists of a Slater determinant of the molecular orbitals for a molecule or, alternatively, the band structure of a solid. Unfortunately, direct solution of this equation is complicated by the electron-electron interactions. Often, it is necessary to introduce a mean-field approximation that neglects the individual dynamical electron-electron correlations but instead treats the electrons as moving in the average field created by the other electrons. Various corrections have been developed to improve upon this approximation [160, 167, 168]. 

Both Monte Carlo and molecular dynamics methods sample directly the phase space of a small but representative component of the crystal, the former by performing stochastic moves through configuration space, the latter by following a specified trajectory according to an equation of motion and chosen initial condition. A typical Hamiltonian for molecular dynamics simulation is [14] ... 

The QM/MM Hamiltonian can be used to cany out Molecular Dynamics simulations of a complex system. In the case of liquid interfaces, the simulation box contains the solute and solvent molecules and one must apply appropriate periodic boundary conditions. Typically, for air-water interface simulations, we use a cubic box with periodic boundary conditions in the X and Y directions, whereas for liquid/liquid interfaces, we use a rectangle cuboid interface with periodic boundary conditions in the three directions. An example of simulation box for a liquid-liquid interface is illustrated in Fig. 11.1. The solute s wave function is computed on the fly at each time step of the simulation using the terms in the whole Hamiltonian that explicitly depend on the solute s electronic coordinates (the Born-Oppenheimer approximation is assumed in this model). To accelerate the convergence of the wavefunction calculation, the initial guess in the SCF iterative procedure is taken from the previous step in the simulation, or better, using an extrapolated density matrix from the last three or four steps [39]. The forces acting on QM nuclei and on MM centers are evaluated analytically, and the classical equations of motion are solved to obtain a set of new atomic positions and velocities. 

The classical microscopic description of molecular processes leads to a mathematical model in terms of Hamiltonian differential equations. In principle, the discretization of such systems permits a simulation of the dynamics. However, as will be worked out below in Section 2, both forward and backward numerical analysis restrict such simulations to only short time spans and to comparatively small discretization steps. Fortunately, most questions of chemical relevance just require the computation of averages of physical observables, of stable conformations or of conformational changes. The computation of averages is usually performed on a statistical physics basis. In the subsequent Section 3 we advocate a new computational approach on the basis of the mathematical theory of dynamical systems we directly solve a... 

---