Branch data Line data Source code
1 : : !!****m* ABINIT/m_tdep_phdos
2 : : !! NAME
3 : : !! m_tdep_phdos
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the routine computing the phonon density of states
7 : : !! and the phonon frequencies and eigenvectors.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2011-2026 ABINIT group (GA,FB,JB)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : module m_tdep_phdos
24 : :
25 : : use defs_basis
26 : : !FB use m_nctk
27 : : use m_errors
28 : : use m_abicore
29 : : use m_phonons
30 : : use m_xmpi
31 : : use m_io_tools
32 : : use m_ifc, only : ifc_type
33 : : use m_crystal, only : crystal_t
34 : : use m_tdep_phi2, only : Eigen_type, tdep_write_yaml, tdep_write_dij
35 : : use m_tdep_qpt, only : Qpoints_type, Qbz_type
36 : : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
37 : :
38 : : implicit none
39 : : !!***
40 : :
41 : : public :: tdep_calc_phdos
42 : :
43 : : contains
44 : :
45 : : !=====================================================================================================
46 : :
47 : : !!****f* ABINIT/m_tdep_phdos/tdep_calc_phdos
48 : : !! NAME
49 : : !! tdep_calc_phdos
50 : : !!
51 : : !! FUNCTION
52 : : !! Compute the phonon density of states (on the q-point grid Invar%ngqpt2, shifted by Invar%q2shft,
53 : : !! with a gaussian smearing) from the IFC, and write it in the _PHDOS file
54 : : !! and, normalized to one, in the _vdos.dat file. The frequency mesh is
55 : : !! enlarged until it contains all the phonon frequencies.
56 : : !! Then, compute the phonon frequencies, eigenvectors and dynamical matrices
57 : : !! along the q-point path (written in the _omega.dat, _eigenvectors.dat,
58 : : !! _dij.dat and YAML files) and on the irreducible q-points of Qbz.
59 : : !!
60 : : !! INPUTS
61 : : !! Invar = Input object containing the input variables.
62 : : !! MPIdata = Info on MPI parallelism.
63 : : !! Ifc = IFC object.
64 : : !! Qbz = Qbz object containing the q-points of the irreducible Brillouin zone.
65 : : !! Qpt = Qpoints object containing the q-point path.
66 : : !!
67 : : !! OUTPUT
68 : : !! PHdos = Phonon density of states (normalized to one).
69 : : !!
70 : : !! SIDE EFFECTS
71 : : !! Crystal = Crystal object describing the unitcell.
72 : : !! Eigen2nd_MP = Eigen object. On exit, contains the phonon frequencies,
73 : : !! eigenvectors and dynamical matrices on the irreducible
74 : : !! q-points of Qbz.
75 : : !! Eigen2nd_path = Eigen object. On exit, contains the phonon frequencies,
76 : : !! eigenvectors and dynamical matrices along the q-point path.
77 : : !!
78 : : !! SOURCE
79 : :
80 [ - + - + : 644 : subroutine tdep_calc_phdos(Invar,MPIdata,Crystal,Ifc,Qbz,Qpt,Eigen2nd_MP,Eigen2nd_path,PHdos)
- + - + -
+ - + - +
- + - + -
+ - + + +
+ + ]
81 : :
82 : : type(atdep_dataset_type),intent(in) :: Invar
83 : : type(phdos_t),intent(out) :: PHdos
84 : : type(ifc_type),intent(in) :: Ifc
85 : : type(crystal_t),intent(inout) :: Crystal
86 : : type(Qbz_type),intent(in) :: Qbz
87 : : type(Qpoints_type),intent(in) :: Qpt
88 : : type(MPI_enreg_type), intent(in) :: MPIdata
89 : : type(Eigen_type),intent(inout) :: Eigen2nd_path
90 : : type(Eigen_type),intent(inout) :: Eigen2nd_MP
91 : :
92 : : integer :: natom_unitcell
93 : : integer :: prtdos,iqpt,iq_ibz,iomega
94 : : integer :: dos_ngqpt(3)
95 : : integer :: count_wminmax(2)
96 : : character (len=fnlen):: phdos_fname
97 : : double precision :: dossmear,integ,domega
98 : : double precision :: dos_qshift(3)
99 [ - + ]: 46 : double precision, allocatable :: displ(:,:)
100 : : character(len=500) :: message
101 : : real(dp) :: wminmax(2)
102 : :
103 : 46 : write(ab_out,*)' '
104 : 46 : write(ab_out,*) '#############################################################################'
105 : 46 : write(ab_out,*) '############## Compute the phonon spectrum and DOS ##########################'
106 : 46 : write(ab_out,*) '#############################################################################'
107 : :
108 : 46 : natom_unitcell = Invar%natom_unitcell
109 : :
110 : : ! Compute the DOS
111 : : ! ===============
112 : 46 : prtdos=1 !Gaussian
113 : : ! prtdos=2 !Tetra
114 : 46 : dossmear=4.5d-6
115 : : ! dossmear=4.5d-5
116 [ + + ]: 184 : dos_qshift(:)=Invar%q2shft(:)
117 [ + + ]: 184 : dos_ngqpt(:)=Invar%ngqpt2(:)
118 : 46 : write(ab_out,'(a)') ' Compute the phonon DOS'
119 : : ! Only 1 shift in q-mesh.
120 : : ! The grid contains Gamma by default (q2shft=0): its three acoustic modes (rigid translations,
121 : : ! at zero frequency) are excluded from the DOS, otherwise their smeared contribution at w~0
122 : : ! makes the mean square displacement and <Omega^(-2)> diverge (see tdep_model_calc_thermo).
123 : 46 : wminmax = zero
124 : 46 : do
125 : : call PHdos%init(Crystal,Ifc,prtdos,Invar%dosdeltae,dossmear,dos_ngqpt,1,dos_qshift, &
126 : 92 : "freq_displ", wminmax, count_wminmax, XMPI_WORLD,prtout=.false.,exclude_gamma_acoustic=.true.)
127 [ + + + + : 229 : if (all(count_wminmax == 0)) exit
+ + ]
128 : 46 : wminmax(1) = wminmax(1) - abs(wminmax(1)) * 0.05
129 : 46 : wminmax(2) = wminmax(2) + abs(wminmax(2)) * 0.05
130 : 46 : call phdos%free()
131 : 46 : write(message, "(a, 2f8.5)")" Initial frequency mesh not large enough. Recomputing PHDOS with wmin, wmax: ",wminmax
132 : 46 : call wrtout(std_out, message)
133 : : end do
134 : :
135 : : ! Print the DOS
136 : : ! =============
137 [ + - ]: 46 : phdos_fname = trim(Invar%output_prefix)//"_PHDOS"
138 [ + - ]: 46 : if (MPIdata%iam_master) then
139 : 46 : write(ab_out,'(a)') ' Writing phonon DOS in file PHDOS'
140 : 46 : call phdos%print(phdos_fname)
141 : : end if
142 : 46 : domega=(Invar%dosdeltae*Ha_meV)
143 : 46 : integ=0.d0
144 [ + + ]: 90332 : do iomega=1,PHdos%nomega
145 : 90332 : integ=integ + domega*PHdos%phdos(iomega)
146 : : end do
147 [ + + ]: 90332 : PHdos%phdos(:)=PHdos%phdos(:)/integ
148 [ + - ]: 46 : if (MPIdata%iam_master) then
149 : 46 : write(ab_out,'(a)') ' Writing vibrational DOS in file vdos.dat'
150 : 46 : open(unit=56,file=trim(Invar%output_prefix)//'_vdos.dat')
151 [ + + ]: 90332 : do iomega=1,PHdos%nomega
152 [ + + ]: 90286 : if (Invar%enunit.eq.0) write(56,'(2(f18.6,1x))') PHdos%omega(iomega)*Ha_eV*1000,PHdos%phdos(iomega)
153 [ + + ]: 90286 : if (Invar%enunit.eq.1) write(56,'(2(f18.6,1x))') PHdos%omega(iomega)*Ha_cmm1 ,PHdos%phdos(iomega)
154 [ + + ]: 90286 : if (Invar%enunit.eq.2) write(56,'(2(f18.6,1x))') PHdos%omega(iomega)*1000 ,PHdos%phdos(iomega)
155 [ + + ]: 90332 : if (Invar%enunit.eq.3) write(56,'(2(f18.6,1x))') PHdos%omega(iomega)*Ha_THz ,PHdos%phdos(iomega)
156 : : end do
157 : 46 : close(56)
158 : : end if
159 : :
160 : : ! Compute the frequencies (2)
161 : : ! =======================
162 [ + - - + : 1372398 : ABI_MALLOC(displ,(2*3*natom_unitcell*3*natom_unitcell,Qpt%nqpt)); displ(:,:)=zero
+ - - + +
- + - + -
- + - + +
+ + + ]
163 [ + + + + : 4783736 : Eigen2nd_path%eigenval=zero ; Eigen2nd_path%eigenvec=zero ; Eigen2nd_path%dynmat=zero
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + ]
164 [ + + ]: 9336 : do iqpt=1,Qpt%nqpt
165 : : call Ifc%fourq(Crystal,Qpt%qpt_red(:,iqpt),Eigen2nd_path%eigenval(:,iqpt),displ(:,iqpt),&
166 : : & out_eigvec=Eigen2nd_path%eigenvec(:,:,:,:,:,iqpt),&
167 : 9336 : & out_d2cart=Eigen2nd_path%dynmat (:,:,:,:,:,iqpt))
168 : : end do
169 [ - + ]: 46 : ABI_FREE(displ)
170 : :
171 : : ! Write the Dij, eigenvalues and eigenvectors (in ASCII and YAML)
172 : : ! ======================================================
173 [ + - ]: 46 : if (MPIdata%iam_master) then
174 : 46 : write(ab_out,'(a)') ' Writing dynamical matrix in file dij.dat'
175 : 46 : write(ab_out,'(a)') ' Writing phonon frequencies in file omega.dat'
176 : 46 : write(ab_out,'(a)') ' Writing phonon eigenvectors in file eigenvectors.dat'
177 : 46 : open(unit=51,file=trim(Invar%output_prefix)//'_eigenvectors.dat')
178 : 46 : open(unit=52,file=trim(Invar%output_prefix)//'_dij.dat')
179 : 46 : open(unit=53,file=trim(Invar%output_prefix)//'_omega.dat')
180 [ + + ]: 46 : if (Invar%enunit.eq.0) write(53,'(a)') '# Phonon frequencies in meV'
181 [ + + ]: 46 : if (Invar%enunit.eq.1) write(53,'(a)') '# Phonon frequencies in cm-1'
182 [ + + ]: 46 : if (Invar%enunit.eq.2) write(53,'(a)') '# Phonon frequencies in mHa'
183 [ + + ]: 46 : if (Invar%enunit.eq.3) write(53,'(a)') '# Phonon frequencies in THz'
184 [ + + ]: 9336 : do iqpt=1,Qpt%nqpt
185 : 9336 : call tdep_write_dij(Eigen2nd_path,iqpt,Invar,Qpt%qpt_red(:,iqpt))
186 : : end do
187 : 46 : close(51)
188 : 46 : close(52)
189 : 46 : close(53)
190 : 46 : call tdep_write_yaml(Eigen2nd_path,Qpt,Invar%output_prefix)
191 : : end if
192 : :
193 : : ! Compute the frequencies (2)
194 : : ! =======================
195 [ - + + - : 253473 : ABI_MALLOC(displ,(2*3*natom_unitcell*3*natom_unitcell,Qbz%nqibz)); displ(:,:)=zero
- + + - +
- + - - +
- + + + +
+ ]
196 [ + + + + : 1658261 : Eigen2nd_MP%eigenval=zero ; Eigen2nd_MP%eigenvec=zero ; Eigen2nd_MP%dynmat=zero
+ + + + +
+ + + + +
+ + + + +
+ + + + +
+ + + + ]
197 [ + + ]: 3135 : do iq_ibz=1,Qbz%nqibz
198 : : call Ifc%fourq(Crystal,Qbz%qibz(:,iq_ibz),Eigen2nd_MP%eigenval(:,iq_ibz),displ(:,iq_ibz),&
199 : : & out_eigvec=Eigen2nd_MP%eigenvec(:,:,:,:,:,iq_ibz),&
200 : 3135 : & out_d2cart=Eigen2nd_MP%dynmat (:,:,:,:,:,iq_ibz))
201 : : end do
202 [ - + ]: 46 : ABI_FREE(displ)
203 : :
204 : 46 : end subroutine tdep_calc_phdos
205 : : !!***
206 : :
207 : : !=====================================================================================================
208 : :
209 : : end module m_tdep_phdos
210 : : !!***
|
