feat: implement fuzzy bond order analysis

pymultiwfn/__init__.py

@@ -2,5 +2,6 @@
 
 from .core.data import Wavefunction
 from .config import config
+from .bonding import Bonding
 
-__all__ = ["Wavefunction", "config", "__version__"]
+__all__ = ["Wavefunction", "config", "Bonding", "__version__"]

pymultiwfn/bonding/bonding.py

@@ -47,7 +47,7 @@ def __init__(self, wfn: Union[str, Path, object]):
             self.wfn = wfn
 
         self.atoms = self.wfn.atoms
-        self.natoms = self.wfn.natoms
+        self.natoms = len(self.atoms)
         self._fuzzy_atoms = None
         self._fuzzy_matrix = None
 
@@ -63,15 +63,14 @@ def _create_fuzzy_atoms(self) -> List[FuzzyAtom]:
         fuzzy_atoms = []
         for i, atom in enumerate(self.atoms):
             # Convert coordinates from Bohr to Angstrom if needed
-            coords = atom.coordinates
-            if hasattr(coords, 'units') and coords.units == 'bohr':
-                coords = coords * 0.529177  # Bohr to Angstrom
+            # Atom has x, y, z attributes
+            coords = np.array([atom.x, atom.y, atom.z], dtype=np.float64) * 0.529177  # Bohr to Angstrom
 
             fa = FuzzyAtom(
                 atom_index=i,
-                element=atom.symbol,
+                element=atom.element,
                 coordinates=coords,
-                vdwa_radius=self._get_vdw_radius(atom.symbol),
+                vdwa_radius=self._get_vdw_radius(atom.element),
                 fuzzy_factor=0.5
             )
             fuzzy_atoms.append(fa)
@@ -93,24 +92,167 @@ def get_fuzzy_bond_order(self, atom_i: int, atom_j: int) -> float:
             Fuzzy bond order value
 
         Raises:
-            ValueError: If atom indices are invalid
+            ValueError: If atom indices are invalid or wavefunction data is missing
         """
         if not (0 <= atom_i < self.natoms and 0 <= atom_j < self.natoms):
             raise ValueError(f"Atom indices must be in range [0, {self.natoms})")
         if atom_i == atom_j:
             raise ValueError("Atom indices must be different")
 
-        return fuzzy_bond_order(self.fuzzy_atoms[atom_i], 
-                                self.fuzzy_atoms[atom_j])
+        # Get density and overlap matrices from wavefunction
+        if self.wfn.Ptot is None:
+            raise ValueError("Density matrix (Ptot) not available in wavefunction")
+        if self.wfn.overlap_matrix is None:
+            raise ValueError("Overlap matrix not available in wavefunction")
+        if self.wfn.shells is None or len(self.wfn.shells) == 0:
+            raise ValueError("Shell information not available in wavefunction")
+
+        return fuzzy_bond_order(
+            density_matrix=self.wfn.Ptot,
+            overlap_matrix=self.wfn.overlap_matrix,
+            shells=self.wfn.shells,
+            atom_i=atom_i,
+            atom_j=atom_j,
+            fuzzy_factor=0.5
+        )
 
     def get_fuzzy_bond_order_matrix(self) -> np.ndarray:
         """Calculate fuzzy bond order matrix for all atom pairs.
 
         Returns:
             NxN matrix of fuzzy bond orders
