Address quality stewardship issues

docs/quality/PERFORMANCE_SIMPLIFICATION_REPORT.md

@@ -0,0 +1,29 @@
+# Performance And Simplification Report
+
+## Scope
+
+This pass reviewed parser selection, element mapping, overlap type conversion,
+and bond-order loops because they were repeatedly called out in the open KISS and
+maintainability issues.
+
+## Changes
+
+- Removed repeated element lookup tables from high-use parsers and routed them
+  through `get_atomic_number`.
+- Replaced repeated atom-pair loops in Mayer and Mulliken bond-order calculations
+  with a shared helper.
+- Replaced the long `_type_to_lmn` branch chain with a constant lookup table.
+- Removed empty placeholder files and unsupported VASP grid parser classes.
+
+## Complexity Impact
+
+The implementation removes more lines than it adds in the touched areas. The
+largest reduction comes from parser element lookup tables and placeholder code,
+while the bond-order refactor reduces four nested-loop copies to one shared pair
+iterator.
+
+## Risk
+
+The changes preserve existing public loader and analyzer names where behavior was
+implemented. Unsupported placeholders were removed only after verifying there
+were no internal references.

docs/quality/PHASE3_YAGNI_REVIEW.md

@@ -0,0 +1,25 @@
+# Phase 3 YAGNI Review
+
+## Decision
+
+Phase 3 advanced bond-order work should remain deferred until there is a
+validated user workflow and reference data for each proposed method.
+
+## Evidence
+
+- Existing open issues prioritized correctness, maintainability, tests, and CI
+  reliability over new analysis breadth.
+- Several modules contained placeholder APIs without working behavior. This pass
+  removed the unreferenced placeholders and kept implemented parsers/analyzers
+  available.
+- The new reference validation tests create a safer baseline for deciding which
+  advanced methods are worth implementing next.
+
+## Roadmap
+
+1. Stabilize Mayer, Mulliken, fuzzy, and delocalization analysis behavior with
+   reference tests.
+2. Add user-facing advanced methods only after their inputs, outputs, and
+   validation references are documented.
+3. Keep experimental code out of default imports until it has a passing test
+   matrix and documented limitations.

docs/quality/QUALITY_IMPROVEMENT_PLAN.md

@@ -0,0 +1,43 @@
+# Quality Improvement Plan
+
+This plan tracks the remaining quality issues opened for PyMultiWFN and ties
+each area to a concrete implementation step.
+
+## Testability
+
+- Added focused regression tests under `tests/quality/` for parser selection,
+  element lookup, overlap type mapping, bond-order helpers, and the analysis base
+  contract.
+- Added `tests/validation/` reference-molecule checks so scientific invariants
+  can run in CI without external data files.
+
+## DRY
+
+- Centralized element-symbol lookup in `pymultiwfn.core.definitions`.
+- Reused one atom-pair iterator for Mayer and Mulliken bond-order loops.
+- Introduced `BaseWavefunctionAnalysis` to standardize analyzer construction.
+
+## Logging
+
+- Replaced parser and core warning/debug `print` calls with module loggers.
+- Debug output now stays silent by default and can be enabled by test or caller
+  logging configuration.
+
+## Security
+
+- Removed unused placeholder modules that presented unsupported behavior as
+  importable APIs.
+- Replaced bare exception handling in parser detection with scoped exceptions and
+  debug diagnostics.
+
+## Configuration
+
+- Kept package-level configuration minimal.
+- Removed unsupported VASP grid formats from the parser factory so the advertised
+  supported-format list matches implemented behavior.
+
+## Next Checks
+
+- Keep Gemini review optional until its required credentials are configured.
+- Keep `test-success`, `docs`, and `quality` as the required branch-protection
+  checks for `main`.

docs/quality/SCIENTIFIC_VALIDATION_FRAMEWORK.md

@@ -0,0 +1,32 @@
+# Scientific Validation Framework
+
+## Reference Set
+
+The CI validation framework currently includes five compact diatomic references:
+
+- H2
+- LiH
+- N2
+- O2
+- F2
+
+These are represented as deterministic synthetic wavefunctions in
+`tests/validation/test_reference_molecule_framework.py`.
+
+## Validated Behavior
+
+- Mayer bond-order matrix shape and symmetry.
+- Expected atom-pair bond order for each reference molecule.
+- Valence diagonal consistency through the shared bond-order helper.
+
+## Why Synthetic References
+
+Several legacy tests depend on example WFN files outside the repository. The new
+validation tests avoid external paths so they can run consistently in local
+checks and GitHub Actions.
+
+## Expansion Criteria
+
+Future reference molecules should include a source calculation, expected values,
+tolerance rationale, and at least one regression test that can run without
+network access.

