Quantum-Classical Computational Molecular Design of
Deuterated High-Efficiency OLED Emitters

$Alq_{3}$ (tris(8-hydroxyquinolinato)aluminium) is a stable metal chelate complex wherein aluminum is bonded to three 8-hydroxyquinoline ligands, as shown in Figure 1. For every ligand, there are six hydrogen atoms which can be replaced by deuterium isotopes [8]. By assuming that every ligand has the same deuterated structure and setting every hydrogen atom and its deuterium as bit values of 1 and 0 respectively, the combinatorial optimization required for discovering a desired deuterated $Alq_{3}$ molecule with high emission quantum efficiency can be viewed as a problem involving searching for an optimal bitstring represented by 6 qubits on a quantum computer.

Figure 2 shows a flowchart illustrating the hybrid quantum-classical procedure combining quantum chemistry, machine learning and quantum optimization to identify deuterated $Alq_{3}$ molecules possessing desired quantum efficiency. This procedure begins with quantum chemistry calculations with a classical computer to obtain the quantum efficiencies of a set of deuterated $Alq_{3}$ molecules. Results obtained from these calculations are used to prepare a training dataset used to train a model to predict the quantum efficiency, and a test dataset used to validate the performance of the trained model by comparing predicted values with those obtained using the test dataset. This procedure is repeated with additional quantum chemistry results added to the training dataset until the predicted accuracy of the trained model achieves a desired threshold. The improved model is used to construct an Ising Hamiltonian which is used to perform unconstrained and constrained optimization using a quantum computer to find the optimally deuterated $Alq_{3}$ derivatives. The following sections describe the details of the different aspects of this procedure involving quantum chemistry calculations, machine learning, and quantum optimizations on quantum devices and simulators.

In order to calculate the FC factors for nonradiative decay from $S_{1}$ to $S_{0}$ of the variously substituted deuterated $Alq_{3}$ molecules, the calculations of harmonic vibrational frequencies and normal modes of the structures with geometries optimized on the $S_{0}$ and $S_{1}$ electronic states are performed with the CAM-B3LYP [19] and (TD)CAM-B3LYP methods [20] and the 6-31G(d,p) basis sets using Gaussian09 [21]. The FC factors for the $S_{1}\rightarrow S_{0}$ nonradiative decay of non-deuterated and deuterated $Alq_{3}$ derivatives were calculated from the $S_{0}$ and $S_{1}$ vibrational frequencies and normal modes using MOMAP version 0.3.001 [22]. These calculations consist of the square of the overlap between the vibrational wavefunctions of the $S_{1}$ and $S_{0}$ states with the adiabatic excitation energy of the $S_{1}$ excited state above the $S_{0}$ vibrational ground state at $300K$.

Using the numbering scheme for the hydrogen atoms in the ligand as shown in Figure 1, every deuterated molecular structure is transformed into a binary feature vector $\vec{x}^{(s)}$ represented by 6 qubits. A dataset is prepared with all the $\vec{x}^{(s)}$ and its corresponding target values of $\vec{y}^{(s)}$ as shown in Table 1. The non-deuterated $Alq_{3}$ molecule (i.e., $\vec{x}^{(s)}=[111111]$) was chosen and a singly deuterated structure was randomly selected as the training data. In addition to this selection, deuterated structures possessing the normal reverse feature vectors (for example, $\vec{x}^{(s)}=[011111]$ and $\vec{x}^{(s)}=[100000]$) of the selected $\vec{x}^{(s)}$ were also used in the training data. The selected training data are used to train a prediction model described by the quadratic unconstrained binary optimization (QUBO) formulation [23]:

where $q_{i}$ is a binary variable which takes either 0 or 1, and coefficients $Q_{i}j$ and $Q_{i}i$ have real values. Here, the second equality holds because $q_{i}^{2}=q_{i}$. The $Q_{i}j$ and $Q_{i}i$ are obtained through the learning process with a Factorization Machine (FM) predictor [24] with model equation:

where $\kappa$ is the dimension of the factorization. Following previous prescriptions [25], $\kappa$ is set to 8 during the learning process to minimize the cost function.

