Add support for Quantum ESPRESSO

[S-Erik]

Feature details

Creating a Hamiltonian and performing VQE calculations starting from a DFT calculation using the Quantum ESPRESSO package. This would also include calculating observables of interest like spin2, spinz or the dipole_moment. Further, active space calculations, like with active_space, would also be very useful, e.g. by using the frozen core approximation.

Implementation

Similar to the molecular_hamiltonian function I suggest a qe_hamiltonian function that loads the output of a DFT calculation done with Quantum ESPRESSO and returns a Hamiltonian object. The Quantum ESPRESSO output is, in its most basic form, just an XML-file that stores general information and a binary file that stores the wavefunction information as complex expansion coefficients. For the implementation, it sould be noted that Quantum ESPRESSO uses plane-waves to describe the Kohn-Sham orbitals instead of Gaussian functions like PySCF. Therefore, methods in PennyLane that rely on contracted Gaussians cannot be used in combination with a qe_hamiltonian. An example would be the dipole_moment which uses the moment_integral function which uses contracted Gaussian function for its calculations. Concrete implementations should be very similar to a related project involving qiskit-nature: Qiskit-Quantum-Espresso-Driver. There, the loading of the Quantum ESPRESSO output files (wfc.py) and the calculation of all relevant Hamiltonian matrix elements (calc_matrix_elements.py for one-electron matrix elments and eri_pair_densities.py for two-electron matrix elements) is already implemented. Defining an active space with the frozen core approximation can be implemented from that in a straightforward manner. For adapting observables like the dipole_moment to the plane-wave basis set, I am eager to hear your suggestions.

Now I present a first draft for the qe_hamiltonian function which loads the Quantum ESPRESSO output and generates a PennyLane Hamiltonian. The following code is based on the Qiskit-Quantum-Espresso-Driver.

def qe_hamiltonian(
    xml_out_file,
    wfc_file,
    mult=1,
    active_electrons=None,
    active_orbitals=None,
    wires=None,
):  # pylint:disable=too-many-arguments, too-many-statements
    r"""Generate the qubit Hamiltonian from a Quantum ESPRESSO DFT simulation.
    The active space is described with the frozen core approximation.
    """
    wires_map = None

    if wires:
        wires_new = qml.qchem.convert._process_wires(wires)
        wires_map = dict(zip(range(len(wires_new)), list(wires_new.labels)))

    wfc_obj = Wfc.from_file(file=wfc_file, output_xml=xml_out_file)

    # For now we only support spin unpolarized calculations, i.e. one wfc file with half the number of electrons
    num_electrons = wfc_obj.num_electrons * 2

    num_orbitals = wfc_obj.num_bands

    # define the active space
    # Same as in pennylane.qchem.openfermion_obs.molecular_hamiltonian
    core, active = qml.qchem.active_space(
        num_electrons, num_orbitals, mult, active_electrons, active_orbitals
    )

    with open(xml_out_file, "r", encoding="utf-8") as file:
        xml_dict = xmltodict.parse(file.read())
    _reference_energy = xml_dict["qes:espresso"]["output"]["total_energy"]["etot"]

    overlaps_ncpp = wfc_obj.get_overlaps()

    if not np.allclose(overlaps_ncpp, np.eye(overlaps_ncpp.shape[0])):
        warnings.warn("The wavefunctions are not orthonormal!")

    p = wfc_obj.k_plus_G  # momenta used for the plane-waves, shape (#waves, 3)
    _c_ip = wfc_obj.evc  # expansion coefficients of the Kohn-Sham orbitals in the plane-wave basis ,shape (#bands, #waves)

    # Load expansion coefficients and occupations of active and core orbitals
    _occupations_active, c_ip_active = wfc_obj.get_orbitals_by_index(active)
    occupations_core, c_ip_core = wfc_obj.get_orbitals_by_index(core)

    # Calculate matrix elements
    # Kinetic energy
    iTj_orbitals = calc_matrix_elements.iTj(p, c_ip_active)

    # Nuclear interaction
    iUj_orbitals = calc_matrix_elements.iUj(
        p, c_ip_active, wfc_obj.atoms, wfc_obj.cell_volume
    )

    # Calculate ERIs (two-electron matrix elements) via pair densities
    two_mo: np.ndarray = (
        eri_pair_densities.eri(c_ip=c_ip_active, b=wfc_obj.b, mill=wfc_obj.mill)
        / wfc_obj.cell_volume
    )

    # Calculate frozen core
    h_pq_core = (
        eri_pair_densities.get_frozen_core_pot(
            c_ip_core=c_ip_core,
            c_ip_active=c_ip_active,
            b=wfc_obj.b,
            mill=wfc_obj.mill,
        )
    ) / wfc_obj.cell_volume

    # Frozen core energy offset
    core_constant = eri_pair_densities.get_frozen_core_energy(
        p=p,
        c_ip_core=c_ip_core,
        b=wfc_obj.b,
        mill=wfc_obj.mill,
        cell_volume=wfc_obj.cell_volume,
        atoms=wfc_obj.atoms,
        occupations_core=occupations_core,
    )

    # Calculate nuclear repulsion
    nucl_repulsion = calc_matrix_elements.nuclear_repulsion_energy(
        c_ip_active, wfc_obj.atoms, wfc_obj.cell_volume
    )

    # Build Hamiltonian
    one_mo = iTj_orbitals - iUj_orbitals + h_pq_core + nucl_repulsion

    hf = qml.qchem.fermionic_observable(core_constant, one_mo, two_mo)

    if qml.operation.active_new_opmath():
        h_pl = qml.jordan_wigner(hf, wire_map=wires_map, tol=1.0e-10).simplify()

    else:
        h_pl = qml.jordan_wigner(
            hf, ps=True, wire_map=wires_map, tol=1.0e-10
        ).hamiltonian()
        h_pl = qml.pauli.utils.simplify(h_pl)

    return h_pl, len(h_pl.wires)

How important would you say this feature is?

1: Not important. Would be nice to have.

Additional information

In general, I would be interested if in the future PennyLane aims to support classical ab-initio DFT codes like Quantum ESPRESSO or VASP, which use plane-wave basis sets, besides PySCF and Psi4, which use Gaussian basis sets. Further it would be interesting to know which observables, that are implemented in Pennylane, depend on the used basis set, e.g. the dipole_moment seems to rely on the usage of Gaussian basis functions.

[CatalinaAlbornoz]

Hi @S-Erik , thank you for your suggestion! I've shared it with our team. I'll let you know their comments in the next few days.

[trbromley]

Hi @S-Erik, thanks for this detailed feature request. Having better PennyLane integration with Quantum ESPRESSO, and in general offering easier ways to create Hamiltonians/observables describing materials is definitely something we'd like to build out in PennyLane.

It's not yet on our near-term roadmap, but it would be valuable to keep discussing to understand better what would be useful to you and most impactful overall. If you're interested in developing something already it might make sense to make a standalone repo with the functionality, allowing us to then think more about how it could fit into PennyLane.