Point charge embedding for ONIOM excited states calculations

The accurate description of electronic excitations in large molecular systems is crucial in biochemical, energy, and materials sciences.^1,2 Even if accurate quantum-mechanical algorithms have become more efficient over the past several years, their application to large molecular structures is still prohibitive. Luckily, this is not always necessary; since in many cases excitations are localized on a particular chromophoric region, which can be partitioned in smaller subunits or layers that can be treated at different levels of theory. Such differential treatment is achieved with hybrid methods, where the core (or model) system, i.e., the region where the excitation is localized, is treated at a high level of theory while the rest is treated at a lower level. The careful combination of both partitioning and levels of theory can indeed provide excitation energies and transition properties of accuracy comparable with that of a full calculation at a high level of theory.

Initially, hybrid schemes combined a quantum mechanics method with a molecular mechanics force field, QM:MM.^3–5 However, nowadays hybrid schemes that combine two QM methods, QM:QM, are also widespread. The QM methods involved can be based on either wave function theory (WFT) or density functional theory (DFT). Among the QM:QM hybrid methods that directly treat the interaction energy between layers, DFT-in-DFT, WFT-in-DFT, and density matrix embedding are very popular, see Refs. 6–11 for an incomplete list. In this work, we use an alternative hybrid approach called ONIOM (Our own N-layered Integrated molecular Orbital molecular Mechanics)^12–18 developed originally by Morokuma and co-workers. Unlike other hybrid approaches, ONIOM is formulated as an extrapolation where the energy of the entire structure computed at a low level of theory is corrected by the difference between calculations at a high and low level of theory on the model system. In the case of a two-layer partitioning, the extrapolated energy is

where RL, MH, and ML refer to real-low, model-high and model-low calculations, respectively, and real refers to the entire system. The ONIOM approach has several advantages: it is very flexible in the choice of methods, and it only requires three truly independent calculations if no electronic embedding is introduced, so that no particular implementation is required. In case of severed bonds at the boundary between the model system and the rest, hydrogen link atoms are usually added for capping.

The ONIOM extrapolation can be also formulated for vertical transition energies^19–31 with QM or MM methods as low level. One of us benchmarked ONIOM(QM:QM) for electronic transition energies and properties against a high level of theory such as the equation of motion coupled cluster singles and doubles (EOM-CCSD) method.^31–33 These studies showed that considerable time savings can be achieved with moderate loss in accuracy. They also showed that the use of link atoms does not affect the accuracy as long as the partitioning is sensible. However, despite these promising results, a bare ONIOM calculation (usually referred as mechanical embedding, ONIOM-ME, since the effect of the environment is only felt on the geometry) neglects the polarization effects of the surrounding on the electron density of the high-level calculation, which may be very important in large structures. Such polarization can be included by embedding the model system with point charges located at the position of the nuclei of the external layer. This approach is referred to as electronic embedding, ONIOM-EE. An extensive investigation of ONIOM-EE has been performed by the Raghavachari group for ground state problems, comparing different types of point charges, e.g., Mulliken and Lowdin,^34,35 and density based embedding.³⁶ They also suggested an approach to balance the charge transfer across the QM:QM boundary.^37,38 One issue with fixed charge embedding is that an unphysical overpolarization of the electronic density may occur at the boundary region. Several techniques have been proposed to alleviate this issue, including the deletion of one-electron integrals, the shifting/redistribution of point charges, or the replacement of point charges with Gaussian charge distributions.^39,40

In this work, we study the performance of the point charge embedding for ONIOM excited state calculations. We compare Mulliken charges and charges fitted to reproduce the molecular electrostatic potential (ESP) of the real-low calculation. For the latter, we also explore various constrained fittings to improve the charge distribution around the link atoms, which is the critical region for overpolarization. The constrained fitting is performed on the atomic point charges of the outer layer in the presence of the fixed charges of both the model system and link atom(s) obtained from the unconstrained fitting. Our results show that the overpolarization can be reduced, and the fitted charges consistently improve the performance of ONIOM-ME.

The paper is organized as follows: Section II reports the methods, while Section III reports the computational details; Section IV collects the results for the test cases we examined and discusses their convergence with respect to reference calculations; a general discussion and conclusions are reported in Section V.

In ONIOM-ME, the excitation energy $ω$ and the oscillator strength f can be computed with a similar partitioning of the structure as in the ground state, where the result of three separate calculations can be combined as^31,32

