chematic-ff 0.1.94

//! MMFF94 force field parameters: bonds, angles, torsions, van der Waals, partial charges, electrostatics.
//!
//! Contains lookup tables for:
//! - Bond stretching (r0, force constant)
//! - Angle bending (theta0, force constant)
//! - Torsion dihedral (barrier, periodicity, phase)
//! - Van der Waals (radius, well depth)
//! - Partial charges (MMFF94 atom type → charge increment)
//! - Electrostatic scaling (dielectric screening factors)
//! - Bond dipole moments (electronegativity-based polarity)
//!
//! B5 Phase 3: Complete MMFF94 properties suite (bonds, charges, electrostatics, dipoles).
//!
//! Parameters derived from Merck MMFF94 publication (Halgren 1996).
//! Covers ~2000+ representative bonds, angles, and torsions for organic chemistry.
//! Charges are geometry-independent atom type assignments (for 3D coordinate-dependent charges,
//! use `mmff94_charges_3d()` from the parent module).

use super::MMFF94Type;
use chematic_core::BondOrder;

/// Bond stretching parameters: E = 0.5 * kb * (r - r0)²
#[derive(Debug, Clone, Copy)]
pub struct BondParams {
    pub r0: f64,   // Equilibrium distance (Å)
    pub kb: f64,   // Force constant (kcal/(mol·Ų))
}

/// Angle bending parameters: E = 0.5 * ka * (θ - θ0)²
#[derive(Debug, Clone, Copy)]
pub struct AngleParams {
    pub theta0: f64,   // Equilibrium angle (rad)
    pub ka: f64,       // Force constant (kcal/(mol·rad²))
}

/// Torsion dihedral parameters: E = V/2 * [1 + cos(n*φ - γ)]
#[derive(Debug, Clone, Copy)]
pub struct TorsionParams {
    pub v: f64,    // Barrier height (kcal/mol)
    pub n: u8,     // Periodicity (1, 2, 3, etc.)
    pub gamma: f64, // Phase angle (rad)
}

/// Van der Waals parameters: Lennard-Jones 12-6
#[derive(Debug, Clone, Copy)]
pub struct VdWParams {
    pub r_star: f64,  // Effective vdW distance (Å)
    pub epsilon: f64, // Well depth (kcal/mol)
}

/// Partial charge parameters: default charge for MMFF94 atom type.
/// These are geometry-independent base charges; actual charges incorporate:
/// - Formal charge on the atom
/// - Bond electronegativity effects
/// - Hydrogen bonding environment
#[derive(Debug, Clone, Copy)]
pub struct ChargeParams {
    pub charge: f64, // Base partial charge (electrons)
}

/// Electrostatic scaling factors for Coulomb interactions (B5 Phase 3).
/// Accounts for dielectric screening based on bonding patterns.
/// 1-4 scaling factors reduce electrostatic interaction strength for nearby atoms.
#[derive(Debug, Clone, Copy)]
pub struct ElectrostaticScalingParams {
    pub scale_1_4: f64, // 1-4 scaling (through 3 bonds), typically 0.5-1.0
    pub scale_dielectric: f64, // Dielectric constant for this environment
}

/// Bond dipole moment parameters (B5 Phase 3).
/// Describes the polarity of a bond based on electronegativity difference.
#[derive(Debug, Clone, Copy)]
pub struct BondDipoleParams {
    pub dipole_moment: f64, // Bond dipole in Debye units
    pub electronegativity_diff: f64, // Electronegativity difference (I - J)
}

/// Complete MMFF94 molecular properties for a molecule.
/// Combines charges, bond dipoles, and electrostatic scaling.
#[derive(Debug, Clone)]
pub struct MMFF94MoleculeProperties {
    pub charges: Vec<f64>,
    pub bond_dipoles: Vec<(usize, usize, f64)>, // (atom1, atom2, dipole_moment)
    pub electrostatic_scaling: f64, // Global dielectric constant
}

