Car-Parrinello Monitor for More Robust Born-Oppenheimer Molecular Dynamics

Abstract

Born-Oppenheimer molecular dynamics (BOMD) is a promising simulation method for exploring the possible reaction pathways of a chemical system, but one significant challenge is the increased difficulty of converging the self-consistent field (SCF) calculation that often accompanies the breaking and forming of chemical bonds. To address this challenge, we developed an enhancement to the BOMD simulation method called the Car-Parrinello monitor (CPMonitor) that uses Car-Parrinello molecular dynamics (CPMD) to recover from SCF convergence failures. CPMonitor works by detecting SCF convergence failures in BOMD and switching to a CPMD Hamiltonian to propagate through the region of configuration space where the SCF calculation is unable to converge, then switching back to BOMD when good convergence behavior is re-established. We present a series of simulation studies that use CPMonitor, including detailed studies of the thermodynamic and dynamical properties of simple systems, as well as ab initio nanoreactor simulations containing transition metal atoms that were previously not possible to simulate using standard BOMD methods. Our studies show that CPMonitor can make BOMD simulations robust to SCF convergence difficulties and improve simulation performance and stability in reaction discovery applications.

Figure 1:

Flowchart illustrating main concepts of CPMonitor. Top row: CPMonitor running in normal mode (green) is the same as normal BOMD except the trajectory is written with a delay of n_S frames. Middle row: At the start of recovery mode (red), the simulation time is set back by n_S frames, the density matrix velocity is initialized, and the simulation continues with a reduced time step. Lower row: Watch mode (purple semicircle) is activated when CPMonitor is in recovery mode and the current time coincides with a BOMD time step. A SCF calculation is started and successful convergence increments a counter. When the number of consecutive successful SCF calculations exceeds a threshold, CPMonitor reverts to normal mode.

Figure 2:

Free energy profile of malonate ion along the hydrogen transfer coordinate.

Figure 3:

Power spectral density of malonate ion computed from energy-conserving trajectories, 100 ps in length, with initial velocities from a 750 K Maxwell-Boltzmann distribution. Error bars (transparency) indicate one standard error computed from 10 semi-independent runs.

Figure 4:

Comparison of CPMonitor simulation trajectory and reference BOMD trajectory with level-shifting. Top panel: HOMO and LUMO energy levels of the reference trajectory, showing the closing of the energy gap at t = 8480 fs. Middle panel: Difference in nuclear gradient between the CPMonitor and reference trajectory (blue, CPMonitor in normal mode; red, CPMonitor in recovery mode). The norm of the nuclear gradient of the reference trajectory (black) is shown for comparison. Bottom panel: RMSD of atomic positions between the CPMonitor and reference trajectory (blue, CPMonitor in normal mode; red, CPMonitor in recovery mode.)

Figure 5:

Distribution of the number of SCF cycles per AIMD time step (macro-iteration). The left column and right column show the same data using linear and log scales on the y-axis respectively. The kernel density estimates in the left column are guides for the eye. The circled bars in the right column indicate SCF calculations that did not converge.

Figure 6:

Center: 3-D rendering of a trajectory snapshot from the CPMonitor nanoreactor simulation. Side panels: Optimized structures of Fe-containing coordination complexes that were encountered over the course of the simulation trajectory.

Figure 7:

Organic molecules and radicals that were found in the nanoreactor simulations, grouped by empirical formula of heavy atoms. Each molecule is color coded by the simulation trajectories that it occurred in. The Venn diagram counts the number of molecules that were found in the different possible combinations of simulations.