ffsim_integration

ffsim Integration Module for QSCI

This module provides integration between ffsim's UCJ/LUCJ ansatz and the QSCI framework. It requires ffsim to be installed: pip install quri-qsci[ffsim]

Key Components: - molecular_systems: Molecular system creation utilities - integration: ffsim UCJ/LUCJ ansatz integration - state_conversion: State format conversion between ffsim and QURI Parts - qsci_interface: High-level QSCI-ffsim integration interface

Modules:

Name	Description
integration	ffsim integration module for UCJ and LUCJ ansatz generation.
molecular_systems	Molecular system setup utilities for ffsim/UCJ ansatz integration.
qsci_interface	Main interface for ffsim/UCJ ansatz integration with QSCI framework.
state_conversion	State conversion bridge between ffsim and QURI Parts.

Classes:

Name	Description
MolecularSystem	Container for molecular system data.
UCJResult	Container for UCJ ansatz optimization results.
ConversionMetrics	Metrics for state conversion quality.
LUCJQSCIResult	Result from LUCJ/UCJ + QSCI workflow.

Functions:

Name	Description
create_h2_molecule	Create H2 molecule system for testing and validation.
create_n2_molecule	Create N2 molecule system for benchmark studies.
get_reference_energies	Get reference energies for common molecules and basis sets.
validate_molecular_system	Validate that a molecular system is properly constructed.
create_ucj_ansatz	Create and optimize UCJ (Unitary Coupled Cluster Jastrow) ansatz.
create_lucj_ansatz	Create and optimize LUCJ (Linear UCJ) ansatz.
ffsim_to_quri_state	Convert ffsim state vector to QURI Parts CircuitQuantumState.
ucj_result_to_quri_state	Convert UCJ/LUCJ result to QURI Parts state with validation.
run_lucj_qsci	Run complete LUCJ/UCJ + QSCI workflow.
print_result_summary	Print a comprehensive summary of LUCJ/UCJ + QSCI results.
run_convergence_study	Run a convergence study with increasing subspace sizes.
benchmark_against_reference	Benchmark LUCJ/UCJ + QSCI results against reference energy.

Attributes:

Name	Type	Description
FFSIM_AVAILABLE

FFSIM_AVAILABLE module-attribute

FFSIM_AVAILABLE = True

MolecularSystem dataclass

MolecularSystem(
    name,
    geometry,
    basis,
    charge,
    spin,
    bond_length,
    mole,
    scf_result,
    active_space,
    mo_integrals,
    hartree_fock_energy,
    fci_energy,
    quri_hamiltonian,
    ffsim_mol_data,
    jw_mapping,
)

Container for molecular system data.

Attributes:

Name	Type	Description
name	str
geometry	str
basis	str
charge	int
spin	int
bond_length	float
mole	Any
scf_result	Any
active_space	Any
mo_integrals	Any
hartree_fock_energy	float
fci_energy	float
quri_hamiltonian	Operator
ffsim_mol_data	Any
jw_mapping	Any

name instance-attribute

name

geometry instance-attribute

geometry

basis instance-attribute

basis

charge instance-attribute

charge

spin instance-attribute

spin

bond_length instance-attribute

bond_length

mole instance-attribute

mole

scf_result instance-attribute

scf_result

active_space instance-attribute

active_space

mo_integrals instance-attribute

mo_integrals

hartree_fock_energy instance-attribute

hartree_fock_energy

fci_energy instance-attribute

fci_energy

quri_hamiltonian instance-attribute

quri_hamiltonian

ffsim_mol_data instance-attribute

ffsim_mol_data

jw_mapping instance-attribute

jw_mapping

UCJResult dataclass

UCJResult(
    ansatz_type,
    optimized_parameters,
    final_energy,
    n_reps,
    state_vector,
    optimization_success,
    n_iterations,
)

Container for UCJ ansatz optimization results.

Attributes:

Name	Type	Description
ansatz_type	str
optimized_parameters	ndarray
final_energy	float
n_reps	int
state_vector	ndarray
optimization_success	bool
n_iterations	int

ansatz_type instance-attribute

ansatz_type

optimized_parameters instance-attribute

optimized_parameters

final_energy instance-attribute

final_energy

n_reps instance-attribute

n_reps

state_vector instance-attribute

state_vector

optimization_success instance-attribute

optimization_success

n_iterations instance-attribute

n_iterations

ConversionMetrics dataclass

ConversionMetrics(
    fidelity,
    probability_overlap,
    state_vector_norm,
    max_probability_diff,
    conversion_method,
)

Metrics for state conversion quality.

Attributes:

Name	Type	Description
fidelity	float
probability_overlap	float
state_vector_norm	float
max_probability_diff	float
conversion_method	str

fidelity instance-attribute

fidelity

probability_overlap instance-attribute

probability_overlap

state_vector_norm instance-attribute

state_vector_norm

max_probability_diff instance-attribute

max_probability_diff

conversion_method instance-attribute

conversion_method

LUCJQSCIResult dataclass

LUCJQSCIResult(
    molecule_name,
    ansatz_type,
    ansatz_energy,
    ansatz_optimization_success,
    qsci_results,
    conversion_metrics,
    hartree_fock_energy,
    fci_energy,
    target_energy,
    total_time,
    ansatz_time,
    conversion_time,
    qsci_time,
)

Result from LUCJ/UCJ + QSCI workflow.

Attributes:

Name	Type	Description
molecule_name	str
ansatz_type	str
ansatz_energy	float
ansatz_optimization_success	bool
qsci_results	Dict[int, QSCIResult]
conversion_metrics	ConversionMetrics
hartree_fock_energy	float
fci_energy	float
target_energy	Optional[float]
total_time	float
ansatz_time	float
conversion_time	float
qsci_time	float
best_qsci_energy	float	Best (lowest) energy from QSCI calculations.
best_subspace_size	int	Subspace size that gave the best energy.
energy_improvement_vs_ansatz	float	Energy improvement of best QSCI result vs ansatz.
energy_error_vs_fci	float	Energy error of best QSCI result vs FCI reference.

molecule_name instance-attribute

molecule_name