However, the partitioning should insure that the excitation of interest is mostly localized on the model region. This can be accomplished by a natural transition orbital (NTO)⁴¹ analysis at the low level of theory on the entire system (the RL calculation).³¹ The extrapolation is accurate when the differences with respect to the target, i.e., $Δω=ωONIOM−ωRH$ and $Δf=fONIOM−fRH$⁠, are small. $ωRH$ and f_RH refer to the real-high transition energy and oscillator strength, i.e., the target values obtained at the high level of theory considering the entire molecular structure. We will use $Δω$ and $Δf$ (whose values are obviously unknown in production calculations) as our metric to assess the accuracy of the various embedding schemes.

A point charge embedding may be also added to the ONIOM extrapolation in order to improve its accuracy, as long as overpolarization is minimized. In fact, overpolarization may be a larger issue for excited states than for the ground state due to the typically larger spatial extent of the excited state electron density. In this work, we test several embedding options based on Mulliken and ESP charges obtained from the real-low calculation. These charges are based on the ground state density and do not change in the excited state part of the calculation. In this respect, this approach is similar to other frozen embedding schemes.⁹ The ESP charges are fitted to minimize the root-mean-square of the difference between the potential generated by the point charges and the potential generated by the QM electron density and nuclei on a grid surrounding the molecule. Once the charges are computed, the interaction with the electrons is added to the one-electron Hamiltonian³⁴ 

where $μ$ and $ν$ are the atomic basis functions of the model-high or model-low calculations, q_A is the embedding charge centered at rA, and d is the damping factor for obtaining a scaled interaction. The interaction between the nuclei and the embedding charges is expressed by a classical Coulomb interaction, and it does not affect the excitation energy.

In order to reduce overpolarization, we propose two other options based on the ESP fitting, adjusted to account for the presence of the link atoms. We will refer to these charges as link atom modified electronic embedding: EE(L). Figures 1 and 2 show the flowcharts of the algorithms involved in the EE(L) fitting. Since the two approaches only differ in the value of the link atom charge(s) used in the fitting, we will focus on the first method, EE(L1), for a detailed explanation and report only the differences for the second one, EE(L2). The first step is to perform a real-low calculation to obtain the unconstrained ESP charges. These charges are separated in three sets: the charges on the atoms in the model system, {q^M}, the charges on the atoms that will be replaced by link atoms, {q^B}, and the charges on the atoms of the rest of the structure, {q^R}. Our goal is to obtain improved values for {q^R} (used in ONIOM as an embedding field) by accounting for the presence of the link atom(s). To this aim, the {q^R} set is optimized by a constrained ESP fitting, where the {q^M} set is kept fixed at the original values. The {q^B} and {q^R} sets change as follows:

1. A ground state calculation is performed for the model system at high and low levels (MH and ML, respectively) in the presence of the embedding charges {q^R}.

2. An ESP fitting is performed to obtain the corresponding new link atom charges: {q^LinkH} and {q^LinkL}, for MH and ML, respectively.

3. A new {q^B} set is then computed from the average of {q^LinkH} and {q^LinkL}.

4. Updated {q^R} charges are obtained with a ESP fitting in the presence of the fixed {q^M} and new {q^B} sets of charges.

5. A new cycle is performed until convergence on the {q^R} set is achieved.

This iterative procedure does not significantly affect the overall cost of the calculation as the bottleneck is still represented by the excited state calculations, which are the same as with no embedding. The link atom positions are the same as that of the atoms being replaced (i.e., no rescaling of the bond length) in order to simplify the iterative fitting procedure. This is not an issue since we showed that the excitation energies are not very sensitive to the link atom bond length.³³ 

The EE(L1) scheme provides the same set of embedding charges for both model system excited state calculations. In the second scheme, EE(L2), the procedure is basically the same except that two different sets of embedding charges are created for the MH and ML calculations. This is accomplished by using the {q^LinkH} and the {q^LinkL} sets separately as {q^B} charges, see Figure 2. The computational cost of the EE(L1) and EE(L2) schemes is similar since the effort for the actual fitting is negligible, and the two schemes converge in a similar number of iterations, 3-5 on average. Once the {q^R} charges are obtained, the ML and MH excited state calculations can be performed in the presence of the embedding field with any standard quantum chemistry package without modification. The RL, MH, and ML excitation energies can then be combined as in Eqs. (2) and (3) to obtain the ONIOM extrapolated transition energy and oscillator strength. The {q^R} charges should be more balanced compared to regular ESP and Mulliken sets because they are adjusted to account for the presence of the link atom(s) at the boundary between layers. This approach does not require any rescaling or shift of the charges in proximity of the link atoms nor at the position of the link atom(s) themselves. For comparison with previous work,⁴² we also consider a third set of {q^R} charges obtained by a constrained ESP fitting where the link atom charges are set to zero. We will refer to these embeddings as EE(L0).

