Electronic Excitations in Complex Molecular Environments: Many-Body Green's Functions Theory in VOTCA-XTP.

Many-body Green's functions theory within the GW approximation and the Bethe-Salpeter Equation (BSE) is implemented in the open-source VOTCA-XTP software, aiming at the calculation of electronically excited states in complex molecular environments. Based on Gaussian-type atomic orbitals and making use of resolution of identity techniques, the code is designed specifically for nonperiodic systems. Application to a small molecule reference set successfully validates the methodology and its implementation for a variety of excitation types covering an energy range from 2 to 8 eV in single molecules. Further, embedding each GW-BSE calculation into an atomistically resolved surrounding, typically obtained from Molecular Dynamics, accounts for effects originating from local fields and polarization. Using aqueous DNA as a prototypical system, different levels of electrostatic coupling between the regions in this GW-BSE/MM setup are demonstrated. Particular attention is paid to charge-transfer (CT) excitations in adenine base pairs. It is found that their energy is extremely sensitive to the specific environment and to polarization effects. The calculated redshift of the CT excitation energy compared to a nucelobase dimer treated in vacuum is of the order of 1 eV, which matches expectations from experimental data. Predicted lowest CT energies are below that of a single nucleobase excitation, indicating the possibility of an initial (fast) decay of such an UV excited state into a binucleobase CT exciton. The results show that VOTCA-XTP's GW-BSE/MM is a powerful tool to study a wide range of types of electronic excitations in complex molecular environments.

Introduction

The functionality of many complex supramolecular assemblies is often defined not only by structural features and dynamics but also by their electronic excitations. Photosynthetic processes or microbial respiratory activity[1,2] are two prominent cases from nature in which complex processes emerge only from the interplay between electronic structure of molecular building blocks, such as a single chromophore, and nano- and mesoscale morphology.[3,4] Designing synthetic materials with similar functionality is attractive for many technological applications, such as in organic electronics,[5,6] thermoelectricity,[7] or in sensing and spectroscopy.[8]

Understanding and controlling this intimate interplay is crucial for a rational design of such materials, and computational studies hold an enormous potential in this context. Explicitly linking electronic and classical degrees of freedom is required, which on the relevant scales can only be achieved by coupled quantum-classical techniques. Approaches are needed that provide a high level of predictability for different types of excitations and that can be applied to realistic system sizes. Density functional theory[9,10] (DFT) is commonly used to determine ground state properties with reasonable accuracy, even for systems of considerable size. Reliable predictions in supramolecular systems are however often sensitive to the choice of an exchange-correlation (xc) functional or the addition of appropriate corrections for van der Waals effects. Interpretation of Kohn–Sham energies as single-particle excitation energies is subject to significant self-interaction errors and a lack of accounting for electronic screening effects as a response to the excitation.[11] Coupled electron–hole excitations, e.g., after photon absorption, are usually determined within the framework of time-dependent DFT (TD-DFT). Quality of results obtained with TD-DFT varies significantly with the excitation’s character. In particular the accurate description of excitons involving extended π-systems or ones with charge transfer character has been notoriously difficult.[11,12] Functionals which separate long-range and short-range interactions[13,14] can alleviate the problem but need to be adjusted to a specific system for good accuracy. The usage of accurate wave function based techniques, like CASPT2 (complete-active-space second-order perturbation theory)[15,16] or CC (coupled cluster) linear response theory using CC3 or CCSDT,[17,18] is often prohibitive due to their computational demands in application to large molecular systems.

Traditionally, many-body Green’s functions theory with the GW approximation and the Bethe-Salpeter Equation (BSE)[19] has been used within the solid state community. Recently, however it has gained attention for the treatment of electronically excited states of molecular systems.[20−23]GW-BSE has been shown to yield accurate descriptions of different types of excitations, such as localized (Frenkel) and bimolecular charge transfer (CT) excitons, on an equal footing.[21,22,24] Its computational cost is comparable to that of TD-DFT, which makes application to molecules or clusters of molecules of technological relevance tractable.[25−27] With the target of studying electronically excited states in complex molecular environments in mind, one of the current challenges is to link GW-BSE to atomistic models of large scale morphologies, thus treating these systems in hybrid quantum-classical (QM/MM) setups.

In this paper, we introduce our Gaussian orbital implementation of GW-BSE and its integration into a polarizable QM/MM workflow in the open-source VOTCA-XTP package. We first benchmark our implementation using the Thiel set[28,29] containing 28 small molecules of different types: unsaturated aliphatic hydrocarbons, aldehydes, ketones, amides, aromatic hydrocarbons, and heterocycles, as well as four nucleobases. The set offers reference data from both experiment and high-order wave function techniques for a representative variety of types of excitations, including π → π* (e.g., ethene), n → π* (e.g., cyclopropene), and σ → π* (e.g., pyrazine) excitations. These different excitations also cover a wide range of energies from 2 to 8 eV. After that, we illustrate the capabilities of the GW-BSE/MM setup by investigating a water-solvated DNA double-strand as a prototypical system.

Photophysical processes triggered by the absorption of ultraviolet solar radiation can cause damage and subsequent mutations in DNA.[30] Due to the complexity of such biological systems, an understanding of the processes involving excited states is extremely challenging. Bimolecular charge transfer excitations between base pairs are considered to play an important role in the excited-state dynamics in DNA and the mechanism by which these dynamics lead to structural or chemical decay and, eventually, gene mutations.[31] One of the proposed processes is the rapid decay of an initially photoexcited π → π* state to a longer-lived CT state, which induces either structural modifications or chemical oxidation/reduction reactions.[32]

As a first step to gain further insight into the exact conditions under which dynamical excited state processes of this kind can occur in DNA, a detailed understanding of the CT state energies is vital. While it is suggested from experiment on water-solvated single-strands of 20 adenine bases (A20) that CT states cause a faint UV absorption at energies below the energy of the UV active π → π* transition at approximately 5 eV,[33] a high-level second-order approximate coupled-cluster method yields CT excitations far above that for an isolated A2 dimer in the gas phase.[34] Due to the long-range electron–hole interaction, CT energies are typically very sensitive to the arrangement of the constituent monomers of the base pair. Optimized gas-phase structures are likely to exhibit different stacking distances and even motifs compared to a real single-strand and lack the effect of the aqueous environment not only on the structure but also on the electrostatic environment entirely.

