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/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/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)