update README

JoelHBierman · Dec 5, 2023 · 8004ef6 · 8004ef6
1 parent 7bc7834
commit 8004ef6
Show file tree

Hide file tree

Showing 11 changed files with 187 additions and 136 deletions.
diff --git a/README.md b/README.md
@@ -1,37 +1,92 @@
 ![ecosystem](https://raw.githubusercontent.com/qiskit-community/ecosystem/main/badges/qiskit-ecosystem_template.svg)
 
-# Qiskit Ecosystem project template
-- This project is a template project with everything need to be add inside the __Qiskit Ecosystem__.
-- The workflows inside this template are the same than the one running adding and weekly checks.
-
-## Structure
-- .pylinrc : config file for python linter
-- .coveragerc : config file for python coverage
-- tox.ini : file for running tox cmd
-- .gitignore : python gitignore content
-- requirements.txt : install for running project
-- requirements-dev.txt : install for dev the project
-- setup.py : python module/package configuration
-- setup.cfg : python module/package configuration
-- ecosystem.json : config file for Ecosystem testing
-- src/ : repository for all the source of the project
-- tests/ : repository for the units tests of the project
-
-## Tests execution
-- Run for style checks `tox -elint`
-- Run for ecosystem.json checks `tox -ejson`
-- Run for tests `tox -epy39`
-- Run coverage `tox -ecoverage`
-- Run black `tox -eblack`
-  - Run black on unitary file `black <file>`
-
-## Actions
-- The worflow tests.yml is here to tests automatically at every push:
-  - units tests
-  - coverage
-  - lint
-  - black
-- The lastest version of the files `.pylinrc` and `.coveragerc` are automatically download during the actions.
-
-### Qiskit actions
-- The workflows qiskit.yml is here to tests the project are still update with the actual and dev version of Qiskit (if Qiskit is not false in `ecosystem.json`).
+# Electronic Structure Orbital Optimization
+
+Orbital optimization refers to a class of methods in quantum chemistry where the single-particle wavefunctions (the orbitals) are mapped via some parameterized transformation onto a new set of orbitals. An eigensolver is then used to solve for the ground and/or excited states in this transformed basis. The solution of this eigensolver is then used to re-initialize the basis transformation to re-optimize the orbitals. This interleaving of orbital optimization and eigensolver subproblems is repeated until some outer loop stopping condition is reached. This repository contains implementations of the OptOrbVQE algorithm (![arxiv2208.14431](https://arxiv.org/abs/2208.14431v2)) and its excited states variants (![arxiv2310.09418](https://arxiv.org/abs/2310.09418)). In these methods, the orbital transformation takes the form of a partial unitary matrix and quantum eigensolvers such as VQE, SSVQE, MCVQE, and VQD are used in between orbital optimization runs. The starting basis is chosen to be one with a large number $M$ of orbitals (e.g. cc-pVTZ, cc-pVQZ, cc-pV5Z, ect...) and compressing the active space to one of size $N < M$. Thus, the orbital optimization searches for an optimal transformation in the space of $M \times N$ real partial unitaries. Currently, orbital optimization supports the VQE and AdaptVQE ground state solvers and the SSVQE, MCVQE, and VQD excited states solvers. Additional eigensolvers such as QPE, QITE, and VarQITE could be added in the future.
+
+## Installation
+
+To install this project, run:
+
+```
+git clone https://github.com/JoelHBierman/electronic-structure-orbital-optimization.git
+cd electronic-structure-orbital-optimization-main
+pip install .
+```
+
+## OptOrbVQE Code Example
+
+This is an example of how a user might use this orbital optimization scheme in combination with VQE to find the orbital-optimized ground state energy.
+
+```
+import numpy as np
+from qiskit_nature.second_q.drivers import PySCFDriver
+from qiskit_nature.units import DistanceUnit
+from qiskit_nature.second_q.mappers import JordanWignerMapper
+from qiskit_algorithms.minimum_eigensolvers import VQE
+from qiskit_algorithms.optimizers import L_BFGS_B, COBYLA
+from qiskit_aer.primitives import Estimator
+from qiskit_nature.second_q.circuit.library import HartreeFock, UCCSD
+
+from electronic_structure_algorithms.orbital_optimization import PartialUnitaryProjectionOptimizer, OptOrbVQE
+
+estimator = Estimator(approximation=True)
+mapper=JordanWignerMapper()
+
+# Initialize an ElectronicStructureProblem with a large basis.
+# Here we use cc-pVTZ.
+driver = PySCFDriver(atom=f'H 0 0 0; H 0 0 {0.735}',
+                     charge=0,
+                     spin=0,
+                     unit=DistanceUnit.ANGSTROM,
+                     basis='cc-pVTZ')
+
+q_molecule = driver.run()
+num_particles = q_molecule.num_particles
+
+# The size of the active space for OptOrbVQE to use.
+num_reduced_spin_orbitals = 4
+
+# Choose an initialization and an ansatz for VQE.
+# Here we initalize the UCCSD ansatz in the Hartree-Fock state.
+HF_state = HartreeFock(qubit_mapper=mapper,
+                       num_spatial_orbitals=int(num_reduced_spin_orbitals/2),
+                       num_particles=num_particles)
+
+ansatz = UCCSD(qubit_mapper=mapper,
+               num_spatial_orbitals=int(num_reduced_spin_orbitals/2),
+               num_particles=num_particles,
+               initial_state=HF_state)
+
+# Initialize the optimizer for the orbital optimization subroutine.
+partial_unitary_optimizer = PartialUnitaryProjectionOptimizer(initial_BBstepsize=10**-3,
+                                                              stopping_tolerance=10**-5,
+                                                              maxiter=10000,
+                                                              gradient_method='autograd')
+
+# Initialize the VQE solver.
+vqe_instance = VQE(ansatz=ansatz,
+                   initial_point=np.zeros(ansatz.num_parameters),
+                   optimizer=L_BFGS_B(maxfun=10**6, maxiter=10**6),
+                   estimator=estimator)
+
+# Initialize the OptOrbVQE solver.
+optorbvqe_instance = OptOrbVQE(num_spin_orbitals=num_reduced_spin_orbitals,
+                               ground_state_solver=vqe_instance,
+                               mapper=mapper,
+                               estimator=estimator,
+                               partial_unitary_optimizer=partial_unitary_optimizer,
+                               problem=q_molecule,
+                               maxiter=20,
+                               wavefuntion_real=True,
+                               spin_conserving=True)
+
+ground_state_energy_result = optorbvqe_instance.compute_minimum_energy()
+print(f'Orbital-optimized ground state energy: {ground_state_energy_result.eigenvalue} Ha (4 spin-orbitals)')
+print(f'Ground state energy in STO-3G basis: {-1.85727503} Ha (4 spin-orbitals)')
+print(f'Ground state energy in cc-pVTZ basis: {-1.89226657} Ha (56 spin-orbitals)')
+
+```
+
+We can see that OptOrbVQE greatly improves upon the VQE ground state energy, despite using the same number of qubits, however there is still room to improve, as evidenced by the fact that the energy in the full cc-pVTZ active space is still much more accurate. The accuracy of OptOrbVQE can be systematically improved
+by increasing the value of `num_reduced_spin_orbitals` at the cost of more qubits being required.
diff --git a/electronic_structure_algorithms/orbital_optimization/base_opt_orb_solver.py b/electronic_structure_algorithms/orbital_optimization/base_opt_orb_solver.py
@@ -17,12 +17,12 @@
 class BaseOptOrbSolver():
 
     def __init__(self,
-        problem: Optional[ElectronicStructureProblem],
-        integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]],
         num_spin_orbitals: int,
         mapper: QubitMapper,
         estimator: BaseEstimator,
         partial_unitary_optimizer: PartialUnitaryProjectionOptimizer,
+        problem: Optional[ElectronicStructureProblem] = None,
+        integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]] = None,
         initial_partial_unitary: Optional[Union[torch.Tensor, np.ndarray]] = None,
         maxiter: int = 10,
         stopping_tolerance: float = 10**-5,
