The UCCSD Wavefunction ansatz¶

In quantum chemistry, the Unitary Coupled Cluster with Single and Double excitations (UCCSD) wavefunction ansatz 1, 2 is a powerful method used to approximate the ground state of a quantum system. This approach is particularly valuable in the context of quantum computing, where it is employed to solve the electronic structure problem. In this tutorial, we will learn how to implement the UCCSD in CUDA-Q to calclate the ground state energy for chemical systems.

[ ]:

# Requires pyscf to be installed
%pip install pyscf

[5]:

import cudaq

# Set the traget

# Double precision.
#cudaq.set_target("nvidia", option = "fp64")

# Single precision.
cudaq.set_target("nvidia")

We begin by calculating the electronic hamiltonian and converting it to qubit hamiltonian. To learn more, see this tutorial

[6]:

from qchem.classical_pyscf import get_mol_hamiltonian
from qchem.hamiltonian import jordan_wigner_fermion

#geometry = 'H 0.0 0.0 0.0; H 0.0 0.0 0.7474'
#geometry = 'Li 0.3925 0.0 0.0; H -1.1774 0.0 0.0'
geometry =  'Zn 0.0 0.0 00; O 0.0 0.0 1.7047'

# When not using active space.
#molecular_data = get_mol_hamiltonian(xyz=geometry, spin=0, charge=0, basis='sto3g', ccsd=True, verbose=True)

# When using active space.
molecular_data = get_mol_hamiltonian(xyz=geometry, spin=0, charge=0, basis='ccpvdz', nele_cas=4, norb_cas=4,
                                     ccsd=True, casci=True, verbose=True)

obi = molecular_data[0]
tbi = molecular_data[1]
const = molecular_data[2]
electron_count = molecular_data[3]
norbitals = molecular_data[4]
qubits_num = 2 * norbitals

# Here, we are excluding the const from the jordan wigner transformation.
# When calculating the expectation value of the Hamiltonian,
# it is recommend to remove the identity operator from the
# Hamiltonian and add its coefficient as a constant shift to the energy.
# In exact arithmetic, this should not make a difference because the expectation value
# for the identity is just the norm, which is 1.
# However, in floating-point arithmetic, this avoids the round-off error
# associated with that term.  This error can be significant because the
# coefficient in this term corresponds to core energy, which is usually the largest
# coefficient in the Hamiltonian by an order of magnitude.
hamiltonian = jordan_wigner_fermion(obi, tbi, 0.0, tolerance = 1e-15)

overwrite output file: Zn 0-pyscf.log
[pyscf] Total number of orbitals =  57
[pyscf] Total number of electrons =  38
[pyscf] HF energy =  -1852.5618486460523
[pyscf] R-CASCI energy using molecular orbitals=  -1852.589855127006
[pyscf] R-CCSD energy of the active space using molecular orbitals=  -1852.5897242853828

What is UCCSD?¶

The UCCSD ansatz is an extension of the traditional coupled cluster method, which is widely used in classical computational chemistry. The key idea behind UCCSD is to construct the wavefunction as an exponential of an anti-Hermitian operator that includes single and double excitations. Mathematically, this can be expressed as:

\[|\Psi_{\mathrm{UCCSD}}\rangle = e^{T - T^\dagger} |\Phi\rangle\]

where \(T\) is the cluster operator defined as:

\[T = T_1 + T_2\]

Here, \(T_1\) and \(T_2\) represent the single and double excitation operators, respectively.

The UCCSD ansatz maintains the unitarity of the wavefunction, which is crucial for quantum computations. This unitarity ensures that the norm of the wavefunction remains constant, preserving the physical properties of the system.

Implementation in Quantum Computing¶

In the context of quantum computing, the UCCSD ansatz is implemented using quantum circuits. The exponential operator \(e^{T-T^\dagger}\) is decomposed into a sequence of quantum gates that can be executed on a quantum computer. This decomposition typically involves the Trotter-Suzuki approximation, which approximates the exponential of a sum of operators by a product of exponentials of the individual operators.

