Trajectory Surface Hopping for a Polarizable Embedding QM/MM Formulation

Abstract

We present the implementation of trajectory surface-hopping nonadiabatic dynamics for a polarizable embedding QM/MM formulation. Time-dependent density functional theory was used at the quantum mechanical level of theory, whereas the molecular mechanics description involved the polarizable AMOEBA force field. This implementation has been obtained by integrating the surface-hopping program Newton-X NS with an interface between the Gaussian 16 and the Tinker suites of codes to calculate QM/AMOEBA energies and forces. The implementation has been tested on a photoinduced electron-driven proton-transfer reaction involving pyrimidine and a hydrogen-bonded water surrounded by a small cluster of water molecules and within a large water droplet.

Figure 1

Schematic representation of Newton-X-Tinker-Gaussian implementation. Newton-X acts as a driver providing at the beginning of each step the coordinates of the system to Tinker that in turn calls Gaussian to compute energies and gradients with TD-DFT/AMOEBA PES. Then, to evaluate the transition probabilities between states, the overlap of wave functions at two subsequent time steps are computed by a program written on purpose (ci ovlp in the scheme, it is provided together with Newton-X) that uses the orbital overlap matrix computed using Gaussian.

Figure 2

Graphical representation of the clusters Pm(H₂O)₄ (left); hydrogen bonding is represented with dashed yellow lines, molecules in the MM region are represented with balls and sticks, while those in the QM region are represented as licorice. Relevant transitions for the EDPT process represented as single excitations between molecular orbitals at the equilibrium geometry of Pm(H₂O)₄. MOs are computed with QM/AMOEBA Hamiltonian and plotted as isosurfaces at +0.01 and −0.01 in orange and blue, respectively.

Figure 3

H–N distance computed along TSH trajectories of QM/TIP3P (left), QM/AMOEBA (center), and Pm(H₂O)₁ (right). Lines corresponding to EDPT events are highlighted; final cross indicates that the trajectory is terminated before 100 fs.

Figure 4

Comparison of a nonreactive (left) and a reactive QM/AMOEBA trajectory (right). (Top) Excitation energies of the low-lying states colored in terms of the intensity of their oscillator strength (see color bar inside the picture); state with green edges is the current state of TSH, vertical lines indicate the hopping events, and gray shadow is the region where the simulation is stopped and considered as relaxed on the GS. (Middle) Same as the top one but the color maps the character of the states. (Bottom) Distance of Pm N atom from the hydrogen atom of the initially H-bonded water molecule.

Figure 5

Energies and CT character of the lowest excited states computed along the QM/AMOEBA reactive trajectory shown in Figure 4 with different representations of environment water molecules: omitted (upper left), LR QM/AMOEBA (upper right), QM/TIP3P (bottom left), and cLR² QM/AMOEBA for each state (bottom right). In the latter case, only the lowest 10 excited state have been analyzed.

Figure 6

Pyrimidine water droplet (top right). At the left, analysis of a reactive trajectory for the droplet. Format is the same as that adopted in Figure 4. Bottom right shows the CT character of the first 30 fs of the same trajectory but considering cLR² AMOEBA.