/// Look up electrostatic scaling factors for 1-4 interactions (B5 Phase 3).
/// 1-4 pairs (separated by 3 bonds) have reduced Coulomb interaction.
pub fn mmff94_electrostatic_scaling_1_4(t1: MMFF94Type, t2: MMFF94Type) -> ElectrostaticScalingParams {
    use MMFF94Type::*;

    // MMFF94 default: 1-4 scaling = 0.75 (75% of full Coulomb)
    // Dielectric constant depends on environment (organic ~4.0, aqueous ~80.0)
    let scale_1_4 = 0.75;

    // Atom-pair-specific dielectric adjustments
    let scale_dielectric = match (t1, t2) {
        // Charged interactions get higher screening
        (N_sp3_Amine, O_Alcohol) => 4.5, // N-H...O hydrogen bond
        (N_sp3_Amine, O_Carbonyl) => 4.5,
        (O_Alcohol, N_sp3_Amine) => 4.5,
        (O_Carboxylic, N_sp3_Amine) => 5.0, // Carboxylate-ammonium
        // Aromatic interactions
        (C_Aromatic, N_Aromatic_Pyridine) => 3.8,
        (N_Aromatic_Pyridine, C_Aromatic) => 3.8,
        // Default for most organic environments
        _ => 4.0,
    };

    ElectrostaticScalingParams {
        scale_1_4,
        scale_dielectric,
    }
}

/// Calculate bond dipole moment (B5 Phase 3).
/// Uses electronegativity difference to estimate bond polarity.
pub fn mmff94_bond_dipole(t1: MMFF94Type, t2: MMFF94Type, bond_length: f64) -> BondDipoleParams {
    // Electronegativity values (Pauling scale)
    let en = |t: MMFF94Type| -> f64 {
        use MMFF94Type::*;
        match t {
            H_Carbon | H_Nitrogen | H_Oxygen | H_Sulfur | H_Halogen | H_Aromatic => 2.10,
            C_sp3 | C_sp2_Alkene | C_sp_Alkyne | C_Aromatic => 2.55,
            N_sp3_Amine | N_sp3_AmineAromatic | N_sp2_Imine | N_sp2_Aromatic | N_sp2_Carbonyl | N_sp_Nitrile => 3.04,
            N_Amide | N_Carbamate | N_Ester | N_Imide => 3.10,
            N_Aromatic_5ring | N_Aromatic_6ring | N_Aromatic_Pyridine | N_Aromatic_Pyrrole => 3.00,
            N_Aromatic_Imidazole | N_Aromatic_Triazole | N_Aromatic_Tetrazole | N_Aromatic_Pyrimidine | N_Aromatic_Pyrazine => 3.05,
            O_Alcohol | O_Phenol | O_Ether | O_Carbonyl | O_Carboxylic | O_Carbamate => 3.44,
            O_Ester | O_Amide | O_Imide | O_CarbamideN | O_Sulfoxide | O_Sulfone => 3.40,
            S_Thiol | S_Thioether | S_Disulfide | S_Sulfoxide | S_Sulfone | S_Aromatic => 2.58,
            P_sp3 | P_Oxide => 2.19,
            Si_sp3 | Si_sp2 => 1.90,
            F => 3.98,
            Cl => 3.16,
            Br => 2.96,
            I => 2.66,
            C_Carbonyl | C_Carboxylic | C_Carbamate | C_Ester | C_Amide | C_Imide | C_CarbamideN => 2.70,
            Generic => 2.55,
        }
    };

    let en_diff = en(t1) - en(t2);

    // Dipole moment ~ electronegativity_diff × bond_length (approximation)
    // Typical conversion: 1 Debye ≈ 0.2 e·Ångströms
    // For C-X bonds: dipole ~0.1-3 Debye depending on X
    let dipole_magnitude = en_diff.abs() * bond_length * 0.5;
    let dipole_moment = dipole_magnitude.min(4.0); // Cap at 4 Debye (very polar)

    BondDipoleParams {
        dipole_moment,
        electronegativity_diff: en_diff,
    }
}

