Polarizable Embedding for Excited-State ... - ACS Publications

Article Cite This: J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

pubs.acs.org/JCTC

Polarizable Embedding for Excited-State Reactions: Dynamically Weighted Polarizable QM/MM Muhammad A. Hagras†,∥ and William J. Glover*,‡,§,∥,⊥ †

NYU Shanghai, 14 East Fourth Street, New York, New York 10003, United States NYU Shanghai, 1555 Century Avenue, Shanghai 200122, China § NYU−ECNU Center for Computational Chemistry at NYU Shanghai, 3663 Zhongshan Road North, Shanghai 200062, China ∥ Department of Chemistry, New York University, New York, New York 10003, United States ⊥ Department of Chemistry, East China Normal University, Shanghai 200062, China ‡

S Supporting Information *

ABSTRACT: Recent interest in polarizable embedding methods for electronic excited states has so far been focused on optical absorption and emission spectra calculations. To explore the suitability of these methods for excited-state reactions, we constructed a simple molecular system with an electronic crossing coupled to a polarizable species: the triatomic LiFBe. We found that current polarizable QM/MM methods inadequately describe the potential energy surfaces in this system, particularly close to the electronic crossing, so we developed a new polarizable embedding method called dynamically weighted polarizable QM/MM. The new method reproduces the potential energy surfaces of LiFBe from fullsystem multireference configuration interaction singles and doubles calculations with near-quantitative accuracy.

1. INTRODUCTION Excited-state dynamics are important in many physical and chemical systems of interest, including the initial events of vision,1−3 photoswitches,4−6 and fluorescent proteins.7−9 Despite the complexity of these systems, they share a common theme that the light-absorbing moiety, the chromophore, is localized to a small region of the system, which naturally lends these systems to theoretical study with embedding methods such as quantum mechanics/molecular mechanics (QM/ MM),10 wherein the chromophore (and perhaps nearby protein residues or solvent molecules)11 are treated at the QM level and the remaining system is treated with classical MM force fields. Traditional QM/MM approaches couple the QM and MM subsystems through electrostatic interactions, assuming fixed charges in the MM environment; however, this scheme neglects electronic polarization of the environment, which can result in errors in electronic excitation energies that exceed 1 eV.12 Errors of this magnitude mean that using nonpolarizable QM/MM could result in predicted state orderings that are wildly incorrect. There has therefore been significant recent interest in extending polarizable QM/MM (QM/MM-Pol) to treat electronic excitations, since this should capture the electronic polarization response of the environment to an electronic excitation in the QM subsystem. In a QM/MM-Pol calculation, typically the MM region is allowed to have induced dipoles that are solved self-consistently with the QM calculation, such that the QM and MM regions mutually polarize each other. QM/MM-Pol methods applied to © XXXX American Chemical Society

ground-state reactions (or reactions on a single excited state) are now quite mature, with the first such calculation dating back to the very first 1976 QM/MM paper;10 however, this is not the case for photochemical reactions involving potential energy curve crossings, for which, to our knowledge, only one study has been published to date.12 A possible reason for the lack of a wider adoption of excited-state QM/MM-Pol methods is that while in single-state reactions it is clear how the QM and MM regions should mutually polarize, this is less clear when electronic states cross. In particular, multiple electronic states, each with different charge distributions, must be considered simultaneously, so the approach of polarizing the MM region to a single electronic state of the QM region is likely inadequate. Of particular concern is the description of adiabatic state crossings, i.e., conical intersections, which are now understood to be vitally important in dictating the mechanism and kinetics of photochemical reactions.13 Of the published ab initio QM/MM-Pol methods for excited states to date, we find that they can be grouped into three approaches. First are polarizable QM/MM methods that treat one state (the reference state) differently from the others. In this approach, the MM-Pol subsystem is polarized selfconsistently with the electric field of just a single reference electronic state, usually taken to be the ground state, and then electronic excitations are computed a number of different ways. Received: January 21, 2018 Published: March 21, 2018 A

DOI: 10.1021/acs.jctc.8b00064 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

We found that none of the previous QM/MM-Pol methods are able to accurately describe simultaneously the PES away from and in the vicinity of the electronic crossing in this system, so we developed a new method that we call dynamically weighted polarizable QM/MM (QM/MM-DW-Pol), which can be considered a hybrid of QM/MM-SS-Pol and QM/MM-SAPol methods. The method is described in detail below, and we show that it accurately reproduces the PES of LiFBe both at and away from the state crossing in this system.

