Quantum Defect Embedding Theory

Mean field calculations

We implement QDET to compute the excitation energies of a negatively charged nitrogen-vacancy defect in diamond [4]. We use DFT to obtain a mean field description of the whole system. The DFT calculations are performed with the QUANTUM ESPRESSO package. This requires downloading pseudopotentials [5] for each atomic species in the system from the QUANTUM ESPRESSO database. To prepare our system, the necessary carbon and nitrogen pseudopotentials can be downloaded by executing the following commands through the terminal or command prompt:

wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/C_ONCV_PBE-1.2.upf
wget -N -q http://www.quantum-simulation.org/potentials/sg15_oncv/upf/N_ONCV_PBE-1.2.upf

Next, we need to create the input file for running QUANTUM ESPRESSO. This contains information about the system and the DFT calculations. More details on how to construct the input file can be found in the QUANTUM ESPRESSO documentation page. For the system taken here, the input file can be downloaded with

wget -N -q https://west-code.org/doc/training/nv_diamond_63/pw.in

DFT calculations can now be initiated using the pw.x executable in WEST, taking pw.in as the input file and directing the output to pw.out. This process is parallelized across 2 cores using mpirun.

mpirun -n 2 pw.x -i pw.in > pw.out

Identify the impurity

Once we have obtained the mean field description, we can identify our impurity by finding the states that are localized around the defect region in real space. To do that, we compute the localization factor \(L_n\) for each state n, defined as:

where \(V\) is the identified volume including the impurity within the supercell volume \(\Omega\) and \(\Psi\) is the wavefunction [4]. We will use the WEST program to compute the localization factor. This requires the westpp.in input file, example for which is shown below. Here, we specify the box parameters within which the localization factor is being computed; the vectors for this box are provided in in atomic units as [x_start, x_end, y_start, y_end, z_start, z_end].

The calculation can now be performed by running the westpp.x executable from WEST using mpirun to parallelize it across two cores.

mpirun -n 2 westpp.x -i westpp.in > westpp.out

This creates the file westpp.json which contains the information we need here. Since computational resources required to run the calculation are large, for the purpose of this tutorial we just download a pre-computed file with:

mkdir -p west.westpp.save
wget -N -q https://west-code.org/doc/training/nv_diamond_63/box_westpp.json -O west.westpp.save/westpp.json

We can plot the computed localization factor for each of the states:

import json
import numpy as np
import matplotlib.pyplot as plt

with open('west.westpp.save/westpp.json','r') as f:
    data = json.load(f)

y = np.array(data['output']['L']['K000001']['local_factor'],dtype='f8')
x = np.array([i+1 for i in range(y.shape[0])])

plt.plot(x,y,'o')
plt.axhline(y=0.08,linestyle='--',color='red')

plt.xlabel('Kohn-Sham orbital index')
plt.ylabel('Localization factor')

plt.show()

From this plot, it is easy to see that the orbitals can be categorized as orbitals with low and high localization factors. For the purpose of defining an impurity, we need highly localized orbitals, so we set a cutoff of \(0.06\) , illustrated by the red dashed line, and choose the orbitals that have a localization factor larger than \(0.06\). We’ll use these orbitals for the calculation of the parameters needed to construct the effective Hamiltonian.

Electronic Integrals

The next step in QDET is to define the effective one-body and two-body integrals for the impurity. The effective two-body integrals \(v^{eff}\) are computed first as matrix elements of the partially screened static Coulomb potential \(W_0^{R}\).

where \(W_0^R\) results from screening the bare Coulomb potential \(v\) with the reduced polarizability \(P_0^R = P - P_{imp}\), where \(P\) is the system’s polarizability and \(P_{imp}\) is the impurity’s polarizability. Since solving the effective Hamiltonian accounts for the exchange and correlation interactions between the active electrons, we remove these interactions from the Kohn-Sham Hamiltonian \(H^{KS}\) to avoid double counting them.

The one-body term \(t^{eff}\) is obtained by subtracting from the Kohn-Sham Hamiltonian the double-counting term accounting for electrostatic and exchange-correlation interactions in the active space.

We use the WEST program to compute these parameters. WEST will first compute the quasiparticle energies, then the partially screened Coulomb potential, and finally the parameters of the effective Hamiltonian. The software offers several execution modes through the wfreq_calculation keyword. We choose XWGQH to ensure the full workflow is executed, which computes both the quasiparticle corrections and the final parameters for the QDET effective Hamiltonian. The input file for such a calculation is shown below:

We can now execute the calculation with:

mpirun -n 2 wfreq.x -i wfreq.in > wfreq.out

This calculation takes some time and requires computational resources, therefore we download a pre-computed output file with

mkdir -p west.wfreq.save
wget -N -q https://west-code.org/doc/training/nv_diamond_63/wfreq.json -O west.wfreq.save/wfreq.json

This output file contains all the information we need to construct the effective Hamiltonian.

Effective Hamiltonian

We now construct the effective Hamiltonian by importing the electron integral results and using WEST:

from westpy.qdet import QDETResult

effective_hamiltonian = QDETResult(filename='west.wfreq.save/wfreq.json')

The effective Hamiltonian can be solved using a high level method such as the full configuration interaction (FCI) algorithm from WEST as:

solution = effective_hamiltonian.solve()

Using solve() prints the excitation energies, spin multiplicity and relative occupation of the active orbitals.

The solution object is a dictionary that contains information about the FCI eigenstates of the system, which includes various excitation energies, spin multiplicities, eigenvectors etc. More importantly, while FCI handles small embedded effective Hamiltonians with ease, it quickly hits a wall with larger impurities. This is precisely where quantum computing steps in, offering the scalability needed to tackle such complex systems. The first step to solving these effective Hamiltonians via quantum algorithms in PennyLane, is to convert them to qubit Hamiltonians.

Quantum Simulation

We now map the effective Hamiltonian to the qubit basis. Note that the two-electron integrals obtained before are represented in chemists’ notation and need to be converted to the physicists’ notation for compatibility with PennyLane. Here’s how to construct the qubit Hamiltonian:

from pennylane.qchem import one_particle, two_particle, observable
import numpy as np

effective_hamiltonian = QDETResult(filename="west.wfreq.save/wfreq.json")

one_e, two_e = effective_hamiltonian.h1e, effective_hamiltonian.eri

t = one_particle(one_e[0])
v = two_particle(np.swapaxes(two_e[0][0], 1, 3))
qubit_op = observable([t, v], mapping="jordan_wigner")
print("Qubit Hamiltonian: ", qubit_op)

We can use this Hamiltonian in a quantum algorithm such as quantum phase estimation (QPE). You can compare the results and verify that the computed energies from quantum algorithm match those that we obtained before.