phonon_modes.F90

Go to the documentation of this file.

!! Copyright (C) 2024 - 2026 M. Lueders
!!
!! This Source Code Form is subject to the terms of the Mozilla Public
!! License, v. 2.0. If a copy of the MPL was not distributed with this
!! file, You can obtain one at https://mozilla.org/MPL/2.0/.
!!
#include "global.h"
 
 
module phonon_modes_oct_m
 
  use, intrinsic :: iso_fortran_env
  use classical_modes_oct_m
  use debug_oct_m
  use global_oct_m
  use io_oct_m
  use namespace_oct_m
  use messages_oct_m
  use parser_oct_m
  use profiling_oct_m
  use unit_system_oct_m
  use wigner_distribution_oct_m
 
  implicit none
 
  private
  public :: &
    phonon_modes_t
 
 
  !
  type, extends(classical_modes_t) :: phonon_modes_t
 
    real(real64), allocatable :: alpha(:)
    real(real64), allocatable :: l_m(:)
    real(real64), allocatable :: prefactor(:)
 
    integer                   :: num_super
 
    type(namespace_t)         :: namespace
 
  contains
 
    procedure :: init => phonon_modes_init 
    procedure :: sample => phonon_modes_sample 
    procedure :: get_displacements => phonon_modes_get_displacements 
    final :: phonon_modes_finalize 
 
  end type phonon_modes_t
 