The most straightforward approach computes electronic excitations with the MM-Pol dipoles frozen to the reference state, which we shall indicate as MM-Fr-Pol.14 This approach can be expected to be accurate only when the charge densities of the reference and excited states are similar, since it completely neglects the polarization response of the MM subsystem to a QM excited state, but these effects can be included perturbatively if needed.14,15 A second common approach is to treat the excitation of the MM-Pol subsystem with linear response (QM/MM-LR-Pol),16,17 which is applicable when the QM electronic excitations are also treated with linear response (e.g., with linear-response time-dependent density functional theory) and has been shown to give good agreement with full-system linear-response excitation energies.17 The disadvantage of reference-state QM/MM-Pol methods is the unequal treatment of electronic states: in particular, one expects artifacts in the potential energy surfaces near electronic crossings between the reference state and other states18a point we explore in this paper. A second method is state-specific polarizable QM/MM (QM/MM-SS-Pol), wherein the polarization of the MM-Pol subsystem is optimized separately to each state, one at a time, in a self-consistent fashion.19,20 This method has the advantage of correctly describing the polarization of the MM-Pol subsystem to each electronic state in a balanced fashion, and it yields good agreement with excitation energies from fullsystem calculations.21 However, in systems with near, or actual, electronic degeneracies, QM/MM-SS-Pol has been reported to have convergence issues, such as root flipping.22 Furthermore, because it is a state-specific approach, one expects that the topology of conical intersections will not be described correctly. Another issue with QM/MM-SS-Pol is that each electronic state arises from a different Hamiltonian, meaning that transition matrix elements between the states are ill-defined,23 although Jacobson and Herbert24 have developed an orthogonalization approach for one-electron Hamiltonians that overcomes this limitation. The third excited-state polarizable QM/MM method, QM/ MM-SA-Pol, uses a state averaging approach, such that the MM-Pol subsystem is optimized self-consistently to a state average of the QM electric fields.22 Although not as widely used as the first two QM/MM methods, QM/MM-SA-Pol has the advantage of ameliorating root-flipping problems, and the electronic states come from a single Hamiltonian, leading to well-defined transition properties. The drawback to the method is that the MM-Pol subsystem does not polarize differently to different electronic states, so computed excitation energies are not in good agreement with full-system calculations; however, Li et al.22 developed a state-specific correction that improves the agreement. As mentioned above, the bulk of the applications of excitedstate polarizable QM/MM have been to vertical electronic excitations, absorption spectra, and emission spectra of solutes in polarizable environments. Before QM/MM-Pol is applied to excited-state dynamics, it is important to benchmark the methods also on more global features of the potential energy surface, in particular their descriptions of electronic state crossings. Thus, in this work we constructed a simple molecular system to test several QM/MM-Pol methods against high-level theory, in particular, multireference configuration interaction singles and doubles (MRCISD). We chose the prototypical charge-transfer (CT) reaction of lithium fluoride interacting with a single polarizable beryllium atom: the triatomic LiFBe.

2. THEORY 2.1. Variational Formulation of State-Specific Polarizable QM/MM-Pol. Before introducing our dynamically weighted polarizable QM/MM method, we first briefly review the working equations for QM/MM-SS-Pol, which in the case of a single electronic state corresponds to traditional QM/MMPol. Without loss of generality, we use the Drude oscillator model of linear polarization, wherein induced dipoles are modeled as pairs of point charges separated by a variable distance.25,26 We chose this model because it provides a particularly simple scheme to couple the QM and MM-Pol subsystems through the Drude point charges. However, the QM/MM-DW-Pol method that we have developed is not limited to the Drude Oscillator Model, and one could go beyond linear polarization by including higher-order induced multipoles as additional variational parameters27 (at the expense of a more complicated QM/MM-Pol coupling). In the Drude oscillator model, one point charge, qα, is fixed to the charge center, while the second charge (called a Drude particle), qα′, is attached the charge center site via a harmonic spring. The total charge and dipole moment of each polarizable site are then given by qα + qα′ and μ = qα′|rα′ − rα|, respectively. The harmonic potential between the charge site and the Drude particle is chosen to give the correct dipole self-energy term: Eself

1 = 2

∑ α′

qα2′ αα ′

|rα ′ − rα|2 (1)

where αα′ is the isotropic polarizability of the site. In QM/MM-Pol, the electrostatic interaction between the QM and MM-Pol systems is introduced in the one-electron Hamiltonian of the QM system, electrons core,QM/MM ‐ Pol Hμν = −⟨μ|

∑ i

⎡ MM q ⎢∑ α + ⎢⎣ α riα

MM ′

∑ α′

qα ′ ⎤ ⎥| ν ⟩ riα ′ ⎥⎦ (2)

where μ and ν are atomic orbitals, such that the QM energy, EQM, includes the electrostatic interaction between the QM and MM systems. Following an additive scheme, the total potential energy of the system is given by Etot = EQM + EMM + EQM/MM + EMM ‐ Pol + Eself

(3)

where EMM is the nonelectrostatic potential energy of the MM subsystem (i.e., bond stretching, bending, and torsion potentials, short-range repulsion, and dispersion interactions), EQM/MM describes the nonelectrostatic interactions between the QM and MM subsystems (i.e., short-range repulsion and dispersion interactions), and EMM‑Pol is the electrostatic potential of the MM subsystem, including induced polarization: B


Article

Journal of Chemical Theory and Computation MM MM

EMM ‐ Pol =

