QHack 2023: 4 week prep challenge, day 6, H-H first excited state

Almost a week passed, it’s time for more tricky tasks! The last and most advanced exercise in the Quantum Chemistry track of the 2022 QHack.

Disclaimer: At the time being, the proposed solution does not produce the correct results.

The task is to obtain the energy of the first excited state of an H-H molecule at a given nuclei distance. This is achieved by firstly calculating the ground state energy, and later modifying the systems Hamiltonian to penalize the actual ground state. Repeating the ground-state calculation with the modified Hamiltonian, you can obtain the first excited state.

The calculation of the ground-state via a VQE routine, is done following the concise Pennylane VQE tutorial, with only minuscule modifications.

import sys
import pennylane as qml
from pennylane import numpy as np
from pennylane import hf


def ground_state_VQE(H):
    """Perform VQE to find the ground state of the H2 Hamiltonian.

    Args:
        - H (qml.Hamiltonian): The Hydrogen (H2) Hamiltonian

    Returns:
        - (float): The ground state energy
        - (np.ndarray): The ground state calculated through your optimization routine
    """

    qubits = 4
    
    dev = qml.device("default.qubit", wires=qubits)
    
    electrons = 2
    hf = qml.qchem.hf_state(electrons, qubits)
    
    def circuit(param, wires):
        qml.BasisState(hf, wires=wires)
        qml.DoubleExcitation(param, wires=range(qubits))
    
    @qml.qnode(dev)
    def cost_fn(param):
        circuit(param, wires=range(qubits))
        return qml.expval(H)
    
    @qml.qnode(dev)
    def ground_state(param, wires):
        circuit(param, wires)
        return qml.state()
    
    opt = qml.GradientDescentOptimizer(stepsize=0.1)
    theta = np.array(0.0, requires_grad=True)
    
    # store the values of the cost function
    energy = [cost_fn(theta)]

    # store the values of the circuit parameter
    angle = [theta]

    max_iterations = 100
    conv_tol = 1e-06

    for n in range(max_iterations):
        theta, prev_energy = opt.step_and_cost(cost_fn, theta)

        energy.append(cost_fn(theta))
        angle.append(theta)

        conv = np.abs(energy[-1] - prev_energy)

        if n % 2 == 0:
            print(f"Step = {n},  Energy = {energy[-1]:.8f} Ha, theta = {theta:.3f}")

        if conv <= conv_tol:
            break

    E0 = energy[-1]
    gstate = ground_state(theta, wires=range(qubits))
    
    return E0, gstate

symbols = ["H", "H"]
coord = 0.6614
geometry = np.array([[0.0, 0.0, -coord], [0.0, 0.0, coord]], requires_grad=False)
mol = hf.Molecule(symbols, geometry)

H = hf.generate_hamiltonian(mol)()
E0, ground_state = ground_state_VQE(H)

The next step is to calculate a modified Hamiltonian.

def create_H1(ground_state, beta, H):
    """Create the H1 matrix, then use `qml.Hermitian(matrix)` to return an observable-form of H1.

    Args:
        - ground_state (np.ndarray): from the ground state VQE calculation
        - beta (float): the prefactor for the ground state projector term
        - H (qml.Hamiltonian): the result of hf.generate_hamiltonian(mol)()

    Returns:
        - (qml.Observable): The result of qml.Hermitian(H1_matrix)
    """

    H1 = qml.matrix(H) + beta * np.outer(ground_state, ground_state)
    qubits = int(np.log2(len(ground_state)))
    return qml.Hermitian(H1, wires=range(qubits))


def excited_state_VQE(H1):
    """Perform VQE using the "excited state" Hamiltonian.

    Args:
        - H1 (qml.Observable): result of create_H1

    Returns:
        - (float): The excited state energy
    """

    E1, estate = ground_state_VQE(H1)
    return E1

beta = 15.0
H1 = create_H1(ground_state, beta, H)
E1 = excited_state_VQE(H1)

print(E1)

While the ground state energy is correctly calculated, some error remains in the calculation of the excited state energy (larger than zero). … I will leave it to rest for the next few days, and maybe revisit it later.

Related