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

This study describes a hybrid quantum-classical computational approach for designing synthesizable deuterated $Alq_3$ emitters possessing desirable emission quantum efficiencies (QEs). This design process has been performed on the tris(8-hydroxyquinolinato) ligands typically bound to aluminum in $Alq_3$. It involves a multi-pronged approach which first utilizes classical quantum chemistry to predict the emission QEs of the $Alq_3$ ligands. These initial results were then used as a machine learning dataset for a factorization machine-based model which was applied to construct an Ising Hamiltonian to predict emission quantum efficiencies on a classical computer. We show that such a factorization machine-based approach can yield accurate property predictions for all 64 deuterated $Alq_3$ emitters with 13 training values. Moreover, another Ising Hamiltonian could be constructed by including synthetic constraints which could be used to perform optimizations on a quantum simulator and device using the variational quantum eigensolver (VQE) and quantum approximate optimization algorithm (QAOA) to discover a molecule possessing the optimal QE and synthetic cost. We observe that both VQE and QAOA calculations can predict the optimal molecule with greater than 0.95 probability on quantum simulators. These probabilities decrease to 0.83 and 0.075 for simulations with VQE and QAOA, respectively, on a quantum device, but these can be improved to 0.90 and 0.084 by mitigating readout error. Application of a binary search routine on quantum devices improves these results to a probability of 0.97 for simulations involving VQE and QAOA.


I. INTRODUCTION
Organic luminescent materials have attracted significant and increasing interest from academia and industry over the past decade, particularly for the fabrication of lightweight and flexible optoelectronic devices which can be used in a variety of applications including OLED displays [1,2], luminescent solar concentrators [3,4] and photocatalysts [5,6]. While possessing many potential advantages, several serious challenges must be overcome before organic luminescent materials can be commercialized.
One challenge involves the need to improve the emission quantum efficiency, Φ, defined as: where κ r and κ nr are the rate constants for all radiative and non-radiative relaxation processes [7]. Since such materials intrinsically suffer from a variety of non-radiative decay modes, one method to achieve a high emission quantum efficiency is to lower κ nr . Based on the Fermi golden rule of time-dependent perturbation theory, the κ nr between two different electronic states with equal * caoch@user.keio.ac.jp † gojones@us.ibm.com ‡ yamamoto@appi.keio.ac.jp spin multiplicity in the weak nonadiabatic coupling case is proportional to the product of an electronic term of the nonadiabatic coupling, and a vibrational term comprising the Franck-Condon (FC) factor defined by the square of the overlap between total vibrational wavefunctions of different electronic states.
Vibrational normal modes with high vibrational frequencies are required to fulfill the energy conservation law. For example, C-H stretching modes (∼3000 cm −1 ) dominate FC factors since they are widely spaced in energy and possess much smaller vibrational quantum numbers than other vibrational normal modes comprising lower electronic states. Therefore, introduction of heavier deuterium isotopes to replace hydrogen atoms in C-H bonds could effectively lower FC factors (κ nr ) since the vibrational frequencies of C-D stretching modes are on the order of 2200 cm −1 . Many experimental studies show that substitution of hydrogen by its deuterium isotope improve the luminescent performance of organic materials [8][9][10][11] and could potentially be used for optoelectronic applications such as OLED devices [8] and image probes [9].
To accelerate the design of deuterated organic emitters with high emission quantum efficiencies, it would be quite useful to predict the positional dependence of deuteron substitutions on the FC factor using ab initio methods, leading to efficiently deuterated molecules possessing high emission QEs. However, such simulations of FC factors could be computationally demanding since an emitter possessing n hydrogen atoms would require eval-uations of 2 n possible deuterated molecules. Moreover, as a practical consideration, increasing the amount of deuterium in molecules significantly increases the cost of chemical synthesis [8,9], thus deuterated molecules produced by this design should not only exhibit high quantum efficiency but must also be synthesizable. Previous investigations have shown that machine learning is orders of magnitude faster than ab initio calculation for predicting desired properties and synthesizability can also be integrated into a training model [12,13], therefore, here we propose the use of a machine learning model to accelerate the discovery of optimally deuterated OLED emitters [14].
Such a machine learning model, trained to predict molecular properties based on molecular structure, can in principle generate data that can be enhanced by hybrid quantum-classical variational quantum optimization algorithms (VQAs) such as the Quantum Approximate Optimization Algorithm (QAOA) [15] and the Variational Quantum Eigensolver (VQE) [16]. These algorithms are suitable for execution on current noisy quantum hardware since they execute a short, parameterized quantum circuit on the devices while parameters are optimized using a classical outer loop to enable the discovery of highquality solutions by the quantum circuit.
In this study, we develop a combined quantum chemistry, machine learning and quantum optimization approach to discover an optimal, synthesizable deuterated molecule of Alq 3 (a well-known OLED emitter [17,18]) possessing high quantum efficiency. Our approach involves a sequence initiated by the evaluation of the Franck-Condon factors of various deuterated Alq 3 derivatives with classical computational chemistry methods. This step is followed by the selection of a suitable prediction model via machine learning performed on classical architecture. The sequence terminates with the use of an Ising Hamiltonian derived from this machine learning approach in both constrained and unconstrained quantum optimization procedures involving the use of the VQE and QAOA algorithms to predict the optimal deuterated molecule.

