compute_reduced_masses Subroutine

public subroutine compute_reduced_masses(eigenvectors, element_numbers, reduced_masses)

Compute reduced masses for each normal mode.

The reduced mass μ_k for mode k is defined as: μ_k = 1 / Σ_i (L_mw_{i,k}² / m_i)

where L_mw is the mass-weighted eigenvector (normalized to 1). This formula arises from the relationship Q_k = Σ_i √m_i * x_i * L_mw_{i,k} and ensures that the harmonic oscillator relation ω² = k/μ holds.

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in) :: eigenvectors(:,:)

Mass-weighted eigenvectors from diagonalization (3N x 3N) Columns are normal modes, assumed normalized (Σ_i L²_{i,k} = 1)
integer, intent(in) :: element_numbers(:)

Atomic numbers for each atom (N atoms)
real(kind=dp), intent(out), allocatable :: reduced_masses(:)

Reduced masses in amu (one per mode)

Calls

proc~~compute_reduced_masses~~CallsGraph proc~compute_reduced_masses compute_reduced_masses proc~element_mass element_mass proc~compute_reduced_masses->proc~element_mass

Called by

proc~~compute_reduced_masses~~CalledByGraph proc~compute_reduced_masses compute_reduced_masses proc~compute_vibrational_analysis compute_vibrational_analysis proc~compute_vibrational_analysis->proc~compute_reduced_masses proc~compute_gmbe compute_gmbe proc~compute_gmbe->proc~compute_vibrational_analysis proc~compute_mbe compute_mbe proc~compute_mbe->proc~compute_vibrational_analysis proc~gmbe_pie_coordinator gmbe_pie_coordinator proc~gmbe_pie_coordinator->proc~compute_vibrational_analysis proc~hessian_coordinator hessian_coordinator proc~hessian_coordinator->proc~compute_vibrational_analysis proc~serial_gmbe_pie_processor serial_gmbe_pie_processor proc~serial_gmbe_pie_processor->proc~compute_vibrational_analysis proc~unfragmented_calculation unfragmented_calculation proc~unfragmented_calculation->proc~compute_vibrational_analysis interface~hessian_coordinator hessian_coordinator interface~hessian_coordinator->proc~hessian_coordinator interface~unfragmented_calculation unfragmented_calculation interface~unfragmented_calculation->proc~unfragmented_calculation proc~global_coordinator global_coordinator proc~global_coordinator->proc~compute_mbe proc~gmbe_run_distributed gmbe_context_t%gmbe_run_distributed proc~gmbe_run_distributed->proc~gmbe_pie_coordinator proc~gmbe_run_serial gmbe_context_t%gmbe_run_serial proc~gmbe_run_serial->proc~serial_gmbe_pie_processor proc~serial_fragment_processor serial_fragment_processor proc~serial_fragment_processor->proc~compute_mbe interface~global_coordinator global_coordinator interface~global_coordinator->proc~global_coordinator interface~serial_fragment_processor serial_fragment_processor interface~serial_fragment_processor->proc~serial_fragment_processor proc~distributed_unfragmented_hessian distributed_unfragmented_hessian proc~distributed_unfragmented_hessian->interface~hessian_coordinator proc~run_unfragmented_calculation run_unfragmented_calculation proc~run_unfragmented_calculation->interface~unfragmented_calculation interface~distributed_unfragmented_hessian distributed_unfragmented_hessian interface~distributed_unfragmented_hessian->proc~distributed_unfragmented_hessian proc~mbe_run_distributed mbe_context_t%mbe_run_distributed proc~mbe_run_distributed->interface~global_coordinator proc~mbe_run_serial mbe_context_t%mbe_run_serial proc~mbe_run_serial->interface~serial_fragment_processor proc~run_calculation run_calculation proc~run_calculation->proc~run_unfragmented_calculation

Variables

Type Visibility Attributes Name Initial
integer, private :: iatom
integer, private :: icoord
integer, private :: idx
integer, private :: k
real(kind=dp), private :: mass
integer, private :: n_atoms
integer, private :: n_coords
real(kind=dp), private :: sum_over_mass

Source Code

   subroutine compute_reduced_masses(eigenvectors, element_numbers, reduced_masses)
      !! Compute reduced masses for each normal mode.
      !!
      !! The reduced mass μ_k for mode k is defined as:
      !!   μ_k = 1 / Σ_i (L_mw_{i,k}² / m_i)
      !!
      !! where L_mw is the mass-weighted eigenvector (normalized to 1).
      !! This formula arises from the relationship Q_k = Σ_i √m_i * x_i * L_mw_{i,k}
      !! and ensures that the harmonic oscillator relation ω² = k/μ holds.
      real(dp), intent(in) :: eigenvectors(:, :)
         !! Mass-weighted eigenvectors from diagonalization (3*N x 3*N)
         !! Columns are normal modes, assumed normalized (Σ_i L²_{i,k} = 1)
      integer, intent(in) :: element_numbers(:)
         !! Atomic numbers for each atom (N atoms)
      real(dp), allocatable, intent(out) :: reduced_masses(:)
         !! Reduced masses in amu (one per mode)

      integer :: n_atoms, n_coords, iatom, icoord, k, idx
      real(dp) :: mass, sum_over_mass

      n_atoms = size(element_numbers)
      n_coords = 3*n_atoms

      allocate (reduced_masses(n_coords))

      ! For each normal mode k
      do k = 1, n_coords
         sum_over_mass = 0.0_dp

         ! Sum over all 3N coordinates: Σ_i (L²_{i,k} / m_i)
         do iatom = 1, n_atoms
            mass = element_mass(element_numbers(iatom))
            do icoord = 1, 3
               idx = 3*(iatom - 1) + icoord
               sum_over_mass = sum_over_mass + eigenvectors(idx, k)**2/mass
            end do
         end do

         ! μ_k = 1 / Σ_i (L²_{i,k} / m_i)
         if (sum_over_mass > 1.0e-14_dp) then
            reduced_masses(k) = 1.0_dp/sum_over_mass
         else
            ! Near-zero contribution (e.g., trans/rot mode) - assign a large mass
            reduced_masses(k) = 1.0e10_dp
         end if
      end do

   end subroutine compute_reduced_masses