BUG: Discrepancy between least squares results #122

# PySCF Calculation
import time
from pyscf import gto, scf
import numpy as np

from fanpy.wfn.cc.ap1rog_generalized_NEW import AP1roGSDGeneralized
from fanpy.ham.restricted_chemical import RestrictedMolecularHamiltonian
from fanpy.tools.sd_list import sd_list
import fanpy.interface as interface


# Record the start time
start_time = time.time()

energy = {}

print('# PySCF calculation for alpha = 20.0...')
mol = gto.M(atom = [['Be', ['0.0', '0.0', '0.0']], ['H', ['0.0', '-0.7', '20.0']], ['H', ['0.0', '0.7', '20.0']]],
            unit = 'B',
            basis = 'sto-6g')

myhf = scf.HF(mol)
myhf.kernel()
new_guess, _, stable, _ = myhf.stability(return_status=True)

if not stable:
    print("HF stability failed, trying Newton method...")
    myhf = myhf.newton().run(new_guess, myhf.mo_occ)
    _, _, stable, _ =  myhf.stability(return_status=True)

if not stable:
    raise RuntimeError("HF stability failed")
# PySCF interface
pyscf_interface = interface.pyscf.PYSCF(myhf)

print("\n>> >> Fanpy")

   
# Number of electrons
print("Number of Electrons: {}".format(pyscf_interface.nelec))

# Number of spin orbitals
print("Number of Spin Orbitals: {}".format(pyscf_interface.nspinorb))

# Nuclear-nuclear repulsion
print("Nuclear-nuclear repulsion: {}".format(pyscf_interface.energy_nuc))

# Fanpy calculation
print('# Fanpy calculation for alpha = 1.0...')

# Initialize wavefunction
wfn = AP1roGSDGeneralized(pyscf_interface.nelec,
                          pyscf_interface.nspinorb,
                          params=None,
                          memory=None,
                          ranks=None,
                          indices=None,
                          refwfn=None,
                          s_type='sen-o',
                          exop_combinations=None)
wfn.assign_params(wfn.params + 0.0001 * 2 * (np.random.rand(*wfn.params.shape) - 0.5))
print('Wavefunction: AP1roGSDGeneralized (sen-o)')
print(f"Number of parameters: {wfn.nparams}")

# Initialize Hamiltonian
ham = RestrictedMolecularHamiltonian(pyscf_interface.one_int, pyscf_interface.two_int, update_prev_params=True)
print('Hamiltonian: RestrictedMolecularHamiltonian')

# Projection space
exc_orders = [1, 2, 3, 4, 5, 6]
pspace = sd_list(pyscf_interface.nelec, pyscf_interface.nspinorb, num_limit=None, exc_orders=exc_orders, spin=0)
nproj = len(pspace)
print("Projection space (orders of excitations): ", exc_orders)
print("Number of projections (CISD): {:}".format(nproj))

print("\n# New interface!")
# Initialize objective
from fanpy.eqn.projected import ProjectedSchrodinger

objective = ProjectedSchrodinger(wfn, ham, energy_type="compute", pspace=pspace)

pyci_interface = interface.pyci.PYCI(objective, pyscf_interface.energy_nuc, legacy_fanci=False)
results = pyci_interface.objective.optimize(
    np.hstack([wfn.params, pyscf_interface.energy_elec]), mode="lstsq", use_jac=False, xtol=1e-3, gtol=1e-3, ftol=1e-3
)
# Record the end time
end_time = time.time()

# Calculate the total time taken
execution_time = end_time - start_time

print(f'Total execution time: {execution_time} seconds')
print("### Final energy: {:}".format(results["energy"] + pyscf_interface.energy_nuc))

print("results[success]: ", results["success"])
print("results[message]: ", results["message"])

results = pyci_interface.objective.optimize(
    np.hstack([wfn.params, pyscf_interface.energy_elec]), mode="lstsq", use_jac=True
)