mirror of https://gitlab.com/QEF/q-e.git
1570 lines
51 KiB
Fortran
1570 lines
51 KiB
Fortran
!
|
|
! Copyright (C) 2020-2020 Quantum ESPRESSO group
|
|
!
|
|
! This file is distributed under the terms of the GNU General Public
|
|
! License. See the file `LICENSE' in the root directory of the
|
|
! present distribution, or http://www.gnu.org/copyleft.gpl.txt .
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
PROGRAM postahc
|
|
!------------------------------------------------------------------------------
|
|
!!
|
|
!! This program
|
|
!! - Reads the electron-phonon quantities calculated by ph.x with the
|
|
!! electron_phonon='ahc' option.
|
|
!! - Calculate the phonon-induced electron self-energy in the full matrix
|
|
!! form at a given temperature.
|
|
!!
|
|
!! Input data (namelist "input") is described in Doc/INPUT_POSTAHC.
|
|
!!
|
|
!------------------------------------------------------------------------------
|
|
USE kinds, ONLY : DP
|
|
USE constants, ONLY : ry_to_kelvin, AMU_RY, RY_TO_CMM1, RYTOEV
|
|
USE mp, ONLY : mp_bcast, mp_sum
|
|
USE mp_world, ONLY : world_comm
|
|
USE mp_global, ONLY : mp_startup, mp_global_end
|
|
USE mp_pools, ONLY : npool, nproc_pool, me_pool, intra_pool_comm
|
|
USE io_global, ONLY : ionode_id, ionode, stdout
|
|
USE io_files, ONLY : tmp_dir, prefix
|
|
USE io_global, ONLY : meta_ionode, meta_ionode_id, qestdin
|
|
USE environment, ONLY : environment_start, environment_end
|
|
USE read_namelists_module, ONLY : check_namelist_read
|
|
USE open_close_input_file, ONLY : open_input_file
|
|
USE run_info, ONLY : title
|
|
USE klist, ONLY : nks, xk
|
|
USE symm_base, ONLY : s, invs, nsym
|
|
USE cell_base, ONLY : at, bg
|
|
USE ions_base, ONLY : ntyp => nsp, nat, ityp
|
|
USE wvfct, ONLY : nbnd
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
INTEGER, PARAMETER :: ntypx = 10
|
|
!! Max number of atom types
|
|
!
|
|
! Input variables
|
|
!
|
|
CHARACTER(LEN=256) :: outdir
|
|
!! Directory containing input, output, and scratch files
|
|
LOGICAL :: skip_upper
|
|
!! Skip calculation of upper Fan and Debye-Waller self-energy
|
|
LOGICAL :: skip_dw
|
|
!! Skip calculation of Debye-Waller self-energy
|
|
LOGICAL :: truncate_fan
|
|
!! Truncate Fan self energy to states inside the AHC window
|
|
LOGICAL :: truncate_dw
|
|
!! Truncate Debye-Waller matrix element to states inside the AHC window
|
|
LOGICAL :: adiabatic
|
|
!! Use adiabatic approximation (i.e. ignore the phonon frequency in the denominator)
|
|
LOGICAL :: use_irr_q
|
|
!! Use q points from the irreducible Brillouin zone.
|
|
INTEGER :: ahc_nbnd
|
|
!! Number of bands for which the electron self-energy is to be computed.
|
|
INTEGER :: ahc_nbndskip
|
|
!! Number of bands to exclude when computing the self-energy. The
|
|
!! self-energy is computed for ibnd from ahc_nbndskip + 1
|
|
!! to ahc_nbndskip + ahc_nbnd.
|
|
CHARACTER(LEN=256) :: filename
|
|
!! Name of binary files
|
|
CHARACTER(LEN=256) :: ahc_dir
|
|
!! Directory where the output binary files for AHC e-ph coupling are located
|
|
CHARACTER(LEN=256) :: flvec
|
|
!! output file for normalized phonon displacements generated by matdyn.x
|
|
!! The normalized phonon displacements are the eigenvectors divided by the
|
|
!! square root of the mass then normalized. As such they are not orthogonal.
|
|
REAL(DP) :: ahc_win_min_eV
|
|
!! Lower bound of AHC window for the lower Fan term.
|
|
REAL(DP) :: ahc_win_max_eV
|
|
!! Upper bound of AHC window for the lower Fan term.
|
|
REAL(DP) :: eta_eV
|
|
!! Infinitesimal to add in the denominator for self-energy, in eV
|
|
REAL(DP) :: temp_kelvin
|
|
!! Temperature at which the electron self-energy is calculated, in Kelvins.
|
|
REAL(DP) :: efermi_eV
|
|
!! Fermi level energy, in eV.
|
|
REAL(DP) :: amass_amu(ntypx)
|
|
!! Mass of atoms one for each type, in atomic mass unit.
|
|
!
|
|
! Local variables
|
|
!
|
|
LOGICAL :: lgamma
|
|
!! .true. if q == Gamma
|
|
LOGICAL :: needwf = .FALSE.
|
|
!! do not need to read wavefunction data
|
|
INTEGER :: ios
|
|
!! io status
|
|
INTEGER :: iat
|
|
!! Counter for atoms
|
|
INTEGER :: it
|
|
!! Counter for types
|
|
INTEGER :: ibnd
|
|
!! Counter for bands
|
|
INTEGER :: jbnd
|
|
!! Counter for bands
|
|
INTEGER :: ik
|
|
!! Counter for k points
|
|
INTEGER :: jk
|
|
!! Counter for k points
|
|
INTEGER :: iq
|
|
!! Counter for q points
|
|
INTEGER :: imode
|
|
!! Counter for modes
|
|
INTEGER :: rest
|
|
!! Temporary variable for distributing q points
|
|
INTEGER :: nq_loc
|
|
!! Number of q points in this core
|
|
INTEGER :: iq_from
|
|
!! Index of first q point in this core
|
|
INTEGER :: iq_to
|
|
!! Index of last q point in this core
|
|
INTEGER :: iun
|
|
!! Unit for reading file
|
|
INTEGER :: recl
|
|
!! Record length for reading file
|
|
INTEGER :: count
|
|
!! Counter for degeneracy
|
|
INTEGER :: nmodes
|
|
!! Number of modes. 3 * nat
|
|
INTEGER :: nq
|
|
!! Number of q points
|
|
INTEGER :: nqs
|
|
!! Number of q points in the star
|
|
INTEGER :: ikirr
|
|
!! Index of irreducible k point
|
|
INTEGER :: nkirr
|
|
!! Number of irreducible k points
|
|
INTEGER :: isk
|
|
!! Index of the symmetry-unfolded k point
|
|
INTEGER :: imq
|
|
!! index of -q in the star of q (0 if not present)
|
|
INTEGER :: isq (48)
|
|
!! index of q in the star of a given sym.op
|
|
REAL(DP) :: temperature
|
|
!! temp_kelvin transformed from Kelvin to Ry
|
|
REAL(DP) :: unorm
|
|
!! Norm of u multiplied with amass
|
|
REAL(DP) :: rval
|
|
!! Temporary real variables
|
|
REAL(DP) :: omega_zero_cutoff = 1.d-4
|
|
!! Cutoff of phonon frequency. Modes with omega smaller than
|
|
!! omega_zero_cutoff is neglected.
|
|
REAL(DP) :: e_degen_cutoff = 2.d-5
|
|
!! degeneracy cutoff. Ignore couping between degenerate states with energy
|
|
!! difference less than e_degen_cutoff.
|
|
REAL(DP) :: deltak(3)
|
|
!! k point vector
|
|
REAL(DP) :: sxq(3, 48)
|
|
!! list of vectors in the star of q
|
|
REAL(DP) :: eta
|
|
!! eta_eV converted to Ry units
|
|
REAL(DP) :: efermi
|
|
!! efermi_eV converted to Ry units
|
|
REAL(DP) :: ahc_win_min
|
|
!! ahc_win_min_eV converted to Ry units
|
|
REAL(DP) :: ahc_win_max
|
|
!! ahc_win_max_eV converted to Ry units
|
|
COMPLEX(DP) :: selfen_avg_temp(5)
|
|
!! Diagonal self-energy averaged over degenerate states
|
|
INTEGER, ALLOCATABLE :: ik_to_ikirr(:)
|
|
!! Index of irreducible k point for each of the k point
|
|
REAL(DP), ALLOCATABLE :: inv_omega(:)
|
|
!! (nmodes) 1 / omega
|
|
REAL(DP), ALLOCATABLE :: occph(:)
|
|
!! (nmodes) Bose-Einstein occupation of phonon
|
|
REAL(DP), ALLOCATABLE :: wtq(:)
|
|
!! (nq) Weight of q points. Set to 1/nq.
|
|
REAL(DP), ALLOCATABLE :: amass(:)
|
|
!! Mass of atoms in Ry
|
|
REAL(DP), ALLOCATABLE :: xq(:, :)
|
|
!! (3, nq) q point vectors in Cartesian basis
|
|
REAL(DP), ALLOCATABLE :: omega(:, :)
|
|
!! (nmodes, nq) Phonon frequency
|
|
REAL(DP), ALLOCATABLE :: etk(:, :)
|
|
!! (nbnd, nks) Energy at k
|
|
REAL(DP), ALLOCATABLE :: etk_all(:, :)
|
|
!! (ahc_nbnd, nks) Energy at k
|
|
COMPLEX(DP), ALLOCATABLE :: u(:, :, :)
|
|
!! (nmodes, nmodes, nq) Phonon modes
|
|
COMPLEX(DP), ALLOCATABLE :: ahc_dw(:, :, :, :, :)
|
|
!! Debye-Waller matrix element
|
|
COMPLEX(DP), ALLOCATABLE :: ahc_dw_trunc(:, :, :, :, :)
|
|
!! Debye-Waller matrix element computed with truncation
|
|
COMPLEX(DP), ALLOCATABLE :: selfen_updw(:, :, :)
|
|
!! Debye-Waller self-energy
|
|
COMPLEX(DP), ALLOCATABLE :: selfen_lodw(:, :, :)
|
|
!! Debye-Waller self-energy
|
|
COMPLEX(DP), ALLOCATABLE :: selfen_upfan(:, :, :)
|
|
!! Upper Fan self-energy
|
|
COMPLEX(DP), ALLOCATABLE :: selfen_lofan(:, :, :)
|
|
!! Lower Fan self-energy
|
|
COMPLEX(DP), ALLOCATABLE :: selfen_tot(:, :, :)
|
|
!! Total self-energy
|
|
COMPLEX(DP), ALLOCATABLE :: selfen_diag(:, :)
|
|
!! Diagonal self-energy
|
|
COMPLEX(DP), ALLOCATABLE :: selfen_diag_avg(:, :, :)
|
|
!! Diagonal self-energy averaged over degenerate states
|
|
!
|
|
LOGICAL, EXTERNAL :: imatches
|
|
CHARACTER(LEN=256), EXTERNAL :: trimcheck
|
|
CHARACTER(len=6), EXTERNAL :: int_to_char
|
|
REAL(DP), EXTERNAL :: wgauss
|
|
!
|
|
! ---------------------------------------------------------------------------
|
|
!
|
|
NAMELIST / input / prefix, outdir, skip_upper, skip_dw, ahc_nbnd, &
|
|
ahc_nbndskip, ahc_dir, flvec, eta_eV, temp_kelvin, efermi_eV, amass_amu, &
|
|
ahc_win_min_eV, ahc_win_max_eV, use_irr_q, adiabatic
|
|
!
|
|
CALL mp_startup()
|
|
CALL environment_start('POSTAHC')
|
|
!
|
|
! ---------------------------------------------------------------------------
|
|
!
|
|
! Reading input
|
|
!
|
|
IF (meta_ionode) THEN
|
|
!
|
|
! ... Input from file (ios=0) or standard input (ios=-1) on unit "qestdin"
|
|
!
|
|
ios = open_input_file ( )
|
|
!
|
|
! ... Read the first line of the input file
|
|
!
|
|
IF ( ios <= 0 ) READ( qestdin, '(A)', IOSTAT = ios ) title
|
|
!
|
|
ENDIF
|
|
!
|
|
CALL mp_bcast(ios, meta_ionode_id, world_comm )
|
|
CALL errore( 'postahc', 'reading input file ', ABS( ios ) )
|
|
CALL mp_bcast(title, meta_ionode_id, world_comm )
|
|
!
|
|
! Rewind the input if the title is actually the beginning of inputph namelist
|
|
!
|
|
IF( imatches("&input", title) ) THEN
|
|
WRITE(stdout,'(6x,a)') "Title line not specified: using 'default'."
|
|
title = 'default'
|
|
IF (meta_ionode) REWIND(qestdin, iostat=ios)
|
|
CALL mp_bcast(ios, meta_ionode_id, world_comm )
|
|
CALL errore('postahc', 'Title line missing from input.', abs(ios))
|
|
ENDIF
|
|
!
|
|
! Set default values for variables in namelist
|
|
!
|
|
prefix = 'pwscf'
|
|
CALL get_environment_variable('ESPRESSO_TMPDIR', outdir)
|
|
IF ( TRIM( outdir ) == ' ' ) outdir = './'
|
|
skip_upper = .FALSE.
|
|
skip_dw = .FALSE.
|
|
use_irr_q = .FALSE.
|
|
ahc_nbnd = -1
|
|
ahc_nbndskip = 0
|
|
ahc_dir = ' '
|
|
flvec = ' '
|
|
eta_eV = -9999.d0
|
|
temp_kelvin = -9999.d0
|
|
efermi_eV = -9999.d0
|
|
amass_amu(:) = -9999.d0
|
|
ahc_win_min_eV = -1.d8
|
|
ahc_win_max_eV = 1.d8
|
|
adiabatic = .FALSE.
|
|
!
|
|
! Read the namelist input
|
|
!
|
|
IF (meta_ionode) THEN
|
|
READ( qestdin, INPUT, IOSTAT = ios )
|
|
ENDIF
|
|
CALL check_namelist_read(ios, qestdin, "input")
|
|
!
|
|
CALL mp_bcast(ios, meta_ionode_id, world_comm )
|
|
IF ( ABS(ios) /= 0 ) THEN
|
|
CALL errore( 'postahc', 'reading input namelist', ABS( ios ) )
|
|
ENDIF
|
|
!
|
|
! Broadcast input arguments
|
|
!
|
|
CALL mp_bcast(prefix, ionode_id, world_comm)
|
|
CALL mp_bcast(outdir, ionode_id, world_comm)
|
|
CALL mp_bcast(skip_upper, ionode_id, world_comm)
|
|
CALL mp_bcast(skip_dw, ionode_id, world_comm)
|
|
CALL mp_bcast(use_irr_q, ionode_id, world_comm)
|
|
CALL mp_bcast(ahc_nbnd, ionode_id, world_comm)
|
|
CALL mp_bcast(ahc_nbndskip, ionode_id, world_comm)
|
|
CALL mp_bcast(ahc_dir, ionode_id, world_comm)
|
|
CALL mp_bcast(flvec, ionode_id, world_comm)
|
|
CALL mp_bcast(eta_eV, ionode_id, world_comm)
|
|
CALL mp_bcast(temp_kelvin, ionode_id, world_comm)
|
|
CALL mp_bcast(efermi_eV, ionode_id, world_comm)
|
|
CALL mp_bcast(amass_amu, ionode_id, world_comm)
|
|
CALL mp_bcast(ahc_win_min_eV, ionode_id, world_comm)
|
|
CALL mp_bcast(ahc_win_max_eV, ionode_id, world_comm)
|
|
CALL mp_bcast(adiabatic, ionode_id, world_comm)
|
|
!
|
|
tmp_dir = trimcheck(outdir)
|
|
ahc_dir = trimcheck(ahc_dir)
|
|
!
|
|
truncate_dw = skip_upper
|
|
truncate_fan = skip_upper
|
|
temperature = temp_kelvin / ry_to_kelvin
|
|
ahc_win_min = ahc_win_min_eV / RYTOEV
|
|
ahc_win_max = ahc_win_max_eV / RYTOEV
|
|
eta = eta_eV / RYTOEV
|
|
efermi = efermi_eV / RYTOEV
|
|
!
|
|
! Read xml data. Do not need wavefunction information.
|
|
!
|
|
needwf = .FALSE.
|
|
CALL read_file_new(needwf)
|
|
!
|
|
! Check input argument validity
|
|
!
|
|
IF (ahc_nbnd < 0) CALL errore('postahc', 'ahc_nbnd must be specified', 1)
|
|
IF (ahc_dir == '') CALL errore('postahc', 'ahc_dir must be specified', 1)
|
|
IF (flvec == '') CALL errore('postahc', 'flvec must be specified', 1)
|
|
IF (eta_eV < -9990.d0) CALL errore('postahc', 'eta must be specified', 1)
|
|
IF (efermi_eV < -9990.d0) CALL errore('postahc', 'efermi must be specified', 1)
|
|
IF (temp_kelvin < -9990.d0) CALL errore('postahc', &
|
|
'temp_kelvin must be specified', 1)
|
|
DO it = 1, ntyp
|
|
IF (amass_amu(it) < -9990.d0) CALL errore('postahc', &
|
|
'amass_amu(it) must be specified for it = 1 to ntyp', 1)
|
|
ENDDO
|
|
!
|
|
! Default parallelization is over q points. Pools not implemented.
|
|
!
|
|
IF (npool > 1) CALL errore('postahc', 'pools not implemented', npool)
|
|
!
|
|
IF (truncate_fan .neqv. truncate_dw) THEN
|
|
WRITE(stdout, '(5x, a)') "WARNING: For double-grid calculations, it is strongly advised"
|
|
WRITE(stdout, '(5x, a)') "to set truncate_fan and truncate_dw to the same value because"
|
|
WRITE(stdout, '(5x, a)') "otherwise the double-grid result converge much slowly."
|
|
WRITE(stdout, *) ""
|
|
ENDIF
|
|
!
|
|
IF (adiabatic) THEN
|
|
WRITE(stdout, '(5x, a)') "WARNING: Using the adiabatic approximation. This approximation"
|
|
WRITE(stdout, '(5x, a)') "is inaccurate and even divergent in some materials."
|
|
WRITE(stdout, '(5x, a)') "(See S. Poncé et al., J. Chem. Phys. 143, 102813 (2015) for details.)"
|
|
WRITE(stdout, '(5x, a)') "One should use adiabatic = .TRUE. (default) to get accurate results."
|
|
WRITE(stdout, *) ""
|
|
ENDIF
|
|
!
|
|
! Setup local variables, unit converstion
|
|
!
|
|
nmodes = 3 * nat
|
|
!
|
|
ALLOCATE(amass(nmodes))
|
|
ALLOCATE(occph(nmodes))
|
|
ALLOCATE(inv_omega(nmodes))
|
|
ALLOCATE(etk_all(nbnd, nks))
|
|
ALLOCATE(etk(ahc_nbnd, nks))
|
|
ALLOCATE(ahc_dw(ahc_nbnd, ahc_nbnd, nmodes, 3, nks))
|
|
ALLOCATE(ahc_dw_trunc(ahc_nbnd, ahc_nbnd, nmodes, 3, nks))
|
|
ALLOCATE(selfen_updw(ahc_nbnd, ahc_nbnd, nks))
|
|
ALLOCATE(selfen_lodw(ahc_nbnd, ahc_nbnd, nks))
|
|
ALLOCATE(selfen_upfan(ahc_nbnd, ahc_nbnd, nks))
|
|
ALLOCATE(selfen_lofan(ahc_nbnd, ahc_nbnd, nks))
|
|
ALLOCATE(selfen_tot(ahc_nbnd, ahc_nbnd, nks))
|
|
!
|
|
etk_all = 0.d0
|
|
etk = 0.d0
|
|
ahc_dw = (0.d0, 0.d0)
|
|
ahc_dw_trunc = (0.d0, 0.d0)
|
|
selfen_updw = (0.d0, 0.d0)
|
|
selfen_lodw = (0.d0, 0.d0)
|
|
selfen_upfan = (0.d0, 0.d0)
|
|
selfen_lofan = (0.d0, 0.d0)
|
|
selfen_tot = (0.d0, 0.d0)
|
|
!
|
|
DO iat = 1, nat
|
|
amass((iat-1) * 3 + 1:(iat-1) * 3 + 3) = amass_amu(ityp(iat)) * AMU_RY
|
|
ENDDO
|
|
!
|
|
! ---------------------------------------------------------------------------
|
|
!
|
|
! Calculate the map of the full k points to the irreducible k points
|
|
!
|
|
IF (use_irr_q) CALL setup_ik_to_ikirr_mapping()
|
|
!
|
|
! ---------------------------------------------------------------------------
|
|
!
|
|
! Read flvec file
|
|
!
|
|
IF (ionode) THEN
|
|
CALL read_flvec(flvec, nmodes, nq, xq, omega, u)
|
|
!
|
|
omega = omega / RY_TO_CMM1
|
|
!
|
|
DO iq = 1, nq
|
|
DO imode = 1, nmodes
|
|
unorm = SUM(CONJG(u(:, imode, iq)) * u(:, imode, iq) * amass(:))
|
|
u(:, imode, iq) = u(:, imode, iq) / SQRT(unorm)
|
|
ENDDO
|
|
ENDDO
|
|
ENDIF ! ionode
|
|
!
|
|
CALL mp_bcast(nq, ionode_id, world_comm)
|
|
!
|
|
IF (.NOT. ionode) THEN
|
|
ALLOCATE(xq(3, nq))
|
|
ALLOCATE(omega(nmodes, nq))
|
|
ALLOCATE(u(nmodes, nmodes, nq))
|
|
ENDIF
|
|
!
|
|
CALL mp_bcast(xq, ionode_id, world_comm)
|
|
CALL mp_bcast(omega, ionode_id, world_comm)
|
|
CALL mp_bcast(u, ionode_id, world_comm)
|
|
!
|
|
ALLOCATE(wtq(nq))
|
|
wtq = 1.d0 / REAL(nq, KIND=DP)
|
|
!
|
|
! ---------------------------------------------------------------------------
|
|
!
|
|
! Distribute q points.
|
|
!
|
|
nq_loc = nq / nproc_pool
|
|
rest = nq - nq_loc * nproc_pool
|
|
IF (me_pool < rest) THEN
|
|
nq_loc = nq_loc + 1
|
|
iq_from = me_pool * nq_loc + 1
|
|
iq_to = iq_from + nq_loc - 1
|
|
ELSE
|
|
iq_from = rest * (nq_loc + 1) + (me_pool - rest) * nq_loc + 1
|
|
iq_to = iq_from + nq_loc - 1
|
|
ENDIF
|
|
WRITE(stdout, '(5x, a, I8, a)') "There are ", nq, " q points in total."
|
|
WRITE(stdout, '(5x, a, I8, a)') "Calculating ", nq_loc, " q points in the root core."
|
|
!
|
|
! ---------------------------------------------------------------------------
|
|
!
|
|
! Calculate weight of the q-points for the irreducible grid case.
|
|
!
|
|
IF (use_irr_q) THEN
|
|
IF (ionode) THEN
|
|
DO iq = 1, nq
|
|
CALL star_q(xq(:, iq), at, bg, nsym, s, invs, nqs, sxq, isq, imq, .FALSE.)
|
|
wtq(iq) = REAL(nqs, DP)
|
|
IF (imq == 0) wtq(iq) = REAL(2 * nqs, DP)
|
|
ENDDO
|
|
WRITE(stdout, '(5x,a)') 'Check whether the sum of the weight is equal&
|
|
& to the number of q points in the full q-point grid.'
|
|
WRITE(stdout, '(5x,a,F15.4)') 'SUM(wtq) = ', SUM(wtq)
|
|
wtq = wtq / SUM(wtq)
|
|
ENDIF
|
|
CALL mp_bcast(wtq, ionode_id, world_comm)
|
|
!
|
|
ENDIF ! use_irr_q
|
|
!
|
|
! ---------------------------------------------------------------------------
|
|
!
|
|
WRITE(stdout, '(5x, a, I8)') "Total number of bands : ", nbnd
|
|
WRITE(stdout, '(5x, "Bands to calculate the self-energy :", I8, " to ", I8)') &
|
|
ahc_nbndskip + 1, ahc_nbndskip + ahc_nbnd
|
|
!
|
|
! ---------------------------------------------------------------------------
|
|
!
|
|
! Read ahc_etk_iq1 (All other q points should contain the same data.)
|
|
!
|
|
filename = TRIM(ahc_dir) // 'ahc_etk_iq1.bin'
|
|
!
|
|
INQUIRE(IOLENGTH=recl) etk_all
|
|
OPEN(NEWUNIT=iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', &
|
|
ACCESS='DIRECT', RECL=recl, IOSTAT=ios)
|
|
IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1)
|
|
READ(iun, REC=1) etk_all
|
|
CLOSE(iun)
|
|
!
|
|
! Skip ahc_nbndskip bands in etk
|
|
!
|
|
etk = etk_all(ahc_nbndskip+1:ahc_nbndskip+ahc_nbnd, :)
|
|
!
|
|
! Read ahc_dw which does not depend on iq.
|
|
!
|
|
IF (ionode) THEN
|
|
!
|
|
CALL compute_ahc_dw_with_truncation(ahc_dw_trunc)
|
|
!
|
|
IF (.NOT. truncate_dw) THEN
|
|
filename = TRIM(ahc_dir) // 'ahc_dw.bin'
|
|
INQUIRE(IOLENGTH=recl) ahc_dw
|
|
OPEN(NEWUNIT=iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', &
|
|
ACCESS='DIRECT', RECL=recl, IOSTAT=ios)
|
|
IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1)
|
|
READ(iun, REC=1) ahc_dw
|
|
CLOSE(iun)
|
|
!
|
|
ahc_dw = ahc_dw - ahc_dw_trunc
|
|
!
|
|
ENDIF
|
|
!
|
|
ENDIF ! ionode
|
|
!
|
|
CALL mp_bcast(ahc_dw, ionode_id, world_comm)
|
|
CALL mp_bcast(ahc_dw_trunc, ionode_id, world_comm)
|
|
!
|
|
! ---------------------------------------------------------------------------
|
|
!
|
|
! Loop over q points and calculate self-energy
|
|
!
|
|
WRITE(stdout, *)
|
|
WRITE(stdout,'(5x,a,I8,a)') 'Calculating electron self-energy. Loop over ', &
|
|
iq_to - iq_from + 1, ' q points'
|
|
!
|
|
DO iq = iq_from, iq_to
|
|
!
|
|
! Print progress
|
|
!
|
|
IF (MOD(iq, 10) == 0) THEN
|
|
WRITE(stdout, '(i8)', ADVANCE='no') iq
|
|
IF(MOD(iq, 100) == 0 ) WRITE(stdout,*)
|
|
FLUSH(stdout)
|
|
ENDIF
|
|
!
|
|
lgamma = .FALSE.
|
|
IF ( ALL( ABS(xq(:, iq)) < 1.d-4 ) ) lgamma = .TRUE.
|
|
!
|
|
! Set Bose-Einstein occupation of phonons and set inv_omega
|
|
!
|
|
DO imode = 1, nmodes
|
|
IF (omega(imode, iq) < omega_zero_cutoff) THEN
|
|
inv_omega(imode) = 0.d0
|
|
occph(imode) = 0.d0
|
|
ELSE
|
|
inv_omega(imode) = 1.d0 / omega(imode, iq)
|
|
IF (temperature < 1.d-4) THEN
|
|
occph = 0.d0
|
|
ELSE
|
|
rval = wgauss( omega(imode, iq) / temperature, -99 )
|
|
occph(imode) = 1.d0 / (4.d0 * rval - 2.d0) - 0.5d0
|
|
ENDIF
|
|
ENDIF
|
|
ENDDO
|
|
!
|
|
IF (.NOT. skip_dw) THEN
|
|
CALL calc_debye_waller(iq, selfen_lodw, .TRUE.)
|
|
IF (.NOT. truncate_dw) CALL calc_debye_waller(iq, selfen_updw, .FALSE.)
|
|
ENDIF
|
|
!
|
|
IF (.NOT. truncate_fan) CALL calc_upper_fan(iq, selfen_upfan)
|
|
!
|
|
CALL calc_lower_fan(iq, selfen_lofan)
|
|
!
|
|
ENDDO ! iq
|
|
!
|
|
! ---------------------------------------------------------------------------
|
|
!
|
|
! Collect self-energy from all q points
|
|
!
|
|
CALL mp_sum(selfen_lofan, intra_pool_comm)
|
|
CALL mp_sum(selfen_upfan, intra_pool_comm)
|
|
CALL mp_sum(selfen_lodw, intra_pool_comm)
|
|
CALL mp_sum(selfen_updw, intra_pool_comm)
|
|
!
|
|
! If the energy is outside the AHC window, set all values to zero because
|
|
! the results (in particular the upper Fan self-energy) are meaningless
|
|
DO ik = 1, nks
|
|
DO ibnd = 1, ahc_nbnd
|
|
IF (etk(ibnd, ik) < ahc_win_min .OR. etk(ibnd, ik) > ahc_win_max) THEN
|
|
selfen_lofan(:, ibnd, ik) = (0.d0, 0.d0)
|
|
selfen_lofan(ibnd, :, ik) = (0.d0, 0.d0)
|
|
selfen_upfan(:, ibnd, ik) = 0.d0
|
|
selfen_upfan(ibnd, :, ik) = 0.d0
|
|
selfen_updw(:, ibnd, ik) = 0.d0
|
|
selfen_updw(ibnd, :, ik) = 0.d0
|
|
selfen_lodw(:, ibnd, ik) = 0.d0
|
|
selfen_lodw(ibnd, :, ik) = 0.d0
|
|
ENDIF
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
selfen_tot = selfen_lodw + selfen_updw + selfen_lofan + selfen_upfan
|
|
!
|
|
! Write self-energy to stdout
|
|
!
|
|
IF (ionode) THEN
|
|
!
|
|
WRITE(stdout, *)
|
|
WRITE(stdout, *)
|
|
IF (skip_dw) THEN
|
|
WRITE(stdout, '(5x,a)') 'skip_dw = .true.: Debye-Waller self-energy &
|
|
&is set to zero'
|
|
ENDIF
|
|
IF (skip_upper) THEN
|
|
WRITE(stdout, '(5x,a)') 'skip_upper = .true.: the upper Fan self-energy is set &
|
|
&to zero, and the'
|
|
WRITE(stdout, '(5x,a)') ' Debye-Waller self-energy is truncated.'
|
|
ENDIF
|
|
!
|
|
WRITE(stdout, *)
|
|
WRITE(stdout, '(5x,a)') 'Diagonal electron self-energy in eV'
|
|
WRITE(stdout, '(5x,a)') 'Self-energy of degenerate states are averaged.'
|
|
WRITE(stdout, '(5x,a)') 'The DW and Upper_Fan terms are real-valued by construction.'
|
|
WRITE(stdout, '(5x,a)') 'For states with energy outside the AHC window, all output are zero.'
|
|
WRITE(stdout, '(5x,a)') 'Total = Upper_DW + Lower_DW + Upper_Fan + Lower_Fan'
|
|
WRITE(stdout, *)
|
|
WRITE(stdout, '(5x,a)') 'Begin postahc output'
|
|
WRITE(stdout, '(5x,a)') ' ik ibnd Re[Total] Upper_DW Lower_DW&
|
|
& Upper_Fan Re[Lower_Fan] Im[Total]'
|
|
!
|
|
ALLOCATE(selfen_diag(ahc_nbnd, 5))
|
|
ALLOCATE(selfen_diag_avg(ahc_nbnd, 5, nks))
|
|
!
|
|
DO ik = 1, nks
|
|
DO ibnd = 1, ahc_nbnd
|
|
selfen_diag(ibnd, 1) = selfen_tot(ibnd, ibnd, ik)
|
|
selfen_diag(ibnd, 2) = selfen_updw(ibnd, ibnd, ik)
|
|
selfen_diag(ibnd, 3) = selfen_lodw(ibnd, ibnd, ik)
|
|
selfen_diag(ibnd, 4) = selfen_upfan(ibnd, ibnd, ik)
|
|
selfen_diag(ibnd, 5) = selfen_lofan(ibnd, ibnd, ik)
|
|
ENDDO
|
|
!
|
|
! Average over degenerate states
|
|
!
|
|
DO ibnd = 1, ahc_nbnd
|
|
!
|
|
selfen_avg_temp = (0.d0, 0.d0)
|
|
count = 0
|
|
!
|
|
DO jbnd = 1, ahc_nbnd
|
|
!
|
|
IF (ABS(etk(ibnd, ik) - etk(jbnd, ik)) < e_degen_cutoff) THEN
|
|
count = count + 1
|
|
selfen_avg_temp = selfen_avg_temp + selfen_diag(jbnd, :)
|
|
ENDIF
|
|
!
|
|
ENDDO
|
|
!
|
|
selfen_diag_avg(ibnd, :, ik) = selfen_avg_temp / REAL(count, DP)
|
|
!
|
|
ENDDO ! ibnd
|
|
ENDDO ! ik
|
|
!
|
|
! Write averaged self-energy
|
|
!
|
|
IF (use_irr_q) THEN
|
|
!
|
|
! Average over symmetry-unfolded k points
|
|
!
|
|
DO ikirr = 1, nkirr
|
|
selfen_diag = (0.d0, 0.d0)
|
|
!
|
|
count = 0
|
|
DO ik = 1, nks
|
|
IF (ik_to_ikirr(ik) == ikirr) THEN
|
|
selfen_diag = selfen_diag + selfen_diag_avg(:, :, ik)
|
|
count = count + 1
|
|
ENDIF
|
|
ENDDO ! ik
|
|
!
|
|
selfen_diag = selfen_diag / REAL(count, KIND=DP)
|
|
!
|
|
DO ibnd = 1, ahc_nbnd
|
|
WRITE(stdout, '(5x, 2I6, 6ES13.4)') ikirr, ibnd, &
|
|
REAL(selfen_diag(ibnd, :)) * RYTOEV, &
|
|
AIMAG(selfen_diag(ibnd, 1)) * RYTOEV
|
|
ENDDO
|
|
!
|
|
ENDDO ! ikirr
|
|
!
|
|
ELSE ! .NOT. use_irr_q
|
|
!
|
|
DO ik = 1, nks
|
|
DO ibnd = 1, ahc_nbnd
|
|
WRITE(stdout, '(5x, 2I6, 6ES13.4)') ik, ibnd, &
|
|
REAL(selfen_diag_avg(ibnd, :, ik)) * RYTOEV, &
|
|
AIMAG(selfen_diag_avg(ibnd, 1, ik)) * RYTOEV
|
|
ENDDO
|
|
ENDDO ! ik
|
|
!
|
|
ENDIF ! use_irr_q
|
|
!
|
|
DEALLOCATE(selfen_diag)
|
|
DEALLOCATE(selfen_diag_avg)
|
|
!
|
|
WRITE(stdout, '(5x,a)') 'End postahc output'
|
|
WRITE(stdout, *)
|
|
!
|
|
IF (use_irr_q) THEN
|
|
!
|
|
WRITE(stdout, '(5x,a)') 'Using q points in the irreducible Brillouin zone.'
|
|
WRITE(stdout, '(5x,a)') 'The off-diagonal self-energy matrix is not computed.'
|
|
!
|
|
ELSE ! use_irr_q
|
|
!
|
|
WRITE(stdout, '(5x,a)') 'Full off-diagonal complex self-energy &
|
|
&matrix is written in files'
|
|
WRITE(stdout, '(5x,a)') 'selfen_real.dat and selfen_imag.dat. &
|
|
&These data can differ from'
|
|
WRITE(stdout, '(5x,a)') 'the output above because the self-energy &
|
|
&of degenerate states are'
|
|
WRITE(stdout, '(5x,a)') 'NOT averaged in the selfen_*.dat output.'
|
|
!
|
|
! Write self-energy to selfen.dat
|
|
!
|
|
OPEN(NEWUNIT=iun, FILE='selfen_real.dat', FORM='FORMATTED')
|
|
WRITE(iun, '(a)') '# Real part of full electron self-energy matrix &
|
|
&Re[sigma(ibnd, jbnd, ik)] in eV'
|
|
WRITE(iun, '(a)') '# ik ibnd jbnd Total Upper_DW Lower_DW&
|
|
& Upper_Fan Lower_Fan'
|
|
DO ik = 1, nks
|
|
DO jbnd = 1, ahc_nbnd
|
|
DO ibnd = 1, ahc_nbnd
|
|
WRITE(iun, '(3I6, 5ES16.7)') ik, ibnd, jbnd, &
|
|
REAL(selfen_tot(ibnd, jbnd, ik)) * RYTOEV, &
|
|
REAL(selfen_updw(ibnd, jbnd, ik)) * RYTOEV, &
|
|
REAL(selfen_lodw(ibnd, jbnd, ik)) * RYTOEV, &
|
|
REAL(selfen_upfan(ibnd, jbnd, ik)) * RYTOEV, &
|
|
REAL(selfen_lofan(ibnd, jbnd, ik)) * RYTOEV
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
CLOSE(iun)
|
|
!
|
|
OPEN(NEWUNIT=iun, FILE='selfen_imag.dat', FORM='FORMATTED')
|
|
WRITE(iun, '(a)') '# Imaginary part of full electron self-energy matrix &
|
|
&Im[sigma(ibnd, jbnd, ik)] in eV'
|
|
WRITE(iun, '(a)') '# ik ibnd jbnd Total Upper_DW Lower_DW&
|
|
& Upper_Fan Lower_Fan'
|
|
DO ik = 1, nks
|
|
DO jbnd = 1, ahc_nbnd
|
|
DO ibnd = 1, ahc_nbnd
|
|
WRITE(iun, '(3I6, 5ES16.7)') ik, ibnd, jbnd, &
|
|
AIMAG(selfen_tot(ibnd, jbnd, ik)) * RYTOEV, &
|
|
AIMAG(selfen_updw(ibnd, jbnd, ik)) * RYTOEV, &
|
|
AIMAG(selfen_lodw(ibnd, jbnd, ik)) * RYTOEV, &
|
|
AIMAG(selfen_upfan(ibnd, jbnd, ik)) * RYTOEV, &
|
|
AIMAG(selfen_lofan(ibnd, jbnd, ik)) * RYTOEV
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
CLOSE(iun)
|
|
!
|
|
ENDIF ! use_irr_q
|
|
!
|
|
ENDIF ! ionode
|
|
!
|
|
! ---------------------------------------------------------------------------
|
|
!
|
|
DEALLOCATE(amass)
|
|
DEALLOCATE(wtq)
|
|
DEALLOCATE(xq)
|
|
DEALLOCATE(omega)
|
|
DEALLOCATE(u)
|
|
DEALLOCATE(occph)
|
|
DEALLOCATE(inv_omega)
|
|
DEALLOCATE(etk_all)
|
|
DEALLOCATE(etk)
|
|
DEALLOCATE(ahc_dw)
|
|
DEALLOCATE(ahc_dw_trunc)
|
|
DEALLOCATE(selfen_updw)
|
|
DEALLOCATE(selfen_lodw)
|
|
DEALLOCATE(selfen_upfan)
|
|
DEALLOCATE(selfen_lofan)
|
|
DEALLOCATE(selfen_tot)
|
|
!
|
|
! ---------------------------------------------------------------------------
|
|
!
|
|
WRITE(stdout, *)
|
|
WRITE(stdout, *)
|
|
CALL print_clock('debye_waller')
|
|
CALL print_clock('lower_fan')
|
|
CALL print_clock('upper_fan')
|
|
!
|
|
1001 CONTINUE
|
|
!
|
|
CALL environment_end('POSTAHC')
|
|
CALL mp_global_end()
|
|
!
|
|
CONTAINS
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
INTEGER FUNCTION get_index_of_k_point(xk_want, xk, nks)
|
|
!------------------------------------------------------------------------------
|
|
!! Find give k vector xk_want from the array of k vectors xk.
|
|
!------------------------------------------------------------------------------
|
|
USE kinds, ONLY : DP
|
|
USE cell_base, ONLY : at
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
INTEGER, INTENT(IN) :: nks
|
|
REAL(DP), INTENT(IN) :: xk_want(3)
|
|
REAL(DP), INTENT(IN) :: xk(3, nks)
|
|
!
|
|
INTEGER :: ik
|
|
!! Index of k point
|
|
LOGICAL :: found
|
|
!! True if the k point in the star is found in the k point list
|
|
!
|
|
DO ik = 1, nks
|
|
!
|
|
! Find ik such that xk(:, ik) = xk_want (mod G)
|
|
deltak(:) = xk(:, ik) - xk_want
|
|
!
|
|
CALL cryst_to_cart(1, deltak, at, -1)
|
|
!
|
|
IF (ALL(ABS(deltak - NINT(deltak)) < 1.d-6)) THEN
|
|
get_index_of_k_point = ik
|
|
RETURN
|
|
ENDIF
|
|
!
|
|
ENDDO
|
|
!
|
|
! If the function has not returned, the k vector is not found. Return -1.
|
|
!
|
|
get_index_of_k_point = -1
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
END FUNCTION get_index_of_k_point
|
|
!------------------------------------------------------------------------------
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
SUBROUTINE setup_ik_to_ikirr_mapping()
|
|
!----------------------------------------------------------------------------
|
|
!! When using irreducible q points, for each k point, its symmetry-equivalent
|
|
!! k points must all be included.
|
|
!! Here, find the list of symmetry-inequivalent (irreducible) k points and
|
|
!! find the mapping from the full k points to the irreducible k points.
|
|
!----------------------------------------------------------------------------
|
|
!
|
|
LOGICAL :: missing
|
|
!! Set to true if any of the symmetry pairs is missing
|
|
REAL(DP) :: xk_crys(3)
|
|
!! k point in crystal coordinates
|
|
!
|
|
ALLOCATE(ik_to_ikirr(nks))
|
|
ik_to_ikirr = -1
|
|
ikirr = 0
|
|
!
|
|
missing = .FALSE.
|
|
!
|
|
IF (ionode) THEN
|
|
!
|
|
DO ik = 1, nks
|
|
! Skip if the irreducible k is already found.
|
|
IF (ik_to_ikirr(ik) /= -1) CYCLE
|
|
!
|
|
ikirr = ikirr + 1
|
|
ik_to_ikirr(ik) = ikirr
|
|
!
|
|
CALL star_q(xk(:, ik), at, bg, nsym, s, invs, nqs, sxq, isq, imq, .FALSE.)
|
|
!
|
|
IF (nqs == 1) CYCLE
|
|
!
|
|
! Find other k points in the star
|
|
!
|
|
DO isk = 2, nqs
|
|
jk = get_index_of_k_point(sxq(:, isk), xk, nks)
|
|
!
|
|
IF (jk == -1) THEN
|
|
xk_crys = sxq(:, isk)
|
|
CALL cryst_to_cart(1, xk_crys, at, -1)
|
|
WRITE(stdout, '(5x,a,3F12.8)') 'missing k point (crystal): ', xk_crys
|
|
missing = .TRUE.
|
|
ELSE
|
|
ik_to_ikirr(jk) = ikirr
|
|
ENDIF
|
|
ENDDO ! isk
|
|
!
|
|
IF (imq == 0) THEN
|
|
DO isk = 1, nqs
|
|
jk = get_index_of_k_point(-sxq(:, isk), xk, nks)
|
|
!
|
|
IF (jk == -1) THEN
|
|
xk_crys = -sxq(:, isk)
|
|
CALL cryst_to_cart(1, xk_crys, at, -1)
|
|
WRITE(stdout, '(5x,a,3F12.8)') 'missing k point (crystal): ', xk_crys
|
|
missing = .TRUE.
|
|
ELSE
|
|
ik_to_ikirr(jk) = ikirr
|
|
ENDIF
|
|
ENDDO ! isk
|
|
ENDIF
|
|
!
|
|
ENDDO ! ik
|
|
!
|
|
ENDIF ! ionode
|
|
!
|
|
CALL mp_bcast(missing, ionode_id, world_comm)
|
|
!
|
|
IF (missing) THEN
|
|
WRITE(stdout, *) ''
|
|
WRITE(stdout, '(5x,a)') 'ERROR: k point in the star not found.'
|
|
WRITE(stdout, '(5x,a)') 'To use use_irr_q = .true., the k point for the NSCF calculation list'
|
|
WRITE(stdout, '(5x,a)') 'must contain all k points in the star once and only once.'
|
|
WRITE(stdout, '(5x,a)') 'See above for the list of missing k points in crystal coordinates.'
|
|
WRITE(stdout, '(5x,a)') 'Add these k points to the input of the NSCF calculation and rerun.'
|
|
CALL errore('postahc', 'k point in the star not found.', 1)
|
|
ENDIF
|
|
!
|
|
CALL mp_bcast(ik_to_ikirr, ionode_id, world_comm)
|
|
nkirr = MAXVAL(ik_to_ikirr)
|
|
!
|
|
! Print the list of irreducible k points
|
|
!
|
|
WRITE(stdout, *) ''
|
|
WRITE(stdout, '(5x, a)') 'List of irreducible k points'
|
|
WRITE(stdout, '(5x, a)') 'ik xk(1:3) (cart. coord. in units 2pi/alat)'
|
|
ikirr = 0
|
|
DO ik = 1, nks
|
|
IF (ik_to_ikirr(ik) == ik) THEN
|
|
ikirr = ikirr + 1
|
|
WRITE(stdout, '(I8,3F12.8)') ikirr, xk(:, ik)
|
|
ENDIF
|
|
ENDDO
|
|
WRITE(stdout, *) ''
|
|
!
|
|
END SUBROUTINE setup_ik_to_ikirr_mapping
|
|
!------------------------------------------------------------------------------
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
SUBROUTINE compute_ahc_dw_with_truncation(ahc_dw_trunc)
|
|
!----------------------------------------------------------------------------
|
|
!! Compute ahc_dw_trunc using truncated sum over states.
|
|
!! Include only the states inside the AHC window.
|
|
!------------------------------------------------------------------------------
|
|
!
|
|
COMPLEX(DP), INTENT(INOUT) :: ahc_dw_trunc(ahc_nbnd, ahc_nbnd, nmodes, 3, nks)
|
|
!! Debye-Waller matrix element
|
|
!
|
|
INTEGER :: ik, ib, jb, pb, idir, imode
|
|
CHARACTER(LEN=256) :: filename
|
|
!! Name of files to read
|
|
COMPLEX(DP), ALLOCATABLE :: ahc_gkk(:, :, :, :)
|
|
!! electron-phonon matrix element
|
|
COMPLEX(DP), ALLOCATABLE :: ahc_p(:, :, :, :)
|
|
!! Momentum matrix elements
|
|
!
|
|
ALLOCATE(ahc_gkk(nbnd, ahc_nbnd, nmodes, nks))
|
|
ALLOCATE(ahc_p(nbnd, ahc_nbnd, 3, nks))
|
|
!
|
|
filename = TRIM(ahc_dir) // 'ahc_gkk_iq1.bin'
|
|
CALL postahc_read_unformatted_file(filename, 1, ahc_gkk)
|
|
!
|
|
filename = TRIM(ahc_dir) // 'ahc_p.bin'
|
|
CALL postahc_read_unformatted_file(filename, 1, ahc_p)
|
|
!
|
|
ahc_dw_trunc = (0.d0, 0.d0)
|
|
DO ik = 1, nks
|
|
DO pb = 1, nbnd
|
|
IF (etk_all(pb, ik) >= ahc_win_min .AND. etk_all(pb, ik) <= ahc_win_max) THEN
|
|
DO ib = 1, ahc_nbnd
|
|
DO jb = 1, ahc_nbnd
|
|
DO idir = 1, 3
|
|
DO imode = 1, nmodes
|
|
ahc_dw_trunc(ib, jb, imode, idir, ik) = ahc_dw_trunc(ib, jb, imode, idir, ik) + &
|
|
+ (0.d0, 1.d0) * CONJG(ahc_gkk(pb, ib, imode, ik)) * ahc_p(pb, jb, idir, ik) &
|
|
- (0.d0, 1.d0) * CONJG(ahc_p(pb, ib, idir, ik)) * ahc_gkk(pb, jb, imode, ik)
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ENDIF
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
DEALLOCATE(ahc_gkk)
|
|
DEALLOCATE(ahc_p)
|
|
!
|
|
END SUBROUTINE compute_ahc_dw_with_truncation
|
|
!------------------------------------------------------------------------------
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
SUBROUTINE calc_debye_waller(iq, selfen_dw, truncate)
|
|
!----------------------------------------------------------------------------
|
|
!!
|
|
!! Compute Debye-Waller self-energy at iq
|
|
!!
|
|
!! Implements Eq.(8) of PHonon/Doc/dfpt_self_energy.pdf
|
|
!!
|
|
!! Here, the "operator-generalized acoustic sum rule" is used to represent
|
|
!! Debye-Waller self-energy as a simple matrix element.
|
|
!! See Eq.(13) of the following reference:
|
|
!! Jae-Mo Lihm and Cheol-Hwan Park, Phys. Rev. B, 101, 121102(R) (2020).
|
|
!!
|
|
!----------------------------------------------------------------------------
|
|
USE kinds, ONLY : DP
|
|
!
|
|
LOGICAL, INTENT(IN) :: truncate
|
|
!! If .TRUE., use ahc_dw_truncated, otherwise, use ahc_dw.
|
|
INTEGER, INTENT(IN) :: iq
|
|
COMPLEX(DP), INTENT(INOUT) :: selfen_dw(ahc_nbnd, ahc_nbnd, nks)
|
|
!
|
|
INTEGER :: imode, jmode, kmode
|
|
!! Counter for modes
|
|
INTEGER :: idir, jdir
|
|
!! Counter for directions
|
|
COMPLEX(DP), ALLOCATABLE :: selfen_dw_iq(:, :, :)
|
|
!! Debye-Waller self-energy at iq
|
|
COMPLEX(DP), ALLOCATABLE :: coeff_dw(:, :, :)
|
|
!! Coefficients for Debye-Waller
|
|
!
|
|
CALL start_clock('debye_waller')
|
|
!
|
|
ALLOCATE(selfen_dw_iq(ahc_nbnd, ahc_nbnd, nks))
|
|
ALLOCATE(coeff_dw(nmodes, 3, nmodes))
|
|
coeff_dw = (0.d0, 0.d0)
|
|
selfen_dw_iq = (0.d0, 0.d0)
|
|
!
|
|
! Set coeff_dw:
|
|
! coeff_dw(3*iat + idir, jdir, kmode)
|
|
! = 0.5 * ( CONJG(u(iat*3+idir, kmode)) * u(iat*3+jdir, kmode) ).real
|
|
! * (occph(kmode) + 0.5) * inv_omega(kmode)
|
|
!
|
|
DO kmode = 1, nmodes
|
|
DO jdir = 1, 3
|
|
DO idir = 1, 3
|
|
DO iat = 1, nat
|
|
imode = 3 * (iat - 1) + idir
|
|
jmode = 3 * (iat - 1) + jdir
|
|
coeff_dw(imode, jdir, kmode) = coeff_dw(imode, jdir, kmode) &
|
|
+ 0.5d0 * CONJG(u(imode, kmode, iq)) * u(jmode, kmode, iq) &
|
|
* (occph(kmode) + 0.5d0) * inv_omega(kmode)
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
coeff_dw = REAL(coeff_dw, KIND=DP)
|
|
!
|
|
DO kmode = 1, nmodes
|
|
DO jdir = 1, 3
|
|
DO imode = 1, nmodes
|
|
IF (truncate) THEN
|
|
selfen_dw_iq = selfen_dw_iq + ahc_dw_trunc(:, :, imode, jdir, :) &
|
|
* coeff_dw(imode, jdir, kmode)
|
|
ELSE
|
|
selfen_dw_iq = selfen_dw_iq + ahc_dw(:, :, imode, jdir, :) &
|
|
* coeff_dw(imode, jdir, kmode)
|
|
ENDIF
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
selfen_dw = selfen_dw + selfen_dw_iq * wtq(iq)
|
|
!
|
|
DEALLOCATE(coeff_dw)
|
|
DEALLOCATE(selfen_dw_iq)
|
|
!
|
|
CALL stop_clock('debye_waller')
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
END SUBROUTINE calc_debye_waller
|
|
!------------------------------------------------------------------------------
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
SUBROUTINE calc_upper_fan(iq, selfen_upfan)
|
|
!----------------------------------------------------------------------------
|
|
!!
|
|
!! Compute upper Fan (high-energy band contribution from Sternheimer
|
|
!! equation) self-energy at iq
|
|
!!
|
|
!! Implements Eq.(11) of PHonon/Doc/dfpt_self_energy.pdf
|
|
!!
|
|
!----------------------------------------------------------------------------
|
|
USE kinds, ONLY : DP
|
|
!
|
|
INTEGER, INTENT(IN) :: iq
|
|
COMPLEX(DP), INTENT(INOUT) :: selfen_upfan(ahc_nbnd, ahc_nbnd, nks)
|
|
!
|
|
CHARACTER(LEN=256) :: filename
|
|
!! Name of files to read
|
|
INTEGER :: iun
|
|
!! Unit for reading file
|
|
INTEGER :: recl
|
|
!! Record length for reading file
|
|
INTEGER :: ibnd, jbnd, ibq
|
|
!! Counter for bands
|
|
INTEGER :: imode, jmode, kmode
|
|
!! Counter for modes
|
|
REAL(DP) :: coeff
|
|
!! real coefficient
|
|
REAL(DP), ALLOCATABLE :: etq(:, :)
|
|
!! (nbnd, nks) Energy at k+q
|
|
COMPLEX(DP), ALLOCATABLE :: ahc_gkk(:, :, :, :)
|
|
!! electron-phonon matrix element
|
|
COMPLEX(DP), ALLOCATABLE :: ahc_gkk_mode(:, :, :, :)
|
|
!! electron-phonon matrix element in mode basis
|
|
COMPLEX(DP), ALLOCATABLE :: selfen_upfan_iq(:, :, :)
|
|
!! Upper Fan self-energy at iq
|
|
COMPLEX(DP), ALLOCATABLE :: ahc_upfan(:, :, :, :, :)
|
|
!! Upper Fan matrix element
|
|
COMPLEX(DP), ALLOCATABLE :: ahc_upfan_mode(:, :, :, :)
|
|
!! Upper Fan matrix element in mode basis
|
|
!
|
|
CALL start_clock('upper_fan')
|
|
!
|
|
ALLOCATE(selfen_upfan_iq(ahc_nbnd, ahc_nbnd, nks))
|
|
ALLOCATE(ahc_upfan(ahc_nbnd, ahc_nbnd, nmodes, nmodes, nks))
|
|
ALLOCATE(ahc_upfan_mode(ahc_nbnd, ahc_nbnd, nmodes, nks))
|
|
ALLOCATE(ahc_gkk(nbnd, ahc_nbnd, nmodes, nks))
|
|
ALLOCATE(ahc_gkk_mode(nbnd, ahc_nbnd, nmodes, nks))
|
|
ALLOCATE(etq(nbnd, nks))
|
|
selfen_upfan_iq = (0.d0, 0.d0)
|
|
ahc_upfan_mode = (0.d0, 0.d0)
|
|
ahc_gkk = (0.d0, 0.d0)
|
|
ahc_gkk_mode = (0.d0, 0.d0)
|
|
!
|
|
! Read files: ahc_etq, ahc_gkk, ahc_upfan_iq*.bin
|
|
!
|
|
filename = TRIM(ahc_dir) // 'ahc_gkk_iq' // TRIM(int_to_char(iq)) // '.bin'
|
|
CALL postahc_read_unformatted_file(filename, 1, ahc_gkk)
|
|
!
|
|
filename = TRIM(ahc_dir) // 'ahc_etq_iq' // TRIM(int_to_char(iq)) // '.bin'
|
|
!
|
|
INQUIRE(IOLENGTH=recl) etq
|
|
OPEN(NEWUNIT=iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', &
|
|
ACCESS='DIRECT', RECL=recl, IOSTAT=ios)
|
|
IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1)
|
|
READ(iun, REC=1) etq
|
|
CLOSE(iun)
|
|
!
|
|
filename = TRIM(ahc_dir) // 'ahc_upfan_iq' // TRIM(int_to_char(iq)) // '.bin'
|
|
!
|
|
INQUIRE(IOLENGTH=recl) ahc_upfan
|
|
OPEN(NEWUNIT=iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', &
|
|
ACCESS='DIRECT', RECL=recl, IOSTAT=ios)
|
|
IF (ios /= 0) CALL errore('postahc', 'Error reading ' // TRIM(filename), 1)
|
|
READ(iun, REC=1) ahc_upfan
|
|
CLOSE(iun)
|
|
!
|
|
! rotate ahc_gkk from Cartesian to eigenmode basis
|
|
!
|
|
DO ik = 1, nks
|
|
DO imode = 1, nmodes
|
|
DO jmode = 1, nmodes
|
|
ahc_gkk_mode(:, :, imode, ik) = ahc_gkk_mode(:, :, imode, ik) &
|
|
+ ahc_gkk(:, :, jmode, ik) * u(jmode, imode, iq)
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
! rotate ahc_upfan from Cartesian to eigenmode basis
|
|
!
|
|
DO ik = 1, nks
|
|
DO imode = 1, nmodes
|
|
DO kmode = 1, nmodes
|
|
DO jmode = 1, nmodes
|
|
ahc_upfan_mode(:, :, imode, ik) = ahc_upfan_mode(:, :, imode, ik) &
|
|
+ ahc_upfan(:, :, jmode, kmode, ik) &
|
|
* CONJG(u(jmode, imode, iq)) * u(kmode, imode, iq)
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
! Add contribution of bands outside the AHC wannier window to the upfan matrix.
|
|
!
|
|
DO ik = 1, nks
|
|
DO imode = 1, nmodes
|
|
DO ibnd = 1, ahc_nbnd
|
|
DO jbnd = 1, ahc_nbnd
|
|
!
|
|
DO ibq = 1, nbnd
|
|
!
|
|
! Skip active states inside the ahc window
|
|
IF (etq(ibq, ik) >= ahc_win_min .AND. etq(ibq, ik) <= ahc_win_max) CYCLE
|
|
!
|
|
IF (ABS(etk(ibnd, ik) - etq(ibq, ik)) > 1.d-8) THEN
|
|
ahc_upfan_mode(ibnd, jbnd, imode, ik) &
|
|
= ahc_upfan_mode(ibnd, jbnd, imode, ik) &
|
|
+ CONJG(ahc_gkk_mode(ibq, ibnd, imode, ik)) &
|
|
* ahc_gkk_mode(ibq, jbnd, imode, ik) &
|
|
/ (etk(ibnd, ik) - etq(ibq, ik))
|
|
ENDIF
|
|
!
|
|
ENDDO ! ibq
|
|
!
|
|
ENDDO ! jbnd
|
|
ENDDO ! ibnd
|
|
ENDDO ! imode
|
|
ENDDO ! ik
|
|
!
|
|
! Compute selfen_upfan_iq
|
|
!
|
|
DO ik = 1, nks
|
|
DO imode = 1, nmodes
|
|
!
|
|
coeff = 0.5d0 * inv_omega(imode) * (occph(imode) + 0.5d0)
|
|
!
|
|
DO jbnd = 1, ahc_nbnd
|
|
DO ibnd = 1, ahc_nbnd
|
|
selfen_upfan_iq(ibnd, jbnd, ik) = selfen_upfan_iq(ibnd, jbnd, ik) &
|
|
+ ( ahc_upfan_mode(ibnd, jbnd, imode, ik) &
|
|
+ CONJG(ahc_upfan_mode(jbnd, ibnd, imode, ik)) ) &
|
|
* coeff
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
selfen_upfan = selfen_upfan + selfen_upfan_iq * wtq(iq)
|
|
!
|
|
DEALLOCATE(selfen_upfan_iq)
|
|
DEALLOCATE(ahc_upfan)
|
|
DEALLOCATE(ahc_upfan_mode)
|
|
DEALLOCATE(ahc_gkk_mode)
|
|
DEALLOCATE(ahc_gkk)
|
|
DEALLOCATE(etq)
|
|
!
|
|
CALL stop_clock('upper_fan')
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
END SUBROUTINE calc_upper_fan
|
|
!------------------------------------------------------------------------------
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
SUBROUTINE calc_lower_fan(iq, selfen_lofan)
|
|
!----------------------------------------------------------------------------
|
|
!!
|
|
!! Compute lower Fan (low-energy band contribution) self-energy at iq
|
|
!!
|
|
!! Implements Eq.(10) of PHonon/Doc/dfpt_self_energy.pdf
|
|
!!
|
|
!----------------------------------------------------------------------------
|
|
USE kinds, ONLY : DP
|
|
!
|
|
INTEGER, INTENT(IN) :: iq
|
|
COMPLEX(DP), INTENT(INOUT) :: selfen_lofan(ahc_nbnd, ahc_nbnd, nks)
|
|
!
|
|
CHARACTER(LEN=256) :: filename
|
|
!! Name of ahc_gkk_iq*.bin file
|
|
INTEGER :: iun
|
|
!! Unit for reading file
|
|
INTEGER :: recl
|
|
!! Record length for reading file
|
|
INTEGER :: ibq, ibk, jbk
|
|
!! Counter for bands
|
|
INTEGER :: imode, jmode
|
|
!! Counter for modes
|
|
REAL(DP) :: sign
|
|
!! coefficients
|
|
COMPLEX(DP) :: de1, de2
|
|
!! coefficients
|
|
REAL(DP), ALLOCATABLE :: etq(:, :)
|
|
!! (nbnd, nks) Energy at k+q
|
|
REAL(DP), ALLOCATABLE :: occq(:, :)
|
|
!! (nbnd, nks) Fermi-Dirac occupation of energy at k+q
|
|
COMPLEX(DP), ALLOCATABLE :: selfen_lofan_iq(:, :, :)
|
|
!! Lower Fan self-energy at iq
|
|
COMPLEX(DP), ALLOCATABLE :: ahc_gkk(:, :, :, :)
|
|
!! electron-phonon matrix element
|
|
COMPLEX(DP), ALLOCATABLE :: ahc_gkk_mode(:, :, :, :)
|
|
!! electron-phonon matrix element in mode basis
|
|
COMPLEX(DP), ALLOCATABLE :: coeff(:, :, :, :)
|
|
!! coefficient for lower Fan self-energy
|
|
COMPLEX(DP), ALLOCATABLE :: coeff1(:, :, :, :)
|
|
!! coefficient for lower Fan self-energy
|
|
COMPLEX(DP), ALLOCATABLE :: coeff2(:, :, :, :)
|
|
!! coefficient for lower Fan self-energy
|
|
!
|
|
CALL start_clock('lower_fan')
|
|
!
|
|
ALLOCATE(selfen_lofan_iq(ahc_nbnd, ahc_nbnd, nks))
|
|
ALLOCATE(ahc_gkk(nbnd, ahc_nbnd, nmodes, nks))
|
|
ALLOCATE(ahc_gkk_mode(nbnd, ahc_nbnd, nmodes, nks))
|
|
ALLOCATE(coeff(nbnd, ahc_nbnd, nmodes, nks))
|
|
ALLOCATE(coeff1(nbnd, ahc_nbnd, nmodes, nks))
|
|
ALLOCATE(coeff2(nbnd, ahc_nbnd, nmodes, nks))
|
|
ALLOCATE(etq(nbnd, nks))
|
|
ALLOCATE(occq(nbnd, nks))
|
|
!
|
|
selfen_lofan_iq = (0.d0, 0.d0)
|
|
ahc_gkk = (0.d0, 0.d0)
|
|
ahc_gkk_mode = (0.d0, 0.d0)
|
|
coeff = (0.d0, 0.d0)
|
|
coeff1 = (0.d0, 0.d0)
|
|
coeff2 = (0.d0, 0.d0)
|
|
etq = 0.d0
|
|
occq = 0.d0
|
|
!
|
|
! Read files: ahc_etq, ahc_gkk
|
|
!
|
|
filename = TRIM(ahc_dir) // 'ahc_gkk_iq' // TRIM(int_to_char(iq)) // '.bin'
|
|
CALL postahc_read_unformatted_file(filename, 1, ahc_gkk)
|
|
!
|
|
filename = TRIM(ahc_dir) // 'ahc_etq_iq' // TRIM(int_to_char(iq)) // '.bin'
|
|
!
|
|
INQUIRE(IOLENGTH=recl) etq
|
|
OPEN(NEWUNIT=iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', &
|
|
ACCESS='DIRECT', RECL=recl, IOSTAT=ios)
|
|
IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1)
|
|
READ(iun, REC=1) etq
|
|
CLOSE(iun)
|
|
!
|
|
! Fermi-Dirac occupation at k+q
|
|
!
|
|
DO ik = 1, nks
|
|
DO ibnd = 1, nbnd
|
|
IF (temperature < 1.d-4) THEN
|
|
IF (etq(ibnd, ik) < efermi) THEN
|
|
occq(ibnd, ik) = 1.0d0
|
|
ELSEIF (etq(ibnd, ik) > efermi) THEN
|
|
occq(ibnd, ik) = 0.0d0
|
|
ELSE
|
|
occq(ibnd, ik) = 0.5d0
|
|
ENDIF
|
|
ELSE
|
|
occq(ibnd, ik) = wgauss( (efermi - etq(ibnd, ik)) / temperature, -99 )
|
|
ENDIF
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
! rotate ahc_gkk from Cartesian to eigenmode basis
|
|
!
|
|
DO ik = 1, nks
|
|
DO imode = 1, nmodes
|
|
DO jmode = 1, nmodes
|
|
ahc_gkk_mode(:, :, imode, ik) = ahc_gkk_mode(:, :, imode, ik) &
|
|
+ ahc_gkk(:, :, jmode, ik) * u(jmode, imode, iq)
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
! Compute coefficients
|
|
!
|
|
! sign = +1 if etk(ibk, ik) > efermi
|
|
! = -1 otherwise
|
|
!
|
|
! coeff1(ibq, ibk, imode, ik)
|
|
! = ( 1 - occq(ibq, ik) + occph(imode) )
|
|
! / ( etk(ibk, ik) - etq(ibq, ik) - omega(imode) + 1j * eta * sign )
|
|
!
|
|
! coeff2(ibq, ibk, imode, ik)
|
|
! = ( occq(ibq, ik) + occph(imode) )
|
|
! / ( etk(ibk, ik) - etq(ibq, ik) + omega(imode) + 1j * eta * sign )
|
|
!
|
|
! coeff = (coeff1 + coeff2) * 0.5 * inv_omega(imode)
|
|
!
|
|
coeff1 = 0.d0
|
|
coeff2 = 0.d0
|
|
!
|
|
DO ik = 1, nks
|
|
DO imode = 1, nmodes
|
|
DO ibk = 1, ahc_nbnd
|
|
!
|
|
sign = 1.d0
|
|
IF (etk(ibk, ik) < efermi) sign = -1.d0
|
|
!
|
|
DO ibq = 1, nbnd
|
|
!
|
|
! Skip states outside the AHC window
|
|
IF (etq(ibq, ik) < ahc_win_min .OR. etq(ibq, ik) > ahc_win_max) CYCLE
|
|
!
|
|
IF (adiabatic) THEN
|
|
! Adiabatic approximation: ignore omega in the denominator
|
|
de1 = CMPLX(etk(ibk, ik) - etq(ibq, ik), eta * sign, KIND=DP)
|
|
de2 = CMPLX(etk(ibk, ik) - etq(ibq, ik), eta * sign, KIND=DP)
|
|
ELSE
|
|
de1 = CMPLX(etk(ibk, ik) - etq(ibq, ik) - omega(imode, iq), eta * sign, KIND=DP)
|
|
de2 = CMPLX(etk(ibk, ik) - etq(ibq, ik) + omega(imode, iq), eta * sign, KIND=DP)
|
|
ENDIF
|
|
!
|
|
coeff1(ibq, ibk, imode, ik) = (1.d0 - occq(ibq, ik) + occph(imode)) / de1
|
|
coeff2(ibq, ibk, imode, ik) = ( occq(ibq, ik) + occph(imode) ) / de2
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
coeff = coeff1 + coeff2
|
|
!
|
|
DO ik = 1, nks
|
|
DO imode = 1, nmodes
|
|
coeff(:, :, imode, ik) = coeff(:, :, imode, ik) * 0.5d0 * inv_omega(imode)
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
! Remove coupling with oneself
|
|
!
|
|
IF (lgamma) THEN
|
|
DO ibk = 1, ahc_nbnd
|
|
coeff(ibk + ahc_nbndskip, ibk, :, :) = (0.d0, 0.d0)
|
|
ENDDO
|
|
ENDIF
|
|
!
|
|
! Remove coupling between degenerate states
|
|
!
|
|
IF (lgamma) THEN
|
|
DO ik = 1, nks
|
|
DO ibk = 1, ahc_nbnd
|
|
DO ibq = 1, nbnd
|
|
IF ( ABS(etk(ibk, ik) - etq(ibq, ik)) < e_degen_cutoff ) THEN
|
|
coeff(ibq, ibk, :, ik) = (0.d0, 0.d0)
|
|
ENDIF
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
ENDIF
|
|
!
|
|
! Compute selfen_lofan_iq
|
|
!
|
|
DO ik = 1, nks
|
|
DO imode = 1, nmodes
|
|
!
|
|
DO jbk = 1, ahc_nbnd
|
|
DO ibk = 1, ahc_nbnd
|
|
DO ibq = 1, nbnd
|
|
!
|
|
selfen_lofan_iq(ibk, jbk, ik) = selfen_lofan_iq(ibk, jbk, ik) &
|
|
+ CONJG(ahc_gkk_mode(ibq, ibk, imode, ik)) &
|
|
* ahc_gkk_mode(ibq, jbk, imode, ik) &
|
|
* ( coeff(ibq, ibk, imode, ik) &
|
|
+ coeff(ibq, jbk, imode, ik) ) * 0.5d0
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
selfen_lofan = selfen_lofan + selfen_lofan_iq * wtq(iq)
|
|
!
|
|
DEALLOCATE(coeff)
|
|
DEALLOCATE(coeff1)
|
|
DEALLOCATE(coeff2)
|
|
DEALLOCATE(ahc_gkk)
|
|
DEALLOCATE(ahc_gkk_mode)
|
|
DEALLOCATE(selfen_lofan_iq)
|
|
DEALLOCATE(etq)
|
|
DEALLOCATE(occq)
|
|
!
|
|
CALL stop_clock('lower_fan')
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
END SUBROUTINE calc_lower_fan
|
|
!------------------------------------------------------------------------------
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
SUBROUTINE read_flvec(flvec, nmodes, nq, xq, omega, u)
|
|
!----------------------------------------------------------------------------
|
|
!!
|
|
!! Read flvec (.modes) file generated by matdyn.x
|
|
!!
|
|
!! Output
|
|
!! - nq : number of q points
|
|
!! - xq : list of q point vectors in Cartesian coordinate
|
|
!! - omega : phonon frequency in Ry
|
|
!! - u : phonon modes, normalized to satisfy u^dagger * M * u = identity
|
|
!! (Eq.(1) of PHonon/Doc/dfpt_self_energy.pdf)
|
|
!!
|
|
!----------------------------------------------------------------------------
|
|
USE kinds, ONLY : DP
|
|
!
|
|
CHARACTER(LEN=256), INTENT(IN) :: flvec
|
|
!! Name of the modes file
|
|
INTEGER, INTENT(IN) :: nmodes
|
|
!! Number of modes
|
|
INTEGER, INTENT(OUT) :: nq
|
|
!! Number of q points
|
|
REAL(DP), ALLOCATABLE, INTENT(INOUT) :: xq(:, :)
|
|
!! q point vectors
|
|
REAL(DP), ALLOCATABLE, INTENT(INOUT) :: omega(:, :)
|
|
!! Phonon frequency
|
|
COMPLEX(DP), ALLOCATABLE, INTENT(INOUT) :: u(:, :, :)
|
|
!! Phonon modes
|
|
!
|
|
INTEGER :: i
|
|
!! dummy variable for reading flvec
|
|
REAL(DP) :: omega_
|
|
!! dummy variable for reading flvec
|
|
INTEGER :: iq
|
|
!! Counter for q points
|
|
INTEGER :: imode
|
|
!! Counter for modes
|
|
INTEGER :: iat
|
|
!! Counter for atoms
|
|
INTEGER :: nat
|
|
!! number of atoms. nmodes / 3
|
|
INTEGER :: iun
|
|
!! Unit for reading flvec
|
|
INTEGER :: ios
|
|
!! io status
|
|
!
|
|
nat = nmodes / 3
|
|
!
|
|
! Find out nq from flvec
|
|
!
|
|
OPEN(NEWUNIT=iun, FILE=TRIM(flvec), STATUS='OLD', FORM='FORMATTED', IOSTAT=ios)
|
|
IF (ios /= 0) CALL errore('postahc', 'problem reading flvec file ' // TRIM(flvec), 1)
|
|
!
|
|
nq = 0
|
|
DO WHILE (.TRUE.)
|
|
! flvec has nmodes * (nat + 1) + 5 lines per q point
|
|
DO i = 1, nmodes * (nat + 1) + 5
|
|
READ(iun, '(a)', END=9000)
|
|
ENDDO
|
|
nq = nq + 1
|
|
ENDDO
|
|
9000 CONTINUE
|
|
!
|
|
CLOSE(iun, IOSTAT=ios)
|
|
IF (ios /= 0) CALL errore('postahc', 'problem closing flvec file ' // TRIM(flvec), 1)
|
|
!
|
|
ALLOCATE(xq(3, nq))
|
|
ALLOCATE(omega(nmodes, nq))
|
|
ALLOCATE(u(nmodes, nmodes, nq))
|
|
!
|
|
! Read xq, omega, and u from flvec
|
|
!
|
|
OPEN(NEWUNIT=iun, FILE=TRIM(flvec), STATUS='OLD', FORM='FORMATTED', IOSTAT=ios)
|
|
IF (ios /= 0) CALL errore('postahc', 'problem reading flvec file ' // TRIM(flvec), 1)
|
|
!
|
|
DO iq = 1, nq
|
|
READ(iun, '(a)')
|
|
READ(iun, '(a)')
|
|
READ(iun, '(1x,6x,3f15.8)') (xq(i, iq), i=1, 3)
|
|
READ(iun, '(a)')
|
|
DO imode = 1, nmodes
|
|
READ(iun, 9010) i, omega_, omega(imode, iq)
|
|
DO iat = 1, nat
|
|
READ(iun, 9020) (u(i, imode, iq), i=(iat-1)*3+1, (iat-1)*3+3)
|
|
ENDDO
|
|
ENDDO
|
|
READ(iun, '(a)')
|
|
ENDDO
|
|
9010 format(5x, 6x, i5, 3x, f15.6, 8x, f15.6, 7x)
|
|
9020 format(1x, 1x, 3(f10.6, 1x, f10.6, 3x), 1x)
|
|
!
|
|
CLOSE(iun, IOSTAT=ios)
|
|
IF (ios /= 0) CALL errore('postahc', 'problem closing flvec file ' // TRIM(flvec), 1)
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
END SUBROUTINE read_flvec
|
|
!------------------------------------------------------------------------------
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
SUBROUTINE postahc_read_unformatted_file(filename, irec, array)
|
|
!----------------------------------------------------------------------------
|
|
!! Read a unformatted file to an array.
|
|
!----------------------------------------------------------------------------
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
CHARACTER(LEN=*), INTENT(IN) :: filename
|
|
!! Name of the unformatted file to read
|
|
INTEGER, INTENT(IN) :: irec
|
|
!! Record index to read
|
|
COMPLEX(DP), INTENT(INOUT) :: array(:, :, :, :)
|
|
!! Output. Data read from the file.
|
|
!
|
|
INTEGER :: iun
|
|
!! Unit for reading file
|
|
INTEGER :: recl
|
|
!! Record length for reading file
|
|
INTEGER :: ios
|
|
!! io status
|
|
!
|
|
INQUIRE(IOLENGTH=recl) array
|
|
OPEN(NEWUNIT=iun, FILE=TRIM(filename), STATUS='OLD', FORM='UNFORMATTED', &
|
|
ACCESS='DIRECT', RECL=recl, IOSTAT=ios)
|
|
IF (ios /= 0) CALL errore('postahc', 'Error opening ' // TRIM(filename), 1)
|
|
READ(iun, REC=irec) array
|
|
CLOSE(iun, STATUS='KEEP')
|
|
END SUBROUTINE
|
|
!------------------------------------------------------------------------------
|
|
!
|
|
!------------------------------------------------------------------------------
|
|
END PROGRAM postahc
|
|
!------------------------------------------------------------------------------
|