contains
 
  !
  subroutine phonon_modes_init(this, namespace, dim_space, num_atoms, periodic, masses)
    class(phonon_modes_t),  intent(inout) :: this
    type(namespace_t),      intent(in)    :: namespace
    integer,                intent(in)    :: dim_space
    integer,                intent(in)    :: num_atoms
    logical,                intent(in)    :: periodic
    real(real64),           intent(in)    :: masses(:)
 
    character(len=256) :: filename, version_string
    character(len=800) :: line
    character(len=128) :: dummy
    integer :: iunit, iatom, idim, i
    integer :: imode, num_modes, file_num_modes, file_num_atoms, num_zero_modes, max_modes
    logical :: eof
 
    real(real64), parameter :: frequency_default_threshold = 1.0e-5_real64
    real(real64), parameter :: mass_rel_tolerance = 1.0e-3_real64
    real(real64) :: frequency_threshold
 
    real(real64) :: eigenvecs(1:dim_space*num_atoms)
    real(real64) :: file_masses(1:num_atoms)
    real(real64) :: alpha, frequency
 
 
    push_sub(phonon_modes_init)
 
    this%namespace = namespace
    this%dim = dim_space * num_atoms
    num_modes = 0
 
    !%Variable PhononModesFile
    !%Type string
    !%Section System
    !%Description
    !% Filename for the phonon modes file, as generated by the oct-phonopy-eigenmodes
    !% utility shipped with Octopus. This file is in plain text with the following format (version 1.0;
    !% lines starting with # are treated as comments and are ignored, except within the
    !% mass and eigenvector value blocks):
    !% <tt>
    !% Version: 1.0
    !% Nmodes: <integer>  # Number of modes in the file
    !% Natoms: <integer>  # Number of atoms in the supercell
    !% Np: <integer>      # Number of primitive cells in the supercell, used to generate the modes
    !% Masses:            # followed by Natoms atomic masses (in AMU, free format)
    !% <float> <float> ...
    !% frequency: <float>           # Angular frequency of first mode (in atomic units)
    !% <float> <float> <float>      # Normal mode eigenvector (for first atom)
    !% <float> <float> <float>      # Normal mode eigenvector (for second atom)
    !% ...                          # repeated for all atoms in the supercell
    !% alpha: <float>               # amplitude weight of the mode: g_m^(-1/2), i.e.
    !%                              # 1 for region A (2q = G) and finite systems, 1/sqrt(2) for region B.
    !%                              # For finite systems this line is optional and defaults to 1.
    !% frequency: <float>           # Frequency of second mode
    !% ...
    !% </tt>
    !% The order of atoms has to be the same as in the input (or geometry) file; the masses
    !% are used to verify consistency between the phonon calculation and this run.
    !%
    !% Note: the full Brillouin zone can be split into 3 regions:
    !%       - region A: (time reversal invariant points): -q = q + G
    !%       - region B: -q + G falls in region C
    !%       - region C: -q + G falls in region B
    !% Eigenvectors of regions B and C can be expressed as real and imaginary components of those of region B only.
    !%End
    call parse_variable(namespace, "PhononModesFile", "NONE", filename)
 
    !%Variable PhononModesZeroThreshold
    !%Type float
    !%Default 1.0e-5
    !%Section System
    !%Description
    !% Frequency threshold (in atomic units), below which phonon modes should be considered as zero frequency modes,
    !% such as translations or rotations, and be ignored for the generation of initial displacements.
    !% This default corresponds to 1e-5 Ha = 0.27 meV = 0.065 THz.
    !%End
    call parse_variable(namespace, "PhononModesZeroThreshold", frequency_default_threshold, frequency_threshold)
 
    if (trim(filename) /= "NONE") then
 
      iunit = io_open(filename, action="read")
 
      call next_data_line(iunit, line, filename, namespace)
      read(line, *) dummy, version_string
      if (trim(dummy) /= 'Version:' .or. trim(version_string) /= '1.0') then
        message(1) = "Unsupported phonon file format in "//trim(filename)//" (expected 'Version: 1.0')."
        message(2) = "Regenerate the file with the current oct-phonopy-eigenmodes utility."
        call messages_fatal(2, namespace=namespace)
      end if
 
      call next_data_line(iunit, line, filename, namespace)
      read(line, *) dummy, file_num_modes
      if (trim(dummy) /= 'Nmodes:') then
        message(1) = "Phonon file ("//trim(filename)//") is ill formatted (expected: 'Nmodes:')."
        call messages_fatal(1, namespace=namespace)
      end if
 
      call next_data_line(iunit, line, filename, namespace)
      read(line, *) dummy, file_num_atoms
      if (trim(dummy) /= 'Natoms:') then
        message(1) = "Phonon file ("//trim(filename)//") is ill formatted (expected 'Natoms:')."
        call messages_fatal(1, namespace=namespace)
      end if
      if (file_num_atoms /= num_atoms) then
        write(message(1), '(A,I0,A,I0,A)') "The phonon file ("//trim(filename)//") was generated for ", &
          file_num_atoms, " atoms, but this system has ", num_atoms, " atoms."
        call messages_fatal(1, namespace=namespace)
      end if
 
      call next_data_line(iunit, line, filename, namespace)
      read(line, *) dummy, this%num_super
      if (trim(dummy) /= 'Np:') then
        message(1) = "Phonon file ("//trim(filename)//") is ill formatted (expected 'Np:')."
        call messages_fatal(1, namespace=namespace)
      end if
      if (this%num_super < 1) then
        write(message(1), '(A,I0,A)') "The phonon file ("//trim(filename)//") declares Np = ", &
          this%num_super, ", but Np must be at least 1."
        call messages_fatal(1, namespace=namespace)
      end if
 
      call next_data_line(iunit, line, filename, namespace)
      read(line, *) dummy
      if (trim(dummy) /= 'Masses:') then
        message(1) = "Phonon file ("//trim(filename)//") is ill formatted (expected 'Masses:')."
        call messages_fatal(1, namespace=namespace)
      end if
      read(iunit, *) file_masses(:)
      do iatom = 1, num_atoms
        if (abs(file_masses(iatom) - masses(iatom)) > mass_rel_tolerance * masses(iatom)) then
          write(message(1), '(A,I0,A)') "Mass mismatch between the phonon file ("//trim(filename)// &
            ") and this run for atom ", iatom, ":"
          write(message(2), '(A,F14.8,A,F14.8,A)') "  file: ", file_masses(iatom), " AMU, system: ", &
            masses(iatom), " AMU."
          message(3) = "The phonon calculation and the MTEF run must use identical masses."
          call messages_fatal(3, namespace=namespace)
        end if
      end do
 
      safe_allocate(this%frequencies(1:file_num_modes))
      safe_allocate(this%alpha(1:file_num_modes))
      safe_allocate(this%eigenvectors(1:this%dim, file_num_modes))
 
      num_zero_modes = 0
 
      do imode=1, file_num_modes
        call next_data_line(iunit, line, filename, namespace)
        read(line, *) dummy, frequency
        if(trim(dummy) /= 'frequency:') then
          message(1) = "Phonon file is ill formatted. ("//trim(filename)//")"
          call messages_fatal(1, namespace=namespace)
        end if
        read(iunit, *) eigenvecs(:)
        if (periodic) then
          call next_data_line(iunit, line, filename, namespace)
          read(line, *) dummy, alpha
          if(trim(dummy) /= 'alpha:') then
            message(1) = "Phonon file is ill formatted. ("//trim(filename)//")"
            call messages_fatal(1, namespace=namespace)
          end if
        else
          ! For finite systems the alpha line is optional and defaults to 1
          ! (all modes are single real modes, g_m = 1).
          alpha = 1.0_real64
          call next_data_line(iunit, line, filename, namespace, eof)
          if (.not. eof) then
            read(line, *) dummy
            if (trim(dummy) == 'alpha:') then
              read(line, *) dummy, alpha
            else
              ! Not an alpha line: push it back for the next mode.
              backspace(iunit)
            end if
          end if
        end if
        if (frequency > frequency_threshold) then
          num_modes = num_modes + 1
          this%frequencies(num_modes) = frequency
          this%eigenvectors(:, num_modes) = eigenvecs(:)
          this%alpha(num_modes) = alpha
        else
          write(message(1), '(A,"Skip zero frequency mode. Omega = ",F15.5)') trim(filename), frequency
          call messages_info(1, namespace=namespace)
          num_zero_modes = num_zero_modes + 1
        end if
      end do
 
      this%num_modes = num_modes
      ! A periodic supercell has the full dim_space*num_atoms modes (of which only the dim_space
      ! acoustic modes at Gamma are zero). For finite systems the bound excludes the translations;
      ! we cannot simply also subtract the dim_space rotational modes, as molecular symmetries
      ! might reduce that.
      if (periodic) then
        max_modes = dim_space*num_atoms
      else
        max_modes = dim_space*(num_atoms-1)
      end if
 
      if (this%num_modes + num_zero_modes > max_modes) then
        message(1) = "The phonon file ("//trim(filename)//") contains too many non-zero frequency modes."
        message(2) = "Most likely the zero frequency modes have not been identified. Try increasing PhononModesZeroThreshold."
        message(3) = "Note that for molecules with rotational symmetry, this warning might be triggered wrongly."
        call messages_warning(3, namespace=namespace)
      end if
 
      ! We can now already precalculate some other quantities, we will use later:
 
      ! l_m = alpha_m * l~_m with l~_m = sqrt(1/(M_0*omega_m)) and M_0 = 1 AMU.
      allocate(this%l_m(1:this%num_modes))
      this%l_m = sqrt(1.0_real64/(unit_amu%factor * this%frequencies(1:this%num_modes))) * this%alpha(1:this%num_modes)
 
      allocate(this%prefactor(1:(dim_space*num_atoms)))
      i = 1
      do iatom = 1, num_atoms
        do idim=1, dim_space
          this%prefactor(i) = 1.0_real64/sqrt(masses(iatom) * this%num_super)
          i=i+1
        end do
      end do
 
 
      call io_close(iunit)
 
    else
 
      this%num_modes = 0
      this%num_super = 1
 
    end if
 
 
    pop_sub(phonon_modes_init)
  end subroutine phonon_modes_init
 
 
  subroutine phonon_modes_get_displacements(this, imode, Q, P, pos_disp_flat, vel_disp_flat)
    class(phonon_modes_t), intent(in)  :: this
    integer,               intent(in)  :: imode
    real(real64),          intent(in)  :: Q
    real(real64),          intent(in)  :: P
    real(real64),          intent(out) :: pos_disp_flat(1:this%dim)
    real(real64),          intent(out) :: vel_disp_flat(1:this%dim)
 
    pos_disp_flat = this%prefactor(:) * this%eigenvectors(:, imode) * q * this%l_m(imode)
    vel_disp_flat = this%prefactor(:) * this%eigenvectors(:, imode) * p * this%l_m(imode) * this%frequencies(imode)
 
  end subroutine phonon_modes_get_displacements
 
 
  !
  subroutine phonon_modes_sample(this, temperature, seed, pos_displacements, vel_displacements)
    class(phonon_modes_t), intent(in)  :: this
    real(real64),          intent(in)  :: temperature
    integer(int64),        intent(in)  :: seed
    real(real64),          intent(out) :: pos_displacements(:, :)
    real(real64),          intent(out) :: vel_displacements(:, :)
 
    type(wigner_distribution_t) :: wigner
    real(real64) :: sigma(1:this%num_modes), mu(1:this%num_modes)
    real(real64) :: Q(1:this%num_modes), P(1:this%num_modes)
    real(real64) :: pos_disp_flat(1:this%dim), vel_disp_flat(1:this%dim)
    real(real64) :: beta_half
    integer :: imode, iatom, dim_space, num_atoms
 
    real(real64), parameter :: low_temperature_tolerance = 1.0e-6_real64
 
    push_sub(phonon_modes_sample)
 
    num_atoms = size(pos_displacements, 2)
    dim_space = size(pos_displacements, 1)
 
    call wigner%init(this%num_modes, seed)
 
    ! The dimensionless mode variables are sampled with variance
    ! sigma^2 = (1 + 2n)/2 = coth(beta*omega/2)/2, which reduces to 1/2 for T -> 0.
    if (temperature < low_temperature_tolerance) then
      sigma = 1.0_real64/sqrt(2.0_real64)
    else
      beta_half = m_one / (2 * p_kb * temperature)
      sigma = 1.0_real64/sqrt(2.0_real64 * tanh(beta_half * this%frequencies(1:this%num_modes)))
    end if
    mu = m_zero
 
 
    q = wigner%get(sigma, mu, wigner_q)
    p = wigner%get(sigma, mu, wigner_p)
 
    pos_displacements = m_zero
    vel_displacements = m_zero
 
    do imode = 1, this%num_modes
      call phonon_modes_get_displacements(this, imode, q(imode), p(imode), pos_disp_flat, vel_disp_flat)
 
      pos_displacements = pos_displacements + reshape(pos_disp_flat, [dim_space, num_atoms])
      vel_displacements = vel_displacements + reshape(vel_disp_flat, [dim_space, num_atoms])
 
    end do
 
 
    call wigner%end()
 
 
    write(message(1), '("Sampled displacements:")')
    call messages_info(1, namespace=this%namespace)
    do iatom = 1, num_atoms
      write(message(1), '(2x,3E15.5)') pos_displacements(1:dim_space, iatom)
      call messages_info(1, namespace=this%namespace)
    end do
 
    write(message(1), '("Sampled velocities:")')
    call messages_info(1, namespace=this%namespace)
    do iatom = 1, num_atoms
      write(message(1), '(2x,3E15.5)') vel_displacements(1:dim_space, iatom)
      call messages_info(1, namespace=this%namespace)
    end do
 
    pop_sub(phonon_modes_sample)
  end subroutine phonon_modes_sample
 
  !
  subroutine next_data_line(iunit, line, filename, namespace, eof)
    integer,           intent(in)  :: iunit
    character(len=*),  intent(out) :: line
    character(len=*),  intent(in)  :: filename
    type(namespace_t), intent(in)  :: namespace
    logical, optional, intent(out) :: eof
 
    integer :: iostat
 
    push_sub(next_data_line)
 
    if (present(eof)) eof = .false.
    do
      read(iunit, '(A)', iostat=iostat) line
      if (iostat /= 0) then
        if (present(eof)) then
          eof = .true.
          line = ""
          pop_sub(next_data_line)
          return
        end if
        message(1) = "Unexpected end of phonon file ("//trim(filename)//")."
        call messages_fatal(1, namespace=namespace)
      end if
      line = adjustl(line)
      if (len_trim(line) == 0) cycle
      if (line(1:1) == '#') cycle
      exit
    end do
 
    pop_sub(next_data_line)
  end subroutine next_data_line
 
  subroutine phonon_modes_finalize(this)
    type(phonon_modes_t), intent(inout) :: this
 
    push_sub(phonon_modes_finalize)
 
    safe_deallocate_a(this%alpha)
 
    this%dim = 0
    this%num_modes = 0
 
    pop_sub(phonon_modes_finalize)
  end subroutine phonon_modes_finalize
 
end module phonon_modes_oct_m