ansatz_type instance-attribute

ansatz_type

ansatz_energy instance-attribute

ansatz_energy

ansatz_optimization_success instance-attribute

ansatz_optimization_success

qsci_results instance-attribute

qsci_results

conversion_metrics instance-attribute

conversion_metrics

hartree_fock_energy instance-attribute

hartree_fock_energy

fci_energy instance-attribute

fci_energy

target_energy instance-attribute

target_energy

total_time instance-attribute

total_time

ansatz_time instance-attribute

ansatz_time

conversion_time instance-attribute

conversion_time

qsci_time instance-attribute

qsci_time

best_qsci_energy property

best_qsci_energy

Best (lowest) energy from QSCI calculations.

best_subspace_size property

best_subspace_size

Subspace size that gave the best energy.

energy_improvement_vs_ansatz property

energy_improvement_vs_ansatz

Energy improvement of best QSCI result vs ansatz.

energy_error_vs_fci property

energy_error_vs_fci

Energy error of best QSCI result vs FCI reference.

create_h2_molecule

create_h2_molecule(
    basis="6-31g", bond_length=0.74, charge=0, spin=0
)

Create H2 molecule system for testing and validation.

Parameters:

Name	Type	Description	Default
basis	str	Basis set (default: 6-31g)	'6-31g'
bond_length	float	H-H bond length in Angstrom (default: 0.74)	0.74
charge	int	Molecular charge (default: 0)	0
spin	int	Spin multiplicity (default: 0 for singlet)	0

Returns:

Type	Description
MolecularSystem	MolecularSystem object with all necessary data

Raises:

Type	Description
ImportError	If ffsim is not installed

Source code in ffsim_integration/molecular_systems.py

def create_h2_molecule(
    basis: str = "6-31g", 
    bond_length: float = 0.74,
    charge: int = 0,
    spin: int = 0
) -> MolecularSystem:
    """Create H2 molecule system for testing and validation.

    Args:
        basis: Basis set (default: 6-31g)
        bond_length: H-H bond length in Angstrom (default: 0.74)
        charge: Molecular charge (default: 0)
        spin: Spin multiplicity (default: 0 for singlet)

    Returns:
        MolecularSystem object with all necessary data

    Raises:
        ImportError: If ffsim is not installed
    """
    _require_ffsim()
    print(f"Creating H2 molecule with {basis} basis, bond length {bond_length} Å")

    # Create geometry string
    geometry = f"H 0 0 0; H 0 0 {bond_length}"

    try:
        # Create PySCF molecule
        geometry_list = [
            ["H", (0.0, 0.0, 0.0)],
            ["H", (0.0, 0.0, bond_length)]
        ]

        mole = gto.M(
            atom=geometry_list,
            basis=basis,
            charge=charge,
            spin=spin,
            verbose=0
        )

        # Perform Hartree-Fock calculation
        mf = scf.RHF(mole)
        mf.verbose = 0
        mf.run()

        if not mf.converged:
            raise RuntimeError("Hartree-Fock calculation did not converge")

        # Generate molecular orbital integrals for QURI Parts
        active_space, mo_eint_set = get_spin_mo_integrals_from_mole(mole, mf.mo_coeff)

        # Create QURI Parts Hamiltonian
        quri_hamiltonian, jw_mapping = get_qubit_mapped_hamiltonian(
            active_space,
            mo_eint_set,
            sz=None,
            fermion_qubit_mapping=jordan_wigner
        )

        # Calculate FCI reference energy
        fci_energy = _calculate_fci_energy(mole, mf)

        # Create ffsim molecular data
        ffsim_mol_data = ffsim.MolecularData.from_scf(mf)

        print(f"✓ H2 system created successfully")
        print(f"  - Active space: {active_space.n_active_ele} electrons, {active_space.n_active_orb} orbitals")
        print(f"  - Qubit count: {2 * active_space.n_active_orb}")
        print(f"  - HF energy: {mf.e_tot:.6f} Ha")
        print(f"  - FCI energy: {fci_energy:.6f} Ha")
        print(f"  - Correlation energy: {mf.e_tot - fci_energy:.6f} Ha")

        return MolecularSystem(
            name="H2",
            geometry=geometry,
            basis=basis,
            charge=charge,
            spin=spin,
            bond_length=bond_length,
            mole=mole,
            scf_result=mf,
            active_space=active_space,
            mo_integrals=mo_eint_set,
            hartree_fock_energy=mf.e_tot,
            fci_energy=fci_energy,
            quri_hamiltonian=quri_hamiltonian,
            ffsim_mol_data=ffsim_mol_data,
            jw_mapping=jw_mapping
        )

    except Exception as e:
        print(f"Error creating H2 molecule: {e}")
        raise

create_n2_molecule

create_n2_molecule(
    basis="6-31g",
    bond_length=1.0,
    charge=0,
    spin=0,
    active_space=None,
)

Create N2 molecule system for benchmark studies.

Parameters:

Name	Type	Description	Default
basis	str	Basis set (default: 6-31g)	'6-31g'
bond_length	float	N-N bond length in Angstrom (default: 1.0)	1.0
charge	int	Molecular charge (default: 0)	0
spin	int	Spin multiplicity (default: 0 for singlet)	0
active_space	Optional[Tuple[int, int]]	(n_electrons, n_orbitals) for active space (None for full space)	None

Returns:

Type	Description
MolecularSystem	MolecularSystem object with all necessary data

Raises:

Type	Description
ImportError	If ffsim is not installed

Source code in ffsim_integration/molecular_systems.py