A. Computational Methodology
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 op- timal 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.

B. Quantum Chemistry and Machine Learning Methodology
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 → 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 x (s) represented by 6 qubits. A dataset is prepared with all the x (s) and its corresponding target values of y (s) as shown in Table 1. The non-deuterated Alq 3 molecule (i.e., 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, x (s) = [011111] and x (s) = [100000]) of the selected 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]: feature vector x ( s) target vector y ( s) 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 2 i = 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 κ is the dimension of the factorization. Following previous prescriptions [25], κ 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 , v l and 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 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 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.

C. Construction of the Constrained Prediction Model
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 β 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 x (s) . The new dataset has a format similar to that shown in Table 1. The x (s) in the new dataset also represent all 64 deuterated molecules and the corresponding target values y (s) are generated by where n i is the number of deuterium atoms (corresponding to bit value 0) in the i-th of x (s) and n 0 is a specific integer value with 0 ≤ n 0 ≤ 6. Thus, the y (s) has a minimum of value 0 when n i = n 0 . Note that in order to make the values in QU BO unconstrained and QU BO penalty possess the same order of magnitude, Q ij and Q ii in QU BO 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 QU BO unconstrained has n 0 deuterium atoms, the value of β 0 in Eq (5) is set to 10.

D. Quantum Optimization Methodology
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 |ψ(θ) is a parameterized Ansätze with U (θ), the unitary operator of the quantum circuit and |ψ 0 , 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 θ as shown in Figure 3(a). The parameter θ is updated so that the mean energy (i.e., E(θ) shown in Eq (9)) decreases in each iteration.
Like VQE, the parameters (γ i ,β i (i = 1, ..., p)) are updated so that the mean energy (given by E(θ) = ψ p ( γ, β)|H(s)|ψ p ( γ, β) ) 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.

E. Application of the binary search algorithm in optimizations on the quantum device
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: 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 0 or 1 i and the probability of the highest bitstring in the original VQE and QAOA simulations to decrease, the values of p 0 or 1 i are larger than the probability of highest bitstring; 2) leaving out qubits with p 0 or 1 i > δ 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.

A. Quantum Chemistry and Machine Learning
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 re-sulted 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. 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 QU BO 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 QU BO 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 QU BO penalty terms.

B. Simulations with the statevector simulator
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 prob-  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 Anstze 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(cd). 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(ab). 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.

C. Simulations on ibm kawasaki
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 proba-bility 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).

IV. CONCLUSIONS
A combined quantum chemistry, machine learning and quantum optimization approach has been developed for the molecular design of high-efficiency, synthesizable deuterated Alq 3 emitters. Quantum chemistry calculations of FC factors provide a consistent picture of the emission quantum efficiency of deuterated Alq 3 emitters. Based in part upon these quantum chemistry studies, an unconstrained Ising Hamiltonian was constructed with a factorization machine (FM) predictor and its ability to reproduce the property of emission quantum efficiency has been validated. Moreover, a synthetic constraint was introduced based on the fact that the synthesizability of the Alq 3 emitters is inversely proportional to the number of deuterium atoms in deuterated Alq 3 emitters. A constrained Ising Hamiltonian was constructed by adding the additional synthetic constraint to the unconstrained Ising Hamiltonia. The ground state of the constrained Ising Hamiltonian results in a deuterated Alq 3 molecule which meets the desiderata of both emission quantum efficiency and synthesizability.
VQE and QAOA calculations with the constrained and unconstrained Ising Hamiltonians are performed on the statevector simulator and IBM quantum devices. Calculations performed with VQE and QAOA algorithms on the statevector simulator accurately converge to the ground state and successfully find the optimal deuterated Alq 3 molecule with probability larger than 0.95. These results show that the VQE and QAOA algorithms can be reliably used to generate lead deuterated Alq 3 molecules possessing desirable emission quantum efficiency and are synthesizable.
Results from VQE and QAOA calculations on IBM quantum devices show that noise significantly limits the accuracy of energy values of converged state. As a consequence, the bitstring probability corresponding to the optimal deuterated Alq 3 molecule decreases to 0.83 and 0.075 for VQE and QAOA calculations, respectively. The inclusion of readout error mitigation in the VQE calculations increases the probability by 0.07, while only a marginal improvement of the probability was observed in QAOA calculations. To obtain more reliable quantum optimization results from simulations on the quantum device, a new scheme which utilizes the binary search algorithm in the VQE and QAOA calculations was developed to systematically search for the correct output binary value on every qubit. Both VQE and QAOA calculations provide the optimal bitstring with a probability of 0.97 when this scheme is used.
Beyond Alq 3 , the application of the binary search algorithm results in reliable VQE and QAOA calculation accuracy on quantum device, highlighting opportunities for the use of quantum device in the discovery of deuterated OLED emitters such as novel iridium and platinum complexes. More broadly, leveraging a combined quantum chemistry, machine learning and quantum optimization approaches into a discovery workflow paves the way