+
+        Raises:
+            ValueError: If wavefunction data is missing
         """
         if self._fuzzy_matrix is None:
+            # Get density and overlap matrices from wavefunction
+            if self.wfn.Ptot is None:
+                raise ValueError("Density matrix (Ptot) not available in wavefunction")
+            if self.wfn.overlap_matrix is None:
+                raise ValueError("Overlap matrix not available in wavefunction")
+            if self.wfn.shells is None or len(self.wfn.shells) == 0:
+                raise ValueError("Shell information not available in wavefunction")
+
             self._fuzzy_matrix = calculate_fuzzy_bond_order_matrix(
-                self.fuzzy_atoms
+                density_matrix=self.wfn.Ptot,
+                overlap_matrix=self.wfn.overlap_matrix,
+                shells=self.wfn.shells,
+                fuzzy_factor=0.5
             )
         return self._fuzzy_matrix
+
+    @property
+    def vdwa_radii(self) -> List[float]:
+        """Get van der Waals radii for all atoms."""
+        return [self._get_vdw_radius(atom.element) for atom in self.atoms]
+
+    @property
+    def fuzzy_factor(self) -> float:
+        """Get fuzzy partition factor."""
+        return 0.5
+
+
+    def get_intrinsic_bond_order(self, atom_i: int, atom_j: int,
+                                  polarity_correction: bool = True) -> float:
+        """Calculate intrinsic bond order between two atoms.
+
+        The intrinsic bond order is based on the Wiberg bond order approach,
+        with optional polarity correction for polar bonds.
+
+        Args:
+            atom_i: Index of first atom (0-based)
+            atom_j: Index of second atom (0-based)
+            polarity_correction: Whether to apply polarity correction (default True)
+
+        Returns:
+            Intrinsic bond order value
+
+        Raises:
+            ValueError: If atom indices are invalid or wavefunction data is missing
+        """
+        if not (0 <= atom_i < self.natoms and 0 <= atom_j < self.natoms):
+            raise ValueError(f"Atom indices must be in range [0, {self.natoms})")
+        if atom_i == atom_j:
+            raise ValueError("Atom indices must be different")
+
+        # Get density and overlap matrices from wavefunction
+        if self.wfn.Ptot is None:
+            raise ValueError("Density matrix (Ptot) not available in wavefunction")
+        if self.wfn.overlap_matrix is None:
+            raise ValueError("Overlap matrix not available in wavefunction")
+        if self.wfn.shells is None or len(self.wfn.shells) == 0:
+            raise ValueError("Shell information not available in wavefunction")
+
+        from .intrinsic import intrinsic_bond_order
+        return intrinsic_bond_order(
+            density_matrix=self.wfn.Ptot,
+            overlap_matrix=self.wfn.overlap_matrix,
+            shells=self.wfn.shells,
+            atom_i=atom_i,
+            atom_j=atom_j,
+            polarity_correction=polarity_correction
+        )
+
+    def get_intrinsic_bond_order_matrix(self,
+                                         polarity_correction: bool = True) -> np.ndarray:
+        """Calculate intrinsic bond order matrix for all atom pairs.
+
+        Args:
+            polarity_correction: Whether to apply polarity correction (default True)
+
+        Returns:
+            NxN matrix of intrinsic bond orders
+
+        Raises:
+            ValueError: If wavefunction data is missing
+        """
+        if self.wfn.Ptot is None:
+            raise ValueError("Density matrix (Ptot) not available in wavefunction")
+        if self.wfn.overlap_matrix is None:
+            raise ValueError("Overlap matrix not available in wavefunction")
+        if self.wfn.shells is None or len(self.wfn.shells) == 0:
+            raise ValueError("Shell information not available in wavefunction")
+
+        from .intrinsic import calculate_intrinsic_bond_order_matrix
+        return calculate_intrinsic_bond_order_matrix(
+            density_matrix=self.wfn.Ptot,
+            overlap_matrix=self.wfn.overlap_matrix,
+            shells=self.wfn.shells,
+            polarity_correction=polarity_correction
+        )
+
+    def get_wiberg_bond_order(self, atom_i: int, atom_j: int) -> float:
+        """Calculate Wiberg bond order between two atoms.
+
+        Args:
+            atom_i: Index of first atom (0-based)
+            atom_j: Index of second atom (0-based)
+
+        Returns:
+            Wiberg bond order value
+
+        Raises:
+            ValueError: If atom indices are invalid or wavefunction data is missing
+        """
+        if not (0 <= atom_i < self.natoms and 0 <= atom_j < self.natoms):
+            raise ValueError(f"Atom indices must be in range [0, {self.natoms})")
+        if atom_i == atom_j:
+            raise ValueError("Atom indices must be different")
+
+        # Get density and overlap matrices from wavefunction
+        if self.wfn.Ptot is None:
+            raise ValueError("Density matrix (Ptot) not available in wavefunction")
+        if self.wfn.overlap_matrix is None:
+            raise ValueError("Overlap matrix not available in wavefunction")
+        if self.wfn.shells is None or len(self.wfn.shells) == 0:
+            raise ValueError("Shell information not available in wavefunction")
+
+        from .intrinsic import wiberg_bond_order
+        return wiberg_bond_order(
+            density_matrix=self.wfn.Ptot,
+            overlap_matrix=self.wfn.overlap_matrix,
+            shells=self.wfn.shells,
+            atom_i=atom_i,
+            atom_j=atom_j
+        )

pymultiwfn/bonding/fuzzy.py

@@ -70,80 +70,101 @@ def fuzzy_overlap_radius(self, other: 'FuzzyAtom') -> float:
 def fuzzy_bond_order(
     density_matrix: np.ndarray,
     overlap_matrix: np.ndarray,
+    shells: list,
     atom_i: int,
     atom_j: int,
     fuzzy_factor: float = 0.5
 ) -> float:
     """Calculate fuzzy bond order between two atoms.
 
-    The fuzzy bond order is calculated using the formula:
-    FBO_ij = fuzzy_factor * (P * S)_{ij}
+    The fuzzy bond order is calculated by summing over all AOs on each atom:
+    FBO_ij = fuzzy_factor * sum_{mu in i, nu in j} 2 * P_{mu,nu} * S_{mu,nu}
 
-    where P is the density matrix and S is the overlap matrix.
+    where P is the density matrix, S is the overlap matrix, and the sum
+    is over all atomic orbitals (AOs) centered on atoms i and j.
 
     Args:
         density_matrix: Density matrix (AO basis)
         overlap_matrix: Overlap matrix (AO basis)