∑∑ α

β>α

qαqβ rαβ

MM ′ MM ′

+

∑ ∑ α′

β ′> α ′

MM MM ′

+

∑∑ α

β ′≠ α

qαqβ ′

wI ({E QM}) =

rαβ ′

qα ′qβ ′ rα ′ β ′

(4)

rα ′ = rα −

(7)

⎛ ⎜ ∇ r ⎜E MM ‐ Pol + qα2′ α′⎝

αα ′

Nstate

⎞

I

⎠

∑ wIEIQM⎟⎟

(8)

As we show in the Supporting Information, the gradient terms in eq 8 are analytic and involve a computational cost similar to that of a QM/MM-SS-Pol optimization, albeit with the additional requirement of solving Nstate QM gradients rather than just the reference-state gradient. In this work, we follow Deskevich et al.30 and choose the weight function in eq 7 to be a hyperbolic secant squared function: f (E) = sech2(βE). This function satisfies the following three conditions: (1) the weights are smooth and analytic functions of energy (which is necessary for smooth and analytic potential energy surfaces); (2) when two states cross, their weights are equal, and therefore, the induced dipoles for the two states are identical at their intersection, so the states arise from the same Hamiltonian; and (3) when the reference state is energetically separated from the other states (by more than 2/ β), its weight approaches unity, and the QM/MM-DW-Pol method reduces to the original QM/MM-SS-Pol approach. Of course, other functional forms would satisfy these requirements; however, we found that the hyperbolic secant squared function gave satisfactory performance, and we leave a more detailed exploration of functional forms to future work. The choice of the parameter β is explored in the Supporting Information, but it is clear from eq 6 that in the limit β → ∞ the method reduces to QM/MM-SS-Pol, while in the limit β → 0 the method reduces to QM/MM-SA-Pol.

The gradients of EMM‑Pol and EQM in eq 5 correspond to electric fields at the Drude particle site from the MM-Pol and QM systems, respectively; the former come trivially from eq 4, and the latter are available from analytical gradient theory and in general require solving coupled-perturbed response equations, which are available at many levels of electronic structure theory, including CIS, TDDFT, EOM-CCSD, CASSCF, and XMSCASPT2. Since both the QM and MM-Pol electric fields depend on the location of all Drude particles, eq 5 is solved self-consistently, typically with a dual-SCF procedure, whereby Drude updates are alternated with electronic structure calculations. For electronic structures where EQM is variational in all wave function parameters, the QM electric fields come exactly from the Hellman−Feynman theorem (i.e., a simple expectation value of the electron−nuclear electrostatic field), and then computational savings can be realized by updating the Drude particles during the QM SCF cycle using the updated one-electron density at each iteration of the SCF.28 An alternative to eq 5, which we use in this paper, is to numerically minimize eq 3 with respect to the Drude particle positions, which ensures that the variational condition is met. 2.2. Dynamically Weighted Polarizable QM/MM. Unfortunately, as mentioned in the Introduction, QM/MMSS-Pol suffers from a number of problems: (1) The induced dipole moments resulting from the optimization procedure differ for each electronic state. This means that each electronic state is an eigenstate of a different Hamiltonian, and one should not expect electronic state crossings to be described correctly. (2) In the vicinity of state crossings there may exist multiple solutions to the induced MM-Pol equations, resulting in artifacts in the PES and/or poor convergence in the QM/MMSS-Pol optimization. To overcome these problems, we have developed a new approach that we call dynamically weighted polarizable QM/MM, or QM/MM-DW-Pol for short. The approach borrows from our recently developed dynamically weighted CASSCF method29 and modifies the objective function (eq 3) to include a weighted average of the individual QM state energies, EQM I :

3. RESULTS AND DISCUSSION To test the performance of QM/MM-Pol methods, we applied them to a simple system involving a single lithium fluoride molecule interacting with a single beryllium atom. This system was chosen for a number of reasons: (1) LiF is a widely studied model system involving a CT state crossing with neutral states; (2) Be is highly polarizable31 and therefore captures a simple polarizable subsystem with a single atom; (3) the system contains few electrons, allowing the use of the MRCISD level of theory; (4) the closed shell of Be means that covalent interactions between the Be and LiF are negligible; and (5) the small number of electrons of Li, F, and Be together with a chosen large separation between Be and LiF (4 Å) ensures minimal dispersion interactions compared with the dominant polarization effects in which we are primarily interested. The geometry of LiF + Be was chosen to be linear, with a variable Li−F bond length denoted as R1. The Be atom was placed on the F side at a distance R2 = 4 Å from the F atom. In sections I and II in the Supporting Information, we show that this arrangement does indeed lead to minimal dispersion and covalent interactions compared with the dominant induced polarization interaction.

Nstate

∑ wIEIQM + EMM + EQM/MM + EMM‐ Pol I

+ Eself

N

QM ∑ J state f (EJQM − E REF )