The CAM-B3LYP hybrid functional⁴³ with the 6-311+G(d,p) basis set was used for both target and model-high calculations. The target calculations are the reference for the comparison with the different extrapolations as they include the entire structure. The CAM-B3LYP functional was chosen for its ability to provide a balanced description of various types of electronic excitations thanks to its range separation.^44–46 Although it would be desirable to use more accurate methods than CAM-B3LYP as reference, the size of the systems we consider would make the target calculations prohibitive. Thus, CAM-B3LYP represents a good compromise between accuracy and cost for the purposes of this work. All vertical transition energies are computed using the usual linear response approach.^47–49 The layer definition follows the standard ONIOM approach, using hydrogen link atoms to cap dangling bonds. Two low-level methods are tested: B3LYP⁵⁰ and Hartree-Fock (HF) with the 6-31+G basis set. Although these two levels of theory are comparable in cost to CAM-B3LYP (if the same basis set is used), we are not interested here in the computational savings. What matters for this study is that the low levels of theory are less accurate than the target.⁴⁴ Savings in computational cost would be large with a wave function method based, for instance, on the coupled cluster theory as high level.³¹ All calculations were performed with a development version of the GAUSSIAN suite of programs.⁵¹ 

We consider six test compounds, shown in Figure 3: azonaphtharylamide (APA), taipinisine (TAI), flavin mononucleotide (FMN), betaine 30 (B30), and the retinal chromophore (RET and RET(AH)). For the latter, we consider two isomers with different protonations: on the Lys residue (RET) and on the Asp residue (RET(AH)). The test compounds were chosen in order to consider chromophores with localized excitations in different molecular and biomolecular environments: a pigment (APA), a natural product (TAI), a highly conjugated probe for solvent polarity (B30), the active site of both the bacterial photoreceptor (FMN) and the bovine retinal photoreceptor (RET). The size of the model systems and the number of link atoms are different across the test set so that we can explore the performance of the extrapolation in different conditions. The figure also shows the model system (ball-and-stick representation) and the rest (tube-frame representation). The geometries for all four molecules as well as all of the excitation energies are reported in the supplementary material. We focus on the first singlet excitation for all compounds, $S0→S1$⁠, except for RET for which we consider the second state, $S0→S2$⁠, as the first excitation is dark. The model system is chosen such that the excitation is mainly localized on this part of the structure. This is confirmed by the NTO⁴¹ analysis shown in Figure 4.

The vertical transition energies are shown in Figure 5, while the oscillator strengths are shown in Figure 6. We also considered scaling the charges (d = 0.8 in Eq. (4)) to test whether further correction for overpolarization was necessary. This is indeed not the case, so we will not discuss these results any further in the main text and report these data in Figures S1 and S2 of the supplementary material. All of the values of $ω$ and f are reported in Tables S6-S12 of the supplementary material.

We will start by discussing the results for the excitation energies for all molecules, followed by a discussion of the results for the oscillator strengths. The first test compound, APA, belongs to an important class of industrial azo pigments.⁵² The naphthalene ring is the central scaffold of the pigment, with the phenyl ring attached via the hydrazone group being less twisted and more conjugated than those attached via the carboxamide groups. The geometry was optimized at the high level of theory. For APA, the target excitation energy is at 3.08 eV. The ONIOM-ME results with B3LYP as low level are already in good agreement with the target, which is overestimated by 0.04 eV. Introducing the embedding with the Mulliken and ESP charges considerably worsens the ONIOM performance, with an error of 0.2 and 0.3 eV, respectively. Introducing the EE(L) charges, the ONIOM results further improve compared to the already good ONIOM-ME scheme. With HF as low level, ONIOM-ME overestimates the target by 0.06 eV. The embedding with the Mulliken and ESP charges overcompensates the correction, leading to an underestimation of the target by a slightly larger magnitude than with the ONIOM-ME scheme. The results with the EE(L) embedding are in excellent agreement with the target, especially with EE(L2).

