May 06, 2024 | 4 min read

\( \newcommand{\ket}[1]{|#1\rangle} \)

Ansatz k-UpCCGSD

Show code cell source Hide code cell source

import qiskit
# from qiskit_ibm_provider import IBMProvider
# from qiskit import IBMQ
# import qiskit.tools.jupyter
from qiskit.tools.parallel import parallel_map
import json
from styles.style import qspain

from qiskit.algorithms.minimum_eigensolvers import VQE
from qiskit.algorithms.optimizers import COBYLA, SPSA, POWELL
from qiskit.primitives import Estimator

from qiskit_nature.second_q.circuit.library import UCC, PUCCD, HartreeFock
from qiskit_nature.second_q.mappers import QubitMapper, JordanWignerMapper
from qiskit_nature.second_q.drivers import PySCFDriver
from qiskit_nature.units import DistanceUnit
from qiskit_nature.second_q.algorithms import GroundStateEigensolver

import pennylane as qml
from pennylane import numpy as np


%matplotlib inline
#IBMQ.load_account()

#%qiskit_backend_overview

Teoría de Coupled Cluster

Para obtener el estado fundamental de una molécula, dado que son sistemas relativamente complejos en los que los diferentes electrones pueden ocupar muchas posiciones en los diferentes orbitales, es prácticamente imposible obtener un resultado exacto. En su lugar, se utilizan métodos aproximados.

El más famoso de estos métodos, así como el punto de partida para empezar a aplicar el resto, es el método de Hartree-Fock. Con él se pueden aproximar las funciones de onda de los electrones del sistema atendiendo a la configuración de los orbitales atómicos o moleculares.

Sin embargo, este método está bastante limitado por varios motivos: no considera la correlación electrónica, hace uso de la aproximación de Born-Oppenheimer para los núcleos, no tiene en cuenta posibles efectos relativistas… Es por todo esto que, tras obtener una primera aproximación del estado fundamental, se utilizan otros métodos para terminar de “pulir” este resultado como Configuration Interaction (CI) o Møller–Plesset perturbation theory (MP) [2].

Es dentro de estos métodos post-Hartree-Fock donde encontramos los métodos de Coupled Cluster [7], [13], [1]. La base del CC es, en resumidas cuentas, considerar de forma ponderada cómo contribuyen los estados con electrones a la energía del estado fundamental molecular. Para ello se añaden términos a la función de onda considerando las excitaciones de un electron de los orbitales libres a los orbitales ocupados (Single excitation), después de dos electrones con todas las combinaciones posibles (Double) y así sucesivamente.

Los métodos CC se usan para crear sistemas correlacionados utilizando un operador de cluster \(\hat{T}\), de forma que el estado fabricado se genera a partir de otro de referencia haciendo

\[ \ket{\psi} = e^{\hat{T}} \ket{\phi_0} \]

siendo \(\ket{\phi_0}\) el estado de referencia en cuestión, típicamente el de Hartree-Fock.

Dentro de estos métodos encontramos los CCSD (Single and Double excitations Coupled Cluster), es decir, métodos que se limitan a considerar excitaciones simples o dobles de los orbitales. Cuando trabajamos con CCSD, el operador de cluster toma la forma

\[ \hat{T} = \hat{T}_1 + \hat{T}_2 \]

siendo

\[ \hat{T}_1 = \sum_{i,a} t_i^a \hat{a}^\dagger_a \hat{a}_i \]

y

\[ \hat{T}_2 = \frac{1}{4} \sum_{i,j ; a,b} t^{ab}_{ij} \hat{a}^\dagger_{a} \hat{a}^\dagger_{b} \hat{a}_j \hat{a}_i \]

Los índices \(i,j\) denotan orbitales ocupados, mientras que los índices \(a,b\) denotan orbitales libres. Los parámetros a optimizar por tanto por el algoritmo VQE son precisamente estos \(t^a_i, t^{ab}_{ij}\).

Los operadores CC presentan un problema a la hora de implementarlos como parte de un algoritmo VQE, ya que la forma de proyectar el Hamiltoniano del sistema sobre los estados fabricados para calcular la energía no es compatible con el método variacional, que es la columna vertebral del VQE. Para solventar este problema tenemos un caso particular de los métodos CC: los métodos UCC (Unitary Coupled Cluster).

La solución al problema variacional es muy simple: la forma de fabricar el ansatz pasa a ser

\[ \ket{\psi} = e^{\hat{T} - \hat{T}^\dagger} \ket{\phi_0} \]

Al igual que con los CC, podemos definir métodos UCCSD, con excitaciones simples y dobles, tomando de nuevo \( \hat{T} = \hat{T}_1 + \hat{T}_2 \). Un ordenador cuántico puede implementar eficientemente estos métodos UCCSD utilizando una trotterización del operador de cluster exponenciado.

Después de esto, podemos buscar otro caso de UCC, el llamado UCCG (Unitary Coupled Cluster with Generalized wave functions), el cual no distingue entre orbitales ocupados y libres. Es por ello que, en el caso UCCGSD, el operador de cluster vuelve a estar definido como \( \hat{T} = \hat{T}_1 + \hat{T}_2 \), pero esta vez

\[ \hat{T}_1 = \frac{1}{2} \sum_{p,q} t^q_p \hat{a}^\dagger_q \hat{a}_p \]

y

\[ \hat{T}_2 = \frac{1}{4} \sum_{p,q,r,s} t^{rs}_{pq} \hat{a}^\dagger_r \hat{a}^\dagger_s \hat{a}_q \hat{a}_p \]

denotando los cuatro índices \(p,q,r,s\) cualquier tipo de orbital sin ocupar su ocupación.

Estos últimos ansätze se pueden implementar en qiskit utilizando qiskit_nature.

De aquí podemos sacar una clase un poco más restringida, la clase UpCCGSD (Unitary Pair Coupled Cluster with Generalized Singles and Doubles product wave functions), que limita las amplitudes de las excitaciones dobles considerando únicamente las que mueven un par de electrones de un orbital espacial a otro. Esto provoca un cambio en el operador \(\hat{T}_2\), que se convierte en

\[ \hat{T}_2 = \sum_{p,q} t^{q_\alpha q_\beta}_{p_\alpha p_\beta} \hat{a}^\dagger_{q_\alpha} \hat{a}^\dagger_{q_\beta} \hat{a}_{p_\beta} \hat{a}_{p_\alpha} \]

Y de aquí por fin llegamos a la clase que queremos estudiar, k-UpCCGSD [10], [8], [11], con \(k \in \mathbb{Z}\). Con este entero \(k\) definimos nuestra función de prueba \(\ket{\psi}\) como

\[ \ket{\psi} = \left( \prod_{\mu = 1}^k e^{\hat{T}^{(\mu)} - \hat{T}^{(\mu) \dagger}} \right) \ket{\phi_0} \]

haciendo que cada \( \hat{T}^{(\mu)} \) tenga su propio conjunto de parámetros \(t^q_p, t^{q_\alpha q_\beta}_{p_\alpha p_\beta}\).

Utilizar k-UpCCGSD

Las mejores alternativas para implementar el ansatz k-UpCCGSD son qiskit, librería con la que ya estamos, familiarizados, y pennylane otra librería enfocada a la Computación Cuántica y, más concretamente, a Quantum Machine Learning.

Como ejemplo prototípico para problemas moleculares, vamos a resolver la molécula de hidrógeno \(H_2\).

Pennylane

# Construimos la molécula de Hidrógeno H2
symbols = ["H", "H"]
coordinates = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 1.400])
H, qubits = qml.qchem.molecular_hamiltonian(symbols, coordinates)