Since eq 7 and therefore also eq 6 depend on the choice of reference state, when constructing the PES of a system, one must select each state as the reference state, one by one, in a series of independent calculations, just as in the QM/MM-SSPol method. For each reference state, the MM dipoles are then optimized to minimize eq 6, which leads to the following modification to eq 5:

where the inequality in the second sum indicates that electrostatic interactions between charge sites and their attached Drude particles are ignored. In QM/MM-SS-Pol, the induced dipoles are found by minimizing Etot with respect to the Drude particle positions, rα′. Applying this condition to eq 3 results in the following set of coupled equations for the Drude particle positions: α rα ′ = rα − α2′ ∇rα′(EMM ‐ Pol + EQM) qα ′ (5)

Etot‐DW =

QM f (EIQM − E REF )

(6)

We choose the weights wI to be normalized functions of the QM energy gaps between the reference state (the state of interest) and the other states: C


Article


(MRCISD/SA-4-CASSCF(6,7)/6-31g,F=6-31+g*) level of theory. Here it is seen that to the left of the avoided crossing (R1 < 6.25 Å), S0 has CT character and its energy in LiFBe is below that of LiF, and to the right of the avoided crossing (R1 > 6.25 Å), S3 takes CT character and its energy in LiFBe is below that of LiF. As a result of the CT stabilization, the avoided crossing in LiFBe moves to R1 = 6.4 Å and the gap at the crossing closes to 0.08 eV. 3.2. QM/MM Results. The need for polarizable QM/MM has been debated now for over two decades, so before testing QM/MM-Pol methods, we first apply standard QM/MM on LiFBe, treating LiF as the QM subsystem and Be as the MM subsystem. Since Be is a neutral atom, the resulting nonpolarizable QM/MM model contains no electrostatic interactions between the MM and QM subsystems, and thus, the only coupling term between these subsystems is EQM/MM, the short-range repulsion and dispersion interactions, often called the mechanical coupling. At a Be distance of R2 = 4 Å, both short-range repulsion and dispersion interactions are negligible, as confirmed in Figure S1, and we set EQM/MM = 0. Thus, the QM/MM results, plotted in Figure 2a, are the same as the isolated LiF results, since there is effectively no interaction between the neutral Be and LiF at the QM/MM level. This leads to poor agreement with the full QM results due to the complete neglect of Be polarization that serves to stabilize the CT state (S0 to the left of the crossing and S3 to the right). This highlights the importance of using a polarizable QM/MM approach for this system. 3.3. QM/MM-Pol Results. To test the QM/MM-Pol schemes discussed in section 2, we now treat LiFBe with QM/MM-Pol, where LiF is the QM subsystem and the Be atom is the MM-Pol subsystem. We modeled Be as a classical Drude oscillator25 (also known as a charge on a spring) with a Drude charge of qα′ = −4e and an isotropic polarizability of 29.86 bohr3, chosen to reproduce the polarizability of Be in LiFBe from a finite field and distributed multipole analysis35 calculation (see the Supporting Information for more details). For a given LiFBe geometry, the induced dipole on Be was found by minimizing the QM/MM-Pol energy (eq 3 or eq 6) with respect to the Drude charge position using the Broyden− Fletcher−Goldfarb−Shanno (BFGS) algorithm36 and centraldifference gradients (with a step of Δrα′ = 0.01 Å) as implemented in CIOpt.36 To avoid getting trapped in a highenergy local minimum, we used two initial guesses in the optimization of the Be Drude position, one corresponding to a zero dipole and the second as the expected dipole moment for a purely ionic Li+F− configuration. Following the variational principle, we selected the lowest-energy solution from the two optimizations. In the following, we test three different QM/MM-Pol schemes, the state-specific (SS) scheme discussed in section 2.1, a state-averaged (SA) scheme that amounts to setting all of the weights to be equal in eq 6, and the dynamically weighted (DW) scheme discussed in section 2.2. 3.3.1. QM/MM-SS-Pol Results. In Figure 2a, we plot the S0 (black curve) and S3 (red curve) PESs of LiFBe as functions of the Li−F distance (R1) for a fixed Be distance (R2 = 4 Å) near the avoided crossing region using the QM/MM-SS-Pol method to compare with the reference full QM calculations (circle symbols). In the SS scheme, the MM-Pol dipoles are optimized to a particular electronic state, i.e., EQM is the QM energy of either S0 or S3. Although the agreement between QM/MM-SSPol and the reference full QM calculations appears to be very

In all of the calculations below we used the MRCISD/SA-4CASSCF(6,7)/6-31g,F=6-31+g* level of theory without symmetry in Molpro 2015.1 update 11. For computational efficiency, polarization and diffuse functions were included only on the F atom (the 6-31g basis for beryllium includes p functions, allowing it to polarize). At dissociation, the (6,7) active space comprises the Li(2s), F(2p), and F(3p) atomic orbitals, which has been shown previously to give a good description of LiF’s electronic structure.32 Figure 1a shows the

