Problem with active space feature in "generate_molecular_hamiltonian" about openfermion-pyscf HOT 4 OPEN

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

from openfermion-pyscf.

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

from openfermion-pyscf.

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

from openfermion-pyscf.

@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.

from openfermion-pyscf.