Extended Lagrangian molecular dynamics on vibronic surfaces in the nuclear-electronic orbital framework

Abstract

Proton transfer is central to many processes of chemical interest. The simulation of proton transfer dynamics requires the inclusion of nuclear quantum effects, such as zero-point energy, nuclear delocalization, and tunneling. Herein, we introduce methods within the nuclear-electronic orbital (NEO) framework, where specified nuclei are treated quantum mechanically on the same level as the electrons, for the simulation of proton transfer dynamics. Specifically, NEO density functional theory is used to treat the transferring protons quantum mechanically, and the other nuclei are propagated classically on the adiabatic vibronic ground-state surface. We formulate a NEO extended Lagrangian molecular dynamics (NEO-ELMD) approach to incorporate the motion of the nuclear basis function centers during such simulations. Density matrix extrapolation and purification are introduced as a means to accelerate the NEO self-consistent field procedure at each time step by reducing the number of iterations required for convergence. We demonstrate the fidelity and efficiency of NEO-ELMD by comparison to related dynamics methods for intramolecular proton transfer in malonaldehyde. We also use these accelerated techniques to simulate the nonequilibrium single and double proton transfer dynamics of proton-coupled electron transfer in much larger benzimidazole-phenol systems. This work provides a foundation for future methodologies to efficiently simulate proton transfer dynamics within the NEO-DFT framework while incorporating nonadiabatic effects between adiabatic vibronic states.

Figures (22)

• Figure 1: Schematic depictions of proton transfers simulated in this work, including the intramolecular proton transfer of malonaldehyde (left), the single proton transfer in E1PT BIP following oxidation (middle), and the double proton transfer relay in E2PT BIP following oxidation (right). Arrows point in the direction of proton transfer. Note that the phenolic proton transfer occurs first in the E2PT system, followed by the imidazolic proton transfer.
• Figure 2: NEO-BOMD malonaldehyde trajectory used as a reference for the NEO-ELMD$(K)$ and CNEO-MD$(K)$ trajectories. (A) Distance between the oxygen donor and proton position operator expectation value (blue) and proton basis function center (light blue) and distance between the oxygen acceptor and the proton position operator expectation value (red) and proton basis function center (light red) along the trajectory. (B) Distance between the donor and acceptor oxygens along the trajectory. (C) Evolution of the proton density (shown in cyan with an isosurface value of 0.04 a.u.) along the trajectory.
• Figure 3: NEO-ELMD$(4)$ and CNEO-MD$(4)$ trajectories compared to the NEO-BOMD trajectory for malonaldehyde. NEO-BOMD data is identical to the data shown in Fig. \ref{['fig:bomd_ref']}. (A) Distance between the oxygen donor and proton position operator expectation value (blue) and distance between the oxygen acceptor and the proton position operator expectation value (red) along the NEO-BOMD (solid lines), NEO-ELMD$(4)$ (dashed lines), and CNEO-MD$(4)$ (dotted lines) trajectories. The corresponding distances for the proton basis function center are shown in light blue and light red dashed lines for the NEO-ELMD(4) trajectory. (B) Distance between the donor and acceptor oxygens along the NEO-BOMD (solid line), NEO-ELMD$(4)$ (dashed line), and CNEO-MD$(4)$ (dotted line) trajectories.
• Figure 4: Number of iterations of the simultaneous DIIS NEO-SCF procedure needed to fully converge the electronic and protonic densities at each time step at varying values of $K$ for the (A) NEO-ELMD$(K)$ and (B) CNEO-MD$(K)$ malonaldehyde trajectories. $K=0$ indicates that the converged densities from the previous time step were used as initial guesses following purification. All orders of extrapolation tested significantly decrease the number of SCF iterations needed to converge for both the NEO-ELMD and CNEO-MD trajectories, often by a factor of two to four.
• Figure 5: Change in extended energy along the NEO-ELMD$(4)$ malonaldehyde trajectories with time steps 0.2 fs (magenta), 0.5 fs (orange), and 1.0 fs (green). As expected, the extended energy is better conserved as the time step decreases.