pymultiwfn/analysis/__init__.py

@@ -0,0 +1,5 @@
+"""Analysis interfaces and implementations."""
+
+from pymultiwfn.analysis.base import BaseWavefunctionAnalysis
+
+__all__ = ["BaseWavefunctionAnalysis"]

pymultiwfn/analysis/base.py

@@ -0,0 +1,21 @@
+"""Shared interfaces for wavefunction analysis classes."""
+
+from __future__ import annotations
+
+from abc import ABC
+
+from pymultiwfn.core.data import Wavefunction
+
+
+class BaseWavefunctionAnalysis(ABC):
+    """Base class for analyses that operate on a Wavefunction."""
+
+    def __init__(self, wavefunction: Wavefunction):
+        self.wavefunction = wavefunction
+        self.wfn = wavefunction
+        self.validate_wavefunction()
+
+    def validate_wavefunction(self) -> None:
+        """Validate required wavefunction data before analysis starts."""
+        if self.wavefunction is None:
+            raise ValueError("A Wavefunction instance is required for analysis")

pymultiwfn/analysis/bonding/bondorder.py

@@ -8,11 +8,40 @@
 - Orbital occupancy-perturbed Mayer bond order
 """
 
+from typing import Dict, List, Optional, Tuple, Union
+
 import numpy as np
-from typing import Dict, List, Tuple, Optional, Union
+
 from ...core.data import Wavefunction
 
 
+def _iter_atom_basis_pairs(
+    atom_to_bfs: Dict[int, List[int]], n_atoms: int
+) -> List[Tuple[int, List[int], int, List[int]]]:
+    """Return atom basis-function pairs shared by bond-order algorithms."""
+    pairs = []
+    for i in range(n_atoms):
+        bfs_i = atom_to_bfs.get(i, [])
+        if not bfs_i:
+            continue
+
+        for j in range(i + 1, n_atoms):
+            bfs_j = atom_to_bfs.get(j, [])
+            if bfs_j:
+                pairs.append((i, bfs_i, j, bfs_j))
+    return pairs
+
+
+def _set_bond_order(matrix: np.ndarray, i: int, j: int, value: float) -> None:
+    matrix[i, j] = value
+    matrix[j, i] = value
+
+
+def _set_valence_diagonal(matrix: np.ndarray) -> None:
+    for i in range(matrix.shape[0]):
+        matrix[i, i] = np.sum(matrix[i, :])
+
+
 def calculate_mayer_bond_order(wfn: Wavefunction) -> Dict[str, np.ndarray]:
     """
     Calculate Mayer bond order matrices using optimized NumPy operations.
@@ -41,6 +70,7 @@ def calculate_mayer_bond_order(wfn: Wavefunction) -> Dict[str, np.ndarray]:
 
     n_atoms = wfn.num_atoms
     atom_to_bfs = wfn.get_atomic_basis_indices()
+    atom_pairs = _iter_atom_basis_pairs(atom_to_bfs, n_atoms)
 
     # Calculate P*S matrices
     PS_total = wfn.Ptot @ wfn.overlap_matrix
@@ -49,29 +79,17 @@ def calculate_mayer_bond_order(wfn: Wavefunction) -> Dict[str, np.ndarray]:
     bnd_total = np.zeros((n_atoms, n_atoms))
 
     # Calculate total bond order using vectorized operations
-    for i in range(n_atoms):
-        bfs_i = atom_to_bfs.get(i, [])
-        if not bfs_i:
-            continue
-
-        for j in range(i + 1, n_atoms):
-            bfs_j = atom_to_bfs.get(j, [])
-            if not bfs_j:
-                continue
-
-            # Extract submatrices and compute bond order using vectorized operations
-            ps_ij = PS_total[np.ix_(bfs_i, bfs_j)]
-            ps_ji = PS_total[np.ix_(bfs_j, bfs_i)]
+    for i, bfs_i, j, bfs_j in atom_pairs:
+        # Extract submatrices and compute bond order using vectorized operations
+        ps_ij = PS_total[np.ix_(bfs_i, bfs_j)]
+        ps_ji = PS_total[np.ix_(bfs_j, bfs_i)]
 
