diff --git a/pymultiwfn/bonding/bonding.py b/pymultiwfn/bonding/bonding.py
index 7fb9040c..73041725 100644
--- a/pymultiwfn/bonding/bonding.py
+++ b/pymultiwfn/bonding/bonding.py
@@ -10,6 +10,7 @@
 
 from ..io import load
 from .fuzzy import FuzzyAtom, fuzzy_bond_order, calculate_fuzzy_bond_order_matrix
+from .intrinsic import intrinsic_bond_order, calculate_intrinsic_bond_order_matrix, IntrinsicBondResult
 
 
 class Bonding:
@@ -114,3 +115,38 @@ def get_fuzzy_bond_order_matrix(self) -> np.ndarray:
                 self.fuzzy_atoms
             )
         return self._fuzzy_matrix
+    
+    def get_intrinsic_bond_order(self, atom_i: int, atom_j: int) -> float:
+        """Calculate intrinsic bond order between two atoms.
+        
+        The intrinsic bond order (IBO) accounts for bond polarity
+        and provides a measure of the intrinsic covalent character
+        of a bond.
+        
+        Args:
+            atom_i: Index of first atom (0-based)
+            atom_j: Index of second atom (0-based)
+            
+        Returns:
+            Intrinsic bond order value
+            
+        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})")
+        if atom_i == atom_j:
+            raise ValueError("Atom indices must be different")
+            
+        return intrinsic_bond_order(self.wfn, atom_i, atom_j)
+    
+    def get_intrinsic_bond_order_matrix(self) -> IntrinsicBondResult:
+        """Calculate intrinsic bond order matrix for all atom pairs.
+        
+        Returns:
+            IntrinsicBondResult containing:
+            - bond_order_matrix: NxN matrix of intrinsic bond orders
+            - polarity_matrix: NxN matrix of bond polarity corrections
+            - wiberg_matrix: NxN matrix of Wiberg bond orders
+        """
+        return calculate_intrinsic_bond_order_matrix(self.wfn)
diff --git a/pymultiwfn/bonding/intrinsic.py b/pymultiwfn/bonding/intrinsic.py
new file mode 100644
index 00000000..2b601703
--- /dev/null
+++ b/pymultiwfn/bonding/intrinsic.py
@@ -0,0 +1,371 @@
+"""Intrinsic Bond Order (IBO) implementation (Issue 3).
+
+This module implements intrinsic bond order analysis with bond polarity
+correction. The intrinsic bond order provides a measure of the intrinsic
+strength of chemical bonds, accounting for ionic character.
+
+The IBO is calculated using the formula:
+    IBO_ij = (PS_ij)^2 * (1 - polarity_correction)
+
+where PS is the density-overlap matrix product and the polarity correction
+accounts for electronegativity differences between atoms.
+
+References:
+- Mayer, I. (1984). Chem. Phys. Lett. 97, 270-274.
+- Wiberg, K. B. (1968). Tetrahedron 24, 1083-1096.
+- Knizhnik, A. et al. (2019). J. Comput. Chem. 40, 1458-1467.
+"""
+
+import numpy as np
+from typing import Dict, List, Tuple, Optional, Union
+from dataclasses import dataclass
+
+
+# Pauling electronegativity values for common elements
+ELECTRONEGATIVITY = {
+    'H': 2.20,  'He': 4.16,
+    'Li': 0.98, 'Be': 1.57, 'B': 2.04,  'C': 2.55,  'N': 3.04,  'O': 3.44,  'F': 3.98,
+    'Ne': 4.79,
+    'Na': 0.93, 'Mg': 1.31, 'Al': 1.61, 'Si': 1.90, 'P': 2.19,  'S': 2.58,  'Cl': 3.16,
+    'Ar': 3.30,
+    'K': 0.82,  'Ca': 1.00, 'Sc': 1.36, 'Ti': 1.54, 'V': 1.63,  'Cr': 1.66,  'Mn': 1.55,
+    'Fe': 1.83, 'Co': 1.88, 'Ni': 1.91, 'Cu': 1.90, 'Zn': 1.65, 'Ga': 1.81,  'Ge': 2.01,
+    'As': 2.18, 'Se': 2.55, 'Br': 2.96, 'Kr': 3.00,
+    'Rb': 0.82, 'Sr': 0.95, 'Y': 1.22,  'Zr': 1.33, 'Nb': 1.6,  'Mo': 2.16,  'Tc': 1.9,
+    'Ru': 2.2,  'Rh': 2.28, 'Pd': 2.20, 'Ag': 1.93, 'Cd': 1.69, 'In': 1.78,  'Sn': 1.96,
+    'Sb': 2.05, 'Te': 2.1,  'I': 2.66,  'Xe': 2.60,
+    'Cs': 0.79, 'Ba': 0.89,
+}
+
+
+@dataclass
+class IntrinsicBondResult:
+    """Results from intrinsic bond order calculation.
+    
+    Attributes:
+        bond_order_matrix: NxN matrix of intrinsic bond orders
+        polarity_matrix: NxN matrix of bond polarity corrections
+        wiberg_matrix: NxN matrix of Wiberg bond orders (for comparison)
+        n_atoms: Number of atoms in the molecule
+    """
+    bond_order_matrix: np.ndarray
+    polarity_matrix: np.ndarray
+    wiberg_matrix: np.ndarray
+    n_atoms: int
+    
+    def get_bond_order(self, atom_i: int, atom_j: int) -> float:
+        """Get intrinsic bond order between two atoms.
+        
+        Args:
+            atom_i: Index of first atom (0-based)
+            atom_j: Index of second atom (0-based)
+            
+        Returns:
+            Intrinsic bond order value
+        """
+        return float(self.bond_order_matrix[atom_i, atom_j])
+    
+    def get_polarity(self, atom_i: int, atom_j: int) -> float:
+        """Get bond polarity correction between two atoms.
+        
+        Args:
+            atom_i: Index of first atom (0-based)
+            atom_j: Index of second atom (0-based)
+            
+        Returns:
+            Bond polarity correction (0-1)
+        """
+        return float(self.polarity_matrix[atom_i, atom_j])
+
+
+def get_electronegativity(element: str) -> float:
+    """Get Pauling electronegativity for an element.
+    
+    Args:
+        element: Atomic symbol (e.g., 'C', 'H', 'O')
+        
+    Returns:
+        Pauling electronegativity value
+        
+    Raises:
+        ValueError: If element is not in the electronegativity table
+    """
+    element = element.capitalize()
+    if element not in ELECTRONEGATIVITY:
+        raise ValueError(f"No electronegativity defined for element: {element}")
+    return ELECTRONEGATIVITY[element]
+
+
+def calculate_bond_polarity(
+    element_i: str,
+    element_j: str,
+    density_matrix: np.ndarray,
+    overlap_matrix: np.ndarray,
+    atom_i_indices: List[int],
+    atom_j_indices: List[int]
+) -> float:
+    """Calculate bond polarity correction between two atoms.
+    
+    The bond polarity correction accounts for the ionic character of a bond
+    based on electronegativity differences and electron density distribution.
+    
+    Args:
+        element_i: Element symbol for atom i
+        element_j: Element symbol for atom j
+        density_matrix: Total density matrix (AO basis)
+        overlap_matrix: Overlap matrix (AO basis)
+        atom_i_indices: List of basis function indices for atom i
+        atom_j_indices: List of basis function indices for atom j
+        
+    Returns:
+        Bond polarity correction factor (0-1)
+        0 = purely covalent, 1 = purely ionic
+    """
+    # Get electronegativity values
+    chi_i = get_electronegativity(element_i)
+    chi_j = get_electronegativity(element_j)
+    
+    # Calculate electronegativity difference
+    delta_chi = abs(chi_i - chi_j)
+    
+    # Calculate Mulliken population on each atom
+    # P_i = sum_{mu in i} (PS)_mu,mu
+    if len(atom_i_indices) > 0 and len(atom_j_indices) > 0:
+        # Extract submatrices
+        PS = density_matrix @ overlap_matrix
+        
+        # Calculate atomic populations
+        pop_i = sum(PS[mu, mu] for mu in atom_i_indices)
+        pop_j = sum(PS[nu, nu] for nu in atom_j_indices)
+        
+        # Calculate charge transfer
+        # Positive value indicates electron transfer from i to j
+        charge_transfer = abs(pop_i - pop_j) / max(pop_i + pop_j, 1e-10)
+    else:
+        charge_transfer = 0.0
+    
+    # Combine electronegativity and charge transfer effects
+    # Pauling's formula: ionic character ≈ 1 - exp(-0.25 * Δχ²)
+    ionic_character = 1.0 - np.exp(-0.25 * delta_chi**2)
+    
+    # Weight by charge transfer
+    polarity = 0.5 * (ionic_character + min(charge_transfer, 1.0))
+    
+    return min(max(polarity, 0.0), 1.0)
+
+
+def calculate_wiberg_bond_order(
+    density_matrix: np.ndarray,
+    overlap_matrix: np.ndarray,
+    atom_i_indices: List[int],
+    atom_j_indices: List[int]
+) -> float:
+    """Calculate Wiberg bond order between two atoms.
+    
+    The Wiberg bond order is defined as:
+        W_ij = sum_{mu in i} sum_{nu in j} (PS)_mu,nu * (PS)_nu,mu
+    
+    Args:
+        density_matrix: Density matrix (AO basis)
+        overlap_matrix: Overlap matrix (AO basis)
+        atom_i_indices: List of basis function indices for atom i
+        atom_j_indices: List of basis function indices for atom j
+        
+    Returns:
+        Wiberg bond order value
+    """
+    if len(atom_i_indices) == 0 or len(atom_j_indices) == 0:
+        return 0.0
+    
+    # Calculate PS matrix
+    PS = density_matrix @ overlap_matrix
+    
+    # Extract submatrices for atoms i and j
+    PS_ij = PS[np.ix_(atom_i_indices, atom_j_indices)]
+    PS_ji = PS[np.ix_(atom_j_indices, atom_i_indices)]
+    
+    # Wiberg bond order: trace(PS_ij @ PS_ji)
+    bond_order = np.trace(PS_ij @ PS_ji)
+    
+    return max(float(bond_order), 0.0)
+
+
+def calculate_intrinsic_bond_order(
+    density_matrix: np.ndarray,
+    overlap_matrix: np.ndarray,
+    element_i: str,
+    element_j: str,
+    atom_i_indices: List[int],
+    atom_j_indices: List[int]
+) -> Tuple[float, float, float]:
+    """Calculate intrinsic bond order between two atoms.
+    
+    The intrinsic bond order is calculated as:
+        IBO_ij = W_ij * (1 - polarity_correction)
+    
+    where W_ij is the Wiberg bond order and the polarity correction
+    accounts for ionic character.
+    
+    Args:
+        density_matrix: Density matrix (AO basis)
+        overlap_matrix: Overlap matrix (AO basis)
+        element_i: Element symbol for atom i
+        element_j: Element symbol for atom j
+        atom_i_indices: List of basis function indices for atom i
+        atom_j_indices: List of basis function indices for atom j
+        
+    Returns:
+        Tuple of (intrinsic_bond_order, polarity, wiberg_bond_order)
+    """
+    # Calculate Wiberg bond order
+    wiberg_bo = calculate_wiberg_bond_order(
+        density_matrix, overlap_matrix, atom_i_indices, atom_j_indices
+    )
+    
+    # Calculate bond polarity
+    polarity = calculate_bond_polarity(
+        element_i, element_j, density_matrix, overlap_matrix,
+        atom_i_indices, atom_j_indices
+    )
+    
+    # Intrinsic bond order = Wiberg BO * (1 - polarity)
+    # This reduces the bond order for polar bonds
+    intrinsic_bo = wiberg_bo * (1.0 - polarity)
+    
+    return intrinsic_bo, polarity, wiberg_bo
+
+
+def calculate_intrinsic_bond_order_matrix(
+    wfn,
+    density_matrix: Optional[np.ndarray] = None
+) -> IntrinsicBondResult:
+    """Calculate intrinsic bond order matrix for all atom pairs.
+    
+    Args:
+        wfn: Wavefunction object containing molecular data
+        density_matrix: Optional density matrix (uses wfn.Ptot if not provided)
+        
+    Returns:
+        IntrinsicBondResult containing bond order matrices
+        
+    Raises:
+        ValueError: If required matrices are not available
+    """
+    # Check required data
+    if wfn.overlap_matrix is None:
+        raise ValueError("Overlap matrix is required for bond order calculation")
+    
+    # Use provided density matrix or get from wavefunction
+    if density_matrix is None:
+        if wfn.Ptot is None:
+            wfn.calculate_density_matrices()
+            if wfn.Ptot is None:
+                raise ValueError("Density matrix could not be calculated")
+        density_matrix = wfn.Ptot
+    
+    n_atoms = wfn.num_atoms
+    
+    # Get atomic basis function indices
+    atom_to_bfs = wfn.get_atomic_basis_indices()
+    
+    # Initialize matrices
+    ibo_matrix = np.zeros((n_atoms, n_atoms))
+    polarity_matrix = np.zeros((n_atoms, n_atoms))
+    wiberg_matrix = np.zeros((n_atoms, n_atoms))
+    
+    # Calculate bond orders for all atom pairs
+    for i in range(n_atoms):
+        bfs_i = atom_to_bfs.get(i, [])
+        
+        for j in range(i + 1, n_atoms):
+            bfs_j = atom_to_bfs.get(j, [])
+            
+            if not bfs_i or not bfs_j:
+                continue
+            
+            # Get element symbols
+            element_i = wfn.atoms[i].element
+            element_j = wfn.atoms[j].element
+            
+            # Calculate intrinsic bond order
+            ibo, polarity, wiberg = calculate_intrinsic_bond_order(
+                density_matrix, wfn.overlap_matrix,
+                element_i, element_j, bfs_i, bfs_j
+            )
+            
+            # Store in matrices (symmetric)
+            ibo_matrix[i, j] = ibo
+            ibo_matrix[j, i] = ibo
+            polarity_matrix[i, j] = polarity
+            polarity_matrix[j, i] = polarity
+            wiberg_matrix[i, j] = wiberg
+            wiberg_matrix[j, i] = wiberg
+    
+    # Set diagonal to zero (no self-bonding)
+    np.fill_diagonal(ibo_matrix, 0.0)
+    np.fill_diagonal(wiberg_matrix, 0.0)
+    
+    return IntrinsicBondResult(
+        bond_order_matrix=ibo_matrix,
+        polarity_matrix=polarity_matrix,
+        wiberg_matrix=wiberg_matrix,
+        n_atoms=n_atoms
+    )
+
+
+def intrinsic_bond_order(
+    wfn,
+    atom_i: int,
+    atom_j: int,
+    density_matrix: Optional[np.ndarray] = None
+) -> float:
+    """Calculate intrinsic bond order between two specific atoms.
+    
+    This is a convenience function for calculating the IBO between
+    a single pair of atoms.
+    
+    Args:
+        wfn: Wavefunction object
+        atom_i: Index of atom i (0-based)
+        atom_j: Index of atom j (0-based)
+        density_matrix: Optional density matrix
+        
+    Returns:
+        Intrinsic bond order value
+        
+    Raises:
+        ValueError: If atom indices are invalid
+    """
+    n_atoms = wfn.num_atoms
+    if not (0 <= atom_i < n_atoms and 0 <= atom_j < n_atoms):
+        raise ValueError(f"Atom indices must be in range [0, {n_atoms})")
+    if atom_i == atom_j:
+        raise ValueError("Atom indices must be different")
+    
+    # Get atomic basis function indices
+    atom_to_bfs = wfn.get_atomic_basis_indices()
+    bfs_i = atom_to_bfs.get(atom_i, [])
+    bfs_j = atom_to_bfs.get(atom_j, [])
+    
+    if not bfs_i or not bfs_j:
+        return 0.0
+    
+    # Use provided density matrix or get from wavefunction
+    if density_matrix is None:
+        if wfn.Ptot is None:
+            wfn.calculate_density_matrices()
+            if wfn.Ptot is None:
+                raise ValueError("Density matrix could not be calculated")
+        density_matrix = wfn.Ptot
+    
+    # Get element symbols
+    element_i = wfn.atoms[atom_i].element
+    element_j = wfn.atoms[atom_j].element
+    
+    # Calculate intrinsic bond order
+    ibo, _, _ = calculate_intrinsic_bond_order(
+        density_matrix, wfn.overlap_matrix,
+        element_i, element_j, bfs_i, bfs_j
+    )
+    
+    return ibo
diff --git a/tests/bonding/test_intrinsic_bond.py b/tests/bonding/test_intrinsic_bond.py
new file mode 100644
index 00000000..049d491b
--- /dev/null
+++ b/tests/bonding/test_intrinsic_bond.py
@@ -0,0 +1,374 @@
+"""Test Intrinsic Bond Order Implementation (Issue 3).
+
+This module contains comprehensive tests for intrinsic bond order analysis.
+Tests cover intrinsic bond strength, bond polarity correction, bond order
+matrix generation, and comparison with Wiberg bond orders.
+"""
+
+import pytest
+import numpy as np
+from pathlib import Path
+
+from pymultiwfn.io import load
+from pymultiwfn.bonding import Bonding
+from pymultiwfn.bonding.intrinsic import (
+    get_electronegativity,
+    calculate_bond_polarity,
+    calculate_wiberg_bond_order,
+    calculate_intrinsic_bond_order,
+    calculate_intrinsic_bond_order_matrix,
+    intrinsic_bond_order,
+    ELECTRONEGATIVITY
+)
+
+
+class TestElectronegativity:
+    """Test electronegativity values and lookups."""
+
+    def test_common_elements(self):
+        """Test electronegativity for common elements."""
+        # Test a few key elements
+        assert abs(get_electronegativity('H') - 2.20) < 0.01
+        assert abs(get_electronegativity('C') - 2.55) < 0.01
+        assert abs(get_electronegativity('N') - 3.04) < 0.01
+        assert abs(get_electronegativity('O') - 3.44) < 0.01
+        assert abs(get_electronegativity('F') - 3.98) < 0.01
+
+    def test_case_insensitive(self):
+        """Test that element lookup is case-insensitive."""
+        assert get_electronegativity('c') == get_electronegativity('C')
+        assert get_electronegativity('o') == get_electronegativity('O')
+
+    def test_invalid_element(self):
+        """Test that invalid elements raise ValueError."""
+        with pytest.raises(ValueError, match="No electronegativity defined"):
+            get_electronegativity('Xx')
+
+
+class TestBondPolarity:
+    """Test bond polarity correction calculations."""
+
+    def test_homodiatomic_polarity(self, sample_wfn):
+        """Test that homodiatomic bonds have low polarity."""
+        wfn = sample_wfn
+        # For H2, polarity should be very low (same element)
+        atom_to_bfs = wfn.get_atomic_basis_indices()
+        
+        if wfn.num_atoms >= 2:
+            polarity = calculate_bond_polarity(
+                'H', 'H',
+                wfn.Ptot, wfn.overlap_matrix,
+                atom_to_bfs[0], atom_to_bfs[1]
+            )
+            assert 0.0 <= polarity < 0.1, \
+                f"H-H polarity should be low, got {polarity:.3f}"
+
+    def test_heterodiatomic_polarity(self, sample_wfn):
+        """Test that heterodiatomic bonds have higher polarity."""
+        wfn = sample_wfn
+        # For bonds between different elements, polarity should be > 0
+        atom_to_bfs = wfn.get_atomic_basis_indices()
+        
+        if wfn.num_atoms >= 2:
+            # Find first pair of different elements
+            for i in range(wfn.num_atoms):
+                for j in range(i + 1, wfn.num_atoms):
+                    elem_i = wfn.atoms[i].element
+                    elem_j = wfn.atoms[j].element
+                    if elem_i != elem_j:
+                        polarity = calculate_bond_polarity(
+                            elem_i, elem_j,
+                            wfn.Ptot, wfn.overlap_matrix,
+                            atom_to_bfs[i], atom_to_bfs[j]
+                        )
+                        assert 0.0 <= polarity <= 1.0, \
+                            f"Polarity should be in [0, 1], got {polarity:.3f}"
+                        return
+
+    def test_polarity_range(self, sample_wfn):
+        """Test that polarity is always in [0, 1]."""
+        wfn = sample_wfn
+        atom_to_bfs = wfn.get_atomic_basis_indices()
+        
+        for i in range(wfn.num_atoms):
+            for j in range(i + 1, wfn.num_atoms):
+                polarity = calculate_bond_polarity(
+                    wfn.atoms[i].element, wfn.atoms[j].element,
+                    wfn.Ptot, wfn.overlap_matrix,
+                    atom_to_bfs[i], atom_to_bfs[j]
+                )
+                assert 0.0 <= polarity <= 1.0, \
+                    f"Polarity {polarity:.3f} not in [0, 1] for {i}-{j}"
+
+
+class TestWibergBondOrder:
+    """Test Wiberg bond order calculations."""
+
+    def test_wiberg_bond_order_positive(self, sample_wfn):
+        """Test that Wiberg bond orders are non-negative."""
+        wfn = sample_wfn
+        atom_to_bfs = wfn.get_atomic_basis_indices()
+        
+        for i in range(wfn.num_atoms):
+            for j in range(i + 1, wfn.num_atoms):
+                wbo = calculate_wiberg_bond_order(
+                    wfn.Ptot, wfn.overlap_matrix,
+                    atom_to_bfs[i], atom_to_bfs[j]
+                )
+                assert wbo >= 0.0, \
+                    f"Wiberg BO should be non-negative, got {wbo:.3f} for {i}-{j}"
+
+    def test_wiberg_symmetric(self, sample_wfn):
+        """Test that Wiberg bond order is symmetric."""
+        wfn = sample_wfn
+        atom_to_bfs = wfn.get_atomic_basis_indices()
+        
+        for i in range(wfn.num_atoms):
+            for j in range(i + 1, wfn.num_atoms):
+                wbo_ij = calculate_wiberg_bond_order(
+                    wfn.Ptot, wfn.overlap_matrix,
+                    atom_to_bfs[i], atom_to_bfs[j]
+                )
+                wbo_ji = calculate_wiberg_bond_order(
+                    wfn.Ptot, wfn.overlap_matrix,
+                    atom_to_bfs[j], atom_to_bfs[i]
+                )
+                assert abs(wbo_ij - wbo_ji) < 1e-10, \
+                    f"Wiberg BO not symmetric: {wbo_ij:.6f} vs {wbo_ji:.6f}"
+
+
+class TestIntrinsicBondOrder:
+    """Test intrinsic bond order calculations."""
+
+    def test_ibo_single_bond(self, h2_fchk):
+        """Test intrinsic bond order for H2 single bond."""
+        bond = Bonding(h2_fchk)
+        ibo = bond.get_intrinsic_bond_order(atom_i=0, atom_j=1)
+
+        # H2 should have bond order close to 1.0
+        assert isinstance(ibo, float), "Bond order should be a float"
+        # IBO should be <= Wiberg BO due to polarity correction
+        assert 0.0 <= ibo <= 1.5, \
+            f"H2 IBO should be in [0, 1.5], got {ibo:.3f}"
+
+    def test_ibo_vs_wiberg(self, sample_wfn):
+        """Test that IBO <= Wiberg BO due to polarity correction."""
+        bond = Bonding(sample_wfn)
+        result = bond.get_intrinsic_bond_order_matrix()
+        
+        ibo_matrix = result.bond_order_matrix
+        wiberg_matrix = result.wiberg_matrix
+        
+        # IBO should be <= Wiberg BO for all bonds
+        for i in range(bond.natoms):
+            for j in range(i + 1, bond.natoms):
+                assert ibo_matrix[i, j] <= wiberg_matrix[i, j] + 1e-6, \
+                    f"IBO {ibo_matrix[i, j]:.3f} > Wiberg {wiberg_matrix[i, j]:.3f} for {i}-{j}"
+
+    def test_ibo_matrix_shape(self, sample_wfn):
+        """Test that IBO matrix has correct shape."""
+        bond = Bonding(sample_wfn)
+        result = bond.get_intrinsic_bond_order_matrix()
+        
+        assert result.bond_order_matrix.shape == (bond.natoms, bond.natoms), \
+            f"Matrix shape {result.bond_order_matrix.shape} != ({bond.natoms}, {bond.natoms})"
+        assert result.polarity_matrix.shape == (bond.natoms, bond.natoms), \
+            f"Polarity matrix shape mismatch"
+        assert result.wiberg_matrix.shape == (bond.natoms, bond.natoms), \
+            f"Wiberg matrix shape mismatch"
+
+    def test_ibo_matrix_symmetric(self, sample_wfn):
+        """Test that IBO matrix is symmetric."""
+        bond = Bonding(sample_wfn)
+        result = bond.get_intrinsic_bond_order_matrix()
+        
+        ibo_matrix = result.bond_order_matrix
+        polarity_matrix = result.polarity_matrix
+        wiberg_matrix = result.wiberg_matrix
+        
+        assert np.allclose(ibo_matrix, ibo_matrix.T), "IBO matrix not symmetric"
+        assert np.allclose(polarity_matrix, polarity_matrix.T), "Polarity matrix not symmetric"
+        assert np.allclose(wiberg_matrix, wiberg_matrix.T), "Wiberg matrix not symmetric"
+
+    def test_ibo_matrix_diagonal_zero(self, sample_wfn):
+        """Test that diagonal elements of IBO matrix are zero."""
+        bond = Bonding(sample_wfn)
+        result = bond.get_intrinsic_bond_order_matrix()
+        
+        assert np.allclose(np.diag(result.bond_order_matrix), 0.0), \
+            "IBO matrix diagonal should be zero"
+        assert np.allclose(np.diag(result.wiberg_matrix), 0.0), \
+            "Wiberg matrix diagonal should be zero"
+
+    def test_atom_indices_validation(self, sample_fchk):
+        """Test validation of atom indices."""
+        bond = Bonding(sample_fchk)
+        wfn = load(sample_fchk)
+
+        # Test out of range atom index
+        with pytest.raises((IndexError, ValueError)):
+            bond.get_intrinsic_bond_order(atom_i=-1, atom_j=0)
+
+        with pytest.raises((IndexError, ValueError)):
+            bond.get_intrinsic_bond_order(atom_i=0, atom_j=wfn.natoms)
+
+
+class TestBondingClassIntegration:
+    """Test integration with Bonding class."""
+
+    def test_bonding_class_has_ibo_methods(self, sample_fchk):
+        """Test that Bonding class has IBO methods."""
+        bond = Bonding(sample_fchk)
+        
+        assert hasattr(bond, 'get_intrinsic_bond_order'), \
+            "Bonding should have get_intrinsic_bond_order method"
+        assert hasattr(bond, 'get_intrinsic_bond_order_matrix'), \
+            "Bonding should have get_intrinsic_bond_order_matrix method"
+
+    def test_bonding_ibo_api(self, sample_fchk):
+        """Test Bonding class IBO API as specified in issue."""
+        bond = Bonding(sample_fchk)
+        
+        # Test the API from the issue description
+        # Note: The issue uses 1-based indexing, but our implementation uses 0-based
+        if bond.natoms >= 2:
+            ibo = bond.get_intrinsic_bond_order(atom_i=0, atom_j=1)
+            assert isinstance(ibo, float), "IBO should be a float"
+            assert ibo >= 0.0, "IBO should be non-negative"
+
+    def test_intrinsic_bond_result_class(self, sample_wfn):
+        """Test IntrinsicBondResult dataclass."""
+        bond = Bonding(sample_wfn)
+        result = bond.get_intrinsic_bond_order_matrix()
+        
+        # Test that result has all required attributes
+        assert hasattr(result, 'bond_order_matrix')
+        assert hasattr(result, 'polarity_matrix')
+        assert hasattr(result, 'wiberg_matrix')
+        assert hasattr(result, 'n_atoms')
+        
+        # Test get_bond_order method
+        if bond.natoms >= 2:
+            ibo = result.get_bond_order(0, 1)
+            assert isinstance(ibo, float)
+            
+        # Test get_polarity method
+        if bond.natoms >= 2:
+            polarity = result.get_polarity(0, 1)
+            assert isinstance(polarity, float)
+            assert 0.0 <= polarity <= 1.0
+
+
+class TestComparisonWithWiberg:
+    """Test comparison between IBO and Wiberg bond orders."""
+
+    def test_ibo_lower_for_polar_bonds(self, water_fchk):
+        """Test that IBO is lower than Wiberg for polar bonds."""
+        if water_fchk is None:
+            pytest.skip("Water test file not available")
+            
+        bond = Bonding(water_fchk)
+        result = bond.get_intrinsic_bond_order_matrix()
+        
+        # For O-H bonds (polar), IBO should be significantly lower than Wiberg
+        # Find O-H bonds (assuming O is atom 0, H are atoms 1 and 2)
+        if bond.natoms >= 3:
+            # Check O-H bonds
+            for h_idx in [1, 2]:
+                ibo = result.bond_order_matrix[0, h_idx]
+                wiberg = result.wiberg_matrix[0, h_idx]
+                
+                # O-H is polar, so IBO should be noticeably lower
+                # (unless the bond is already very weak)
+                if wiberg > 0.5:
+                    assert ibo < wiberg, \
+                        f"O-H IBO {ibo:.3f} should be < Wiberg {wiberg:.3f}"
+
+    def test_ibo_similar_for_nonpolar_bonds(self, sample_wfn):
+        """Test that IBO is similar to Wiberg for nonpolar bonds."""
+        bond = Bonding(sample_wfn)
+        result = bond.get_intrinsic_bond_order_matrix()
+        
+        # For homodiatomic bonds, IBO should be very close to Wiberg
+        # (polarity correction should be minimal)
+        for i in range(bond.natoms):
+            for j in range(i + 1, bond.natoms):
+                elem_i = bond.wfn.atoms[i].element
+                elem_j = bond.wfn.atoms[j].element
+                
+                if elem_i == elem_j:  # Same element
+                    ibo = result.bond_order_matrix[i, j]
+                    wiberg = result.wiberg_matrix[i, j]
+                    polarity = result.polarity_matrix[i, j]
+                    
+                    # Polarity should be low for same element
+                    assert polarity < 0.1, \
+                        f"Polarity {polarity:.3f} too high for {elem_i}-{elem_j}"
+                    
+                    # IBO should be close to Wiberg
+                    if wiberg > 0.1:  # Only check if bond exists
+                        relative_diff = abs(ibo - wiberg) / wiberg
+                        assert relative_diff < 0.1, \
+                            f"IBO {ibo:.3f} differs too much from Wiberg {wiberg:.3f}"
+
+
+# Pytest fixtures for test data
+@pytest.fixture
+def sample_fchk():
+    """Return path to a sample FCHK file."""
+    test_data_dir = Path(__file__).parent.parent.parent / 'test_data' / 'fchk'
+    sample_file = test_data_dir / 'sample.fchk'
+    
+    if not sample_file.exists():
+        # Try alternative locations
+        alternatives = [
+            Path('/home/yhm/software/PyMultiWFN/test_data/fchk/sample.fchk'),
+            Path('/home/yhm/software/PyMultiWFN/test_files/sample.fchk'),
+        ]
+        for alt in alternatives:
+            if alt.exists():
+                return str(alt)
+        
+        pytest.skip(f"Sample FCHK file not found at {sample_file}")
+    
+    return str(sample_file)
+
+
+@pytest.fixture
+def sample_wfn(sample_fchk):
+    """Return a Wavefunction object for testing."""
+    from pymultiwfn.io import load
+    wfn = load(sample_fchk)
+    
+    # Ensure density matrices are calculated
+    if wfn.Ptot is None:
+        wfn.calculate_density_matrices()
+    
+    return wfn
+
+
+@pytest.fixture
+def h2_fchk():
+    """Return path to H2 FCHK file."""
+    test_data_dir = Path(__file__).parent.parent.parent / 'test_data' / 'fchk'
+    h2_file = test_data_dir / 'H2.fchk'
+    
+    if not h2_file.exists():
+        # Try alternative location
+        alt = Path('/home/yhm/software/PyMultiWFN/test_data/fchk/H2.fchk')
+        if alt.exists():
+            return str(alt)
+        pytest.skip(f"H2 FCHK file not found at {h2_file}")
+    
+    return str(h2_file)
+
+
+@pytest.fixture
+def water_fchk():
+    """Return path to H2O FCHK file."""
+    test_data_dir = Path(__file__).parent.parent.parent / 'test_data' / 'fchk'
+    water_file = test_data_dir / 'H2O.fchk'
+    
+    if not water_file.exists():
+        pytest.skip(f"H2O FCHK file not found at {water_file}")
+    return str(water_file)