BUG: constraints in interface

[kzsigmond]

Describe the problem:

The constraints are not passed correctly between fanpy and PyCI. Pyci requires constraints as {"const 1": [f, dfdx], ...}, while fanpy requires constraints to be a list of BaseSchrodinger instances. This results in a bug if we want to use constraints in FanPT.

Code that reproduces the bug:

Step 1: update fanpt init

In fanpt init change the pyci class assignment to:

        # Convert ProjectedSchrodinger to PyCI objective interface
        constraints = {}
        constraints[f"<\\psi_{{{fanpy_objective.refwfn}}}|\\Psi> - v_{{{fanpy_objective.refwfn}}}"] = fanpy_objective.constraints[0]
        constraints["Energy"] = fanpy_objective.constraints[1]
        fanci_interface = fanpy.interface.pyci.PYCI(fanpy_objective, energy_nuc, legacy_fanci=legacy_fanci, constraints=constraints)

Note: this is technically not the correct datatype for constraints, but the docstring for PyCI class does not mention any data types and doesn't perform any checks.

Step 2: Run the following script to get the traceback error:

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

import sys
from fanpy.wfn.ci.cisd import CISD
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 = 2.5...')
mol = gto.M(atom = [['H', ['0.0', '0.0', '0.0']], ['H', ['0.0', '0.0', 2.5]], ['H', ['0.0', '0.0', 5.0]], ['H', ['0.0', '0.0', 7.5]], ['H', ['0.0', '0.0', 10.0]], ['H', ['0.0', '0.0', 12.5]], ['H', ['0.0', '0.0', 15.0]], ['H', ['0.0', '0.0', 17.5]]],
            unit = 'angstrom',
            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 = 2.5...')

# Initialize wavefunction
wfn = CISD(pyscf_interface.nelec, pyscf_interface.nspinorb)
wfn.assign_params(wfn.params + 0.0001 * 2 * (np.random.rand(*wfn.params.shape) - 0.5))
print('Wavefunction: CISD')
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, 7, 8]
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
from fanpy.eqn.constraints.norm import NormConstraint
from fanpy.eqn.constraints.energy import EnergyConstraint

norm_const = NormConstraint(wfn)
hf_electronic_energy = pyscf_interface.energy_elec
if hf_electronic_energy < 0:
    ref_e = 2*hf_electronic_energy
else:
    ref_e = - hf_electronic_energy
energy_const = EnergyConstraint(wfn, ham, ref_energy=ref_e)

objective = ProjectedSchrodinger(wfn, ham, energy_type="compute", pspace=pspace, constraints=[norm_const, energy_const])

from fanpy.fanpt import FANPT

fanpt = FANPT(objective, pyscf_interface.energy_nuc, 
        energy_active=False, final_order=1,
              steps = 10, legacy_fanci=False)

results = fanpt.optimize(
    guess_params = wfn.params,
    guess_energy = pyscf_interface.energy_elec,
    mode='lstsq', 
    use_jac=True, 
    xtol=1.0e-8,
    ftol=1.0e-8, 
    gtol=1.0e-8,
    max_nfev=wfn.nparams,
    verbose=2)

# 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))

Traceback error:

Traceback (most recent call last):
  File "/home/quintana/Documents/code/Fanpy/tests/test_calc.py", line 92, in <module>
    results = fanpt.optimize(
  File "/home/quintana/Documents/code/Fanpy/fanpy/fanpt/fanpt.py", line 274, in optimize
    self.fanci_interface.update_objective(self.ham0)
  File "/home/quintana/Documents/code/Fanpy/fanpy/interface/pyci.py", line 331, in update_objective
    new_fanpy_objective = fanpy_objective_class(
  File "/home/quintana/Documents/code/Fanpy/fanpy/eqn/projected.py", line 267, in __init__
    self.assign_constraints(constraints)
  File "/home/quintana/Documents/code/Fanpy/fanpy/eqn/projected.py", line 398, in assign_constraints
    raise TypeError(
TypeError: Constraints must be given as a BaseSchrodinger instance or list/tuple of BaseSchrodinger instances.

Package versions:

python: 3.10.16
fanpy: N/A
numpy: 1.26.4

Additonal comments

No response