diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 4f5120c..cba01e0 100644 --- a/.pre-commit-config.yaml +++ b/.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 diff --git a/doc/source/api-reference/ansatz.rst b/doc/source/api-reference/ansatz.rst index daea608..ece5596 100644 --- a/doc/source/api-reference/ansatz.rst +++ b/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 diff --git a/doc/source/conf.py b/doc/source/conf.py index 768a7bd..e84d6fc 100644 --- a/doc/source/conf.py +++ b/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 diff --git a/examples/br_example.py b/examples/br_example.py new file mode 100644 index 0000000..065f362 --- /dev/null +++ b/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) diff --git a/examples/h3p.xyz b/examples/h3p.xyz new file mode 100644 index 0000000..6353a63 --- /dev/null +++ b/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 diff --git a/examples/ucc_example1.py b/examples/ucc_example1.py index 4529b30..bb453f8 100644 --- a/examples/ucc_example1.py +++ b/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 diff --git a/examples/ucc_example2.py b/examples/ucc_example2.py index 67ba4cf..4241876 100644 --- a/examples/ucc_example2.py +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index 824591b..49a9f3e 100644 --- a/pyproject.toml +++ b/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"] diff --git a/src/qibochem/ansatz/basis_rotation.py b/src/qibochem/ansatz/basis_rotation.py new file mode 100644 index 0000000..4a35dd7 --- /dev/null +++ b/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 diff --git a/src/qibochem/ansatz/hardware_efficient.py b/src/qibochem/ansatz/hardware_efficient.py index b290226..c84c96a 100644 --- a/src/qibochem/ansatz/hardware_efficient.py +++ b/src/qibochem/ansatz/hardware_efficient.py @@ -1,26 +1,31 @@ from qibo import gates -def hea_su2(nlayers, nqubits): - """ - implements a hardware-efficient ansatz SU2 - inputs: +def hea(n_layers, n_qubits, parameter_gates=["RY", "RZ"], coupling_gates="CZ"): + """Builds the generalized hardware-efficient ansatz, i.e. the choice of rotation and entangling gates can be + manually defined + + Args: + n_layers: Number of layers of rotation and entangling gates + n_qubits: Number of qubits in the quantum circuit + parameter_gates: List of single-qubit rotation gates to be used in the ansatz. The gates should be given as + strings representing valid ``Qibo`` one-qubit gates. Default: ``["RY", "RZ"]`` + coupling_gates: String representing the two-qubit entangling gate to be used in the ansatz; should be a + valid two-qubit ``Qibo`` gate. Default: ``"CZ"`` + + Returns: + List of gates corresponding to the hardware-efficient ansatz """ gate_list = [] - for ilayer in range(nlayers): - for iqubit in range(nqubits): - gate_list.append(gates.RY(iqubit, theta=0)) - for iqubit in range(nqubits): - gate_list.append(gates.RZ(iqubit, theta=0)) - # entanglement - for iqubit in range(nqubits - 1): - gate_list.append(gates.CNOT(iqubit, iqubit + 1)) - gate_list.append(gates.CNOT(nqubits - 1, 0)) + for ilayer in range(n_layers): + for rgate in parameter_gates: + for iqubit in range(n_qubits): + gate_list.append(getattr(gates, rgate)(iqubit, theta=0)) - for iqubit in range(nqubits): - gate_list.append(gates.RY(iqubit, theta=0)) - for iqubit in range(nqubits): - gate_list.append(gates.RZ(iqubit, theta=0)) + # entanglement + for iqubit in range(n_qubits - 1): + gate_list.append(getattr(gates, coupling_gates)(iqubit, iqubit + 1)) + gate_list.append(getattr(gates, coupling_gates)(n_qubits - 1, 0)) return gate_list diff --git a/src/qibochem/ansatz/hf_reference.py b/src/qibochem/ansatz/hf_reference.py index 8e9d846..5683a50 100644 --- a/src/qibochem/ansatz/hf_reference.py +++ b/src/qibochem/ansatz/hf_reference.py @@ -51,16 +51,15 @@ def bk_matrix(n: int): # Main function def hf_circuit(n_qubits, n_electrons, ferm_qubit_map=None): - """Circuit for a Hartree-Fock reference + """Circuit to prepare a Hartree-Fock state Args: - n_qubits: Number of qubits in circuit - n_electrons: Number of electrons in molecular system - ferm_qubit_map: Fermion to qubit mapping; Jordan-Wigner ('jw') or Brayvi-Kitaev ('bk') - Defaults to 'jw' + n_qubits: Number of qubits in the quantum circuit + n_electrons: Number of electrons in the molecular system + ferm_qubit_map: Fermion to qubit map. Must be either Jordan-Wigner (``jw``) or Brayvi-Kitaev (``bk``). Default value is ``jw``. Returns: - Qibo Circuit initialized in a HF reference state + Qibo ``Circuit`` initialized in a HF reference state """ # Which fermion-to-qubit map to use if ferm_qubit_map is None: diff --git a/src/qibochem/ansatz/ucc.py b/src/qibochem/ansatz/ucc.py index a62f776..7196269 100644 --- a/src/qibochem/ansatz/ucc.py +++ b/src/qibochem/ansatz/ucc.py @@ -58,19 +58,21 @@ def expi_pauli(n_qubits, theta, pauli_string): def ucc_circuit(n_qubits, excitation, theta=0.0, trotter_steps=1, ferm_qubit_map=None, coeffs=None): - """ - Build a circuit corresponding to the unitary coupled-cluster ansatz for only one excitation + r""" + Circuit corresponding to the unitary coupled-cluster ansatz for a single excitation Args: - n_qubits: No. of qubits in the quantum circuit - excitation: List of orbitals involved in the excitation; must have even number of elements. - E.g. [0, 1, 2, 3] represents the excitation of electrons in orbitals (0, 1) to (2, 3) + n_qubits: Number of qubits in the quantum circuit + excitation: Iterable of orbitals involved in the excitation; must have an even number of elements + E.g. ``[0, 1, 2, 3]`` represents the excitation of electrons in orbitals ``(0, 1)`` to ``(2, 3)`` theta: UCC parameter. Defaults to 0.0 - trotter_steps: number of Trotter steps; i.e. number of times the UCC ansatz is applied with - theta=theta/trotter_steps - ferm_qubit_map: fermion-to-qubit transformation. Default is Jordan-Wigner (jw) + trotter_steps: Number of Trotter steps; i.e. number of times the UCC ansatz is applied with \theta = \theta / trotter_steps + ferm_qubit_map: Fermion-to-qubit transformation. Default is Jordan-Wigner (``jw``). coeffs: List to hold the coefficients for the rotation parameter in each Pauli string. May be useful in running the VQE. WARNING: Will be modified in this function + + Returns: + Qibo ``Circuit`` corresponding to a single UCC excitation """ # Check size of orbitals input n_orbitals = len(excitation) diff --git a/src/qibochem/driver/molecule.py b/src/qibochem/driver/molecule.py index 48dfe61..fca61fd 100644 --- a/src/qibochem/driver/molecule.py +++ b/src/qibochem/driver/molecule.py @@ -20,35 +20,19 @@ class Molecule: """ Class representing a single molecule - :ivar geometry(list): molecular coordinates in OpenFermion format - e.g. [('H', (0.0, 0.0, 0.0)), ('H', (0.0, 0.0, 0.7))] - :ivar charge(int): net charge of molecule - :ivar multiplicity(int): spin multiplicity of molecule - :ivar basis(str): atomic orbital basis set to be used for PySCF calculation - :ivar xyz_file(str): xyz file containing molecular coordinates in chemical xyz format, - comment line should include charge, multiplicity values - :ivar active(list): iterable representing the set of MOs to be included in - quantum simulation, e.g. list(range(3,6)) for active space of orbitals 3,4,5. + Args: + geometry (list): Molecular coordinates in OpenFermion format, e.g. ``[('H', (0.0, 0.0, 0.0)), ('H', (0.0, 0.0, 0.7))]`` + charge (int): Net charge of molecule + multiplicity (int): Spin multiplicity of molecule + basis (str): Atomic orbital basis set, used for the PySCF/PSI4 calculations. Default: "STO-3G" (minimal basis) + xyz_file (str): .xyz file containing the molecular coordinates. The comment line should follow + "{charge} {multiplicity}" + active: Iterable representing the set of MOs to be included in the quantum simulation + e.g. ``list(range(3,6))`` for an active space with orbitals 3, 4 and 5. """ def __init__(self, geometry=None, charge=0, multiplicity=1, basis=None, xyz_file=None, active=None): - """ - Args: - geometry (list): Molecular coordinates in OpenFermion format - e.g. [('H', (0.0, 0.0, 0.0)), ('H', (0.0, 0.0, 0.7))] - charge (int): Net charge of molecule - multiplicity (int): Spin multiplicity of molecule - basis (str): Atomic orbital basis set, used for the PySCF/PSI4 calculations - xyz_file (str): .xyz file containing the molecular coordinates - Comment line should follow "{charge} {multiplicity}" - active: Iterable representing the set of MOs to be included in the quantum simulation - e.g. list(range(3,6)) for orbitals 3, 4, 5 active space. - - Example: - .. testcode:: - TODO - """ # Basic properties # Define using the function arguments if xyz_file not given if xyz_file is None: @@ -58,7 +42,7 @@ def __init__(self, geometry=None, charge=0, multiplicity=1, basis=None, xyz_file else: # Check if xyz_file exists, then fill in the Molecule attributes assert Path(f"{xyz_file}").exists(), f"{xyz_file} not found!" - self.process_xyz_file(xyz_file) + self._process_xyz_file(xyz_file) if basis is None: # Default bais is STO-3G self.basis = "sto-3g" @@ -96,7 +80,7 @@ def __init__(self, geometry=None, charge=0, multiplicity=1, basis=None, xyz_file self.n_active_e = None self.n_active_orbs = None # Number of spin-orbitals in the active space - def process_xyz_file(self, xyz_file): + def _process_xyz_file(self, xyz_file): """ Reads a .xyz file to obtain and set the molecular coordinates (in OpenFermion format), charge, and multiplicity @@ -128,7 +112,7 @@ def run_pyscf(self, max_scf_cycles=50): molecular integrals Args: - basis: Atomic orbital basis set + max_scf_cycles: Maximum number of SCF cycles in PySCF """ import pyscf @@ -183,8 +167,8 @@ def run_psi4(self, output=None): molecular integrals Args: - output: Name of PSI4 output file. None suppresses the output on non-Windows systems, - and uses 'psi4_output.dat' otherwise + output: Name of PSI4 output file. ``None`` suppresses the output on non-Windows systems, + and uses ``psi4_output.dat`` otherwise """ import psi4 # pylint: disable=import-error @@ -251,7 +235,7 @@ def run_psi4(self, output=None): self.da = self.ca.T @ self.overlap @ self.pa @ self.overlap @ self.ca # HF embedding functions - def inactive_fock_matrix(self, frozen): + def _inactive_fock_matrix(self, frozen): """ Returns the full inactive Fock matrix @@ -273,8 +257,8 @@ def inactive_fock_matrix(self, frozen): def hf_embedding(self, active=None, frozen=None): """ - Turns on HF embedding for a given active/frozen space, i.e. - fills the class attributes: inactive_energy, embed_oei, and embed_tei + Turns on HF embedding for a given active/frozen space, and fills in the class attributes: + ``inactive_energy``, ``embed_oei``, and ``embed_tei``. Args: active: Iterable representing the active-space for quantum simulation @@ -301,7 +285,7 @@ def hf_embedding(self, active=None, frozen=None): ) # Build the inactive Fock matrix first - inactive_fock = self.inactive_fock_matrix(frozen) + inactive_fock = self._inactive_fock_matrix(frozen) # Calculate the inactive Fock energy # Only want frozen part of original OEI and inactive Fock matrix @@ -331,22 +315,15 @@ def hamiltonian( Builds a molecular Hamiltonian using the one-/two- electron integrals Args: - ham_type: Format of molecular Hamiltonian returned - - ("f", "ferm"): OpenFermion FermionOperator - - ("q", "qubit"): OpenFermion QubitOperator - - ("s", "sym"): Qibo SymbolicHamiltonian (default) - - oei: 1-electron integrals. Default: self.oei (MO basis) - - tei: 2-electron integrals in 2ndQ notation. Default: self.tei (MO basis) - + ham_type: Format of molecular Hamiltonian returned. The available options are: + ``("f", "ferm")``: OpenFermion ``FermionOperator``, + ``("q", "qubit")``: OpenFermion ``QubitOperator``, or + ``("s", "sym")``: Qibo ``SymbolicHamiltonian`` (default) + oei: 1-electron integrals. Default: ``self.oei`` (MO basis) + tei: 2-electron integrals in 2ndQ notation. Default: ``self.tei`` (MO basis) constant: For inactive Fock energy if embedding used. Default: 0.0 - ferm_qubit_map: Which fermion to qubit transformation to use. - Must be either "jw" (default) or "bk" + Must be either ``jw`` (default) or ``bk`` Returns: Molecular Hamiltonian in the format of choice @@ -389,11 +366,11 @@ def hamiltonian( def eigenvalues(hamiltonian): """ Finds the lowest 6 exact eigenvalues of the molecular Hamiltonian - Note: Use the .eigenvalues() method for a Qibo SymbolicHamiltonian object + Note: Uses the ``eigenvalues()`` class method for a Qibo ``SymbolicHamiltonian`` object Args: - hamiltonian: Molecular Hamiltonian, given as a FermionOperator, QubitOperator, or - SymbolicHamiltonian (not recommended) + hamiltonian: Molecular Hamiltonian, given as a ``FermionOperator``, ``QubitOperator``, or + ``SymbolicHamiltonian`` (not recommended) """ if isinstance(hamiltonian, (openfermion.FermionOperator, openfermion.QubitOperator)): from scipy.sparse import linalg diff --git a/tests/test_basis_rotation.py b/tests/test_basis_rotation.py new file mode 100644 index 0000000..c1de41c --- /dev/null +++ b/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) diff --git a/tests/test_hardware_efficient.py b/tests/test_hardware_efficient.py new file mode 100644 index 0000000..d202904 --- /dev/null +++ b/tests/test_hardware_efficient.py @@ -0,0 +1,54 @@ +import numpy as np +import pytest +from qibo import gates, models +from scipy.optimize import minimize + +from qibochem.ansatz import hardware_efficient +from qibochem.driver.molecule import Molecule +from qibochem.measurement.expectation import expectation + + +def test_hea_ansatz(): + mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.74804))]) + mol.run_pyscf() + mol_classical_hf_energy = mol.e_hf + mol_sym_ham = mol.hamiltonian("s") + + nlayers = 1 + nqubits = mol.nso + ntheta = 2 * nqubits * nlayers + theta = np.zeros(ntheta) + + hea_ansatz = hardware_efficient.hea(nlayers, nqubits) + qc = models.Circuit(nqubits) + qc.add(gates.X(_i) for _i in range(sum(mol.nelec))) + qc.add(hea_ansatz) + qc.set_parameters(theta) + + hf_energy = expectation(qc, mol_sym_ham) + assert mol_classical_hf_energy == pytest.approx(hf_energy) + + +def test_vqe_hea_ansatz(): + def test_vqe_hea_ansatz_cost(parameters, circuit, hamiltonian): + circuit.set_parameters(parameters) + return expectation(circuit, hamiltonian) + + mol = Molecule([("H", (0.0, 0.0, 0.0)), ("H", (0.0, 0.0, 0.75))]) + mol.run_pyscf() + mol_classical_hf_energy = mol.e_hf + mol_sym_ham = mol.hamiltonian("s") + + nlayers = 3 + nqubits = mol.nso + ntheta = 2 * nqubits * nlayers + theta = np.full(ntheta, np.pi / 4) + + hea_ansatz = hardware_efficient.hea(nlayers, nqubits, ["RY", "RZ"], "CNOT") + qc = models.Circuit(nqubits) + qc.add(gates.X(_i) for _i in range(sum(mol.nelec))) + qc.add(hea_ansatz) + qc.set_parameters(theta) + + vqe_object = minimize(test_vqe_hea_ansatz_cost, theta, args=(qc, mol_sym_ham), method="Powell") + assert vqe_object.fun == pytest.approx(-1.1371170019838128)