From 0e4b4c260943ed3608b0f87142d751e1d00bcf0a Mon Sep 17 00:00:00 2001 From: Krisztina Zsigmond <144738961+kzsigmond@users.noreply.github.com> Date: Thu, 30 Apr 2026 14:51:48 -0400 Subject: [PATCH 1/7] chore: remove psi4 wrapper --- fanpy/tools/wrapper/psi4.py | 65 -------------------------------- tests/test_tools_wrapper_psi4.py | 32 ---------------- 2 files changed, 97 deletions(-) delete mode 100644 fanpy/tools/wrapper/psi4.py delete mode 100644 tests/test_tools_wrapper_psi4.py diff --git a/fanpy/tools/wrapper/psi4.py b/fanpy/tools/wrapper/psi4.py deleted file mode 100644 index eeadec3e..00000000 --- a/fanpy/tools/wrapper/psi4.py +++ /dev/null @@ -1,65 +0,0 @@ -"""Script for utilizing Psi4.""" -import numpy as np - -import psi4 - - -def hartreefock(xyz_file, basis, is_unrestricted=False): # pylint: disable=W0613 - """Run HF using PySCF. - - Parameters - ---------- - xyz_file : str - XYZ file location. - Units are in Angstrom. - basis : str - Basis set available in PySCF. - is_unrestricted : bool - Flag to run unrestricted HF. - Default is restricted HF. - - Returns - ------- - result : dict - "hf_energy" - The electronic energy. - "nuc_nuc" - The nuclear repulsion energy. - "one_int" - The tuple of the one-electron interal. - "two_int" - The tuple of the two-electron integral in Physicist's notation. - - Raises - ------ - ValueError - If given xyz file does not exist. - NotImplementedError - If calculation is unrestricted or generalized. - - """ - # pylint: disable=E1101,C0103 - # get coordinates - with open(xyz_file, "r") as fh: - lines = [i.strip() for i in fh.readlines()[2:]] - geometry = "\n".join(lines + ["symmetry c1"]) - - mol = psi4.geometry(geometry) - psi4.set_options({"basis": basis, "reference": "rhf", "scf_type": "direct"}) - hf_energy, hf_wfn = psi4.energy("scf", return_wfn=True) - C = hf_wfn.Ca() # noqa: N806 - npC = C.to_array() # noqa: N806 - - mints = psi4.core.MintsHelper(hf_wfn.basisset()) - one_int = np.asarray(mints.ao_kinetic()) + np.asarray(mints.ao_potential()) - one_int = np.einsum("uj,vi,uv", npC, npC, one_int) - two_int = np.asarray(mints.mo_eri(C, C, C, C)) - two_int = two_int.swapaxes(1, 2) - - result = { - "hf_energy": hf_energy, - "nuc_nuc": mol.nuclear_repulsion_energy(), - "one_int": one_int, - "two_int": two_int, - } - return result diff --git a/tests/test_tools_wrapper_psi4.py b/tests/test_tools_wrapper_psi4.py deleted file mode 100644 index 9a9becab..00000000 --- a/tests/test_tools_wrapper_psi4.py +++ /dev/null @@ -1,32 +0,0 @@ -"""Test fanpy.tools.wrapper.psi4.""" -from fanpy.ham.restricted_chemical import RestrictedMolecularHamiltonian -from fanpy.solver.ci import brute -from fanpy.wfn.ci.fci import FCI - -import pytest - -from utils import find_datafile - - -def test_generate_fci_cimatrix_h2_631gdp(): - """Test PySCF's FCI calculation against H2 FCI 6-31G** data from Gaussian. - - HF energy: -1.13126983927 - FCI energy: -1.1651487496 - - """ - pytest.importorskip("psi4") - from fanpy.tools.wrapper.psi4 import hartreefock - - hf_data = hartreefock(find_datafile("data/data_h2.xyz"), "6-31gss") - - nelec = 2 - nuc_nuc = hf_data["nuc_nuc"] - one_int = hf_data["one_int"] - two_int = hf_data["two_int"] - - wfn = FCI(nelec, one_int.shape[0] * 2) - ham = RestrictedMolecularHamiltonian(one_int, two_int) - results = brute(wfn, ham) - energy = results["energy"] - assert abs(energy + nuc_nuc - (-1.1651486697)) < 1e-7 From c54c6ca61931817154a5eab1e4129a6e880d2e4c Mon Sep 17 00:00:00 2001 From: Krisztina Zsigmond <144738961+kzsigmond@users.noreply.github.com> Date: Thu, 30 Apr 2026 14:56:33 -0400 Subject: [PATCH 2/7] chore: remove python wrapper --- fanpy/tools/wrapper/python_wrapper.py | 132 -------------------------- tests/test_wrapper_python_wrapper.py | 73 +------------- 2 files changed, 1 insertion(+), 204 deletions(-) delete mode 100644 fanpy/tools/wrapper/python_wrapper.py diff --git a/fanpy/tools/wrapper/python_wrapper.py b/fanpy/tools/wrapper/python_wrapper.py deleted file mode 100644 index bd4d372b..00000000 --- a/fanpy/tools/wrapper/python_wrapper.py +++ /dev/null @@ -1,132 +0,0 @@ -"""Wraps appropriate Python versions around current Python version (3.6+). - -Because HORTON does not support many versions of Python, some sort of workaround is needed. Since -only certain arrays are needed from HORTON, these arrays are extracted and saved to disk as a `.npy` -file. Then, user can load these files to have access to these numbers. - -Notes ------ -This is only a temporary solution. We can make things back-compatible to certain versions of HORTON, -or we can wait until these modules catch up to the latest Python versions. For now, this module can -act as a temporary hack to access these modules. - -""" -import os -import sys -from subprocess import call # nosec:B404 - -import numpy as np - -DIRNAME = os.path.dirname(os.path.abspath(__file__)) - - -def generate_hartreefock_results( - calctype, - energies_name="energies.npy", - oneint_name="oneint.npy", - twoint_name="twoint.npy", - remove_npyfiles=False, - **kwargs -): # pragma: no cover - """Extract results from a Hartree Fock calculation. - - Parameters - ---------- - calctype : {'horton_hartreefock.py', 'horton_gaussian_fchk.py'} - Name of the python script to be used. - energies_name : {str, 'energies.npy} - Name of the file to be generated that contains the electronic and nuclear-nuclear repulsion - energy. - First entry is the electronic energy. - Second entry is the nuclear-nuclear repulsion energy. - oneint_name : {str, 'oneint.npy} - Name of the file to be generated that contains the one electron integrals. - If two-dimensional matrix, then the orbitals are restricted. - If three-dimensional matrix, then the orbitals are unrestricted. - twoint_name : {str, 'twoint.npy} - Name of the file to be generated that contains the two electron integrals. - If four-dimensional matrix, then the orbitals are restricted. - If five-dimensional matrix, then the orbitals are unrestricted. - remove_npyfiles : bool - Option to remove generated numpy files. - True will remove numpy files. - kwargs - Keyword arguments for the script. - - Returns - ------- - el_energy : float - Electronic energy. - nuc_nuc_energy : float - Nuclear-nuclear repulsion energy. - oneint : {np.ndarray, tuple of np.ndarray} - One electron integrals. - If numpy array, then orbitals are restricted. - If tuple of numpy arrays, then orbitals are unrestricted. - twoint : {np.ndarray, tuple of np.ndarray} - Two electron integrals. - If numpy array, then orbitals are restricted. - If tuple of numpy arrays, then orbitals are unrestricted. - - Raises - ------ - ValueError - If calctype is not one of 'horton_hartreefock.py' or 'horton_gaussian_fchk.py' - - Note - ---- - Since HORTON is available only for Python 2.7, it is not compatible with this module. This - function executes the given scripts using the Python interpreter that is compatible with HORTOn, - provided via an environment variable `HORTONPYTHON`. However, this function is a security risk - since we do not actually check that the file provided by the environment variable is a Python - interpreter. This means that **it is up to the user to specify the correct path for the Python - interpreter**. - - """ - # get python interpreter - try: - if calctype in ["horton_hartreefock.py", "horton_gaussian_fchk.py"]: - python_name = os.environ["HORTONPYTHON"] - else: - raise ValueError( - "The calctype must be one of 'horton_hartreefock.py' or 'horton_gaussian_fchk.py'." - ) - except KeyError: - python_name = sys.executable - if not os.path.isfile(python_name): - raise FileNotFoundError( - "The file provided in HORTONPYTHON environment variable is not a python executable." - ) - # FIXME: I can't think of a way to make sure that the python_name is a python interpreter. - - # turn keywords to pair of key and value - kwargs = [str(i) for item in kwargs.items() for i in item] - # call script with appropriate python - # NOTE: this is a possible security risk since we don't check that python_name is actually a - # python interpreter. However, it is up to the user to make sure that their environment variable - # is set up properly. Ideally, we shouldn't even have to call a python outside the current - # python, but here we are. - call( # nosec: B603 - [ - python_name, - os.path.join(DIRNAME, calctype), - energies_name, - oneint_name, - twoint_name, - *kwargs, - ] - ) - el_energy, nuc_nuc_energy = np.load(energies_name) - oneint = np.load(oneint_name) - if oneint.ndim == 3: - oneint = tuple(oneint) - twoint = np.load(twoint_name) - if twoint.ndim == 5: - twoint = tuple(twoint) - - if remove_npyfiles: - os.remove(energies_name) - os.remove(oneint_name) - os.remove(twoint_name) - - return el_energy, nuc_nuc_energy, oneint, twoint diff --git a/tests/test_wrapper_python_wrapper.py b/tests/test_wrapper_python_wrapper.py index fbf6ede3..e62bd63c 100644 --- a/tests/test_wrapper_python_wrapper.py +++ b/tests/test_wrapper_python_wrapper.py @@ -3,14 +3,8 @@ import sys from subprocess import call -from fanpy.tools.wrapper.python_wrapper import generate_hartreefock_results - import numpy as np -import pytest - -from utils import find_datafile - def check_data_h2_rhf_sto6g(el_energy, nuc_nuc_energy, one_int, two_int): """Check data for h2 rhf sto6g calculation.""" @@ -123,69 +117,4 @@ def check_dependency(dependency): # is set up properly. Ideally, we shouldn't even have to call a python outside the current # python, but here we are. exit_code = call([python_name, "-c", "import {0}".format(dependency)]) - return exit_code == 0 - - -def test_generate_hartreefock_results_error(): - """Test python_wrapper.generate_hartreefock_results. - - Check if it raises correct error. - - """ - with pytest.raises(ValueError): - generate_hartreefock_results("sdf") - os.environ["HORTONPYTHON"] = "Asfasdfsadf" - with pytest.raises(FileNotFoundError): - generate_hartreefock_results("horton_hartreefock.py") - - -@pytest.mark.skipif( - not check_dependency("horton"), reason="HORTON is not available or HORTONPATH is not set." -) -def test_horton_hartreefock_h2_rhf_sto6g(): - """Test HORTON"s hartreefock against H2 HF STO-6G data from Gaussian.""" - hf_data = generate_hartreefock_results( - "horton_hartreefock.py", - energies_name="energies.npy", - oneint_name="oneint.npy", - twoint_name="twoint.npy", - remove_npyfiles=True, - fn=find_datafile("data/data_h2.xyz"), - basis="sto-6g", - nelec=2, - ) - check_data_h2_rhf_sto6g(*hf_data) - - -@pytest.mark.skipif( - not check_dependency("horton"), reason="HORTON is not available or HORTONPATH is not set." -) -def test_horton_gaussian_fchk_h2_rhf_sto6g(): - """Test HORTON"s gaussian_fchk against H2 HF STO-6G data from Gaussian.""" - fchk_data = generate_hartreefock_results( - "horton_gaussian_fchk.py", - energies_name="energies.npy", - oneint_name="oneint.npy", - twoint_name="twoint.npy", - remove_npyfiles=True, - fchk_file=find_datafile("data/data_h2_hf_sto6g.fchk"), - horton_internal=False, - ) - check_data_h2_rhf_sto6g(*fchk_data) - - -@pytest.mark.skipif( - not check_dependency("horton"), reason="HORTON is not available or HORTONPATH is not set." -) -def test_gaussian_fchk_h2_uhf_sto6g(): - """Test HORTON"s gaussian_fchk against H2 UHF STO-6G data from Gaussian.""" - fchk_data = generate_hartreefock_results( - "horton_gaussian_fchk.py", - energies_name="energies.npy", - oneint_name="oneint.npy", - twoint_name="twoint.npy", - remove_npyfiles=True, - fchk_file=find_datafile("data/data_h2_uhf_sto6g.fchk"), - horton_internal=False, - ) - check_data_h2_uhf_sto6g(*fchk_data) + return exit_code == 0 \ No newline at end of file From be2c25e11127dedc79bf1ea7db81437334250b47 Mon Sep 17 00:00:00 2001 From: Krisztina Zsigmond <144738961+kzsigmond@users.noreply.github.com> Date: Thu, 30 Apr 2026 14:58:50 -0400 Subject: [PATCH 3/7] chore: remove horton wrappers --- fanpy/tools/wrapper/horton_gaussian_fchk.py | 118 ------------- fanpy/tools/wrapper/horton_hartreefock.py | 180 -------------------- 2 files changed, 298 deletions(-) delete mode 100644 fanpy/tools/wrapper/horton_gaussian_fchk.py delete mode 100644 fanpy/tools/wrapper/horton_hartreefock.py diff --git a/fanpy/tools/wrapper/horton_gaussian_fchk.py b/fanpy/tools/wrapper/horton_gaussian_fchk.py deleted file mode 100644 index b25c0ef2..00000000 --- a/fanpy/tools/wrapper/horton_gaussian_fchk.py +++ /dev/null @@ -1,118 +0,0 @@ -"""Script for generating one- and two-electron integrals from Gaussian using HORTON. - -Functions ---------- -gaussian_fchk(fchk_file, horton_internal=False, compute_nuc=True) - Extract appropriate information from a Gaussian FCHK file. - -""" -# pylint: disable=W0212,C0103 -import sys - -from horton import IOData, compute_nucnuc - -import numpy as np - - -# FIXME: remove horton_internal -def gaussian_fchk(fchk_file, horton_internal=False): - """Extract the appropriate data from Gaussian fchk file (using HORTON). - - Parameters - ---------- - fchk_file : str - Formatted chk file. - horton_internal : bool - Flag to return horton_internal variables. - - Returns - ------- - result : dict - "el_energy" - Electronic energy. - "nuc_nuc_energy" - Nuclear nuclear repulsion energy. - "one_int" - Tuple of the one-electron Hamiltonian. - "two_int" - Tuple of the two-electron Hamiltonian. - - Raises - ------ - NotImplementedError - If `horton_internal` is True. - - """ - mol = IOData.from_file(fchk_file) - - # for spin orbitals - exps = [mol.exp_alpha] - if hasattr(mol, "exp_beta"): - exps.append(mol.exp_beta) - - obasis = mol.obasis - kin = obasis.compute_kinetic(mol.lf)._array - na = obasis.compute_nuclear_attraction(mol.coordinates, mol.pseudo_numbers, mol.lf)._array - - one_ab = kin - one_ab += na - two_ab = obasis.compute_electron_repulsion(mol.lf)._array - - # compute nuclear nuclear repulsion - nuc_nuc = compute_nucnuc(mol.coordinates, mol.pseudo_numbers) - - # for spin orbitals - one_mo = [] - two_mo = [] - for i, exp_i in enumerate(exps): - for j, exp_j in enumerate(exps[i:]): - j += i - temp = np.einsum("sd,pqrs->pqrd", exp_j.coeffs, two_ab, casting="no", order="C") - temp = np.einsum("rc,pqrd->pqcd", exp_i.coeffs, temp, casting="no", order="C") - temp = np.einsum("qb,pqcd->pbcd", exp_j.coeffs, temp, casting="no", order="C") - temp = np.einsum("pa,pbcd->abcd", exp_i.coeffs, temp, casting="no", order="C") - two_mo.append(temp) - one_mo.append(exp_i.coeffs.T.dot(one_ab).dot(exp_i.coeffs)) - - output = { - "el_energy": mol.energy - nuc_nuc, - "nuc_nuc_energy": nuc_nuc, - "one_int": tuple(one_mo), - "two_int": tuple(two_mo), - "t_ab_mo": tuple(i.coeffs for i in exps), - } - if horton_internal: - raise NotImplementedError( - "horton_internal storage is unsupported until it become Python " "3.6 compatible." - ) - # output["horton_internal"] = { - # "mol": mol, - # "lf": mol.lf, - # "one": one_mo, - # "two": two_mo, - # "orb": exps, - # } - return output - - -if __name__ == "__main__": - # extract keyword from command line - kwargs = {key: val for key, val in zip(sys.argv[5::2], sys.argv[6::2])} # pylint: disable=R1721 - # change data types - if "horton_internal" in kwargs: - kwargs["horton_internal"] = kwargs["horton_internal"] == "True" - - data = gaussian_fchk(**kwargs) - np.save(sys.argv[1], [data["el_energy"], data["nuc_nuc_energy"]]) - if len(data["one_int"]) == 1: - np.save(sys.argv[2], data["one_int"][0]) - else: - np.save(sys.argv[2], data["one_int"]) - if len(data["two_int"]) == 1: - np.save(sys.argv[3], data["two_int"][0]) - else: - np.save(sys.argv[3], data["two_int"]) - if len(data["t_ab_mo"]) == 1: - np.save(sys.argv[4], data["t_ab_mo"][0]) - else: - np.save(sys.argv[4], data["t_ab_mo"]) diff --git a/fanpy/tools/wrapper/horton_hartreefock.py b/fanpy/tools/wrapper/horton_hartreefock.py deleted file mode 100644 index a4887568..00000000 --- a/fanpy/tools/wrapper/horton_hartreefock.py +++ /dev/null @@ -1,180 +0,0 @@ -"""Script for generating one and two electron integrals using HORTON. - -Functions ---------- -hartreefock(xyz_file, basis, is_unrestricted=False) - Runs HF in HORTON. - -""" -# pylint: disable=W0212,C0103,W0621 -import sys - -from horton import ( - AufbauOccModel, - DenseLinalgFactory, - EDIIS2SCFSolver, - IOData, - PlainSCFSolver, - RDirectTerm, - REffHam, - RExchangeTerm, - RTwoIndexTerm, - compute_nucnuc, - get_gobasis, - guess_core_hamiltonian, - transform_integrals, -) - -import numpy as np - - -def hartreefock( - fn=None, - basis=None, - nelec=None, - solver=EDIIS2SCFSolver, - tol=1.0e-12, - horton_internal=False, - **kwargs -): - """Run a HF calculation using HORTON. - - Parameters - ---------- - fn : str - File name. - Supports XYZ file (in Angstrom). - basis : str - Basis set name. - nelec : int - Number of electrons. - solver : {PlainSCFSolver, ODASSCFSolver, CDIISSCFSolver, EDIISSCFSolver, EDIIS2SCFSolver} - HORTON's SCF solver. - tol : float - The convergence threshold for the wavefunction. - horton_internal : bool - Flag for returning HORTON's internal data. - kwargs - Other optional arguments for the SCF Solver. - - Returns - ------- - result : dict - "el_energy" - Electronic energy. - "nuc_nuc_energy" - Nuclear repulsion energy. - "one_int" - Tuple of the one-electron Hamiltonian. - "two_int" - Tuple of the two-electron Hamiltonian. - - Raises - ------ - ValueError - If solver is not callable. - NotImplemntedError - If `horton_internal` is True. - - Notes - ----- - While HORTON does support unrestricted calculations, the wrapper does not support it. - - """ - # Initialize molecule and basis set from specified file - if isinstance(fn, IOData): - mol = fn - else: - mol = IOData.from_file(fn) - obasis = get_gobasis(mol.coordinates, mol.numbers, basis) - npair = nelec // 2 - - # Fill in orbital expansion and overlap - occ_model = AufbauOccModel(npair) - lf = DenseLinalgFactory(obasis.nbasis) - orb = lf.create_expansion(obasis.nbasis) - olp = obasis.compute_overlap(lf) - - # Construct Hamiltonian and density matrix - kin = obasis.compute_kinetic(lf) - na = obasis.compute_nuclear_attraction(mol.coordinates, mol.pseudo_numbers, lf) - two = obasis.compute_electron_repulsion(lf) - external = {"nn": compute_nucnuc(mol.coordinates, mol.pseudo_numbers)} - terms = [ - RTwoIndexTerm(kin, "kin"), - RDirectTerm(two, "hartree"), - RExchangeTerm(two, "x_hf"), - RTwoIndexTerm(na, "ne"), - ] - ham = REffHam(terms, external) - guess_core_hamiltonian(olp, kin, na, orb) - - # Do Hartree-Fock SCF - if solver is PlainSCFSolver: - PlainSCFSolver(tol)(ham, lf, olp, occ_model, orb, **kwargs) - elif hasattr(solver, "__call__"): - occ_model.assign(orb) - dm = orb.to_dm() - solver(tol)(ham, lf, olp, occ_model, dm, **kwargs) - # transform back to orbitals - fock_alpha = lf.create_two_index() - ham.compute_fock(fock_alpha) - orb.from_fock_and_dm(fock_alpha, dm, olp) - else: - raise ValueError("Given solver, {0}, is not callable".format(solver)) - energy = ham.cache["energy"] - - # Transform one- and two- electron integrals into MO basis - one = kin - one.iadd(na) - one_mo, two_mo = transform_integrals(one, two, "tensordot", orb) - one_mo = tuple(i._array for i in one_mo) - two_mo = tuple(i._array for i in two_mo) - - output = { - # Collect HF energy and integrals - "el_energy": energy - external["nn"], - "nuc_nuc_energy": external["nn"], - "one_int": one_mo, - "two_int": two_mo, - } - - if horton_internal: - output["horton_internal"] = { - "lf": lf, - "occ_model": occ_model, - "one": one, - "two": two, - "orb": [orb], - "olp": olp, - } - raise NotImplementedError( - "horton_internal storage is unsupported until it become Python " "3.6 compatible." - ) - - return output - - -if __name__ == "__main__": - # extract keyword from command line - kwargs = {key: val for key, val in zip(sys.argv[4::2], sys.argv[5::2])} # pylint: disable=R1721 - # change data types - if "nelec" in kwargs: - kwargs["nelec"] = int(kwargs["nelec"]) - if "tol" in kwargs: - kwargs["tol"] = float(kwargs["tol"]) - if "horton_internal" in kwargs: - kwargs["horton_internal"] = bool(kwargs["horton_internal"]) - if "solver" in kwargs: - kwargs["solver"] = locals()[kwargs["solver"]] - - data = hartreefock(**kwargs) - np.save(sys.argv[1], [data["el_energy"], data["nuc_nuc_energy"]]) - if len(data["one_int"]) == 1: - np.save(sys.argv[2], data["one_int"][0]) - else: - np.save(sys.argv[2], data["one_int"]) - if len(data["two_int"]) == 1: - np.save(sys.argv[3], data["two_int"][0]) - else: - np.save(sys.argv[3], data["two_int"]) From 36a4da3f35cca59ff2b5e7096fa59e9757010b81 Mon Sep 17 00:00:00 2001 From: Krisztina Zsigmond <144738961+kzsigmond@users.noreply.github.com> Date: Thu, 30 Apr 2026 14:59:41 -0400 Subject: [PATCH 4/7] chore: remove trial loc script this script does not have any functions or classes that fanpy relies on. It also does not use any fanpy functionality. It is likely a remnant of past development that should not be part of Fanpy. --- fanpy/tools/wrapper/trial_loc.py | 132 ------------------------------- 1 file changed, 132 deletions(-) delete mode 100644 fanpy/tools/wrapper/trial_loc.py diff --git a/fanpy/tools/wrapper/trial_loc.py b/fanpy/tools/wrapper/trial_loc.py deleted file mode 100644 index 20c36b71..00000000 --- a/fanpy/tools/wrapper/trial_loc.py +++ /dev/null @@ -1,132 +0,0 @@ -from functools import reduce -import numpy as np -import scipy.linalg -import os -import re -from pyscf import gto -from pyscf import scf -from pyscf import __config__ -from pyscf.lo.iao import reference_mol -from fanpy.interface.pyscf_tools import convert_gbs_nwchem - - -xyz_file = "./system.xyz" -basis_file = "../basis.gbs" -system_inds = [0, 1, 2] -nelec = 6 - -# check xyz file -cwd = os.path.dirname(__file__) -if os.path.isfile(os.path.join(cwd, xyz_file)): - xyz_file = os.path.join(cwd, xyz_file) -elif not os.path.isfile(xyz_file): # pragma: no branch - raise ValueError("Given xyz_file does not exist") - -# get coordinates -with open(xyz_file, "r") as f: - lines = [i.strip() for i in f.readlines()[2:]] - atoms = ";".join(lines) - -# get mol -if os.path.splitext(basis_file)[1] == ".gbs": - basis = convert_gbs_nwchem(basis_file) -else: - basis = basis_file -mol = gto.M( - atom=atoms, basis={i: j for i, j in basis.items() if i + ' ' in atoms}, parse_arg=False, - unit="bohr" -) - -# get hf -hf = scf.RHF(mol) -# run hf -hf.scf() -mo_coeff = hf.mo_coeff - - -MINAO = getattr(__config__, 'lo_iao_minao', 'minao') - -pmol = reference_mol(mol, MINAO) -orig_s12 = gto.mole.intor_cross('int1e_ovlp', mol, pmol) - -minao_labels = pmol.ao_labels() -minao_inds = np.array([system_inds[int(re.search(r'^(\d+)\s+', label).group(1))] for label in minao_labels]) -print(minao_labels) -print(minao_inds) - -coeff_mo_lo = [] -for sub_mo_coeff in [mo_coeff[:, hf.mo_occ > 0], mo_coeff[:, hf.mo_occ == 0]]: - s12 = sub_mo_coeff.T.dot(orig_s12) - system_s12 = [] - for i in range(max(system_inds) + 1): - system_s12.append(s12[:, minao_inds == i]) - - system_mo_lo_T = [[] for i in range(max(system_inds) + 1)] - counter = 0 - cum_transform = np.identity(s12.shape[0]) - while counter < sub_mo_coeff.shape[1]: - system_u = [] - system_s = [] - print('x'*99) - for s12 in system_s12: - u, s, vdag = np.linalg.svd(s12) - system_u.append(u.T) - system_s.append(s) - # TEST: negative singular value? - assert np.all(s > 0) - print(s12) - print(s) - # find which system had the largest singular value - max_system_ind = np.argmax([max(s) for s in system_s]) - # find largest singular value - max_sigma_ind = np.argmax(system_s[max_system_ind]) - # add corresponding left singular vector - system_mo_lo_T[max_system_ind].append(system_u[max_system_ind][max_sigma_ind].dot(cum_transform)) - # update overlap matrix (remove singular vector) - trunc_transform = np.delete(system_u[max_system_ind], max_sigma_ind, axis=0) - assert np.allclose(trunc_transform.dot(trunc_transform.T), np.identity(trunc_transform.shape[0])) - for i in range(len(system_s12)): - system_s12[i] = trunc_transform.dot(system_s12[i]) - # update transformation matrix - cum_transform = trunc_transform.dot(cum_transform) - # increment - counter += 1 - - lo_inds = [[i] * len(rows) for i, rows in enumerate(system_mo_lo_T)] - lo_inds = np.array([j for i in lo_inds for j in i]) - print(lo_inds) - system_mo_lo_T = np.vstack([np.vstack(rows) for rows in system_mo_lo_T if rows]) - np.set_printoptions(linewidth=200) - print(system_mo_lo_T) - - # check orthogonalization - s1 = mol.intor_symmetric('int1e_ovlp') - transform = system_mo_lo_T.dot(sub_mo_coeff.T) - olp = transform.dot(s1).dot(transform.T) - print(olp) - assert np.allclose(olp, np.identity(transform.shape[0])) - - # check span - s12 = gto.mole.intor_cross('int1e_ovlp', mol, pmol) - s12 = sub_mo_coeff.T.dot(s12) - s12 = system_mo_lo_T.dot(s12) - for i in range(max(system_inds) + 1): - # how much is inside - system_s12 = s12[lo_inds == i][:, minao_inds == i] - print(np.sum(np.diag(system_s12.T.dot(system_s12)))) - # how much is outside - system_s12 = s12[lo_inds == i][:, minao_inds != i] - print(np.sum(np.diag(system_s12.T.dot(system_s12)))) - - coeff_mo_lo.append(system_mo_lo_T.T) - -from pyscf.tools import molden -# visualize -temp = np.zeros((coeff_mo_lo[0].shape[0] + coeff_mo_lo[1].shape[0], coeff_mo_lo[0].shape[1] + coeff_mo_lo[1].shape[1])) -temp[np.where(hf.mo_occ > 0)[0][:, None], np.where(np.arange(temp.shape[1]) < nelec // 2)[0][None, :]] = coeff_mo_lo[0] -temp[np.where(hf.mo_occ == 0)[0][:, None], np.where(np.arange(temp.shape[1]) >= nelec // 2)[0][None, :]] = coeff_mo_lo[1] -coeff_mo_lo = scipy.linalg.block_diag(*coeff_mo_lo) -assert np.allclose(coeff_mo_lo, temp) -t_ab_lo = mo_coeff.dot(coeff_mo_lo) -molden.from_mo(mol, 'trial.molden', t_ab_lo) - From 585467c5a25ff1e38eb4bca7779448d1cd6544ab Mon Sep 17 00:00:00 2001 From: Krisztina Zsigmond <144738961+kzsigmond@users.noreply.github.com> Date: Thu, 30 Apr 2026 15:02:13 -0400 Subject: [PATCH 5/7] chore: remove init of wrappers --- fanpy/tools/wrapper/__init__.py | 1 - 1 file changed, 1 deletion(-) delete mode 100644 fanpy/tools/wrapper/__init__.py diff --git a/fanpy/tools/wrapper/__init__.py b/fanpy/tools/wrapper/__init__.py deleted file mode 100644 index 702ba85f..00000000 --- a/fanpy/tools/wrapper/__init__.py +++ /dev/null @@ -1 +0,0 @@ -"""Module for wrapping other quantum chemistry python modules.""" From 6ea6b766f810ce975e198f151f96291fb634e20e Mon Sep 17 00:00:00 2001 From: Krisztina Zsigmond <144738961+kzsigmond@users.noreply.github.com> Date: Thu, 30 Apr 2026 15:08:03 -0400 Subject: [PATCH 6/7] chore: rename python wrapper test. This only houses utility functions used in the PySCF interface --- tests/{test_wrapper_python_wrapper.py => interface_utils.py} | 0 tests/test_interface_pyscf.py | 2 +- tests/test_interface_pyscf_tools.py | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) rename tests/{test_wrapper_python_wrapper.py => interface_utils.py} (100%) diff --git a/tests/test_wrapper_python_wrapper.py b/tests/interface_utils.py similarity index 100% rename from tests/test_wrapper_python_wrapper.py rename to tests/interface_utils.py diff --git a/tests/test_interface_pyscf.py b/tests/test_interface_pyscf.py index 3a520a50..0e8dac9b 100644 --- a/tests/test_interface_pyscf.py +++ b/tests/test_interface_pyscf.py @@ -4,7 +4,7 @@ from pyscf import gto, scf from fanpy.interface.pyscf import PYSCF -from test_wrapper_python_wrapper import check_data_h2_rhf_sto6g +from tests.interface_utils import check_data_h2_rhf_sto6g def test_hf_data_h2_sto6g(): """Make sure computed integrals, and extracted HF energies are the same as Gaussian.""" diff --git a/tests/test_interface_pyscf_tools.py b/tests/test_interface_pyscf_tools.py index 6ac0b66a..4a888f2b 100644 --- a/tests/test_interface_pyscf_tools.py +++ b/tests/test_interface_pyscf_tools.py @@ -5,7 +5,7 @@ from fanpy.interface.pyscf_tools import hartreefock, fci_cimatrix, localize -from test_wrapper_python_wrapper import check_data_h2_rhf_sto6g, check_data_lih_rhf_sto6g +from tests.interface_utils import check_data_h2_rhf_sto6g, check_data_lih_rhf_sto6g from utils import find_datafile From 2bcf2852761ff10ed6087deff9629464aa53e938 Mon Sep 17 00:00:00 2001 From: Krisztina Zsigmond <144738961+kzsigmond@users.noreply.github.com> Date: Thu, 30 Apr 2026 15:15:38 -0400 Subject: [PATCH 7/7] fix: interface util imports --- tests/test_interface_pyscf.py | 2 +- tests/test_interface_pyscf_tools.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_interface_pyscf.py b/tests/test_interface_pyscf.py index 0e8dac9b..ba5a73f2 100644 --- a/tests/test_interface_pyscf.py +++ b/tests/test_interface_pyscf.py @@ -4,7 +4,7 @@ from pyscf import gto, scf from fanpy.interface.pyscf import PYSCF -from tests.interface_utils import check_data_h2_rhf_sto6g +from interface_utils import check_data_h2_rhf_sto6g def test_hf_data_h2_sto6g(): """Make sure computed integrals, and extracted HF energies are the same as Gaussian.""" diff --git a/tests/test_interface_pyscf_tools.py b/tests/test_interface_pyscf_tools.py index 4a888f2b..62fdac23 100644 --- a/tests/test_interface_pyscf_tools.py +++ b/tests/test_interface_pyscf_tools.py @@ -5,7 +5,7 @@ from fanpy.interface.pyscf_tools import hartreefock, fci_cimatrix, localize -from tests.interface_utils import check_data_h2_rhf_sto6g, check_data_lih_rhf_sto6g +from interface_utils import check_data_h2_rhf_sto6g, check_data_lih_rhf_sto6g from utils import find_datafile