compute_moments_of_inertia Subroutine

   subroutine compute_moments_of_inertia(coords, atomic_numbers, n_atoms, &
                                         center_of_mass, moments, principal_axes, is_linear, total_mass)
      !! Compute the principal moments of inertia and detect linear molecules.
      !!
      !! Calculates the center of mass, inertia tensor, and diagonalizes to get
      !! principal moments. A molecule is considered linear if one moment is
      !! essentially zero (< LINEAR_THRESHOLD).
      real(dp), intent(in) :: coords(:, :)       !! Atomic coordinates (3, n_atoms) in Bohr
      integer, intent(in) :: atomic_numbers(:)   !! Atomic numbers
      integer, intent(in) :: n_atoms             !! Number of atoms
      real(dp), intent(out) :: center_of_mass(3)  !! Center of mass in Angstrom
      real(dp), intent(out) :: moments(3)        !! Principal moments in amu*Angstrom^2
      real(dp), intent(out) :: principal_axes(3, 3)  !! Principal axis vectors (columns)
      logical, intent(out) :: is_linear          !! True if molecule is linear
      real(dp), intent(out) :: total_mass        !! Total mass in amu

      real(dp) :: coords_ang(3, n_atoms)
      real(dp) :: masses(n_atoms)
      real(dp) :: rel_coords(3, n_atoms)
      real(dp) :: inertia_tensor(3, 3)
      integer :: i, info
      real(dp) :: mass_i, x, y, z

      ! Convert coordinates to Angstrom
      coords_ang = coords*BOHR_TO_ANGSTROM

      ! Get atomic masses
      total_mass = 0.0_dp
      do i = 1, n_atoms
         masses(i) = element_mass(atomic_numbers(i))
         total_mass = total_mass + masses(i)
      end do

      ! Compute center of mass
      center_of_mass = 0.0_dp
      do i = 1, n_atoms
         center_of_mass(:) = center_of_mass(:) + masses(i)*coords_ang(:, i)
      end do
      center_of_mass = center_of_mass/total_mass

      ! Compute coordinates relative to center of mass
      do i = 1, n_atoms
         rel_coords(:, i) = coords_ang(:, i) - center_of_mass(:)
      end do

      ! Build inertia tensor
      inertia_tensor = 0.0_dp
      do i = 1, n_atoms
         mass_i = masses(i)
         x = rel_coords(1, i)
         y = rel_coords(2, i)
         z = rel_coords(3, i)

         ! Diagonal elements: I_xx = sum(m * (y^2 + z^2)), etc.
         inertia_tensor(1, 1) = inertia_tensor(1, 1) + mass_i*(y*y + z*z)
         inertia_tensor(2, 2) = inertia_tensor(2, 2) + mass_i*(x*x + z*z)
         inertia_tensor(3, 3) = inertia_tensor(3, 3) + mass_i*(x*x + y*y)

         ! Off-diagonal elements: I_xy = -sum(m * x * y), etc.
         inertia_tensor(1, 2) = inertia_tensor(1, 2) - mass_i*x*y
         inertia_tensor(1, 3) = inertia_tensor(1, 3) - mass_i*x*z
         inertia_tensor(2, 3) = inertia_tensor(2, 3) - mass_i*y*z
      end do

      ! Symmetrize
      inertia_tensor(2, 1) = inertia_tensor(1, 2)
      inertia_tensor(3, 1) = inertia_tensor(1, 3)
      inertia_tensor(3, 2) = inertia_tensor(2, 3)

      ! Diagonalize to get principal moments
      principal_axes = inertia_tensor
      call pic_syev(principal_axes, moments, 'V', 'U', info)

      if (info /= 0) then
         call logger%warning("Failed to diagonalize inertia tensor, info = "// &
                             trim(adjustl(to_char(info))))
         moments = 0.0_dp
         is_linear = .false.
         return
      end if

      ! Check for linear molecule: one moment should be ~0
      ! Moments are returned in ascending order
      is_linear = (moments(1) < LINEAR_THRESHOLD)

   end subroutine compute_moments_of_inertia