Basis rotation circuit

examples/br_example.py

@@ -0,0 +1,82 @@
+"""
+Example of the basis rotation circuit with H3+ molecule. Starts with the guess wave function from the core Hamiltonian,
+    and runs the VQE to obtain the HF energy.
+"""
+
+
+import numpy as np
+from qibo.optimizers import optimize
+from scipy.optimize import minimize
+
+from qibochem.ansatz.basis_rotation import br_circuit
+from qibochem.ansatz.hf_reference import hf_circuit
+from qibochem.driver.molecule import Molecule
+from qibochem.measurement.expectation import expectation
+
+# Define molecule and populate
+mol = Molecule(xyz_file="h3p.xyz")
+try:
+    mol.run_pyscf()
+except ModuleNotFoundError:
+    mol.run_psi4()
+
+
+# Diagonalize H_core to use as the guess wave function
+
+# First, symmetric orthogonalization of overlap matrix S using np:
+u, s, vh = np.linalg.svd(mol.overlap)
+A = u @ np.diag(s ** (-0.5)) @ vh
+
+# Transform guess Fock matrix formed with H_core
+F_p = A.dot(mol.hcore).dot(A)
+# Diagonalize F_p for eigenvalues and eigenvectors
+_e, C_p = np.linalg.eigh(F_p)
+# Transform C_p back into AO basis
+C = A.dot(C_p)
+# Form OEI/TEI with the (guess) MO coefficients
+oei = np.einsum("pQ, pP -> PQ", np.einsum("pq, qQ -> pQ", mol.hcore, C, optimize=True), C, optimize=True)
+# TEI code from https://pycrawfordprogproj.readthedocs.io/en/latest/Project_04/Project_04.html
+tei = np.einsum("up, vq, uvkl, kr, ls -> prsq", C, C, mol.aoeri, C, C, optimize=True)
+
+# Molecular Hamiltonian with the guess OEI/TEI
+hamiltonian = mol.hamiltonian(oei=oei, tei=tei)
+
+# Check that the hamiltonian with a HF reference ansatz doesn't yield the correct HF energy
+circuit = hf_circuit(mol.nso, sum(mol.nelec))
+print(
+    f"\nElectronic energy: {expectation(circuit, hamiltonian):.8f} (From the H_core guess, should be > actual HF energy)"
+)
+
+
+def electronic_energy(parameters):
+    """
+    Loss function (Electronic energy) for the basis rotation ansatz
+    """
+    circuit = hf_circuit(mol.nso, sum(mol.nelec))
+    circuit += br_circuit(mol.nso, parameters, sum(mol.nelec))
+
+    return expectation(circuit, hamiltonian)
+
+
+# Build  a basis rotation_circuit
+params = np.random.rand(sum(mol.nelec) * (mol.nso - sum(mol.nelec)))  # n_occ * n_virt
+
+best, params, extra = optimize(electronic_energy, params)
+
+print("\nResults using Qibo optimize:")
+print(f" HF energy: {mol.e_hf:.8f}")
+print(f"VQE energy: {best:.8f} (Basis rotation ansatz)")
+# print()
+# print("Optimized parameters:", params)
+
+
+params = np.random.rand(sum(mol.nelec) * (mol.nso - sum(mol.nelec)))  # n_occ * n_virt
+
+result = minimize(electronic_energy, params)
+best, params = result.fun, result.x
+
+print("\nResults using scipy.optimize:")
+print(f" HF energy: {mol.e_hf:.8f}")
+print(f"VQE energy: {best:.8f} (Basis rotation ansatz)")
+# print()
+# print("Optimized parameters:", params)

examples/h3p.xyz

@@ -0,0 +1,5 @@
+3
+1 1
+H 0.0 0.0 0.0
+H 0.0 0.0 1.0
+H 0.0 0.0 2.0

examples/ucc_example1.py

@@ -9,6 +9,7 @@
 from qibochem.ansatz.hf_reference import hf_circuit
 from qibochem.ansatz.ucc import ucc_circuit
 from qibochem.driver.molecule import Molecule
+from qibochem.measurement.expectation import expectation
 
 # Define molecule and populate
 mol = Molecule(xyz_file="lih.xyz")
@@ -89,7 +90,7 @@ def electronic_energy(parameters):
     all_parameters = [_x for _param in all_parameters for _x in _param]
     circuit.set_parameters(all_parameters)
 
-    return mol.expectation(circuit, hamiltonian)
+    return expectation(circuit, hamiltonian)
 
 
 # Reference energy

examples/ucc_example2.py