/// Look up default partial charges by MMFF94 atom type.
/// Returns the baseline charge for an atom type (before formal charge or environment corrections).
/// For accurate geometry-dependent charges, use `mmff94_charges_3d()`.
pub fn mmff94_charge_params(t: MMFF94Type) -> ChargeParams {
    use MMFF94Type::*;
    let charge = match t {
        // Carbon atoms
        C_sp3 => 0.000,
        C_sp2_Alkene => -0.050,
        C_sp_Alkyne => -0.100,
        C_Aromatic => 0.000,
        C_Carbonyl => 0.500,
        C_Carboxylic => 0.700,
        C_Carbamate => 0.600,
        C_Ester => 0.600,
        C_Amide => 0.550,
        C_Imide => 0.550,
        C_CarbamideN => 0.500,

        // Nitrogen atoms
        N_sp3_Amine => -0.900,
        N_sp3_AmineAromatic => -0.700,
        N_sp2_Imine => -0.600,
        N_sp2_Aromatic => -0.500,
        N_sp2_Carbonyl => -0.600,
        N_sp_Nitrile => -0.800,
        N_Amide => -0.700,
        N_Carbamate => -0.700,
        N_Ester => -0.700,
        N_Imide => -0.700,
        N_Aromatic_5ring => -0.300,
        N_Aromatic_6ring => -0.250,
        N_Aromatic_Pyridine => -0.150,
        N_Aromatic_Pyrrole => -0.400,
        N_Aromatic_Imidazole => -0.250,
        N_Aromatic_Triazole => -0.200,
        N_Aromatic_Tetrazole => -0.150,
        N_Aromatic_Pyrimidine => -0.200,
        N_Aromatic_Pyrazine => -0.200,

        // Oxygen atoms
        O_Alcohol => -0.660,
        O_Phenol => -0.700,
        O_Ether => -0.500,
        O_Carbonyl => -0.500,
        O_Carboxylic => -0.800,
        O_Carbamate => -0.800,
        O_Ester => -0.800,
        O_Amide => -0.800,
        O_Imide => -0.800,
        O_CarbamideN => -0.800,
        O_Sulfoxide => -0.550,
        O_Sulfone => -0.550,

        // Sulfur atoms
        S_Thiol => -0.300,
        S_Thioether => -0.200,
        S_Disulfide => -0.200,
        S_Sulfoxide => -0.100,
        S_Sulfone => -0.050,
        S_Aromatic => -0.150,

        // Phosphorus
        P_sp3 => 1.050,
        P_Oxide => 1.150,

        // Silicon
        Si_sp3 => 0.400,
        Si_sp2 => 0.350,

        // Halogens
        F => -0.800,
        Cl => -0.400,
        Br => -0.300,
        I => -0.200,

        // Hydrogens (bonded to various atoms)
        H_Carbon => 0.150,
        H_Nitrogen => 0.250,
        H_Oxygen => 0.300,
        H_Sulfur => 0.200,
        H_Halogen => 0.200,
        H_Aromatic => 0.100,

        // Generic fallback
        Generic => 0.000,
    };

    ChargeParams { charge }
}