Figure 1. Potential energy surfaces of LiF and LiFBe at the MRCI/SA4-CASSCF(6,7)/6-31g,F=6-31+g* level of theory. (a) PESs of LiF (S0, dashed black curve; S1, dot-dashed blue curve; S2, dotted green curve; S3, dashed red curve), indicating the avoided crossing between the S1 and S3 states, which dissociate to the Li0F0 and Li+F− limits, respectively. (b) Zoomed view of the avoided crossing region. The stabilizing effect of including a single Be atom on states with CT character is shown by plotting the PESs of LiFBe (solid curves) with the Be atom fixed at R2 = 4 Å from the fluorine atom in a linear arrangement.

calculated potential energy surfaces (PESs) of the four lowest states of the isolated LiF system, displaying an avoided crossing between the 1Σ+g ground state (S0) and the 2Σ+g excited state (S3) at R1 = 6.25 Å with an energy gap of 0.1 eV.33 The 1Πu (S1, S2) states are shown for reference in Figure 1a but have been removed from all subsequent figures for clarity. In this and all of the subsequent figures, the zero of energy is taken to be the ground-state energy of the system with all atoms at dissociation.34 3.1. Reference Full QM Results. Addition of a Be atom at R2 = 4 Å from the LiF molecule results in an expected stabilization of states with CT character relative to the neutral states. This is shown in Figure 1b, which compares the avoided crossing regions for LiF and LiFBe, both treated at the full QM D


Article


moment is plotted for each state, as shown in Figure 2b. It can be seen that the dipole moment of S0 (solid black curve) changes discontinuously at R1 = 6.4 Å, which can be understood as the character of S0 changing discontinuously from CT at R1 < 6.4 Å to neutral at R1 > 6.4 Å. The presence of a derivative discontinuity in the PES is unphysical since S0 and S3 should exhibit an avoided crossing rather than a conical intersection because of the linear geometry of LiFBe. Furthermore, the presence of discontinuities in the forces makes QM/MM-SS-Pol unsuitable for ab initio molecular dynamics. The issue of the discontinuous change in character of the electronic states in QM/MM-SS-Pol is explored in more detail in Figure 2c, which plots the energy of S0 as a function of the induced Be dipole for several Li−F distances. Here it is apparent that near the avoided crossing between S0 and S3, the QM/MM-SS-Pol energy exhibits two minima: one at μ ≈ −0.025 e Å with neutral S0 character and one at μ ≈ −0.175 e Å with CT S0 character. Since in the QM/MM-SS-Pol method the MM dipoles are optimized such that they minimize the QM/MM-SS-Pol energy (eq 3), the presence of multiple minima in the objective function means that there are multiple solutions for the MM dipoles. This explains why root flipping was observed using the traditional dipole update procedure (eq 5). Furthermore, as a function of increasing Li−F distance, the energy difference between the minima in Figure 2c is seen to decrease until they become degenerate at the crossing point (R1 = 6.3925 Å), after which the minimum at μ ≈ −0.025 e Å becomes the variational solution. This behavior, akin to a firstorder phase transition, explains the derivative discontinuity in the S0 energy as a function of Li−F distance. It is important to note that this electronic discontinuity occurs in S0. This means that reference-state QM/MM-Pol methods, including QM/ MM-Fr-Pol and QM/MM-LR-Pol, will also display this discontinuity, ruling out their use for excited-state dynamics involving photodeactivation to the ground state. 3.3.2. QM/MM-SA-Pol Results. Others have noted the problem of root flipping in QM/MM-SS-Pol, and one suggestion to alleviate the problem is to use state averaging,22 which we denote as QM/MM-SA-Pol, such that the induced MM dipoles are found by minimizing the QM/MM-Pol energy averaged over the electronic states of interest. This can be realized in eq 6 by setting the weights, wI, to be equal. We explore this approach in Figure 3a, which plots the S0 and S3 PESs of LiFBe with QM/MM-SA-Pol (dashed curves) compared with the reference full QM results (circles). Here it can be seen that QM/MM-SA-Pol appears to solve the problem of multiple solutions, and both S0 and S3 have continuous derivatives, as also confirmed by their continuous dipole moments (dashed curves in Figure 3b). However, away from the S0−S3 crossing, the QM/MM-SA-Pol energies are seen to be higher than the full QM PES for states of CT character (S0 for R1 < 6 Å and S3 for R1 > 6.6 Å). This overestimation of the CT energy is a result of the optimization of the induced dipole on Be to the average of the electric fields of S0−S3. For R1 > 6.6 Å, S3 is the CT state and S0−S2 are neutral states, so because of the state averaging Be experiences an electric field with approximately only one-quarter the magnitude it should have for S3 and is therefore considerably underpolarized. This means that, in general, QM/MM-SA-Pol energetics will be inaccurate even away from electronic state crossings. Li and co-workers introduced a state-specific correction to approximately address this issue;22 however,