Based on the experimental evidence, one should expect a redshift on the order of 1 eV for CT energies of DNA in an aqueous versus a vacuum environment. Inclusion of a polarizable environment using a Polarizable Continuum Model (PCM) on top of TD-DFT[35] fails to reproduce this observation, with the redshift being reported to be as small as 0.1 eV, even using functionals specifically containing long-range interactions. The use of model structures consisting of idealized molecular dimers and the lack of explicit local fields from a molecular environment comprising, e.g., a charged backbone, water solvent, and ions, and the quality of the functional used in the TD-DFT calculation to describe CT excitons are likely origins of this discrepancy. To reliably distinguish the different effects of the aqueous environment and to quantify how they affect the character of CT excitations contributing to the observed redshift, it is important to consider a realistic morphology of aqueous DNA and treat it with an accurate set of techniques. The GW-BSE formalism has been reported to yield very accurate predictions of CT excitation energies in prototypical small-molecule dimers. Subsequently, Yin et al.[36] studied small complexes of adenine dimers and water (A2-(H20)). Geometries of the complex were obtained from Classical Molecular Dynamics, while GW-BSE was used to evaluate the excitation energies. It was found that CT energies are strongly affected by the dipole electric fields in the first hydration shell around the A2, giving rise to an overall energetic shift to below that of the π → π* transition in single adenine, much more in line with the experimental observation.

Inspired by these results, we will use GW-BSE within the QM/MM framework discussed above on aqueous DNA. We go beyond the model used by Yin et al.,[36] and instead of a single hydrated adenine dimer, we consider a full double-strand of DNA solvated in explicit water. This will allow us to study, among other things, the effects of realistic stacking, the differences between intra- and interstrand charge transfer excitations, and the explicit effect of the DNA backbone. Given the sensitivity of the CT excitations to water hydration, inclusion of these structural parameters can have a substantial influence since the electrostatic environment will be vastly different from the idealized situation of a hydrated dimer. We will also be able to study dimers formed by different types of nucleobases on equal footing. We will further consider different embedding variants, including vacuum calculations as well as GW-BSE/MM calculations with static and additional polarizable interactions. Analyzing these different scenarios makes it possible to disentangle effects of geometric structure of the dimer, local electric fields or the structure of the environment, and electronic polarization.

This paper is organized as follows: In section , the theoretical concept of GW-BSE for the calculation of one- and two-particle excitations is briefly outlined, as well as the idea of coupling to a classical atomistic environment. Details of the implementation in VOTCA-XTP are given in section . The results for the Thiel set and the aqueous DNA system are discussed in section . A brief summary concludes the paper.

Electronically Excited States via GW-BSE: A Theoretical Framework

In the following section, major concepts behind GW-BSE are summarized. In order to keep the notation simple, we restrict the discussion to a spin singlet, closed shell system of N electrons. Its ground state, |N,0⟩, can be calculated using DFT by solving the Kohn–Sham (KS) equations[19]in which T0 is the kinetic energy, Vext is an external potential, VHartree is the Hartree potential, and Vxc is the exchange-correlation potential.

One-Particle Excitations

Particle-like excitations, quasiparticles (QP), in which an electron is added (N → N+1) to or removed (N → N–1) from |N,0⟩, are described by the one-body Green’s function[37,38]

In this notation, time and space variables are combined into a single variable (e.g. (r1,t1 ≡ 1)), T is the time ordering operator, and ψ and ψ† are the annhilation and the creation electron field operators, respectively.

G1 obeys a Dyson-type equation of motion, which in spectral representation iswhere H0 = T0 + Vext + VHartree is the DFT Hamiltonian in the Hartree approximation, while the exchange-correlation effects are described by the electron self-energy operator Σ(E). Knowing its exact form is crucial to solve the problem. It can be shown that eq is part of a closed set of coupled equations, known as Hedin equations.[39,40] DFT, for instance, can be seen as an approximated solution for the excited electrons problem, in which Σ ∼ Vxc. The related Green’s function G0, solution of eq , is .

An approximated solution beyond DFT of this system is given by the GW approximation, in whichwhere W = ϵ–1v is the screened Coulomb interaction, and v(r,r′) = 1/|r–r′| is the bare Coulomb interaction, respectively. The inverse dielectric function, ϵ–1, is calculated in the random-phase approximation (RPA).[41] Within this GW approximation, eq is transformed into a Dyson equation of motion for the quasiparticles[42,43]where εQP are the one-particle excitation energies of the system (i.e., the QP electron and holes states), and |ϕQP⟩ are the quasiparticle wave functions.

In practice, these quasiparticle wave functions are expanded in terms of the KS states according to . Assuming that |ϕQP⟩ ≈ |ϕKS⟩, the quasiparticle energies can be obtained perturbatively as

Both the correction term Δε and the nonlocal, energy-dependent microscopic dielectric function calculated within the RPA depend on εQP. Solutions to eq therefore in general need to be found self-consistently. It can be avoided by setting εQP = εKS in the evaluation of Σ and W, typically called the G0W0 approximation. To improve on this one shot approach, the G0W0 results can be used as the first step of an iterative evaluation of Σ(εQP), called GW0. In the evGW procedure, the quasiparticle energies are additionally updated in the RPA calculation, see also eq and Figure , until eigenvalue (ev) self-consistency.

Figure 1

GW-BSE workflow as implemented in VOTCA-XTP. The inner self-consistency loop corresponds to the GW0 algorithm, and the outer convergence loop, which requires the recalculation of the RPA, is the evGW.

The determination of εQP via eq typically holds if the off-diagonal elements of the self-energy, i.e., ⟨ϕKS|Σ(E)|ϕKS⟩, are small. Otherwise, expressing the QP wave functions as a linear combination of KS states needs to be fully taken into account. Quasiparticle wave functions and energies can then be obtained by diagonalizing the energy dependent QP Hamiltonian

The GW approach in which also the resulting quasiparticle wave functions, and not only the energies, are fed back into the RPA eq is also referred to as self-consistent scGW.

Two-Particle Excitations

Neutral excitations, in which the number of electrons remains the same but they assume an excited configuration S (|N,0⟩ → |N,S⟩), can be described based on the two-particle Green’s function.[38] It can be obtained solving a Dyson-like equation of motion, known as the Bethe-Salpeter Equation (BSE).[40] Defining the electron–hole correlation function aswhere the second term represents the independent movement of two particles (i.e electron and hole) as a product of single-particle Green’s functions and G2 as the two-particle Green’s function. The BSE readswith K(35,46) being the interaction kernel, and L0(12,1′2′) = G1(1,1′)G1(2,2′) is the two-particle noninteracting correlation function.

Under the assumption of optical excitations, which involve the simultaneous creation and annihilation of quasiparticles, we can reduce the four time variables to two. Due to time homogeneity L(12,1′2′) can be reduced to L(12,1′2′;ω), with the indices only representing position. The kernel K is given by the functional derivative of the full self-energy with respect to noninteracting quasiparticles. Using the GW approximation and assuming δW/δG1 ≈ 0, e.g., the screening is not influenced by the excitation, it can be shown thatK is called the direct interaction and originates from the screened interaction W between electron and hole and is responsible for the binding in the electron hole pair. K originates from the unscreened interaction ν and is responsible for the singlet–triplet splitting. It is denoted exchange interaction.

