spintronics 0.3.0

//! Orbital Torque Module
//!
//! This module implements orbital torques on magnetization, which arise when
//! orbital angular momentum currents generated by the Orbital Hall Effect
//! are absorbed by a ferromagnetic layer.
//!
//! ## Physics Background
//!
//! Orbital torques (OT) provide an alternative mechanism to spin-orbit torques
//! (SOT) for current-driven magnetization manipulation. The key advantage is
//! that orbital currents can be generated in light metals without strong SOC,
//! and then converted to effective torques through orbital-to-spin conversion
//! at interfaces.
//!
//! The orbital torque on magnetization has two components:
//!
//! - **Damping-like** (anti-damping): tau_DL ~ m x (m x sigma)
//!   Drives magnetization switching
//! - **Field-like**: tau_FL ~ m x sigma
//!   Acts as an effective field
//!
//! ## Orbital Rashba Effect
//!
//! At interfaces lacking inversion symmetry, the orbital angular momentum
//! can split in momentum space (orbital Rashba effect), analogous to the
//! spin Rashba effect but for orbital degrees of freedom. This leads to
//! the orbital Edelstein effect: charge current generates orbital accumulation
//! at the interface.
//!
//! ## Key References
//!
//! - D. Go and H.-W. Lee, "Orbital torque: Torque generation by orbital
//!   current injection", Phys. Rev. Research 2, 013177 (2020)
//! - S. Lee et al., "Efficient conversion of orbital Hall current to spin
//!   current for spin-orbit torque switching", Commun. Phys. 4, 234 (2021)
//! - D. Go et al., "Theory of current-induced angular momentum transfer
//!   dynamics in spin-orbit coupled systems", Phys. Rev. Research 2, 033401 (2020)
//! - V. Adventures et al., "Orbital Rashba effect in surface oxidized Cu",
//!   Phys. Rev. B 105, L201407 (2022)

use std::fmt;

#[cfg(feature = "serde")]
use serde::{Deserialize, Serialize};

use super::orbital_hall::{OrbitalHallMaterial, OrbitalToSpinConverter};
use crate::constants::{E_CHARGE, GAMMA, HBAR};
use crate::error::{Error, Result};
use crate::vector3::Vector3;

// =============================================================================
// Orbital Torque Calculator
// =============================================================================

/// Orbital torque on magnetization from orbital current injection
///
/// Computes the torque exerted on a ferromagnetic layer by orbital angular
/// momentum currents generated via the Orbital Hall Effect and converted
/// to spin at the interface.
///
/// The total effective field from orbital torque:
///
/// H_OT = (hbar / (2 * e * M_s * d)) * eta_LS * sigma_OH * E
///
/// which decomposes into damping-like and field-like components.
///
/// # Example
///
/// ```
/// use spintronics::orbitronics::orbital_torque::OrbitalTorque;
/// use spintronics::orbitronics::orbital_hall::{OrbitalHallMaterial, OrbitalToSpinConverter};
/// use spintronics::vector3::Vector3;
///
/// let torque = OrbitalTorque::new(
///     OrbitalHallMaterial::chromium(),
///     OrbitalToSpinConverter::cr_pt(),
///     1.0e6,   // M_s = 1 MA/m
///     2.0e-9,  // d = 2 nm FM thickness
///     0.8,     // damping-like efficiency
///     0.2,     // field-like efficiency
/// ).expect("valid parameters");
///
/// let m = Vector3::new(0.0, 0.0, 1.0); // perpendicular magnetization
/// let j_charge = 1.0e11; // A/m^2
/// let current_dir = Vector3::new(1.0, 0.0, 0.0);
///
/// let (tau_dl, tau_fl) = torque.compute_torques(m, j_charge, current_dir);
/// assert!(tau_dl.magnitude() > 0.0);
/// ```
#[derive(Debug, Clone)]
pub struct OrbitalTorque {
    /// OHE material providing orbital current
    pub material: OrbitalHallMaterial,
    /// Orbital-to-spin converter at interface
    pub converter: OrbitalToSpinConverter,
    /// Saturation magnetization of FM layer \[A/m\]
    pub ms: f64,
    /// Ferromagnetic layer thickness \[m\]
    pub thickness_fm: f64,
    /// Damping-like torque efficiency factor (0 to 1)
    pub damping_like_efficiency: f64,
    /// Field-like torque efficiency factor (0 to 1)
    pub field_like_efficiency: f64,
}