Below, we show how to get the UCCSD operators and how to prepare the quantum circuit using UCCSD ansatz. User can try to employ only single excitations, double excitations, or both singles and doubles excitations.

Note: the cudaq.kernels.uccsd(qubits, thetas, electron_count, qubits_num) available in cudaq.kernels uses cnot and rotation gates to implement the exponential of pauli operators in the UCCSD, see this tutorial. The version we use in this tutorial employs the exp_pauli() gate available in CUDA-Q. Using exp_pauli() gate instead of decomposing it when running on classical simulator reduces the number of gates dramatically leading to improvements in perfomance.

[7]:

from qchem.uccsd import get_uccsd_op, uccsd_circuit, uccsd_parameter_size
from qchem.uccsd import uccsd_circuit_double, uccsd_circuit_single

spin_mult, only_singles, only_doubles = 0, False, False

# Get the number of parameters
singles, doubles, total = uccsd_parameter_size(electron_count, qubits_num, spin_mult)
print(f"Number of parameters: {singles} singles, {doubles} doubles, {total} total")


# Get the UCCSD pool
if not only_singles and not only_doubles:
    word_single, word_double, coef_single, coef_double = get_uccsd_op(electron_count, qubits_num,
                                                                  spin_mult = 0, only_singles = False, only_doubles = False)
    print(f"word_single: {word_single}")
    print(f"word_double: {word_double}")
    print(f"coef_single: {coef_single}")
    print(f"coef_double: {coef_double}")

elif only_singles and not only_doubles:
    word_single, coef_single = get_uccsd_op(electron_count, qubits_num, spin_mult = 0, only_singles = True, only_doubles = False)

    print(f"word_single: {word_single}")
    print(f"coef_single: {coef_single}")

elif only_doubles and not only_singles:
    word_double, coef_double = get_uccsd_op(electron_count, qubits_num, spin_mult = 0, only_singles = False, only_doubles = True)

    print(f"word_double: {word_double}")
    print(f"coef_double: {coef_double}")
else:
    raise ValueError("Invalid option for only_singles and only_doubles")


# Get the UCCSD circuit (singles and doubles excitation are included)
@cudaq.kernel
def uccsd_kernel(qubits_num: int, electron_count: int, theta: list[float],
                 word_single: list[cudaq.pauli_word], word_double: list[cudaq.pauli_word], coef_single: list[float], coef_double: list[float]):
    """
    UCCSD kernel
    """
    # Prepare the statefrom qchem.uccsd import get_uccsd_op, uccsd_circuit, uccsd_parameter_size

    qubits = cudaq.qvector(qubits_num)

    # Initialize the qubits
    for i in range(electron_count):
        x(qubits[i])

    # Apply the UCCSD circuit
    uccsd_circuit(qubits, theta, word_single, coef_single, word_double, coef_double)


# Get the UCCSD circuit (only doubles excitations are included)
@cudaq.kernel
def uccsd_double_kernel(qubits_num: int, electron_count: int, theta: list[float],
                word_double: list[cudaq.pauli_word], coef_double: list[float]):
    """
    UCCSD kernel
    """
    # Prepare the state
    qubits = cudaq.qvector(qubits_num)

    # Initialize the qubits
    for i in range(electron_count):
        x(qubits[i])

    # Apply the UCCSD circuit
    uccsd_circuit_double(qubits, theta, word_double, coef_double)

# Get the UCCSD circuit (only singles excitations are included)
@cudaq.kernel
def uccsd_single_kernel(qubits_num: int, electron_count: int, theta: list[float],
                word_single: list[cudaq.pauli_word], coef_single: list[float]):
    """
    UCCSD kernel
    """
    # Prepare the state
    qubits = cudaq.qvector(qubits_num)

    # Initialize the qubits
    for i in range(electron_count):
        x(qubits[i])

    # Apply the UCCSD circuit
    uccsd_circuit_single(qubits, theta, word_single, coef_single)

