John M. Millam

Ab initio molecular dynamics: Propagating the density matrix with Gaussian orbitals

We propose and implement an alternative approach to the original Car–Parrinello method where the density matrix elements (instead of the molecular orbitals) are propagated together with the nuclear degrees of freedom. Our new approach has the advantage of leading to an O(N) computational scheme in the large system limit. Our implementation is based on atom-centered Gaussian orbitals, which are especially suited to deal effectively with general molecular systems. The methodology is illustrated by applications to the three-body dissociation of triazine and to the dynamics of a cluster of a chloride ion with 25 water molecules.

Ab initio molecular dynamics: Propagating the density matrix with Gaussian orbitals. III. Comparison with Born-Oppenheimer dynamics

In a recently developed approach to ab initio molecular dynamics (ADMP), we used an extended Lagrangian to propagate the density matrix in a basis of atom centered Gaussian functions. Results of trajectory calculations obtained by this method are compared with the Born–Oppenheimer approach (BO), in which the density is converged at each step rather than propagated. For NaCl, the vibrational frequency with ADMP is found to be independent of the fictitious electronic mass and to be equal to the BO trajectory result. For the photodissociation of formaldehyde, H2CO→H2+CO, and the three body dissociation of glyoxal, C2H2O2→H2+2CO, very good agreement is found between the Born–Oppenheimer trajectories and the extended Lagrangian approach in terms of the rotational and vibrational energy distributions of the products. A 1.2 ps simulation of the dynamics of chloride ion in a cluster of 25 water molecules was used as a third test case. The Fourier transform of the velocity–velocity autocorrelation function showed ...

Linear scaling conjugate gradient density matrix search as an alternative to diagonalization for first principles electronic structure calculations

Advances in the computation of the Coulomb, exchange, and correlation contributions to Gaussian-based Hartree–Fock and density functional theory Hamiltonians have demonstrated near-linear scaling with molecular size for these steps. These advances leave the O(N3) diagonalization bottleneck as the rate determining step for very large systems. In this work, a conjugate gradient density matrix search (CG-DMS) method has been successfully extended and computationally implemented for use with first principles calculations. A Cholesky decomposition of the overlap matrix and its inverse is used to transform to and back from an orthonormal basis, which can be formed in near-linear time for sparse systems. Linear scaling of CPU time for the density matrix search and crossover of CPU time with diagonalization is demonstrated for polyglycine chains containing up to 493 atoms and water clusters up to 900 atoms.

Ab initio molecular dynamics: Propagating the density matrix with Gaussian orbitals. II. Generalizations based on mass-weighting, idempotency, energy conservation and choice of initial conditions

A generalization is presented here for a newly developed approach to ab initio molecular dynamics, where the density matrix is propagated with Gaussian orbitals. Including a tensorial fictitious mass facilitates the use of larger time steps for the dynamics process. A rigorous analysis of energy conservation is presented and used to control the deviation of the fictitious dynamics trajectory from the corresponding Born–Oppenheimer dynamics trajectory. These generalizations are tested for the case of the Cl−(H2O)25 cluster. It is found that, even with hydrogen atoms present in the system, no thermostats are necessary to control the exchange of energy between the nuclear and the fictitious electronic degrees of freedom.

Ab initio classical trajectories on the Born-Oppenheimer surface: Hessian-based integrators using fifth-order polynomial and rational function fits

Classical trajectories can be computed directly from electronic structure calculations without constructing a global potential-energy surface. When the potential energy and its derivatives are needed during the integration of the classical equations of motion, they are calculated by electronic structure methods. In the Born–Oppenheimer approach the wave function is converged rather than propagated to generate a more accurate potential-energy surface. If analytic second derivatives (Hessians) can be computed, steps of moderate size can be taken by integrating the equations of motion on a local quadratic approximation to the surface (a second-order algorithm). A more accurate integration method is described that uses a second-order predictor step on a local quadratic surface, followed by a corrector step on a better local surface fitted to the energies, gradients, and Hessians computed at the beginning and end points of the predictor step. The electronic structure work per step is the same as the second-ord...

Semiempirical methods with conjugate gradient density matrix search to replace diagonalization for molecular systems containing thousands of atoms