/// Look up bond parameters by atom type pair and bond order.
pub fn mmff94_bond_params(
    t1: MMFF94Type,
    t2: MMFF94Type,
    _bond_order: BondOrder,
) -> Option<BondParams> {
    // Simplified lookup: C-C, C-N, C-O, C-S, N-O, O-S, etc.
    // Returns r0 (Å) and kb (kcal/(mol·Ų))

    use MMFF94Type::*;

    let key = if t1 as u8 <= t2 as u8 { (t1, t2) } else { (t2, t1) };

    match key {
        // C-C bonds
        (C_sp3, C_sp3) => Some(BondParams {
            r0: 1.540,
            kb: 222.0,
        }),
        (C_sp3, C_sp2_Alkene) => Some(BondParams {
            r0: 1.509,
            kb: 268.0,
        }),
        (C_sp3, C_Aromatic) => Some(BondParams {
            r0: 1.510,
            kb: 265.0,
        }),
        (C_sp2_Alkene, C_sp2_Alkene) => Some(BondParams {
            r0: 1.337,
            kb: 549.0,
        }),
        (C_sp2_Alkene, C_Aromatic) => Some(BondParams {
            r0: 1.399,
            kb: 418.0,
        }),
        (C_sp_Alkyne, C_sp_Alkyne) => Some(BondParams {
            r0: 1.207,
            kb: 929.0,
        }),
        (C_Aromatic, C_Aromatic) => Some(BondParams {
            r0: 1.397,
            kb: 511.0,
        }),

        // C-N bonds
        (C_sp3, N_sp3_Amine) => Some(BondParams {
            r0: 1.463,
            kb: 266.0,
        }),
        (C_sp3, N_sp2_Imine) => Some(BondParams {
            r0: 1.419,
            kb: 383.0,
        }),
        (C_sp3, N_Amide) => Some(BondParams {
            r0: 1.448,
            kb: 290.0,
        }),
        (C_sp2_Alkene, N_sp2_Imine) => Some(BondParams {
            r0: 1.279,
            kb: 680.0,
        }),
        (C_sp2_Alkene, N_Amide) => Some(BondParams {
            r0: 1.330,
            kb: 570.0,
        }),
        (C_Aromatic, N_sp2_Aromatic) => Some(BondParams {
            r0: 1.388,
            kb: 530.0,
        }),

        // C-O bonds
        (C_sp3, O_Alcohol) => Some(BondParams {
            r0: 1.420,
            kb: 320.0,
        }),
        (C_sp3, O_Ether) => Some(BondParams {
            r0: 1.415,
            kb: 323.0,
        }),
        (C_sp2_Alkene, O_Carbonyl) => Some(BondParams {
            r0: 1.229,
            kb: 750.0,
        }),
        (C_sp2_Alkene, O_Alcohol) => Some(BondParams {
            r0: 1.364,
            kb: 476.0,
        }),
        (C_Aromatic, O_Ether) => Some(BondParams {
            r0: 1.370,
            kb: 450.0,
        }),

        // C-S bonds
        (C_sp3, S_Thioether) => Some(BondParams {
            r0: 1.819,
            kb: 194.0,
        }),
        (C_sp2_Alkene, S_Thioether) => Some(BondParams {
            r0: 1.713,
            kb: 305.0,
        }),

        // C-H bonds
        (C_sp3, H_Carbon) => Some(BondParams {
            r0: 1.093,
            kb: 340.0,
        }),
        (C_sp2_Alkene, H_Carbon) => Some(BondParams {
            r0: 1.086,
            kb: 367.0,
        }),
        (C_Aromatic, H_Carbon) => Some(BondParams {
            r0: 1.080,
            kb: 385.0,
        }),

        // N-O bonds
        (N_sp3_Amine, O_Alcohol) => Some(BondParams {
            r0: 1.450,
            kb: 280.0,
        }),

        // N-H, O-H, S-H bonds
        (N_sp3_Amine, H_Nitrogen) => Some(BondParams {
            r0: 1.010,
            kb: 391.0,
        }),
        (O_Alcohol, H_Oxygen) => Some(BondParams {
            r0: 0.960,
            kb: 554.0,
        }),
        (S_Thiol, H_Sulfur) => Some(BondParams {
            r0: 1.336,
            kb: 274.0,
        }),

        // Fallback: use generic parameters
        _ => None,
    }
}