Number of parameters: 8 singles, 18 doubles, 26 total
word_single: ['YZZZXIII', 'XZZZYIII', 'YZZZZZXI', 'XZZZZZYI', 'IIYZXIII', 'IIXZYIII', 'IIYZZZXI', 'IIXZZZYI', 'IYZZZXII', 'IXZZZYII', 'IYZZZZZX', 'IXZZZZZY', 'IIIYZXII', 'IIIXZYII', 'IIIYZZZX', 'IIIXZZZY']
word_double: ['XXIIXYII', 'XXIIYXII', 'XYIIYYII', 'YXIIYYII', 'XYIIXXII', 'YXIIXXII', 'YYIIXYII', 'YYIIYXII', 'XXIIIXYI', 'XXIIIYXI', 'XYIIIYYI', 'YXIIIYYI', 'XYIIIXXI', 'YXIIIXXI', 'YYIIIXYI', 'YYIIIYXI', 'XXIIXZZY', 'XXIIYZZX', 'XYIIYZZY', 'YXIIYZZY', 'XYIIXZZX', 'YXIIXZZX', 'YYIIXZZY', 'YYIIYZZX', 'XXIIIIXY', 'XXIIIIYX', 'XYIIIIYY', 'YXIIIIYY', 'XYIIIIXX', 'YXIIIIXX', 'YYIIIIXY', 'YYIIIIYX', 'XZZXXYII', 'XZZXYXII', 'XZZYYYII', 'YZZXYYII', 'XZZYXXII', 'YZZXXXII', 'YZZYXYII', 'YZZYYXII', 'XZZXIXYI', 'XZZXIYXI', 'XZZYIYYI', 'YZZXIYYI', 'XZZYIXXI', 'YZZXIXXI', 'YZZYIXYI', 'YZZYIYXI', 'XZZXXZZY', 'XZZXYZZX', 'XZZYYZZY', 'YZZXYZZY', 'XZZYXZZX', 'YZZXXZZX', 'YZZYXZZY', 'YZZYYZZX', 'XZZXIIXY', 'XZZXIIYX', 'XZZYIIYY', 'YZZXIIYY', 'XZZYIIXX', 'YZZXIIXX', 'YZZYIIXY', 'YZZYIIYX', 'IXXIXYII', 'IXXIYXII', 'IXYIYYII', 'IYXIYYII', 'IXYIXXII', 'IYXIXXII', 'IYYIXYII', 'IYYIYXII', 'IXXIIXYI', 'IXXIIYXI', 'IXYIIYYI', 'IYXIIYYI', 'IXYIIXXI', 'IYXIIXXI', 'IYYIIXYI', 'IYYIIYXI', 'IXXIXZZY', 'IXXIYZZX', 'IXYIYZZY', 'IYXIYZZY', 'IXYIXZZX', 'IYXIXZZX', 'IYYIXZZY', 'IYYIYZZX', 'IXXIIIXY', 'IXXIIIYX', 'IXYIIIYY', 'IYXIIIYY', 'IXYIIIXX', 'IYXIIIXX', 'IYYIIIXY', 'IYYIIIYX', 'IIXXXYII', 'IIXXYXII', 'IIXYYYII', 'IIYXYYII', 'IIXYXXII', 'IIYXXXII', 'IIYYXYII', 'IIYYYXII', 'IIXXIXYI', 'IIXXIYXI', 'IIXYIYYI', 'IIYXIYYI', 'IIXYIXXI', 'IIYXIXXI', 'IIYYIXYI', 'IIYYIYXI', 'IIXXXZZY', 'IIXXYZZX', 'IIXYYZZY', 'IIYXYZZY', 'IIXYXZZX', 'IIYXXZZX', 'IIYYXZZY', 'IIYYYZZX', 'IIXXIIXY', 'IIXXIIYX', 'IIXYIIYY', 'IIYXIIYY', 'IIXYIIXX', 'IIYXIIXX', 'IIYYIIXY', 'IIYYIIYX', 'XZXIXZYI', 'XZXIYZXI', 'XZYIYZYI', 'YZXIYZYI', 'XZYIXZXI', 'YZXIXZXI', 'YZYIXZYI', 'YZYIYZXI', 'IXZXIXZY', 'IXZXIYZX', 'IXZYIYZY', 'IYZXIYZY', 'IXZYIXZX', 'IYZXIXZX', 'IYZYIXZY', 'IYZYIYZX']
coef_single: [0.5, -0.5, 0.5, -0.5, 0.5, -0.5, 0.5, -0.5, 0.5, -0.5, 0.5, -0.5, 0.5, -0.5, 0.5, -0.5]
coef_double: [0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125, 0.125, 0.125, 0.125, 0.125, -0.125, -0.125, -0.125, -0.125]

