From 04ca4da51b4a651bb6e37e29ed2b148dc4f8a3a4 Mon Sep 17 00:00:00 2001 From: Haoming Yan Date: Tue, 3 Mar 2026 07:25:22 +0800 Subject: [PATCH] fix: Implement Delocalization Index (DI) - Add delocalization.py module with 2-center and 3-center DI calculations - Implement aromaticity indices (AI, PDI, FLU) from DI matrices - Add bond classification based on DI values - Create comprehensive test suite with 35 passing tests - Update bonding __init__.py to export new DI classes and functions Fixes newtontech/PyMultiWFN#4 --- pymultiwfn/bonding/__init__.py | 20 + pymultiwfn/bonding/delocalization.py | 617 +++++++++++++++++++++ tests/bonding/test_delocalization_index.py | 551 ++++++++++++++++++ 3 files changed, 1188 insertions(+) create mode 100644 pymultiwfn/bonding/delocalization.py create mode 100644 tests/bonding/test_delocalization_index.py diff --git a/pymultiwfn/bonding/__init__.py b/pymultiwfn/bonding/__init__.py index e165f383..2d9a36b7 100644 --- a/pymultiwfn/bonding/__init__.py +++ b/pymultiwfn/bonding/__init__.py @@ -8,10 +8,30 @@ from .fuzzy import FuzzyAtom, fuzzy_bond_order, calculate_fuzzy_bond_order_matrix from .bonding import Bonding +from .delocalization import ( + DelocalizationIndex, + DelocalizationResult, + delocalization_index, + three_center_delocalization_index, + calculate_di_matrix, + classify_bond_from_di, + calculate_aromaticity_index, + calculate_pdi, + calculate_flu, +) __all__ = [ 'Bonding', 'FuzzyAtom', 'fuzzy_bond_order', 'calculate_fuzzy_bond_order_matrix', + 'DelocalizationIndex', + 'DelocalizationResult', + 'delocalization_index', + 'three_center_delocalization_index', + 'calculate_di_matrix', + 'classify_bond_from_di', + 'calculate_aromaticity_index', + 'calculate_pdi', + 'calculate_flu', ] diff --git a/pymultiwfn/bonding/delocalization.py b/pymultiwfn/bonding/delocalization.py new file mode 100644 index 00000000..bcdb8ba4 --- /dev/null +++ b/pymultiwfn/bonding/delocalization.py @@ -0,0 +1,617 @@ +"""Delocalization Index Implementation (Issue 4). + +This module implements electron pair delocalization index analysis for +investigating electron sharing and aromaticity. The delocalization index (DI) +measures the number of electron pairs shared between atoms, providing +insights into bond character and aromatic systems. + +References: +- Bader, R.F.W. & Stephens, M.E. (1975). J. Am. Chem. Soc. 97, 7391-7399. +- Angyan, J.G. et al. (2006). J. Chem. Sci. 118, 159-169. +- Matito, E. et al. (2005). J. Phys. Chem. A 109, 5600-5609. +""" + +import numpy as np +from typing import Union, Tuple, List, Optional, Dict +from dataclasses import dataclass +from pathlib import Path + + +@dataclass +class DelocalizationResult: + """Container for delocalization index results. + + Attributes: + atom_i: Index of first atom (0-based) + atom_j: Index of second atom (0-based) + di_value: Delocalization index value in electrons + bond_type: Classification of bond type (single, double, aromatic, etc.) + """ + atom_i: int + atom_j: int + di_value: float + bond_type: str = "unknown" + + def __repr__(self) -> str: + return f"DI({self.atom_i},{self.atom_j}) = {self.di_value:.4f} e" + + +def delocalization_index( + density_matrix: np.ndarray, + overlap_matrix: np.ndarray, + atom_i_indices: List[int], + atom_j_indices: List[int], +) -> float: + """Calculate the 2-center delocalization index between two atoms. + + The delocalization index δ(A,B) measures the number of electron pairs + shared between atoms A and B. It is calculated using: + + δ(A,B) = 2 * Σ_μ∈A Σ_ν∈B |P_μν * S_μν| + + where P is the density matrix and S is the overlap matrix. + + Args: + density_matrix: Total density matrix (AO basis), shape (n_ao, n_ao) + overlap_matrix: Overlap matrix (AO basis), shape (n_ao, n_ao) + atom_i_indices: List of basis function indices belonging to atom i + atom_j_indices: List of basis function indices belonging to atom j + + Returns: + Delocalization index value in electrons + + Raises: + ValueError: If matrices have incompatible shapes or indices are empty + """ + # Validate inputs + if density_matrix.shape != overlap_matrix.shape: + raise ValueError( + f"Density matrix shape {density_matrix.shape} must match " + f"overlap matrix shape {overlap_matrix.shape}" + ) + + if not atom_i_indices or not atom_j_indices: + return 0.0 + + n_ao = density_matrix.shape[0] + + # Validate indices + for idx in atom_i_indices + atom_j_indices: + if idx < 0 or idx >= n_ao: + raise ValueError(f"Basis function index {idx} out of range [0, {n_ao})") + + # Calculate DI using the formula: δ(A,B) = 2 * Σ_μ∈A Σ_ν∈B |P_μν * S_μν| + di_value = 0.0 + for mu in atom_i_indices: + for nu in atom_j_indices: + # Use absolute value to ensure positive DI + di_value += abs(density_matrix[mu, nu] * overlap_matrix[mu, nu]) + + # Multiply by 2 for electron pairs + di_value *= 2.0 + + return float(di_value) + + +def three_center_delocalization_index( + density_matrix: np.ndarray, + overlap_matrix: np.ndarray, + atom_a_indices: List[int], + atom_b_indices: List[int], + atom_c_indices: List[int], +) -> float: + """Calculate the 3-center delocalization index. + + The 3-center DI measures electron sharing among three atoms, + important for analyzing multi-center bonds and aromaticity. + + δ(A,B,C) = 2 * Σ_μ∈A Σ_ν∈B Σ_λ∈C |P_μνλ * S_μνλ| + + For practical implementation, we use the approximation: + δ(A,B,C) ≈ √[δ(A,B) * δ(A,C) * δ(B,C)] / 2 + + Args: + density_matrix: Total density matrix (AO basis) + overlap_matrix: Overlap matrix (AO basis) + atom_a_indices: List of basis function indices for atom A + atom_b_indices: List of basis function indices for atom B + atom_c_indices: List of basis function indices for atom C + + Returns: + 3-center delocalization index value in electrons + """ + # Calculate 2-center DIs for all pairs + di_ab = delocalization_index( + density_matrix, overlap_matrix, atom_a_indices, atom_b_indices + ) + di_ac = delocalization_index( + density_matrix, overlap_matrix, atom_a_indices, atom_c_indices + ) + di_bc = delocalization_index( + density_matrix, overlap_matrix, atom_b_indices, atom_c_indices + ) + + # Use geometric mean approximation for 3-center DI + # δ(A,B,C) ≈ √[δ(A,B) * δ(A,C) * δ(B,C)] / 2 + if di_ab > 0 and di_ac > 0 and di_bc > 0: + di_3c = np.sqrt(di_ab * di_ac * di_bc) / 2.0 + else: + di_3c = 0.0 + + return float(di_3c) + + +def calculate_di_matrix( + density_matrix: np.ndarray, + overlap_matrix: np.ndarray, + atom_basis_indices: Dict[int, List[int]], + natoms: int, +) -> np.ndarray: + """Calculate the full delocalization index matrix. + + Args: + density_matrix: Total density matrix (AO basis) + overlap_matrix: Overlap matrix (AO basis) + atom_basis_indices: Dictionary mapping atom index to basis function indices + natoms: Number of atoms in the system + + Returns: + Symmetric natoms x natoms matrix of delocalization indices + """ + di_matrix = np.zeros((natoms, natoms)) + + for i in range(natoms): + for j in range(i + 1, natoms): + indices_i = atom_basis_indices.get(i, []) + indices_j = atom_basis_indices.get(j, []) + + if indices_i and indices_j: + di_val = delocalization_index( + density_matrix, overlap_matrix, indices_i, indices_j + ) + di_matrix[i, j] = di_val + di_matrix[j, i] = di_val + + return di_matrix + + +def classify_bond_from_di(di_value: float) -> str: + """Classify bond type based on delocalization index value. + + Args: + di_value: Delocalization index in electrons + + Returns: + Bond type string: 'single', 'double', 'triple', 'aromatic', or 'weak' + """ + if di_value < 0.3: + return "weak" + elif di_value < 0.8: + return "single" + elif di_value < 1.3: + return "aromatic" + elif di_value < 1.8: + return "double" + elif di_value < 2.5: + return "triple" + else: + return "very_strong" + + +def calculate_aromaticity_index( + di_matrix: np.ndarray, + ring_atoms: List[int], +) -> float: + """Calculate aromaticity index from DI matrix for a ring system. + + The aromaticity index based on DI (also related to PDI - Para + Delocalization Index) measures electron delocalization in rings. + + For a ring with n atoms, we compute: + AI = (1/n) * Σ_ δ(i,j) / δ_ref + + where δ_ref is the reference DI for a perfect aromatic bond (~1.4 e). + + Args: + di_matrix: Full delocalization index matrix + ring_atoms: List of atom indices forming the ring (0-based) + + Returns: + Aromaticity index (1.0 = perfect aromatic, <1.0 = less aromatic) + """ + if len(ring_atoms) < 3: + return 0.0 + + # Reference DI for perfect aromatic bond (e.g., benzene) + di_ref = 1.4 # electrons + + # Calculate average DI for adjacent atoms in ring + n_atoms = len(ring_atoms) + total_di = 0.0 + + for i in range(n_atoms): + # Adjacent atoms in ring (with periodic boundary) + j = (i + 1) % n_atoms + atom_i = ring_atoms[i] + atom_j = ring_atoms[j] + total_di += di_matrix[atom_i, atom_j] + + # Average DI per bond + avg_di = total_di / n_atoms + + # Normalize to aromaticity index + aromaticity_index = avg_di / di_ref + + return float(aromaticity_index) + + +def calculate_pdi( + di_matrix: np.ndarray, + ring_atoms: List[int], +) -> float: + """Calculate Para-Delocalization Index (PDI) for a six-membered ring. + + PDI measures electron delocalization between para positions in + six-membered rings, useful for aromaticity analysis. + + PDI = (δ(1,4) + δ(2,5) + δ(3,6)) / 3 + + Reference: Matito, E. et al. (2005). J. Phys. Chem. A 109, 5600. + + Args: + di_matrix: Full delocalization index matrix + ring_atoms: List of 6 atom indices forming the ring (0-based) + + Returns: + PDI value in electrons + + Raises: + ValueError: If ring doesn't have exactly 6 atoms + """ + if len(ring_atoms) != 6: + raise ValueError(f"PDI requires exactly 6 ring atoms, got {len(ring_atoms)}") + + # Para positions: (1,4), (2,5), (3,6) - 0-indexed: (0,3), (1,4), (2,5) + pdi = ( + di_matrix[ring_atoms[0], ring_atoms[3]] + + di_matrix[ring_atoms[1], ring_atoms[4]] + + di_matrix[ring_atoms[2], ring_atoms[5]] + ) / 3.0 + + return float(pdi) + + +def calculate_flu( + di_matrix: np.ndarray, + ring_atoms: List[int], + natoms_total: int = None, +) -> float: + """Calculate Aromatic Fluctuation Index (FLU) for a ring. + + FLU measures the fluctuation of electron delocalization relative + to a reference aromatic system. Lower FLU indicates more aromatic. + + FLU = (1/n) * Σ_i [(δ_ref - δ_i) / δ_ref]^2 + + where δ_i is the sum of DIs for atom i with its neighbors in the ring. + + Reference: Matito, E. et al. (2006). Chem. Phys. Lett. 420, 287. + + Args: + di_matrix: Full delocalization index matrix + ring_atoms: List of atom indices forming the ring (0-based) + natoms_total: Total number of atoms (unused, kept for API compatibility) + + Returns: + FLU value (0 = perfectly aromatic, higher = less aromatic) + """ + if len(ring_atoms) < 3: + return 1.0 + + n_ring = len(ring_atoms) + + # Reference DI sum for aromatic atom (benzene reference ~2.8 e total) + # Each C in benzene has DI ~1.4 with each of 2 neighbors + di_ref_sum = 2.8 + + flu = 0.0 + for i in range(n_ring): + # Get DI with two neighbors in ring + j_prev = ring_atoms[(i - 1) % n_ring] + j_next = ring_atoms[(i + 1) % n_ring] + atom_i = ring_atoms[i] + + di_sum = di_matrix[atom_i, j_prev] + di_matrix[atom_i, j_next] + + # FLU contribution + if di_ref_sum > 0: + flu += ((di_ref_sum - di_sum) / di_ref_sum) ** 2 + + flu /= n_ring + + return float(flu) + + +class DelocalizationIndex: + """Main class for delocalization index analysis. + + Provides methods for calculating 2-center and 3-center delocalization + indices, aromaticity indices, and bond classification. + + Attributes: + wfn: Wavefunction object containing molecular data + natoms: Number of atoms + density_matrix: Total density matrix (AO basis) + overlap_matrix: Overlap matrix (AO basis) + atom_basis_indices: Mapping from atom index to basis function indices + + Example: + >>> from pymultiwfn.bonding import DelocalizationIndex + >>> di = DelocalizationIndex(wfn) + >>> di_val = di.get_delocalization_index(atom_i=0, atom_j=1) + >>> print(f"DI: {di_val:.3f} e") + """ + + def __init__(self, wfn: Union[str, Path, object]): + """Initialize DelocalizationIndex analysis. + + Args: + wfn: Path to wavefunction file or Wavefunction object + + Raises: + FileNotFoundError: If wavefunction file doesn't exist + ValueError: If wavefunction is invalid or missing required data + """ + # Import here to avoid circular imports + from ..io import load + + if isinstance(wfn, (str, Path)): + self.wfn = load(wfn) + else: + self.wfn = wfn + + self.natoms = self.wfn.num_atoms + self.atoms = self.wfn.atoms + + # Get density and overlap matrices + self._validate_wavefunction() + self.density_matrix = self._get_density_matrix() + self.overlap_matrix = self._get_overlap_matrix() + self.atom_basis_indices = self._get_atom_basis_indices() + + # Cached DI matrix + self._di_matrix = None + + def _validate_wavefunction(self) -> None: + """Validate that wavefunction has required data.""" + if self.wfn is None: + raise ValueError("Wavefunction object is None") + + if self.natoms == 0: + raise ValueError("Wavefunction has no atoms") + + def _get_density_matrix(self) -> np.ndarray: + """Get or calculate the total density matrix.""" + if self.wfn.Ptot is not None: + return self.wfn.Ptot + + # Calculate if not present + if hasattr(self.wfn, 'calculate_density_matrices'): + self.wfn.calculate_density_matrices() + if self.wfn.Ptot is not None: + return self.wfn.Ptot + + # Fallback: construct from coefficients and occupations + if self.wfn.coefficients is not None and self.wfn.occupations is not None: + # P = C * occ * C^T (for MOs as rows) + C = self.wfn.coefficients + occ = self.wfn.occupations + P = np.einsum('oi,oj->ij', C * occ[:, np.newaxis], C) + return P + + raise ValueError("Cannot construct density matrix from wavefunction") + + def _get_overlap_matrix(self) -> np.ndarray: + """Get or calculate the overlap matrix.""" + if self.wfn.overlap_matrix is not None: + return self.wfn.overlap_matrix + + # Calculate if method exists + if hasattr(self.wfn, 'calculate_overlap_matrix'): + self.wfn.calculate_overlap_matrix() + if self.wfn.overlap_matrix is not None: + return self.wfn.overlap_matrix + + # Fallback: use identity matrix (crude approximation) + n_basis = self.wfn.num_basis + return np.eye(n_basis) + + def _get_atom_basis_indices(self) -> Dict[int, List[int]]: + """Get mapping from atom index to basis function indices.""" + if hasattr(self.wfn, 'get_atomic_basis_indices'): + return self.wfn.get_atomic_basis_indices() + + # Fallback: assign basis functions equally (approximation) + n_basis = self.wfn.num_basis + indices = {} + basis_per_atom = max(1, n_basis // self.natoms) + + for i in range(self.natoms): + start = i * basis_per_atom + end = min((i + 1) * basis_per_atom, n_basis) + indices[i] = list(range(start, end)) + + # Handle remaining basis functions + remaining = n_basis - self.natoms * basis_per_atom + if remaining > 0: + for i in range(remaining): + atom_idx = i % self.natoms + bf_idx = self.natoms * basis_per_atom + i + indices[atom_idx].append(bf_idx) + + return indices + + def get_delocalization_index(self, atom_i: int, atom_j: int) -> float: + """Calculate 2-center delocalization index between two atoms. + + Args: + atom_i: Index of first atom (0-based) + atom_j: Index of second atom (0-based) + + Returns: + Delocalization index value in electrons + + Raises: + ValueError: If atom indices are invalid + """ + 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}), " + f"got {atom_i} and {atom_j}" + ) + + if atom_i == atom_j: + return 0.0 + + indices_i = self.atom_basis_indices.get(atom_i, []) + indices_j = self.atom_basis_indices.get(atom_j, []) + + return delocalization_index( + self.density_matrix, self.overlap_matrix, indices_i, indices_j + ) + + def get_three_center_di( + self, atom_a: int, atom_b: int, atom_c: int + ) -> float: + """Calculate 3-center delocalization index. + + Args: + atom_a: Index of first atom (0-based) + atom_b: Index of second atom (0-based) + atom_c: Index of third atom (0-based) + + Returns: + 3-center delocalization index value in electrons + + Raises: + ValueError: If atom indices are invalid or not distinct + """ + if len({atom_a, atom_b, atom_c}) != 3: + raise ValueError("Three center DI requires three distinct atoms") + + for atom in [atom_a, atom_b, atom_c]: + if not (0 <= atom < self.natoms): + raise ValueError( + f"Atom index {atom} out of range [0, {self.natoms})" + ) + + indices_a = self.atom_basis_indices.get(atom_a, []) + indices_b = self.atom_basis_indices.get(atom_b, []) + indices_c = self.atom_basis_indices.get(atom_c, []) + + return three_center_delocalization_index( + self.density_matrix, + self.overlap_matrix, + indices_a, + indices_b, + indices_c, + ) + + def get_di_matrix(self) -> np.ndarray: + """Calculate full delocalization index matrix. + + Returns: + Symmetric natoms x natoms matrix of delocalization indices + """ + if self._di_matrix is None: + self._di_matrix = calculate_di_matrix( + self.density_matrix, + self.overlap_matrix, + self.atom_basis_indices, + self.natoms, + ) + return self._di_matrix + + def get_bond_type(self, atom_i: int, atom_j: int) -> str: + """Classify bond type based on delocalization index. + + Args: + atom_i: Index of first atom (0-based) + atom_j: Index of second atom (0-based) + + Returns: + Bond type string + """ + di_val = self.get_delocalization_index(atom_i, atom_j) + return classify_bond_from_di(di_val) + + def get_aromaticity_index(self, ring_atoms: List[int]) -> float: + """Calculate aromaticity index for a ring. + + Args: + ring_atoms: List of atom indices forming the ring (0-based) + + Returns: + Aromaticity index (1.0 = perfect aromatic) + """ + di_matrix = self.get_di_matrix() + return calculate_aromaticity_index(di_matrix, ring_atoms) + + def get_pdi(self, ring_atoms: List[int]) -> float: + """Calculate Para-Delocalization Index for a six-membered ring. + + Args: + ring_atoms: List of 6 atom indices (0-based) + + Returns: + PDI value in electrons + + Raises: + ValueError: If ring doesn't have exactly 6 atoms + """ + di_matrix = self.get_di_matrix() + return calculate_pdi(di_matrix, ring_atoms) + + def get_flu(self, ring_atoms: List[int]) -> float: + """Calculate Aromatic Fluctuation Index for a ring. + + Args: + ring_atoms: List of atom indices forming the ring (0-based) + + Returns: + FLU value (lower = more aromatic) + """ + di_matrix = self.get_di_matrix() + return calculate_flu(di_matrix, ring_atoms) + + def is_aromatic_ring(self, ring_atoms: List[int], threshold: float = 0.85) -> bool: + """Determine if a ring is aromatic based on aromaticity index. + + Args: + ring_atoms: List of atom indices forming the ring (0-based) + threshold: Aromaticity threshold (default 0.85) + + Returns: + True if ring is aromatic, False otherwise + """ + ai = self.get_aromaticity_index(ring_atoms) + return ai >= threshold + + def get_all_dihedral_dis(self, ring_atoms: List[int]) -> Dict[Tuple[int, int], float]: + """Get all delocalization indices for atom pairs in a ring. + + Args: + ring_atoms: List of atom indices forming the ring (0-based) + + Returns: + Dictionary mapping (atom_i, atom_j) pairs to DI values + """ + di_matrix = self.get_di_matrix() + result = {} + + for i, atom_i in enumerate(ring_atoms): + for j, atom_j in enumerate(ring_atoms): + if i < j: + result[(atom_i, atom_j)] = di_matrix[atom_i, atom_j] + + return result diff --git a/tests/bonding/test_delocalization_index.py b/tests/bonding/test_delocalization_index.py new file mode 100644 index 00000000..4b661a26 --- /dev/null +++ b/tests/bonding/test_delocalization_index.py @@ -0,0 +1,551 @@ +"""Test Delocalization Index Implementation (Issue 4). + +This module contains comprehensive tests for delocalization index analysis. +Tests cover 2-center DI, 3-center DI, aromaticity indices, and bond +classification. +""" + +import pytest +import numpy as np +from pathlib import Path + +from pymultiwfn.io import load +from pymultiwfn.bonding.delocalization import ( + delocalization_index, + three_center_delocalization_index, + calculate_di_matrix, + classify_bond_from_di, + calculate_aromaticity_index, + calculate_pdi, + calculate_flu, + DelocalizationIndex, + DelocalizationResult, +) + + +class TestDelocalizationIndexFunction: + """Test the core delocalization_index function.""" + + def test_basic_di_calculation(self): + """Test basic DI calculation with simple matrices.""" + # Simple 2x2 case: two atoms with 1 basis function each + P = np.array([[1.0, 0.5], [0.5, 1.0]]) + S = np.array([[1.0, 0.75], [0.75, 1.0]]) + + di = delocalization_index(P, S, [0], [1]) + + # DI = 2 * |P_01 * S_01| = 2 * |0.5 * 0.75| = 0.75 + assert isinstance(di, float), "DI should be a float" + assert di > 0, "DI should be positive for bonded atoms" + assert abs(di - 0.75) < 0.01, f"Expected DI ~0.75, got {di:.4f}" + + def test_di_with_multiple_basis_functions(self): + """Test DI with multiple basis functions per atom.""" + # 4 basis functions: 2 for atom A, 2 for atom B + P = np.array([ + [1.0, 0.3, 0.5, 0.2], + [0.3, 1.0, 0.2, 0.5], + [0.5, 0.2, 1.0, 0.3], + [0.2, 0.5, 0.3, 1.0], + ]) + S = np.eye(4) # Orthogonal basis + S[0, 2] = S[2, 0] = 0.7 + S[1, 3] = S[3, 1] = 0.6 + + di = delocalization_index(P, S, [0, 1], [2, 3]) + + assert isinstance(di, float), "DI should be a float" + assert di >= 0, "DI should be non-negative" + + def test_di_zero_for_non_bonded(self): + """Test DI is zero when there's no overlap.""" + P = np.array([[1.0, 0.0], [0.0, 1.0]]) + S = np.eye(2) # No off-diagonal overlap + + di = delocalization_index(P, S, [0], [1]) + + assert di == 0.0, "DI should be zero for non-interacting atoms" + + def test_di_empty_indices(self): + """Test DI with empty basis function indices.""" + P = np.eye(2) + S = np.eye(2) + + # Empty list for one atom + di = delocalization_index(P, S, [], [0]) + assert di == 0.0, "DI should be zero with empty indices" + + di = delocalization_index(P, S, [0], []) + assert di == 0.0, "DI should be zero with empty indices" + + def test_di_matrix_shape_validation(self): + """Test that mismatched matrix shapes raise error.""" + P = np.eye(3) + S = np.eye(2) + + with pytest.raises(ValueError, match="shape"): + delocalization_index(P, S, [0], [1]) + + def test_di_index_out_of_range(self): + """Test that out-of-range indices raise error.""" + P = np.eye(2) + S = np.eye(2) + + with pytest.raises(ValueError, match="out of range"): + delocalization_index(P, S, [0], [5]) + + +class TestThreeCenterDI: + """Test 3-center delocalization index.""" + + def test_three_center_di_basic(self): + """Test basic 3-center DI calculation.""" + # Create 3x3 matrices for 3 atoms + P = np.array([ + [1.0, 0.5, 0.3], + [0.5, 1.0, 0.5], + [0.3, 0.5, 1.0], + ]) + S = np.array([ + [1.0, 0.7, 0.5], + [0.7, 1.0, 0.7], + [0.5, 0.7, 1.0], + ]) + + di_3c = three_center_delocalization_index(P, S, [0], [1], [2]) + + assert isinstance(di_3c, float), "3-center DI should be a float" + assert di_3c >= 0, "3-center DI should be non-negative" + + def test_three_center_di_symmetric(self): + """Test 3-center DI is symmetric in atom ordering.""" + P = np.array([ + [1.0, 0.5, 0.3], + [0.5, 1.0, 0.5], + [0.3, 0.5, 1.0], + ]) + S = np.array([ + [1.0, 0.7, 0.5], + [0.7, 1.0, 0.7], + [0.5, 0.7, 1.0], + ]) + + di_3c_abc = three_center_delocalization_index(P, S, [0], [1], [2]) + di_3c_bca = three_center_delocalization_index(P, S, [1], [2], [0]) + di_3c_cab = three_center_delocalization_index(P, S, [2], [0], [1]) + + # All permutations should give same result + assert abs(di_3c_abc - di_3c_bca) < 1e-10, "3C DI should be symmetric" + assert abs(di_3c_abc - di_3c_cab) < 1e-10, "3C DI should be symmetric" + + +class TestDIMatrix: + """Test DI matrix calculation.""" + + def test_di_matrix_shape(self): + """Test DI matrix has correct shape.""" + # 4 basis functions for 2 atoms + P = np.eye(4) + S = np.eye(4) + for i in range(2): + for j in range(2, 4): + P[i, j] = P[j, i] = 0.5 + S[i, j] = S[j, i] = 0.7 + + atom_indices = {0: [0, 1], 1: [2, 3]} + + di_matrix = calculate_di_matrix(P, S, atom_indices, 2) + + assert di_matrix.shape == (2, 2), f"Expected (2,2), got {di_matrix.shape}" + + def test_di_matrix_symmetric(self): + """Test DI matrix is symmetric.""" + P = np.array([ + [1.0, 0.5, 0.3, 0.2], + [0.5, 1.0, 0.2, 0.3], + [0.3, 0.2, 1.0, 0.5], + [0.2, 0.3, 0.5, 1.0], + ]) + S = np.array([ + [1.0, 0.8, 0.6, 0.4], + [0.8, 1.0, 0.4, 0.6], + [0.6, 0.4, 1.0, 0.8], + [0.4, 0.6, 0.8, 1.0], + ]) + atom_indices = {0: [0, 1], 1: [2, 3]} + + di_matrix = calculate_di_matrix(P, S, atom_indices, 2) + + assert np.allclose(di_matrix, di_matrix.T), "DI matrix should be symmetric" + + def test_di_matrix_zero_diagonal(self): + """Test DI matrix has zero diagonal (no self-delocalization).""" + P = np.eye(4) + S = np.eye(4) + atom_indices = {0: [0, 1], 1: [2, 3]} + + di_matrix = calculate_di_matrix(P, S, atom_indices, 2) + + assert np.allclose(np.diag(di_matrix), 0), "Diagonal should be zero" + + +class TestBondClassification: + """Test bond type classification from DI values.""" + + def test_weak_bond(self): + """Test classification of weak bonds.""" + assert classify_bond_from_di(0.1) == "weak" + assert classify_bond_from_di(0.2) == "weak" + + def test_single_bond(self): + """Test classification of single bonds.""" + assert classify_bond_from_di(0.5) == "single" + assert classify_bond_from_di(0.7) == "single" + + def test_aromatic_bond(self): + """Test classification of aromatic bonds.""" + assert classify_bond_from_di(1.0) == "aromatic" + assert classify_bond_from_di(1.2) == "aromatic" + + def test_double_bond(self): + """Test classification of double bonds.""" + assert classify_bond_from_di(1.5) == "double" + assert classify_bond_from_di(1.7) == "double" + + def test_triple_bond(self): + """Test classification of triple bonds.""" + assert classify_bond_from_di(2.0) == "triple" + assert classify_bond_from_di(2.3) == "triple" + + def test_very_strong_bond(self): + """Test classification of very strong bonds.""" + assert classify_bond_from_di(3.0) == "very_strong" + assert classify_bond_from_di(5.0) == "very_strong" + + +class TestAromaticityIndices: + """Test aromaticity index calculations.""" + + def test_aromaticity_index_perfect_aromatic(self): + """Test aromaticity index for perfect aromatic system.""" + # Create ideal aromatic DI matrix (benzene-like) + di_matrix = np.zeros((6, 6)) + # Adjacent bonds have DI = 1.4 (perfect aromatic) + for i in range(6): + j = (i + 1) % 6 + di_matrix[i, j] = di_matrix[j, i] = 1.4 + + ring = [0, 1, 2, 3, 4, 5] + ai = calculate_aromaticity_index(di_matrix, ring) + + # Should be close to 1.0 for perfect aromatic + assert abs(ai - 1.0) < 0.01, f"Expected AI ~1.0, got {ai:.4f}" + + def test_aromaticity_index_non_aromatic(self): + """Test aromaticity index for non-aromatic system.""" + # Create non-aromatic DI matrix (localized bonds) + di_matrix = np.zeros((6, 6)) + # Alternating single/double bonds + for i in range(6): + j = (i + 1) % 6 + di_matrix[i, j] = di_matrix[j, i] = 1.0 if i % 2 == 0 else 1.8 + + ring = [0, 1, 2, 3, 4, 5] + ai = calculate_aromaticity_index(di_matrix, ring) + + # Should be different from 1.0 + assert ai != 1.0, "Non-aromatic should have AI != 1.0" + + def test_pdi_benzene_like(self): + """Test PDI for benzene-like system.""" + di_matrix = np.zeros((6, 6)) + # Adjacent: DI = 1.4, para: DI = 0.1 + for i in range(6): + j_adj = (i + 1) % 6 + di_matrix[i, j_adj] = di_matrix[j_adj, i] = 1.4 + + j_para = (i + 3) % 6 + di_matrix[i, j_para] = di_matrix[j_para, i] = 0.1 + + ring = [0, 1, 2, 3, 4, 5] + pdi = calculate_pdi(di_matrix, ring) + + # PDI should be around 0.1 for this setup + assert 0 < pdi < 0.5, f"Expected small PDI, got {pdi:.4f}" + + def test_pdi_requires_six_atoms(self): + """Test PDI raises error for non-6-membered rings.""" + di_matrix = np.zeros((5, 5)) + ring = [0, 1, 2, 3, 4] + + with pytest.raises(ValueError, match="6"): + calculate_pdi(di_matrix, ring) + + def test_flu_perfect_aromatic(self): + """Test FLU for perfect aromatic system.""" + di_matrix = np.zeros((6, 6)) + # Each atom has DI = 1.4 with each neighbor + for i in range(6): + j_prev = (i - 1) % 6 + j_next = (i + 1) % 6 + di_matrix[i, j_prev] = di_matrix[j_prev, i] = 1.4 + di_matrix[i, j_next] = di_matrix[j_next, i] = 1.4 + + ring = [0, 1, 2, 3, 4, 5] + flu = calculate_flu(di_matrix, ring) + + # FLU should be close to 0 for perfect aromatic + assert flu < 0.1, f"Expected FLU ~0, got {flu:.4f}" + + def test_flu_non_aromatic(self): + """Test FLU for non-aromatic system.""" + di_matrix = np.zeros((6, 6)) + # Localized bonds: alternating strong/weak + for i in range(6): + j_prev = (i - 1) % 6 + j_next = (i + 1) % 6 + # Alternating bond strengths + strength = 1.8 if i % 2 == 0 else 1.0 + di_matrix[i, j_prev] = di_matrix[j_prev, i] = strength + di_matrix[i, j_next] = di_matrix[j_next, i] = strength + + ring = [0, 1, 2, 3, 4, 5] + flu = calculate_flu(di_matrix, ring) + + # FLU should be positive for non-aromatic + assert flu > 0, "FLU should be positive for non-aromatic" + + +class TestDelocalizationIndexClass: + """Test the DelocalizationIndex class.""" + + def test_initialization_with_wavefunction(self, simple_h2_wfn): + """Test initialization with a wavefunction object.""" + di_analyzer = DelocalizationIndex(simple_h2_wfn) + + assert di_analyzer.natoms == 2, "Should have 2 atoms" + assert di_analyzer.density_matrix is not None, "Should have density matrix" + assert di_analyzer.overlap_matrix is not None, "Should have overlap matrix" + + def test_get_delocalization_index(self, simple_h2_wfn): + """Test getting DI between two atoms.""" + di_analyzer = DelocalizationIndex(simple_h2_wfn) + + di_val = di_analyzer.get_delocalization_index(0, 1) + + assert isinstance(di_val, float), "DI should be a float" + assert di_val >= 0, "DI should be non-negative" + + def test_get_di_matrix(self, simple_h2_wfn): + """Test getting full DI matrix.""" + di_analyzer = DelocalizationIndex(simple_h2_wfn) + + di_matrix = di_analyzer.get_di_matrix() + + assert di_matrix.shape == (2, 2), f"Expected (2,2), got {di_matrix.shape}" + assert np.allclose(di_matrix, di_matrix.T), "DI matrix should be symmetric" + + def test_get_three_center_di(self, simple_h3_wfn): + """Test 3-center DI calculation.""" + di_analyzer = DelocalizationIndex(simple_h3_wfn) + + di_3c = di_analyzer.get_three_center_di(0, 1, 2) + + assert isinstance(di_3c, float), "3-center DI should be a float" + assert di_3c >= 0, "3-center DI should be non-negative" + + def test_get_bond_type(self, simple_h2_wfn): + """Test bond type classification.""" + di_analyzer = DelocalizationIndex(simple_h2_wfn) + + bond_type = di_analyzer.get_bond_type(0, 1) + + assert isinstance(bond_type, str), "Bond type should be a string" + assert bond_type in ["weak", "single", "aromatic", "double", "triple", "very_strong"] + + def test_invalid_atom_indices(self, simple_h2_wfn): + """Test error handling for invalid atom indices.""" + di_analyzer = DelocalizationIndex(simple_h2_wfn) + + with pytest.raises(ValueError): + di_analyzer.get_delocalization_index(0, 5) + + def test_same_atom_di(self, simple_h2_wfn): + """Test DI for same atom is zero.""" + di_analyzer = DelocalizationIndex(simple_h2_wfn) + + di_val = di_analyzer.get_delocalization_index(0, 0) + + assert di_val == 0.0, "DI for same atom should be zero" + + def test_is_aromatic_ring(self, benzene_wfn): + """Test aromatic ring detection.""" + di_analyzer = DelocalizationIndex(benzene_wfn) + ring = [0, 1, 2, 3, 4, 5] + + is_aromatic = di_analyzer.is_aromatic_ring(ring) + + assert isinstance(is_aromatic, bool), "Should return boolean" + + def test_get_all_ring_dis(self, benzene_wfn): + """Test getting all DIs in a ring.""" + di_analyzer = DelocalizationIndex(benzene_wfn) + ring = [0, 1, 2, 3, 4, 5] + + all_dis = di_analyzer.get_all_dihedral_dis(ring) + + assert isinstance(all_dis, dict), "Should return dictionary" + assert len(all_dis) == 15, "6-membered ring has 15 unique pairs" + + +class TestDelocalizationResult: + """Test DelocalizationResult dataclass.""" + + def test_result_creation(self): + """Test creating a DelocalizationResult.""" + result = DelocalizationResult(atom_i=0, atom_j=1, di_value=1.4) + + assert result.atom_i == 0 + assert result.atom_j == 1 + assert result.di_value == 1.4 + assert result.bond_type == "unknown" + + def test_result_with_bond_type(self): + """Test creating a DelocalizationResult with bond type.""" + result = DelocalizationResult( + atom_i=0, atom_j=1, di_value=1.4, bond_type="aromatic" + ) + + assert result.bond_type == "aromatic" + + def test_result_repr(self): + """Test string representation of result.""" + result = DelocalizationResult(atom_i=0, atom_j=1, di_value=1.4) + repr_str = repr(result) + + assert "DI" in repr_str + assert "0" in repr_str + assert "1" in repr_str + + +# Pytest fixtures for test data +@pytest.fixture +def simple_h2_wfn(): + """Create a simple H2 wavefunction for testing.""" + from pymultiwfn.core.data import Atom, Shell, Wavefunction + + atoms = [ + Atom(element="H", index=1, x=0.0, y=0.0, z=-0.7, charge=1.0), + Atom(element="H", index=1, x=0.0, y=0.0, z=0.7, charge=1.0), + ] + + shells = [ + Shell(type=0, center_idx=0, exponents=np.array([1.0]), coefficients=np.array([1.0])), + Shell(type=0, center_idx=1, exponents=np.array([1.0]), coefficients=np.array([1.0])), + ] + + coeff = 1.0 / np.sqrt(2) + wfn = Wavefunction( + atoms=atoms, + num_electrons=2.0, + charge=0, + multiplicity=1, + num_basis=2, + num_atomic_orbitals=2, + num_primitives=2, + num_shells=2, + shells=shells, + occupations=np.array([2.0, 0.0]), + coefficients=np.array([[coeff, coeff], [coeff, -coeff]]), + overlap_matrix=np.array([[1.0, 0.75], [0.75, 1.0]]), + ) + wfn.calculate_density_matrices() + + return wfn + + +@pytest.fixture +def simple_h3_wfn(): + """Create a simple H3 wavefunction for 3-center testing.""" + from pymultiwfn.core.data import Atom, Shell, Wavefunction + + atoms = [ + Atom(element="H", index=1, x=0.0, y=0.0, z=-1.0, charge=1.0), + Atom(element="H", index=1, x=0.0, y=0.0, z=0.0, charge=1.0), + Atom(element="H", index=1, x=0.0, y=0.0, z=1.0, charge=1.0), + ] + + shells = [ + Shell(type=0, center_idx=0, exponents=np.array([1.0]), coefficients=np.array([1.0])), + Shell(type=0, center_idx=1, exponents=np.array([1.0]), coefficients=np.array([1.0])), + Shell(type=0, center_idx=2, exponents=np.array([1.0]), coefficients=np.array([1.0])), + ] + + wfn = Wavefunction( + atoms=atoms, + num_electrons=3.0, + charge=0, + multiplicity=2, + num_basis=3, + num_atomic_orbitals=3, + num_primitives=3, + num_shells=3, + shells=shells, + occupations=np.array([2.0, 1.0, 0.0]), + coefficients=np.eye(3), + overlap_matrix=np.array([[1.0, 0.6, 0.3], [0.6, 1.0, 0.6], [0.3, 0.6, 1.0]]), + ) + wfn.calculate_density_matrices() + + return wfn + + +@pytest.fixture +def benzene_wfn(): + """Create a simplified benzene wavefunction for aromaticity testing.""" + from pymultiwfn.core.data import Atom, Shell, Wavefunction + + # Create 6 carbon atoms in a hexagon + atoms = [] + for i in range(6): + angle = i * np.pi / 3 + x = 1.4 * np.cos(angle) + y = 1.4 * np.sin(angle) + z = 0.0 + atoms.append(Atom(element="C", index=6, x=x, y=y, z=z, charge=6.0)) + + # Simplified basis: 1 function per atom + shells = [] + for i in range(6): + shells.append( + Shell(type=0, center_idx=i, exponents=np.array([1.0]), coefficients=np.array([1.0])) + ) + + # Create density and overlap matrices for aromatic system + n_basis = 6 + P = np.zeros((n_basis, n_basis)) + S = np.eye(n_basis) + + # Set up aromatic-like bonding pattern + for i in range(6): + j = (i + 1) % 6 + P[i, j] = P[j, i] = 0.5 # Bond order ~1.5 + S[i, j] = S[j, i] = 0.7 # Overlap between neighbors + + wfn = Wavefunction( + atoms=atoms, + num_electrons=30.0, # 6 carbons * 4 valence electrons (simplified) + charge=0, + multiplicity=1, + num_basis=n_basis, + num_atomic_orbitals=n_basis, + num_primitives=n_basis, + num_shells=n_basis, + shells=shells, + occupations=np.ones(n_basis), + coefficients=np.eye(n_basis), + overlap_matrix=S, + Ptot=P, + ) + + return wfn