Some issues when migrating vqe_h2o.ipynb to use CO2

[buttercutter]

I am trying to migrate vqe_h2o.ipynb to use CO2 carbon dioxide molecule.

1. For print(fh.terms[((0, 1), (0, 0))]) which is the two-body coefficients, it gives -36.01932840424451. When I try to access the one-body coefficients using the following code snippet, nothing is printed out. Could you advise how to decide which orbitals to be removed for CO2 ? The video comes with a notebook code which no longer runs since it is using deprecated qiskit-chemistry package. Besides, when I checked https://qiskit.org/documentation/nature/tutorials/05_problem_transformers.html#The-ActiveSpaceTransformer , code box [11] still need to be manually hand-picking the active orbitals. Correct me if I miss anything though

one_body_terms = fh.terms
for term, coefficient in one_body_terms.items():
    if len(term) == 2:
        # Only consider terms with one creation and one annihilation operator
        i, j = term
        if i[0] == "c" and j[0] == "c":
            print("Term: c^\dagger_{} c_{} Coefficient: {}".format(i[1], j[1], coefficient))

2. Besides, the kernel crashes when running ma = tc.quantum.PauliStringSum2COO_numpy(lsa, wa)

[liwt31]

1. fh is an openfermion object. Please refer to the official document for the data structure. If you wish to print one-body integrals, the following code works

one_body_terms = fh.terms
for term, coefficient in one_body_terms.items():
    if len(term) == 2:
        # Only consider terms with one creation and one annihilation operator
        i, j = term
        print("Term: c^\dagger_{} c_{} Coefficient: {}".format(i[0], j[0], coefficient))

As for the orbital index to be removed, this usually involves point-group symmetry considerations in quantum chemistry and is different from case to case. This is a relevant paper. 2. The error is likely due to the large memory required to construct the Hamiltonian matrix. In general CO2 is a very large system for quantum computing, please consider using a smaller system such as LiH.

[buttercutter]

@liwt31 Why did you remove the line if i[0] == "c" and j[0] == "c": ?

[liwt31]

@liwt31 Why did you remove the line if i[0] == "c" and j[0] == "c": ?

i, j are tuples of integers. No string is contained.

[buttercutter]

@liwt31 I suppose that I cannot use one_body_terms = fh.terms to do point group symmetry for determining orbital removal index ?

Do I need to resort to using Q-Chem or ORCA ?

I noticed that there is an external user library that does point group symmetry, but I am not familiar with it.

[liwt31]

@liwt31 I suppose that I cannot use one_body_terms = fh.terms to do point group symmetry for determining orbital removal index ?

No, I believe more advanced tools are necessary.

Do I need to resort to using Q-Chem or ORCA ?

I noticed that there is an external user library that does point group symmetry, but I am not familiar with it.

To the best of my knowledge currently there's no readily available tool for this, unfortunately. And I don't believe tensorcircuit is the correct place to implement it. IMHO the only viable way is to read the paper I mentioned above along with several important references in the paper, understand the theory, and then implement the reduction by yourself.

[buttercutter]

@liwt31 @gigaBrainIO

Abby Mitchell from qiskit team told me to check the paper reference found in https://qiskit.org/documentation/nature/stubs/qiskit_nature.second_q.transformers.ActiveSpaceTransformer.html

What do you think ?

[liwt31]

This is certainly a viable method. If you do not care so much about accuracy, you can directly remove orbitals with the highest and lowest energy and retain the ones that are closest to HOMO/LUMO orbitals

[buttercutter]

The relevant code snippet within ActiveSpaceTransformer is as follows, however I am not sure how it works.

    def _determine_active_space(
        self, grouped_property: GroupedElectronicProperty
    ) -> Tuple[List[int], List[int]]:
        """Determines the active and inactive orbital indices.
        Args:
            grouped_property: the `GroupedElectronicProperty` to be transformed.
        Returns:
            The list of active and inactive orbital indices.
        """
        particle_number = grouped_property.get_property(ParticleNumber)
        if isinstance(self._num_electrons, tuple):
            num_alpha, num_beta = self._num_electrons
        elif isinstance(self._num_electrons, (int, np.integer)):
            num_alpha = num_beta = self._num_electrons // 2

        # compute number of inactive electrons
        nelec_total = particle_number._num_alpha + particle_number._num_beta
        nelec_inactive = nelec_total - num_alpha - num_beta

        self._validate_num_electrons(nelec_inactive)
        self._validate_num_orbitals(nelec_inactive, particle_number)

        # determine active and inactive orbital indices
        if self._active_orbitals is None:
            norbs_inactive = nelec_inactive // 2
            inactive_orbs_idxs = list(range(norbs_inactive))
            active_orbs_idxs = list(
                range(norbs_inactive, norbs_inactive + self._num_molecular_orbitals)
            )
        else:
            active_orbs_idxs = self._active_orbitals
            inactive_orbs_idxs = [
                o
                for o in range(nelec_total // 2)
                if o not in self._active_orbitals and self._mo_occ_total[o] > 0
            ]

        return (active_orbs_idxs, inactive_orbs_idxs)

[buttercutter]

directly remove orbitals with the highest and lowest energy and retain the ones that are closest to HOMO/LUMO orbitals

You mean Gouterman orbitals ?

PS: I need some theoretical study on LUMO and LUMO+1 Gouterman orbitals, as well as HOMO-LUMO gaps

[buttercutter]

@liwt31 Thanks for your suggestion, I have the LUMO and HOMO molecular orbitals now for Fe2+ metallic ion for carbon capture MOFs.