impl OrbitalTorque {
    /// Create a new OrbitalTorque calculator
    ///
    /// # Arguments
    /// * `material` - Orbital Hall material
    /// * `converter` - Orbital-to-spin converter
    /// * `ms` - Saturation magnetization \[A/m\]
    /// * `thickness_fm` - FM layer thickness \[m\]
    /// * `damping_like_efficiency` - DL torque efficiency (0 to 1)
    /// * `field_like_efficiency` - FL torque efficiency (0 to 1)
    ///
    /// # Errors
    /// Returns error if physical parameters are invalid
    pub fn new(
        material: OrbitalHallMaterial,
        converter: OrbitalToSpinConverter,
        ms: f64,
        thickness_fm: f64,
        damping_like_efficiency: f64,
        field_like_efficiency: f64,
    ) -> Result<Self> {
        if ms <= 0.0 {
            return Err(Error::InvalidParameter {
                param: "ms".to_string(),
                reason: "saturation magnetization must be positive".to_string(),
            });
        }
        if thickness_fm <= 0.0 {
            return Err(Error::InvalidParameter {
                param: "thickness_fm".to_string(),
                reason: "FM thickness must be positive".to_string(),
            });
        }
        if !(0.0..=1.0).contains(&damping_like_efficiency) {
            return Err(Error::InvalidParameter {
                param: "damping_like_efficiency".to_string(),
                reason: "must be between 0 and 1".to_string(),
            });
        }
        if !(0.0..=1.0).contains(&field_like_efficiency) {
            return Err(Error::InvalidParameter {
                param: "field_like_efficiency".to_string(),
                reason: "must be between 0 and 1".to_string(),
            });
        }
        Ok(Self {
            material,
            converter,
            ms,
            thickness_fm,
            damping_like_efficiency,
            field_like_efficiency,
        })
    }

    /// Create Cr/Pt/CoFeB orbital torque system
    ///
    /// Cr generates giant orbital current via OHE, Pt converts to spin current,
    /// which then exerts torque on the CoFeB ferromagnetic layer.
    pub fn cr_pt_cofeb() -> Self {
        Self {
            material: OrbitalHallMaterial::chromium(),
            converter: OrbitalToSpinConverter::cr_pt(),
            ms: 1.0e6,
            thickness_fm: 2.0e-9,
            damping_like_efficiency: 0.8,
            field_like_efficiency: 0.2,
        }
    }

    /// Create Ti/Pt/CoFeB orbital torque system
    pub fn ti_pt_cofeb() -> Self {
        Self {
            material: OrbitalHallMaterial::titanium(),
            converter: OrbitalToSpinConverter::ti_pt(),
            ms: 1.0e6,
            thickness_fm: 2.0e-9,
            damping_like_efficiency: 0.75,
            field_like_efficiency: 0.25,
        }
    }

    /// Create Cu/Pt/CoFeB orbital torque system
    pub fn cu_pt_cofeb() -> Self {
        Self {
            material: OrbitalHallMaterial::copper(),
            converter: OrbitalToSpinConverter::cu_pt(),
            ms: 1.0e6,
            thickness_fm: 2.0e-9,
            damping_like_efficiency: 0.7,
            field_like_efficiency: 0.15,
        }
    }

    /// Compute the orbital torque prefactor
    ///
    /// Returns (gamma / (M_s * d)) * eta_LS * theta_OH in appropriate units.
    /// This prefactor multiplies the charge current density to give the torque
    /// per unit area.
    ///
    /// The full prefactor: (hbar / (2 * e)) * (gamma / (M_s * d)) * eta_LS * theta_OH
    fn torque_prefactor(&self) -> f64 {
        let hbar_over_2e = HBAR / (2.0 * E_CHARGE);
        let eta_ls = self.converter.conversion_efficiency;
        let theta_oh = self.material.orbital_hall_angle;
        hbar_over_2e * GAMMA * eta_ls * theta_oh / (self.ms * self.thickness_fm)
    }

    /// Compute damping-like and field-like orbital torques
    ///
    /// # Arguments
    /// * `m` - Normalized magnetization vector
    /// * `j_charge` - Charge current density [A/m^2]
    /// * `current_direction` - Unit vector of current flow
    ///
    /// # Returns
    /// (tau_DL, tau_FL) tuple of torque vectors \[rad/s\]
    ///
    /// The damping-like torque: tau_DL = prefactor * j * DL_eff * m x (m x sigma)
    /// The field-like torque: tau_FL = prefactor * j * FL_eff * m x sigma
    pub fn compute_torques(
        &self,
        m: Vector3<f64>,
        j_charge: f64,
        current_direction: Vector3<f64>,
    ) -> (Vector3<f64>, Vector3<f64>) {
        // Spin polarization from OHE: perpendicular to both current and normal
        let normal = Vector3::unit_z();
        let sigma = current_direction.cross(&normal);
        // Normalize sigma if non-zero
        let sigma_mag = sigma.magnitude();
        let sigma = if sigma_mag > 1e-30 {
            sigma * (1.0 / sigma_mag)
        } else {
            Vector3::zero()
        };

        let prefactor = self.torque_prefactor() * j_charge;

        // Damping-like: tau_DL = prefactor * DL_eff * m x (m x sigma)
        let m_cross_sigma = m.cross(&sigma);
        let tau_dl = m.cross(&m_cross_sigma) * (prefactor * self.damping_like_efficiency);

        // Field-like: tau_FL = prefactor * FL_eff * m x sigma
        let tau_fl = m_cross_sigma * (prefactor * self.field_like_efficiency);

        (tau_dl, tau_fl)
    }