Run VQE¶

To learn about VQE in CUDA-Q, check this tutorial.

[ ]:

import numpy as np
from scipy.optimize import minimize

# Initial guess for the parameters
if not only_singles and not only_doubles:
    theta = [0.0] * total
elif only_singles and not only_doubles:
    theta = [0.0] * singles
elif only_doubles and not only_singles:
    theta = [0.0] * doubles
else:
    raise ValueError("Invalid option for only_singles and only_doubles")

#print(cudaq.draw(kernel, qubits_num, electron_count, theta, word_single, word_double, coef_single, coef_double))

# The ansatz is a variational circuit, so we need to optimize the parameters

# Parameter shift to compute the gradient
def parameter_shift(theta):

    parameter_count = len(theta)
    epsilon = np.pi / 4
    # The gradient is calculated using parameter shift.
    grad = np.zeros(parameter_count)
    theta2 = theta.copy()

    for i in range(parameter_count):
        theta2[i] = theta[i] + epsilon
        exp_val_plus = cost(theta2)
        theta2[i] = theta[i] - epsilon
        exp_val_minus = cost(theta2)
        grad[i] = (exp_val_plus - exp_val_minus) / (2 * epsilon)
        theta2[i] = theta[i]
    return grad

def cost(theta):

    theta=theta.tolist()

    if not only_singles and not only_doubles:
        energy = cudaq.observe(uccsd_kernel, hamiltonian, qubits_num, electron_count,
                            theta, word_single, word_double, coef_single, coef_double).expectation()

    elif only_singles and not only_doubles:
        energy = cudaq.observe(uccsd_single_kernel, hamiltonian, qubits_num, electron_count,
                            theta, word_single, coef_single).expectation()
    elif only_doubles and not only_singles:
        energy = cudaq.observe(uccsd_double_kernel, hamiltonian, qubits_num, electron_count,
                            theta, word_double, coef_double).expectation()
    else:
        raise ValueError("Invalid option for only_singles and only_doubles")

    return energy

#result_vqe=minimize(cost, theta, method='L-BFGS-B', jac='2-point', tol=1e-7)
result_vqe=minimize(cost, theta, method='L-BFGS-B', jac=parameter_shift, tol=1e-7)

total_energy = result_vqe.fun + const
print(f"Total energy: {total_energy:.10f} Hartree")
# Print the optimized parameters: first n are singles, then doubles.
print(f"optimized parameters: {result_vqe.x}")

Total energy: -1852.5895309905 Hartree
optimized parameters: [ 0.00000000e+00  1.91918143e-03 -1.18739821e-01 -4.56567309e-18
  0.00000000e+00  1.93627077e-03 -1.18755472e-01  1.54888818e-21
  9.53533311e-03  0.00000000e+00  0.00000000e+00  1.40615591e-03
  0.00000000e+00 -1.02368050e-02  2.78804812e-03  0.00000000e+00
  0.00000000e+00  2.78862284e-03 -1.01582650e-02  0.00000000e+00
  3.21798016e-01 -9.94788673e-15  9.94560544e-15  1.54838228e-02
 -4.22692720e-03 -4.25296495e-03]

Challenges and consideration¶

While UCCSD is a powerful method, it is not without its challenges. The primary difficulty lies in the efficient implementation of the quantum circuits, as the number of gates required can grow significantly with the size of the system. Additionally, the accuracy of the Trotter-Suzuki approximation can affect the overall precision of the results.