Note that although the unconstrained QUBO model can also be determined by the regression of $Q_{ij}$ and $Q_{ii}$ in equation (1), the advantage of using the FM predictor is that it can reliably fit the parameters (i.e. $w_{i}$, $\vec{v}_{l}$ and $\vec{v}_{j}$ with a small training dataset and prevent the trained model from overfitting as previously shown [25]. This prediction model is termed the unconstrained QUBO model. For a thorough discussion on the learning ability of FM, we refer the reader to previous studies [24, 26, 27].

The performance of the unconstrained QUBO model is evaluated by checking the correlation (i.e., the square of the sample correlation coefficient denoted by $R^{2}$) between the FC values of all deuterated molecular structures not selected by the unconstrained QUBO model and the quantum chemistry calculations. If $R^{2}$ is less than a threshold value of 0.95, a randomly chosen $\vec{x}^{(s)}$ possessing one deuterium atom and its reverse feature vector is added to the training data to retrain the unconstrained QUBO model. The procedure repeats by choosing $\vec{x}^{(s)}$ while gradually increasing the number of deuterium atoms (i.e., the bit value of 0) in the feature vector, until $R^{2}$ surpasses the threshold value.

After creating an unconstrained QUBO model to accurately predict the FC values, we then build a new constrained QUBO model which includes a penalty term added to the unconstrained QUBO model:

where $\beta_{0}$ is a weighting parameter . The penalty QUBO model is also obtained using the FM predictor with a new dataset corresponding to the feature vector $\vec{x}^{(s)}$. The new dataset has a format similar to that shown in Table 1. The $\vec{x}^{(s)}$ in the new dataset also represent all 64 deuterated molecules and the corresponding target values $\vec{y}^{(s)}$ are generated by

where $n_{i}$ is the number of deuterium atoms (corresponding to bit value 0) in the i-th of $\vec{x}^{(s)}$ and $n_{0}$ is a specific integer value with $0\leq n_{0}\leq 6$. Thus, the $\vec{y}^{(s)}$ has a minimum of value 0 when $n_{i}=n_{0}$. Note that in order to make the values in $QUBO_{unconstrained}$ and $QUBO_{penalty}$ possess the same order of magnitude, $Q_{ij}$ and $Q_{ii}$ in $QUBO_{unconstrained}$ are scaled by dividing the minimum value of $Q_{ij}$ and $Q_{ii}$ into all $Q_{ij}$ and $Q_{ii}$. Note that in order to guarantee the minimum of $QUBO_{unconstrained}$ has $n_{0}$ deuterium atoms, the value of $\beta_{0}$ in Eq (5) is set to 10.

To perform quantum optimizations on quantum simulators and devices using the VQE and QAOA algorithms, both unconstrained and constrained QUBO models are converted into a spin-based Ising model Hamiltonian (denoted $s_{i}$) which takes a value of +1 or -1 [28]:

where $s_{i}$ is the Pauli Z operator acting on the i-th qubit, and the relation between $q_{i}$ in eq (1) and $s_{i}$ is given by Equation 8,

which indicates how the classical parameter $q_{i}$ is replaced by the quantum parameter $s_{i}$.

In the VQE method, the ground state energy of the $H(s)$ is calculated by minimizing the mean energy

where $|\psi(\theta)\rangle$ is a parameterized Ansätze

with $U(\theta)$, the unitary operator of the quantum circuit and $|\psi_{0}\rangle$, the initial state.

In the QAOA method, the parameterized Ansätze state is constructed from $H(s)$ with single qubit Pauli $X$ operators as:

where $p$ is an integer parameter representing the number of repetitions of the unitary operators.

The $R_{y}$ heuristic Ansätze, which is suitable for calculations on quantum devices, was used in the VQE simulations with a circuit depth equal to 1. The 6-qubit circuit featuring the $R_{y}$ Ansätze has 5 nearest-neighbor CNOT gates and 12 optimization parameters $\theta$ as shown in Figure 3(a). The parameter $\theta$ is updated so that the mean energy (i.e., $E(\theta)$ shown in Eq (9)) decreases in each iteration.

Like VQE, the parameters ($\gamma_{i}$,$\beta_{i}(i=1,...,p)$) are updated so that the mean energy (given by $E(\theta)=\langle{\psi_{p}(\vec{\gamma},\vec{\beta)}|H(s)|\psi_{p}(\vec{\gamma},\vec{\beta)}}\rangle$) decreases in each iteration. An example QAOA circuit with $p=3$ is shown in Figure 3(b) which requires 90 CNOT gates and 6 optimization parameters. The advantage of using the QAOA Ansätze is that adiabatic dynamics can be used to effectively find the ground state of the Hamiltonian to solve the combinatorial optimization problem [15, 29].

All quantum optimization calculations were performed with the statevector simulator contained in the Aer module of Qiskit 0.25 [30], as well as on the $ibm\_kawasaki$ 27-qubit quantum device. The statevector simulator uses linear algebra operations to compute the Ansätze state and expectation values exactly. The Constrained Optimization BY Linear Approximation (COBYLA) optimizer [31] was used to update parameters for calculations on both simulator and quantum device.

The readout error mitigation routine in Qiskit was used to mitigate errors occurring when computing ground state energies. In detail, the 64 computational basis states of the 6-qubit system are measured for the purpose of readout error mitigation, then the vectors of probability for the 64 basis states are created for every measurement, and the combination of all 64 vectors resulted in the formation of a $64x64$ matrix from which a calibration matrix was generated by the least-squares fitting. The calibration matrix was applied to the measurements obtained from VQE and QAOA calculations to correct the mean energy value. The parameter values in the VQE algorithm with the $R_{y}$ Ansätze and QAOA are then calculated using the COBYLA optimizer. Six linearly connected qubits were selected to reduce the influence of CNOT errors on the accuracy of results computed on the quantum device. Figure 4 show the architecture of the $ibm\_kawasaki$ device and the quantum circuits comprising the 6 chosen qubits (i.e. [q0,q1,q2,q3,q5,q8]) implementing VQE with the $R_{y}$ Ansätze and QAOA with $p=3$.

Since one of the computational basis states must comprise the optimal bitstring, we can view searching for the optimal bitstring as a problem involving determining the binary values of every qubit. This inspired us to suggest a new approach, illustrated in Figure 5, which applies the binary search algorithm to VQE and QAOA calculations to improve the accuracy on a quantum device. In this approach, VQE and QAOA calculations are first performed on all the 6 qubits with a loose optimization convergence criterion. Using the computed results, the probability of the output of 0 or 1 for every qubit is then calculated as:

For qubits with $p_{i}^{\textrm{0 or 1}}$ larger than a threshold value denoted by $\delta$, we determine the binary values of the qubits and construct a new Ising model Hamiltonian by omitting the terms of these related qubits. With the new Ising model Hamiltonian in hand, VQE and QAOA calculations are repeated for the remaining qubits with $p_{i}^{\textrm{0 or 1}}\leq\delta$ until the binary values of all the 6 qubits are determined. We expect that this approach is more effective for finding the optimal bitstring on a quantum device than the original VQE and QAOA approach due to the following reasons: 1) although noise from the quantum device causes both the accuracy of $p_{i}^{\textrm{0 or 1}}$ and the probability of the highest bitstring in the original VQE and QAOA simulations to decrease, the values of $p_{i}^{\textrm{0 or 1}}$ are larger than the probability of highest bitstring; 2) leaving out qubits with $p_{i}^{\textrm{0 or 1}}>\delta$ causes the number of qubits required for VQE and QAOA calculations to be reduced and consequently the binary values of some qubits are determined with smaller-depth quantum circuits than the original circuits for VQE and QAOA calculations.