def create_n2_molecule(
    basis: str = "6-31g",
    bond_length: float = 1.0,
    charge: int = 0,
    spin: int = 0,
    active_space: Optional[Tuple[int, int]] = None
) -> MolecularSystem:
    """Create N2 molecule system for benchmark studies.

    Args:
        basis: Basis set (default: 6-31g) 
        bond_length: N-N bond length in Angstrom (default: 1.0)
        charge: Molecular charge (default: 0)
        spin: Spin multiplicity (default: 0 for singlet)
        active_space: (n_electrons, n_orbitals) for active space (None for full space)

    Returns:
        MolecularSystem object with all necessary data

    Raises:
        ImportError: If ffsim is not installed
    """
    _require_ffsim()
    print(f"Creating N2 molecule with {basis} basis, bond length {bond_length} Å")

    # Create geometry string
    geometry = f"N 0 0 0; N 0 0 {bond_length}"

    try:
        # Create PySCF molecule
        geometry_list = [
            ["N", (0.0, 0.0, 0.0)],
            ["N", (0.0, 0.0, bond_length)]
        ]

        mole = gto.M(
            atom=geometry_list,
            basis=basis,
            charge=charge,
            spin=spin,
            verbose=0
        )

        # Perform Hartree-Fock calculation
        mf = scf.RHF(mole)
        mf.verbose = 0
        mf.run()

        if not mf.converged:
            raise RuntimeError("Hartree-Fock calculation did not converge")

        # Generate molecular orbital integrals for QURI Parts
        active_space, mo_eint_set = get_spin_mo_integrals_from_mole(mole, mf.mo_coeff)

        # Create QURI Parts Hamiltonian
        quri_hamiltonian, jw_mapping = get_qubit_mapped_hamiltonian(
            active_space,
            mo_eint_set,
            sz=None,
            fermion_qubit_mapping=jordan_wigner
        )

        # Calculate FCI reference energy
        fci_energy = _calculate_fci_energy(mole, mf)

        # Create ffsim molecular data
        ffsim_mol_data = ffsim.MolecularData.from_scf(mf)

        print(f"✓ N2 system created successfully")
        print(f"  - Active space: {active_space.n_active_ele} electrons, {active_space.n_active_orb} orbitals")
        print(f"  - Qubit count: {2 * active_space.n_active_orb}")
        print(f"  - HF energy: {mf.e_tot:.6f} Ha")
        print(f"  - FCI energy: {fci_energy:.6f} Ha")
        print(f"  - Correlation energy: {mf.e_tot - fci_energy:.6f} Ha")
        print(f"  - ffsim molecule details:")
        print(f"    - norb: {ffsim_mol_data.norb}")
        print(f"    - nelec: {ffsim_mol_data.nelec}")
        print(f"    - Expected ffsim Fock space: C({ffsim_mol_data.norb},{ffsim_mol_data.nelec[0]}) * C({ffsim_mol_data.norb},{ffsim_mol_data.nelec[1]})")

        from scipy.special import comb
        expected_ffsim_dim = int(comb(ffsim_mol_data.norb, ffsim_mol_data.nelec[0]) * comb(ffsim_mol_data.norb, ffsim_mol_data.nelec[1]))
        print(f"    - Expected ffsim dimension: {expected_ffsim_dim}")

        return MolecularSystem(
            name="N2",
            geometry=geometry,
            basis=basis,
            charge=charge,
            spin=spin,
            bond_length=bond_length,
            mole=mole,
            scf_result=mf,
            active_space=active_space,
            mo_integrals=mo_eint_set,
            hartree_fock_energy=mf.e_tot,
            fci_energy=fci_energy,
            quri_hamiltonian=quri_hamiltonian,
            ffsim_mol_data=ffsim_mol_data,
            jw_mapping=jw_mapping
        )

    except Exception as e:
        print(f"Error creating N2 molecule: {e}")
        raise

get_reference_energies

get_reference_energies(molecule, basis='6-31g')

Get reference energies for common molecules and basis sets.

Parameters:

Name	Type	Description	Default
molecule	str	Molecule name ("H2" or "N2")	required
basis	str	Basis set name	'6-31g'

Returns:

Type	Description
Dict[str, float]	Dictionary with reference energies in Hartree

Source code in ffsim_integration/molecular_systems.py

def get_reference_energies(molecule: str, basis: str = "6-31g") -> Dict[str, float]:
    """Get reference energies for common molecules and basis sets.

    Args:
        molecule: Molecule name ("H2" or "N2")
        basis: Basis set name

    Returns:
        Dictionary with reference energies in Hartree
    """
    # Reference energies from literature/high-level calculations
    references = {
        "H2": {
            "6-31g": {
                "hartree_fock": -1.123,  # Approximate HF energy
                "fci": -1.135,           # Approximate FCI energy
                "experimental": -1.139   # Approximate experimental
            },
            "sto-3g": {
                "hartree_fock": -1.117,
                "fci": -1.127,
                "experimental": -1.139
            }
        },
        "N2": {
            "6-31g": {
                "hartree_fock": -108.9,  # Approximate HF energy
                "fci": -109.1,           # Approximate FCI energy  
                "experimental": -109.3   # Target benchmark energy
            },
            "sto-3g": {
                "hartree_fock": -107.5,
                "fci": -107.7,
                "experimental": -109.3
            }
        }
    }

    molecule_upper = molecule.upper()
    basis_lower = basis.lower()

    if molecule_upper in references and basis_lower in references[molecule_upper]:
        return references[molecule_upper][basis_lower]
    else:
        print(f"Warning: No reference energies available for {molecule} with {basis} basis")
        return {"hartree_fock": 0.0, "fci": 0.0, "experimental": 0.0}

validate_molecular_system

validate_molecular_system(system)

Validate that a molecular system is properly constructed.

Parameters:

Name	Type	Description	Default
system	MolecularSystem	MolecularSystem to validate	required

Returns:

Type	Description
bool	True if system is valid, False otherwise

Source code in ffsim_integration/molecular_systems.py

def validate_molecular_system(system: MolecularSystem) -> bool:
    """Validate that a molecular system is properly constructed.

    Args:
        system: MolecularSystem to validate

    Returns:
        True if system is valid, False otherwise
    """
    try:
        # Check basic attributes
        assert system.mole is not None
        assert system.scf_result is not None
        assert system.scf_result.converged

        # Check active space
        assert system.active_space is not None
        assert system.active_space.n_active_ele > 0
        assert system.active_space.n_active_orb > 0

        # Check Hamiltonians
        assert system.quri_hamiltonian is not None
        assert system.ffsim_mol_data is not None

        # Check energy ordering (FCI should be lower than HF)
        assert system.fci_energy <= system.hartree_fock_energy

        print(f"✓ {system.name} molecular system validation passed")
        return True

    except AssertionError as e:
        print(f"✗ {system.name} molecular system validation failed: {e}")
        return False
    except Exception as e:
        print(f"✗ {system.name} molecular system validation error: {e}")
        return False

