$Δ$ADAPT-VQE: Toward Accurate Calculation of Excitation Energies on Quantum Computers for BODIPY Molecules With Application in Photodynamic Therapy
Anton Nykänen Leander Thiessen Elsi-Mari Borrelli Vijay Krishna Stefan Knecht Fabijan Pavošević























Quantum chemistry simulations offer a cost-effective way for computational design of BODIPY photosensitizers with potential use in photodynamic therapy (PDT). However, accurate predictions of photophysical properties, such as excitation energies, pose a challenge for the popular time-dependent density functional theory (TDDFT) and equation-of-motion coupled cluster with singles and doubles (EOM-CCSD) methods. By contrast, reliable descriptions can be achieved by multi-reference quantum chemistry methods, though unfortunately, their computational cost grows exponentially with the number of correlated electrons. Alternatively, quantum computing holds a great potential for exact simulation of photophysical properties in a computationally more efficient way. To this end, we introduce the state-specific ΔUCCSD-VQE (unitary coupled cluster with singles and doubles variational quantum eigensolver) and ΔADAPT-VQE methods in which the electronically excited state is calculated via a non-Aufbau electronic configuration. The accuracy and capability of the developed methods are  assessed against experimentally determined excitation energies for six BODIPY-derivatives. We show that the proposed methods predict accurate vertical excitation energies that are not only in good agreement with experimental reference data but also outperform popular quantum chemistry methods, such as TDDFT and EOM-CCSD. Spurred by its impressive performance and simplicity, we are confident that ΔADAPT will emerge as the method of choice for guiding the rational design of photosensitizers for PDT and photocatalysis in the era of near-term quantum computing.



Fast photon-mediated entanglement of continuously-cooled trapped ions 
 for quantum networking
    Christopher Monroe
    April 29, 2024
==============================================================================================





§ INTRODUCTION

Boron-dipyrromethene (also known as BODIPY) and its derivatives have emerged as an important class of organic fluorescent dyes characterized by high molecular absorption coefficients, high fluorescent quantum yields, and excellent thermal and photochemical stabilities <cit.>. These favorable photophysical properties can be further tuned by appropriate structural modifications of the BODIPY core <cit.>, making them versatile compounds for medical imaging <cit.>, labelling <cit.>, photoelectrochemistry <cit.>, artificial photosynthesis <cit.>, optoelectronics <cit.>, and
photovoltaics <cit.>. BODIPYs can also be used in photodynamic therapy (PDT) as a photosensitizer (PS) <cit.>. PDT is an emerging procedure for noninvasive cancer treatment (including bladder, prostate, lung, breast, skin, and esophagus cancers) in which the PS upon light activation transfers electrons (Type I PDT mechanism) or energy (Type II PDT mechanism) to the surrounding environment to generate cytotoxic reactive species that cause cellular death (apoptosis and necrosis) <cit.>. The Type II PDT mechanism <cit.> mainly occurs through the following three steps: (i) upon light absorption, the PS is excited to the first excited state S_1; (ii) a subsequent population of the first triplet (T_1) state via non-radiative intersystem crossing (ISC); (iii) energy transfer from PS in T_1 state to cellular triplet molecular oxygen to form reactive oxygen species (ROS). Unlike Type I, the Type II mechanism maintains the chemical form of the PS which allows for multiple therapeutic cycles during a single procedure. 