Figure 2. Performance of QM/MM and QM/MM-SS-Pol on LiFBe with Be treated as an MM atom. (a) PESs of S0 (black curve) and S3 (red curve) computed with QM/MM (dashed curves) and QM/MMSS-Pol (solid curves) compared with the reference full QM calculations from Figure 1b (black and red circles for S0 and S3, respectively). (b) Total dipole moment of the system for S0 and S3 computed with QM/MM-SS-Pol compared with full QM (same color scheme as in (a)), where it is seen that the dipole moment of S0 exhibits a discontinuity (indicated by a dotted line) at the crossing with S3 despite the use of increments in R1 as small as 0.0075 Å. (c) QM/MM-SS-Pol energy of S0 as a function of the induced dipole on Be for five Li−F distances evenly spaced between R1 = 6.3775 Å and R1 = 6.4075 Å that bracket the S0−S3 crossing. Here, multiple minima are seen, which explains why the QM/MM-SS-Pol procedure, which optimizes the Be dipole to minimize the S0 energy, results in the dipole discontinuity seen in (b).

good, a number of problems were noted. First, when we tried to use the traditional dipole update scheme of eq 5, we found that near the avoided crossing, the dipole optimization failed to converge because of root flipping (the character of the reference electronic state would flip from CT to neutral and back at each iteration). Only by using the BFGS algorithm could we get stable convergence. Second, the S0 PES appears to have a derivative discontinuity at the crossing point at R1 = 6.4 Å. This problem is more apparent when the total dipole E


Article


and state-averaged methods by optimizing the MM dipoles to a weighted average of electronic-state fields (eq 8) with weights that depend on the electronic energy gaps. Figure 3a,b shows that QM/MM-DW-Pol with β = 4 (eV)−1 gives very good agreement with the reference full QM calculations on LiFBe, both in the vicinity of and away from the S0−S3 crossing. Figure 3c,d illustrates how the method works: at the S0−S3 crossing, the weights of S0, S1, S2, and S3 in the objective energy function (eq 6) are equal, so the method reproduces the QM/MM-SA-Pol results, which also agree reasonably well with the reference full QM energies and dipoles. Away from the crossing, for R1 < 6 Å, the weight of S0 in the S0 reference calculation (Figure 3c, solid black curve) approaches unity, so the method reproduces the QM/MM-SSPol results, which are in good agreement with the reference full QM results. Likewise, for R1 > 6.5 Å, the weight of S3 in the S3 reference calculation (Figure 3d, dotted red curve) approaches unity, giving good agreement with the full QM S3 PES. It should be noted that the weight of S0 in the S0 reference calculation does not approach unity for R1 > 6.5 Å. This is desired behavior, since S0, S1, and S2 are degenerate for dissociated Li−F and therefore have equal weights in QM/ MM-DW-Pol. Finally, the effect of averaging over state energies in the objective function is to remove the presence of multiple minima, as illustrated in Figure 3e, which shows that eq 8 has only a single solution for each Li−F distance.

4. CONCLUSIONS In this work, we explored the accuracy of several excited-state polarizable QM/MM methods on a simple model system of lithium fluoride interacting with a polarizable beryllium atom. We found that none of the previous methods performed well over the entire PES: either they displayed artifacts in the vicinity of the electronic state crossing (QM/MM-SS-Pol) or they had inaccurate excitation energies away from the crossing (QM/MM-SA-Pol). This motivated us to develop a new method, called dynamically weighted polarizable QM/MM, that is a hybrid of QM/MM-SS-Pol and QM/MM-SA-Pol. By treating the interaction of the QM and MM-Pol subsystems for different QM states in an equal fashion at state crossings, the qualitative features of the PES near the state crossing are captured. When the electronic states are energetically separated, the method reduces to state-specific QM/MM Pol, and the PES is in quantitative agreement with high-level QM calculations on the entire system. Although we have demonstrated our method on a relatively simple QM system at the MRCISD level and a single polarizable MM atom, it is trivial to generalize the approach to other electronic structure methods such as TDDFT and to extend the MM region to many coupled polarizable particles. Given the similarities of the objective function underlying the QM/MM-DW-Pol method to that underlying traditional QM/ MM-SS-Pol, the scaling of the method with system size will also be the same, allowing us to perform ab initio QM/MM-DWPol excited-state dynamics of complex systems. In addition, given the implementational similarities between polarizable QM/MM and polarizable continuum models (PCMs), it is expected that dynamic weighting can also be applied to excitedstate PCM calculations. Work in these directions is underway.

Figure 3. Performance of QM/MM-SA-Pol and QM/MM-DW-Pol on LiFBe with Be treated as an MM-Pol atom. (a, b) Same as Figure 2a,b, but comparing QM/MM-SA-Pol (dashed curves) and QM/MM-DWPol (solid curves) to full QM results (circles). (c, d) Weights for QM/ MM-DW-Pol for S0 and S3 reference state calculations, respectively. (e) Same as Figure 2c but plotting the objective energy function for QM/MM-DW-Pol (eq 6) as a function of induced dipole on Be for five Li−F distances evenly spaced between R1 = 6.2 and 6.6 Å that bracket the S0−S3 crossing. Here only a single minimum is seen for each Li−F distance, explaining why QM/MM-DW-Pol yields smooth PESs.