@@ -12,6 +12,7 @@
 from qibochem.ansatz.hf_reference import hf_circuit
 from qibochem.ansatz.ucc import ucc_ansatz
 from qibochem.driver.molecule import Molecule
+from qibochem.measurement.expectation import expectation
 
 # Define molecule and populate
 mol = Molecule(xyz_file="lih.xyz")
@@ -47,7 +48,7 @@ def electronic_energy(parameters):
     """
     circuit.set_parameters(parameters)
 
-    return mol.expectation(circuit, hamiltonian)
+    return expectation(circuit, hamiltonian)
 
 
 # Reference energy

src/qibochem/ansatz/basis_rotation.py

@@ -0,0 +1,143 @@
+"""
+Circuit representing an unitary rotation of the molecular (spin-)orbital basis set
+"""
+
+import numpy as np
+from openfermion.linalg.givens_rotations import givens_decomposition
+from qibo import gates, models
+from scipy.linalg import expm
+
+# Helper functions
+
+
+def unitary_rot_matrix(parameters, occ_orbitals, virt_orbitals):
+    r"""
+    Returns the unitary rotation matrix U = exp(\kappa) mixing the occupied and virtual orbitals,
+         where \kappa is an anti-Hermitian matrix
+
+    Args:
+        parameters: List/array of rotation parameters. dim = len(occ_orbitals)*len(virt_orbitals)
+        occ_orbitals: Iterable of occupied orbitals
+        virt_orbitals: Iterable of virtual orbitals
+    """
+    n_orbitals = len(occ_orbitals) + len(virt_orbitals)
+
+    kappa = np.zeros((n_orbitals, n_orbitals))
+    # Get all possible (occ, virt) pairs
+    orbital_pairs = [(_o, _v) for _o in occ_orbitals for _v in virt_orbitals]
+    # occ/occ and virt/virt orbital mixing doesn't change the energy, so can ignore
+    for _i, (_occ, _virt) in enumerate(orbital_pairs):
+        kappa[_occ, _virt] = parameters[_i]
+        kappa[_virt, _occ] = -parameters[_i]
+    return expm(kappa)
+
+
+def swap_matrices(permutations, n_qubits):
+    r"""
+    Build matrices that permute (swap) the columns of a given matrix
+
+    Args:
+        permutations [List(tuple, ), ]: List of lists of 2-tuples permuting two consecutive numbers
+            e.g. [[(0, 1), (2, 3)], [(1, 2)], ...]
+        n_qubits: Dimension of matrix (\exp(\kappa) matrix) to be operated on i.e. # of qubits
+
+    Returns:
+        List of matrices that carry out \exp(swap) for each pair in permutations
+    """
+    exp_swaps = []
+    _temp = np.array([[-1.0, 1.0], [1.0, -1.0]])
+    for pairs in permutations:
+        gen = np.zeros((n_qubits, n_qubits))
+        for pair in pairs:
+            # Add zeros to pad out the (2, 2) matrix to (n_qubits, n_qubits) with 0.
+            gen += np.pad(_temp, ((pair[0], n_qubits - pair[1] - 1),), "constant")
+        # Take matrix exponential, remove all entries close to zero, and convert to real matrix
+        exp_mat = expm(-0.5j * np.pi * gen)
+        exp_mat[np.abs(exp_mat) < 1e-10] = 0.0
+        exp_mat = np.real_if_close(exp_mat)
+        # Add exp_mat to exp_swaps once cleaned up
+        exp_swaps.append(exp_mat)
+    return exp_swaps
+
+
+def givens_rotation_parameters(n_qubits, exp_k, n_occ):
+    r"""
+    Get the parameters of the Givens rotation matrix obtained after decomposing \exp(\kappa)
+
+    Args:
+        n_qubits: Number of qubits (i.e. spin-orbitals)
+        exp_k: Unitary rotation matrix
+        n_occ: Number of occupied orbitals
+    """
+    # List of 2-tuples to link all qubits via swapping
+    swap_pairs = [
+        [(_i, _i + 1) for _i in range(_d % 2, n_qubits - 1, 2)]
+        # Only 2 possible ways of swapping: [(0, 1), ... ] or [(1, 2), ...]
+        for _d in range(2)
+    ]
+    # Get their corresponding matrices
+    swap_unitaries = swap_matrices(swap_pairs, n_qubits)
+
+    # Create a copy of exp_k for rotating
+    exp_k_shift = np.copy(exp_k)
+    # Apply the column-swapping transformations
+    for u_rot in swap_unitaries:
+        exp_k_shift = u_rot @ exp_k_shift
+    # Only want the n_occ by n_qubits (n_orb) sub-matrix
+    qmatrix = exp_k_shift.T[:n_occ, :]
+
+    # Apply Givens decomposition of the submatrix to get a list of individual Givens rotations
+    # (orb1, orb2, theta, phi), where orb1 and orb2 are rotated by theta(real) and phi(imag)
+    qdecomp, _, _ = givens_decomposition(qmatrix)
+    # Lastly, reverse the list of Givens rotations (because U = G_k ... G_2 G_1)
+    return list(reversed(qdecomp))
+
+
+def givens_rotation_gate(n_qubits, orb1, orb2, theta):
+    """
+    Circuit corresponding to a Givens rotation between two qubits (spin-orbitals)
+
+    Args:
+        n_qubits: Number of qubits in the quantum circuit
+        orb1, orb2: Orbitals used for the Givens rotation
+        theta: Rotation parameter
+
+    Returns:
+        circuit: Qibo Circuit object representing a Givens rotation between orb1 and orb2
+    """
+    # Define 2x2 sqrt(iSwap) matrix
+    iswap_mat = np.array([[1.0, 1.0j], [1.0j, 1.0]]) / np.sqrt(2.0)
+    # Build and add the gates to circuit
+    circuit = models.Circuit(n_qubits)
+    circuit.add(gates.GeneralizedfSim(orb1, orb2, iswap_mat, 0.0, trainable=False))
+    circuit.add(gates.RZ(orb1, -theta))
+    circuit.add(gates.RZ(orb2, theta + np.pi))
+    circuit.add(gates.GeneralizedfSim(orb1, orb2, iswap_mat, 0.0, trainable=False))
+    circuit.add(gates.RZ(orb1, np.pi, trainable=False))
+    return circuit
+
+
+def br_circuit(n_qubits, parameters, n_occ):
+    r"""
+    Google's basis rotation circuit between the occupied/virtual orbitals. Forms the \exp(\kappa) matrix, decomposes
+        it into Givens rotations, and sets the circuit parameters based on the Givens rotation decomposition
+
+    Args:
+        n_qubits: Number of qubits in the quantum circuit
+        parameters: Rotation parameters for \exp(\kappa)
+        n_occ: Number of occupied orbitals
+    """
+    # Unitary rotation matrix \exp(\kappa)
+    exp_k = unitary_rot_matrix(parameters, range(n_occ), range(n_occ, n_qubits))
+    # List of Givens rotation parameters
+    g_rotation_parameters = givens_rotation_parameters(n_qubits, exp_k, n_occ)
+
+    # Start building the circuit
+    circuit = models.Circuit(n_qubits)
+    # Add the Givens rotation circuits
+    for g_rot_parameter in g_rotation_parameters:
+        for orb1, orb2, theta, _phi in g_rot_parameter:
+            if not np.allclose(_phi, 0.0):
+                raise ValueError("Unitary rotation is not real")
+            circuit += givens_rotation_gate(n_qubits, orb1, orb2, theta)
+    return circuit

tests/test_basis_rotation.py

@@ -0,0 +1,42 @@
+"""
+Test for basis rotation ansatz
+"""
+
+import numpy as np
+import pytest
+from qibo import Circuit, gates
+from qibo.optimizers import optimize
+
+from qibochem.ansatz.basis_rotation import br_circuit
+from qibochem.driver.molecule import Molecule
+from qibochem.measurement.expectation import expectation
+
+
+def test_br_ansatz():
+    """Test of basis rotation ansatz against hardcoded HF energies"""
+    h2_ref_energy = -1.117349035
+
+    mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.7))])
+    mol.run_pyscf(max_scf_cycles=1)
+    # Use an un-converged wave function to build the Hamiltonian
+    mol_sym_ham = mol.hamiltonian()
+
+    # Define quantum circuit
+    circuit = Circuit(mol.nso)
+    circuit.add(gates.X(_i) for _i in range(sum(mol.nelec)))
+
+    # Add basis rotation ansatz
+    # Initialize all at zero
+    parameters = np.zeros(sum(mol.nelec) * (mol.nso - sum(mol.nelec)))  # n_occ * n_virt
+    circuit += br_circuit(mol.nso, parameters, sum(mol.nelec))
+
+    def electronic_energy(parameters):
+        """
+        Loss function (Electronic energy) for the basis rotation ansatz
+        """
+        circuit.set_parameters(parameters)
+        return expectation(circuit, mol_sym_ham)
+
+    hf_energy, parameters, _extra = optimize(electronic_energy, parameters)
+
+    assert hf_energy == pytest.approx(h2_ref_energy)