For BODIPYs to be useful in the PDT, they need to meet certain minimum criteria <cit.>, that are, non-toxic PS at the ground state S_0 <cit.>, a high and sharp absorption maximum within the tissue transparency window (750-900 nm), <cit.> and a large S_1/T_1 spin-orbit coupling (SOC) to allow for an efficient formation of T_1 state with a high enough energy (>0.98 eV <cit.>) to generate ROS. To design functionalized BODIPYs with optimal photophysical properties mentioned above, quantum chemistry simulations play an important role <cit.>. From a theoretical standpoint, a computational aided design entails the calculation of singlet (S_1) and triplet (T_1) excitation energies (and potentially higher-excited states) which are important for the prediction of absorption spectra, fluorescent quantum yield <cit.>, SOC <cit.>, and triplet lifetime. Because the S_1 is essential for many of these photophysical properties, the focus of the present work is primarily on the calculation of vertical excitation energies (S_1). Unfortunately, the most popular and efficient quantum chemistry method for excited states calculations, namely time-dependent density functional theory (TDDFT), is not able to accurately predict the photophysical properties of BODIPY dyes <cit.>. The predictive power of TDDFT strongly depends on the selected exchange-correlation functional exhibiting  errors in calculated S_1 excitation energies of 0.3-0.6 eV compared to experimentally determined values for a representative set of 17 chemically diverse BODIPY derivatives <cit.>. Further improvements are obtained by employing various single-reference wave function approaches, such as equation-of-motion coupled cluster with singles and doubles (EOM-CCSD), second-order coupled cluster (CC2), scaled-opposite spin configuration interactions with singles and perturbative doubles (SOS-CIS(D)), and symmetry adapted cluster configuration interaction (SAC-CI), 
that yield on average an error on the order of 0.3 eV for calculated S_1 vertical excitation energies compared to experiment <cit.>. The most accurate vertical excitation energies are obtained by employing the multi-reference complete-active space second-order perturbation theory (CASPT2) method with errors less than 0.2 eV relative to experiment <cit.>. Despite its success, the accuracy of CASPT2 strongly depends on the active space size <cit.> and on the value of the shift parameter used to alleviate the intruder-state problem <cit.>. Moreover, since the number of configurations grows exponentially with increasing active space sizes, CASPT2 remains – for practical purposes – unsuitable for black-box applications to large molecular systems <cit.>. 

By contrast, recent advances of quantum computing technology and quantum algorithms <cit.> have opened a path to cope with the cost of exponential scaling. Therefore, quantum computing holds great promise for exact simulation of molecular processes in excited states involved in PDT. The most popular quantum computing algorithm for the currently available noisy intermediate-scale quantum (NISQ) devices <cit.>, is the hybrid quantum-classical variational quantum eigensolver (VQE) <cit.>, which employs the quantum device for handling only classically intractable parts of the computation. It is used to find the ground state energy of a given molecular electronic Hamiltonian. The VQE algorithm has also been extended for excited state calculations, for example, by means of the quantum equation-of-motion (qEOM) approach <cit.> or a quantum subspace expansion (QSE) framework <cit.> to name just a few. The overall accuracy and efficiency of any VQE ansatz is determined by the choice of the wave function parametrization, and in its original implementation <cit.>, it exploited a fixed unitary coupled cluster with singles and doubles (UCCSD) ansatz <cit.> that demonstrates good performance for selected multi-reference problems <cit.>. However, a genuine implementation of UCCSD-VQE on NISQ devices requires very deep circuits, limiting its applicability only to the smallest chemical systems <cit.>. Alternatively, instead of employing a fixed UCCSD ansatz, the subsequently proposed  
ADAPT-VQE (Adaptive Derivative-Assembled Problem-Tailored) <cit.> algorithm grows the ansatz iteratively, leading to significantly shallower circuits than those of UCCSD-VQE <cit.>.

To lay the foundation for accurate simulations of vertical excitation energies for BODIPY molecules on the NISQ devices, in this work we introduce ΔUCCSD-VQE and ΔADAPT-VQE, respectively, which are extensions of the classical ΔSCF <cit.> (self-consistent field) and ΔCCSD <cit.> methods within a quantum computing framework. In the remainder of the text, we will refer to these methods simply as ΔUCCSD and ΔADAPT. Both ΔUCCSD and ΔADAPT are state-specific approaches in which the excited state is calculated via a non-Aufbau electronic configuration. We assess the performance and accuracy of ΔUCCSD and ΔADAPT by comparing with  experimentally determined S_1 excitation energies of six BODIPY systems (see Fig. <ref>A). Moreover, we discuss quantum resource estimates for an implementation of the two approaches on near-term quantum devices. The developments and numerical performance analysis emphasise the efficiency and accuracy of ΔUCCSD and ΔADAPT geared toward the simulation of molecular excitation energies on the contemporary quantum computers. Moreover, this work lays the foundation for the development of other quantum-hardware efficient methods applicable to the calculation of excited-state energies and spectroscopic properties. 





§ THEORY

The VQE <cit.> is a quantum computing algorithm used for finding the ground state energy of a molecular system defined by the electronic molecular Hamiltonian

    Ĥ = h^p_q a^q_p + 1/2g^pq_rsa^rs_pq