The ability of the model obtained from the trained FM method to accurately predict FC factors was determined by comparing results predicted by the unconstrained QUBO model with results obtained from ab initio quantum chemistry calculations, as shown in Figure 6. With the smallest training dataset containing 3 deuterated molecules, the unconstrained QUBO model yields FC factors that are 100-1000 larger than results obtained from quantum chemistry results and the correlation between FC factors obtained with the unconstrained QUBO model and the ab initio prediction is 0.032. Increasing the number of deuterated molecules in the training dataset to 11, improves the order of the predicted FC values from the unconstrained QUBO model but the predicted FC values are still approximately an order of magnitude higher than quantum chemistry predictions. Moreover, the use of this larger training dataset results in only slightly better correlation between the QUBO model and ab initio prediction with an $R^{2}$ of 0.11. Gratifyingly, the use of 13 deuterated molecules as a training dataset resulted in excellent agreement between the unconstrained QUBO model predictions and the quantum chemistry calculations and improved the $R^{2}$ to 0.99. These results signify that machine learning can accurately predict the FC factors of $2^{n}$ deuterated molecules from training data derived from order n quantum chemistry calculations.

These 13 deuterated $Alq_{3}$ derivatives comprise the non-deuterated molecules with feature vector $\vec{x}^{(s)}=[111111]$ as well as the deuterated molecules containing only one deuterium atom or one hydrogen atom i.e., the one-hot ([000001], [000010], etc.) and one-cold ([111110], [111101], etc.) feature vectors. Because training an unconstrained QUBO model can be viewed as determining the values of entries in a 6x6 matrix, the values of 6 diagonal entries in the matrix are mainly determined by the 6 one-hot vectors in $\vec{x}^{(s)}$. Similarly, the 6 one-cold vectors and [111111] in $\vec{x}^{(s)}$ are used to determine the values of off-diagonal entries. This means that exclusion of some deuterated molecules in the training data will cause the learning process to be more difficult when trying to accurately determine the values in the matrix and emphasizes the need for an accurate unconstrained QUBO model.