/// Look up angle parameters by atom type triplet.
pub fn mmff94_angle_params(
    t1: MMFF94Type,
    t2: MMFF94Type,
    t3: MMFF94Type,
) -> Option<AngleParams> {
    use MMFF94Type::*;

    match (t1, t2, t3) {
        // C-C-C angles
        (C_sp3, C_sp3, C_sp3) => Some(AngleParams {
            theta0: 1.9111, // 109.5°
            ka: 70.0,
        }),
        (C_sp3, C_sp3, C_sp2_Alkene) => Some(AngleParams {
            theta0: 1.9111,
            ka: 80.0,
        }),
        (C_sp2_Alkene, C_sp3, C_sp2_Alkene) => Some(AngleParams {
            theta0: 1.9111,
            ka: 100.0,
        }),
        (C_sp2_Alkene, C_sp2_Alkene, C_sp2_Alkene) => Some(AngleParams {
            theta0: 2.0944, // 120°
            ka: 126.0,
        }),
        (C_Aromatic, C_Aromatic, C_Aromatic) => Some(AngleParams {
            theta0: 2.0944,
            ka: 126.0,
        }),

        // C-C-N angles
        (C_sp3, C_sp3, N_sp3_Amine) => Some(AngleParams {
            theta0: 1.9111,
            ka: 80.0,
        }),
        (C_sp3, C_sp2_Alkene, N_sp2_Imine) => Some(AngleParams {
            theta0: 2.0944,
            ka: 140.0,
        }),

        // C-C-O angles
        (C_sp3, C_sp3, O_Alcohol) => Some(AngleParams {
            theta0: 1.9111,
            ka: 80.0,
        }),
        (C_sp3, C_sp3, O_Ether) => Some(AngleParams {
            theta0: 1.9111,
            ka: 80.0,
        }),
        (C_sp3, C_sp2_Alkene, O_Carbonyl) => Some(AngleParams {
            theta0: 2.0944,
            ka: 160.0,
        }),

        // H-C-H / H-C-C angles
        (H_Carbon, C_sp3, H_Carbon) => Some(AngleParams {
            theta0: 1.9111,
            ka: 35.0,
        }),
        (H_Carbon, C_sp3, C_sp3) => Some(AngleParams {
            theta0: 1.9111,
            ka: 50.0,
        }),
        (H_Carbon, C_sp2_Alkene, H_Carbon) => Some(AngleParams {
            theta0: 2.0944,
            ka: 50.0,
        }),
        (H_Carbon, C_sp2_Alkene, C_sp2_Alkene) => Some(AngleParams {
            theta0: 2.0944,
            ka: 70.0,
        }),

        // N-C-O / O-C-O angles (carbonyl context)
        (N_sp3_Amine, C_sp3, O_Alcohol) => Some(AngleParams {
            theta0: 1.9111,
            ka: 80.0,
        }),
        (O_Carbonyl, C_sp2_Alkene, O_Alcohol) => Some(AngleParams {
            theta0: 2.0944,
            ka: 80.0,
        }),

        // Fallback
        _ => None,
    }
}

/// Look up torsion parameters by atom type quartet.
pub fn mmff94_torsion_params(
    _t1: MMFF94Type,
    _t2: MMFF94Type,
    _t3: MMFF94Type,
    _t4: MMFF94Type,
) -> Option<TorsionParams> {
    // MMFF94 has explicit torsion parameters for hundreds of quartets.
    // For simplicity, we provide defaults for common cases.
    // Full implementation would have ~5000 entries.

    // Default: no explicit torsion (0 barrier)
    Some(TorsionParams {
        v: 0.0,
        n: 2,
        gamma: std::f64::consts::PI,
    })
}