L0 can be written, assuming that G1 is fully given by electron and hole quasiparticles of the system, as a combination of independent excitations. In position space it readswhere v runs over all occupied (hole) states, and c runs over all unoccupied (electron) states.

Defining the electron–hole amplitude asallows eq to be rewritten aswhere S labels the two-particle excitation, and Ω is the corresponding excitation energy. To practically evaluate the BSE in eq , all quantities in eq and eq are expressed in terms of a basis set of single-particle electron and hole states. That is, introducingtransforms the BSE into an eigenvalue problem of the formor in the block matrix form:

The matrix elements Hres and K are given bywhere

Here it is assumed that the dynamic properties of W(ω) are negligible and the static approximation ω = 0 is used, which reduces the computational cost significantly. This is only valid if Ω – (ε–ε) ≪ ω, where ω is the plasmon frequency, which determines the screening properties.

For many systems the off-diagonal blocks K in eq are small and can be neglected. This leads to the Tamm-Dancoff approximation (TDA)[44]and the resulting electron–hole amplitude

This approximation halves the size of the BSE matrix. Additionally, it helps to reduce triplet instabilities,[45] but especially for small molecules the error from neglecting the coupling between resonant and antiresonant part can be significant.[20]

The spin structure of the BSE solutions depends on the spin–orbit coupling. If the ground state is a spin singlet state and spin–orbit coupling is small, the Hamiltonian decouples into singlet and triplet class solutions, with HsingletBSE = D + K + 2K and HtripletBSE = D + K. If spin–orbit coupling is large, the BSE Hamiltonian must be evaluated using the full spin structure. More complex spin contributions also arise for open shell systems, where the ground state is not a singlet.

GW-BSE/MM

Excitation energies in complex molecular environments can be obtained via a QM/MM procedure.[46−50] This method relies on treating the active subpart of the system quantum-mechanically (QM) while embedding it into an environment described at molecular mechanics resolution (MM). Recently, discrete static point charge models of a molecular environment have been introduced to plane-wave implementations of GW-BSE.[50,51] Such static classical models do not include environment polarization effects. Furthermore, the use of plane waves and the implied artificial periodicity for the study of isolated systems, such as molecules, is considered less efficient than the use of localized basis sets. Gaussian type orbital implementations of GW-BSE, on the other hand, have recently been coupled to continuum polarization models,[52] which lack explicit local electric fields from a complex molecular environment. Inclusion of polarization effects has further been reported for the GW formalism using polarizable point charge models.[53] In our QM/MM scheme, we employ a distributed atomic multipole representation, which allows for a general treatment of static electric fields and polarization effects on equal footing. The QM subpart can be treated, for instance, at the GW-BSE level, and molecules in the MM region are represented by static atomic multipole moments Q(54) where t indicates the multipole rank and a indicates the associated atom in the molecule. Additionally, each atom can be assigned a polarizability α which determines the induced moments ΔQ due to the field generated by moment t′ of atom a′. The classical total energy of a system in the state (s) (i.e. ground or excited) composed of A molecules is given by the sum of the external (electrostatic) and internal (polarization) contribution[54]where interactions between the multipole moment Q and Q are described by the interaction tensor T. Eq follows a variational principle with respect to the induced moments, and a self-consistent procedure iteratively updating ΔQ is required to find the minimum energy. Induced interactions are modified using Thole’s damping functions[55,56] to avoid overpolarization. Unlike ref (53) in which environment polarization effects are included explicitly in the GW equations as an additional screening contribution, we employ a double-level self-consistency cycle. At iteration step m, the potential generated by the static and induced moments of the MM region acting on the QM region is added to the external potential in eq 1, and a self-consistent converged QM calculation is performed yielding an electron density ρQM for the QM part. Since the excited state density for state S is not directly accessible in GW-BSE, we calculate it as ρ(r) = ρDFT(r) + ρe(r) – ρh(r), with the hole and electron contribution of the exciton to the density obtained from the electron–hole wave function according towhere χ(r,r) is the two-particle amplitude, introduced in eq . At the moment an equivalent of relaxed excited state densities as in ref (57) is not accessible as analytic gradients are not implemented, and their proper definition is considered among the theoretic challenges in GW-BSE.[58]

The energy of the QM region is (DFT for ground s = n, DFT+GW-BSE for excited s = x states) EQM = EDFT + δΩ. Once ρQM is obtained, an effective multipole representation {Q̃} is used in the next evaluation of the MM energy in eq . Since the QM electron density already contains the polarization response to the outside field, no atomic polarizabilities are added to the QM region representation in this step. These effective multipoles are thus used to determine (self-consistently) new induced dipoles in the MM region using eq , treating the whole system classically.

Obtaining the total energy at step m for the coupled QM/MM system requires the subtraction of the interaction energy of the QM charge distribution with the field generated by the total MM multipoles, already included in EQM. In this way, double counting is avoided. The total energy at step m is thus

The whole procedure is repeated until the change of total energy ΔEQMMM = |EQMMM – EQMMM|, as well as those of the individual contributions, is smaller than 10–4 eV. As stated before, this procedure is valid for sytems in the ground or excited state: calculating separately EQMMM in both cases and subtracting the two will give the excitation energy in the polarizable environment. This procedure assumes that the states of interest and, in particular, their localization characteristics on the QM cluster are easily identifiable.

The explicit state dependence of the coupled QM/MM system introduces another difficulty, in particular when excited states via GW-BSE are calculated. The solution of the BSE yields a spectrum of excitations, which are ordered according to their energy. These states can be energetically separated or very close, depending on the specific system. As a consequence, the index of a specific excitation of interest can vary for different external potentials at the individual steps of the QM/MM self-consistency procedure. It is therefore important to be able to identify the electronically excited state of interest during the calculation. In practice, a filtering of the total spectrum is employed which selects states according to some predefined property. Currently, the selectable properties are the oscillator strength f for optically active excitations and the amount of charge transferred (Δq) from one fragment to another for charge transfer states. For such filtering criteria to be applicable, it is implicitly assumed that the overall characteristics of the excited states do not change significantly during the QM/MM calculation.

Implementation

The theoretical concepts outlined in the previous section are implemented in the open source VOTCA-XTP software, available at github.com/votca. In the following, we briefly describe the most important implementational details as they pertain to GW-BSE.