Choosing a constrained QUBO model to find an optimally deuterated molecule possessing high desired quantum efficiency requires the constrained QUBO model to be in good agreement with quantum chemistry predictions of FC factors for all deuterated structures. Moreover, the deuterated structure corresponding to the global minimum of the constrained QUBO model must contain the expected number of deuterium atoms.

Figure 7 shows a representative plot of values obtained for an unconstrained QUBO model in comparison with a constrained QUBO model built using equations (5) and (6) with $n_{0}=3$ for the set of deuterated molecules. These results show that the $QUBO_{penalty}$ term has no effect on deuterated molecules possessing 3 deuterium atoms i.e., the values predicted by unconstrained and constrained QUBO models are the same for these molecules and are in good agreement with the FC factors obtained from quantum chemistry calculations. Consequently, the constrained QUBO model can reliably search the global minimum for $Alq_{3}$ derivatives possessing 3 deuterium atoms. On the other hand, for deuterated $Alq_{3}$ derivatives with more or less than 3 deuterium atoms, inclusion of the $QUBO_{penalty}$ term in the constrained QUBO model predicts values that are 40% to 860% larger than the values of the unconstrained QUBO model. Similar results were also observed for other constrained QUBO models penalized with $n_{0}=2,4,or5$. These results indicate that constrained QUBO models can be used to evaluate additional desirable traits, such as stability and low toxicity, by including adequate $QUBO_{penalty}$ terms.

The search for an optimal deuterated molecule can be mapped to a quantum optimization problem whose ground state corresponds to the optimal deuterated molecule by converting the unconstrained and constrained QUBO models described in the preceding section to the unconstrained and constrained Ising Hamiltonian models using equations (7) and (8). The ground state energy and the corresponding bitstrings appearing with highest probability, calculated by the exact eigensolver (i.e., an exact diagonalization of the Ising model Hamiltonian), VQE and QAOA algorithms with the statevector simulator for the unconstrained and constrained Ising models are shown in Table 2.

VQE with the heuristic $R_{y}$ Ansätze reproduces the ground state calculated with the exact eigensolver using an unconstrained Ising model Hamiltonian. The [000000] bitstring, corresponding to the fully deuterated $Alq_{3}$ ligand, is obtained with probability of 1.00, indicating that quantum optimization with VQE successfully finds the exact solution for the unconstrained Ising model Hamiltonian. While QAOA converges to bitstrings corresponding to the fully deuterated $Alq_{3}$ derivative, the probabilities for this selection only range from 0.44 to 0.97 increasing with the values of $p$ and not quite matching the probability of 1.00 observed with VQE. Similarly, the converged energy obtained from QAOA calculations decreases with increasing values of $p$. These improved behaviors are presumably caused by the increased accuracy of the description of the time-evolution operator with the increasing integer value of the parameters of $p$ [15].