    /// Compute effective damping-like field from orbital torque
    ///
    /// H_DL = (hbar/(2e)) * (theta_OH * eta_LS * j / (M_s * d)) * DL_eff
    ///
    /// # Arguments
    /// * `j_charge` - Charge current density [A/m^2]
    /// * `m` - Normalized magnetization vector
    /// * `current_direction` - Unit vector of current flow
    ///
    /// # Returns
    /// Damping-like effective field \[A/m\]
    #[inline]
    pub fn damping_like_field(
        &self,
        j_charge: f64,
        m: Vector3<f64>,
        current_direction: Vector3<f64>,
    ) -> Vector3<f64> {
        let normal = Vector3::unit_z();
        let sigma = current_direction.cross(&normal);
        let sigma_mag = sigma.magnitude();
        let sigma = if sigma_mag > 1e-30 {
            sigma * (1.0 / sigma_mag)
        } else {
            Vector3::zero()
        };

        let hbar_over_2e = HBAR / (2.0 * E_CHARGE);
        let eta_ls = self.converter.conversion_efficiency;
        let theta_oh = self.material.orbital_hall_angle;
        let h_prefactor =
            hbar_over_2e * j_charge * eta_ls * theta_oh * self.damping_like_efficiency
                / (self.ms * self.thickness_fm);

        // H_DL = prefactor * m x (m x sigma)
        let m_cross_sigma = m.cross(&sigma);
        m.cross(&m_cross_sigma) * h_prefactor
    }

    /// Compute effective field-like field from orbital torque
    ///
    /// # Arguments
    /// * `j_charge` - Charge current density [A/m^2]
    /// * `m` - Normalized magnetization vector
    /// * `current_direction` - Unit vector of current flow
    ///
    /// # Returns
    /// Field-like effective field \[A/m\]
    #[inline]
    pub fn field_like_field(
        &self,
        j_charge: f64,
        _m: Vector3<f64>,
        current_direction: Vector3<f64>,
    ) -> Vector3<f64> {
        let normal = Vector3::unit_z();
        let sigma = current_direction.cross(&normal);
        let sigma_mag = sigma.magnitude();
        let sigma = if sigma_mag > 1e-30 {
            sigma * (1.0 / sigma_mag)
        } else {
            Vector3::zero()
        };

        let hbar_over_2e = HBAR / (2.0 * E_CHARGE);
        let eta_ls = self.converter.conversion_efficiency;
        let theta_oh = self.material.orbital_hall_angle;
        let h_prefactor = hbar_over_2e * j_charge * eta_ls * theta_oh * self.field_like_efficiency
            / (self.ms * self.thickness_fm);

        // H_FL = prefactor * sigma (effective field along spin polarization)
        // The field-like torque tau_FL = -gamma * m x H_FL
        // With tau_FL = prefactor * m x sigma, we get H_FL proportional to sigma
        sigma * h_prefactor
    }

    /// Compare orbital torque magnitude with conventional SOT from Pt
    ///
    /// Returns the ratio |tau_OT| / |tau_SOT_Pt| for the same current density.
    /// Values > 1 indicate orbital torque exceeds conventional Pt SOT.
    ///
    /// # Arguments
    /// * `j_charge` - Charge current density [A/m^2]
    /// * `m` - Normalized magnetization
    /// * `current_direction` - Current flow direction
    pub fn ot_to_sot_ratio(
        &self,
        j_charge: f64,
        m: Vector3<f64>,
        current_direction: Vector3<f64>,
    ) -> f64 {
        let (tau_dl_ot, _) = self.compute_torques(m, j_charge, current_direction);

        // Pt SOT reference: theta_SH = 0.07, same geometry
        let theta_sh_pt: f64 = 0.07;
        let transparency_pt: f64 = 0.5;
        let thickness_pt: f64 = 5.0e-9;
        let lambda_sd_pt: f64 = 1.5e-9;
        let tanh_factor = (thickness_pt / lambda_sd_pt).tanh();
        let theta_eff_pt = theta_sh_pt * tanh_factor * transparency_pt;

        let hbar_over_2e = HBAR / (2.0 * E_CHARGE);
        let sot_prefactor =
            hbar_over_2e * j_charge * theta_eff_pt * GAMMA / (self.ms * self.thickness_fm);

        // SOT damping-like field direction is m x (m x sigma)
        let normal = Vector3::unit_z();
        let sigma = current_direction.cross(&normal);
        let sigma_mag = sigma.magnitude();
        let sigma = if sigma_mag > 1e-30 {
            sigma * (1.0 / sigma_mag)
        } else {
            return 0.0;
        };
        let m_cross_sigma = m.cross(&sigma);
        let tau_dl_sot = m.cross(&m_cross_sigma) * sot_prefactor;

        let sot_mag = tau_dl_sot.magnitude();
        if sot_mag > 1e-30 {
            tau_dl_ot.magnitude() / sot_mag
        } else {
            0.0
        }
    }

