TABLE OF CONTENTS
ABINIT/mkph_linwid [ Functions ]
NAME
mkph_linwid
FUNCTION
Calculate the phonon linewidths on a trajectory in q space
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (MVer) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
Cryst<crystal_t>=Info on the unit cell and symmetries. Ifc<ifc_type>=Object containing the interatomic force constants. elph_ds = datastructure with phonon matrix elements nqpath = dimension of qpath_vertices qpath_vertices = vertices of reciprocal space trajectory
OUTPUT
SIDE EFFECTS
PARENTS
elphon
CHILDREN
ftgam,ftgam_init,gam_mult_displ,ifc_fourq,make_path,phdispl_cart2red wrap2_pmhalf,wrtout,zgemm
SOURCE
37 #if defined HAVE_CONFIG_H 38 #include "config.h" 39 #endif 40 41 #include "abi_common.h" 42 43 subroutine mkph_linwid(Cryst,ifc,elph_ds,nqpath,qpath_vertices) 44 45 use defs_basis 46 use defs_elphon 47 use m_profiling_abi 48 use m_errors 49 use m_bz_mesh 50 use m_ifc 51 52 use m_io_tools, only : open_file 53 use m_geometry, only : phdispl_cart2red 54 use m_dynmat, only : ftgam_init, ftgam 55 use m_crystal, only : crystal_t 56 use m_numeric_tools, only : wrap2_pmhalf 57 58 !This section has been created automatically by the script Abilint (TD). 59 !Do not modify the following lines by hand. 60 #undef ABI_FUNC 61 #define ABI_FUNC 'mkph_linwid' 62 use interfaces_14_hidewrite 63 use interfaces_77_ddb, except_this_one => mkph_linwid 64 !End of the abilint section 65 66 implicit none 67 68 !Arguments ------------------------------------ 69 !scalars 70 integer,intent(in) :: nqpath 71 type(crystal_t),intent(in) :: Cryst 72 type(ifc_type),intent(in) :: ifc 73 type(elph_type),intent(inout) :: elph_ds 74 !arrays 75 real(dp),intent(in) :: qpath_vertices(3,nqpath) 76 77 !Local variables------------------------------- 78 !scalars 79 integer :: ibranch,natom,ii,indx,ipoint,nbranch,nqbz,nsppol,nrpt 80 integer :: isppol,jbranch,qtor,unit_bs,unit_lambda,unit_lwd,npt_tot 81 real(dp) :: diagerr,res 82 character(len=500) :: msg 83 character(len=fnlen) :: fname,base_name 84 !arrays 85 integer :: ndiv(nqpath-1) 86 integer, allocatable :: indxprtqpt(:) 87 real(dp),parameter :: c0(2)=(/0._dp,0._dp/),c1(2)=(/1._dp,0._dp/) 88 real(dp) :: displ_cart(2,3*Cryst%natom,3*Cryst%natom) 89 real(dp) :: displ_red(2,3*Cryst%natom,3*Cryst%natom) 90 real(dp) :: eigval(3*Cryst%natom) 91 real(dp) :: gam_now(2,(3*Cryst%natom)**2) 92 real(dp) :: imeigval(3*Cryst%natom) 93 real(dp) :: lambda(3*Cryst%natom) 94 real(dp) :: pheigvec(2*3*Cryst%natom*3*Cryst%natom),phfrq_tmp(3*Cryst%natom) 95 real(dp) :: qpt(3),redkpt(3) 96 real(dp) :: tmpgam1(2,3*Cryst%natom,3*Cryst%natom) 97 real(dp) :: tmpgam2(2,3*Cryst%natom,3*Cryst%natom) 98 real(dp), allocatable :: coskr(:,:), sinkr(:,:),finepath(:,:) 99 100 ! ********************************************************************* 101 102 DBG_ENTER("COLL") 103 104 natom = Cryst%natom 105 nbranch = elph_ds%nbranch 106 nsppol = elph_ds%nsppol 107 base_name = elph_ds%elph_base_name 108 nrpt = ifc%nrpt 109 110 !=================================================================== 111 !Definition of the q path along which ph linwid will be interpolated 112 !=================================================================== 113 call make_path(nqpath,qpath_vertices,Cryst%gmet,'G',20,ndiv,npt_tot,finepath) 114 ABI_ALLOCATE(indxprtqpt,(npt_tot)) 115 indxprtqpt = 0 116 117 !========================================================== 118 !Open _LWD file and write header 119 !========================================================== 120 fname=trim(base_name) // '_LWD' 121 if (open_file(fname,msg,newunit=unit_lwd,status="unknown") /= 0) then 122 MSG_ERROR(msg) 123 end if 124 125 write (unit_lwd,'(a)') '#' 126 write (unit_lwd,'(a)') '# ABINIT package : Phonon linewidth file' 127 write (unit_lwd,'(a)') '#' 128 write (unit_lwd,'(a,i10,a)') '# Phonon linewidths calculated on ',npt_tot,' points along the qpath' 129 write (unit_lwd,'(a)') '# Description of the Q-path :' 130 write (unit_lwd, '(a,i10)') '# Number of line segments = ',nqpath-1 131 write (unit_lwd,'(a)') '# Vertices of the Q-path and corresponding index = ' 132 133 indx=1 134 indxprtqpt(1) = 1 135 indxprtqpt(npt_tot) = 1 136 137 do ii=1,nqpath 138 write (unit_lwd,'(a,3(e16.6,1x),i8)')'# ',qpath_vertices(:,ii),indx 139 if (ii<nqpath) then 140 indx=indx+ndiv(ii) 141 indxprtqpt(indx) = 1 142 end if 143 end do 144 145 write (unit_lwd,'(a)')'#' 146 147 !========================================================== 148 !Open _BST file and write header 149 !========================================================== 150 fname=trim(base_name) // '_BST' 151 if (open_file(fname,msg,newunit=unit_bs,status="unknown") /= 0) then 152 MSG_ERROR(msg) 153 end if 154 155 write (unit_bs, '(a)') '#' 156 write (unit_bs, '(a)') '# ABINIT package : Phonon band structure file' 157 write (unit_bs, '(a)') '#' 158 write (unit_bs, '(a,i10,a)')'# Phonon BS calculated on ', npt_tot,' points along the qpath' 159 write (unit_bs, '(a,i10)') '# Number of line segments = ', nqpath-1 160 indx=1 161 do ii=1,nqpath 162 write (unit_bs,'(a,3(E16.6,1x),i8)')'# ',qpath_vertices(:,ii),indx 163 if (ii<nqpath) indx=indx+ndiv(ii) 164 end do 165 write (unit_bs,'(a)')'#' 166 167 !MG20060606 168 !========================================================== 169 !open _LAMBDA file and write header 170 !contains \omega(q,n) and \lambda(q,n) and can be plotted using xmgrace 171 !========================================================== 172 fname=trim(base_name) // '_LAMBDA' 173 if (open_file(fname,msg,newunit=unit_lambda,status="unknown") /= 0) then 174 MSG_ERROR(msg) 175 end if 176 177 write (unit_lambda,'(a)') '#' 178 write (unit_lambda,'(a)') '# ABINIT package : Lambda file' 179 write (unit_lambda,'(a)') '#' 180 write (unit_lambda,'(a,i10,a)')'# Lambda(q,nu) calculated on ',npt_tot,' Q-points' 181 write (unit_lambda,'(a)') '# Description of the Q-path :' 182 write (unit_lambda,'(a,i10)') '# Number of line segments = ',nqpath-1 183 write (unit_lambda,'(a)') '# Vertices of the Q-path and corresponding index = ' 184 185 indx=1 186 do ii=1,nqpath 187 write (unit_lambda,'(a,3(E16.6,1x),i8)')'# ',qpath_vertices(:,ii),indx 188 if (ii<nqpath) indx=indx+ndiv(ii) 189 end do 190 write (unit_lambda,'(a)')'#' 191 write (unit_lambda,'(a)')'# index frequency lambda(q,n) frequency lambda(q,n) .... lambda_tot' 192 write (unit_lambda,'(a)')'#' 193 194 !real space to q space 195 qtor=0 196 197 !initialize the maximum phonon frequency 198 elph_ds%omega_min = zero 199 elph_ds%omega_max = zero 200 201 ABI_ALLOCATE(coskr, (npt_tot,nrpt)) 202 ABI_ALLOCATE(sinkr, (npt_tot,nrpt)) 203 call ftgam_init(ifc%gprim, npt_tot, nrpt, finepath, ifc%rpt, coskr, sinkr) 204 205 write (std_out,*) ' mkph_linwid : shape(elph_ds%gamma_qpt) = ',shape(elph_ds%gamma_qpt) 206 nqbz = SIZE(elph_ds%gamma_qpt,DIM=4) 207 write(std_out,*) " nqbz = SIZE(elph_ds%gamma_qpt,DIM=4) = ",nqbz 208 ! 209 !Big do loop over spin polarizations 210 !could put in locally, so phonon stuff is not done twice... 211 ! 212 do isppol=1,nsppol 213 indx=1 214 215 ! Output to the main output file 216 write(msg,'(a,a)')ch10,& 217 & ' Output of the linewidths for the first point of each segment. Linewidths are given in Hartree.' 218 call wrtout(std_out,msg,'COLL') 219 call wrtout(ab_out,msg,'COLL') 220 221 write (std_out,*) ' mkph_linwid : elph_ds%ep_scalprod = ', elph_ds%ep_scalprod 222 223 qtor = 0 224 225 ! Interpolation along specified path in q space 226 do ipoint=1,npt_tot 227 228 ! Get qpoint along the path from qpath_vertices 229 qpt(:) = finepath(:,ipoint) 230 231 call wrap2_pmhalf(qpt(1),redkpt(1),res) 232 call wrap2_pmhalf(qpt(2),redkpt(2),res) 233 call wrap2_pmhalf(qpt(3),redkpt(3),res) 234 qpt(:) = redkpt(:) 235 ! 236 ! This reduced version of ftgkk supposes the kpoints have been integrated 237 ! in integrate_gamma. Do FT from real-space gamma grid to 1 qpt. 238 call ftgam(ifc%wghatm,gam_now,elph_ds%gamma_rpt(:,:,isppol,:),natom,1,ifc%nrpt,qtor, & 239 & coskr(ipoint,:), sinkr(ipoint,:)) 240 ! 241 ! get phonon freqs and eigenvectors anyway 242 ! 243 call ifc_fourq(ifc,cryst,qpt,phfrq_tmp,displ_cart,out_eigvec=pheigvec) 244 ! 245 ! additional frequency factor for some cases 246 ! 247 ! If the matrices do not contain the scalar product with the displ_cart vectors yet do it now 248 if (elph_ds%ep_scalprod == 0) then 249 250 call phdispl_cart2red(natom,Cryst%gprimd,displ_cart,displ_red) 251 252 tmpgam2 = reshape (gam_now, (/2,nbranch,nbranch/)) 253 call gam_mult_displ(nbranch, displ_red, tmpgam2, tmpgam1) 254 255 do jbranch=1,nbranch 256 eigval(jbranch) = tmpgam1(1, jbranch, jbranch) 257 imeigval(jbranch) = tmpgam1(2, jbranch, jbranch) 258 259 if (abs(imeigval(jbranch)) > tol8) then 260 write (msg,'(a,i0,a,es16.8)')' imaginary values for branch = ',jbranch,' imeigval = ',imeigval(jbranch) 261 MSG_WARNING(msg) 262 end if 263 end do 264 265 else if (elph_ds%ep_scalprod == 1) then 266 ! 267 ! Diagonalize gamma matrix at qpoint (complex matrix). 268 ! MJV NOTE: gam_now is recast implicitly here to matrix 269 call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, c1, gam_now, 3*natom,& 270 & pheigvec, 3*natom, c0, tmpgam1, 3*natom) 271 272 call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, c1, pheigvec, 3*natom,& 273 & tmpgam1, 3*natom, c0, tmpgam2, 3*natom) 274 275 diagerr = zero 276 do ibranch=1,nbranch 277 278 eigval(ibranch) = tmpgam2(1,ibranch,ibranch) 279 280 do jbranch=1,ibranch-1 281 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))+abs(tmpgam2(2,jbranch,ibranch)) 282 end do 283 do jbranch=ibranch+1,nbranch 284 diagerr = diagerr + abs(tmpgam2(1,jbranch,ibranch))+abs(tmpgam2(2,jbranch,ibranch)) 285 end do 286 diagerr = diagerr + abs(tmpgam2(2,ibranch,ibranch)) 287 end do 288 289 if (diagerr > tol12) then 290 write (msg,'(a,es14.6)')' Numerical error in diagonalization of gamma with phon eigenvectors: ', diagerr 291 MSG_WARNING(msg) 292 end if 293 294 else 295 write (msg,'(a,i0)')' Wrong value for elph_ds%ep_scalprod = ',elph_ds%ep_scalprod 296 MSG_BUG(msg) 297 end if ! end elph_ds%ep_scalprod if 298 ! 299 ! ========================================================== 300 ! write data to files for each q point 301 ! ========================================================== 302 write (unit_lwd,'(i5)', advance='no') indx 303 do ii=1, nbranch 304 write (unit_lwd,'(E16.5)',advance='no') eigval(ii) 305 end do 306 write (unit_lwd,*) 307 308 ! only print phonon BS for isppol 1: independent of electron spins 309 if (isppol==1) then 310 write (unit_bs,'(i5)', advance='no') indx 311 do ii=1, nbranch 312 write (unit_bs,'(E16.5)',advance='no') phfrq_tmp(ii) 313 end do 314 write (unit_bs,*) 315 end if 316 317 write (unit_lambda,'(i5)', advance='no') indx 318 do ii=1,nbranch 319 lambda(ii)=zero 320 if (abs(phfrq_tmp(ii)) > tol10) lambda(ii)=eigval(ii)/(pi*elph_ds%n0(isppol)*phfrq_tmp(ii)**2) 321 write (unit_lambda,'(es16.8)',advance='no')phfrq_tmp(ii),lambda(ii) 322 end do 323 write (unit_lambda,'(es16.8)',advance='no') sum(lambda) 324 write (unit_lambda,*) 325 326 ! MG NOTE: I wrote a piece of code to output all these quantities using units 327 ! chosen by the user, maybe in version 5.2? 328 ! In this version the output of lambda(q,\nu) has been added 329 330 ! Output to the main output file, for first point in segment 331 if(indxprtqpt(ipoint)==1)then 332 write(msg,'(a,a,3es16.6,a,i4,a,a)')ch10,& 333 & ' Q point =',qpt(:),' isppol = ',isppol,ch10,& 334 & ' Mode number Frequency (Ha) Linewidth (Ha) Lambda(q,n)' 335 call wrtout(std_out,msg,'COLL') 336 call wrtout(ab_out,msg,'COLL') 337 do ii=1,nbranch 338 write(msg,'(i8,es20.6,2es16.6)' )ii,phfrq_tmp(ii),eigval(ii),lambda(ii) 339 call wrtout(std_out,msg,'COLL') 340 call wrtout(ab_out,msg,'COLL') 341 end do 342 end if 343 344 ! find max/min phonon frequency along path chosen 345 ! presumed to be representative of full BZ to within 10 percent 346 elph_ds%omega_min = min(elph_ds%omega_min,1.1_dp*phfrq_tmp(1)) 347 elph_ds%omega_max = max(elph_ds%omega_max,1.1_dp*phfrq_tmp(nbranch)) 348 349 indx = indx+1 350 end do ! end ipoint do 351 352 ! add blank lines to output files between sppol 353 write(msg,'(a)' ) '' 354 call wrtout(unit_lwd,msg,'COLL') 355 call wrtout(unit_lambda,msg,'COLL') 356 call wrtout(std_out,msg,'COLL') 357 call wrtout(ab_out,msg,'COLL') 358 end do ! isppol 359 360 ABI_DEALLOCATE(coskr) 361 ABI_DEALLOCATE(sinkr) 362 363 close(unit=unit_lwd) 364 close(unit=unit_bs) 365 close(unit=unit_lambda) 366 367 ABI_DEALLOCATE(finepath) 368 ABI_DEALLOCATE(indxprtqpt) 369 370 write(std_out,*) ' elph_linwid : omega_min, omega_max = ',elph_ds%omega_min, elph_ds%omega_max 371 372 DBG_EXIT("COLL") 373 374 end subroutine mkph_linwid