Results produced with the constrained Ising model Hamiltonian are largely like those found with the unconstrained Ising model Hamiltonian: VQE produces the same accuracy and QAOA results also improve with increasing values of $p$. However, it was observed that to have a result comparable with the exact solution, $p=4$ was needed to converge the ground state with a probability of 0.95 with QAOA. Such behavior may be due to the energy differences between the unconstrained and constrained Ising model Hamiltonians shown in Figure 8. While 4 excited states are obtained for the unconstrained Ising model Hamiltonian, 19 excited states are obtained for the constrained Ising model Hamiltonian, and the gaps between the ground state and the first excited state for the unconstrained and constrained Ising model Hamiltonians are 1.00 and 0.40, respectively. Because more excited states and a smaller energy band gap are obtained for the constrained Ising model Hamiltonian, the optimization is more sensitive to the influence of the excited states, and it is more difficult to obtain the ground state for the QAOA calculation using the constrained Ising model Hamiltonian with $p=3$ than a similar calculation using the unconstrained Ising model Hamiltonian. Nevertheless, since the most probable bitstring from all the QAOA calculation is the same as the optimal bitstring obtained via the exact eigensolver, these results can be useful when searching for find the optimal deuterated molecule.

Table 2 also shows the number of CNOT gates in the quantum circuits and the number of optimization parameters used in the VQE and QAOA algorithms. This data was used to determine whether the optimization algorithms could be implemented on the $ibm\_kawasaki$ device for simulations involving 6 qubits. The quantum circuits for VQE calculations involving the use of the $R_{y}$ Ansätze possess 5 CNOT gates, whereas 30 CNOT gates are required for QAOA simulations with $p=1$ and this linearly increases with the value of $p$. Note that since the $R_{y}$ heuristic Ansätze employs nearest-neighbor coupling between the qubits in the so-called light-cone structure, by using the linearly connected 6-qubit subsection of $ibm\_kawasaki$, the number of CNOT gates would not change for when the circuit is implemented. In contrast, since QAOA does not have such a convenient coupling structure, an overhead of 3x 4x more qubits would be needed for implementing the circuit on the quantum device. Since the CNOT error rate can negatively affect the accuracy of the computed results, simulation with the VQE Ansätze on $ibm\_kawasaki$ is expected to provide a higher accuracy than provided by QAOA.

VQE requires 12 optimization parameters, while only 2 parameters are needed for QAOA with $p=1$; this value linearly increases with the value of $p$, and thus even with $p=4$ the number of optimization parameters is 8 which is still smaller than the number of parameters required for VQE. QAOA is expected to search these parameters more capably on a classical computer than VQE since a smaller number of optimization parameters infers better trainability and, thereby, a reduction in the number of iterations required for convergence of the energy.

The energy minimization and the success probability of the converged state in the VQE and QAOA calculations using the constrained Ising Hamiltonian with and without the readout error mitigations are shown in Figure 9(c-d). The deviation of ground state energy and the highest bitstring probability in VQE calculations on the quantum device from predictions obtained with the quantum simulator are very similar to VQE calculations using the unconstrained Ising Hamiltonian as shown in Figure 9(a-b). On the other hand, the deviations observed for QAOA calculations using the constrained Ising Hamiltonian are much larger than those with the unconstrained Ising Hamiltonian. For example, the energy deviations from the QAOA calculations using the constrained and unconstrained Ising Hamiltonian with readout error mitigation are 10.18 and 4.64, respectively, and their corresponding highest probabilities are 0.08 and 0.14. The lower accuracy of QAOA with the constrained Ising Hamiltonian may be due to the energy differences of the Hamiltonians as shown in Figure 8 and discussed in the preceding section. Overall, that fact that VQE can predict the optimal bitstring with probability larger than 0.9 on the quantum device is quite satisfactory. Although the success probability obtained from QAOA calculation is one order lower than those from VQE, the predicted results can still be useful for finding the optimally deuterated $Alq_{3}$ derivative.

The five highest probability bitstrings were also investigated for VQE and QAOA calculations and are shown in Figures 10(b,d). Bitstrings with probabilities ranging from the 2nd highest to 5th highest occurrences mostly differ from the dominant bitstring by only a single qubit. For example, as shown in Figure 9(b) for VQE calculation with the unconstrained Ising Hamiltonian the highest probability bitstring comprises the [000000] configuration and the next four bitstrings comprise the [000010], [001000], [010000], [100000] configurations. Consequently, the probability of the correct value 0 on qubits q1, q2, q3, q4 and q6 is above 0.96. The value is much higher than the probability of the dominant bitstring (i.e. [000000]) which is about 0.90. Moreover, we also found that the probability of the correct output on every qubit is quite different. For example, among the top five highest probability bitstrings, the probability of the correct value $0$ on qubit q5 is about 0.90 which is much lower than the probability on qubits q1, q2, q3, q4 and q6. This is a consequence of the fact that the impact of noise on each of these qubits is variable.