In this equation, a^q_p = a_q^†a_p and a^rs_pq = a_r^†a_s^†a_qa_p are the one- and two-electron excitation operators (expressed in terms of fermionic creation/annihilation, a^†/a, operators), respectively, whereas h^p_q=⟨ q |ĥ|p ⟩ and g^pq_rs=⟨ rs|pq ⟩ are one- and two-electron integrals in the spin-orbital basis, respectively. Througout this work, summation over repeating indices is implied. Moreover, indices denoted by p,q,r,s,... correspond to general spin-orbitals, whereas indices denoted by i,j,k,l,... and a,b,c,d,... correspond to occupied and unoccupied (virtual) spin-orbitals, respectively. The ground state VQE energy, E_VQE, is obtained by minimizing the expectation value of the molecular Hamiltonian (defined in Eq. (<ref>)) 

    E_VQE=θmin⟨Ψ(θ)|Ĥ|Ψ(θ)⟩

with respect to parameters θ of some trial wave function |Ψ(θ)⟩.

In UCCSD-VQE, the wave function reads as

    |Ψ_UCCSD⟩=e^T̂-T̂^†|0⟩

where T̂=T̂_1+T̂_2=θ^i_aa_i^a+1/4θ^ij_aba_ij^ab is the cluster operator with singles and doubles and  |0⟩ denotes a (single-determinantal) reference wave function, commonly coinciding with the Hartree-Fock (HF) determinant. UCCSD-VQE can be implemented with a polynomial number of one- and two-qubit quantum gates on NISQ devices, although in practice, even for small systems their number is too high, rendering this method impractical for larger molecular applications <cit.>. To push toward larger scale simulations, the ADAPT-VQE algorithm adaptively grows the wave function from a predefined set of the generalized fermionic pool of operators {a_p^q-a^p_q,a_pq^rs-a^pq_rs} <cit.>. At each iteration of the ADAPT-VQE procedure, the operator with the highest electronic-energy gradient is added to augment the current trial wave function. Consequently, favorable features of the ADAPT-VQE algorithm are its accuracy paired with a resulting shallow circuit that allows for quantum chemistry simulations on NISQ devices. <cit.>

As previously stated, the VQE algorithm is commonly used to find the ground state energy of a molecular system. In this work, we lift this limitation and introduce ΔUCCSD and ΔADAPT, respectively, that enable us to locate a higher energy solution. Prediction of the excitation energies with the ΔUCCSD or ΔADAPT methods, are obtained from two separate calculations; the ground state and the excited state energy calculations. The ground state energy is obtained by performing the usual UCCSD-VQE or ADAPT-VQE calculation with the Aufbau reference configuration. For obtaining the excited state energy, in the first step, the state-specific molecular orbitals (MO) of the non-Aufbau HF reference configuration are optimized via the SCF procedure (also known as ΔSCF <cit.>). Such non-Aufbau reference configuration is constructed by placing one or more electrons in higher lying orbitals instead of the lowest lying orbitals (Aufbau). In the second step, these new set of orbitals are utilized to construct the state-specific molecular Hamiltonian defined in Eq. (<ref>). Finally, using the Hamiltonian from the previous step within the VQE algorithm provides the excited state UCCSD-VQE or ADAPT-VQE energies. Once the ground and excited state UCCSD-VQE or ADAPT-VQE energies are obtained, their difference defines the ΔUCCSD or ΔADAPT excitation energy, respectively. The schematic representation of the overall procedure is given in Fig. <ref>B.

To reach quantitative predictions with quantum-chemical calculations, it is necessary to employ large and flexible basis sets. Unfortunately, such large basis sets drastically increase the number of qubits required for quantum simulations. The frozen natural orbital (FNO) approximation offers a robust path for reducing the number of qubits without sacrificing the accuracy <cit.>. Within the FNO approximation, the FNO orbitals (defined as the eigenvectors of the virtual-virtual block of the one-particle density matrix, γ_a^b=⟨Ψ|a_a^b|Ψ⟩) are used to span the unoccupied orbital space. Moreover, their importance is ranked by their eigenvalues which correspond to the natural occupation numbers. Therefore, FNOs with larger natural occupation numbers will have a greater impact on the correlation energy than those with the smaller natural occupation numbers. Because the FNOs provide a more suitable basis for spanning the unoccupied orbitals compared to HF, for the same number of FNO and HF unoccupied orbitals, the FNO approach recovers significantly more correlation energy. In this work, the FNOs are obtained from the density matrix that employ the low-cost first-order Møller–Plesset (MP1) wave function, |Ψ_MP1^(1)⟩=1/4(g^ij_ab-g^ij_ba)/(ϵ_i+ϵ_j-ϵ_a-ϵ_b)a_ij^ab|0⟩, where the ϵ's are the HF orbital energies. Lastly, for recovering the correlation energy that is outside of the reduced set of FNO orbitals, the missing correlation energy due to truncation is accounted by the following correction Δ E_MP2=E_MP2^MO-E_MP2^FNO <cit.> that is added to the final energy. In this equation, E_MP2^MO and E_MP2^FNO are the energies calculated with the second-order Møller–Plesset (MP2) method in the untruncated MO basis and in the truncated FNO basis, respectively.