create_ucj_ansatz

create_ucj_ansatz(
    mol_system,
    n_reps=1,
    optimization_method="BFGS",
    max_iterations=100,
)

Create and optimize UCJ (Unitary Coupled Cluster Jastrow) ansatz.

Parameters:

Name	Type	Description	Default
mol_system	MolecularSystem	MolecularSystem containing molecular data	required
n_reps	int	Number of repetitions in the ansatz circuit	1
optimization_method	str	Optimization method for parameter optimization	'BFGS'
max_iterations	int	Maximum optimization iterations	100

Returns:

Type	Description
UCJResult	UCJResult with optimized parameters and state vector

Raises:

Type	Description
ImportError	If ffsim is not installed

Source code in ffsim_integration/integration.py

def create_ucj_ansatz(
    mol_system: MolecularSystem,
    n_reps: int = 1,
    optimization_method: str = "BFGS",
    max_iterations: int = 100
) -> UCJResult:
    """Create and optimize UCJ (Unitary Coupled Cluster Jastrow) ansatz.

    Args:
        mol_system: MolecularSystem containing molecular data
        n_reps: Number of repetitions in the ansatz circuit
        optimization_method: Optimization method for parameter optimization
        max_iterations: Maximum optimization iterations

    Returns:
        UCJResult with optimized parameters and state vector

    Raises:
        ImportError: If ffsim is not installed
    """
    _require_ffsim()
    print(f"Creating UCJ ansatz with {n_reps} repetitions...")

    # Extract molecular data
    mol_data = mol_system.ffsim_mol_data
    norb = mol_data.norb
    nelec = mol_data.nelec

    print(f"  - Norb: {norb}, Nelec: {nelec}")

    # Define interaction pairs for UCJ ansatz (use all possible interactions)
    # For UCJ, we can use None to allow all interactions or define comprehensive pairs
    interaction_pairs = None  # This allows all possible interactions for UCJ

    print(f"  - Interaction pairs: All possible (UCJ)")

    # Create reference state (Hartree-Fock)
    reference_state = ffsim.hartree_fock_state(norb, nelec)

    # Determine correct parameter count by directly extracting from error message
    def get_param_count(norb, n_reps, interaction_pairs=None):
        # Try with a small number first to get the error message
        try:
            test_params = np.zeros(1)
            ffsim.UCJOpSpinBalanced.from_parameters(
                test_params, norb=norb, n_reps=n_reps, interaction_pairs=interaction_pairs
            )
            return 1  # If it works with 1, return 1
        except Exception as e:
            # Extract expected parameter count from error message
            if "Expected" in str(e) and "but got" in str(e):
                try:
                    expected_str = str(e).split("Expected ")[1].split(" but got")[0]
                    expected_count = int(expected_str)
                    print(f"  - ffsim expects {expected_count} parameters for norb={norb}, n_reps={n_reps}")
                    return expected_count
                except:
                    pass

        # Fallback: use empirical formula or brute force with smaller range
        return 210 if interaction_pairs is None else 119

    n_params = get_param_count(norb, n_reps, interaction_pairs)
    initial_params = np.random.random(n_params) * 0.1

    print(f"  - Number of parameters: {n_params}")

    def objective_function(params):
        """Objective function for optimization."""
        try:
            # Create ansatz operator using from_parameters method
            ansatz_op = ffsim.UCJOpSpinBalanced.from_parameters(
                params,
                norb=norb,
                n_reps=n_reps,
                interaction_pairs=interaction_pairs
            )

            # Apply ansatz to reference state
            state = ffsim.apply_unitary(reference_state, ansatz_op, norb=norb, nelec=nelec)

            # Calculate energy expectation value using linear operator
            hamiltonian_linop = ffsim.linear_operator(mol_data.hamiltonian, norb=norb, nelec=nelec)
            energy = np.real(np.conj(state) @ hamiltonian_linop @ state)

            return energy

        except Exception as e:
            print(f"Warning: Energy evaluation failed: {e}")
            return 1e6  # Return high energy on failure

    # Optimize parameters
    print("  - Starting parameter optimization...")

    try:
        opt_result = scipy.optimize.minimize(
            objective_function,
            initial_params,
            method=optimization_method,
            options={'maxiter': max_iterations}
        )

        optimization_success = opt_result.success
        optimized_params = opt_result.x
        final_energy = opt_result.fun
        n_iterations = opt_result.nit

        print(f"  - Optimization {'succeeded' if optimization_success else 'failed'}")
        print(f"  - Final energy: {final_energy:.6f} Ha")
        print(f"  - Iterations: {n_iterations}")

        # Generate final optimized state
        ansatz_op = ffsim.UCJOpSpinBalanced.from_parameters(
            optimized_params,
            norb=norb,
            n_reps=n_reps,
            interaction_pairs=interaction_pairs
        )
        final_state = ffsim.apply_unitary(reference_state, ansatz_op, norb=norb, nelec=nelec)

        return UCJResult(
            ansatz_type="UCJ",
            optimized_parameters=optimized_params,
            final_energy=final_energy,
            n_reps=n_reps,
            state_vector=final_state,
            optimization_success=optimization_success,
            n_iterations=n_iterations
        )

    except Exception as e:
        print(f"Error in UCJ optimization: {e}")
        # Return fallback result with HF state
        return UCJResult(
            ansatz_type="UCJ",
            optimized_parameters=initial_params,
            final_energy=mol_system.hartree_fock_energy,
            n_reps=n_reps,
            state_vector=reference_state,
            optimization_success=False,
            n_iterations=0
        )

create_lucj_ansatz

create_lucj_ansatz(
    mol_system,
    n_reps=1,
    optimization_method="BFGS",
    max_iterations=100,
)

Create and optimize LUCJ (Linear UCJ) ansatz.

Parameters:

