qml.qchem

Overview

The quantum chemistry module provides the functionality to perform Hartree-Fock calculations and construct observables such as molecular Hamiltonians as well as dipole moment, spin and particle number observables.

This subpackage provides the functionality to perform quantum chemistry calculations.

Functions

active_space(electrons, orbitals[, mult, …])	Build the active space for a given number of active electrons and active orbitals.
atom_basis_data(name, atom)	Generate default basis set parameters for an atom.
attraction_integral(r, basis_a, basis_b)	Return a function that computes the nuclear attraction integral for two contracted Gaussian functions.
attraction_matrix(basis_functions, charges, r)	Return a function that computes the electron-nuclear attraction matrix for a given set of basis functions.
clifford(generators, paulixops)	Compute a Clifford operator from a set of generators and Pauli-X operators.
contracted_norm(l, alpha, a)	Compute the normalization constant for a contracted Gaussian function.
core_matrix(basis_functions, charges, r)	Return a function that computes the core matrix for a given set of basis functions.
decompose(hf_file[, mapping, core, active])	Decomposes the molecular Hamiltonian into a linear combination of Pauli operators using OpenFermion tools.
diff_hamiltonian(mol[, cutoff, core, active])	Return a function that computes the qubit hamiltonian.
dipole_integrals(mol[, core, active])	Return a function that computes the dipole moment integrals over the molecular orbitals.
dipole_moment(mol[, cutoff, core, active])	Return a function that computes the qubit dipole moment observable.
dipole_of(symbols, coordinates[, name, …])	Computes the electric dipole moment operator in the Pauli basis.
electron_integrals(mol[, core, active])	Return a function that computes the one- and two-electron integrals in the molecular orbital basis.
electron_repulsion(la, lb, lc, ld, ra, rb, …)	Compute the electron-electron repulsion integral between four primitive Gaussian functions.
excitations(electrons, orbitals[, delta_sz])	Generate single and double excitations from a Hartree-Fock reference state.
excitations_to_wires(singles, doubles[, wires])	Map the indices representing the single and double excitations generated with the function excitations() to the wires that the Unitary Coupled-Cluster (UCCSD) template will act on.
expansion(la, lb, ra, rb, alpha, beta, t)	Compute Hermite Gaussian expansion coefficients recursively for two Gaussian functions.
fermionic_dipole(mol[, cutoff, core, active])	Return a function that builds the fermionic dipole moment observable.
fermionic_hamiltonian(mol[, cutoff, core, …])	Return a function that computes the fermionic hamiltonian.
fermionic_observable(constant[, one, two, …])	Create a fermionic observable from molecular orbital integrals.
gaussian_kinetic(la, lb, ra, rb, alpha, beta)	Compute the kinetic integral for two primitive Gaussian functions.
gaussian_moment(li, lj, ri, rj, alpha, beta, …)	Compute the one-dimensional multipole moment integral for two primitive Gaussian functions.
gaussian_overlap(la, lb, ra, rb, alpha, beta)	Compute overlap integral for two primitive Gaussian functions.
hermite_moment(alpha, beta, t, order, r)	Compute the Hermite moment integral recursively.
hf_energy(mol)	Return a function that computes the Hartree-Fock energy.
hf_state(electrons, orbitals)	Generate the occupation-number vector representing the Hartree-Fock state.
import_operator(qubit_observable[, format, …])	Convert an external operator to a Pennylane operator.
jordan_wigner(op)	Convert a fermionic operator to a qubit operator using the Jordan-Wigner mapping.
kinetic_integral(basis_a, basis_b)	Return a function that computes the kinetic integral for two contracted Gaussian functions.
kinetic_matrix(basis_functions)	Return a function that computes the kinetic matrix for a given set of basis functions.
meanfield(symbols, coordinates[, name, …])	Generates a file from which the mean field electronic structure of the molecule can be retrieved.
mol_basis_data(name, symbols)	Generates default basis set parameters for a molecule.
mol_density_matrix(n_electron, c)	Compute the molecular density matrix.
molecular_hamiltonian(symbols, coordinates)	Generate the qubit Hamiltonian of a molecule.
moment_integral(basis_a, basis_b, order, idx)	Return a function that computes the multipole moment integral for two contracted Gaussians.
moment_matrix(basis_functions, order, idx)	Return a function that computes the multipole moment matrix for a set of basis functions.
nuclear_attraction(la, lb, ra, rb, alpha, …)	Compute nuclear attraction integral between primitive Gaussian functions.
nuclear_energy(charges, r)	Return a function that computes the nuclear-repulsion energy.
observable(fermion_ops[, init_term, …])	Builds the fermionic many-body observable whose expectation value can be measured in PennyLane.
one_particle(matrix_elements[, core, …])	Generates the FermionOperator representing a given one-particle operator required to build many-body qubit observables.
optimal_sector(qubit_op, generators, …)	Get the optimal sector which contains the ground state.
overlap_integral(basis_a, basis_b)	Return a function that computes the overlap integral for two contracted Gaussian functions.
overlap_matrix(basis_functions)	Return a function that computes the overlap matrix for a given set of basis functions.
particle_number(orbitals)	Compute the particle number observable \(\hat{N}=\sum_\alpha \hat{n}_\alpha\) in the Pauli basis.
paulix_ops(generators, num_qubits)	Generate the single qubit Pauli-X operators \(\sigma^{x}_{i}\) for each symmetry \(\tau_j\), such that it anti-commutes with \(\tau_j\) and commutes with all others symmetries \(\tau_{k\neq j}\).
primitive_norm(l, alpha)	Compute the normalization constant for a primitive Gaussian function.
qubit_observable(o_ferm[, cutoff])	Convert a fermionic observable to a PennyLane qubit observable.
read_structure(filepath[, outpath])	Read the structure of the polyatomic system from a file and returns a list with the symbols of the atoms in the molecule and a 1D array with their positions \([x_1, y_1, z_1, x_2, y_2, z_2, \dots]\) in atomic units (Bohr radius = 1).
repulsion_integral(basis_a, basis_b, …)	Return a function that computes the electron-electron repulsion integral for four contracted Gaussian functions.
repulsion_tensor(basis_functions)	Return a function that computes the electron repulsion tensor for a given set of basis functions.
scf(mol[, n_steps, tol])	Return a function that performs the self-consistent-field calculations.
simplify(h[, cutoff])	Add together identical terms in the Hamiltonian.
spin2(electrons, orbitals)	Compute the total spin observable \(\hat{S}^2\).
spinz(orbitals)	Computes the total spin projection observable \(\hat{S}_z\).
symmetry_generators(h)	Compute the generators \(\{\tau_1, \ldots, \tau_k\}\) for a Hamiltonian over the binary field \(\mathbb{Z}_2\).
taper(h, generators, paulixops, paulix_sector)	Transform a Hamiltonian with a Clifford operator and taper qubits.
taper_hf(generators, paulixops, …)	Transform a Hartree-Fock state with a Clifford operator and taper qubits.
two_particle(matrix_elements[, core, …])	Generates the FermionOperator representing a given two-particle operator required to build many-body qubit observables.