Finding the solutions to the quasi-particle equations eq , and subsequently to the BSE as in eq , requires converged Kohn–Sham molecular orbitals, their energies, and the contribution of Vxc to them as a starting point. VOTCA-XTP can read this information from standard packages using Gaussian-type orbitals (GTOs) as basis functions {ψ(r)} to expressand currently provides interfaces to GAUSSIAN,[59] ORCA,[60] and NWChem.[61] The modular nature of the interfaces allows for straightforward extension to other packages, provided information about the atomic orbital function order and input/output files is available. Matrix elements ⟨ϕKS|Vxc|ϕKS⟩ needed in eq are numerically integrated using spherical Lebedev and radial Euler-Maclaurin grids as used in NWChem,[61] with XC functionals provided by the LibXC library.[62]

As an alternative, VOTCA-XTP also contains a minimal implementation of DFT with GTOs, which is currently limited to closed shell systems. One- and two-electron integrals are computed with modified recursive algorithms.[63,64] Initial guesses can either be constructed solving the Hamiltonian of a noninteracting system or by superposition of atomic densities.[65] Convergence acceleration can be achieved by mixing techniques using an approximate energy functional (ADIIS)[66] or the commutator of Fock and density matrix (DIIS).[67]

The most time-consuming step in a DFT calculation is commonly the computation of the electron–electron interaction integralswhere D̲ is the density matrix andare four-center integrals of Gaussian basis functions. VOTCA-XTP makes use of the RI-V approximation, in which the introduction of an auxiliary basis set with functions {ξν(r)} allows eq to be rewritten as[68]Here, (ν|μ)−1 is the inverse of the two-center repulsion matrixand (ij|ν) is the three-center repulsion matrix

The RI-V approximation, sometimes also referred to as density-fitting, reduces the scaling from N4 to N3, where N is the number of basis functions, and is particularly useful in application to large systems. It is also an integral part in the implementation of GW-BSE.

Solving the quasiparticle equations eq or eq relies on the determination of the self-energy Σ with the help of the dielectric function within the RPA. To avoid numerical instabilities in the calculation of the long-range part of ϵ, we introduce a symmetrized (with respect to r1 and r2) form of the Coulomb interaction ṽ(r1,r2) = π–3/2/|r1–r2|2, which is related via the convolution v(r1,r2) = ∫ṽ(r1,r′)ṽ(r′,r2)d3r′ to the actual Coulomb interaction. This can be shown by making use of the convolution theorem of Fourier transforms, which yields . With being the Fourier transformed of ṽ(r1,r2), it follows that v = δ4π/|G|2, which is the known Fourier transform of v(r1,r2) (see, e.g., the Appendix of ref (69) for details). The associated symmetrized dielectric function reads ϵ̃ = 1 – ṽPṽ or, equivalently, ϵ̃ = ṽ–1ϵṽ. Due to the simultaneous symmetrization of the Coulomb interaction and the dielectric function, the screened Coulomb interaction is obtained via a double convolution as

From eq it is apparent that the use of the symmetrized form of the Coulomb interaction does not change the screening behavior.

Using the RPA[41,69] for the polarizability (P = iG0G0) and the resolution of identity in eq , we can write ϵ̃μν in the auxiliary basis explicitly aswhere and (ν|μ)−1/2 is the matrix square root of the inverse of eq . With the definition of mixed molecular-atomic three-center Coulomb integralsone obtains the expression of the expectation values of the self-energy Σ(E) with respect to two Kohn–Sham orbitalswhere the plus (minus) sign in the ±iη term in the denominator is used if the Kohn–Sham orbital l is occupied (empty).

VOTCA-XTP employs a generalized plasmon pole model (PPM) as outlined in ref (69) to perform the frequency integration. This model allows for a quick evaluation of the integral in eq , but at the same time turns the self-energy into a real operator.[70] The PPM was chosen with the application to complex molecular systems of considerable size, e.g., with relevance to organic electronics such as polymer–fullerene clusters, in mind. The particular model used in this work has been successfully applied to determine quasiparticle and optical excitations in bulk semiconductor and insulator crystals,[71,72] their surfaces,[69,73] defect levels,[74] inorganic clusters,[70] and polymers,[27,75,76] as well as inorganic and organic molecules.[20,22,26,36,43,77] Explicit integration of the complex integral using (partially) analytic techniques[78−81] is planned for future versions. The respective matrix elements of the electron–hole interaction kernel in eq can be analogously expressed using eq and eq .

As mentioned in section , the use of the Kohn–Sham energies εKS as in eq and eq corresponds to the so-called single-shot G0W0 approximation. In VOTCA-XTP, the partial self-consistent GW0 and evGW schemes are available. Both schemes converge the quasiparticle energies εQP but not the quasi-particle states ϕQP. Fully self-consistent scGW, which diagonalizes HQP (eq ) in every iteration is currently not implemented. The respective steps of the GW-BSE calculation are depicted in Figure .

The general GW-BSE/MM procedure, together with the use of the PPM setting VOTCA-XTP apart from other (closed source) GW codes such as Turbomole[82] or Fiesta,[21] builds on the classical trajectory parsers of VOTCA-CSG[83] and the system partitioning functionality and electrostatic and polarizable interactions and potentials in the VOTCA-CTP library.[84] In addition to the core functionality described in this paper, VOTCA-XTP also contains visualization tools as well as modules for Mulliken[85] or Löwdin[86] population analysis, CHELPG[87] partial charge fitting for ground, excited, and transition densities with optional constraints, and numerical excited state gradients and geometry optimization. It provides methods for the calculation of nonadiabatic coupling elements for electrons, holes,[88] and singlet or triplet excitons[89] and links to a rate-based model of excited state dynamics using kinetic Monte Carlo techniques.

Results

In this section, VOTCA-XTP’s GW-BSE implementation is first tested using a small molecule reference set, also known as the Thiel set. After that, we consider as a prototypical complex molecular system double-stranded DNA with specific focus on the effects of local-electric fields and environment polarization on charge transfer excitations.

Single Molecule Data: Thiel Set

For benchmarking, the following procedure has been used. First, the ground state geometries of the molecules have been optimized using DFT with the hybrid PBE0 functional[90] at three different levels of theory, including all-electron (AE) calculations with the aug-cc-pVTZ and cc-pVTZ basis sets,[91] respectively, as well as calculations making use of effective core potentials and an associated basis set[92] that has been augmented by a single shell of polarization functions taken from the 6-311G** basis.[93] Due to the significantly reduced computational requirements, the latter case can be considered a minimal setup and is further referred to as min-ubecppol. Optimized auxiliary basis sets for (aug-)cc-pVTZ[94,95] taken from the Basis Set Exchange[96] have been used in the resolution-of-identity steps. For the min-ubecppol basis, we constructed an auxiliary basis using the technique employed in the SAPT code.[97−99] For all cases, we have compared the obtained results to those from calculations using large auxiliary bases created with the AutoAux functionality[100] available in Orca[60] and found agreement within a few 10 meV. A full list of size of the corresponding basis sets and auxiliary basis sets is given in Table S1 of the Supporting Information.