§ RESULTS AND DISCUSSION

All of the reported results are obtained on the geometries of the six BODIPY systems compiled in Fig. <ref> as obtained from Ref. . The TDDFT, EOM-CCSD, ΔSCF, and ΔCCSD calculations were performed with the  quantum chemistry software <cit.> employing the cc-pVDZ basis set <cit.>. The one-particle (h) and two-particle (g) molecular integrals that enter Hamiltonian defined by Eq. (<ref>) were obtained with the in-house modified version of  library <cit.>. In all calculations, effects due to solvation are accounted for by the conductor-like polarizable continuum model (CPCM) <cit.>. The state-specific SCF procedure (ΔSCF) utilizes the unrestricted orbitals and the maximum overlap method <cit.> for avoiding the variational collapse to the ground state. Fermionic operators were mapped into qubit space with the Jordan-Wigner mapper <cit.>. The ΔUCCSD and ΔADAPT calculations are performed within  <cit.> and Algorithmiq’s software framework , respectively, employing a noise-free quantum statevector simulator model. The reported ΔUCCSD results are performed with (4,4), (6,6), and (8,8) active spaces, where the first value represents the number of electrons, and the second, the combined number of active occupied and natural (unoccupied) orbitals. The reported ΔADAPT results are performed with (6,6) and (8,8) active spaces. The reported ADAPT results utilize the procedure in which operators are selected from the operator pool consisting of majorana excitations (fermionic excitations split into invidiual Pauli strings, but without removing the Z-chains as usually done when using the Jordan-Wigner mapping) <cit.>. The number of two-qubit controlled NOT (CNOT) gates was estimated using the  transpile pass <cit.> for all-to-all qubit connectivity with optimization level set to 3.



To test the accuracy of the developed ΔUCCSD and ΔADAPT quantum computing methods, we have calculated the vertical excitation energies (S_0 →S_1) for a set of six BODIPY systems (Fig. <ref>A) and compared the predicted excitation energies to the experimentally determined values (experimental data are obtained from Refs. ). Figure <ref> shows the individual absolute errors (AE) and mean absolute errors (MAE) of the calculated vertical excitation energies from experimental data, as obtained with different quantum chemical methods. Numerical values of the excitation energies calculated with different methods along with MAE and standard deviation (STDEV) are collected in Table S1 of the Supporting Information (SI). The MAE for the TDDFT calculations with three different exchange-correlation functionals, BLYP <cit.> (generalized gradient approximation), B3LYP <cit.> (hybrid), and CAM-B3LYP <cit.> (range-separated hybrid)
are 0.199 eV, 0.368 eV, and 0.392 eV, respectively. The same trend of performance for different functionals is also reported on a larger set of BODIPY systems <cit.>. Overall, we observe that the TDDFT method with different functionals overestimate the experimental excitation energies by on average ∼0.3 eV which is in accordance with previous findings <cit.>. While, BLYP exhibits the smallest MAE between these functionals, it has the largest deviation between the minimum and maximum value of the AE as shown in Fig. <ref>. This is also supported by the largest value of STDEV as given in Table S1. Among all studied methods herein, the largest MAE of 0.498 eV was observed for the EOM-CCSD method. However, unlike the TDDFT results, all of the EOM-CCSD excitation energy errors are closely clustered around the MAE. This is also evident from a very small STDEV of 0.015 eV and small deviation between the maximum and minimum AE. The ΔHF method exhibits a very large discrepancy from experiments with MAE of 0.496 eV, which is mainly due to lack of electron correlation. Inclusion of electron correlation effects by moving to a CCSD ansatz within the ΔCCSD method yields the best agreement with experiments for which MAE is only 0.054 eV. In passing we note that due to the non-unitary form of the ansatz, CCSD is unsuitable for implementation on quantum hardware.