Classes

BasisFunction(l, alpha, coeff, r)	Create a basis function object.
Molecule(symbols, coordinates[, charge, …])	Create a molecule object that stores molecular information and default basis set parameters.

Differentiable Hartree-Fock

The differentiable Hartree-Fock (HF) solver allows performing differentiable HF calculations and computing exact gradients with respect to molecular geometry, basis set, and circuit parameters simultaneously using the techniques of automatic differentiation available in Autograd. This makes the solver more versatile and robust compared to non-differentiable tools that mainly rely on numerical approaches for computing gradients, which can lead to inaccuracies and instability. Additionally, optimizing the basis set parameters allows reaching lower ground-state energies without increasing the size of the basis set. Overall, the solver allows users to execute end-to-end differentiable algorithms for quantum chemistry.

The differentiable HF solver computes the integrals over basis functions, constructs the relevant matrices, and performs self-consistent-field iterations to obtain a set of optimized molecular orbital coefficients. These coefficients and the computed integrals over basis functions are used to construct the one- and two-body electron integrals in the molecular orbital basis, which can be used to generate differentiable second-quantized Hamiltonians and dipole moments in the fermionic and qubit basis.

The following code shows the construction of the Hamiltonian for the hydrogen molecule where the geometry of the molecule and the basis set parameters are all differentiable.

import pennylane as qml
from pennylane import numpy as np

symbols = ["H", "H"]
# This initial geometry is suboptimal and will be optimized by the algorithm
geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]], requires_grad=True)

# The exponents and contraction coefficients of the Gaussian basis functions
alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
                  [3.42525091, 0.62391373, 0.1688554]], requires_grad = True)
coeff = np.array([[0.15432897, 0.53532814, 0.44463454],
                  [0.15432897, 0.53532814, 0.44463454]], requires_grad = True)

We then construct the Hamiltonian.

args = [geometry, alpha, coeff] # initial values of the differentiable parameters

hamiltonian, qubits = qml.qchem.molecular_hamiltonian(symbols, geometry, alpha=alpha, coeff=coeff, args=args)

>>> print(hamiltonian)
  (-0.35968235922631075) [I0]
+ (-0.11496335836166222) [Z2]
+ (-0.11496335836166222) [Z3]
+ (0.13082414303722753) [Z1]
+ (0.13082414303722759) [Z0]
+ (0.1031689825163302) [Z0 Z2]
+ (0.1031689825163302) [Z1 Z3]
+ (0.15329959994281844) [Z0 Z3]
+ (0.15329959994281844) [Z1 Z2]
+ (0.15405495855815063) [Z0 Z1]
+ (0.1609686663987323) [Z2 Z3]
+ (-0.05013061742648825) [Y0 Y1 X2 X3]
+ (-0.05013061742648825) [X0 X1 Y2 Y3]
+ (0.05013061742648825) [Y0 X1 X2 Y3]
+ (0.05013061742648825) [X0 Y1 Y2 X3]