For the optimized geometries, excited state energies are determined within the GW-BSE formalism making use of the full BSE (eq ) on top of evGW self-consistent quasi-particle energies using the procedure outlined in section , in which all GW energies are converged to 10–5 Hartree. Transitions between all occupied and empty states, with their total number determined by the respective basis set sizes as in Table S1 of the Supporting Information, are taken into account in the calculation of the dielectric screening in the RPA. This choice is conservatively large, since including about 10 times as many empty as occupied states has typically shown to be sufficient to yield converged low energy excitations (see also Figure S1 in the Supporting Information). Similarly, quasi-particle corrections are determined for all available states, which are then used to construct the basis of product states for the expansion of the electron–hole wave functions in the BSE. For example in the smallest system, ethene, our calculations with the aug-cc-pVTZ basis include 8 occupied and 130 empty states, leading to 1040 transitions in the RPA and for the BSE product basis. For naphthalene, inclusion of 34 occupied and 610 empty states amounts to 20740 RPA transitions/BSE product functions.

From the resulting set of excitations for the respective molecules, the excitations with optical activity are identified, and their energies are compared to the ones obtained from experiment, as summarized in Figure . The four different categories of small molecules are represented by differently colored symbols (see caption for details). For the aug-cc-pVTZ basis set that contains additional diffuse functions, the results depicted in Figure (a) indicate a very good agreement with the reference data with a RMSD of 0.24 eV. The largest deviation is found for cyclopropene, whose excitation is reported to be at 7.19 eV in experiment compared to 6.38 eV in our GW-BSE calculation. Such a deviation is, however, not unique to our implementation. In ref (101), a GW-BSE excitation energy of 6.14 eV was reported, which is very close to the value of 6.18 eV obtained by TD-DFT with the PBE0 functional. Even the Theoretical Best Estimate based on high-order wave function methods of 6.65 eV shows a similar deviation. We note that the difference of some of our GW-BSE results from those in ref (101) is likely an effect of the different treatment of the frequency dependence of the dielectric functions (PPM vs complex contour integration). Overall we find a mean absolute error of 0.14 eV between our PPM approach and the literature results. For the moderately sized nucleobases, for which one would expect the PPM to be a better approximation, this error is as small as 0.03 eV. Such an error is negligible compared to the effects of the molecular environment on the excitation energies, which can be on the order of 1 eV (see Section ). Smaller deviations could also be attributed to more subtle variations in the computational protocol, such as in the stabilization of near linear dependencies in the basis sets and auxiliary basis sets. In general, a more direct comparison between the various theoretical approaches is made somewhat difficult by the fact that molecular geometries have been optimized at different levels of theory and therefore can distort the picture slightly.

Figure 2

Comparison of calculated lowest singlet excitation energies with experimental data (in eV) for the 28 small molecules in Thiel’s set. Ground state DFT calculations including geometry optimizations have been performed on an all-electron (AE) level with the (a) aug-cc-pVTZ and (b) cc-pVTZ basis sets, as well as (c) employing effective core potentials and the min-ubecppol basis set, respectively. The PBE0 functional was used. The same fitting auxiliary basis functions have been used for both DFT and GW-BSE stages. Data is given as follows: for nucleobases by green squares, for unsaturated aliphatic hydrocarbons by red up-triangles, formaldehydes, ketones, and amides by ocher diamonds, and for aromatic hydrocarbons by blue down-triangles. Panel (d) shows the relative error between the smaller AE/cc-pVTZ (green filled circles) and ECP/min-ubecppol (open blue circles) calculations as compared to the more complete AE/aug-cc-pVTZ as a function of energy. The green dashed (blue dotted) lines indicate the mean error of (3.2 ± 1.0)% and (6.3 ± 3.1)%, respectively.

To scrutinize whether our results are affected by the choice of functional in the underlying ground state DFT calculation, we have computed the respective excitation energies also with the gradient-corrected PBE[102] functional instead of its hybrid variant PBE0. The full data for both G0W0 and evGW variants are given in Table S2 in the Supporting Information. Inclusion of quasi-particle energy self-consistency reduces the mean-absolute error between the PBE0 and PBE functionals from 0.087 ± 0.053 eV to 0.052 ± 0.028 eV. The largest difference on the G0W0 level is 0.18 eV for formaldehyde, compared to only 0.02 eV with evGW. Overall, we note only a very weak starting point dependence, in particular for evGW.

Using diffuse basis functions in quantum-chemical calculations is typically associated with significant computational costs due to increased number of functions not only in the basis set itself but also in the auxiliary basis sets for RI. Concomitantly, one occasionally encounters problems with linear dependencies in the basis sets that require careful treatment. In this situation, it is desirable to avoid such diffuse functions, especially in applications to larger molecules. In Figure (b), the GW-BSE results obtained with the cc-pVTZ basis set show overall an excellent agreement with the experimental reference. On average, the RMSD of 0.28 eV is as expected larger than that for the aug-cc-pVTZ basis. This is illustrated in Figure (d), in which the relative deviation of the excitation energies (in %, indicated by green filled circles) obtained with cc-pVTZ from those obtained with the more complete aug-cc-pVTZ basis sets are shown depending on the absolute aug-cc-pVTZ energies. It can clearly be seen that on the energy range covered by the test set, the relative deviation varies between 1% and 9%, yielding a mean relative error of 3.2% with standard deviation of 1.0%. More importantly, however, the average run time is reduced to (25.2 ± 6.5)%.

While neglecting diffuse functions already massively reduces computational costs with only minimal loss of overall accuracy and reliability, all-electron calculations explicitly include the typically inert core electrons, such as the two electrons in the 1s shell of carbon. It is therefore possible to simply exclude them from the active space of product functions. However, the presence of such explicit core electrons requires the use of normal and auxiliary basis sets with strongly localized functions in the DFT ground state calculation underlying the GW-BSE formalism.

