diff --git a/.ruff.toml b/.ruff.toml index c79a2d33..e5bb3848 100644 --- a/.ruff.toml +++ b/.ruff.toml @@ -1,6 +1,7 @@ [lint] # see the rules [here](https://docs.astral.sh/ruff/rules/) select = ["E", "F", "I", "NPY", "PL", "ARG"] +exclude = ["tests/fragmentation_test.py"] ignore = [ "S101", # https://docs.astral.sh/ruff/rules/assert/ diff --git a/docs/source/conf.py b/docs/source/conf.py index 6e2813c8..3c366a5b 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -58,6 +58,7 @@ "python": ("https://docs.python.org/3", None), "pyscf": ("https://pyscf.org/", None), "h5py": ("https://docs.h5py.org/en/stable/", None), + "networkx": ("https://networkx.org/documentation/stable/", None), } diff --git a/example/molbe_ppp.py b/example/molbe_ppp.py index 0af4f253..55251ee3 100644 --- a/example/molbe_ppp.py +++ b/example/molbe_ppp.py @@ -32,7 +32,7 @@ mf.kernel() # Define fragments; use IAO scheme with 'sto-3g' as the minimal basis set -fobj = fragpart(be_type="be2", mol=mol, valence_basis="sto-3g", frozen_core=True) +fobj = fragpart(be_type="be2", mol=mol, iao_valence_basis="sto-3g", frozen_core=True) # Initialize BE mybe = BE(mf, fobj, lo_method="iao") diff --git a/mypy.ini b/mypy.ini index 9fb77851..ea7bff8c 100644 --- a/mypy.ini +++ b/mypy.ini @@ -18,7 +18,7 @@ disallow_untyped_defs = False check_untyped_defs = False -[mypy-tests.chem_dm_kBE_test,tests.chempot_molBE_test,tests.dm_molBE_test,tests.dmrg_molBE_test,tests.eri_onthefly_test,tests.hf-in-hf_BE_test,tests.kbe_polyacetylene_test,tests.molbe_h8_test,tests.molbe_io_fcidump_test,tests.molbe_octane_get_rdms_test,tests.molbe_oneshot_rbe_hcore_test,tests.molbe_oneshot_rbe_qmmm-fromchk_test,tests.ube-oneshot_test] +[mypy-tests.fragmentation_test,tests.chem_dm_kBE_test,tests.chempot_molBE_test,tests.dm_molBE_test,tests.dmrg_molBE_test,tests.eri_onthefly_test,tests.hf-in-hf_BE_test,tests.kbe_polyacetylene_test,tests.molbe_h8_test,tests.molbe_io_fcidump_test,tests.molbe_octane_get_rdms_test,tests.molbe_oneshot_rbe_hcore_test,tests.molbe_oneshot_rbe_qmmm-fromchk_test,tests.ube-oneshot_test] disallow_untyped_defs = False check_untyped_defs = False diff --git a/setup.py b/setup.py index 428febd2..7bef7d49 100644 --- a/setup.py +++ b/setup.py @@ -22,6 +22,7 @@ "numpy>=1.22.0", "scipy>=1.7.0", "pyscf>=2.0.0", + "networkx", "matplotlib", "libdmet @ git+https://github.com/gkclab/libdmet_preview.git", "attrs", diff --git a/src/quemb/kbe/autofrag.py b/src/quemb/kbe/autofrag.py index fc62dcb3..62bedd53 100644 --- a/src/quemb/kbe/autofrag.py +++ b/src/quemb/kbe/autofrag.py @@ -241,7 +241,7 @@ def autogen( nx=False, ny=False, nz=False, - valence_basis=None, + iao_valence_basis=None, interlayer=False, print_frags=True, ): @@ -274,7 +274,7 @@ def autogen( write_geom : bool, optional Whether to write a 'fragment.xyz' file which contains all the fragments in Cartesian coordinates. Defaults to False. - valence_basis : str, optional + iao_valence_basis : str, optional Name of minimal basis set for IAO scheme. 'sto-3g' is sufficient for most cases. Defaults to None. valence_only : bool, optional @@ -1922,11 +1922,11 @@ def autogen( ) w.close() - pao = valence_basis is not None + pao = iao_valence_basis is not None if pao: cell2 = cell.copy() - cell2.basis = valence_basis + cell2.basis = iao_valence_basis cell2.build() bas2list = cell2.aoslice_by_atom() diff --git a/src/quemb/kbe/fragment.py b/src/quemb/kbe/fragment.py index 167d95c2..b27d7e0a 100644 --- a/src/quemb/kbe/fragment.py +++ b/src/quemb/kbe/fragment.py @@ -23,7 +23,7 @@ def __init__( ny=False, nz=False, kpt=None, - valence_basis=None, + iao_valence_basis=None, be_type="be2", mol=None, frozen_core=False, @@ -56,7 +56,7 @@ def __init__( mol : pyscf.pbc.gto.cell.Cell pyscf.pbc.gto.cell.Cell object. This is required for the options, 'autogen', and 'chain' as frag_type. - valence_basis: str + iao_valence_basis: str Name of minimal basis set for IAO scheme. 'sto-3g' suffice for most cases. frozen_core: bool Whether to invoke frozen core approximation. This is set to False by default @@ -90,7 +90,7 @@ def __init__( self.frozen_core = frozen_core self.self_match = self_match self.allcen = allcen - self.valence_basis = valence_basis + self.iao_valence_basis = iao_valence_basis self.kpt = kpt self.molecule = False # remove this @@ -117,7 +117,7 @@ def __init__( kpt, be_type=be_type, frozen_core=frozen_core, - valence_basis=valence_basis, + iao_valence_basis=iao_valence_basis, unitcell=unitcell, nx=nx, ny=ny, diff --git a/src/quemb/kbe/lo.py b/src/quemb/kbe/lo.py index 4fca07b4..03e7c228 100644 --- a/src/quemb/kbe/lo.py +++ b/src/quemb/kbe/lo.py @@ -45,7 +45,7 @@ class Mixin_k_Localize: def localize( self, lo_method, - valence_basis="sto-3g", + iao_valence_basis="sto-3g", core_basis="sto-3g", iao_wannier=True, iao_val_core=True, @@ -60,7 +60,7 @@ def localize( lo_method : str Localization method in quantum chemistry. 'lowdin', 'boys','iao', and 'wannier' are supported. - valence_basis : str + iao_valence_basis : str Name of valence basis set for IAO scheme. 'sto-3g' suffice for most cases. core_basis : str Name of core basis set for IAO scheme. 'sto-3g' suffice for most cases. @@ -128,7 +128,7 @@ def localize( elif lo_method == "iao": if not iao_val_core or not self.frozen_core: Co = self.C[:, :, : self.Nocc].copy() - S12, S2 = get_xovlp_k(self.cell, self.kpts, basis=valence_basis) + S12, S2 = get_xovlp_k(self.cell, self.kpts, basis=iao_valence_basis) ciao_ = get_iao_k(Co, S12, self.S, S2=S2) # tmp - aos are not rearrange and so below is not necessary @@ -145,7 +145,7 @@ def localize( # Cpao = get_pao_k(Ciao, self.S, S12, S2, self.cell) # get_pao_native_k returns symm orthogonalized orbitals cpao_ = get_pao_native_k( - Ciao_, self.S, self.cell, valence_basis, self.kpts + Ciao_, self.S, self.cell, iao_valence_basis, self.kpts ) nk, nao, nlo = cpao_.shape @@ -192,7 +192,7 @@ def localize( # Begin valence s12_val_, s2_val = get_xovlp_k( - self.cell, self.kpts, basis=valence_basis + self.cell, self.kpts, basis=iao_valence_basis ) C_nocore = self.C[:, :, self.ncore :].copy() C_nocore_occ_ = C_nocore[:, :, : self.Nocc].copy() @@ -235,7 +235,12 @@ def localize( ) cpao_ = get_pao_native_k( - c_core_val, self.S, self.cell, valence_basis, self.kpts, ortho=True + c_core_val, + self.S, + self.cell, + iao_valence_basis, + self.kpts, + ortho=True, ) nk, nao, nlo = cpao_.shape Cpao_ = zeros((nk, nao, nlo), dtype=complex128) diff --git a/src/quemb/kbe/lo_k.py b/src/quemb/kbe/lo_k.py index b174f9e7..5f555324 100644 --- a/src/quemb/kbe/lo_k.py +++ b/src/quemb/kbe/lo_k.py @@ -192,7 +192,7 @@ def get_pao_k(Ciao, S, S12): return asarray(Cpao) -def get_pao_native_k(Ciao, S, mol, valence_basis, ortho=True): +def get_pao_native_k(Ciao, S, mol, iao_valence_basis, ortho=True): """ Parameters @@ -203,7 +203,7 @@ def get_pao_native_k(Ciao, S, mol, valence_basis, ortho=True): ao ovlp matrix mol : mol object - valence_basis: + iao_valence_basis: basis used for valence orbitals Returns @@ -215,7 +215,7 @@ def get_pao_native_k(Ciao, S, mol, valence_basis, ortho=True): # Form a mol object with the valence basis for the ao_labels mol_alt = mol.copy() - mol_alt.basis = valence_basis + mol_alt.basis = iao_valence_basis mol_alt.build() full_ao_labels = mol.ao_labels() diff --git a/src/quemb/kbe/pbe.py b/src/quemb/kbe/pbe.py index 583cbc64..d08fed13 100644 --- a/src/quemb/kbe/pbe.py +++ b/src/quemb/kbe/pbe.py @@ -263,7 +263,7 @@ def __init__( # Localize orbitals self.localize( lo_method, - valence_basis=fobj.valence_basis, + iao_valence_basis=fobj.iao_valence_basis, iao_wannier=iao_wannier, iao_val_core=iao_val_core, ) diff --git a/src/quemb/molbe/autofrag.py b/src/quemb/molbe/autofrag.py index 9b52f4fb..4e65dfc5 100644 --- a/src/quemb/molbe/autofrag.py +++ b/src/quemb/molbe/autofrag.py @@ -1,18 +1,298 @@ -# Author: Oinam Romesh Meitei +# Author: Oinam Romesh Meitei, Shaun Weatherly +from typing import Any +import networkx as nx +import numpy as np +from attrs import define +from networkx import single_source_all_shortest_paths # type: ignore[attr-defined] from numpy.linalg import norm +from pyscf import gto from quemb.molbe.helper import get_core from quemb.shared.helper import unused +@define +class FragmentMap: + """Dataclass for fragment bookkeeping. + + Parameters + ---------- + fsites: + List whose entries are tuples containing all AO indices for a fragment. + fs: + List whose entries are tuples of tuples, containing AO indices per atom + per fragment. + edge: + List whose entries are tuples of tuples, containing edge AO + indices per atom (inner tuple) per fragment (outer tuple). + center: + List whose entries are tuples of tuples, containing all fragment AO + indices per atom (inner tuple) and per fragment (outer tuple). + centerf_idx: + List whose entries are tuples containing the relative index of all + center sites within a fragment (ie, with respect to fsites). + ebe_weights: + Weights determining the energy contributions from each center site + (ie, with respect to centerf_idx). + sites: + List whose entries are tuples containing all AO indices per atom + (excluding frozen core indices, if applicable). + dnames: + List of strings giving fragment data names. Useful for bookkeeping and + for constructing fragment scratch directories. + adjacency_mat: + The adjacency matrix for all sites (atoms) in the system. + adjacency_graph: + The adjacency graph corresponding to `adjacency_mat`. + """ + + fsites: list[tuple[int, ...]] + fs: list[tuple[tuple[int, ...], ...]] + edge: list[tuple[tuple[int, ...], ...]] + center: list[tuple[int, ...]] + centerf_idx: list[tuple[int, ...]] + ebe_weights: list[tuple] + sites: list[tuple] + dnames: list + center_atoms: list[tuple[str, ...]] + edge_atoms: list[tuple[str, ...]] + adjacency_mat: np.ndarray + adjacency_graph: nx.Graph + + def remove_nonnunique_frags(self, natm: int) -> None: + """Remove all fragments which are strict subsets of another. + + Remove all fragments whose AO indices can be identified as subsets of + another fragment's. The center site for the removed frag is then + added to that of the superset. Because doing so will necessarily + change the definition of fragments, we repeat it up to `natm` times + such that all fragments are guaranteed to be distinct sets. + another fragment's. The center site for the removed frag is then + added to that of the superset. Because doing so will necessarily + change the definition of fragments, we repeat it up to `natm` times + such that all fragments are guaranteed to be distinct sets. + """ + for _ in range(0, natm): + for adx, basa in enumerate(self.fsites): + for bdx, basb in enumerate(self.fsites): + if adx == bdx: + pass + elif set(basb).issubset(set(basa)): + tmp = set(self.center[adx] + self.center[bdx]) + self.center[adx] = tuple(tmp) + del self.center[bdx] + del self.fsites[bdx] + del self.fs[bdx] + + return None + + +def euclidean_norm( + i_coord: np.ndarray, + j_coord: np.ndarray, +) -> np.floating[Any]: + return norm(np.asarray(i_coord - j_coord)) + + +def graphgen( + mol: gto.Mole, + be_type: str = "BE2", + frozen_core: bool = True, + remove_nonunique_frags: bool = True, + frag_prefix: str = "f", + connectivity: str = "euclidean", + iao_valence_basis: str | None = None, +) -> FragmentMap: + """Generate fragments via adjacency graph. + + Generalizes the BEn fragmentation scheme to arbitrary fragment sizes using a + graph theoretic heuristic. In brief: atoms are assigned to nodes in an + adjacency graph and edges are weighted by some distance metric. For a given + fragment center site, Dijkstra's algorithm is used to find the shortest path + from that center to its neighbors. The number of nodes visited on that shortest + path determines the degree of separation of the corresponding neighbor. I.e., + all atoms whose shortest paths from the center site visit at most 1 node must + be direct neighbors to the center site, which gives BE2-type fragments; all + atoms whose shortest paths visit at most 2 nodes must then be second-order + neighbors, hence BE3; and so on. + + Currently does not support periodic calculations. + + Parameters + ---------- + mol : + The molecule object. + be_type : + The order of nearest neighbors (with respect to the center atom) + included in a fragment. Supports all 'BEn', with 'n' in - + [1, 2, 3, 4, 5, 6, 7, 8, 9] having been tested. + frozen_core: + Whether to exclude core AO indices from the fragmentation process. + True by default. + remove_nonunique_frags: + Whether to remove fragments which are strict subsets of another + fragment in the system. True by default. + frag_prefix: + Prefix to be appended to the fragment datanames. Useful for managing + fragment scratch directories. + connectivity: + Keyword string specifying the distance metric to be used for edge + weights in the fragment adjacency graph. Currently supports "euclidean" + (which uses the square of the distance between atoms in real + space to determine connectivity within a fragment.) + + Returns + ------- + FragmentMap : + FragmentMap mapping various fragment components to AO indices, data names, + and other info. + """ + assert mol is not None + if iao_valence_basis is not None: + raise NotImplementedError("IAOs not yet implemented for graphgen.") + + fragment_type_order = int(be_type[-1]) + natm = mol.natm + + adx_map = { + adx: { + "bas": bas, + "label": mol.atom_symbol(adx), + "coord": mol.atom_coord(adx), + "shortest_paths": dict(), + } + for adx, bas in enumerate(mol.aoslice_by_atom()) + } + + fragment_map = FragmentMap( + fsites=(list(tuple())), + fs=list(tuple(tuple())), + edge=list(tuple(tuple())), + center=list(tuple()), + centerf_idx=list(tuple()), + ebe_weights=list(tuple()), + sites=list(tuple()), + dnames=list(), + center_atoms=list(), + edge_atoms=list(), + adjacency_mat=np.zeros((natm, natm), np.float64), + adjacency_graph=nx.Graph(), + ) + fragment_map.adjacency_graph.add_nodes_from(adx_map) + + _core_offset = 0 + for adx, map in adx_map.items(): + start_ = map["bas"][2] + stop_ = map["bas"][3] + if frozen_core: + _, _, core_list = get_core(mol) + start_ -= _core_offset + ncore_ = int(core_list[adx]) + stop_ -= _core_offset + ncore_ + _core_offset += ncore_ + fragment_map.sites.append(tuple([i for i in range(start_, stop_)])) + else: + fragment_map.sites.append(tuple([i for i in range(start_, stop_)])) + + if connectivity.lower() in ["euclidean_distance", "euclidean"]: + # Begin by constructing the adjacency matrix and adjacency graph + # for the system. Each node corresponds to an atom, such that each + # pair of nodes can be assigned an edge weighted by the square of + # their distance in real space. + for adx in range(natm): + for bdx in range(adx + 1, natm): + dr = ( + euclidean_norm( + adx_map[adx]["coord"], + adx_map[bdx]["coord"], + ) + ** 2 + ) + fragment_map.adjacency_mat[adx, bdx] = dr + fragment_map.adjacency_graph.add_edge(adx, bdx, weight=dr) + + # For a given center site (adx), find the set of shortest + # paths to all other sites. The number of nodes visited + # on that path gives the degree of separation of the + # sites. + for adx, map in adx_map.items(): + fragment_map.center_atoms.append(tuple()) + fsites_temp = fragment_map.sites[adx] + fs_temp = [] + fs_temp.append(fragment_map.sites[adx]) + map["shortest_paths"] = dict( + single_source_all_shortest_paths( + fragment_map.adjacency_graph, + source=adx, + weight=lambda a, b, _: ( + fragment_map.adjacency_graph[a][b]["weight"] + ), + method="dijkstra", + ) + ) + + # If the degree of separation is smaller than the *n* + # in your fragment type, BE*n*, then that site is appended to + # the set of fragment sites for adx. + for bdx, path in map["shortest_paths"].items(): + if 0 < (len(path[0]) - 1) < fragment_type_order: + fsites_temp = tuple(fsites_temp + fragment_map.sites[bdx]) + fs_temp.append(tuple(fragment_map.sites[bdx])) + + fragment_map.fsites.append(tuple(fsites_temp)) + fragment_map.fs.append(tuple(fs_temp)) + fragment_map.center.append(tuple(fragment_map.sites[adx])) + + elif connectivity.lower() in ["resistance_distance", "resistance"]: + raise NotImplementedError("Work in progress...") + + elif connectivity.lower() in ["entanglement"]: + raise NotImplementedError("Work in progress...") + + else: + raise AttributeError(f"Connectivity metric not recognized: '{connectivity}'") + + if remove_nonunique_frags: + fragment_map.remove_nonnunique_frags(natm) + + # Define the 'edges' for fragment A as the intersect of its sites + # with the set of all center sites outside of A: + for adx, fs in enumerate(fragment_map.fs): + edge: set[tuple] = set() + for bdx, center in enumerate(fragment_map.center): + if adx == bdx: + pass + else: + overlap = set(fs).intersection(set((center,))) + if overlap: + edge = edge.union(overlap) + fragment_map.edge.append(tuple(edge)) + + # Update relative center site indices (centerf_idx) and weights + # for center site contributions to the energy (ebe_weights): + for adx, center in enumerate(fragment_map.center): + centerf_idx = tuple( + set([fragment_map.fsites[adx].index(cdx) for cdx in center]) + ) + ebe_weight = (1.0, tuple(centerf_idx)) + fragment_map.centerf_idx.append(centerf_idx) + fragment_map.ebe_weights.append(ebe_weight) + + # Finally, set fragment data names for scratch and bookkeeping: + for adx, _ in enumerate(fragment_map.fs): + fragment_map.dnames.append(str(frag_prefix) + str(adx)) + + return fragment_map + + def autogen( mol, frozen_core=True, be_type="be2", write_geom=False, - valence_basis=None, + iao_valence_basis=None, valence_only=False, print_frags=True, ): @@ -41,7 +321,7 @@ def autogen( write_geom : bool, optional Whether to write a 'fragment.xyz' file which contains all the fragments in Cartesian coordinates. Defaults to False. - valence_basis : str, optional + iao_valence_basis : str, optional Name of minimal basis set for IAO scheme. 'sto-3g' is sufficient for most cases. Defaults to None. valence_only : bool, optional @@ -84,7 +364,7 @@ def autogen( cell = mol.copy() else: cell = mol.copy() - cell.basis = valence_basis + cell.basis = iao_valence_basis cell.build() ncore, no_core_idx, core_list = get_core(cell) @@ -278,11 +558,11 @@ def autogen( w.close() # Prepare for PAO basis if requested - pao = bool(valence_basis and not valence_only) + pao = bool(iao_valence_basis and not valence_only) if pao: cell2 = cell.copy() - cell2.basis = valence_basis + cell2.basis = iao_valence_basis cell2.build() bas2list = cell2.aoslice_by_atom() diff --git a/src/quemb/molbe/fragment.py b/src/quemb/molbe/fragment.py index 36f1441a..c6b05c97 100644 --- a/src/quemb/molbe/fragment.py +++ b/src/quemb/molbe/fragment.py @@ -1,7 +1,7 @@ # Author: Oinam Romesh Meitei -from quemb.molbe.autofrag import autogen +from quemb.molbe.autofrag import autogen, graphgen from quemb.molbe.helper import get_core from quemb.molbe.lchain import chain as _ext_chain from quemb.shared.helper import copy_docstring @@ -31,7 +31,7 @@ class fragpart: mol : pyscf.gto.mole.Mole This is required for the options, 'autogen' and 'chain' as frag_type. - valence_basis: str + iao_valence_basis: str Name of minimal basis set for IAO scheme. 'sto-3g' suffice for most cases. valence_only: bool If this option is set to True, all calculation will be performed in @@ -50,7 +50,7 @@ def __init__( self, frag_type="autogen", closed=False, - valence_basis=None, + iao_valence_basis=None, valence_only=False, print_frags=True, write_geom=False, @@ -71,7 +71,7 @@ def __init__( self.centerf_idx = [] self.be_type = be_type self.frozen_core = frozen_core - self.valence_basis = valence_basis + self.iao_valence_basis = iao_valence_basis self.valence_only = valence_only # Initialize class attributes necessary for mixed-basis BE @@ -82,30 +82,45 @@ def __init__( # Check for frozen core approximation if frozen_core: - self.ncore, self.no_core_idx, self.core_list = get_core(mol) + self.ncore, self.no_core_idx, self.core_list = get_core(self.mol) + + if frag_type != "hchain_simple" and self.mol is None: + raise ValueError("Provide pyscf gto.M object in fragpart() and restart!") # Check type of fragmentation function if frag_type == "hchain_simple": # This is an experimental feature. self.hchain_simple() + elif frag_type == "chain": - if mol is None: - raise ValueError( - "Provide pyscf gto.M object in fragpart() and restart!" - ) - self.chain(mol, frozen_core=frozen_core, closed=closed) - elif frag_type == "autogen": - if mol is None: - raise ValueError( - "Provide pyscf gto.M object in fragpart() and restart!" - ) + self.chain(self.mol, frozen_core=frozen_core, closed=closed) + elif frag_type == "graphgen": + fragment_map = graphgen( + mol=self.mol.copy(), + be_type=be_type, + frozen_core=frozen_core, + remove_nonunique_frags=True, + frag_prefix="f", + connectivity="euclidean", + iao_valence_basis=iao_valence_basis, + ) + + self.fsites = fragment_map.fsites + self.edge = fragment_map.edge + self.center = fragment_map.center + # self.edge_idx = fragment_map["edge"] + self.centerf_idx = fragment_map.centerf_idx + self.ebe_weight = fragment_map.ebe_weights + self.Nfrag = len(self.fsites) + + elif frag_type == "autogen": fgs = autogen( mol, be_type=be_type, frozen_core=frozen_core, write_geom=write_geom, - valence_basis=valence_basis, + iao_valence_basis=iao_valence_basis, valence_only=valence_only, print_frags=print_frags, ) @@ -124,6 +139,7 @@ def __init__( self.add_center_atom, ) = fgs self.Nfrag = len(self.fsites) + else: raise ValueError(f"Fragmentation type = {frag_type} not implemented!") diff --git a/src/quemb/molbe/lo.py b/src/quemb/molbe/lo.py index 6ab763e0..d6dbe717 100644 --- a/src/quemb/molbe/lo.py +++ b/src/quemb/molbe/lo.py @@ -177,7 +177,9 @@ def get_pao(Ciao: Matrix, S: Matrix, S12: Matrix) -> Matrix: return cano_orth(Cpao_redundant, ovlp=S) -def get_pao_native(Ciao: Matrix, S: Matrix, mol: Mole, valence_basis: str) -> Matrix: +def get_pao_native( + Ciao: Matrix, S: Matrix, mol: Mole, iao_valence_basis: str +) -> Matrix: """ Parameters @@ -188,7 +190,7 @@ def get_pao_native(Ciao: Matrix, S: Matrix, mol: Mole, valence_basis: str) -> Ma ao ovlp matrix mol: mol object - valence_basis: + iao_valence_basis: basis used for valence orbitals Returns ------- @@ -200,7 +202,7 @@ def get_pao_native(Ciao: Matrix, S: Matrix, mol: Mole, valence_basis: str) -> Ma # Form a mol object with the valence basis for the ao_labels mol_alt = mol.copy() - mol_alt.basis = valence_basis + mol_alt.basis = iao_valence_basis mol_alt.build() full_ao_labels = mol.ao_labels() @@ -256,7 +258,7 @@ class MixinLocalize: def localize( self, lo_method, - valence_basis="sto-3g", + iao_valence_basis="sto-3g", hstack=False, pop_method=None, init_guess=None, @@ -273,7 +275,7 @@ def localize( lo_method : str Localization method in quantum chemistry. 'lowdin', 'boys', and 'iao' are supported. - valence_basis : str + iao_valence_basis : str Name of minimal basis set for IAO scheme. 'sto-3g' suffice for most cases. valence_only : bool If this option is set to True, all calculation will be performed in the @@ -360,7 +362,7 @@ def localize( W_ = C_ @ W_ self.W = get_loc( - self.mol, W_, "PM", pop_method=pop_method, init_guess=init_guess + self.mf.mol, W_, "PM", pop_method=pop_method, init_guess=init_guess ) if not self.frozen_core: @@ -374,7 +376,7 @@ def localize( # Occupied mo_coeff (with core) Co = self.C[:, : self.Nocc] # Get necessary overlaps, second arg is IAO basis - S12, S2 = get_xovlp(self.mol, basis=valence_basis) + S12, S2 = get_xovlp(self.mol, basis=iao_valence_basis) # Use these to get IAOs Ciao = get_iao(Co, S12, self.S, S2=S2) @@ -384,7 +386,10 @@ def localize( Cpao = get_pao(Ciao, self.S, S12) elif loc_type.upper() == "SO": Cpao = get_pao_native( - Ciao, self.S, self.mol, valence_basis=valence_basis + Ciao, + self.S, + self.mol, + iao_valence_basis=iao_valence_basis, ) # rearrange by atom diff --git a/src/quemb/molbe/mbe.py b/src/quemb/molbe/mbe.py index 73c77b30..520d3ba6 100644 --- a/src/quemb/molbe/mbe.py +++ b/src/quemb/molbe/mbe.py @@ -211,7 +211,7 @@ def __init__( self.localize( lo_method, pop_method=pop_method, - valence_basis=fobj.valence_basis, + iao_valence_basis=fobj.iao_valence_basis, valence_only=fobj.valence_only, ) @@ -219,7 +219,7 @@ def __init__( self.Ciao_pao = self.localize( lo_method, pop_method=pop_method, - valence_basis=fobj.valence_basis, + iao_valence_basis=fobj.iao_valence_basis, hstack=True, valence_only=False, nosave=True, diff --git a/src/quemb/molbe/solver.py b/src/quemb/molbe/solver.py index 636675c2..8462c4a0 100644 --- a/src/quemb/molbe/solver.py +++ b/src/quemb/molbe/solver.py @@ -908,8 +908,8 @@ def solve_block2( mc.fcisolver.twodot_to_onedot = DMRG_args.twodot_to_onedot mc.fcisolver.maxIter = DMRG_args.max_iter mc.fcisolver.block_extra_keyword = DMRG_args.block_extra_keyword - mc.fcisolver.scratchDirectory = str(frag_scratch) - mc.fcisolver.runtimeDir = str(frag_scratch) + mc.fcisolver.scratchDirectory = frag_scratch.path + mc.fcisolver.runtimeDir = frag_scratch.path mc.fcisolver.memory = DMRG_args.max_mem os.chdir(frag_scratch) diff --git a/src/quemb/molbe/ube.py b/src/quemb/molbe/ube.py index 83d24f10..5fd3153b 100644 --- a/src/quemb/molbe/ube.py +++ b/src/quemb/molbe/ube.py @@ -141,7 +141,7 @@ def __init__( self.localize( lo_method, - valence_basis=fobj.valence_basis, + iao_valence_basis=fobj.iao_valence_basis, valence_only=fobj.valence_only, ) diff --git a/tests/fragmentation_test.py b/tests/fragmentation_test.py new file mode 100644 index 00000000..a5f31962 --- /dev/null +++ b/tests/fragmentation_test.py @@ -0,0 +1,1419 @@ +""" +Tests for the fragmentation modules. + +Author(s): Shaun Weatherly +""" + +import os +import unittest + +from pyscf import gto, scf + +from quemb.molbe import BE, fragpart + + +class TestBE_Fragmentation(unittest.TestCase): + def test_autogen_h_linear_be1(self): + mol = gto.M() + mol.atom = [["H", (0.0, 0.0, i)] for i in range(8)] + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + + target = { + "fsites": [[0], [1], [2], [3], [4], [5], [6], [7]], + "edge": [], + "center": [], + "centerf_idx": [[0], [0], [0], [0], [0], [0], [0], [0]], + "ebe_weight": [ + [1.0, [0]], + [1.0, [0]], + [1.0, [0]], + [1.0, [0]], + [1.0, [0]], + [1.0, [0]], + [1.0, [0]], + [1.0, [0]], + ], + } + + self.run_indices_test( + mf, + "be1", + "autogen_h_linear_be1", + "autogen", + target, + ) + + def test_autogen_h_linear_be2(self): + mol = gto.M() + mol.atom = [["H", (0.0, 0.0, i)] for i in range(8)] + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + + target = { + "fsites": [ + [1, 0, 2], + [2, 1, 3], + [3, 2, 4], + [4, 3, 5], + [5, 4, 6], + [6, 7, 5], + ], + "edge": [[[2]], [[1], [3]], [[2], [4]], [[3], [5]], [[4], [6]], [[5]]], + "center": [[1], [0, 2], [1, 3], [2, 4], [3, 5], [4]], + "centerf_idx": [[0], [0], [0], [0], [0], [0]], + "ebe_weight": [ + [1.0, [0, 1]], + [1.0, [0]], + [1.0, [0]], + [1.0, [0]], + [1.0, [0]], + [1.0, [0, 1]], + ], + } + + self.run_indices_test( + mf, + "be2", + "autogen_h_linear_be2", + "autogen", + target, + ) + + def test_autogen_h_linear_be3(self): + mol = gto.M() + mol.atom = [["H", (0.0, 0.0, i)] for i in range(8)] + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + + target = { + "fsites": [ + [2, 0, 1, 3, 4], + [3, 2, 1, 4, 5], + [4, 3, 2, 5, 6], + [5, 6, 7, 4, 3], + ], + "edge": [ + [[3], [4]], + [[2], [1], [4], [5]], + [[3], [2], [5], [6]], + [[4], [3]], + ], + "center": [[1, 2], [0, 0, 2, 3], [1, 0, 3, 3], [2, 1]], + "centerf_idx": [[0], [0], [0], [0]], + "ebe_weight": [[1.0, [0, 1, 2]], [1.0, [0]], [1.0, [0]], [1.0, [0, 1, 2]]], + } + + self.run_indices_test( + mf, + "be3", + "autogen_h_linear_be3", + "autogen", + target, + ) + + def test_autogen_octane_be1(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/octane.xyz") + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + + target = { + "fsites": [ + [0, 1, 2, 3, 4, 11, 13], + [5, 6, 7, 8, 9, 10, 12], + [14, 15, 16, 17, 18, 24, 26], + [19, 20, 21, 22, 23, 25, 27], + [28, 29, 30, 31, 32, 38, 40], + [33, 34, 35, 36, 37, 39, 41], + [42, 43, 44, 45, 46, 52, 54, 57], + [47, 48, 49, 50, 51, 53, 55, 56], + ], + "edge": [], + "center": [], + "centerf_idx": [ + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6, 7], + [0, 1, 2, 3, 4, 5, 6, 7], + ], + "ebe_weight": [ + [1.0, [0, 1, 2, 3, 4, 5, 6]], + [1.0, [0, 1, 2, 3, 4, 5, 6]], + [1.0, [0, 1, 2, 3, 4, 5, 6]], + [1.0, [0, 1, 2, 3, 4, 5, 6]], + [1.0, [0, 1, 2, 3, 4, 5, 6]], + [1.0, [0, 1, 2, 3, 4, 5, 6]], + [1.0, [0, 1, 2, 3, 4, 5, 6, 7]], + [1.0, [0, 1, 2, 3, 4, 5, 6, 7]], + ], + } + + self.run_indices_test( + mf, + "be1", + "autogen_octane_be1", + "autogen", + target, + ) + + def test_autogen_octane_be2(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/octane.xyz") + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + + target = { + "fsites": [ + [ + 0, + 1, + 2, + 3, + 4, + 11, + 13, + 5, + 6, + 7, + 8, + 9, + 10, + 12, + 19, + 20, + 21, + 22, + 23, + 25, + 27, + ], + [ + 5, + 6, + 7, + 8, + 9, + 10, + 12, + 0, + 1, + 2, + 3, + 4, + 11, + 13, + 14, + 15, + 16, + 17, + 18, + 24, + 26, + ], + [ + 14, + 15, + 16, + 17, + 18, + 24, + 26, + 5, + 6, + 7, + 8, + 9, + 10, + 12, + 28, + 29, + 30, + 31, + 32, + 38, + 40, + ], + [ + 19, + 20, + 21, + 22, + 23, + 25, + 27, + 0, + 1, + 2, + 3, + 4, + 11, + 13, + 33, + 34, + 35, + 36, + 37, + 39, + 41, + ], + [ + 28, + 29, + 30, + 31, + 32, + 38, + 40, + 42, + 43, + 44, + 45, + 46, + 52, + 54, + 57, + 14, + 15, + 16, + 17, + 18, + 24, + 26, + ], + [ + 33, + 34, + 35, + 36, + 37, + 39, + 41, + 47, + 48, + 49, + 50, + 51, + 53, + 55, + 56, + 19, + 20, + 21, + 22, + 23, + 25, + 27, + ], + ], + "edge": [ + [[5, 6, 7, 8, 9, 10, 12], [19, 20, 21, 22, 23, 25, 27]], + [[0, 1, 2, 3, 4, 11, 13], [14, 15, 16, 17, 18, 24, 26]], + [[5, 6, 7, 8, 9, 10, 12], [28, 29, 30, 31, 32, 38, 40]], + [[0, 1, 2, 3, 4, 11, 13], [33, 34, 35, 36, 37, 39, 41]], + [[14, 15, 16, 17, 18, 24, 26]], + [[19, 20, 21, 22, 23, 25, 27]], + ], + "center": [[1, 3], [0, 2], [1, 4], [0, 5], [2], [3]], + "centerf_idx": [ + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + ], + "ebe_weight": [ + [1.0, [0, 1, 2, 3, 4, 5, 6]], + [1.0, [0, 1, 2, 3, 4, 5, 6]], + [1.0, [0, 1, 2, 3, 4, 5, 6]], + [1.0, [0, 1, 2, 3, 4, 5, 6]], + [1.0, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]], + [1.0, [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]], + ], + } + + self.run_indices_test( + mf, + "be2", + "autogen_octane_be2", + "autogen", + target, + ) + + def test_autogen_octane_be3(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/octane.xyz") + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + + target = { + "fsites": [ + [ + 0, + 1, + 2, + 3, + 4, + 11, + 13, + 5, + 6, + 7, + 8, + 9, + 10, + 12, + 14, + 15, + 16, + 17, + 18, + 24, + 26, + 19, + 20, + 21, + 22, + 23, + 25, + 27, + 33, + 34, + 35, + 36, + 37, + 39, + 41, + ], + [ + 5, + 6, + 7, + 8, + 9, + 10, + 12, + 0, + 1, + 2, + 3, + 4, + 11, + 13, + 19, + 20, + 21, + 22, + 23, + 25, + 27, + 14, + 15, + 16, + 17, + 18, + 24, + 26, + 28, + 29, + 30, + 31, + 32, + 38, + 40, + ], + [ + 14, + 15, + 16, + 17, + 18, + 24, + 26, + 28, + 29, + 30, + 31, + 32, + 38, + 40, + 42, + 43, + 44, + 45, + 46, + 52, + 54, + 57, + 5, + 6, + 7, + 8, + 9, + 10, + 12, + 0, + 1, + 2, + 3, + 4, + 11, + 13, + ], + [ + 19, + 20, + 21, + 22, + 23, + 25, + 27, + 33, + 34, + 35, + 36, + 37, + 39, + 41, + 47, + 48, + 49, + 50, + 51, + 53, + 55, + 56, + 0, + 1, + 2, + 3, + 4, + 11, + 13, + 5, + 6, + 7, + 8, + 9, + 10, + 12, + ], + ], + "edge": [ + [ + [5, 6, 7, 8, 9, 10, 12], + [14, 15, 16, 17, 18, 24, 26], + [19, 20, 21, 22, 23, 25, 27], + [33, 34, 35, 36, 37, 39, 41], + ], + [ + [0, 1, 2, 3, 4, 11, 13], + [19, 20, 21, 22, 23, 25, 27], + [14, 15, 16, 17, 18, 24, 26], + [28, 29, 30, 31, 32, 38, 40], + ], + [[5, 6, 7, 8, 9, 10, 12], [0, 1, 2, 3, 4, 11, 13]], + [[0, 1, 2, 3, 4, 11, 13], [5, 6, 7, 8, 9, 10, 12]], + ], + "center": [[1, 2, 3, 3], [0, 3, 2, 2], [1, 0], [0, 1]], + "centerf_idx": [ + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + [0, 1, 2, 3, 4, 5, 6], + ], + "ebe_weight": [ + [1.0, [0, 1, 2, 3, 4, 5, 6]], + [1.0, [0, 1, 2, 3, 4, 5, 6]], + [ + 1.0, + [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20, + 21, + ], + ], + [ + 1.0, + [ + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20, + 21, + ], + ], + ], + } + + self.run_indices_test( + mf, + "be3", + "autogen_octane_be3", + "autogen", + target, + ) + + def test_graphgen_h_linear_be1(self): + mol = gto.M() + mol.atom = [["H", (0.0, 0.0, i)] for i in range(8)] + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + + target = { + "fsites": [(0,), (1,), (2,), (3,), (4,), (5,), (6,), (7,)], + "edge": [(), (), (), (), (), (), (), ()], + "center": [(0,), (1,), (2,), (3,), (4,), (5,), (6,), (7,)], + "centerf_idx": [(0,), (0,), (0,), (0,), (0,), (0,), (0,), (0,)], + "ebe_weight": [ + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + ], + } + + self.run_indices_test( + mf, + "be1", + "graphgen_h_linear_be1", + "graphgen", + target, + ) + + def test_graphgen_h_linear_be2(self): + mol = gto.M() + mol.atom = [["H", (0.0, 0.0, i)] for i in range(8)] + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + + target = { + "fsites": [ + (1, 0, 2), + (2, 1, 3), + (3, 2, 4), + (4, 3, 5), + (5, 4, 6), + (6, 5, 7), + ], + "edge": [((2,),), ((3,),), ((2,), (4,)), ((3,), (5,)), ((4,),), ((5,),)], + "center": [(0, 1), (2,), (3,), (4,), (5,), (6, 7)], + "centerf_idx": [(0, 1), (0,), (0,), (0,), (0,), (0, 2)], + "ebe_weight": [ + (1.0, (0, 1)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0, 2)), + ], + } + + self.run_indices_test( + mf, + "be2", + "graphgen_h_linear_be2", + "graphgen", + target, + ) + + def test_graphgen_h_linear_be3(self): + mol = gto.M() + mol.atom = [["H", (0.0, 0.0, i)] for i in range(8)] + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + + target = { + "fsites": [ + (2, 0, 1, 3, 4), + (3, 1, 2, 4, 5), + (4, 2, 3, 5, 6), + (5, 3, 4, 6, 7), + ], + "edge": [((3,), (4,)), ((4,),), ((3,),), ((3,), (4,))], + "center": [(0, 1, 2), (3,), (4,), (5, 6, 7)], + "centerf_idx": [(0, 1, 2), (0,), (0,), (0, 3, 4)], + "ebe_weight": [ + (1.0, (0, 1, 2)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0, 3, 4)), + ], + } + + self.run_indices_test( + mf, + "be3", + "graphgen_h_linear_be3", + "graphgen", + target, + ) + + def test_graphgen_octane_be1(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/octane.xyz") + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + + target = { + "fsites": [ + (0, 1, 2, 3, 4), + (5, 6, 7, 8, 9), + (10,), + (11,), + (12,), + (13,), + (14, 15, 16, 17, 18), + (19, 20, 21, 22, 23), + (24,), + (25,), + (26,), + (27,), + (28, 29, 30, 31, 32), + (33, 34, 35, 36, 37), + (38,), + (39,), + (40,), + (41,), + (42, 43, 44, 45, 46), + (47, 48, 49, 50, 51), + (52,), + (53,), + (54,), + (55,), + (56,), + (57,), + ], + "edge": [ + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + (), + ], + "center": [ + (0, 1, 2, 3, 4), + (5, 6, 7, 8, 9), + (10,), + (11,), + (12,), + (13,), + (14, 15, 16, 17, 18), + (19, 20, 21, 22, 23), + (24,), + (25,), + (26,), + (27,), + (28, 29, 30, 31, 32), + (33, 34, 35, 36, 37), + (38,), + (39,), + (40,), + (41,), + (42, 43, 44, 45, 46), + (47, 48, 49, 50, 51), + (52,), + (53,), + (54,), + (55,), + (56,), + (57,), + ], + "centerf_idx": [ + (0, 1, 2, 3, 4), + (0, 1, 2, 3, 4), + (0,), + (0,), + (0,), + (0,), + (0, 1, 2, 3, 4), + (0, 1, 2, 3, 4), + (0,), + (0,), + (0,), + (0,), + (0, 1, 2, 3, 4), + (0, 1, 2, 3, 4), + (0,), + (0,), + (0,), + (0,), + (0, 1, 2, 3, 4), + (0, 1, 2, 3, 4), + (0,), + (0,), + (0,), + (0,), + (0,), + (0,), + ], + "ebe_weight": [ + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + ], + } + + self.run_indices_test( + mf, + "be1", + "graphgen_octane_be1", + "graphgen", + target, + ) + + def test_graphgen_octane_be2(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/octane.xyz") + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + + target = { + "fsites": [ + (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 13, 19, 20, 21, 22, 23), + (5, 6, 7, 8, 9, 0, 1, 2, 3, 4, 10, 12, 14, 15, 16, 17, 18), + (10, 5, 6, 7, 8, 9, 27, 40), + (11, 0, 1, 2, 3, 4, 26, 41), + (12, 5, 6, 7, 8, 9, 25, 38), + (13, 0, 1, 2, 3, 4, 24, 39), + (14, 15, 16, 17, 18, 5, 6, 7, 8, 9, 24, 26, 28, 29, 30, 31, 32), + (19, 20, 21, 22, 23, 0, 1, 2, 3, 4, 25, 27, 33, 34, 35, 36, 37), + (24, 13, 14, 15, 16, 17, 18, 52), + (25, 12, 19, 20, 21, 22, 23, 53), + (26, 11, 14, 15, 16, 17, 18, 54), + (27, 10, 19, 20, 21, 22, 23, 55), + (28, 29, 30, 31, 32, 14, 15, 16, 17, 18, 38, 40, 42, 43, 44, 45, 46), + (33, 34, 35, 36, 37, 19, 20, 21, 22, 23, 39, 41, 47, 48, 49, 50, 51), + (38, 12, 28, 29, 30, 31, 32), + (39, 13, 33, 34, 35, 36, 37), + (40, 10, 28, 29, 30, 31, 32), + (41, 11, 33, 34, 35, 36, 37), + (42, 43, 44, 45, 46, 28, 29, 30, 31, 32, 52, 54, 57), + (47, 48, 49, 50, 51, 33, 34, 35, 36, 37, 53, 55, 56), + (52, 24, 42, 43, 44, 45, 46), + (53, 25, 47, 48, 49, 50, 51), + (54, 26, 42, 43, 44, 45, 46), + (55, 27, 47, 48, 49, 50, 51), + ], + "edge": [ + ((5, 6, 7, 8, 9), (13,), (19, 20, 21, 22, 23), (11,)), + ((12,), (0, 1, 2, 3, 4), (10,), (14, 15, 16, 17, 18)), + ((5, 6, 7, 8, 9), (40,), (27,)), + ((41,), (0, 1, 2, 3, 4), (26,)), + ((5, 6, 7, 8, 9), (25,), (38,)), + ((24,), (0, 1, 2, 3, 4), (39,)), + ((5, 6, 7, 8, 9), (28, 29, 30, 31, 32), (24,), (26,)), + ((0, 1, 2, 3, 4), (25,), (27,), (33, 34, 35, 36, 37)), + ((52,), (13,), (14, 15, 16, 17, 18)), + ((12,), (53,), (19, 20, 21, 22, 23)), + ((11,), (54,), (14, 15, 16, 17, 18)), + ((55,), (19, 20, 21, 22, 23), (10,)), + ((40,), (38,), (14, 15, 16, 17, 18)), + ((41,), (19, 20, 21, 22, 23), (39,)), + ((28, 29, 30, 31, 32), (12,)), + ((13,), (33, 34, 35, 36, 37)), + ((28, 29, 30, 31, 32), (10,)), + ((11,), (33, 34, 35, 36, 37)), + ((28, 29, 30, 31, 32), (52,), (54,)), + ((53,), (55,), (33, 34, 35, 36, 37)), + ((24,),), + ((25,),), + ((26,),), + ((27,),), + ], + "center": [ + (0, 1, 2, 3, 4), + (5, 6, 7, 8, 9), + (10,), + (11,), + (12,), + (13,), + (14, 15, 16, 17, 18), + (19, 20, 21, 22, 23), + (24,), + (25,), + (26,), + (27,), + (28, 29, 30, 31, 32), + (33, 34, 35, 36, 37), + (38,), + (39,), + (40,), + (41,), + (42, 43, 44, 45, 46, 57), + (47, 48, 49, 50, 51, 56), + (52,), + (53,), + (54,), + (55,), + ], + "centerf_idx": [ + (0, 1, 2, 3, 4), + (0, 1, 2, 3, 4), + (0,), + (0,), + (0,), + (0,), + (0, 1, 2, 3, 4), + (0, 1, 2, 3, 4), + (0,), + (0,), + (0,), + (0,), + (0, 1, 2, 3, 4), + (0, 1, 2, 3, 4), + (0,), + (0,), + (0,), + (0,), + (0, 1, 2, 3, 4, 12), + (0, 1, 2, 3, 4, 12), + (0,), + (0,), + (0,), + (0,), + ], + "ebe_weight": [ + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0, 1, 2, 3, 4, 12)), + (1.0, (0, 1, 2, 3, 4, 12)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + ], + } + + self.run_indices_test( + mf, + "be2", + "graphgen_octane_be2", + "graphgen", + target, + ) + + def test_graphgen_octane_be3(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/octane.xyz") + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + + target = { + "fsites": [ + ( + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20, + 21, + 22, + 23, + 25, + 27, + 33, + 34, + 35, + 36, + 37, + ), + ( + 5, + 6, + 7, + 8, + 9, + 0, + 1, + 2, + 3, + 4, + 10, + 11, + 12, + 13, + 14, + 15, + 16, + 17, + 18, + 19, + 20, + 21, + 22, + 23, + 24, + 26, + 28, + 29, + 30, + 31, + 32, + ), + (10, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 12, 14, 15, 16, 17, 18, 27, 40), + (11, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 13, 19, 20, 21, 22, 23, 26, 41), + (12, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 14, 15, 16, 17, 18, 25, 38), + (13, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 19, 20, 21, 22, 23, 24, 39), + ( + 14, + 15, + 16, + 17, + 18, + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 10, + 12, + 24, + 26, + 28, + 29, + 30, + 31, + 32, + 38, + 40, + 42, + 43, + 44, + 45, + 46, + ), + ( + 19, + 20, + 21, + 22, + 23, + 0, + 1, + 2, + 3, + 4, + 5, + 6, + 7, + 8, + 9, + 11, + 13, + 25, + 27, + 33, + 34, + 35, + 36, + 37, + 39, + 41, + 47, + 48, + 49, + 50, + 51, + ), + (24, 5, 6, 7, 8, 9, 13, 14, 15, 16, 17, 18, 26, 28, 29, 30, 31, 32, 52), + (25, 0, 1, 2, 3, 4, 12, 19, 20, 21, 22, 23, 27, 33, 34, 35, 36, 37, 53), + (26, 5, 6, 7, 8, 9, 11, 14, 15, 16, 17, 18, 24, 28, 29, 30, 31, 32, 54), + (27, 0, 1, 2, 3, 4, 10, 19, 20, 21, 22, 23, 25, 33, 34, 35, 36, 37, 55), + ( + 28, + 29, + 30, + 31, + 32, + 5, + 6, + 7, + 8, + 9, + 14, + 15, + 16, + 17, + 18, + 24, + 26, + 38, + 40, + 42, + 43, + 44, + 45, + 46, + 52, + 54, + 57, + ), + ( + 33, + 34, + 35, + 36, + 37, + 0, + 1, + 2, + 3, + 4, + 19, + 20, + 21, + 22, + 23, + 25, + 27, + 39, + 41, + 47, + 48, + 49, + 50, + 51, + 53, + 55, + 56, + ), + ], + "edge": [ + ((5, 6, 7, 8, 9), (12,), (10,), (11,), (13,), (25,), (27,)), + ((12,), (26,), (10,), (24,), (11,), (13,), (0, 1, 2, 3, 4)), + ((5, 6, 7, 8, 9), (12,), (0, 1, 2, 3, 4), (27,)), + ((5, 6, 7, 8, 9), (13,), (0, 1, 2, 3, 4), (26,)), + ((5, 6, 7, 8, 9), (0, 1, 2, 3, 4), (25,), (10,)), + ((5, 6, 7, 8, 9), (24,), (0, 1, 2, 3, 4), (11,)), + ((5, 6, 7, 8, 9), (12,), (26,), (10,), (24,), (0, 1, 2, 3, 4)), + ((5, 6, 7, 8, 9), (11,), (13,), (0, 1, 2, 3, 4), (25,), (27,)), + ((5, 6, 7, 8, 9), (13,), (26,)), + ((12,), (0, 1, 2, 3, 4), (27,)), + ((5, 6, 7, 8, 9), (24,), (11,)), + ((0, 1, 2, 3, 4), (25,), (10,)), + ((5, 6, 7, 8, 9), (24,), (26,)), + ((0, 1, 2, 3, 4), (25,), (27,)), + ], + "center": [ + (0, 1, 2, 3, 4), + (5, 6, 7, 8, 9), + (10,), + (11,), + (12,), + (13,), + (38, 40, 14, 15, 16, 17, 18), + (39, 41, 19, 20, 21, 22, 23), + (24,), + (25,), + (26,), + (27,), + (32, 42, 43, 44, 45, 46, 52, 54, 57, 28, 29, 30, 31), + (33, 34, 35, 36, 37, 47, 48, 49, 50, 51, 53, 55, 56), + ], + "centerf_idx": [ + (0, 1, 2, 3, 4), + (0, 1, 2, 3, 4), + (0,), + (0,), + (0,), + (0,), + (0, 1, 2, 3, 4, 24, 25), + (0, 1, 2, 3, 4, 24, 25), + (0,), + (0,), + (0,), + (0,), + (0, 1, 2, 3, 4, 19, 20, 21, 22, 23, 24, 25, 26), + (0, 1, 2, 3, 4, 19, 20, 21, 22, 23, 24, 25, 26), + ], + "ebe_weight": [ + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0, 1, 2, 3, 4)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0, 1, 2, 3, 4, 24, 25)), + (1.0, (0, 1, 2, 3, 4, 24, 25)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0,)), + (1.0, (0, 1, 2, 3, 4, 19, 20, 21, 22, 23, 24, 25, 26)), + (1.0, (0, 1, 2, 3, 4, 19, 20, 21, 22, 23, 24, 25, 26)), + ], + } + + self.run_indices_test( + mf, + "be3", + "graphgen_octane_be3", + "graphgen", + target, + ) + + def test_graphgen_autogen_h_linear_be2(self): + mol = gto.M() + mol.atom = [["H", (0.0, 0.0, i)] for i in range(8)] + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + mf.kernel() + target = -0.13198886164212092 + + self.run_energies_test( + mf, + "be2", + "energy_graphgen_autogen_h_linear_be2", + target, + delta=1e-2, + ) + + def test_graphgen_autogen_octane_be2(self): + mol = gto.M() + mol.atom = os.path.join(os.path.dirname(__file__), "xyz/octane.xyz") + mol.basis = "sto-3g" + mol.charge = 0.0 + mol.spin = 0.0 + mol.build() + + mf = scf.RHF(mol) + mf.kernel() + target = -0.5499456086311243 + + self.run_energies_test( + mf, + "be2", + "energy_graphgen_autogen_octane_be2", + target, + delta=1e-2, + ) + + def run_energies_test( + self, + mf, + be_type, + test_name, + target, + delta, + ): + Es = {"target": target} + for frag_type in ["autogen", "graphgen"]: + fobj = fragpart(frag_type=frag_type, be_type=be_type, mol=mf.mol) + mbe = BE(mf, fobj) + mbe.oneshot(solver="CCSD") + Es.update({frag_type: mbe.ebe_tot}) + + for frag_type_A, E_A in Es.items(): + for frag_type_B, E_B in Es.items(): + self.assertAlmostEqual( + float(E_A), + float(E_B), + msg=f"{test_name}: BE Correlation Energy (oneshot) for " + + frag_type_A + + " does not match" + + frag_type_B + + f" ({E_A} != {E_B}) ", + delta=delta, + ) + + def run_indices_test( + self, + mf, + be_type, + test_name, + frag_type, + target, + ): + fobj = fragpart(frag_type=frag_type, be_type=be_type, mol=mf.mol) + try: + assert fobj.fsites == target["fsites"] + assert fobj.edge == target["edge"] + assert fobj.center == target["center"] + assert fobj.centerf_idx == target["centerf_idx"] + assert fobj.ebe_weight == target["ebe_weight"] + except AssertionError as e: + print(f"Fragmentation test failed at {test_name}:") + raise e + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/static_analysis_requirements.txt b/tests/static_analysis_requirements.txt index 2fef4fb8..ff9812ef 100644 --- a/tests/static_analysis_requirements.txt +++ b/tests/static_analysis_requirements.txt @@ -5,4 +5,5 @@ pylint mypy scipy-stubs types-PyYAML +types-networkx pytest