Help with refactoring / migrating OO-VQE

@mrossinek May I know where this energy library originates from ?

Traceback (most recent call last):
  Cell In[29], line 39
    oovqe_result = oovqe.solve(driver)
  Cell In[28], line 463 in solve
    self._initialize_additional_parameters(driver)
  Cell In[28], line 376 in _initialize_additional_parameters
    self._set_operator_and_vqe(driver)
  Cell In[28], line 319 in _set_operator_and_vqe
    operator, aux_operators = self._do_transform(self._qmolecule, self._transformation)
  Cell In[28], line 99 in _do_transform
    self._hf_energy = qmolecule.hf_energy
AttributeError: 'ElectronicStructureProblem' object has no attribute 'hf_energy'

    def _set_operator_and_vqe(self, driver: BaseDriver):
        """ Initializes the operators using provided driver of qmolecule."""

        if not isinstance(self._transformation, FermionicOp):
            raise QiskitError('OrbitalOptimizationVQE requires a FermionicTransformation.')
        if not isinstance(driver, PySCFDriver):
            raise QiskitError('OrbitalOptimizationVQE only works with PySCF Drivers.')

        if self._qmolecule is None:
            # in future, self._transformation.transform should return also qmolecule
            # to avoid running the driver twice
            self._qmolecule = driver.run()
            operator, aux_operators = self._do_transform(self._qmolecule, self._transformation)
        else:
            operator, aux_operators = self._do_transform(self._qmolecule, self._transformation)

Hi @buttercutter, please excuse the late reply.