    /// Estimate critical switching current density for orbital torque switching
    ///
    /// For perpendicular magnetization switching:
    /// j_c = (2e/hbar) * (M_s * d / (eta_LS * theta_OH * DL_eff)) * (H_k + M_s)
    ///
    /// # Arguments
    /// * `h_k` - Perpendicular anisotropy field \[A/m\]
    ///
    /// # Returns
    /// Critical current density [A/m^2]
    pub fn critical_switching_current(&self, h_k: f64) -> f64 {
        let eta_ls = self.converter.conversion_efficiency;
        let theta_oh = self.material.orbital_hall_angle;
        let effective_angle = (eta_ls * theta_oh * self.damping_like_efficiency).abs();

        if effective_angle < 1e-30 {
            return f64::MAX;
        }

        let h_eff = h_k + self.ms;
        (2.0 * E_CHARGE / HBAR) * (self.ms * self.thickness_fm / effective_angle) * h_eff
    }

    /// Compare switching efficiency of OT vs SOT (Pt)
    ///
    /// Returns j_c(SOT_Pt) / j_c(OT). Values > 1 mean OT switches at lower current.
    ///
    /// # Arguments
    /// * `h_k` - Perpendicular anisotropy field \[A/m\]
    pub fn switching_efficiency_vs_sot(&self, h_k: f64) -> f64 {
        let j_c_ot = self.critical_switching_current(h_k);

        // Pt SOT critical current with same FM parameters
        let theta_sh_pt: f64 = 0.07;
        let transparency_pt: f64 = 0.5;
        let thickness_pt: f64 = 5.0e-9;
        let lambda_sd_pt: f64 = 1.5e-9;
        let tanh_factor = (thickness_pt / lambda_sd_pt).tanh();
        let theta_eff_pt = (theta_sh_pt * tanh_factor * transparency_pt).abs();

        let h_eff = h_k + self.ms;
        let j_c_sot =
            (2.0 * E_CHARGE / HBAR) * (self.ms * self.thickness_fm / theta_eff_pt) * h_eff;

        if j_c_ot > 1e-30 {
            j_c_sot / j_c_ot
        } else {
            0.0
        }
    }
}

impl fmt::Display for OrbitalTorque {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            f,
            "OrbitalTorque({}/{}, M_s={:.2e} A/m, d={:.1} nm, DL={:.2}, FL={:.2})",
            self.material.name,
            self.converter.interface_type,
            self.ms,
            self.thickness_fm * 1e9,
            self.damping_like_efficiency,
            self.field_like_efficiency
        )
    }
}

// =============================================================================
// Orbital Rashba Effect
// =============================================================================

/// Orbital Rashba effect at interfaces
///
/// At interfaces lacking inversion symmetry, the orbital angular momentum
/// splits in momentum space, creating an orbital texture. This orbital
/// Rashba effect can generate orbital accumulation from charge current
/// flow (orbital Edelstein effect).
///
/// The orbital Rashba parameter alpha_OR characterizes the strength of
/// the orbital splitting, analogous to the spin Rashba parameter but
/// typically much larger.
///
/// # Physics
///
/// The orbital Rashba Hamiltonian: H_OR = alpha_OR * (k x z_hat) . L
///
/// where L is the orbital angular momentum operator, k is the crystal
/// momentum, and z_hat is the interface normal.
///
/// The orbital Edelstein effect generates orbital accumulation:
/// ⟨L⟩ = alpha_OR * tau * j / (e * hbar)
///
/// where tau is the momentum relaxation time.
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
pub struct OrbitalRashba {
    /// Orbital Rashba parameter \[eV*Angstrom\]
    pub alpha_or: f64,
    /// Interface description
    pub interface: &'static str,
    /// Momentum relaxation time \[s\]
    pub relaxation_time: f64,
    /// Effective orbital quantum number (L_eff, typically 1 for p-orbitals, 2 for d-orbitals)
    pub l_effective: f64,
}