To avoid the expensive calculation of these core states altogether, effective core potentials can be used in combination with the min-ubecppol basis set. In Figure (c), the obtained excitation energies are shown compared to the experimental reference. The overall RMSD of 0.42 eV, while slightly larger than that recorded for aug-cc-pVTZ and cc-pVTZ, respectively, is still very good. One can observe a general tendency for the ECP/min-ubecppol combination to overestimate the measured data. This is also apparent considering the relative deviations from aug-cc-pVTZ shown as open circles in Figure (d). Interestingly, the relative deviation varies between 1% and 10%, only slightly larger than for cc-pVTZ. However, the mean error is larger and amounts to (6.7 ± 2.0)%, which can be considered acceptable, in particular when one takes into account that the computational cost is reduced to as much as (6.3 ± 3.1)% as compared to aug-cc-pVTZ. These numbers highlight that the use of the minimal ECP/min-ubecppol variant offers a great compromise between accuracy and computational cost, which make it particularly attractive for the application to large, relevant molecular systems.

For completeness, a comparison of the electronic excitation energies obtained with GW-BSE to the Theoretical Best Estimate (TBE) clearly reveals that all three basis set variants considered in this work exhibit a very satisfying agreement with the high-order reference. The data and a figure are available in Figure S2 and Table S3 in the Supporting Information.

Additional savings can in principle be achieved by resorting to the Tamm-Dancoff Approximation (TDA), in which as explained in section the resonant-antiresonant coupling terms are neglected in the Bethe-Salpeter Equation. The dimension of the matrix system is reduced by a factor of 2 which directly translates into significant numerical gains. This omission of the corresponding coupling terms in the BSE can reduce the associated energies by several 0.1 eV, depending on the size of the π-conjugated system. The smaller the π-system, the stronger the effect. For the relatively small molecules in the Thiel test set, it is therefore expected that the TDA deviations will be noticeable.

Also for the Thiel set, the TDA energies are typically larger than those from the full BSE, see Figure S3 and Table S3 in the Supporting Information. Also the size dependence is clearly visible. The strongest effects can be seen for ethene (C2H4), the molecule with the smallest π system. For the aug-cc-pVTZ basis, TDA yields an excitation energy of 8.04 eV as compared to 7.51 eV obtained by the full BSE formalism. Resonant-antiresonant coupling accounts for as much as 0.53 eV in this case. In contrast, for a larger molecule such as adenine, the effect is reduced to just 0.02 eV. These results illustrate that the TDA can be a useful approximation depending on the specific system of interest and should therefore be carefully evaluated.

With these promising conclusions regarding the application to single small molecule systems at hand, the following section will focus on the integration of GW-BSE in coupled quantum-classical QM/MM setups for complex molecular environments.

Charge Transfer Excitations in Double-Stranded Aqueous DNA

To obtain the atomistic structural information, an exemplary DNA double strand with 23 base pairs in the sequence shown in Figure was prepared. This double strand was solvated by 42216 water molecules and 44 sodium counterions. For this system, in the following referred to as aqDNA, classical molecular dynamics simulations were performed using the AMBER99 force field[103] for DNA and sodium and the SPC/E water model.[104] Geometric mixing rules [ and ] for Lennard-Jones (LJ) diameters (σ) and LJ energies (ϵ) were used for atoms of different species.[105−107] Nonbonded interactions between atom pairs within a molecule separated by one or two bonds were excluded. Interaction was reduced by a factor of 1/2 for atoms separated by three bonds and more. Simulations were run using GROMACS version 5.[108] A 0.9 nm cutoff was employed for the real space part of electrostatics and Lennard-Jones interactions. The long-range electrostatics were calculated using particle-mesh Ewald (PME)[109,110] with the reciprocal-space interactions evaluated on a 0.16 grid with cubic interpolation of order 4. First, the system was energy minimized using the steepest descents algorithm. Then, 10 ns simulations in constant particle number, volume, and temperature (NVT) ensemble at 300 K were performed using the stochastic velocity rescaling thermostat[111] with time constant of 0.1 ps. The velocity-Verlet algorithm[112] was employed to integrate the equations of motions with a 2 fs time step. The simulation box size was (12 × 12 × 8) nm3. Simulations were then continued in constant particle number, pressure, and temperature (NpT) ensemble at 300 K and 1 bar controlled by the Parrinello–Rahman[113] barostat with a coupling time constant of 2.0 ps. Molecular visualizations were done using Visual Molecular Dynamics (VMD) software.[114]

Figure 3

DNA strand sequence.

Figure illustrates the partitioning of the MD system of the solvated DNA double strand into QM and MM regions. A single nucleobase or a base pair is chosen as the QM region, while the rest of the system that is within a certain distance is assigned to the MM region. We differentiate between two distinct MM regions, here referred to as MM0 and MM1. In MM0, both static and polarizable effects are taken into account, while in MM1, only static multipoles are considered. In this particular case, we restrict the static multipoles to point charges[87] and the induced moments to dipoles.

Figure 4

Schematic representation of aqDNA and separation into MM0 and MM1 for an adenine nucleobase. The QM region is seen in the small inset.

For the parametrization of the polarizable model used in the coupled QM/MM calculations, atomic partial charges and molecular polarizability tensors were determined for the nucleobases and for water based on DFT calculations using the PBE0 functional and the cc-pVTZ basis set. Classical atomic polarizabilities were then optimized to reproduce the molecular polarizable volume of the DFT reference calculation. For the DNA backbone, partial charges were taken from the force field used in the MD simulation, and the default atomic polarizabilities in the AMOEBA force field[115] were employed. Either a single nucleobase or a pair of nucleobases is chosen as the QM region in the QM/MM setup. As this region is covalently bonded to the MM region, the bond to the frontier atom was truncated and saturated with a hydrogen atom. All residues within a closest contact distance of 4.3 nm to the molecules defining the QM region were assigned to the MM region. When polarized QM/MM calculations were performed, polarization effects were included for all residues within a closest contact distance of 2.0 nm.

From the simulated DNA structure, neighboring nucleobases with separation less than 1 nm were defined as pairs (yielding 59 pairs in total) between which CT excitations are calculated. These include both intra- and interstrand excitations. Due to the presence of the four nucleobases adenine (A), guanine (G), cytosine (C), and thymine (T) in the present system of aqDNA, 10 different types of dimers can be formed.

First, we compare the results obtained using QM calculation of gas-phase dimers with those obtained using QM/MM with only static classical interactions. Figure shows the distribution of CT exciton energies for both cases. We refer to these distributions and their Gaussian broadened guide-to-the-eye as CT Density of States (DOS). CT DOS for dimers in vacuum are represented by blue bars, while red bars indicate results for QM/MM dimers embedded in a static background. The inset labels show the base pair combination and in parentheses the total number of CT states found in vacuum and static QM/MM, respectively. In all cases, an excitation was labeled as a CT state if the charge transfer between the two nucleobases exceeded 0.5 e.

Figure 5