Next, we discuss the performance of the quantum ΔUCCSD and ΔADAPT methods proposed in this work. As shown in Fig. <ref>, in case of ΔUCCSD, the MAE systematically improves with increase of the active space size such that the ΔUCCSD(8,8) method outperforms classical TDDFT and EOM-CCSD methods yielding MAE of 0.131 eV. Moreover, as shown in Table S1, the ΔUCCSD method provides nearly identical excitation energies to that of ΔCCSD with the same active space size. This indicates that in the large active space limit, ΔUCCSD will have very similar accuracy to ΔCCSD. In addition, due to the variational nature of ΔUCCSD, it is expected to outperform ΔCCSD in situations where CCSD fails due to its nonvariational nature of the resulting wave function <cit.>. Finally, the ΔADAPT method shows nearly identical performance compared to ΔUCCSD. We note that the reported ΔADAPT excitation energies are calculated as the energy difference between excited state and ground state wave functions that converged to within 1e-2 Hartree accuracy relative to the exact energy. Because of their unitary nature, both ΔUCCSD and ΔADAPT are implementable on a quantum device. However, the implementation of individual states for ΔUCCSD(6,6) and ΔUCCSD(8,8) would require 9200 and 39568 CNOTs, respectively, rendering this simulations practically impossible on current near-term quantum devices. As already mentioned, we note that these estimates assume an all-to-all connectivity map of physical qubits as implemented in . By contrast, ΔADAPT(6,6) and ΔADAPT(8,8) require on average (average between ground state and excited state ADAPT ansatz) merely 70 CNOTs and 583 CNOTs, respectively (for more information regarding the CNOT count of each individual state see Table S2). Importantly, the latter CNOT counts are well within contemporary quantum hardware capabilities and work on experimental hardware-based quantum simulations are ongoing in our laboratories. We would like to stress out that these active space sizes are tractable with the conventional multi-reference methods, however the simulations performed herein provides us with a valuable resources estimates and their scaling behavior as we increase the active space size, which is an important step for enabling simulations on the NISQ devices.



§ CONCLUSION

In this work, we have introduced the ΔUCCSD and ΔADAPT methods for the calculation of excitation energies on noisy near-term quantum devices within a VQE algorithm formulation. The proposed approaches are non-Aufbau methods for which we employ in the first step non-Aufbau Hartree-Fock electronic references optimized by means of an SCF procedure <cit.>. Electronic correlation effects are subsequently accounted for by UCCSD-VQE and ADAPT-VQE, respectively. To  recover a significant portion of the total correlation energy, we additionally make use of the FNO approach. Moreover, solvation effects are taken into account with the CPCM implicit model. The developed methods were then employed to calculate vertical excitation energies (S_0 →S_1) of six representative and chemically diverse BODIPY molecules with potential application in PDT. We compared computationally predicted excitation energies to corresponding experimentally determined values. The latter strikingly reveals that the proposed ΔUCCSD and ΔADAPT approaches  outperform the popular classical TDDFT and EOM-CCSD methods. We also demonstrate that ΔADAPT can be implemented on NISQ devices by requiring only a modest number of CNOTs that is well within contemporary quantum hardware capabilities. 

The accuracy of both methods can be further improved by systematically increasing the active space size, thus approaching a <0.1 eV accuracy of predicted vertical excitation energy compared to experiment. Additional enhancement of the method compared to experiments (especially conducted in polar solvents that form hydrogen bonds with a solute) can be achieved by taking into account the explicit solvent effects instead of the implicit one, as employed herein, albeit at a higher computational cost due to statistical sampling <cit.>. The use of spin-pure reference configurations instead of the unrestricted one could further improve the obtained results <cit.>, and the work along those lines is in progress within our team.

In view of the accuracy and simplicity of ΔADAPT, we believe that this method will become the method of choice for the calculation of molecular excitation energies on NISQ devices. Following the growing interest in the development of near-term quantum algorithms for excited state quantum chemistry applications, the ΔADAPT method will open new frontiers and guide an effective design of molecular photosensitizers equipped with the optimal properties for usage in photodynamic therapy <cit.> and photocatalysis <cit.>. Due to its variational nature, ΔADAPT will be particularly useful for NISQ simulations of strongly correlated transition metal-containing complexes that are at the core of the vast majority of  prominent photosensitizers <cit.>. 


Work on “Quantum Computing for Photon-Drug Interactions in Cancer Prevention and Treatment” is supported by Wellcome Leap as part of the Q4Bio Program.




§ AUTHOR CONTRIBUTIONS

F.P. proposed the methods, wrote the code and the first draft. F.P. and V.K. prepared the figures. F.P. and A.N. obtained, evaluated, and analyzed the data. All authors participated in editing the manuscript.



Supporting Information Available: The supporting information includes: experimental and calculated excitation energies for six BODIPY systems; number of two-qubit CNOT gates for ADAPT ansatz.



Conflict of interest

The authors declare no conflict of interest.

1