# Obtenemos el estado de referencia de Hartree-Fock
electrons = 2
ref_state = qml.qchem.hf_state(electrons, qubits)

Qiskit

# Construimos la molécula de Hidrógeno H2
driver = PySCFDriver(
    atom="H 0 0 0; H 0 0 0.741",
    basis="sto3g",
    charge=0,
    spin=0,
    unit=DistanceUnit.ANGSTROM,
    # El método de aproximación por defecto es Hartree-Fock: method=MethodType.RHF
)

es_problem = driver.run()

# Utilizamos el mapping de Jordan-Wigner
mapper = JordanWignerMapper()

# Obtenemos el estado de referencia de Hartree-Fock
ref_state = HartreeFock(es_problem.num_spatial_orbitals, es_problem.num_particles, mapper)

Un detalle importante a tener en cuenta a la hora de trabajar con estas librerías es que, mientras que en qiskit podemos decidir con qué unidades trabajamos, en pennylane trabajamos con unidades atómicas.En estas unidades, la distancia de enlace entre los átomos de hidrógeno de la molécula es \( 0.741 \cdot 10^{-10} \text{m} =1.40 a_0\)

Pennylane

# En pennylane tenemos que definir el circuito sobre el que trabajará el ansatz
dev = qml.device('default.qubit', wires=qubits)

# Definimos el ansatz k-UpCCGSD
# weights son los parámetros t^q_p, t^{q_\alpha q_\beta}_{p_\alpha p_\beta}
# wires son los qubits (o en general los sistemas) sobre los que se define el ansatz.
# delta_sz especifica las reglas de selección para la proyección del espín en dirección z 
#     en los orbitales de las excitaciones simples.