Simulations using VQE and QAOA with $p=3$ using the unconstrained Ising model Hamiltonian were performed on $ibm\_kawasaki$. The iteration steps (i.e., the number of updates of the parameters) versus the energy values, as well as the five greatest probability bitstrings of the converged state, are shown in Figure 9(a-b). To mitigate the measurement errors from the quantum device, we have applied readout error mitigation to ground states computed with VQE and QAOA. VQE without readout error mitigation yields energies that are about 0.43 higher than the results from the exact values calculated on the simulator. The bitstring appearing with the highest probability is the same as that obtained on the simulator but the probability is lowered by 0.18. The difference of the energy and the success probability are improved to about 0.21 and 0.10, which highlights the importance of the inclusion of readout error mitigation in VQE calculations on the quantum device. Although the most probable bitstring obtained via QAOA is the same as that predicted by VQE, the deviations of the success probability and the converged state from those obtained with QAOA on the simulator are approximately an order of magnitude larger than those obtained with VQE. This is caused by the fact that QAOA involves many more gate operations than VQE as shown in Figure 4. For example, there are 5 and 156 CNOT gates, which significantly limits the accuracy of quantum algorithms, in the VQE and QAOA quantum circuits, respectively. We also found that the application of readout error mitigation in QAOA results only in marginal improvement of the ground state energy and the value of highest probability of bitstring. These results suggest that errors such as decoherence noise is the main cause for the deviation of QAOA results from the ideal.

The preceding results inspired us to use the binary search algorithm in the VQE and QAOA calculations to reduce the effect of noise from the quantum device by systematically searching the correct output binary value on every qubit as shown in Figure 5, instead of directly searching for the optimal bitstring. The results of VQE and QAOA calculations using the binary search for the unconstrained Ising Hamiltonian are shown in Figure 10(a-b).

Applying the binary search algorithm to determine binary values of target qubits improves the accuracy of ground state energies compared to both the original VQE and QAOA calculations. The original VQE and QAOA calculations using the unconstrained Ising Hamiltonian yields energies that are about 0.21 and 4.64 higher than the exact values; this gap can be improved to about 0.08 and 0.01 by using the binary search algorithm. Notably, the improvement of the ground state prediction observed for QAOA is more significant than observed for VQE. This is due to the fact that the decreased number of target qubits causes the number of gate operations to decrease by a factor of 10-100 in the QAOA calculations whereas with VQE this decrease is only by a factor of 2-3. However, it must be noted that many more iterations are required to converge QAOA combined with binary search in comparison to VQE. Because the QAOA calculation for 6 qubits is much more sensitive to the effect of noise from the quantum device than the VQE calculation, the first application of binary search in the QAOA calculation can only leave 1 qubit out of 6 qubits whereas 3 qubits can be leaved out in the VQE calculations. Consequently, the VQE calculation can converge to the ground state using two binary searches, whereas 4 binary searches are needed for the QAOA calculation, as shown in Figure 10(a). The improvement of the values of the highest probability for the VQE and QAOA calculations with the binary search is also observed. Without the binary search, the values of the highest probability are 0.90 and 0.14 for the original VQE and QAOA calculations, respectively. Calculations using the binary search improves these values to 0.97 and 0.99 respectively, signifying the effectiveness of using the binary search to find the optimal bitstring on quantum device. Moreover, we note that the value of the highest probability from the QAOA calculations using the binary search on the quantum device is somewhat higher than the values obtained on the statevector simulator, as shown in Figure 10(b) and Table 2. This is caused by the fact that the number of terms in the Ising Hamiltonian decreases when the binary search algorithm is applied to QAOA simulations. Consequently, the energy eigenstates have smaller variance, and the converged state becomes closer to the exact ground state. Similar results are also found in the VQE and QAOA results using the constrained Ising Hamiltonian as shown in Figure 10(c-d).