Conventional semiempirical methods using diagonalization are not practical for calculations on molecular systems containing more than a few hundred atoms because of O(N3) time and O(N2) memory requirements, where N is the number of atoms. Currently, the time dominating step is diagonalization of the Fock matrix. This paper demonstrates how O(N3) diagonalization and O(N2) memory requirements are eliminated by using a conjugate gradient search for the density matrix with sparse matrix techniques. Our method makes high accuracy energy calculations on molecules containing thousands of atoms possible on the typical workstation. Benchmark examples are presented on polyglycine chains (20000 atoms), water clusters (up to 1800 atoms), and nucleic acids (up to 6304 atoms).

Ab initio classical trajectories on the Born–Oppenheimer surface: Updating methods for Hessian-based integrators

For the integration of the classical equations of motion in the Born–Oppenheimer approach, each time the energy and gradient of the potential energy surface are needed, a properly converged wave function is calculated. If Hessians (second derivatives) can be calculated, significantly larger steps can be taken in the numerical integration of the equations of motion without loss of accuracy. Even larger steps can be taken with a Hessian-based predictor–corrector algorithm. Since updated Hessians are used successfully in quasi-Newton methods for geometry optimization, it should be possible to improve the performance of trajectory calculations using updated Hessians. The Murtagh–Sargent (MS) update, the Powell-symmetric–Broyden (PSB) update and Bofill’s update (a weighted combination of MS and PSB) were tested, and Bofill’s update was found to be the best. Slightly smaller step sizes were needed with Hessian updating to maintain good conservation of the energy, but this was more than compensated by the reduction in total computational cost. An overall factor of 3 in speed-up was obtained for trajectories of systems containing 4 to 6 heavy atoms computed at the HF/3-21G level.

Ab initio molecular dynamics studies of the photodissociation of formaldehyde, H2CO→H2+CO: Direct classical trajectory calculations by MP2 and density functional theory

The dynamics of H2CO→H2+CO photodissociation have been studied by classical trajectory calculations at the MP2/6-311G(d,p), B3LYP/6-311G(d,p), and VSXC/6-311G(d,p) levels of theory. The trajectories were calculated directly from the electronic structure computations without first fitting a global potential energy surface. A Hessian based method with updating was used to integrate the trajectories. The translational energy distribution of the products is in better agreement with experiment than the previous Hartree–Fock direct trajectory calculations, since the MP2 and density functional methods reproduce the barrier height better. The MP2 and density functional calculations give very good descriptions of the product rotational state distributions and the CO vibrational state populations. The MP2 method yields a very good representation of the H2 vibrational state populations, whereas the density functional methods give too little H2 vibrational excitation and Hartree–Fock produces too much. This can be at...

Ab initio molecular dynamics: Propagating the density matrix with gaussian orbitals. IV. Formal analysis of the deviations from born-oppenheimer dynamics

In the context of the recently developed Atom-centered Density Matrix Propagation (ADMP) approach to ab initio molecular dynamics, a formal analysis of the deviations from the Born—Oppenheimer surface is conducted. These deviations depend on the fictitious mass and on the magnitude of the commutator of the Fock and density matrices. These quantities are found to be closely interrelated and the choice of the fictitious mass provides a lower bound on the deviations from the Born—Oppenheimer surface. The relations are illustrated with an example calculation for the Cl−(H2O)25 cluster. We also show that there exists a direct one-to-one correspondence between approximate Born—Oppenheimer dynamics, where SCF convergence is restricted by a chosen threshold value for the commutator of the Fock and density matrices, and extended Lagrangian dynamics performed using a finite value for the fictitious mass. The analysis is extended to the nuclear forces used in the ADMP approximation. The forces are shown to be more general than those standardly used in Born—Oppenheimer dynamics, with the addition terms in the nuclear forces depending on the commutator mentioned above.

Density matrix search using direct inversion in the iterative subspace as a linear scaling alternative to diagonalization in electronic structure calculations

For electronic structure calculations on large systems, solving the self-consistent-field (SCF) equations is one of the key bottlenecks. Density matrix search methods provide an efficient linear scaling approach for circumventing the traditional O(N3) diagonalization procedure for solving the SCF equations. The conjugate gradient density matrix search (CG-DMS) method is a successful implementation of this approach. An alternative density matrix search method, QN–DMS, employs direct inversion in the iterative subspace using a quasi-Newton (QN) step to estimate the error vector. For linear polyglycine chains of 10–100 residues, the present approach scales linearly and is 30% faster than CG-DMS. For clusters of up to 300 water molecules, this method shows smoother and efficient convergence, and displays nearly linear scaling.