impl OrbitalRashba {
    /// Create a new OrbitalRashba system with validation
    ///
    /// # Arguments
    /// * `alpha_or` - Orbital Rashba parameter \[eV*Angstrom\]
    /// * `interface` - Interface description
    /// * `relaxation_time` - Momentum relaxation time \[s\]
    /// * `l_effective` - Effective orbital quantum number
    ///
    /// # Errors
    /// Returns error for invalid physical parameters
    pub fn new(
        alpha_or: f64,
        interface: &'static str,
        relaxation_time: f64,
        l_effective: f64,
    ) -> Result<Self> {
        if relaxation_time <= 0.0 {
            return Err(Error::InvalidParameter {
                param: "relaxation_time".to_string(),
                reason: "must be positive".to_string(),
            });
        }
        if l_effective < 0.0 {
            return Err(Error::InvalidParameter {
                param: "l_effective".to_string(),
                reason: "must be non-negative".to_string(),
            });
        }
        Ok(Self {
            alpha_or,
            interface,
            relaxation_time,
            l_effective,
        })
    }

    /// Cu/oxide interface orbital Rashba system
    ///
    /// Cu surfaces with oxide overlayers exhibit large orbital Rashba
    /// splitting due to symmetry breaking at the interface.
    /// alpha_OR ~ 1-3 eV*Angstrom, much larger than typical spin Rashba.
    pub fn cu_oxide() -> Self {
        Self {
            alpha_or: 2.0,
            interface: "Cu/oxide",
            relaxation_time: 1.0e-14,
            l_effective: 2.0,
        }
    }

    /// Ag/Bi interface: strong orbital Rashba from heavy element
    pub fn ag_bi() -> Self {
        Self {
            alpha_or: 3.5,
            interface: "Ag/Bi",
            relaxation_time: 0.8e-14,
            l_effective: 1.0,
        }
    }

    /// Al/oxide interface: light metal orbital Rashba
    pub fn al_oxide() -> Self {
        Self {
            alpha_or: 1.0,
            interface: "Al/oxide",
            relaxation_time: 0.5e-14,
            l_effective: 1.0,
        }
    }

    /// Compute the orbital Edelstein effect: orbital accumulation from current
    ///
    /// The orbital accumulation per unit current density is:
    /// ⟨L⟩ / j = alpha_OR * tau / (e * hbar)
    ///
    /// This is in units of [hbar * m / A], representing the orbital
    /// accumulation (in hbar) per unit current density.
    ///
    /// # Returns
    /// Orbital accumulation response coefficient [s/(kg*m)]
    pub fn orbital_edelstein_coefficient(&self) -> f64 {
        // alpha_OR is in eV*Angstrom, convert to J*m:
        // 1 eV = 1.602e-19 J, 1 Angstrom = 1e-10 m
        let alpha_si = self.alpha_or * E_CHARGE * 1e-10;
        alpha_si * self.relaxation_time / (E_CHARGE * HBAR)
    }

    /// Compute orbital accumulation vector from charge current
    ///
    /// ⟨L⟩ = coefficient * (j x z_hat)
    ///
    /// The orbital accumulation is perpendicular to both current flow
    /// and the interface normal.
    ///
    /// # Arguments
    /// * `j_charge` - Charge current density [A/m^2]
    /// * `current_direction` - Unit vector of current flow direction
    ///
    /// # Returns
    /// Orbital accumulation vector [hbar per unit area]
    pub fn orbital_accumulation(
        &self,
        j_charge: f64,
        current_direction: Vector3<f64>,
    ) -> Vector3<f64> {
        let coeff = self.orbital_edelstein_coefficient();
        let z_hat = Vector3::unit_z();
        let j_vec = current_direction * j_charge;
        j_vec.cross(&z_hat) * coeff
    }

    /// Compute orbital Rashba splitting energy at a given k-point
    ///
    /// Delta E = alpha_OR * |k| for the splitting between orbital branches.
    ///
    /// # Arguments
    /// * `k_magnitude` - Crystal momentum magnitude \[1/m\]
    ///
    /// # Returns
    /// Energy splitting \[J\]
    pub fn orbital_splitting(&self, k_magnitude: f64) -> f64 {
        // alpha_OR in eV*Angstrom, convert to J*m
        let alpha_si = self.alpha_or * E_CHARGE * 1e-10;
        alpha_si * k_magnitude
    }

    /// Compute orbital Rashba splitting at the Fermi wavevector
    ///
    /// Uses a typical Fermi wavevector for the given interface material.
    ///
    /// # Arguments
    /// * `k_fermi` - Fermi wavevector \[1/m\]
    ///
    /// # Returns
    /// Splitting energy \[eV\]
    pub fn splitting_at_fermi(&self, k_fermi: f64) -> f64 {
        let splitting_j = self.orbital_splitting(k_fermi);
        splitting_j / E_CHARGE
    }

