Problem with active space feature in "generate_molecular_hamiltonian"

[hrgrimsl]

I am trying to get frozen cores working with an OpenFermion-PySCF workflow, and am having issues. As an example, here I am freezing two electrons (i.e., one core orbital, i.e. 2 qubits.) When I set N_e = 6 and N_qubits = 12, i.e. the full active space, I get matching results between CASCI in PySCF and diagonalization of the Hamiltonian that OF-PSCF gives me. In this example, I use N_e = 4 and N_qubits = 10, and get the wrong reference and CASCI energies from the OF-PSCF Hamiltonian.

def test_adapt_vqe():
"""Test ADAPT on H6."""
if os.path.exists('test') == False:
os.makedirs('test')

geom = 'H 0 0 0; H 0 0 1; H 0 0 2; H 0 0 3; H 0 0 4; H 0 0 5'

basis = "sto-3g"
reference = "rhf"
N_e = 4
N_qubits = 10
multiplicity = 1

#E_nuc, H_core, g, D, C, hf_energy = get_integrals(geom, basis, reference, chkfile = "scr.chk", read = False, active = (int(N_qubits/2), N_e))

#PySCF SCF and CASCI Calculations
mol = gto.M(atom = geom, basis = basis, verbose = False)
mol.build()
mf = scf.RHF(mol)
mf.conv_tol_grad = 1e-8
mf.max_cycle = 10000
mf.init_guess = 'atom'
print("PySCF HF Energy:")
print(mf.kernel())
mycas = mcscf.CASCI(mf, int(N_qubits/2), N_e)
casci = mycas.kernel(verbose=False)
print("PYSCF CASCI:")
print(casci[0])


H = generate_molecular_hamiltonian(geom, basis, multiplicity, n_active_electrons = N_e, n_active_orbitals = int(N_qubits/2))
H = of.linalg.get_sparse_operator(H).real

refstr = ''
for i in range(N_e, N_qubits):
    refstr += '0'
for i in range(0, N_e):
    refstr += '1'
ref = np.zeros(2**N_qubits)
ref[int(refstr,2)] = 1
ref = scipy.sparse.csc_matrix(ref).T

print("OpenFermion-PySCF HF Energy:")
print((ref.T@H@ref)[0,0])
print("OpenFermion-PySCF CASCI Energy:")
w, v = scipy.sparse.linalg.eigsh(H, which = "SA")
print(w[0])
exit()

Output:

PySCF HF Energy:
-3.1355322139663198
PYSCF CASCI:
-3.1981030445646623
OpenFermion-PySCF HF Energy:
-1.3744190019034384
OpenFermion-PySCF CASCI Energy:
-3.1614199163979224

[ncrubin]

Does the ground state wavefunction have the correct number of electrons in the case that fails?

[hrgrimsl]

No, when I use a number operator on the ground-state wavefunction of the frozen core Hamiltonian, I get 6 electrons instead of 4.

[tsubasa-iino]

@hrgrimsl

There seem to be two problems in your code.

The first is how you specify the geometry.
PySCF accepts both "string" and "list" as the geometry specification. (here)
However, In "MolecularData" class of OpenFermion, the geometry is assumed to be a list.
If you provide a string, MolecularData's "n_electrons" will be set to 0. (here)
The function "generate_molecular_hamiltonian" in OpenFermion-PySCF uses this "n_electrons" to generate occupied_indices and active_indices (here), and calls the OpenFermion function "get_molecular_hamiltonian", resulting in the improper active space.
In the case of your code, active_indices=[-2, -1, 0, 1, 2], occupied_indices=[].
Note that if you did not set the active space in the "generate_molecular_hamiltonian" function, the "get_molecular_hamiltonian" function seems to work correctly. (here)

The second problem is that the "refstr" (Hartree-Fock state basis) is incorrect.
In the case of your code, the "refstr" is "0000001111", but the correct "refstr" is "1111000000".

The correct code and output are shown below.

from pyscf import gto, scf, mcscf
from openfermionpyscf import generate_molecular_hamiltonian
from openfermion.linalg import get_sparse_operator
import numpy as np
import scipy

# geom = "H 0 0 0; H 0 0 1; H 0 0 2; H 0 0 3; H 0 0 4; H 0 0 5"
geom = [["H", (0., 0., 0.)],
        ["H", (0., 0., 1.)],
        ["H", (0., 0., 2.)],
        ["H", (0., 0., 3.)],
        ["H", (0., 0., 4.)],
        ["H", (0., 0., 5.)]]

basis = "sto-3g"
reference = "rhf"
N_e = 4
N_qubits = 10
multiplicity = 1

# PySCF SCF and CASCI Calculations
mol = gto.M(atom=geom, basis=basis, verbose=False)
mol.build()
mf = scf.RHF(mol)
mf.conv_tol_grad = 1e-8
mf.max_cycle = 10000
mf.init_guess = "atom"
print("PySCF HF Energy:")
print(mf.kernel())
mycas = mcscf.CASCI(mf, int(N_qubits/2), N_e)
casci = mycas.kernel(verbose=False)
print("PYSCF CASCI:")
print(casci[0])

H = generate_molecular_hamiltonian(geom, 
                                   basis,
                                   multiplicity, 
                                   n_active_electrons=N_e,
                                   n_active_orbitals=int(N_qubits/2))
H = get_sparse_operator(H).real

refstr = ""
# for i in range(N_e, N_qubits):
#     refstr += "0"
# for i in range(0, N_e):
#     refstr += "1"
for i in range(0, N_e):
    refstr += "1"
for i in range(N_e, N_qubits):
    refstr += "0"
ref = np.zeros(2**N_qubits)
ref[int(refstr, 2)] = 1
ref = scipy.sparse.csc_matrix(ref).T

print("OpenFermion-PySCF HF Energy:")
print((ref.T@H@ref)[0,0])
print("OpenFermion-PySCF CASCI Energy:")
w, v = scipy.sparse.linalg.eigsh(H, which = "SA")
print(w[0])
exit()

Output:

PySCF HF Energy:
-3.135532213966319
PYSCF CASCI:
-3.1981030445646605
OpenFermion-PySCF HF Energy:
-3.1355322139663198
OpenFermion-PySCF CASCI Energy:
-3.1981030446364223

It should be noted that the lowest eigenvalue obtained by Hamiltonian diagonalization does not necessarily correspond to the CASCI energy. This is because Hamiltonian diagonalization is performed by linear combinations of all electron configurations, so the obtained eigenstate do not necessarily preserve the number of electrons, Sz, etc.

[hrgrimsl]

Sorry for the delay, but I just confirmed this works. Thanks for the help!