Name	Type	Description	Default
mol_system	MolecularSystem	MolecularSystem containing molecular data	required
n_reps	int	Number of repetitions in the ansatz circuit	1
optimization_method	str	Optimization method for parameter optimization	'BFGS'
max_iterations	int	Maximum optimization iterations	100

Returns:

Type	Description
UCJResult	UCJResult with optimized parameters and state vector

Raises:

Type	Description
ImportError	If ffsim is not installed

Source code in ffsim_integration/integration.py

def create_lucj_ansatz(
    mol_system: MolecularSystem,
    n_reps: int = 1,
    optimization_method: str = "BFGS",
    max_iterations: int = 100
) -> UCJResult:
    """Create and optimize LUCJ (Linear UCJ) ansatz.

    Args:
        mol_system: MolecularSystem containing molecular data
        n_reps: Number of repetitions in the ansatz circuit
        optimization_method: Optimization method for parameter optimization
        max_iterations: Maximum optimization iterations

    Returns:
        UCJResult with optimized parameters and state vector

    Raises:
        ImportError: If ffsim is not installed
    """
    _require_ffsim()
    print(f"Creating LUCJ ansatz with {n_reps} repetitions...")

    # Extract molecular data
    mol_data = mol_system.ffsim_mol_data
    norb = mol_data.norb
    nelec = mol_data.nelec

    print(f"  - Norb: {norb}, Nelec: {nelec}")

    # Define interaction pairs for LUCJ ansatz
    # For larger systems, limit interaction pairs to reduce parameter count dramatically
    if norb > 6:  # Large system - use very minimal interactions for manageable parameters
        pairs_aa = [(p, p + 1) for p in range(min(2, norb - 1))]  # Only 2 α-α pairs
        pairs_ab = [(p, p) for p in range(min(2, norb))]          # Only 2 α-β pairs
        print(f"  - Using very reduced interaction pairs for large system (norb={norb})")
        print(f"  - This should give approximately 10-20 parameters (much more manageable)")
    else:  # Small system - use all interactions
        pairs_aa = [(p, p + 1) for p in range(norb - 1)]  # Alpha-alpha pairs
        pairs_ab = [(p, p) for p in range(norb)]          # Alpha-beta pairs

    interaction_pairs = (pairs_aa, pairs_ab)

    print(f"  - Alpha-alpha pairs: {pairs_aa}")
    print(f"  - Alpha-beta pairs: {pairs_ab}")

    # Create reference state (Hartree-Fock)
    reference_state = ffsim.hartree_fock_state(norb, nelec)

    # Determine correct parameter count by directly extracting from error message
    def get_param_count(norb, n_reps, interaction_pairs=None):
        # Try with a small number first to get the error message
        try:
            test_params = np.zeros(1)
            ffsim.UCJOpSpinBalanced.from_parameters(
                test_params, norb=norb, n_reps=n_reps, interaction_pairs=interaction_pairs
            )
            return 1  # If it works with 1, return 1
        except Exception as e:
            # Extract expected parameter count from error message
            if "Expected" in str(e) and "but got" in str(e):
                try:
                    expected_str = str(e).split("Expected ")[1].split(" but got")[0]
                    expected_count = int(expected_str)
                    print(f"  - ffsim expects {expected_count} parameters for norb={norb}, n_reps={n_reps}")
                    return expected_count
                except:
                    pass

        # Fallback: use empirical formula or brute force with smaller range
        return 119 if interaction_pairs is not None else 210

    n_params = get_param_count(norb, n_reps, interaction_pairs)
    initial_params = np.random.random(n_params) * 0.1

    print(f"  - Number of parameters: {n_params}")

    def objective_function(params):
        """Objective function for optimization."""
        try:
            # Create ansatz operator using UCJ with locality constraints (LUCJ)
            ansatz_op = ffsim.UCJOpSpinBalanced.from_parameters(
                params,
                norb=norb,
                n_reps=n_reps,
                interaction_pairs=interaction_pairs  # This makes it LUCJ (local UCJ)
            )

            # Apply ansatz to reference state
            state = ffsim.apply_unitary(reference_state, ansatz_op, norb=norb, nelec=nelec)

            # Calculate energy expectation value using linear operator
            hamiltonian_linop = ffsim.linear_operator(mol_data.hamiltonian, norb=norb, nelec=nelec)
            energy = np.real(np.conj(state) @ hamiltonian_linop @ state)

            return energy

        except Exception as e:
            print(f"Warning: Energy evaluation failed: {e}")
            return 1e6  # Return high energy on failure

    # Optimize parameters
    print("  - Starting parameter optimization...")

    try:
        opt_result = scipy.optimize.minimize(
            objective_function,
            initial_params,
            method=optimization_method,
            options={'maxiter': max_iterations}
        )

        optimization_success = opt_result.success
        optimized_params = opt_result.x
        final_energy = opt_result.fun
        n_iterations = opt_result.nit

        print(f"  - Optimization {'succeeded' if optimization_success else 'failed'}")
        print(f"  - Final energy: {final_energy:.6f} Ha")
        print(f"  - Iterations: {n_iterations}")

        # Generate final optimized state
        ansatz_op = ffsim.UCJOpSpinBalanced.from_parameters(
            optimized_params,
            norb=norb,
            n_reps=n_reps,
            interaction_pairs=interaction_pairs  # This makes it LUCJ (local UCJ)
        )
        final_state = ffsim.apply_unitary(reference_state, ansatz_op, norb=norb, nelec=nelec)

        return UCJResult(
            ansatz_type="LUCJ",
            optimized_parameters=optimized_params,
            final_energy=final_energy,
            n_reps=n_reps,
            state_vector=final_state,
            optimization_success=optimization_success,
            n_iterations=n_iterations
        )

    except Exception as e:
        print(f"Error in LUCJ optimization: {e}")
        # Return fallback result with HF state
        return UCJResult(
            ansatz_type="LUCJ",
            optimized_parameters=initial_params,
            final_energy=mol_system.hartree_fock_energy,
            n_reps=n_reps,
            state_vector=reference_state,
            optimization_success=False,
            n_iterations=0
        )

ffsim_to_quri_state