Density of states (DOS) for charge transfer (CT) excitations in aqDNA as obtained from dimers in vacuum (blue bars) and QM/MM embedded in a static background of point charges (red bars), respectively. The individual panels show different base pair combinations, in which neighboring nucleobases within a closest contact distance of less than 1 nm are considered as pairs. Due to the specific sequence of the model strand used in this work, different numbers of pairs are found for each combination. The inset labels indicate both the type of combination and in parentheses the total number of CT states found in vacuum and static QM/MM. A cutoff of 4.3 nm was used for the atomistic electrostatic embedding.

A general observation is that the total number of CT states found in the covered energy region of 5 to 9 eV is always larger in the QM/MM case. This observation can be attributed to two effects. First, some of the CT states that fall outside of the energy interval in the gas-phase calculation get pushed down in energy to values below 9 eV in static QM/MM. Second, some of the CT states change their character by embedding in the static background.

A more detailed analysis of the changes in distributions, in particular in the low energy regions, shows no universal behavior. In some cases such as for the adenine dimers (A-A), some individual excitations demonstrate lower energies in static QM/MM than in the gas-phase. While not resolved in Figure , the lowest energy CT excitation at about 5.35 eV is an intrastrand adenine dimer of the kind previously discussed by Yin et al.[36] in a more idealized structure. We will scrutinize the properties of this particular excitation in more detail below.

In contrast to the behavior of A-A pairs, dimers formed from two cytosine bases exhibit CT excitons at higher energy than in the respective gas-phase calculation, irrespective of whether it is an intra- or interstrand excitation, cf. Figure .

Given the nonuniversal behavior observed upon inclusion of a static environment, we limit the following discussion to only the lowest energy CT excitation in each of the 59 pairs. The aim is to understand the additional influence of polarization in the GW-BSE/MM calculations. In Figure , CT excitation energies resulting from both static (open symbols) and polarized (closed symbols) calculations are shown against the respective vacuum energy, also resolving intrastrand (circles) and interstrand (squares) excitations. As in the static case, no general trend can be discerned. CT excitation energies are both lowered and raised due to the presence of the environment. There appears to be a tendency that the lower-energy interstrand CTs up to an energy of 7 eV are all resulting at about 0.5 eV higher energies in the static case.

Figure 6

Comparison of CT excitation energies (in eV) calculated in static (open symbols) and polarizable (filled symbols) QM/MM setups with vacuum QM results. Interstrand (intrastrand) excitations are represented by green squares (red circles). The gray shaded areas indicate the range of single nucleobase UV absorption energies of adenine.

Taking polarization effects into account universally lowers the energy, not only with respect to the static QM/MM results but also most importantly with respect to the vacuum calculation. On average, we observe a redshift of the interstrand CT energies by (−0.83 ± 0.5) eV, while intrastrand CTs are red-shifted by (−1.15 ± 0.6) eV, compared to respective vacuum results. Notably, these redshifts are on the order of the redshift observed in experiment. Also, the CT excitation with the lowest energy of 4.81 eV is found for a A2 dimer in the chain.

In addition to the individual CT energies, the gray shaded areas in Figure indicate the energy range in which single adenine nuclobases absorb UV light, according to gas-phase and QM/MM calculations. While not shown here explicitly, the inclusion of a polarizable environment does not affect the energetic properties of these localized Frenkel excitons perceptively, with absorption predicted to be in the range (5.12 ± 0.02) eV. The lowest energies of CT excitations found in our data set are approximately 0.3 eV below this absorption energy, indicating that the decay of the UV excitation to a CT excited state is energetically possible, as speculated.

Due to this energetic situation, it is worthwhile to analyze the A2 CT exciton in further detail and to illustrate how the atomistic environment affects not only its energy but also its electron–hole wave function. To this end, we show in Figure the distributions of electron and hole densities on the A2 dimer for (a) vacuum QM, (b) static QM/MM, and (c) polarized QM/MM, respectively. The associated excitation energies and effective charge transfer are indicated below. As discussed before, for the vacuum case the CT energy of 5.78 eV is several 0.1 eV above the energy of the UV active excitation. The amount of charge transferred in the CT state is only 0.6 e, with the hole contribution on the lower nucleobase (AL) and the electron contribution on the upper one (AU). Upon inclusion of the static environment, the energy of this excitation is lowered by 0.44 to 5.34 eV, while the amount of charge transferred between the two adenines remains at 0.6 e. Despite this similarity, the characteristic of the excitation is changed significantly, as can be seen in Figure (b). The localization of electron and hole contribution in the excitation is inverted. Including polarization effects, the general character of the CT excitation remains unaffected, i.e., the hole is localized on AU and the electron on AL. Most notably, however, the excitation exhibits integer charge transfer character in this situation.

Figure 7

Isosurfaces (±2 × 10–3 e/Å3) of differential electron densities of the lowest energy adenine dimer resulting from (a) a gas-phase (vacuum) calculation, (b) a QM/MM calculation with static environment, and (c) a QM/MM calculation with polarizable environment. Red color corresponds to negative values (hole density), and blue color corresponds to positive values (electron density).

The observation that the nature of the CT excitation can be affected dramatically by the complex molecular environment can be attributed to a combination of a shift of energy levels and changed composition of transitions. We analyze the quasi-particle energy levels obtained at the GW step of the respective calculations. Figure shows the energies of two highest occupied and two lowest empty quasi-particle levels for vacuum, static, and polarized calculations. Note that for an easier comparison, the zero of the energy scale has been set to the center of the HOMO–LUMO gap in all individual cases. The spatial distribution of all quasi-particle wave functions has been inspected and is indicated by the horizontal lines’ color. Brown (dark green) lines indicate states that are localized on AL (AU). In addition, the vertical lines show the contributions of the quasi-particle transitions to the respective CT excitations in Figure , with the weights given in the inset.

Figure 8

Quasi-particle energy levels (eV) for HOMO–1, HOMO, LUMO, and LUMO+1 resulting from (a) a gas-phase (vacuum) calculation, (b) a QM/MM calculation with static environment, and (c) a QM/MM calculation with polarizable environment. The color of horizontal lines indicates the localization of the quasi-particle states on either of the two nucleobases. Brown (dark green) represents localization on AL (AU). For HOMO–1 and HOMO in the vacuum case, the quasi-particle states are distributed over the whole base pair, which is noted as a dashed line. Vertical arrows show the dominant transitions forming the CT excitation.

Figure 9

Effective charge transfer character (in e) in the CT excitations as a function of center-of-mass distance of the involved monomers (in nm). Results for intra- and interstrand excitations are compared for the three different calculation setups: vacuum, static QM/MM, and polarized QM/MM.