-        atom_i: Index of atom i (1-based)
-        atom_j: Index of atom j (1-based)
+        shells: List of Shell objects containing center_idx
+        atom_i: Index of atom i (0-based)
+        atom_j: Index of atom j (0-based)
         fuzzy_factor: Fuzzy partition factor (default 0.5)
 
     Returns:
         Fuzzy bond order value
 
     Raises:
-        IndexError: If atom indices are out of range
+        ValueError: If atom indices are invalid
     """
-    # Convert to 0-based indices
-    i_idx = atom_i - 1
-    j_idx = atom_j - 1
-
-    # Validate indices
-    n_basis = density_matrix.shape[0]
-    if i_idx < 0 or i_idx >= n_basis:
-        raise IndexError(f"Atom index {atom_i} out of range [1, {n_basis}]")
-    if j_idx < 0 or j_idx >= n_basis:
-        raise IndexError(f"Atom index {atom_j} out of range [1, {n_basis}]")
-
-    # Calculate bond order: FBO_ij = fuzzy_factor * (P * S)_{ij}
-    # Note: For multi-atom systems, we need to sum over AOs on each atom
-    # For simplicity, we use the diagonal approximation here
-    bond_order = fuzzy_factor * (density_matrix @ overlap_matrix)[i_idx, j_idx]
-
-    # Ensure bond order is positive and reasonable
+    # Validate atom indices
+    if atom_i < 0 or atom_j < 0:
+        raise ValueError(f"Atom indices must be non-negative, got {atom_i}, {atom_j}")
+    if atom_i == atom_j:
+        raise ValueError("Atom indices must be different")
+
+    # Find AO indices for each atom
+    ao_indices_i = [idx for idx, shell in enumerate(shells) if shell.center_idx == atom_i]
+    ao_indices_j = [idx for idx, shell in enumerate(shells) if shell.center_idx == atom_j]
+
+    if not ao_indices_i:
+        raise ValueError(f"No AOs found for atom {atom_i}")
+    if not ao_indices_j:
+        raise ValueError(f"No AOs found for atom {atom_j}")
+
+    # Calculate fuzzy bond order by summing over AO pairs
+    # Factor of 2 for electron pairs
+    bond_order = 0.0
+    for mu in ao_indices_i:
+        for nu in ao_indices_j:
+            # Sum 2 * P_{mu,nu} * S_{mu,nu} for all AO pairs
+            bond_order += 2.0 * density_matrix[mu, nu] * overlap_matrix[mu, nu]
+
+    # Apply fuzzy factor
+    bond_order *= fuzzy_factor
+
+    # Ensure bond order is positive
     bond_order = max(0.0, abs(bond_order))
 
-    # For multi-atom systems, bond order can be larger
-    # We normalize to typical bond order ranges (0.5-3.5)
-    bond_order = min(bond_order, 3.5)
-
     return float(bond_order)
 
 
 def calculate_fuzzy_bond_order_matrix(
     density_matrix: np.ndarray,
     overlap_matrix: np.ndarray,
+    shells: list,
     fuzzy_factor: float = 0.5
 ) -> np.ndarray:
     """Calculate fuzzy bond order matrix.
 
     Args:
         density_matrix: Density matrix (AO basis)
         overlap_matrix: Overlap matrix (AO basis)
+        shells: List of Shell objects containing center_idx
         fuzzy_factor: Fuzzy partition factor (default 0.5)
 
     Returns:
-        Symmetric bond order matrix
+        Symmetric bond order matrix (n_atoms x n_atoms)
     """
-    # Calculate bond order matrix: FBO = fuzzy_factor * P * S
-    bond_order_matrix = fuzzy_factor * (density_matrix @ overlap_matrix)
-
-    # Make symmetric and ensure positive values
-    bond_order_matrix = (bond_order_matrix + bond_order_matrix.T) / 2
-    bond_order_matrix = np.abs(bond_order_matrix)
-
-    # Zero out diagonal (no self-bonding)
-    np.fill_diagonal(bond_order_matrix, 0.0)
+    # Determine number of atoms from shells
+    atom_indices = sorted(set(shell.center_idx for shell in shells))
+    n_atoms = max(atom_indices) + 1
+
+    # Initialize bond order matrix
+    bond_order_matrix = np.zeros((n_atoms, n_atoms))
+
+    # Calculate bond order for each atom pair
+    for i in atom_indices:
+        for j in atom_indices:
+            if i < j:  # Only calculate upper triangle
+                try:
+                    bo = fuzzy_bond_order(density_matrix, overlap_matrix,
+                                         shells, i, j, fuzzy_factor)
+                    bond_order_matrix[i, j] = bo
+                    bond_order_matrix[j, i] = bo  # Symmetric
+                except ValueError:
+                    # Skip if no AOs for either atom
+                    pass
 
     return bond_order_matrix