ffsim_to_quri_state(
    state_vector,
    n_qubits,
    threshold=1e-10,
    method="sampling_circuit",
)

Convert ffsim state vector to QURI Parts CircuitQuantumState.

This is a key function that bridges ffsim output with QSCI input requirements.

Parameters:

Name	Type	Description	Default
state_vector	ndarray	State vector from ffsim (normalized)	required
n_qubits	int	Number of qubits in the system	required
threshold	float	Threshold for significant amplitudes	1e-10
method	str	Conversion method ("sampling_circuit" or "superposition")	'sampling_circuit'

Returns:

Type	Description
CircuitQuantumState	CircuitQuantumState compatible with QSCI algorithms

Source code in ffsim_integration/state_conversion.py

def ffsim_to_quri_state(
    state_vector: np.ndarray, 
    n_qubits: int,
    threshold: float = 1e-10,
    method: str = "sampling_circuit"
) -> CircuitQuantumState:
    """Convert ffsim state vector to QURI Parts CircuitQuantumState.

    This is a key function that bridges ffsim output with QSCI input requirements.

    Args:
        state_vector: State vector from ffsim (normalized)
        n_qubits: Number of qubits in the system
        threshold: Threshold for significant amplitudes
        method: Conversion method ("sampling_circuit" or "superposition")

    Returns:
        CircuitQuantumState compatible with QSCI algorithms
    """
    print(f"Converting ffsim state vector to QURI Parts format...")
    print(f"  - State vector dimension: {len(state_vector)}")
    print(f"  - Number of qubits: {n_qubits}")
    print(f"  - Expected qubit dimension: 2^{n_qubits} = {2**n_qubits}")
    print(f"  - Conversion method: {method}")

    # FIXED: Handle fermionic Fock space representation correctly
    import math
    from scipy.special import comb

    # ffsim uses fermionic Fock space with dimension C(norb, n_alpha) * C(norb, n_beta)
    # NOT the full qubit Hilbert space 2^n_qubits
    expected_qubit_dim = 2**n_qubits
    actual_state_dim = len(state_vector)

    print(f"  - ffsim fermionic state dimension: {actual_state_dim}")
    print(f"  - Expected qubit space dimension: {expected_qubit_dim}")

    if actual_state_dim == expected_qubit_dim:
        print(f"  - Direct qubit representation detected")
    else:
        print(f"  - Fermionic Fock space representation detected")
        print(f"  - Will map fermionic basis states to computational basis states")

    # Normalize state vector
    norm = np.linalg.norm(state_vector)
    if norm == 0:
        raise ValueError("State vector has zero norm")

    normalized_state = state_vector / norm
    print(f"  - State vector norm: {norm:.6f}")

    if method == "sampling_circuit":
        return _create_sampling_circuit(normalized_state, n_qubits, threshold)
    elif method == "superposition":
        return _create_superposition_state(normalized_state, n_qubits, threshold)
    else:
        raise ValueError(f"Unknown conversion method: {method}")

ucj_result_to_quri_state

ucj_result_to_quri_state(
    ucj_result,
    n_qubits,
    conversion_method="sampling_circuit",
)

Convert UCJ/LUCJ result to QURI Parts state with validation.

This is a convenience function that combines conversion and validation.

Parameters:

Name	Type	Description	Default
ucj_result	UCJResult	Result from UCJ/LUCJ optimization	required
n_qubits	int	Number of qubits in the system	required
conversion_method	str	Method for state conversion	'sampling_circuit'

Returns:

Type	Description
Tuple[CircuitQuantumState, ConversionMetrics]	Tuple of (converted_state, conversion_metrics)

Source code in ffsim_integration/state_conversion.py

def ucj_result_to_quri_state(
    ucj_result: UCJResult,
    n_qubits: int,
    conversion_method: str = "sampling_circuit"
) -> Tuple[CircuitQuantumState, ConversionMetrics]:
    """Convert UCJ/LUCJ result to QURI Parts state with validation.

    This is a convenience function that combines conversion and validation.

    Args:
        ucj_result: Result from UCJ/LUCJ optimization
        n_qubits: Number of qubits in the system
        conversion_method: Method for state conversion

    Returns:
        Tuple of (converted_state, conversion_metrics)
    """
    print(f"Converting {ucj_result.ansatz_type} result to QURI Parts state...")

    # Convert state
    quri_state = ffsim_to_quri_state(
        ucj_result.state_vector, 
        n_qubits, 
        method=conversion_method
    )

    # Validate conversion
    metrics = validate_state_conversion(
        ucj_result.state_vector, 
        quri_state, 
        n_shots=5000
    )

    print(f"✓ {ucj_result.ansatz_type} state conversion completed")
    print(f"  - Conversion fidelity: {metrics.fidelity:.4f}")

    return quri_state, metrics

run_lucj_qsci

run_lucj_qsci(
    molecule,
    ansatz_type="LUCJ",
    subspace_sizes=[50, 100, 150],
    basis="sto-3g",
    bond_length=None,
    n_reps=1,
    max_optimization_iterations=50,
    total_shots=5000,
    conversion_method="sampling_circuit",
    active_space=None,
    use_homo_lumo=True,
)

Run complete LUCJ/UCJ + QSCI workflow.

Parameters:

Name	Type	Description	Default
molecule	str	Molecule name ("H2" or "N2")	required
ansatz_type	str	"UCJ" or "LUCJ"	'LUCJ'
subspace_sizes	List[int]	List of subspace dimensions to test	[50, 100, 150]
basis	str	Basis set for quantum chemistry	'sto-3g'
bond_length	Optional[float]	Bond length (None for default)	None
n_reps	int	Number of ansatz repetitions	1
max_optimization_iterations	int	Max iterations for ansatz optimization	50
total_shots	int	Total measurement shots for QSCI	5000
conversion_method	str	State conversion method	'sampling_circuit'
active_space	Optional[Tuple[int, int]]	(n_electrons, n_orbitals) for N2 active space (None for full space)	None
use_homo_lumo	bool	If True, use HOMO-LUMO focused active space for N2 (more efficient)	True

Returns:

Type	Description
LUCJQSCIResult	LUCJQSCIResult with complete workflow results

Source code in ffsim_integration/qsci_interface.py

def run_lucj_qsci(
    molecule: str,
    ansatz_type: str = "LUCJ",
    subspace_sizes: List[int] = [50, 100, 150],
    basis: str = "sto-3g",
    bond_length: Optional[float] = None,
    n_reps: int = 1,
    max_optimization_iterations: int = 50,
    total_shots: int = 5000,
    conversion_method: str = "sampling_circuit",
    active_space: Optional[Tuple[int, int]] = None,
    use_homo_lumo: bool = True
) -> LUCJQSCIResult:
    """Run complete LUCJ/UCJ + QSCI workflow.

    Args:
        molecule: Molecule name ("H2" or "N2")
        ansatz_type: "UCJ" or "LUCJ" 
        subspace_sizes: List of subspace dimensions to test
        basis: Basis set for quantum chemistry
        bond_length: Bond length (None for default)
        n_reps: Number of ansatz repetitions
        max_optimization_iterations: Max iterations for ansatz optimization
        total_shots: Total measurement shots for QSCI
        conversion_method: State conversion method
        active_space: (n_electrons, n_orbitals) for N2 active space (None for full space)
        use_homo_lumo: If True, use HOMO-LUMO focused active space for N2 (more efficient)

    Returns:
        LUCJQSCIResult with complete workflow results
    """
    print(f"=" * 60)
    print(f"Running {ansatz_type} + QSCI workflow for {molecule}")
    print(f"=" * 60)

    start_time = time.time()

    # 1. Create molecular system
    print(f"Step 1: Creating {molecule} molecular system...")
    ansatz_start = time.time()

    if molecule.upper() == "H2":
        bond_length = bond_length or 0.74
        mol_system = create_h2_molecule(basis=basis, bond_length=bond_length)
    elif molecule.upper() == "N2":
        bond_length = bond_length or 1.0
        if use_homo_lumo:
            # Use HOMO-LUMO focused active space for efficient optimization
            mol_system = create_n2_homo_lumo_molecule(
                basis=basis,
                bond_length=bond_length,
                active_space_size=4
            )
        elif active_space is not None:
            n_frozen = active_space[0] if len(active_space) > 0 else 2
            mol_system = create_n2_active_space_molecule(
                n_frozen=n_frozen,
                basis=basis,
                bond_length=bond_length
            )
        else:
            mol_system = create_n2_molecule(basis=basis, bond_length=bond_length)
    else:
        raise ValueError(f"Unsupported molecule: {molecule}")

    # 2. Generate ansatz state
    print(f"\nStep 2: Generating {ansatz_type} ansatz...")

    if ansatz_type.upper() == "UCJ":
        ansatz_result = create_ucj_ansatz(
            mol_system, 
            n_reps=n_reps, 
            max_iterations=max_optimization_iterations
        )
    elif ansatz_type.upper() == "LUCJ":
        ansatz_result = create_lucj_ansatz(
            mol_system, 
            n_reps=n_reps, 
            max_iterations=max_optimization_iterations
        )
    else:
        raise ValueError(f"Unsupported ansatz type: {ansatz_type}")

    ansatz_time = time.time() - ansatz_start

    # 3. Convert to QURI Parts format
    print(f"\nStep 3: Converting {ansatz_type} state to QURI Parts format...")
    conversion_start = time.time()

    # Jordan-Wigner mapping: each spatial orbital needs 2 qubits (alpha and beta spin)
    # This gives us the full computational basis space needed for QSCI
    n_qubits = 2 * mol_system.ffsim_mol_data.norb
    quri_state, conversion_metrics = ucj_result_to_quri_state(
        ansatz_result, n_qubits, conversion_method=conversion_method
    )

    conversion_time = time.time() - conversion_start

    # 4. Run QSCI with different subspace sizes
    print(f"\nStep 4: Running VanillaQSCI with different subspace sizes...")
    qsci_start = time.time()

    qsci_results = {}

    for subspace_size in subspace_sizes:
        print(f"\n  Running QSCI with subspace size {subspace_size}...")

        # Create QSCI algorithm
        qsci_algo = VanillaQSCI(
            hamiltonian=mol_system.quri_hamiltonian,
            sampler=None,  # Will be created automatically
            num_states_pick_out=subspace_size
        )

        # Create sampler
        from quri_parts.qulacs.sampler import create_qulacs_vector_concurrent_sampler
        qsci_algo.sampler = create_qulacs_vector_concurrent_sampler()

        # Run QSCI
        try:
            qsci_result = qsci_algo.run(
                input_states=[quri_state],
                total_shots=total_shots,
                start_time=time.time()
            )
            qsci_results[subspace_size] = qsci_result

            print(f"    ✓ QSCI R={subspace_size}: Energy = {qsci_result.ground_state_energy:.6f} Ha")

        except Exception as e:
            print(f"    ✗ QSCI R={subspace_size} failed: {e}")
            continue

    qsci_time = time.time() - qsci_start
    total_time = time.time() - start_time

    # 5. Create result summary
    print(f"\nStep 5: Summarizing results...")

    # Determine target energy for benchmark
    target_energy = None
    if molecule.upper() == "N2" and basis == "6-31g":
        target_energy = -109.0  # Target benchmark energy

    result = LUCJQSCIResult(
        molecule_name=molecule.upper(),
        ansatz_type=ansatz_type.upper(),
        ansatz_energy=ansatz_result.final_energy,
        ansatz_optimization_success=ansatz_result.optimization_success,
        qsci_results=qsci_results,
        conversion_metrics=conversion_metrics,
        hartree_fock_energy=mol_system.hartree_fock_energy,
        fci_energy=mol_system.fci_energy,
        target_energy=target_energy,
        total_time=total_time,
        ansatz_time=ansatz_time,
        conversion_time=conversion_time,
        qsci_time=qsci_time
    )

    # Print summary
    print_result_summary(result)

    return result

print_result_summary

print_result_summary(result)

Print a comprehensive summary of LUCJ/UCJ + QSCI results.

Source code in ffsim_integration/qsci_interface.py