-            # Mayer bond order formula: trace(ps_ij @ ps_ji)
-            # This is equivalent to sum_{a in A} sum_{b in B} (PS)_{ab} (PS)_{ba}
-            accum = np.trace(ps_ij @ ps_ji)
-            bnd_total[i, j] = accum
-            bnd_total[j, i] = accum  # Symmetric matrix
+        # Mayer bond order formula: trace(ps_ij @ ps_ji)
+        # This is equivalent to sum_{a in A} sum_{b in B} (PS)_{ab} (PS)_{ba}
+        _set_bond_order(bnd_total, i, j, np.trace(ps_ij @ ps_ji))
 
     # Set diagonal elements as sum of row elements (Mayer valence)
-    for i in range(n_atoms):
-        bnd_total[i, i] = np.sum(bnd_total[i, :])
+    _set_valence_diagonal(bnd_total)
 
     result = {"total": bnd_total}
 
@@ -88,38 +106,23 @@ def calculate_mayer_bond_order(wfn: Wavefunction) -> Dict[str, np.ndarray]:
         bnd_alpha = np.zeros((n_atoms, n_atoms))
         bnd_beta = np.zeros((n_atoms, n_atoms))
 
-        for i in range(n_atoms):
-            bfs_i = atom_to_bfs.get(i, [])
-            if not bfs_i:
-                continue
-
-            for j in range(i + 1, n_atoms):
-                bfs_j = atom_to_bfs.get(j, [])
-                if not bfs_j:
-                    continue
-
-                # Mayer bond order formula for alpha and beta
-                ps_alpha_ij = PS_alpha[np.ix_(bfs_i, bfs_j)]
-                ps_alpha_ji = PS_alpha[np.ix_(bfs_j, bfs_i)]
-                accum_alpha = np.trace(ps_alpha_ij @ ps_alpha_ji)
+        for i, bfs_i, j, bfs_j in atom_pairs:
+            # Mayer bond order formula for alpha and beta
+            ps_alpha_ij = PS_alpha[np.ix_(bfs_i, bfs_j)]
+            ps_alpha_ji = PS_alpha[np.ix_(bfs_j, bfs_i)]
+            ps_beta_ij = PS_beta[np.ix_(bfs_i, bfs_j)]
+            ps_beta_ji = PS_beta[np.ix_(bfs_j, bfs_i)]
 
-                ps_beta_ij = PS_beta[np.ix_(bfs_i, bfs_j)]
-                ps_beta_ji = PS_beta[np.ix_(bfs_j, bfs_i)]
-                accum_beta = np.trace(ps_beta_ij @ ps_beta_ji)
-
-                bnd_alpha[i, j] = accum_alpha
-                bnd_alpha[j, i] = accum_alpha
-                bnd_beta[i, j] = accum_beta
-                bnd_beta[j, i] = accum_beta
+            _set_bond_order(bnd_alpha, i, j, np.trace(ps_alpha_ij @ ps_alpha_ji))
+            _set_bond_order(bnd_beta, i, j, np.trace(ps_beta_ij @ ps_beta_ji))
 
         # Scale by 2 for unrestricted (following Multiwfn convention)
         bnd_alpha = 2 * bnd_alpha
         bnd_beta = 2 * bnd_beta
 
         # Set diagonal elements
-        for i in range(n_atoms):
-            bnd_alpha[i, i] = np.sum(bnd_alpha[i, :])
-            bnd_beta[i, i] = np.sum(bnd_beta[i, :])
+        _set_valence_diagonal(bnd_alpha)
+        _set_valence_diagonal(bnd_beta)
 
         result["alpha"] = bnd_alpha
         result["beta"] = bnd_beta
@@ -152,6 +155,7 @@ def calculate_mulliken_bond_order(wfn: Wavefunction) -> Dict[str, np.ndarray]:
 
     n_atoms = wfn.num_atoms
     atom_to_bfs = wfn.get_atomic_basis_indices()
+    atom_pairs = _iter_atom_basis_pairs(atom_to_bfs, n_atoms)
 
     # Initialize bond order matrices
     bnd_total = np.zeros((n_atoms, n_atoms))
