phonon_modes.F90

Go to the documentation of this file.

!! Copyright (C) 2024 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
 
  implicit none
 
  private
  public :: &
    phonon_modes_t
 
 
  !
  type, extends(classical_modes_t) :: phonon_modes_t
 
    real(real64), allocatable :: alpha(:)
    integer                   :: num_super
 
  contains
 
    procedure :: init => phonon_modes_init 
    final :: phonon_modes_finalize 
 
  end type phonon_modes_t
 
contains
 
  !
  subroutine phonon_modes_init(this, namespace, dim_space, num_atoms, periodic)
    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
 
    character(len=256) :: filename, version_string
    character(len=128) :: dummy
    integer :: iunit
    integer :: imode, num_modes, file_num_modes, num_zero_modes
 
    real(real64), parameter :: frequency_default_threshold = 1.0e-3_real64
    real(real64) :: frequency_threshold
 
    real(real64) :: eigenvecs(1:dim_space*num_atoms)
    real(real64) :: alpha, frequency
 
    push_sub(phonon_modes_init)
 
    this%dim = dim_space * num_atoms
    num_modes = 0
 
    !%Variable PhononModesFile
    !%Type string
    !%Section System
    !%Description
    !% Filename for the phonon modes file. This file is in plain text with the following format:
    !% <tt>
    !% Version: <string> # Version number (currently ignored)
    !% Nmodes: <integer> # Number of modes in the file
    !% Np: <integer>     # Number of primitive cells in the supercell, used to generate the modes
    !% # Comment line
    !% frequency: <float>           # 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 unit cell
    !% frequency: <float>           # Frequency of second mode
    !% ...
    !% </tt>
    !% Note that the order of atoms has to be the same as in the input (or geometry) file.
    !%End
    call parse_variable(namespace, "PhononModesFile", "NONE", filename)
 
    !%Variable PhononModesZeroThreshold
    !%Type float
    !%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.
    !%End
    call parse_variable(namespace, "PhononModesZeroThreshold", frequency_default_threshold, frequency_threshold)
 
    if (trim(filename) /= "NONE") then
 
      iunit = io_open(filename, action="read")
 
      read(iunit, *) version_string
      ! TODO: check whether the version is current.
 
      read(iunit, *) dummy, file_num_modes
      if (trim(dummy) /= 'Nmodes:') then
        message(1) = "Phonon file ("//trim(filename)//") is ill formatted."
        call messages_fatal(1, namespace=namespace)
      end if
 
      read(iunit, *) dummy, this%num_super
      if (trim(dummy) /= 'Np:') then
        message(1) = "Phonon file ("//trim(filename)//") is ill formatted."
        call messages_fatal(1, namespace=namespace)
      end if
 
      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
        read(iunit, *) dummy
        read(iunit, *) 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
          read(iunit, *) 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
          alpha = 1.0
        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
      ! TODO: This should be adjusted for finite systems. Note, however, that we cannot simply add the dim_space rotational
      !       modes, as molecular symmetries might reduce that.
 
      if (this%num_modes + num_zero_modes > dim_space*(num_atoms-1)) then
        message(1) = "The phonon file ("//trim(filename)//") contains too many 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
 
 
      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_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