The generated Hamiltonian can be used in a circuit where the molecular geometry, the basis set parameters, and the circuit parameters are optimized simultaneously. Further information about molecular geometry optimization with PennyLane is provided in this paper and this demo.

dev = qml.device("default.qubit", wires=4)
hf_state = np.array([1, 1, 0, 0])
params = [np.array([0.0], requires_grad=True)] # initial values of the circuit parameters

def generate_circuit(mol):
    @qml.qnode(dev)
    def circuit(*args):
        qml.BasisState(hf_state, wires=[0, 1, 2, 3])
        qml.DoubleExcitation(*args[0][0], wires=[0, 1, 2, 3])
        return qml.expval(qml.qchem.molecular_hamiltonian(mol.symbols, mol.coordinates, alpha=mol.alpha, coeff=mol.coeff, args=args[1:])[0])
    return circuit

Now that the circuit is defined, we can create a geometry and parameter optimization loop. For convenience, we create a molecule object that stores the molecular parameters.

for n in range(21):

    mol = qml.qchem.Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
    args = [params, geometry, alpha, coeff] # initial values of the differentiable parameters

    # compute gradients with respect to the circuit parameters and update the parameters
    g_params = qml.grad(generate_circuit(mol), argnum = 0)(*args)
    params = params - 0.25 * g_params[0]

    # compute gradients with respect to the nuclear coordinates and update geometry
    g_coor = qml.grad(generate_circuit(mol), argnum = 1)(*args)
    geometry = geometry - 0.5 * g_coor

    # compute gradients with respect to the Gaussian exponents and update the exponents
    g_alpha = qml.grad(generate_circuit(mol), argnum = 2)(*args)
    alpha = alpha - 0.25 * g_alpha

    # compute gradients with respect to the Gaussian contraction coefficients and update them
    g_coeff = qml.grad(generate_circuit(mol), argnum = 3)(*args)
    coeff = coeff - 0.25 * g_coeff

    if n%5 == 0:
        print(f'Step: {n}, Energy: {generate_circuit(mol)(*args)}, Maximum Absolute Force: {abs(g_coor).max()}')

Running this optimization, we get the following output in atomic units:

Step: 0, Energy: -1.0491709019856188, Maximum Absolute Force: 0.1580194718925249
Step: 5, Energy: -1.1349862621177522, Maximum Absolute Force: 0.037660768852544046
Step: 10, Energy: -1.1399960666483346, Maximum Absolute Force: 0.005175323916673413
Step: 15, Energy: -1.140321384816611, Maximum Absolute Force: 0.0004138319900744425
Step: 20, Energy: -1.1403680839339787, Maximum Absolute Force: 8.223248376348913e-06

Note that the computed energy is lower than the Full-CI energy, -1.1373060483 Ha, obtained without optimizing the basis set parameters.

The components of the HF solver can also be differentiated individually. For instance, the overlap integral can be differentiated with respect to the basis set parameters as follows

symbols = ["H", "H"]
geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]], requires_grad=False)
alpha = np.array([[3.42525091, 0.62391373, 0.1688554],
                [3.42525091, 0.62391373, 0.1688554]], requires_grad = True)
coeff = np.array([[0.15432897, 0.53532814, 0.44463454],
                [0.15432897, 0.53532814, 0.44463454]], requires_grad = True)

mol = qml.qchem.Molecule(symbols, geometry, alpha=alpha, coeff=coeff)
args = [alpha, coeff]

a = mol.basis_set[0]
b = mol.basis_set[1]

g_alpha = qml.grad(qml.qchem.overlap_integral(a, b), argnum = 0)(*args)
g_coeff = qml.grad(qml.qchem.overlap_integral(a, b), argnum = 1)(*args)

>>> print(g_alpha)
[[ 0.00169332 -0.14826928 -0.37296693]
 [ 0.00169332 -0.14826928 -0.37296693]]

OpenFermion-PySCF backend

The molecular_hamiltonian() function can be also used to construct the molecular Hamiltonian with a non-differentiable backend that uses the OpenFermion-PySCF plugin interfaced with the electronic structure package PySCF. The non-differentiable backend can be selected by setting method='pyscf' in molecular_hamiltonian():

import pennylane as qml
from pennylane import numpy as np

symbols = ["H", "H"]
geometry = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 2.0]])
hamiltonian, qubits = qml.qchem.molecular_hamiltonian(
    symbols,
    geometry,
    charge=0,
    mult=1,
    basis='sto-3g',
    method='pyscf'
)

The non-differentiable backend requires the OpenFermion-PySCF plugin to be installed by the user with

pip install openfermionpyscf