The second test compound, TAI, is a recently isolated in-dole alkaloid from the extract of the flowering plant Malayan Tabernaemontana.⁵³ Its structural model is taken from the supplementary material of Ref. 53, and the target excitation energy is 4.80 eV. The B3LYP results without embedding underestimate the target by more than 0.2 eV. In this case, introducing the embedding with the Mulliken and ESP charges improves the ONIOM performance considerably, reducing the error to 0.1 and 0.02 eV, respectively. The error with the EE(L) schemes is larger than with the standard charges, and it is quite different between the two sets. $Δω<0$ with EE(L1) and $>0$ with EE(L2), and of the same order of magnitude than the ONIOM-ME error. Using HF as low level provides a smaller error for the ONIOM-ME scheme, $∼0.1$ eV. The error with the Mulliken charges is very large and negative, $|Δω|>0.5$ eV. With ESP, the target is still largely underestimated, with $|Δω|>0.5$ eV. The use of either EE(L1) or EE(L2) schemes brings the ONIOM results extremely close to the target. The EE(L0) results underestimate the target, with an error of the order of 0.1 eV in the case of HF and 0.2 eV in the case of B3LYP.

The third test compound, FMN, is the active site of a bacterial photoreceptor. This system was recently studied theoretically with QM:MM methods⁵⁴ due to its potential technological applications. The structure is obtained from experimental X-ray crystal data (PDB ID 2HFN)⁵⁵ considering 188 atoms for the target calculation. For FMN, the target excitation energy is 2.73 eV. With B3LYP as low level, all ONIOM results are in basically perfect agreement with the target independently of the choice of embedding. This is an important result as it shows no overpolarization effects. However, ONIOM-ME overestimates the target by 0.15 eV with HF as low level. Using Mulliken charges as electronic embedding improves the agreement with the target, $|Δω|≃0.1$ eV. With ESP, the target is underestimated by a similar amount as with the Mulliken charges. The agreement with the target is perfect with all EE(L) schemes and both low levels of theory.

The fourth test compound, RET, is the active site of a bovine retinal photoreceptor. This system was already investigated with QM/MM methods^27,56 because of the unique spectral properties of this chromophore. The retinal chromophore is covalently bound to a Lys residue interacting via hydrogen bonding with a nearby Asp residue. The structure is taken from the experimental X-Ray data (PDB ID 1U19), considering 310 atoms. For RET, the target excitation energy is 2.43 eV. The B3LYP results without embedding are in perfect agreement with the target, while HF underestimates the target by less than 0.1 eV. Introducing the embedding with the Mulliken and ESP charges worsens the agreement for B3LYP and slightly improves HF. The error with the EE(L) schemes is smaller than with the standard charges. We obtain a perfect agreement with EE(L1) and EE(L0), while a slight overestimation is obtained with EE(L2).

In order to further investigate the performance of the various embedding schemes, we moved the H⁺ from the Lys to the nearby Asp residue, changing the charge of the model system from +1 to 0 (RET(AH)). In this case, the target excitation energy is 2.79 eV. The B3LYP results without embedding underestimate the target by more than 0.1 eV. Introducing the embedding with the Mulliken or ESP charges worsens the agreement, and the results now overestimate the target by more than 0.2 eV. The error with the EE(L2) scheme is very small (0.02 eV), while with EE(L0) and EE(L1) the error is of similar magnitude as that obtained with ONIOM-ME. With HF as low level, ONIOM-ME and all EE(L) schemes show perfect agreement with the target, while Mulliken and ESP underestimate it by 0.1-0.15 eV.

The last test compound, B30, is a large conjugated molecule and a well known probe for solvent polarity because of the sensitivity of its photophysical properties to the nature of the solvent.^57,58 The geometry was optimized at the high level of theory. The target excitation is at 1.73 eV. B30 is challenging for several reasons: the excitation is charge transfer in nature, see Figure 4, thus problematic for a global functional like B3LYP; in the embedded calculations, the fixed point charges replace only a part of the phenyl substituents (one of the C centers is replaced by a link atom, and therefore it is not a part of the embedding), thus the description of the electrostatic effect of each aromatic group may be unbalanced. The results with B3LYP as low level show that ONIOM-ME underestimates the target by about 0.2 eV. Introducing the embedding does not improve the results as Mulliken and ESP increase the underestimation, while EE(L1) and EE(L2) overestimate the target by about 0.5 eV. Mulliken and ESP are close to each other, as are EE(L1) and EE(L2). The situation is mostly reversed with HF as the ONIOM-ME results largely overestimate the target by about 0.3 eV. The results with Mulliken and ESP overestimate the target even more than without embedding. On the other hand, the EE(L1) and EE(L2) schemes overcorrect in the opposite direction, now underestimating the target, although the magnitude of the error (⁠$∼0.2$ eV) is reduced compared to the ME case. Surprisingly, for this molecule, the use of the EE(L0) charges leads to a decisive improvement with respect to ONIOM-ME with both B3LYP and HF. In this case, setting the link atom charges to zero in the constrained fitting helps in balancing the description of the electrostatic field generated by the aromatic substituents.