    /// Compare orbital Rashba with typical spin Rashba
    ///
    /// Returns alpha_OR / alpha_SR ratio. The orbital Rashba parameter
    /// is typically 10-100x larger than the spin Rashba parameter for
    /// the same interface.
    ///
    /// # Arguments
    /// * `alpha_spin_rashba` - Spin Rashba parameter \[eV*Angstrom\]
    pub fn orbital_to_spin_rashba_ratio(&self, alpha_spin_rashba: f64) -> f64 {
        if alpha_spin_rashba.abs() > 1e-30 {
            self.alpha_or / alpha_spin_rashba
        } else {
            f64::MAX
        }
    }
}

impl fmt::Display for OrbitalRashba {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(
            f,
            "OrbitalRashba({}: alpha_OR={:.2} eV*A, tau={:.2e} s, L_eff={:.1})",
            self.interface, self.alpha_or, self.relaxation_time, self.l_effective
        )
    }
}

// =============================================================================
// Predefined orbital torque systems
// =============================================================================

/// Create Cr/Pt/CoFeB orbital torque system
///
/// This system demonstrates giant orbital torque from Cr's OHE,
/// efficiently converted at the Pt interface to exert torque on CoFeB.
pub fn cr_pt_cofeb_system() -> OrbitalTorque {
    OrbitalTorque::cr_pt_cofeb()
}

/// Create Ti/Pt/CoFeB orbital torque system
pub fn ti_pt_cofeb_system() -> OrbitalTorque {
    OrbitalTorque::ti_pt_cofeb()
}