def print_result_summary(result: LUCJQSCIResult):
    """Print a comprehensive summary of LUCJ/UCJ + QSCI results."""
    print(f"\n" + "=" * 60)
    print(f"WORKFLOW SUMMARY: {result.ansatz_type} + QSCI for {result.molecule_name}")
    print(f"=" * 60)

    print(f"\n📊 ENERGY RESULTS:")
    print(f"  Hartree-Fock energy:     {result.hartree_fock_energy:.6f} Ha")
    print(f"  FCI reference energy:    {result.fci_energy:.6f} Ha")
    print(f"  {result.ansatz_type} ansatz energy:      {result.ansatz_energy:.6f} Ha")
    print(f"  Best QSCI energy:        {result.best_qsci_energy:.6f} Ha")

    if result.target_energy:
        print(f"  Target benchmark energy: {result.target_energy:.6f} Ha")
        error_vs_target = result.best_qsci_energy - result.target_energy
        print(f"  Error vs target:         {error_vs_target:.6f} Ha ({error_vs_target*1000:.1f} mHa)")

    print(f"\n📈 PERFORMANCE METRICS:")
    print(f"  {result.ansatz_type} vs HF improvement:  {result.hartree_fock_energy - result.ansatz_energy:.6f} Ha")
    print(f"  QSCI vs {result.ansatz_type} improvement: {result.energy_improvement_vs_ansatz:.6f} Ha")
    print(f"  Error vs FCI:            {result.energy_error_vs_fci:.6f} Ha")
    print(f"  Best subspace size:      R = {result.best_subspace_size}")

    print(f"\n🔄 CONVERSION QUALITY:")
    print(f"  State conversion fidelity: {result.conversion_metrics.fidelity:.4f}")
    print(f"  Probability overlap:       {result.conversion_metrics.probability_overlap:.4f}")

    print(f"\n⏱️  TIMING BREAKDOWN:")
    print(f"  {result.ansatz_type} optimization time:  {result.ansatz_time:.2f} seconds")
    print(f"  State conversion time:    {result.conversion_time:.2f} seconds")
    print(f"  QSCI calculation time:    {result.qsci_time:.2f} seconds")
    print(f"  Total workflow time:      {result.total_time:.2f} seconds")

    print(f"\n📊 SUBSPACE CONVERGENCE:")
    if result.qsci_results:
        for size in sorted(result.qsci_results.keys()):
            energy = result.qsci_results[size].ground_state_energy
            error = energy - result.fci_energy
            print(f"  R={size:3d}: {energy:.6f} Ha (error: {error:.6f} Ha)")

    print(f"\n" + "=" * 60)

run_convergence_study

run_convergence_study(
    molecule,
    ansatz_type="LUCJ",
    max_subspace=200,
    step_size=25,
    **kwargs
)

Run a convergence study with increasing subspace sizes.

Parameters:

Name	Type	Description	Default
molecule	str	Molecule name	required
ansatz_type	str	Ansatz type	'LUCJ'
max_subspace	int	Maximum subspace size to test	200
step_size	int	Step size for subspace sizes	25
**kwargs		Additional arguments for run_lucj_qsci	{}

Returns:

Type	Description
LUCJQSCIResult	LUCJQSCIResult with convergence data

Source code in ffsim_integration/qsci_interface.py

def run_convergence_study(
    molecule: str,
    ansatz_type: str = "LUCJ",
    max_subspace: int = 200,
    step_size: int = 25,
    **kwargs
) -> LUCJQSCIResult:
    """Run a convergence study with increasing subspace sizes.

    Args:
        molecule: Molecule name
        ansatz_type: Ansatz type
        max_subspace: Maximum subspace size to test
        step_size: Step size for subspace sizes
        **kwargs: Additional arguments for run_lucj_qsci

    Returns:
        LUCJQSCIResult with convergence data
    """
    subspace_sizes = list(range(25, max_subspace + 1, step_size))

    print(f"Running convergence study for {ansatz_type} + QSCI")
    print(f"Subspace sizes: {subspace_sizes}")

    return run_lucj_qsci(
        molecule=molecule,
        ansatz_type=ansatz_type,
        subspace_sizes=subspace_sizes,
        **kwargs
    )

benchmark_against_reference

benchmark_against_reference(
    result, reference_energy, tolerance=0.001
)

Benchmark LUCJ/UCJ + QSCI results against reference energy.

Parameters:

Name	Type	Description	Default
result	LUCJQSCIResult	LUCJQSCIResult to benchmark	required
reference_energy	float	Reference energy in Hartree	required
tolerance	float	Energy tolerance in Hartree	0.001

Returns:

Type	Description
Dict[str, bool]	Dictionary with benchmark results

Source code in ffsim_integration/qsci_interface.py

def benchmark_against_reference(
    result: LUCJQSCIResult, 
    reference_energy: float,
    tolerance: float = 0.001
) -> Dict[str, bool]:
    """Benchmark LUCJ/UCJ + QSCI results against reference energy.

    Args:
        result: LUCJQSCIResult to benchmark
        reference_energy: Reference energy in Hartree
        tolerance: Energy tolerance in Hartree

    Returns:
        Dictionary with benchmark results
    """
    best_energy = result.best_qsci_energy
    energy_error = abs(best_energy - reference_energy)

    benchmark = {
        "energy_within_tolerance": energy_error < tolerance,
        "convergence_demonstrated": len(result.qsci_results) > 1,
        "ansatz_improves_hf": result.ansatz_energy < result.hartree_fock_energy,
        "qsci_improves_ansatz": result.best_qsci_energy < result.ansatz_energy,
        "reasonable_conversion_fidelity": result.conversion_metrics.fidelity > 0.5
    }

    print(f"\nBENCHMARK RESULTS vs {reference_energy:.6f} Ha:")
    for criterion, passed in benchmark.items():
        status = "✓" if passed else "✗"
        print(f"  {status} {criterion.replace('_', ' ').title()}")

    all_passed = all(benchmark.values())
    print(f"\nOverall benchmark: {'✓ PASSED' if all_passed else '✗ FAILED'}")

    return benchmark