@@ -35,14 +35,14 @@ def __init__(self,
         """
         
         Args:
-            problem: An ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained.
-            integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first
-                entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in
-                the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`.
             num_spin_orbitals: The number of spin-orbitals to use for the active space.
             mapper: A QubitMapper to use for the RDM calculations.
             partial_unitary_optimizer: An instance of :class:`PartialUnitaryProjectionOptimizer` to use for
                 the basis optimization subroutine.
+            problem: An ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained.
+            integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first
+                entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in
+                the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`.
             initial_partial_unitary: The initial guess for the orbital rotation matrix. If ``None``, then a permutation matrix
                 selecting the spatial orbitals with the lowest energy will be generated.
             maxiter: The maximum number of outerloop iterations. (The number of times the wavefunction optimization is run.)

diff --git a/electronic_structure_algorithms/orbital_optimization/opt_orb_adapt_vqe.py b/electronic_structure_algorithms/orbital_optimization/opt_orb_adapt_vqe.py
@@ -13,13 +13,13 @@
 class OptOrbAdaptVQE(OptOrbMinimumEigensolver):
 
     def __init__(self,
-        problem: Optional[ElectronicStructureProblem],
-        integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]],
         num_spin_orbitals: int,
         ground_state_solver: MinimumEigensolver,
         mapper: QubitMapper,
         estimator: BaseEstimator,
         partial_unitary_optimizer: PartialUnitaryProjectionOptimizer,
+        problem: Optional[ElectronicStructureProblem] = None,
+        integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]] = None,
         initial_partial_unitary: Optional[Union[torch.Tensor, np.ndarray]] = None,
         maxiter: int = 10,
         stopping_tolerance: float = 10**-5,
@@ -32,14 +32,14 @@ def __init__(self,
         """
         
         Args:
-            problem: The ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained.
-            integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first
-                entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in
-                the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`.
             num_spin_orbitals: The number of spin-orbitals to use for the active space.
             ground_state_solver: An instance of VQE to use for the wavefunction optimization.
             mapper: A QubitMapper to use for the RDM calculations.
             partial_unitary_optimizer: An instance of PartialUnitaryProjectionOptimizer to use for the basis optimization.
+            problem: The ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained.
+            integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first
+                entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in
+                the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`.
             initial_partial_unitary: The initial guess for the orbital rotation matrix. If ``None``, then a permutation matrix
                 selecting the spatial orbitals with the lowest energy will be generated.
             maxiter: The maximum number of outerloop iterations. (The number of times the wavefunction optimization is run.)
@@ -63,13 +63,13 @@ def __init__(self,
 
         """
 
-        super().__init__(problem=problem,
-                         integral_tensors=integral_tensors,
-                         num_spin_orbitals=num_spin_orbitals,
+        super().__init__(num_spin_orbitals=num_spin_orbitals,
                          ground_state_solver=ground_state_solver,
                          mapper=mapper,
                          estimator=estimator,
                          partial_unitary_optimizer=partial_unitary_optimizer,
+                         problem=problem,
+                         integral_tensors=integral_tensors,
                          initial_partial_unitary=initial_partial_unitary,
                          maxiter=maxiter,
                          stopping_tolerance=stopping_tolerance,

diff --git a/electronic_structure_algorithms/orbital_optimization/opt_orb_eigensolver.py b/electronic_structure_algorithms/orbital_optimization/opt_orb_eigensolver.py
@@ -14,13 +14,13 @@
 class OptOrbEigensolver(BaseOptOrbSolver):
 
     def __init__(self,
-        problem: Optional[ElectronicStructureProblem],
-        integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]],
         num_spin_orbitals: int,
         excited_states_solver: Eigensolver,
         mapper: QubitMapper,
         estimator: BaseEstimator,
         partial_unitary_optimizer: PartialUnitaryProjectionOptimizer,
+        problem: Optional[ElectronicStructureProblem] = None,
+        integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]] = None,
         initial_partial_unitary: Optional[Union[torch.Tensor, np.ndarray]] = None,
         maxiter: int = 10,
         stopping_tolerance: float = 10**-5,
@@ -34,14 +34,15 @@ def __init__(self,
         """
         
         Args:
+            num_spin_orbitals: The number of spin-orbitals to use for the active space.
+            excited_states_solver: An instance of :class:`Eigensolver` for the excited states solver subroutine.
+            mapper: A QubitMapper to use for the RDM calculations.
+            partial_unitary_optimizer: An instance of PartialUnitaryProjectionOptimizer to use for the basis optimization.
             problem: The ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained.
             integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first
                 entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in
                 the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`.
             num_spin_orbitals: The number of spin-orbitals to use for the active space.
-            excited_states_solver: An instance of :class:`Eigensolver` for the excited states solver subroutine.
-            mapper: A QubitMapper to use for the RDM calculations.
-            partial_unitary_optimizer: An instance of PartialUnitaryProjectionOptimizer to use for the basis optimization.
             initial_partial_unitary: The initial guess for the orbital rotation matrix. If ``None``, then a permutation matrix
                 selecting the spatial orbitals with the lowest energy will be generated.
             maxiter: The maximum number of outerloop iterations. (The number of times the wavefunction optimization is run.)
@@ -63,12 +64,12 @@ def __init__(self,
                 objective function whose length is equal to the number of excited states.
 
         """
-        super().__init__(problem=problem,
-                         integral_tensors=integral_tensors,
-                         num_spin_orbitals=num_spin_orbitals,
+        super().__init__(num_spin_orbitals=num_spin_orbitals,
                          mapper=mapper,
                          estimator=estimator,
                          partial_unitary_optimizer=partial_unitary_optimizer,
+                         problem=problem,
+                         integral_tensors=integral_tensors,
                          initial_partial_unitary=initial_partial_unitary,
                          maxiter=maxiter,
                          stopping_tolerance=stopping_tolerance,

diff --git a/electronic_structure_algorithms/orbital_optimization/opt_orb_mcvqe.py b/electronic_structure_algorithms/orbital_optimization/opt_orb_mcvqe.py
@@ -13,13 +13,13 @@
 class OptOrbMCVQE(OptOrbEigensolver):
 
     def __init__(self,
-        problem: ElectronicStructureProblem,
-        integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]],
         num_spin_orbitals: int,
         excited_states_solver: Eigensolver,
         mapper: QubitMapper,
         estimator: BaseEstimator,
         partial_unitary_optimizer: PartialUnitaryProjectionOptimizer,
+        problem: ElectronicStructureProblem = None,
+        integral_tensors: Optional[Union[tuple[torch.Tensor, torch.Tensor], tuple[np.ndarray, np.ndarray]]] = None,
         initial_partial_unitary: Optional[Union[torch.Tensor, np.ndarray]] = None,
         maxiter: int = 10,
         stopping_tolerance: float = 10**-5,
@@ -33,14 +33,14 @@ def __init__(self,
         """
         
         Args:
-            problem: The ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained.
-            integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first
-                entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in
-                the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`.
             num_spin_orbitals: The number of spin-orbitals to use for the active space.
             excited_states_solver: An instance of VQE to use for the wavefunction optimization.
             mapper: A QubitMapper to use for the RDM calculations.
             partial_unitary_optimizer: An instance of PartialUnitaryProjectionOptimizer to use for the basis optimization.
+            problem: The ElectronicStructureProblem from which molecule information such as one and two body integrals is obtained.
+            integral_tensors: A tuple storing the one and two body integrals of the full orbital space. The first
+                entry stores the one body integrals and the second entry stores the two body integrals in the Physicist's notation in
+                the dense spin-orbital representation, represented as either :class:`torch.Tensor` or :class:`np.ndarray`.
             initial_partial_unitary: The initial guess for the orbital rotation matrix. If ``None``, then a permutation matrix
                 selecting the spatial orbitals with the lowest energy will be generated.
             maxiter: The maximum number of outerloop iterations. (The number of times the wavefunction optimization is run.)
@@ -64,13 +64,13 @@ def __init__(self,
 
         """
 
-        super().__init__(problem=problem,
-                         integral_tensors=integral_tensors,
-                         num_spin_orbitals=num_spin_orbitals,
+        super().__init__(num_spin_orbitals=num_spin_orbitals,
                          excited_states_solver=excited_states_solver,
                          mapper=mapper,
                          estimator=estimator,
                          partial_unitary_optimizer=partial_unitary_optimizer,
+                         problem=problem,
+                         integral_tensors=integral_tensors,
                          initial_partial_unitary=initial_partial_unitary,
                          maxiter=maxiter,
                          stopping_tolerance=stopping_tolerance,