compute_force_constants Subroutine

public subroutine compute_force_constants(eigenvalues, reduced_masses, force_constants, force_constants_mdyne)

Compute force constants for each normal mode.

From the harmonic oscillator relation: ω² = k/μ → k = ω² × μ = eigenvalue × μ

Returns force constants in both atomic units (Hartree/Bohr²) and mdyne/Å.

Arguments

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

Eigenvalues from mass-weighted Hessian diagonalization (1/amu)
real(kind=dp), intent(in) :: reduced_masses(:)

Reduced masses in amu
real(kind=dp), intent(out), allocatable :: force_constants(:)

Force constants in atomic units (Hartree/Bohr²)
real(kind=dp), intent(out), optional, allocatable :: force_constants_mdyne(:)

Force constants in mdyne/Å (common experimental unit)

Called by

proc~~compute_force_constants~~CalledByGraph proc~compute_force_constants compute_force_constants proc~compute_vibrational_analysis compute_vibrational_analysis proc~compute_vibrational_analysis->proc~compute_force_constants 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 :: k
integer, private :: n_modes

Source Code

   subroutine compute_force_constants(eigenvalues, reduced_masses, force_constants, &
                                      force_constants_mdyne)
      !! Compute force constants for each normal mode.
      !!
      !! From the harmonic oscillator relation:
      !!   ω² = k/μ  →  k = ω² × μ = eigenvalue × μ
      !!
      !! Returns force constants in both atomic units (Hartree/Bohr²) and mdyne/Å.
      real(dp), intent(in) :: eigenvalues(:)
         !! Eigenvalues from mass-weighted Hessian diagonalization (1/amu)
      real(dp), intent(in) :: reduced_masses(:)
         !! Reduced masses in amu
      real(dp), allocatable, intent(out) :: force_constants(:)
         !! Force constants in atomic units (Hartree/Bohr²)
      real(dp), allocatable, intent(out), optional :: force_constants_mdyne(:)
         !! Force constants in mdyne/Å (common experimental unit)

      integer :: n_modes, k

      n_modes = size(eigenvalues)
      allocate (force_constants(n_modes))

      ! k = eigenvalue × μ (eigenvalue has units Hartree/(Bohr²·amu), μ in amu)
      ! So force_constant has units Hartree/Bohr²
      do k = 1, n_modes
         if (eigenvalues(k) >= 0.0_dp) then
            force_constants(k) = eigenvalues(k)*reduced_masses(k)
         else
            ! Imaginary frequency mode - report absolute value
            force_constants(k) = -abs(eigenvalues(k))*reduced_masses(k)
         end if
      end do

      ! Optionally convert to mdyne/Å
      if (present(force_constants_mdyne)) then
         allocate (force_constants_mdyne(n_modes))
         force_constants_mdyne = force_constants*AU_TO_MDYNE_ANG
      end if

   end subroutine compute_force_constants