@@ -169,26 +173,15 @@ def calculate_mulliken_bond_order(wfn: Wavefunction) -> Dict[str, np.ndarray]:
             )
 
         # Calculate bond orders using vectorized operations
-        for i in range(n_atoms):
-            bfs_i = atom_to_bfs.get(i, [])
-            if not bfs_i:
-                continue
-
-            for j in range(i + 1, n_atoms):
-                bfs_j = atom_to_bfs.get(j, [])
-                if not bfs_j:
-                    continue
-
-                # Vectorized sum over basis function pairs
-                bnd_total[i, j] = np.sum(PS_total[np.ix_(bfs_i, bfs_j)])
-                bnd_total[j, i] = bnd_total[i, j]  # Symmetric matrix
+        for i, bfs_i, j, bfs_j in atom_pairs:
+            # Vectorized sum over basis function pairs
+            _set_bond_order(bnd_total, i, j, np.sum(PS_total[np.ix_(bfs_i, bfs_j)]))
 
         # Scale by 2 for closed-shell
         bnd_total = 2 * bnd_total
 
         # Set diagonal elements
-        for i in range(n_atoms):
-            bnd_total[i, i] = np.sum(bnd_total[i, :])
+        _set_valence_diagonal(bnd_total)
 
     else:
         # Unrestricted case
@@ -204,30 +197,18 @@ def calculate_mulliken_bond_order(wfn: Wavefunction) -> Dict[str, np.ndarray]:
         bnd_beta = np.zeros((n_atoms, n_atoms))
 
         # Calculate alpha and beta bond orders using vectorized operations
-        for i in range(n_atoms):
-            bfs_i = atom_to_bfs.get(i, [])
-            if not bfs_i:
-                continue
-
-            for j in range(i + 1, n_atoms):
-                bfs_j = atom_to_bfs.get(j, [])
-                if not bfs_j:
-                    continue
-
-                # Vectorized calculations for alpha and beta
-                bnd_alpha[i, j] = np.sum(PS_alpha[np.ix_(bfs_i, bfs_j)])
-                bnd_alpha[j, i] = bnd_alpha[i, j]
-                bnd_beta[i, j] = np.sum(PS_beta[np.ix_(bfs_i, bfs_j)])
-                bnd_beta[j, i] = bnd_beta[i, j]
+        for i, bfs_i, j, bfs_j in atom_pairs:
+            # Vectorized calculations for alpha and beta
+            _set_bond_order(bnd_alpha, i, j, np.sum(PS_alpha[np.ix_(bfs_i, bfs_j)]))
+            _set_bond_order(bnd_beta, i, j, np.sum(PS_beta[np.ix_(bfs_i, bfs_j)]))
 
         # Scale by 2 for unrestricted
         bnd_alpha = 2 * bnd_alpha
         bnd_beta = 2 * bnd_beta
 
         # Set diagonal elements
-        for i in range(n_atoms):
-            bnd_alpha[i, i] = np.sum(bnd_alpha[i, :])
-            bnd_beta[i, i] = np.sum(bnd_beta[i, :])
+        _set_valence_diagonal(bnd_alpha)
+        _set_valence_diagonal(bnd_beta)
 
         bnd_total = bnd_alpha + bnd_beta
 

pymultiwfn/analysis/orbitals/__init__.py

@@ -4,6 +4,7 @@
 
 import numpy as np
 
+from pymultiwfn.analysis.base import BaseWavefunctionAnalysis
 from pymultiwfn.core.data import Wavefunction
 from pymultiwfn.math.basis import evaluate_basis
 
@@ -12,11 +13,14 @@
 # Import localization methods
 
 
-class OrbitalAnalyzer:
+class OrbitalAnalyzer(BaseWavefunctionAnalysis):
     def __init__(self, wavefunction: Wavefunction):
-        if wavefunction.coefficients is None:
+        super().__init__(wavefunction)
+
+    def validate_wavefunction(self) -> None:
+        super().validate_wavefunction()
+        if self.wfn.coefficients is None:
             raise ValueError("Wavefunction object must contain MO coefficients.")
-        self.wfn = wavefunction
 
     def get_orbital_properties(self, mo_idx: int, is_beta: bool = False) -> dict:
         """