In the case of the vacuum calculation on the adenine dimer taken from the MD snapshot, it turns out that the two occupied levels cannot be uniquely assigned to either of the two nucleobases. Instead, the quasi-particle states delocalize over the dimer, however not at equal distribution. Note, though, that they are only separated by 0.13 eV in energy. To make this also visually clear, the two levels are shown as dashed lines in Figure . As can be seen from the two arrows, the CT excitation in this environment-free QM calculation is composed of HOMO–1 → LUMO and HOMO → LUMO transitions with nearly equal weight. The fact that the combined weight is only 70% emphasizes that even more quasi-particle transitions play a significant role here. It is likely that this is directly linked to the delocalized nature of the occupied states. Taken as a whole, the hole contribution of the CT, arising in large parts from the HOMO and HOMO–1 states, is consequently localized on AL. For the two unoccupied levels shown here, no strong delocalization over the dimer can be identified. Since the LUMO is localized on AU, also the electron density in the CT state is found on this nucleobase.

Turning now toward the results obtained from calculations performed in the static QM/MM setup, one can spot significant changes as compared to the vacuum only calculation. First, all quasi-particle states around the HOMO–LUMO gap are localized on either of the two nucleobases of the excimer. In the occupied manifold, one can now assign the HOMO to be uniquely localized on AU and HOMO–1 on AL. As a consequence the energetic separation is more pronounced, amounting to 0.62 eV. At the same time the two unoccupied states change character. While also localized on either of the two nucleobases in the vacuum calculation, one finds that the specific localization site is switched. The LUMO is now localized on AL, and LUMO+1 is localized on AU. Combined with the fact that the dominant transition in the CT excitation is a HOMO to LUMO transition from AU to AL with a weight of approximately 60%, see Figure , the localization behavior of hole and electron densities is inverted as compared to the vacuum case. The total transferred charge however remains 0.6 eV, which can be attributed to the additional transitions that collectively contribute to 40% of the excited state.

We note in passing that the HOMO–LUMO gap is also slightly reduced by embedding in a static molecular environment, namely from 9.07 to 8.64 eV. A word of caution: The fact that the reduction by 0.43 eV of this gap is numerically similar to the lowering of the CT excitation energy by 0.44 eV is likely coincidental. Typically, a change in localization of the contributing quasi-particle states leads to a very different composition of the effective electron–hole interaction that determines the exciton binding energy and, concomitantly, the excitation energy.

From the quasi-particle levels in the polarized QM/MM calculation as shown on the right-hand side of Figure , one can see that the environment polarization response modifies this picture even more. First of all, now one finds the two occupied states shown being localized on the upper adenine nucleobase, and the two unoccupied states being localized on the lower one. The HOMO–LUMO gap is further reduced to 7.52 eV, and the energetic separation of the occupied and unoccupied levels is increased. Most remarkable is now that the CT excitation is in this case given as a pure HOMO to LUMO transition. The hole and electron contributions are fully localized on AU and AL, respectively, corresponding to integer charge transfer.

The above detailed analysis of the characteristics of the quasi-particle and CT excited states for the minimum energy CT found in our data set clearly reveals that the resulting excitation energies in complex molecular environments obtained from QM/MM calculations are a result of an intricate interplay of several effects. In particular, modifications on the nature of the quasi-particle states are significant since their localization/delocalization characteristics have a profound and direct effect on the two-particle excitations. This interplay cannot be captured by adding a perturbative energy correction due to the environment to a vacuum QM calculation.

To scrutinize whether the change of effective charge transfer in the CT excitation observed for the intrastrand adenine dimer observed above is a more general effect of embedding into a static and/or polarizable molecular environment, we show in Figure the calculated amount of transferred charge as a function of center-of-mass distance for the various calculation setups. We differentiate also between intra- and interstrand excitations.

It can be seen for the excimers with the closest intermolecular separation between 3 and 4 Å, which are exclusively intrastrand excitations, vacuum calculations yield only partial charge transfer upon excitation between 0.5 e and 0.9 e. The same holds for interstrand dimers with distances of 5 and 6 Å. All these short distance dimers are essentially neighboring molecules whose electron density can spatially overlap and the associated interaction yielding (partially) delocalized quasi-particle states. For all dimers with separation larger than 0.7 nm center-of-mass distance, i.e., second-nearest neighbors, such a direct interaction is not possible. In the case of intrastrand excitations, it means that in a stack of three bases (base trimer), only the outer two nucleobases are treated quantum-mechanically, while the center one is part of the polarizable MM region. This is strictly speaking a fairly strong approximation. When base stacking interactions are strong, the purely classical treatment cannot cover possible effects of forming delocalized states and the associated partial charge transfer. Also, such explicit base pair interactions might affect the CT excitation energies directly. A possible pathway to cover such effects is to treat the full base trimer quantum-mechanically and embed this in a classical environment. However, this case goes beyond the scope of this work and is left for future studies.

We focus in the following on the short-distance excimers. When the molecular environment is taken into account, the static-only interactions (open symbols in Figure ) affect the amount of effectively transferred charge roughly in the same fashion as observed for the A2 system with minimal CT excitation energy discussed above. In some cases, one can note a change of this effective charge by up to 0.3 e. However, at least for the first shell of intra- and interstrand dimers, there is no observable integer charge transfer state.

Only upon adding environment polarization effects (filled symbols in Figure ), most of the CT states are approaching such an integer CT character. It stands to reason that remnant delocalization for quasi-particle states is responsible for that.

Summary

In this paper, the Gaussian-orbital based implementation of many-body Green’s functions theory within the GW approximation and the Bethe-Salpeter Equation (BSE) in the open-source VOTCA-XTP software has been introduced. Application to the standard small molecule Thiel set has been used to benchmark the obtained excitation energies. The results are in very good agreement with the experimental reference for a variety of excitation types and an energy range from 2 to 8 eV, validating both the methodology and its implementation.

It has further been demonstrated how coupling GW-BSE to a classical atomistic environment in QM/MM schemes allows studying electronic excitations in complex molecular environments, here in prototypical aqueous DNA. It is found that charge transfer excitations are extremely sensitive to the specific environment. For the lowest energy CT excitations in an intrastrand adenine dimer, the approach predicts energies below that of the UV active single nucleobase excitation. This has a tremendous impact on the possibility of an initial (fast) decay of such an UV excited state into a binucleobase CT exciton, which is considered one of the pathways for UV-induced DNA damage. The calculated redshift of the CT excitation energy compared to a nucelobase dimer treated only in vacuum is of the order of 1 eV, which matches expectations from experimental data. The GW-BSE/MM methodology used here allows for gaining very detailed insight into the mechanisms leading to the observed energies. It is possible to disentangle the effects of the different levels of the explicit molecular environment on single-particle and two-particle excitations. Incorporating GW-BSE into the presented QM/MM setup is therefore an extremely powerful tool to study a wide range of types of electronic excitations in complex molecular environments.