The algorithms module of Qiskit is currently undergoing some changes (see also https://github.com/Qiskit/RFCs/pull/44) which leaves the state of Qiskit Nature's algorithms still a bit uncertain. That is why your timing to work on this is not ideal and why I had not replied earlier.

I also saw your question here (https://github.com/qiskit-community/qiskit-nature/pull/629#issuecomment-1616626612) and just noticed that I never replied. In fact, you can do orbital rotation with the PySCF plugin but you need to use the CASSCF solver rather than the CASCI one (which is the default in the demo). The unittests also use CASSCF though whose API is basically identical to that of the CASCI solver (this latter one does not perform orbital rotation as you had noticed).

The existence of this CASSCF-integration via PySCF is also the main reason I do not see an advantage of an ooVQE implementation.

That said, let me at least reply to some questions you raised here.

I need some technical help on migrating _do_transform() code since it involves dipole and spin orbital

I cannot access this Google Colab link. If you have a specific code question, I suggest you upload a minimal Gist (not the massive Jupyter notebook of your first link). Alternatively, you can also embed a small Python script directly in a comment here.

Relevant issues : #629 (review) , FermionicTransformation class replacement by QubitConverter class

The QubitConverter recently got deprecated in favor of the refactored QubitMapper classes (see this issue https://github.com/qiskit-community/qiskit-nature/issues/967 and the various PRs linking to it). Furthermore, the _do_transform method which you have been seeing refers back to a FermionicTransformation as it existed back in Qiskit Aqua (see here https://github.com/qiskit-community/qiskit-aqua/blob/c1564af8792c6664670807614a378147fd04d28f/qiskit/chemistry/transformations/fermionic_transformation.py). This object was refactored and removed prior to Qiskit Nature's first release (0.1). It used to have the job of mapping (now done via the QubitMapper objects), active space reduction (now done via the ActiveSpaceTransformer or FreezeCoreTransformer) and particle hole transformation (not implemented in Nature, see #48). You will need to investigate the original code of ooVQE in detail to figure out how to replace this objects usage. It was built around the QMolecule which no longer exists (see #148 and the various PRs linking to it).

May I know where this energy library originates from ?

What you are looking at here is code belonging to the Psi4Driver (it reads this template file). As a result, this energy function comes from Psi4. That said, your ooVQE implementation should not have to rely on any driver-related code. These drivers are in my opinion a legacy interface of Qiskit Nature to classical computational chemistry codes which I would like to replace with plugins (in the style of https://github.com/qiskit-community/qiskit-nature-pyscf) in the long run. As of now, they still serve a nice and easy entry point for educational and testing purposes, but plugins are superior in their functionality and flexibility.

The ooVQE implementation should follow the interface set forth by the GroundStateSolver (https://github.com/qiskit-community/qiskit-nature/blob/c8207e39f33ea7a692d1c3ce0d2509d124f3bd80/qiskit_nature/second_q/algorithms/ground_state_solvers/ground_state_solver.py). Thus, all you will need is an ElectronicStructureProblem (https://github.com/qiskit-community/qiskit-nature/blob/c8207e39f33ea7a692d1c3ce0d2509d124f3bd80/qiskit_nature/second_q/problems/electronic_structure_problem.py) as the input to .solve() which you can obtain from a driver as the output of its .run() method.

I will deal with ooVQE issue with do_transform() later after the following runtime error is resolved.

The following code is the use of CASSCF library to rotate the orbital, feel free to correct me if the code is wrong / misses anything.

I am having trouble loading the resulting orbital-rotated hamiltonian into qiskit-nature. Please advise.

"""# **Orbital rotation using pySCF CASSCF library**

In the context of the CASSCF (Complete Active Space Self-Consistent Field) method, "rotate_orb_cc" and "rotate_mo" are two different functions used to rotate the molecular orbitals (MOs) during the CASSCF optimization process.

"rotate_orb_cc" stands for "rotate orbitals using coupled-cluster theory". This function performs a rotation of the MOs using a coupled-cluster approach, which involves solving for the amplitudes of the excitation operators that describe the correlation effects beyond the Hartree-Fock approximation. This method is generally more computationally expensive but can lead to more accurate results.

On the other hand, "rotate_mo" simply rotates the MOs using a diagonalization method. This method is less computationally expensive but may be less accurate than the coupled-cluster approach.

In summary, the main difference between "rotate_orb_cc" and "rotate_mo" is the method used to rotate the MOs during the CASSCF optimization process. "rotate_orb_cc" uses a coupled-cluster approach, while "rotate_mo" uses a simpler diagonalization method.

Credit : https://github.com/qiskit-community/qiskit-nature/issues/1219
"""

from pyscf import gto, scf, mcscf

# Set up molecule
mol = gto.M(atom='H 0 0 0; F 0 0 1.1', basis='cc-pvdz')

# SCF calculation
mf = scf.RHF(mol)
mf.kernel()

# CASSCF calculation
mc = mcscf.CASSCF(mf, 4, 4)
mc.kernel()

"""After CASSCF calculation, get the arguments needed

rotate_orb_cc() needs:

      mo: MO coefficients
      fcivec: CI vector
      fcasdm1: Spin-separated 1-particle RDM
      fcasdm2: Spin-separated 2-particle RDM
      eris: transformed electron repulsion integrals

which are generated from the CASSCF calculation. Passing those in allows rotate_orb_cc() to do the coupled-cluster orbital rotation.
"""

mo = mc.mo_coeff
fcivec = mc.ci
fcasdm1, fcasdm2 = mc.fcisolver.make_rdm12(fcivec, mc.ncas, mc.nelecas)
eris = mc.ao2mo(mo)

# Print MO coefficients before rotation
print('MO coefficients before:')
print(mc.mo_coeff)

# Rotate orbitals using coupled-cluster (rotate_orb_cc)
mc.rotate_orb_cc(mo, fcivec, fcasdm1, fcasdm2, eris)

# Print MO coefficients after rotation
print('MO coefficients after:')
print(mc.mo_coeff)

"""rotate_mo() needs:

      mc.mo_coeff: Current MO coefficients
      u: Rotation matrix
"""

# CASSCF calculation
mc = mcscf.CASSCF(mf, 4, 4)
mc.kernel()

# Get MO coefficients
mo_coeff = mc.mo_coeff

# Make CI-driven 1-RDM and 2-RDM
fcivec = mc.ci
casdm1, casdm2 = mc.fcisolver.make_rdm12(fcivec, mc.ncas, mc.nelecas)

# Pass as a tuple
rdm = (casdm1, casdm2)

# Compute orbital rotation gradient
g_orb = mc.get_grad(mo_coeff, rdm)
dx = g_orb[0]

# Get U matrix with the gradient
u = mc.update_rotate_matrix(dx)

# Rotate orbitals using diagonalization (rotate_mo)
mc.rotate_mo(mc.mo_coeff, u)

"""Here is one way to take the rotated CASSCF orbitals from pySCF and create a Hamiltonian file for use in Qiskit Nature:

1. After rotating the orbitals in pySCF, export the one- and two-electron integrals:

"""

# In pySCF

h1 = mc.get_hcore() # One-electron integrals
h2 = mc.get_h2eff(mo_coeff) # Two-electron integrals

"""2. Convert the integrals to Numpy arrays:"""

h1 = np.array(h1)
h2 = np.array(h2)

# Convert numpy arrays to lists
h1 = h1.tolist()
h2 = h2.tolist()

"""3. Create a Python dictionary containing the Hamiltonian parameters:"""

ham_dict = {
  'one_body_integrals': h1,
  'two_body_integrals': h2
}

"""4. Save the dictionary to a Hamiltonian file:"""

from qiskit_nature.second_q.hamiltonians import Hamiltonian

from qiskit_nature.second_q.operators import FermionicOp
from qiskit_nature.second_q.mappers import JordanWignerMapper

class MolecularHamiltonian(Hamiltonian):
    def __init__(self, one_body_integrals, two_body_integrals, basis):
        self.one_body_integrals = one_body_integrals
        self.two_body_integrals = two_body_integrals
        self.basis = basis

    def second_q_op(self):
        fo = FermionicOp()
        for orbital1, integral in self.one_body_integrals.items():
            fo += integral * (f'^{orbital1} a_{orbital1}')
        for orbital1, orbital2, integral in self.two_body_integrals.items():
            fo += integral * (f'^{orbital1} ^{orbital2} a_{orbital1} a_{orbital2}')
        # Convert FermionicOp to SparseLabelOp using Jordan-Wigner transform
        return JordanWignerMapper().transform(fo)

    @property
    def register_length(self):
        max_orbital = max(orbital for orbital in self.one_body_integrals.keys())
        max_orbital = max(max_orbital, *self.two_body_integrals.keys())
        # Add 1 to account for zero orbital
        return max_orbital + 1

# Build Hamiltonian from dict
hamiltonian = MolecularHamiltonian(
   one_body_integrals=ham_dict['one_body_integrals'],
   two_body_integrals=ham_dict['two_body_integrals'],
   basis='sto-3g'
)

"""5. In Qiskit Nature, load the Hamiltonian file:"""

from qiskit_nature.second_q.drivers import PySCFDriver

driver = PySCFDriver(atom='H 0 0 0; F 0 0 1.1',
                     basis='sto-3g',
                     charge=0,
                     spin=0,
                     hamiltonian=hamiltonian)

# Set Hamiltonian property
driver.set_property(Energies([hamiltonian]))

If you are initially working with PySCF directly, the PySCFDriver is not the object you will want to use to interface with Qiskit. It is written to run its own separate calculation using PySCF and, as such, does not provide you with the means to load an arbitrary Hamiltonian.

That is exactly, what the plugin https://github.com/qiskit-community/qiskit-nature-pyscf has been developed for. You can see an example on how this works with CASSCF here: https://github.com/qiskit-community/qiskit-nature-pyscf/blob/4969f3b250dd02b9986d7ff6ea97d65ed7e0de9a/test/test_qiskit_solver.py#L74-L93

Regarding your errors:

The error you see states that hamiltonian is not a valid keyword argument to the PySCFDriver class. This is exactly due to my earlier statement.
Also the set_property which you are trying to use does not exist. Again, this points to the same problem: you are trying to use the PySCFDriver in a wrong way.

If all you want is to use a different set of orbitals (the fully rotated once from the result of a prior CASSCF calculation) you can do something similar to how it is done in the plugin repo: https://github.com/qiskit-community/qiskit-nature-pyscf/blob/4969f3b250dd02b9986d7ff6ea97d65ed7e0de9a/qiskit_nature_pyscf/qiskit_solver.py#L156-L183

I would like to take this moment to ask whether you are developing a purely Qiskit Nature-based implementation of the ooVQE algorithm which does not require PySCF? If not and you are only looking for support, I suggest that the quantum computing stackexchange is a better place to get support. The Github issue tracker is not a support system and should be used for bug reports, feature requests, or code design discussions.

I can help with developing a purely Qiskit Nature-based implementation of the ooVQE algorithm which does not require PySCF.

Let me know when I should start doing so for _do_transform() given all the transitions of class and API changes.

@mrossinek Looking at nvidia's cuda-quantum, their quantum chemistry code is still using pySCF library.

Feel free to correct me if wrong / miss anything.

Let me know when I should start doing so for _do_transform() given all the transitions of class and API changes.

No more transitions or class API changes are planned at the moment, so this should not block you any longer.

Looking at nvidia's cuda-quantum, their quantum chemistry code is still using pySCF library.

I do not see how this is relevant. If you want an ooVQE implementation to be merged into Qiskit Nature it should be an implementation of the GroundStateSolver interface. This, by design, needs to (and can be) independent from any classical backend.

Due to current priorities, I do not currently know of anyone who will be able or has the time to help you in case you get stuck. You can try taking problems either to stackoverflow (as mentioned in my previous comment) or reply here, but I cannot guarantee that you will get help. Just so that you are aware of this, in case I go radio silent, that will be due to time and priority constraints.

qiskit-community / qiskit-nature

Help with refactoring / migrating OO-VQE #1219