given the issues with state-specific methods already discussed, we seek an alternative approach. 3.3.3. QM/MM-DW-Pol Results. The fact that QM/MM-SSPol provides a good description of PESs away from electronic state crossings but performs poorly in the vicinity of crossings, while the reverse is true for QM/MM-SA-Pol, suggests that a hybrid approach that interpolates between these two methods based on the electronic energy gaps should be fruitful. We thus developed the QM/MM-DW-Pol method described in section 2.2, which achieves an interpolation between the state-specific F


Article


■

(12) Koch, D. M.; Peslherbe, G. H. Importance of Polarization in Quantum Mechanics/Molecular Mechanics Descriptions of Electronic Excited States: NaI(H2O)n Photodissociation Dynamics as a Case Study. J. Phys. Chem. B 2008, 112, 636−649. (13) Levine, B. G.; Martínez, T. J. Isomerization Through Conical Intersections. Annu. Rev. Phys. Chem. 2007, 58, 613−634. (14) Arora, P.; Slipchenko, L. V.; Webb, S. P.; DeFusco, A.; Gordon, M. S. Solvent-Induced Frequency Shifts: Configuration Interaction Singles Combined with the Effective Fragment Potential Method. J. Phys. Chem. A 2010, 114, 6742−6750. (15) Thompson, M. A. QM/MMpol: A Consistent Model for Solute/Solvent Polarization. Application to the Aqueous Solvation and Spectroscopy of Formaldehyde, Acetaldehyde, and Acetone. J. Phys. Chem. 1996, 100, 14492−14507. (16) Kongsted, J.; Osted, A.; Mikkelsen, K. V.; Christiansen, O. Linear response functions for coupled cluster/molecular mechanics including polarization interactions. J. Chem. Phys. 2003, 118, 1620− 1633. (17) Curutchet, C.; Muñoz-Losa, A.; Monti, S.; Kongsted, J.; Scholes, G. D.; Mennucci, B. Electronic Energy Transfer in Condensed Phase Studied by a Polarizable QM/MM Model. J. Chem. Theory Comput. 2009, 5, 1838−1848. (18) Glover, W. J. Communication: Smoothing out excited-state dynamics: Analytical gradients for dynamically weighted complete active space self-consistent field. J. Chem. Phys. 2014, 141, 171102. (19) Improta, R.; Barone, V.; Scalmani, G.; Frisch, M. J. A statespecific polarizable continuum model time dependent density functional theory method for excited state calculations in solution. J. Chem. Phys. 2006, 125, 054103. (20) Zeng, Q.; Liang, W. Analytic energy gradient of excited electronic state within TDDFT/MMpol framework: Benchmark tests and parallel implementation. J. Chem. Phys. 2015, 143, 134104. (21) Daday, C.; Curutchet, C.; Sinicropi, A.; Mennucci, B.; Filippi, C. Chromophore−Protein Coupling beyond Nonpolarizable Models: Understanding Absorption in Green Fluorescent Protein. J. Chem. Theory Comput. 2015, 11, 4825−4839. (22) Li, Q.; Mennucci, B.; Robb, M. A.; Blancafort, L.; Curutchet, C. Polarizable QM/MM Multiconfiguration Self-Consistent Field Approach with State-Specific Corrections: Environment Effects on Cytosine Absorption Spectrum. J. Chem. Theory Comput. 2015, 11, 1674−1682. (23) DeFusco, A.; Minezawa, N.; Slipchenko, L. V.; Zahariev, F.; Gordon, M. S. Modeling Solvent Effects on Electronic Excited States. J. Phys. Chem. Lett. 2011, 2, 2184−2192. (24) Jacobson, L. D.; Herbert, J. M. A Simple Algorithm for Determining Orthogonal, Self-Consistent Excited-State Wave Functions for a State-Specific Hamiltonian: Application to the Optical Spectrum of the Aqueous Electron. J. Chem. Theory Comput. 2011, 7, 2085−2093. (25) Drude, P.; Mann, C.; Millikan, R. The Theory of Optics; Longmans, Green: New York, 1902. (26) Lamoureux, G.; Roux, B. Modeling induced polarization with classical Drude oscillators: Theory and molecular dynamics simulation algorithm. J. Chem. Phys. 2003, 119, 3025−3039. (27) Wilson, M.; Madden, P. A.; Costa-Cabral, B. J. Quadrupole Polarization in Simulations of Ionic Systems: Application to AgCl. J. Phys. Chem. 1996, 100, 1227−1237. (28) Lu, Z.; Zhang, Y. Interfacing ab initio Quantum Mechanical Method with Classical Drude Osillator Polarizable Model for Molecular Dynamics Simulation of Chemical Reactions. J. Chem. Theory Comput. 2008, 4, 1237−48. (29) Chibani, S.; Le Guennic, B.; Charaf-Eddin, A.; Maury, O.; Andraud, C.; Jacquemin, D. On the Computation of Adiabatic Energies in Aza-Boron-Dipyrromethene Dyes. J. Chem. Theory Comput. 2012, 8, 3303−3313. (30) Deskevich, M. P.; Nesbitt, D. J.; Werner, H.-J. Dynamically weighted multiconfiguration self-consistent field: Multistate calculations for F + H2O → HF + OH reaction paths. J. Chem. Phys. 2004, 120, 7281−7289.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b00064. Estimate of the magnitudes of dispersion and induced polarization interactions in LiFBe, potential energy surfaces of LiFBe as a function of Be distance, derivation of QM/MM-DW-Pol analytical gradients, and sensitivity of QM/MM-DW-Pol to the width parameter (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Muhammad A. Hagras: 0000-0002-6062-2544 William J. Glover: 0000-0002-2908-5680 Funding

This work was supported by an NYU Global Seed Grant for Collaborative Research. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Profs. Yingkai Zhang and Todd Martinez for helpful comments. W.J.G. thanks NYU Shanghai for startup funds.

■

REFERENCES

(1) Warshel, A. Bicycle-pedal model for the first step in the vision process. Nature 1976, 260, 679−683. (2) Ben-Nun, M.; Molnar, F.; Lu, H.; Phillips, J. C.; Martínez, T. J.; Schulten, K. Quantum dynamics of the femtosecond photoisomerization of retinal in bacteriorhodopsin. Faraday Discuss. 1998, 110, 447− 462. (3) Kukura, P.; McCamant, D. W.; Yoon, S.; Wandschneider, D. B.; Mathies, R. A. Structural observation of the primary isomerization in vision with femtosecond-stimulated Raman. Science 2005, 310, 1006− 1009. (4) The Nobel Prize in Chemistry 2016 - Advanced Information. http://www.nobelprize.org/nobel_prizes/chemistry/laureates/2016/ advanced.html (accessed Jan 21, 2018). (5) Martínez, T. J. Insights for Light-Driven Molecular Devices from Ab Initio Multiple Spawning Excited-State Dynamics of Organic and Biological Chromophores. Acc. Chem. Res. 2006, 39, 119−126. (6) Böckmann, M.; Doltsinis, N. L.; Marx, D. Enhanced photoswitching of bridged azobenzene studied by nonadiabatic ab initio simulation. J. Chem. Phys. 2012, 137, 22A505. (7) Tsien, R. Y. The Green Fluorescent Protein. Annu. Rev. Biochem. 1998, 67, 509−544. (8) Olsen, S.; Lamothe, K.; Martínez, T. J. Protonic Gating of Excited-State Twisting and Charge Localization in GFP Chromophores: A Mechanistic Hypothesis for Reversible Photoswitching. J. Am. Chem. Soc. 2010, 132, 1192−1193. (9) Park, J. W.; Rhee, Y. M. Electric Field Keeps Chromophore Planar and Produces High Yield Fluorescence in Green Fluorescent Protein. J. Am. Chem. Soc. 2016, 138, 13619−13629. (10) Warshel, A.; Levitt, M. Theoretical studies of enzymic reactions: Dielectric, electrostatic and steric stabilization of the carbonium ion in the reaction of lysozyme. J. Mol. Biol. 1976, 103, 227−249. (11) Isborn, C. M.; Götz, A. W.; Clark, M. A.; Walker, R. C.; Martínez, T. J. Electronic Absorption Spectra from MM and ab Initio QM/MM Molecular Dynamics: Environmental Effects on the Absorption Spectrum of Photoactive Yellow Protein. J. Chem. Theory Comput. 2012, 8, 5092−5106. G


Article

Journal of Chemical Theory and Computation (31) Stiehler, J.; Hinze, J. Calculation of static polarizabilities and hyperpolarizabilities for the atoms He through Kr with a numerical RHF method. J. Phys. B: At., Mol. Opt. Phys. 1995, 28, 4055. (32) Bauschlicher, C. W.; Langhoff, S. R. Full configurationinteraction study of the ionic−neutral curve crossing in LiF. J. Chem. Phys. 1988, 89, 4246−4254. (33) Our calculated avoided crossing distance of 6.3 Å is somewhat shorter than the established benchmark value of 6.7 Å, mainly because of our use of a relatively small basis set. (34) Since MRCISD is not size-consistent, for systems involving just the LiF subsystem treated at the QM level we also included a QM Be atom at R2 = 100 Å. This allowed a direct comparison to the full QM LiFBe results by correcting for non-size-consistency. (35) Stone, A. J. Distributed multipole analysis, or how to describe a molecular charge distribution. Chem. Phys. Lett. 1981, 83, 233−239. (36) Fletcher, R. Practical Methods of Optimization, 2nd ed.; John Wiley & Sons: New York, 1987.

H


Polarizable Embedding for Excited-State ... - ACS Publications

Recommend Documents