The data for the oscillator strength, reported in Figure 6, show a trend that is similar to that of the excitation energy: the EE(L) schemes tend to improve upon the ONIOM-ME results (or leave them unaltered when the latter are already good), while the Mulliken and regular ESP embedding schemes have a more erratic behavior. For APA with B3LYP, ME overestimates the target while most embedding schemes underestimate it (the only exception is ESP); however, all embedding schemes are closer to the target than ONIOM-ME. With HF as low level, all schemes are in very good agreement with the target. For TAI and both low levels of theory, basically all schemes underestimate the target and behave approximately the same. The EE(L1) and EE(L2) schemes are in perfect agreement with the target in the HF calculations. For FMN and B3LYP, the ME scheme overestimates the target, while all embedding schemes are on target. With HF as low level, ME and Mulliken overestimate while ESP underestimates the target; all EE(L) schemes provide a nearly perfect performance. For RET, all methods overestimate the target by about the same amount (ESP is slightly better) with B3LYP and underestimate it with HF. For RET(AH), ONIOM-ME is basically on target with both low level methods, and EE(L1) performs almost as well. EE(L2) and EE(L0) perform very well with HF and slightly worse with B3LYP. On the other hand, Mulliken and ESP severely underestimate f with both low levels. Finally, all methods perform comparatively well for B30 with both low levels of theory.

In this work, we investigate the use of fixed point charges as electronic embedding in ONIOM excited state calculations. We compare several charge sets: Mulliken, ESP, EE(L1), EE(L2), and EE(L0). The last three sets are obtained as a constrained fitting of the ESP of the real-low calculation when the charge on the link atom(s) is taken into account (see Section II).

In order to discuss the overall performance of the various schemes, we consider averages of the relative error in the form of root mean square (RMS) and mean absolute error (MAE), shown in Figure 7. These data do not include the results for B30, which we will discuss separately. Both RMS and MAE show the same trends: embedding with Mulliken and ESP charges on average increases the error compared to no embedding for both the excitation energy and oscillator strength. The EE(L) constrained fitting schemes provide reduced errors compared to no embedding. For the excitation energy, this improvement is considerable with EE(L2) with both B3LYP and HF as low levels. For EE(L1) and EE(L0), the improvement is minimal with B3LYP and substantial for HF; in fact EE(L1) works better for HF than EE(L2). For the oscillator strength, EE(L2) performs on average slightly worse than no embedding with B3LYP, while EE(L1) and EE(L0) are marginally better. For HF, all EE(L) schemes outperform the ME scheme, and EE(L2) is the best among them. The B30 case is a cautionary tale that embedding cannot always be used to improve the ONIOM results, except for the surprising performance of the EE(L0) scheme. The reason is related to where the charges are positioned. The ONIOM partitioning is reasonable as only single C−C bonds are severed, and indeed ME results are rather good with DFT. However, the embedding charges try to reproduce the electronic distribution of aromatic rings while replacing only five of the six carbon centers in the ring (the sixth C center is replaced by a link atom). In this situation, all embedding schemes overpolarize the model system electron density; thus deteriorating the ONIOM performance.

These results show that the electronic embedding with fixed point charges can indeed improve the performance of ONIOM for excited states. The EE(L) schemes that we suggest seem to accomplish this in most circumstances, especially with two separate sets of charges, one for ML and one for MH: EE(L2). Nevertheless, we do not find a clear preferred choice for all cases as sometimes even Mulliken and ESP charges perform better than the EE(L) charges (e.g., for TAI with B3LYP). This is likely due to the lack of flexibility of ground state fixed point charges, which cannot reproduce any polarization response of the embedding field to the electronic excitation in the core region. Other terms such as exchange and repulsion can also affect the electronic excitation. We expect that including such terms will considerably improve the performance of ONIOM-EE, and we will present this development in future work. This future development will also include the implementation of energy gradients, which may probably follow a similar route as for the ground state.^34,35

See supplementary material for the geometries of test molecules, and all of the excitation energies.

Support from startup funds provided by the University of Kansas is gratefully acknowledged. J.B. also thanks the NSF Research Experiences for Undergraduates (REU) Grant No. 1263259.