From dfbf1235a2105427a6984c311416f1559ddcc24d Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Sat, 17 Jan 2026 16:35:47 +0000 Subject: [PATCH 1/3] Add FCHK connectivity (bonds) parsing support - Add 'bonds' to @document_load_one optional attributes - Add MxBond, NBond, IBond, RBond label patterns - Add _load_connectivity() to parse bond data - Add _dump_connectivity() to write bond data - Add 6 tests for connectivity loading/dumping Fixes #399 --- iodata/formats/fchk.py | 159 ++++++++++++++++++++++++++++++++++++++- iodata/test/test_fchk.py | 73 ++++++++++++++++++ 2 files changed, 231 insertions(+), 1 deletion(-) diff --git a/iodata/formats/fchk.py b/iodata/formats/fchk.py index 6249349d2..477ed50f5 100644 --- a/iodata/formats/fchk.py +++ b/iodata/formats/fchk.py @@ -76,7 +76,17 @@ "run_type", "title", ], - ["energy", "atfrozen", "atgradient", "athessian", "atmasses", "one_rdms", "extra", "moments"], + [ + "bonds", + "energy", + "atfrozen", + "atgradient", + "athessian", + "atmasses", + "one_rdms", + "extra", + "moments", + ], ) def load_one(lit: LineIterator) -> dict: """Do not edit this docstring. It will be overwritten.""" @@ -125,6 +135,10 @@ def load_one(lit: LineIterator) -> dict: "Cartesian Gradient", "Cartesian Force Constants", "MicOpt", + "MxBond", + "NBond", + "IBond", + "RBond", ], ) @@ -278,6 +292,11 @@ def load_one(lit: LineIterator) -> dict: if atcharges: result["atcharges"] = atcharges + # F) Load connectivity (bonds) + bonds = _load_connectivity(fchk, len(result["atnums"])) + if bonds is not None: + result["bonds"] = bonds + return result @@ -468,6 +487,139 @@ def _load_fchk_field(lit: LineIterator, label_patterns: list[str]) -> tuple[str, return label, value +def _load_connectivity(fchk: dict, natoms: int) -> NDArray[int] | None: + """Load connectivity (bond) information from FCHK data. + + Parameters + ---------- + fchk + Dictionary containing parsed FCHK fields. + natoms + Number of atoms in the molecule. + + Returns + ------- + bonds + An (nbond, 3) integer array with [atom1, atom2, bond_type] for each bond, + or None if connectivity information is not present or invalid. + Atom indices are 0-based. Bond types follow the convention in + ``iodata.periodic.bond2num`` (1=single, 2=double, 3=triple). + + Notes + ----- + The FCHK format stores connectivity using four fields: + - MxBond: Maximum number of bonds per atom (scalar) + - NBond: Number of bonds for each atom (array of length natoms) + - IBond: Bonded atom indices, 1-based (array of length natoms * mxbond) + - RBond: Bond orders (array of length natoms * mxbond) + + """ + # Check if MxBond is present and valid + mxbond = fchk.get("MxBond") + if mxbond is None or mxbond < 1: + return None + + # Get required arrays + nbond_arr = fchk.get("NBond") + ibond_arr = fchk.get("IBond") + + # NBond and IBond are required if MxBond is set + if nbond_arr is None or ibond_arr is None: + return None + + # Check array sizes match expected dimensions + if len(nbond_arr) != natoms: + return None + if len(ibond_arr) != natoms * mxbond: + return None + + # Get optional RBond array (bond orders) + rbond_arr = fchk.get("RBond") + if rbond_arr is not None and len(rbond_arr) != natoms * mxbond: + rbond_arr = None + + # Reshape IBond and RBond to (natoms, mxbond) + ibond = ibond_arr.reshape(natoms, mxbond) + rbond = rbond_arr.reshape(natoms, mxbond) if rbond_arr is not None else None + + # Build bond list, avoiding duplicates + # Each bond appears twice in FCHK (once for each atom in the pair) + # We only add bonds where atom_i < atom_j to avoid duplicates + bonds = [] + for iatom in range(natoms): + num_bonds = int(nbond_arr[iatom]) + for j in range(min(num_bonds, mxbond)): + partner = int(ibond[iatom, j]) - 1 # Convert from 1-based to 0-based + if partner < 0 or partner >= natoms: + continue + if partner > iatom: # Only add each bond once + bond_order = round(rbond[iatom, j]) if rbond is not None else 1 + # FCHK uses 0=no bond, 1=single, 2=double, 3=triple + # This matches IOData's bond type convention + if bond_order >= 1: + bonds.append([iatom, partner, bond_order]) + + if len(bonds) == 0: + return None + + return np.array(bonds, dtype=int) + + +def _dump_connectivity(bonds: NDArray[int], natom: int, f: TextIO): + """Dump connectivity (bond) information to FCHK file. + + Parameters + ---------- + bonds + An (nbond, 3) array with [atom1, atom2, bond_type] for each bond. + Atom indices are 0-based. + natom + Number of atoms in the molecule. + f + The file object to write to. + + """ + if len(bonds) == 0: + return + + # Calculate MxBond (maximum bonds per atom) + bond_counts = np.zeros(natom, dtype=int) + for iatom, jatom, _ in bonds: + bond_counts[iatom] += 1 + bond_counts[jatom] += 1 + mxbond = int(bond_counts.max()) + + if mxbond == 0: + return + + # Build NBond, IBond, RBond arrays + nbond = bond_counts + ibond = np.zeros((natom, mxbond), dtype=int) + rbond = np.zeros((natom, mxbond), dtype=float) + + # Track how many bonds we've added for each atom + bond_idx = np.zeros(natom, dtype=int) + + for iatom, jatom, bond_order in bonds: + # Add bond from iatom's perspective + idx_i = bond_idx[iatom] + ibond[iatom, idx_i] = jatom + 1 # Convert to 1-based + rbond[iatom, idx_i] = float(bond_order) + bond_idx[iatom] += 1 + + # Add bond from jatom's perspective + idx_j = bond_idx[jatom] + ibond[jatom, idx_j] = iatom + 1 # Convert to 1-based + rbond[jatom, idx_j] = float(bond_order) + bond_idx[jatom] += 1 + + # Write to file + _dump_integer_scalars("MxBond", mxbond, f) + _dump_integer_arrays("NBond", nbond, f) + _dump_integer_arrays("IBond", ibond.flatten(), f) + _dump_real_arrays("RBond", rbond.flatten(), f) + + def _load_dm(label: str, fchk: dict, result: dict, key: str): """Load a density matrix from the FCHK file if present. @@ -608,6 +760,7 @@ def prepare_dump(data: IOData, allow_changes: bool, filename: str) -> IOData: "atgradient", "athessian", "atmasses", + "bonds", "charge", "energy", "lot", @@ -771,6 +924,10 @@ def dump_one(f: TextIO, data: IOData): arr = data.athessian[np.tril_indices(data.athessian.shape[0])] _dump_real_arrays("Cartesian Force Constants", arr, f) + # write connectivity (bonds) + if data.bonds is not None and len(data.bonds) > 0: + _dump_connectivity(data.bonds, data.natom, f) + # write moments if (1, "c") in data.moments: _dump_real_arrays("Dipole Moment", data.moments[(1, "c")], f) diff --git a/iodata/test/test_fchk.py b/iodata/test/test_fchk.py index d912b5aae..4914ac7eb 100644 --- a/iodata/test/test_fchk.py +++ b/iodata/test/test_fchk.py @@ -742,3 +742,76 @@ def test_methanol_g16_scan(): -115.43621498, ] ) + + +def test_load_fchk_connectivity_h2o(): + """Test parsing connectivity from H2O FCHK file.""" + mol = load_fchk_helper("h2o_sto3g.fchk") + # H2O should have 2 bonds: O-H and O-H + assert mol.bonds is not None + assert len(mol.bonds) == 2 + # Check atom indices (0-based): O(0) bonded to H(1) and H(2) + assert_equal(mol.bonds[:, :2], [[0, 1], [0, 2]]) + # Check bond types: all single bonds (type 1) + assert_equal(mol.bonds[:, 2], [1, 1]) + + +def test_load_fchk_connectivity_peroxide(): + """Test parsing connectivity from peroxide FCHK file.""" + mol = load_fchk_helper("peroxide_opt.fchk") + # H2O2 (peroxide) has 3 bonds: O-O, O-H, O-H + assert mol.bonds is not None + assert len(mol.bonds) == 3 + # All should be single bonds + assert_equal(mol.bonds[:, 2], [1, 1, 1]) + + +def test_load_fchk_connectivity_li2(): + """Test parsing connectivity from Li2 FCHK file.""" + mol = load_fchk_helper("li2_g09_nbasis_indep.fchk") + # Li2 should have 1 bond: Li-Li + assert mol.bonds is not None + assert len(mol.bonds) == 1 + assert_equal(mol.bonds[0, :2], [0, 1]) + assert mol.bonds[0, 2] == 3 # Triple bond (as stored in the file) + + +def test_load_fchk_no_connectivity(): + """Test FCHK file with zero bonds (NBond all zeros).""" + # water_atcharges.fchk has MxBond=1 but NBond=[0,0,0] - no actual bonds + mol = load_fchk_helper("water_atcharges.fchk") + # Should return None when no actual bonds are present + assert mol.bonds is None + + +def test_dump_load_connectivity_roundtrip(tmpdir): + """Test that connectivity survives a dump/load roundtrip.""" + mol1 = load_fchk_helper("h2o_sto3g.fchk") + assert mol1.bonds is not None + original_bonds = mol1.bonds.copy() + + # Dump to temporary file + fn_tmp = os.path.join(tmpdir, "h2o_roundtrip.fchk") + dump_one(mol1, fn_tmp) + + # Load back + mol2 = load_one(fn_tmp) + + # Verify bonds match + assert mol2.bonds is not None + assert_equal(mol2.bonds, original_bonds) + + +def test_dump_load_connectivity_peroxide_roundtrip(tmpdir): + """Test connectivity roundtrip for peroxide.""" + mol1 = load_fchk_helper("peroxide_opt.fchk") + assert mol1.bonds is not None + + fn_tmp = os.path.join(tmpdir, "peroxide_roundtrip.fchk") + dump_one(mol1, fn_tmp) + mol2 = load_one(fn_tmp) + + assert mol2.bonds is not None + assert len(mol2.bonds) == len(mol1.bonds) + # Bond order should be preserved + assert_equal(mol1.bonds[:, 2], mol2.bonds[:, 2]) From d914929e996f88b418f3aa7e3c27a562e00d4558 Mon Sep 17 00:00:00 2001 From: Amrinder Singh Date: Sun, 18 Jan 2026 12:46:12 +0530 Subject: [PATCH 2/3] Update iodata/test/test_fchk.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- iodata/test/test_fchk.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/iodata/test/test_fchk.py b/iodata/test/test_fchk.py index 4914ac7eb..1411ae377 100644 --- a/iodata/test/test_fchk.py +++ b/iodata/test/test_fchk.py @@ -773,7 +773,8 @@ def test_load_fchk_connectivity_li2(): assert mol.bonds is not None assert len(mol.bonds) == 1 assert_equal(mol.bonds[0, :2], [0, 1]) - assert mol.bonds[0, 2] == 3 # Triple bond (as stored in the file) + # Bond order 3 as stored in the test file; this test verifies parsing, not chemical accuracy. + assert mol.bonds[0, 2] == 3 def test_load_fchk_no_connectivity(): From 656ef1abb3db0116c7654ffa249f00470c919f92 Mon Sep 17 00:00:00 2001 From: aamrindersingh Date: Sun, 18 Jan 2026 07:43:39 +0000 Subject: [PATCH 3/3] vectorize connectivity functions and improve error handling --- iodata/formats/fchk.py | 111 +++++++++++++++++++++++------------------ 1 file changed, 63 insertions(+), 48 deletions(-) diff --git a/iodata/formats/fchk.py b/iodata/formats/fchk.py index 477ed50f5..f1d333c46 100644 --- a/iodata/formats/fchk.py +++ b/iodata/formats/fchk.py @@ -293,7 +293,7 @@ def load_one(lit: LineIterator) -> dict: result["atcharges"] = atcharges # F) Load connectivity (bonds) - bonds = _load_connectivity(fchk, len(result["atnums"])) + bonds = _load_connectivity(fchk, len(result["atnums"]), lit) if bonds is not None: result["bonds"] = bonds @@ -487,7 +487,7 @@ def _load_fchk_field(lit: LineIterator, label_patterns: list[str]) -> tuple[str, return label, value -def _load_connectivity(fchk: dict, natoms: int) -> NDArray[int] | None: +def _load_connectivity(fchk: dict, natoms: int, lit: LineIterator) -> NDArray[int] | None: """Load connectivity (bond) information from FCHK data. Parameters @@ -496,12 +496,14 @@ def _load_connectivity(fchk: dict, natoms: int) -> NDArray[int] | None: Dictionary containing parsed FCHK fields. natoms Number of atoms in the molecule. + lit + The line iterator for error reporting. Returns ------- bonds An (nbond, 3) integer array with [atom1, atom2, bond_type] for each bond, - or None if connectivity information is not present or invalid. + or None if connectivity information is not present. Atom indices are 0-based. Bond types follow the convention in ``iodata.periodic.bond2num`` (1=single, 2=double, 3=triple). @@ -525,44 +527,47 @@ def _load_connectivity(fchk: dict, natoms: int) -> NDArray[int] | None: # NBond and IBond are required if MxBond is set if nbond_arr is None or ibond_arr is None: - return None + raise LoadError("MxBond is set but NBond or IBond sections are missing.", lit.filename) # Check array sizes match expected dimensions if len(nbond_arr) != natoms: - return None + raise LoadError( + f"NBond array size {len(nbond_arr)} does not match number of atoms {natoms}.", + lit.filename, + ) if len(ibond_arr) != natoms * mxbond: - return None + raise LoadError( + f"IBond array size {len(ibond_arr)} does not match expected {natoms * mxbond}.", + lit.filename, + ) - # Get optional RBond array (bond orders) rbond_arr = fchk.get("RBond") if rbond_arr is not None and len(rbond_arr) != natoms * mxbond: rbond_arr = None - # Reshape IBond and RBond to (natoms, mxbond) ibond = ibond_arr.reshape(natoms, mxbond) rbond = rbond_arr.reshape(natoms, mxbond) if rbond_arr is not None else None - # Build bond list, avoiding duplicates - # Each bond appears twice in FCHK (once for each atom in the pair) - # We only add bonds where atom_i < atom_j to avoid duplicates - bonds = [] - for iatom in range(natoms): - num_bonds = int(nbond_arr[iatom]) - for j in range(min(num_bonds, mxbond)): - partner = int(ibond[iatom, j]) - 1 # Convert from 1-based to 0-based - if partner < 0 or partner >= natoms: - continue - if partner > iatom: # Only add each bond once - bond_order = round(rbond[iatom, j]) if rbond is not None else 1 - # FCHK uses 0=no bond, 1=single, 2=double, 3=triple - # This matches IOData's bond type convention - if bond_order >= 1: - bonds.append([iatom, partner, bond_order]) + # Vectorized bond extraction + nb = np.asarray(nbond_arr, dtype=np.int64) + partner = np.asarray(ibond, dtype=np.int64) - 1 # 0-based indexing + i = np.arange(natoms)[:, None] + j = np.arange(mxbond)[None, :] - if len(bonds) == 0: + present = j < nb[:, None] + valid = present & (partner >= 0) & (partner < natoms) & (partner > i) + + if rbond is None: + bo = np.ones_like(partner, dtype=np.int64) + else: + bo = np.rint(np.asarray(rbond)).astype(np.int64) + valid &= bo >= 1 + + ii, jj = np.nonzero(valid) + if len(ii) == 0: return None - return np.array(bonds, dtype=int) + return np.stack((ii, partner[ii, jj], bo[ii, jj]), axis=1) def _dump_connectivity(bonds: NDArray[int], natom: int, f: TextIO): @@ -582,38 +587,48 @@ def _dump_connectivity(bonds: NDArray[int], natom: int, f: TextIO): if len(bonds) == 0: return - # Calculate MxBond (maximum bonds per atom) - bond_counts = np.zeros(natom, dtype=int) - for iatom, jatom, _ in bonds: - bond_counts[iatom] += 1 - bond_counts[jatom] += 1 - mxbond = int(bond_counts.max()) + atom1 = bonds[:, 0] + atom2 = bonds[:, 1] + bond_order = bonds[:, 2] + + # Calculate bond counts per atom + nbond = np.zeros(natom, dtype=int) + np.add.at(nbond, atom1, 1) + np.add.at(nbond, atom2, 1) + mxbond = int(nbond.max()) if mxbond == 0: return - # Build NBond, IBond, RBond arrays - nbond = bond_counts ibond = np.zeros((natom, mxbond), dtype=int) rbond = np.zeros((natom, mxbond), dtype=float) - # Track how many bonds we've added for each atom - bond_idx = np.zeros(natom, dtype=int) + # Compute per-atom bond indices + order1 = np.argsort(atom1, kind="stable") + sorted_atom1 = atom1[order1] + idx1 = np.zeros(len(atom1), dtype=int) + idx1[order1] = ( + np.arange(len(atom1)) + - np.searchsorted(sorted_atom1, sorted_atom1, side="left")[np.argsort(order1)] + ) - for iatom, jatom, bond_order in bonds: - # Add bond from iatom's perspective - idx_i = bond_idx[iatom] - ibond[iatom, idx_i] = jatom + 1 # Convert to 1-based - rbond[iatom, idx_i] = float(bond_order) - bond_idx[iatom] += 1 + order2 = np.argsort(atom2, kind="stable") + sorted_atom2 = atom2[order2] + idx2_base = np.zeros(len(atom2), dtype=int) + idx2_base[order2] = ( + np.arange(len(atom2)) + - np.searchsorted(sorted_atom2, sorted_atom2, side="left")[np.argsort(order2)] + ) - # Add bond from jatom's perspective - idx_j = bond_idx[jatom] - ibond[jatom, idx_j] = iatom + 1 # Convert to 1-based - rbond[jatom, idx_j] = float(bond_order) - bond_idx[jatom] += 1 + atom1_counts = np.zeros(natom, dtype=int) + np.add.at(atom1_counts, atom1, 1) + idx2 = idx2_base + atom1_counts[atom2] - # Write to file + # Fill arrays with 1-based indexing for FCHK format + ibond[atom1, idx1] = atom2 + 1 + rbond[atom1, idx1] = bond_order.astype(float) + ibond[atom2, idx2] = atom1 + 1 + rbond[atom2, idx2] = bond_order.astype(float) _dump_integer_scalars("MxBond", mxbond, f) _dump_integer_arrays("NBond", nbond, f) _dump_integer_arrays("IBond", ibond.flatten(), f)