// =============================================================================
// Tests
// =============================================================================

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

    #[test]
    fn test_orbital_torque_magnitude() {
        let ot = OrbitalTorque::cr_pt_cofeb();
        let m = Vector3::new(0.0, 0.0, 1.0); // perpendicular
        let j = 1.0e11; // 10^11 A/m^2
        let dir = Vector3::new(1.0, 0.0, 0.0);

        let (tau_dl, tau_fl) = ot.compute_torques(m, j, dir);

        // Both torques should be non-zero
        assert!(tau_dl.magnitude() > 0.0, "DL torque should be non-zero");
        assert!(tau_fl.magnitude() > 0.0, "FL torque should be non-zero");

        // DL should be larger than FL (DL_eff=0.8 > FL_eff=0.2)
        assert!(
            tau_dl.magnitude() > tau_fl.magnitude(),
            "DL torque ({:.2e}) should exceed FL torque ({:.2e})",
            tau_dl.magnitude(),
            tau_fl.magnitude()
        );
    }

    #[test]
    fn test_orbital_torque_vs_sot_comparison() {
        let ot = OrbitalTorque::cr_pt_cofeb();
        let m = Vector3::new(0.1, 0.0, 1.0).normalize();
        let j = 1.0e11;
        let dir = Vector3::new(1.0, 0.0, 0.0);

        let ratio = ot.ot_to_sot_ratio(j, m, dir);
        // Cr/Pt OT should be comparable to or exceed Pt SOT
        // eta_LS * theta_OH * DL_eff = 0.5 * 0.065 * 0.8 = 0.026
        // Pt SOT effective: theta_SH * tanh * T = 0.07 * tanh(5/1.5) * 0.5 ~ 0.035
        // So ratio ~ 0.026/0.035 ~ 0.74
        assert!(ratio > 0.0, "OT/SOT ratio should be positive: {}", ratio);
        assert!(ratio.is_finite(), "OT/SOT ratio should be finite");
    }

    #[test]
    fn test_orbital_torque_perpendicular_to_m() {
        let ot = OrbitalTorque::cr_pt_cofeb();
        let m = Vector3::new(0.1, 0.2, 0.97).normalize();
        let j = 1.0e11;
        let dir = Vector3::new(1.0, 0.0, 0.0);

        let (tau_dl, _) = ot.compute_torques(m, j, dir);

        // DL torque = m x (m x sigma) should be perpendicular to m
        // since m x (m x sigma) lies in the plane of m and sigma, perpendicular to m
        let dot_m = tau_dl.dot(&m);
        assert!(
            dot_m.abs() < 1e-6,
            "DL torque should be perpendicular to m, dot = {}",
            dot_m
        );
    }

    #[test]
    fn test_light_metal_enhancement() {
        // Compare Cr/Pt system with Cu/Pt system
        let cr_ot = OrbitalTorque::cr_pt_cofeb();
        let cu_ot = OrbitalTorque::cu_pt_cofeb();

        let m = Vector3::new(0.0, 0.0, 1.0);
        let j = 1.0e11;
        let dir = Vector3::new(1.0, 0.0, 0.0);

        let (tau_cr, _) = cr_ot.compute_torques(m, j, dir);
        let (tau_cu, _) = cu_ot.compute_torques(m, j, dir);

        // Cr should produce larger torque than Cu due to larger sigma_OH
        // Cr: eta_LS=0.5, theta_OH=0.065, DL=0.8 => product = 0.026
        // Cu: eta_LS=0.6, theta_OH=0.003, DL=0.7 => product = 0.00126
        assert!(
            tau_cr.magnitude() > tau_cu.magnitude(),
            "Cr should produce larger OT than Cu: Cr={:.2e}, Cu={:.2e}",
            tau_cr.magnitude(),
            tau_cu.magnitude()
        );
    }

    #[test]
    fn test_critical_switching_current() {
        let ot = OrbitalTorque::cr_pt_cofeb();
        let h_k = 5.0e4; // 50 kA/m

        let j_c = ot.critical_switching_current(h_k);

        // Should be positive and finite
        assert!(j_c > 0.0, "Critical current should be positive");
        assert!(j_c.is_finite(), "Critical current should be finite");
        // Should be in physically reasonable range
        assert!(
            j_c > 1.0e10,
            "Critical current should be > 10^10 A/m^2, got {:.2e}",
            j_c
        );
    }

    #[test]
    fn test_switching_efficiency_comparison() {
        let ot = OrbitalTorque::cr_pt_cofeb();
        let h_k = 5.0e4;

        let ratio = ot.switching_efficiency_vs_sot(h_k);

        // The ratio should be positive and finite
        assert!(ratio > 0.0, "Switching efficiency ratio should be positive");
        assert!(
            ratio.is_finite(),
            "Switching efficiency ratio should be finite"
        );
    }

    #[test]
    fn test_damping_like_field_magnitude() {
        let ot = OrbitalTorque::cr_pt_cofeb();
        let m = Vector3::new(0.0, 0.1, 1.0).normalize();
        let j = 1.0e11;
        let dir = Vector3::new(1.0, 0.0, 0.0);

        let h_dl = ot.damping_like_field(j, m, dir);

        // Should be finite and non-zero
        assert!(h_dl.magnitude() > 0.0, "DL field should be non-zero");
        assert!(h_dl.magnitude().is_finite(), "DL field should be finite");
    }

    #[test]
    fn test_orbital_rashba_splitting() {
        let or = OrbitalRashba::cu_oxide();

        // Typical Fermi wavevector for Cu: k_F ~ 1.36e10 /m
        let k_fermi = 1.36e10;
        let splitting_ev = or.splitting_at_fermi(k_fermi);

        // Orbital Rashba splitting should be in meV to eV range
        assert!(splitting_ev > 0.0, "Splitting should be positive");
        assert!(
            splitting_ev < 10.0,
            "Splitting should be < 10 eV for physical reasonableness"
        );

        // For Cu/oxide with alpha_OR=2 eV*A:
        // Delta E = 2 * 1.36e10 * 1e-10 eV = 2 * 1.36 = 2.72 eV
        let expected = 2.0 * 1.36e10 * 1e-10;
        let rel_err = ((splitting_ev - expected) / expected).abs();
        assert!(
            rel_err < 1e-6,
            "Splitting mismatch: got {}, expected {}",
            splitting_ev,
            expected
        );
    }

    #[test]
    fn test_orbital_rashba_vs_spin_rashba() {
        let or = OrbitalRashba::cu_oxide();
        // Typical spin Rashba for Cu/oxide: ~0.01-0.1 eV*Angstrom
        let alpha_sr = 0.05;
        let ratio = or.orbital_to_spin_rashba_ratio(alpha_sr);

        // Orbital Rashba should be much larger (10-100x)
        assert!(
            ratio > 10.0,
            "Orbital Rashba should be >> spin Rashba: ratio = {}",
            ratio
        );
    }

    #[test]
    fn test_orbital_edelstein_accumulation() {
        let or = OrbitalRashba::cu_oxide();
        let j = 1.0e11;
        let dir = Vector3::new(1.0, 0.0, 0.0);

        let accumulation = or.orbital_accumulation(j, dir);

        // Should be transverse to current and interface normal
        // Current along x => accumulation along -y (j x z = -y)
        assert!(
            accumulation.y < 0.0,
            "Accumulation should be along -y for current along x"
        );
        assert!(
            accumulation.x.abs() < 1e-30,
            "Accumulation x should be zero"
        );
        assert!(
            accumulation.z.abs() < 1e-30,
            "Accumulation z should be zero"
        );
    }

    #[test]
    fn test_orbital_torque_parameter_validation() {
        // Valid parameters
        let valid = OrbitalTorque::new(
            OrbitalHallMaterial::chromium(),
            OrbitalToSpinConverter::cr_pt(),
            1.0e6,
            2.0e-9,
            0.8,
            0.2,
        );
        assert!(valid.is_ok(), "Valid parameters should be accepted");

        // Invalid: zero Ms
        let bad_ms = OrbitalTorque::new(
            OrbitalHallMaterial::chromium(),
            OrbitalToSpinConverter::cr_pt(),
            0.0,
            2.0e-9,
            0.8,
            0.2,
        );
        assert!(bad_ms.is_err(), "Zero Ms should be rejected");

        // Invalid: negative thickness
        let bad_d = OrbitalTorque::new(
            OrbitalHallMaterial::chromium(),
            OrbitalToSpinConverter::cr_pt(),
            1.0e6,
            -1.0e-9,
            0.8,
            0.2,
        );
        assert!(bad_d.is_err(), "Negative thickness should be rejected");

        // Invalid: DL efficiency > 1
        let bad_dl = OrbitalTorque::new(
            OrbitalHallMaterial::chromium(),
            OrbitalToSpinConverter::cr_pt(),
            1.0e6,
            2.0e-9,
            1.5,
            0.2,
        );
        assert!(bad_dl.is_err(), "DL efficiency > 1 should be rejected");
    }

    #[test]
    fn test_orbital_rashba_parameter_validation() {
        // Valid
        let valid = OrbitalRashba::new(2.0, "test", 1.0e-14, 2.0);
        assert!(valid.is_ok(), "Valid Rashba parameters should be accepted");

        // Invalid: zero relaxation time
        let bad_tau = OrbitalRashba::new(2.0, "test", 0.0, 2.0);
        assert!(bad_tau.is_err(), "Zero relaxation time should be rejected");

        // Invalid: negative L_eff
        let bad_l = OrbitalRashba::new(2.0, "test", 1.0e-14, -1.0);
        assert!(bad_l.is_err(), "Negative L_eff should be rejected");
    }

    #[test]
    fn test_current_direction_dependence_torque() {
        let ot = OrbitalTorque::cr_pt_cofeb();
        let m = Vector3::new(0.0, 0.0, 1.0);
        let j = 1.0e11;

        // Current along x
        let dir_x = Vector3::new(1.0, 0.0, 0.0);
        let (tau_x, _) = ot.compute_torques(m, j, dir_x);

        // Current along y
        let dir_y = Vector3::new(0.0, 1.0, 0.0);
        let (tau_y, _) = ot.compute_torques(m, j, dir_y);

        // Both should have same magnitude (by symmetry)
        let rel_diff = ((tau_x.magnitude() - tau_y.magnitude()) / tau_x.magnitude()).abs();
        assert!(
            rel_diff < 1e-6,
            "Torque magnitude should be same for x and y current: x={:.2e}, y={:.2e}",
            tau_x.magnitude(),
            tau_y.magnitude()
        );

        // But different directions
        // For current along x: sigma ~ -y, so tau_dl ~ m x (m x (-y))
        // For m = z: m x (-y) = z x (-y) = x, then m x x = z x x = -y
        // => tau_dl_x ~ -y direction
        // For current along y: sigma ~ x, so tau_dl ~ m x (m x x)
        // m x x = z x x = -y, then m x (-y) = z x (-y) = x
        // Wait, m x (m x sigma): first m x sigma, then m x (that)
        // For m=z, sigma=-y: m x sigma = z x (-y) = x, then m x x = z x x = -y => tau_dl along -y? No.
        // Actually m x (m x sigma) = m(m.sigma) - sigma(m.m) = -sigma for unit m perp to sigma
        // So for sigma=-y: tau_dl ~ -(-y) = +y
        // For sigma=+x: tau_dl ~ -(+x) = -x
        // Directions should be orthogonal
        let dot = tau_x.dot(&tau_y);
        assert!(
            dot.abs() / (tau_x.magnitude() * tau_y.magnitude()) < 1e-6,
            "Torques for x and y current should be orthogonal"
        );

        // Current along z should give zero torque (z x z = 0)
        let dir_z = Vector3::new(0.0, 0.0, 1.0);
        let (tau_z, _) = ot.compute_torques(m, j, dir_z);
        assert!(
            tau_z.magnitude() < 1e-20,
            "Current along z should give zero torque"
        );
    }

    #[test]
    fn test_display_formats() {
        let ot = OrbitalTorque::cr_pt_cofeb();
        let display = format!("{}", ot);
        assert!(
            display.contains("Cr"),
            "Display should contain material name"
        );

        let or = OrbitalRashba::cu_oxide();
        let display = format!("{}", or);
        assert!(
            display.contains("Cu/oxide"),
            "Display should contain interface"
        );
    }
}