k = 3

@qml.qnode(dev)
def ansatz(weights):
    qml.kUpCCGSD(weights, wires=[0, 1, 2, 3],
                    k=k, delta_sz=0, init_state=ref_state)
    return qml.expval(H)

shape = qml.kUpCCGSD.shape(k=k,
                    n_wires=qubits, delta_sz=0)

Qiskit

# Preparamos el ansatz k-UpCCGSD

k = 3

# Para preparar este ansatz partimos de una función que ejecuta un ansatz UpCCD
# Podemos añadirle las k repeticiones con la opción reps
# Podemos añadirle las excitaciones S con la opción include_singles
# Podemos añadirle las excitaciones generalizadas con la opción generalized

ansatz = PUCCD(
    num_spatial_orbitals=es_problem.num_spatial_orbitals, 
    num_particles=es_problem.num_particles, 
    qubit_mapper=mapper, 
    include_singles=(True,True), 
    generalized=True, 
    reps=k
)

Lo siguiente que tenemos que hacer es escoger el optimizador clásico.

Pennylane

# Escogemos un optimizador
opt = qml.GradientDescentOptimizer(stepsize=0.4)
# opt = qml.GradientDescentOptimizer(stepsize=0.04)
# opt = qml.QNGOptimizer(stepsize=0.01)
# opt = qml.RMSPropOptimizer(stepsize=0.01)
# opt = qml.AdamOptimizer(stepsize=0.01)
# opt = qml.AdagradOptimizer(stepsize=0.01)
# opt = qml.SPSAOptimizer(maxiter=250) 

Qiskit

# Escogemos un optimizador
opt = COBYLA(maxiter=1000)
# opt = POWELL()

Y ya podemos resolver el problema VQE

Pennylane

# Inicializamos los primeros parámetros
weights = np.random.random(size=shape)
energy = [ansatz(weights)]

# E iniciamos el proceso iterativo
angle = [weights]
max_iterations = 250
conv_tol = 1e-08 # Escogemos una tolerancia en la que queremos que se detenga el algoritmo
for n in range(max_iterations):
    weights, prev_energy = opt.step_and_cost(ansatz, weights)
    energy.append(ansatz(weights))
    angle.append(weights)
    conv = np.abs(energy[-1] - prev_energy)
    # if n % 4 == 0:
        # print(f"Iteración: {n},  Energía = {energy[-1]:.8f} Ha")
    if conv <= conv_tol:
        break

print("\n" f"Energía del estado fundamental (obtenida) = {energy[-1]:.8f} Ha")
print("\n" f"Valores óptimos del circuito = {angle[-1]}")

Qiskit

# Ejecutamos el algoritmo
vqe_solver = VQE(Estimator(), ansatz, opt)

calc = GroundStateEigensolver(mapper, vqe_solver)

res = calc.solve(es_problem)
best = res.total_energies
print('E_G =',best)

Energía del estado fundamental (obtenida) = -1.13727592 Ha

Valores óptimos del circuito = [[0.64602866 0.34079656 0.42082063 0.16922154 0.62086335 0.3463154 ]
 [0.2275569  0.04592067 0.12976618 0.68158533 0.29590652 0.35328037]
 [0.39530577 0.64276052 0.72391752 0.18231534 0.23546012 0.22749385]]

Mostramos aquí algunas tablas del valor al que converge la energía utilizando diferentes optimizadores y valores de \(k\).

k	E (Ha) GradDesc	E (Ha) SPSA	E (Ha) ADAM	E (Ha) RMSProp	E (Ha) Adargrad	E (Ha) QNG
3	-1.13727592	-1.13575585	-1.13727562	-1.13628251	-1.1372758	-1.13727591
4	-1.13727592	-1.1355345	-1.13727573	-1.13727594	-1.13727586	-1.13727586
5	-1.13727585	-1.13606927	-1.13726938	-1.13465174	-1.1372759	-1.13727585
6	-1.13726475	-1.13358301	-1.13727593	-1.13727594	-1.1372757	-1.13727587
7	-1.13727592	-1.13560152	-1.13727591	-1.13339437	-1.13727577	-1.1372759

---

Authors:

Irais Bautista (CESGA), Sergio Martínez (BIFi-UNIZAR), Jaime Scharfhausen (UAM) y Alejandro Jaramillo (CSUC)

---

License: Licencia Creative Commons Atribución-CompartirIgual 4.0 Internacional.

This work has been financially supported by the Ministry for Digital Transformation and of Civil Service of the Spanish Government through the QUANTUM ENIA project call - Quantum Spain project, and by the European Union through the Recovery, Transformation and Resilience Plan - NextGenerationEU within the framework of the Digital Spain 2026 Agenda.