merge with main

.pre-commit-config.yaml

@@ -2,14 +2,14 @@ ci:
   autofix_prs: true
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.4.0
+    rev: v4.5.0
     hooks:
       - id: trailing-whitespace
       - id: end-of-file-fixer
       - id: check-yaml
       - id: debug-statements
   - repo: https://github.com/psf/black
-    rev: 23.7.0
+    rev: 23.10.1
     hooks:
       - id: black
         args:
@@ -20,6 +20,6 @@ repos:
       - id: isort
         args: ["--profile", "black"]
   - repo: https://github.com/asottile/pyupgrade
-    rev: v3.10.1
+    rev: v3.15.0
     hooks:
       - id: pyupgrade

doc/source/api-reference/ansatz.rst

@@ -1,4 +1,21 @@
 Ansatz
 ======
 
-This section covers the API reference for the ansatz module.
+This section covers the API reference for various chemistry ansatzes.
+
+Hardware-efficient
+------------------
+
+.. autofunction:: qibochem.ansatz.hardware_efficient.hea
+
+
+Hartree-Fock
+------------
+
+.. autofunction:: qibochem.ansatz.hf_reference.hf_circuit
+
+
+Unitary Coupled Cluster
+-----------------------
+
+.. autofunction:: qibochem.ansatz.ucc.ucc_circuit

doc/source/conf.py

@@ -15,6 +15,10 @@
 from recommonmark.transform import AutoStructify
 
 sys.path.insert(0, os.path.abspath(".."))
+<<<<<<< HEAD
+=======
+
+>>>>>>> main
 import qibochem
 
 # -- Project information -----------------------------------------------------
@@ -26,16 +30,19 @@
 
 # The full version, including alpha/beta/rc tags
 release = qibochem.__version__
-
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration
 # https://stackoverflow.com/questions/56336234/build-fail-sphinx-error-contents-rst-not-found
 master_doc = "index"
 
+master_doc = "index"
+
 extensions = [
     "sphinx.ext.autodoc",
+    "sphinx.ext.napoleon",
     "sphinx.ext.doctest",
     "sphinx.ext.coverage",
+    "recommonmark",
     "sphinxcontrib.katex",
     "sphinx.ext.napoleon",
     "sphinx.ext.intersphinx",
@@ -131,5 +138,4 @@ def setup(app):
     app.add_transform(AutoStructify)
     app.add_css_file("css/style.css")
 
-
 html_show_sourcelink = False

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

pyproject.toml

@@ -20,6 +20,7 @@ qibo = ">=0.2.0"
 openfermion = "^1.5"
 pyscf = "2.3.0"
 scipy = "^1.10.1"
+furo = "^2023.09.10"
 
 [build-system]
 requires = ["poetry-core"]

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