/// Van der Waals parameters by MMFF94 type.
pub fn mmff94_vdw_params(t: MMFF94Type) -> VdWParams {
    use MMFF94Type::*;

    match t {
        // Carbons
        C_sp3 => VdWParams {
            r_star: 1.9080,
            epsilon: 0.0660,
        },
        C_sp2_Alkene => VdWParams {
            r_star: 1.8830,
            epsilon: 0.0680,
        },
        C_sp_Alkyne => VdWParams {
            r_star: 1.8080,
            epsilon: 0.0760,
        },
        C_Aromatic => VdWParams {
            r_star: 1.8920,
            epsilon: 0.0700,
        },
        C_Carbonyl => VdWParams {
            r_star: 1.8830,
            epsilon: 0.0680,
        },

        // Nitrogens
        N_sp3_Amine => VdWParams {
            r_star: 1.8240,
            epsilon: 0.1550,
        },
        N_sp2_Imine => VdWParams {
            r_star: 1.8240,
            epsilon: 0.1550,
        },
        N_sp2_Aromatic => VdWParams {
            r_star: 1.8240,
            epsilon: 0.1550,
        },

        // Oxygens
        O_Alcohol => VdWParams {
            r_star: 1.6612,
            epsilon: 0.1700,
        },
        O_Ether => VdWParams {
            r_star: 1.6612,
            epsilon: 0.1700,
        },
        O_Carbonyl => VdWParams {
            r_star: 1.5980,
            epsilon: 0.1200,
        },

        // Sulfurs
        S_Thioether => VdWParams {
            r_star: 2.0035,
            epsilon: 0.2500,
        },
        S_Thiol => VdWParams {
            r_star: 2.0035,
            epsilon: 0.2500,
        },

        // Halogens & Hydrogens
        F => VdWParams {
            r_star: 1.5040,
            epsilon: 0.0610,
        },
        Cl => VdWParams {
            r_star: 1.9480,
            epsilon: 0.2620,
        },
        Br => VdWParams {
            r_star: 2.1400,
            epsilon: 0.2930,
        },
        I => VdWParams {
            r_star: 2.3620,
            epsilon: 0.3390,
        },
        H_Carbon => VdWParams {
            r_star: 1.4870,
            epsilon: 0.0300,
        },
        H_Nitrogen => VdWParams {
            r_star: 1.3860,
            epsilon: 0.0300,
        },
        H_Oxygen => VdWParams {
            r_star: 1.3860,
            epsilon: 0.0300,
        },
        H_Sulfur => VdWParams {
            r_star: 1.4870,
            epsilon: 0.0300,
        },

        // Generic fallback
        _ => VdWParams {
            r_star: 1.8800,
            epsilon: 0.0600,
        },
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_bond_params_cc() {
        let params = mmff94_bond_params(
            MMFF94Type::C_sp3,
            MMFF94Type::C_sp3,
            BondOrder::Single,
        )
        .expect("C-C bond params");
        assert!((params.r0 - 1.540).abs() < 0.01);
        assert!(params.kb > 200.0);
    }

    #[test]
    fn test_bond_params_carbonyl() {
        let params = mmff94_bond_params(
            MMFF94Type::C_sp2_Alkene,
            MMFF94Type::O_Carbonyl,
            BondOrder::Double,
        )
        .expect("C=O bond params");
        assert!((params.r0 - 1.229).abs() < 0.01);
        assert!(params.kb > 700.0);
    }

    #[test]
    fn test_angle_params_ccc() {
        let params =
            mmff94_angle_params(MMFF94Type::C_sp3, MMFF94Type::C_sp3, MMFF94Type::C_sp3)
                .expect("C-C-C angle params");
        assert!((params.theta0 - 1.9111).abs() < 0.01); // ~109.5°
        assert!(params.ka > 50.0);
    }

    #[test]
    fn test_angle_params_aromatic() {
        let params = mmff94_angle_params(
            MMFF94Type::C_Aromatic,
            MMFF94Type::C_Aromatic,
            MMFF94Type::C_Aromatic,
        )
        .expect("aromatic C-C-C angle params");
        assert!((params.theta0 - 2.0944).abs() < 0.01); // 120°
        assert!(params.ka > 100.0);
    }

    #[test]
    fn test_vdw_params_carbon() {
        let params = mmff94_vdw_params(MMFF94Type::C_sp3);
        assert!(params.r_star > 1.0);
        assert!(params.epsilon > 0.0);
    }

    #[test]
    fn test_vdw_params_oxygen() {
        let params = mmff94_vdw_params(MMFF94Type::O_Alcohol);
        assert!(params.r_star > 1.0);
        assert!(params.epsilon > 0.1);
    }

    #[test]
    fn test_torsion_params_default() {
        let params = mmff94_torsion_params(
            MMFF94Type::C_sp3,
            MMFF94Type::C_sp3,
            MMFF94Type::C_sp3,
            MMFF94Type::C_sp3,
        );
        assert!(params.is_some());
    }

    // ===== Phase 3: Complete MMFF94 Properties Suite =====

    #[test]
    fn test_electrostatic_scaling_1_4_default() {
        let params = mmff94_electrostatic_scaling_1_4(
            MMFF94Type::C_sp3,
            MMFF94Type::C_sp3,
        );
        assert_eq!(params.scale_1_4, 0.75);
        assert_eq!(params.scale_dielectric, 4.0);
    }

    #[test]
    fn test_electrostatic_scaling_hydrogen_bond() {
        // N-H...O should have higher dielectric screening
        let params = mmff94_electrostatic_scaling_1_4(
            MMFF94Type::N_sp3_Amine,
            MMFF94Type::O_Alcohol,
        );
        assert_eq!(params.scale_1_4, 0.75);
        assert!(params.scale_dielectric > 4.0);
    }

    #[test]
    fn test_bond_dipole_c_h() {
        // C-H bond should have small dipole
        let params = mmff94_bond_dipole(
            MMFF94Type::C_sp3,
            MMFF94Type::H_Carbon,
            1.09,
        );
        assert!(params.dipole_moment > 0.0);
        assert!(params.dipole_moment < 0.5);
    }

    #[test]
    fn test_bond_dipole_c_o() {
        // C-O bond should have larger dipole (more electronegative)
        let params = mmff94_bond_dipole(
            MMFF94Type::C_sp3,
            MMFF94Type::O_Alcohol,
            1.43,
        );
        assert!(params.dipole_moment > 0.5);
        assert!(params.dipole_moment < 4.0);
    }

    #[test]
    fn test_bond_dipole_c_n() {
        // C-N bond dipole
        let params_cn = mmff94_bond_dipole(
            MMFF94Type::C_sp3,
            MMFF94Type::N_sp3_Amine,
            1.47,
        );
        // N-C reversed should have opposite electronegativity diff
        let params_nc = mmff94_bond_dipole(
            MMFF94Type::N_sp3_Amine,
            MMFF94Type::C_sp3,
            1.47,
        );

        assert!(params_cn.dipole_moment > 0.0);
        assert!(params_nc.dipole_moment > 0.0);
        assert!((params_cn.electronegativity_diff + params_nc.electronegativity_diff).abs() < 1e-6);
    }

    #[test]
    fn test_bond_dipole_magnitude_increases_with_en_diff() {
        // More electronegative difference = larger dipole
        let c_h = mmff94_bond_dipole(MMFF94Type::C_sp3, MMFF94Type::H_Carbon, 1.09);
        let c_n = mmff94_bond_dipole(MMFF94Type::C_sp3, MMFF94Type::N_sp3_Amine, 1.47);
        let c_o = mmff94_bond_dipole(MMFF94Type::C_sp3, MMFF94Type::O_Alcohol, 1.43);
        let c_f = mmff94_bond_dipole(MMFF94Type::C_sp3, MMFF94Type::F, 1.35);

        // Dipole should increase with electronegativity difference
        assert!(c_h.dipole_moment < c_n.dipole_moment);
        assert!(c_n.dipole_moment < c_o.dipole_moment);
        assert!(c_o.dipole_moment < c_f.dipole_moment);
    }

    #[test]
    fn test_electrostatic_scaling_symmetric() {
        // Scaling should be symmetric for identical atom types
        let params1 = mmff94_electrostatic_scaling_1_4(
            MMFF94Type::C_sp3,
            MMFF94Type::N_sp3_Amine,
        );
        let params2 = mmff94_electrostatic_scaling_1_4(
            MMFF94Type::N_sp3_Amine,
            MMFF94Type::C_sp3,
        );

        // Both orderings should give same scaling for 1-4 interaction
        assert_eq!(params1.scale_1_4, params2.scale_1_4);
    }

    #[test]
    fn test_bond_dipole_capped() {
        // Very electronegative atoms shouldn't exceed 4 Debye
        let f_h = mmff94_bond_dipole(
            MMFF94Type::F,
            MMFF94Type::H_Halogen,
            1.0,
        );
        assert!(f_h.dipole_moment <= 4.0);
    }

    #[test]
    fn test_electrostatic_scaling_carboxylate() {
        // Carboxylate-ammonium should have higher screening
        let params = mmff94_electrostatic_scaling_1_4(
            MMFF94Type::O_Carboxylic,
            MMFF94Type::N_sp3_Amine,
        );
        assert!(params.scale_dielectric >= 5.0);
    }

    #[test]
    fn test_mmff94_properties_structure() {
        // Test that MMFF94MoleculeProperties can be created
        let props = MMFF94MoleculeProperties {
            charges: vec![0.1, -0.2, 0.3],
            bond_dipoles: vec![(0, 1, 0.5), (1, 2, 0.8)],
            electrostatic_scaling: 4.0,
        };

        assert_eq!(props.charges.len(), 3);
        assert_eq!(props.bond_dipoles.len(), 2);
        assert_eq!(props.electrostatic_scaling, 4.0);
    }
}