TABLE OF CONTENTS
ABINIT/m_thmeig [ Modules ]
NAME
m_thmeig
FUNCTION
Calculate thermal corrections to the eigenvalues.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (PB, XG, GA) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_thmeig 23 24 use defs_basis 25 use m_abicore 26 use m_tetrahedron 27 use m_htetra 28 use m_errors 29 use m_ddb 30 use m_ddb_hdr 31 use m_xmpi 32 use m_sort 33 34 use m_geometry, only : mkrdim, xred2xcart, metric 35 use m_symfind, only : symfind, symlatt 36 use m_symtk, only : mati3inv, matr3inv, symatm 37 use m_crystal, only : crystal_t 38 use m_io_tools, only : open_file 39 use m_dynmat, only : asria_corr, dfpt_phfrq 40 use m_anaddb_dataset, only : anaddb_dataset_type 41 use m_pawtab, only : pawtab_type,pawtab_nullify,pawtab_free 42 use m_kpts, only : getkgrid 43 44 implicit none 45 46 private
m_thmeig/outg2f [ Functions ]
[ Top ] [ m_thmeig ] [ Functions ]
NAME
outg2f
FUNCTION
Output g2f function to file. FIXME: Paul, please explain what g2f is. Probably a variant on the Eliashberg spectral function a2F
INPUTS
OUTPUT
only write
SOURCE
1097 subroutine outg2f(deltaene,enemin,enemax,filnam,g2f,g2fsmear,kpnt,mband,nene,nkpt,nqpt,ntetra,telphint,unit_g2f) 1098 1099 !Arguments ------------------------------------ 1100 !scalars 1101 integer,intent(in) :: mband,nene,nkpt,nqpt,ntetra,telphint,unit_g2f 1102 character(len=fnlen),intent(in) :: filnam 1103 real(dp) :: deltaene,enemin,enemax,g2fsmear 1104 !arrays 1105 real(dp) :: g2f(mband,nkpt,nene),kpnt(3,nkpt,nqpt) 1106 1107 !Local variables------------------------------- 1108 !scalars 1109 integer :: iband,ikpt,iomega,iost 1110 real(dp) :: omega 1111 character(len=fnlen) :: outfile 1112 character(len=500) :: message 1113 !arrays 1114 1115 ! ************************************************************************* 1116 1117 !output the g2f 1118 outfile = trim(filnam) // '_G2F' 1119 open (unit=unit_g2f,file=outfile,status='unknown',iostat=iost) 1120 if (iost /= 0) then 1121 write (message,'(3a)')' thmeig : ERROR- opening file ',trim(outfile),' as new' 1122 ABI_ERROR(message) 1123 end if 1124 1125 write(std_out,*) ' g2f function' 1126 write (unit_g2f,'(a)') '#' 1127 write (unit_g2f,'(a)') '# ABINIT package : g2f file' 1128 write (unit_g2f,'(a)') '#' 1129 write (unit_g2f,'(a,I10)') '# number of qpoints integrated over : ', nqpt 1130 write (unit_g2f,'(a,I10)') '# number of energy points : ', nene 1131 write (unit_g2f,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', enemin, & 1132 & ' Ha and omega_max = ', enemax, ' Ha' 1133 if(telphint==1)then 1134 write (unit_g2f,'(a,E16.6)') '# and the smearing width for gaussians is ', g2fsmear 1135 write (unit_g2f,'(a)') '#' 1136 end if 1137 if(telphint==0)then 1138 write (unit_g2f,'(a,I10)') '# number of tetrahedrons', ntetra 1139 write (unit_g2f,'(a)') '#' 1140 end if 1141 1142 !Write only the a2f function for the first K point 1143 !ikpt=1 1144 do ikpt=1,nkpt 1145 write(unit_g2f,'(a,3es16.8)')' Kpt :', kpnt(:,ikpt,1) 1146 do iband=1,mband 1147 write(unit_g2f,*) 'band :', iband 1148 omega = enemin 1149 do iomega=1,nene 1150 write (unit_g2f,*) omega*Ha_eV*1000, g2f(iband, ikpt,iomega) 1151 omega=omega+deltaene 1152 end do 1153 end do 1154 end do 1155 1156 close (unit=unit_g2f) 1157 1158 end subroutine outg2f
m_thmeig/outphdos [ Functions ]
[ Top ] [ m_thmeig ] [ Functions ]
NAME
outphdos
FUNCTION
Print out phonon density of states
INPUTS
deltaene = step on energy/frequency grid, in Hartree dos_phon = phonon DOS calculated on a grid enemin = minimal frequency enemax = maximal frequency filnam = file name for output to disk g2fsmear = smearing width nene = number of points on energy axis nqpt = number of q-points ntetra = number of tetrahedra, if tetrahedron interpolation is used telphint = flag for el-phonon interpolation method (to indicate Gaussian or tetrahedron integration) unit_phdos = unit for phonon DOS output
OUTPUT
only write
SIDE EFFECTS
NOTES
FIXME overcomplete inputs. Eliminate unit_phdos (just filnam) and deltaene (gotten from max-min/nene)
SOURCE
1010 subroutine outphdos(deltaene,dos_phon,enemin,enemax,filnam,g2fsmear,nene,nqpt,ntetra,telphint,unit_phdos) 1011 1012 !Arguments ------------------------------------ 1013 !scalars 1014 integer,intent(in) :: nene,nqpt,ntetra,telphint,unit_phdos 1015 character(len=fnlen),intent(in) :: filnam 1016 real(dp) :: deltaene,enemin,enemax,g2fsmear 1017 !arrays 1018 real(dp) :: dos_phon(nene) 1019 1020 !Local variables------------------------------- 1021 !scalars 1022 integer :: iomega,iost,step10 1023 real(dp) :: dos_effective,omega 1024 character(len=fnlen) :: outfile 1025 character(len=500) :: message 1026 !arrays 1027 1028 ! ************************************************************************* 1029 1030 outfile = trim(filnam) // '_PDS' 1031 write(message, '(3a)')ch10,& 1032 & ' Will write phonon DOS in file ',trim(outfile) 1033 call wrtout(ab_out,message,'COLL') 1034 call wrtout(std_out,message,'COLL') 1035 1036 write(message, '(4a)')ch10,& 1037 & ' For checking purposes, write ten values in the present file.',ch10,& 1038 & ' Index Energy (in Ha) DOS ' 1039 call wrtout(ab_out,message,'COLL') 1040 call wrtout(std_out,message,'COLL') 1041 1042 open (unit=unit_phdos,file=outfile,status='replace',iostat=iost) 1043 if (iost /= 0) then 1044 write (message,'(3a)')' Opening file ',trim(outfile),' as new' 1045 ABI_ERROR(message) 1046 end if 1047 1048 write (unit_phdos,'(a)') '#' 1049 write (unit_phdos,'(a)') '# ABINIT package : phonon DOS file' 1050 write (unit_phdos,'(a)') '#' 1051 write (unit_phdos,'(a,i10)') '# Number of Qpoints integrated over : ', nqpt 1052 write (unit_phdos,'(a,i10)') '# Number of energy points : ', nene 1053 write (unit_phdos,'(a,es16.6,a,es16.6,a)') '# between omega_min = ', enemin, & 1054 & ' Ha and omega_max = ', enemax, ' Ha' 1055 if(telphint==1)then 1056 write (unit_phdos,'(a,es16.6)') '# The smearing width for gaussians is ', g2fsmear 1057 end if 1058 if(telphint==0)then 1059 write (unit_phdos,'(a,i10)') '# Number of tetrahedrons', ntetra 1060 end if 1061 write (unit_phdos,'(a)') '#' 1062 write (unit_phdos,'(a)') '# Index Energy (in Ha) DOS ' 1063 1064 omega = enemin 1065 do iomega=1,nene 1066 dos_effective=dos_phon(iomega) 1067 if(abs(dos_effective)<tol16)then 1068 dos_effective=zero 1069 end if 1070 step10=nene/10 1071 if(mod(iomega,step10)==1)write (std_out,'(i10,es18.6,es18.6)')iomega, omega, dos_effective 1072 if(mod(iomega,step10)==1)write (ab_out,'(i10,es18.6,es18.6)')iomega, omega, dos_effective 1073 write (unit_phdos, '(i10,es18.6,es18.6)')iomega, omega, dos_effective 1074 omega=omega+deltaene 1075 end do 1076 1077 close (unit=unit_phdos) 1078 1079 end subroutine outphdos
m_thmeig/thmeig [ Functions ]
[ Top ] [ m_thmeig ] [ Functions ]
NAME
thmeig
FUNCTION
This routine calculates the thermal corrections to the eigenvalues. The output is this quantity for the input k point.
INPUTS
elph_base_name = root filename for outputs eig2_filnam = name of the eig2 database file comm=MPI communicator
OUTPUT
SOURCE
72 subroutine thmeig(inp, ddb, crystal, elph_base_name, eig2_filnam, iout, natom, mpert, msize, d2asr, comm) 73 74 !Arguments ------------------------------------ 75 !scalars 76 integer,intent(inout) :: natom 77 integer,intent(in) :: mpert,msize 78 integer,intent(in) :: comm 79 character(len=*),intent(in) :: elph_base_name, eig2_filnam 80 integer,intent(in) :: iout 81 type(crystal_t), intent(inout) :: crystal 82 type(anaddb_dataset_type),intent(inout) :: inp 83 type(ddb_type),intent(inout) :: ddb 84 !arrays 85 real(dp),intent(inout) :: d2asr(2,3,natom,3,natom) 86 87 88 !Local variables------------------------------- 89 !scalars 90 integer,parameter :: msppol=2,master=0,bcorr0=0 91 integer :: msym 92 integer :: nkpt,mband,ntypat 93 integer :: usepaw,natifc 94 integer :: nsym,occopt,nblok2 95 integer :: ntemper,telphint,thmflag 96 integer :: brav,chksymbreak,found,gqpt,iatom1,iatom2,iband,iblok,iblok2,idir1,idir2,ii,ikpt,ilatt,imod,index 97 integer :: iomega,iqpt,iqpt1,iqpt2,iqpt2_previous,iqpt3,iscf_fake,itemper 98 integer :: mpert_eig2,msize2,nene,ng2f,nqshft,nsym_new,unit_g2f,nqpt,nqpt_computed,qptopt,rftyp 99 !integer :: mqpt,nqpt2,option 100 integer :: unit_phdos,unitout 101 integer :: isym 102 integer :: nptsym,use_inversion 103 integer :: ierr 104 real(dp) :: ucvol 105 real(dp) :: g2fsmear,temperinc,tempermin 106 real(dp) :: bosein,deltaene,det,domega,enemax,enemin,fact2i,fact2r,factr 107 real(dp) :: gaussfactor,gaussprefactor,gaussval,invdet,omega,omega_max,omega_min,qnrm,qptrlen 108 real(dp) :: rcvol,tmp,tol,vec1i,vec1r,vec2i,vec2r,veci,vecr,xx 109 real(dp) :: tolsym,tolsym8 !new 110 character(len=500) :: message 111 character(len=fnlen) :: outfile 112 type(ddb_type) :: ddb_eig2 113 type(ddb_hdr_type) :: ddb_hdr 114 !arrays 115 ! FIXME now these must be allocated 116 integer :: ngqpt(9),qptrlatt(3,3),rfelfd(4),rfphon(4),rfstrs(4),vacuum(3) 117 integer :: bravais(11) 118 integer,allocatable :: typat(:),atifc(:) 119 integer,allocatable :: symrel(:,:,:),symrec(:,:,:) 120 integer,allocatable :: indsym(:,:,:) 121 integer,allocatable :: indqpt(:) 122 integer,allocatable :: symafm(:),symafm_new(:) 123 integer,allocatable :: carflg_eig2(:,:,:,:) 124 integer,allocatable :: ptsymrel(:,:,:),symrel_new(:,:,:) 125 !integer,allocatable :: symrec_new(:,:,:) 126 real(dp) :: rprim(3,3),gprim(3,3),rmet(3,3),gmet(3,3) 127 real(dp) :: acell(3) 128 real(dp) :: diff_qpt(3) 129 real(dp) :: gprimd(3,3),mesh(3,3) 130 real(dp) :: qlatt(3,3),qphnrm(3),qpt_search(3,3) 131 real(dp) :: rprimd(3,3),shiftq(3,MAX_NSHIFTK),tempqlatt(3) 132 real(dp) :: dummy(0),dummy2(0,0) 133 real(dp),allocatable :: xcart(:,:),xred(:,:) 134 real(dp),allocatable :: amu(:),zion(:) 135 real(dp),allocatable :: tnons(:,:) 136 real(dp),allocatable :: deigi(:,:), deigr(:,:), multi(:,:), multr(:,:) 137 real(dp),allocatable :: dwtermi(:,:), dwtermr(:,:) 138 real(dp),allocatable :: slope(:,:,:),thmeigen(:,:,:),zeropoint(:,:,:) 139 real(dp),allocatable :: displ(:) 140 real(dp),allocatable :: dos_phon(:),dtweightde(:,:),d2cart(:,:) 141 real(dp),allocatable :: eigvec(:,:,:,:),eigval(:,:),g2f(:,:,:),intweight(:,:,:) 142 real(dp),allocatable :: indtweightde(:,:,:),tmpg2f(:,:,:),tmpphondos(:),total_dos(:),tweight(:,:) 143 real(dp),allocatable :: phfreq(:,:) 144 real(dp),allocatable :: eig2dGamma(:,:,:,:),kpnt(:,:,:) 145 real(dp),allocatable :: dedni(:,:,:,:),dednr(:,:,:,:) 146 real(dp),allocatable :: eigen_in(:) 147 real(dp),allocatable :: qpt_full(:,:),qptnrm(:) 148 real(dp),allocatable :: spqpt(:,:),tnons_new(:,:),spinat(:,:) 149 real(dp),allocatable :: wghtq(:) 150 151 type(t_tetrahedron) :: tetrahedra 152 !type(htetra_t) :: tetrahedra 153 character(len=80) :: errstr 154 155 ! ********************************************************************* 156 157 ! Only master works for the time being 158 if (xmpi_comm_rank(comm) /= master) return 159 160 write(message,'(83a)') ch10,('=',ii=1,80),ch10,& 161 & ' Computation of the electron-phonon changes to the electronic eigenenergies ' 162 call wrtout(ab_out,message,'COLL') 163 call wrtout(std_out,message,'COLL') 164 165 166 !========================================================================= 167 !0) Initializations 168 !========================================================================= 169 170 171 g2fsmear = inp%a2fsmear 172 173 telphint = inp%telphint 174 temperinc = inp%temperinc 175 tempermin = inp%tempermin 176 thmflag = inp%thmflag 177 178 ntemper = inp%ntemper 179 natifc = inp%natifc 180 181 !Open Derivative DataBase then r/w Derivative DataBase preliminary information. 182 183 write(std_out, '(a)' ) '- thmeig: Initialize the second-order electron-phonon file with name :' 184 write(std_out, '(a,a)' )'- ',trim(eig2_filnam) 185 186 call ddb_hdr%open_read(eig2_filnam, xmpi_comm_self) 187 188 mband = ddb_hdr%mband * ddb_hdr%nsppol 189 nkpt = ddb_hdr%nkpt 190 ntypat = ddb_hdr%ntypat 191 192 msym = ddb_hdr%msym 193 nblok2 = ddb_hdr%nblok 194 nsym = ddb_hdr%nsym 195 occopt = ddb%occopt 196 usepaw = ddb_hdr%usepaw 197 198 ABI_MALLOC(typat, (natom)) 199 ABI_MALLOC(atifc, (natom)) 200 ABI_MALLOC(zion, (ntypat)) 201 ABI_MALLOC(amu, (ntypat)) 202 203 ABI_MALLOC(xcart,(3,natom)) 204 ABI_MALLOC(xred,(3,natom)) 205 206 ABI_MALLOC(symafm, (nsym)) 207 ABI_MALLOC(spinat,(3,natom)) 208 209 ABI_MALLOC(symrel, (3,3,nsym)) 210 ABI_MALLOC(symrec, (3,3,nsym)) 211 ABI_MALLOC(tnons, (3,nsym)) 212 ABI_MALLOC(indsym, (4,nsym,natom)) 213 214 ABI_MALLOC(deigi, (mband,nkpt)) 215 ABI_MALLOC(deigr, (mband,nkpt)) 216 ABI_MALLOC(dwtermi, (mband,nkpt)) 217 ABI_MALLOC(dwtermr, (mband,nkpt)) 218 ABI_MALLOC(multi, (mband,nkpt)) 219 ABI_MALLOC(multr, (mband,nkpt)) 220 ABI_MALLOC(slope, (2,mband,nkpt)) 221 ABI_MALLOC(thmeigen, (2,mband,nkpt)) 222 ABI_MALLOC(zeropoint, (2,mband,nkpt)) 223 224 !At present, only atom-type perturbations are allowed for eig2 type matrix elements. 225 mpert_eig2=natom 226 msize2=3*mpert_eig2*3*mpert_eig2 227 228 ddb_eig2%nsppol = ddb_hdr%nsppol 229 call ddb_eig2%malloc(msize2,nblok2,natom,ntypat,mpert_eig2,nkpt,mband) 230 231 ABI_MALLOC(eig2dGamma,(2,msize2,mband,nkpt)) 232 233 ABI_MALLOC(eigvec,(2,3,natom,3*natom)) 234 ABI_MALLOC(phfreq,(3*natom,ddb%nblok)) 235 236 atifc = inp%atifc 237 238 !amu = ddb%amu 239 amu(:) = ddb_hdr%amu(1:ntypat) 240 typat(:) = ddb_hdr%typat(1:natom) 241 zion(:) = ddb_hdr%zion(1:ntypat) 242 symrel(:,:,1:nsym) = ddb_hdr%symrel(:,:,1:nsym) 243 tnons(:,1:nsym) = ddb_hdr%tnons(:,1:nsym) 244 245 xred(:,:) = ddb_hdr%xred(:,:) 246 247 symafm(:) = ddb_hdr%symafm(1:nsym) 248 spinat(:,:) = ddb_hdr%spinat(:,1:natom) 249 250 !symrel = ddb_hdr%symrel ! out 251 !tnons = ddb_hdr%tnons ! out 252 253 !acell = ddb%acell 254 !natom = ddb_hdr%natom 255 acell = ddb_hdr%acell 256 rprim = ddb_hdr%rprim 257 258 !Compute different matrices in real and reciprocal space, also 259 !checks whether ucvol is positive. 260 call mkrdim(acell,rprim,rprimd) 261 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 262 263 !Obtain reciprocal space primitive transl g from inverse trans of r 264 call matr3inv(rprim,gprim) 265 266 !Generate atom positions in cartesian coordinates 267 call xred2xcart(natom,rprimd,xcart,xred) 268 269 !Transposed inversion of the symmetry matrices, for use in 270 !the reciprocal space 271 do isym=1,nsym 272 call mati3inv(symrel(:,:,isym),symrec(:,:,isym)) 273 end do 274 275 !SYMATM generates for all the atoms and all the symmetries, the atom 276 !on which the referenced one is sent and also the translation bringing 277 !back this atom to the referenced unit cell 278 tolsym8=tol8 279 call symatm(indsym,natom,nsym,symrec(:,:,1:nsym),tnons(:,1:nsym),tolsym8,typat,xred) 280 281 !Check the correctness of some input parameters, 282 !and perform small treatment if needed. 283 call chkin9(atifc,natifc,natom) 284 285 eig2dGamma(:,:,:,:)=zero 286 287 ABI_MALLOC(carflg_eig2,(3,mpert_eig2,3,mpert_eig2)) 288 ABI_MALLOC(kpnt,(3,nkpt,1)) 289 290 ! Copy a bunch of stuff back into crystal (to retain old behavior) 291 ! TODO comment these: doesnt make a difference 292 crystal%xcart = xcart 293 crystal%ucvol = ucvol 294 crystal%zion = zion 295 crystal%gmet = gmet 296 crystal%rmet = rmet 297 crystal%nsym = nsym 298 crystal%symrel = symrel 299 crystal%symrec = symrec 300 crystal%tnons = tnons 301 crystal%indsym = indsym 302 303 304 !========================================================================= 305 !1) Take care of the Gamma point for thmflag=3, 5 or 7 306 !========================================================================= 307 308 if(thmflag==3 .or. thmflag==5 .or. thmflag==7) then 309 found=0 310 do iblok2=1,nblok2 311 312 call ddb_eig2%read_d2eig(ddb_hdr, iblok2, iblok2) 313 314 qnrm = ddb_eig2%qpt(1,iblok2)*ddb_eig2%qpt(1,iblok2)+ & 315 & ddb_eig2%qpt(2,iblok2)*ddb_eig2%qpt(2,iblok2)+ & 316 & ddb_eig2%qpt(3,iblok2)*ddb_eig2%qpt(3,iblok2) 317 if(qnrm < DDB_QTOL) then 318 eig2dGamma(:,:,:,:) = ddb_eig2%eig2dval(:,:,:,:) 319 gqpt=iblok2 320 write(std_out,*)'-thmeig: found Gamma point in EIG2 DDB, blok number ',iblok2 321 found=1 322 exit 323 end if 324 end do 325 326 if(found==0)then 327 write(message,'(a,i3,2a)')& 328 & 'Was unable to find the blok for Gamma point in EIG2 DDB file, while thmflag= ',thmflag,ch10,& 329 & 'Action: compute the contribution from Gamma, and merge it in your EIG2 DDB file.' 330 ABI_ERROR(message) 331 end if 332 333 ! Put eig2dGamma in cartesian coordinates 334 call carttransf(ddb_eig2%flg,eig2dGamma,carflg_eig2,gprimd,gqpt,mband,& 335 & mpert_eig2,msize2,natom,nblok2,nkpt,rprimd) 336 337 end if 338 339 call ddb_hdr%close() 340 341 !========================================================================= 342 !2) Calculation of dE(n,k)/dn(Q,j) : consider all q and modes 343 !========================================================================= 344 345 if(thmflag==3 .or. thmflag==4)then 346 347 348 ! Use the first list of q wavevectors 349 nqpt=inp%nph1l 350 ABI_MALLOC(spqpt,(3,nqpt)) 351 do iqpt=1,inp%nph1l 352 spqpt(:,iqpt)=inp%qph1l(:,iqpt)/inp%qnrml1(iqpt) 353 end do 354 ABI_MALLOC(wghtq,(nqpt)) 355 wghtq(:)=one/nqpt 356 357 else if(thmflag>=5 .and. thmflag<=8)then 358 359 ! Generates the q point grid 360 ngqpt(1:3)=inp%ngqpt(1:3) 361 nqshft=inp%nqshft 362 qptrlatt(:,:)=0 363 qptrlatt(1,1)=ngqpt(1) 364 qptrlatt(2,2)=ngqpt(2) 365 qptrlatt(3,3)=ngqpt(3) 366 367 ABI_MALLOC(ptsymrel,(3,3,msym)) 368 ABI_MALLOC(symafm_new,(msym)) 369 ABI_MALLOC(symrel_new,(3,3,msym)) 370 ABI_MALLOC(tnons_new,(3,msym)) 371 if(thmflag==7 .or. thmflag==8) then 372 ! Re-generate symmetry operations from the lattice and atomic coordinates 373 tolsym=tol8 374 call symlatt(bravais,std_out,msym,nptsym,ptsymrel,rprimd,tolsym) 375 use_inversion=1 376 call symfind(gprimd,msym,natom,nptsym,1,nsym_new,0,& 377 & ptsymrel,spinat,symafm_new,symrel_new,tnons_new,tolsym,typat,use_inversion,xred) 378 write(std_out,*)' thmeig : found ',nsym_new,' symmetries ',ch10 379 qptopt=1 380 else 381 nsym_new=1 382 symrel_new(:,:,1)=0 ; symrel_new(1,1,1)=1 ; symrel_new(2,2,1)=1 ; symrel_new(3,3,1)=1 383 tnons_new(:,1)=zero 384 symafm_new(1)=1 385 qptopt=3 386 end if 387 388 ! GA: This is useless 389 brav=inp%brav 390 if(abs(brav)/=1)then 391 message = ' The possibility to have abs(brav)/=1 for thmeig was disabled.' 392 ABI_ERROR(message) 393 end if 394 395 ! Prepare to compute the q-point grid in the ZB or IZB 396 iscf_fake=5 ! Need the weights 397 chksymbreak=0 398 vacuum=0 399 shiftq(:,1:nqshft)=inp%q1shft(:,1:nqshft) 400 ! Compute the final number of q points 401 call getkgrid(chksymbreak,0,iscf_fake,dummy2,qptopt,qptrlatt,qptrlen,& 402 & nsym_new,0,nqpt,nqshft,nsym_new,rprimd,& 403 & shiftq,symafm_new,symrel_new,vacuum,dummy) 404 ABI_MALLOC(spqpt,(3,nqpt)) 405 ABI_MALLOC(wghtq,(nqpt)) 406 call getkgrid(chksymbreak,iout,iscf_fake,spqpt,qptopt,qptrlatt,qptrlen,& 407 & nsym_new,nqpt,nqpt_computed,nqshft,nsym_new,rprimd,& 408 & shiftq,symafm_new,symrel_new,vacuum,wghtq) 409 410 ABI_FREE(ptsymrel) 411 ABI_FREE(symafm_new) 412 ABI_FREE(symrel_new) 413 ABI_FREE(tnons_new) 414 415 end if 416 417 call ddb_hdr%free() 418 419 420 write(message,'(a,a)')ch10,' thmeig : list of q wavevectors, with integration weights ' 421 call wrtout(ab_out,message,'COLL') 422 call wrtout(std_out,message,'COLL') 423 do iqpt=1,nqpt 424 write(message,'(i6,3es16.6,es20.6)')iqpt,spqpt(:,iqpt),wghtq(iqpt) 425 call wrtout(ab_out,message,'COLL') 426 call wrtout(std_out,message,'COLL') 427 end do 428 429 if(.not.allocated(indqpt))then 430 ABI_MALLOC(indqpt,(nqpt)) 431 end if 432 ABI_MALLOC(dedni,(mband,nkpt,3*natom,nqpt)) 433 ABI_MALLOC(dednr,(mband,nkpt,3*natom,nqpt)) 434 ABI_MALLOC(eigen_in,(nqpt)) 435 ABI_MALLOC(qpt_full,(3,nqpt)) 436 ABI_MALLOC(qptnrm,(nqpt)) 437 438 dednr(:,:,:,:) = zero 439 dedni(:,:,:,:) = zero 440 441 !!Prepare the reading of the EIG2 files 442 call ddb_hdr%open_read(eig2_filnam, xmpi_comm_self, msym=msym) 443 call ddb_hdr%free() 444 445 !iqpt2 will be the index of the q point bloks inside the EIG2 file 446 iqpt2=0 447 448 !Sum on all phonon wavevectors and modes 449 do iqpt=1,nqpt 450 451 ! Finding the target wavevector in DDB file 452 qpt_search(:,:)=0.0d0 453 qpt_search(:,1)=spqpt(:,iqpt) 454 qphnrm(:)=one 455 rfphon(1:2)=1 456 ! NOTE : at present, no LO-TO splitting included !!! 457 rfelfd(1:2)=0 458 rfstrs(1:2)=0 459 rftyp=1 460 461 write(std_out,'(a,3es16.6)' )' Looking for spqpt=',qpt_search(:,1) 462 463 call ddb%get_block(iblok,qpt_search,qphnrm,rfphon,rfelfd,rfstrs,rftyp) 464 465 if(iblok==0) then 466 write(message,'(a,3es16.6,2a)')& 467 & 'Was unable to find in DDB file, the blok for point ',spqpt(:,iqpt),ch10,& 468 & 'Action: compute the contribution from this point, and merge it in your DDB file.' 469 ABI_ERROR(message) 470 end if 471 472 ABI_MALLOC(d2cart,(2,msize)) 473 ! Copy the dynamical matrix in d2cart 474 d2cart(:,1:msize)=ddb%val(:,:,iblok) 475 476 ! Eventually impose the acoustic sum rule based on previously calculated d2asr 477 !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart) 478 if (inp%asr==1 .or. inp%asr==2 .or. inp%asr==5) then 479 call asria_corr(inp%asr,d2asr,d2cart,mpert,natom) 480 end if 481 482 ! Calculation of the eigenvectors and eigenvalues 483 ! of the dynamical matrix 484 ABI_MALLOC(displ,(2*3*natom*3*natom)) 485 ABI_MALLOC(eigval,(3,natom)) 486 call dfpt_phfrq(amu,displ,d2cart,eigval,eigvec,indsym,& 487 & mpert,msym,natom,nsym,ntypat,phfreq(:,iqpt),qphnrm(1),spqpt(:,iqpt),rprimd,inp%symdynmat,& 488 & symrel,symafm,typat,ucvol) 489 ABI_FREE(displ) 490 ABI_FREE(eigval) 491 ABI_FREE(d2cart) 492 493 494 ! Read the next bloks to find the next q point. 495 found=0 ; iqpt2_previous=iqpt2 496 do while (iqpt2<nblok2) 497 iqpt2=iqpt2+1 498 call ddb_eig2%read_d2eig(ddb_hdr, iqpt2, iqpt2) 499 diff_qpt(:)=ddb_eig2%qpt(1:3,iqpt2)/ddb_eig2%nrm(1,iqpt2)-spqpt(:,iqpt) 500 if(diff_qpt(1)**2+diff_qpt(2)**2+diff_qpt(3)**2 < DDB_QTOL )then 501 found=1 502 exit 503 end if 504 end do 505 506 ! Usually, the q points come in the right order. However, this is not always the case... 507 if(found==0)then 508 509 ! If the EIG2 database file has to be read again, close it, then search for the right q point, 510 ! from the beginning of the file 511 call ddb_hdr%close() 512 513 call ddb_hdr%open_read(eig2_filnam, xmpi_comm_self, msym=msym) 514 515 ! And examine again the EIG2 file. Still, not beyond the previously examined value. 516 found=0 517 do iqpt2=1,iqpt2_previous 518 call ddb_eig2%read_d2eig(ddb_hdr, iqpt2, iqpt2) 519 diff_qpt(:)=ddb_eig2%qpt(1:3,iqpt2)/ddb_eig2%nrm(1,iqpt2)-spqpt(:,iqpt) 520 if(diff_qpt(1)**2+diff_qpt(2)**2+diff_qpt(3)**2 < DDB_QTOL )then 521 found=1 522 exit 523 end if 524 end do 525 526 if(found==0)then 527 write(message,'(a,3es16.6,2a)')& 528 & 'Was unable to find in EIG2 DDB file, the blok for point ',spqpt(:,iqpt),ch10,& 529 & 'Action: compute the contribution from this point, and merge it in your EIG2 DDB file.' 530 ABI_ERROR(message) 531 end if 532 533 call ddb_hdr%free() 534 535 end if 536 537 ! Put eig2dval in cartesian coordinates 538 call carttransf(ddb_eig2%flg,ddb_eig2%eig2dval,carflg_eig2,gprimd,iqpt,mband,& 539 & mpert_eig2,msize2,natom,nblok2,nkpt,rprimd) 540 541 do imod=1,3*natom 542 543 ! Calculate the derivative 544 deigr(:,:) = zero 545 deigi(:,:) = zero 546 dwtermr(:,:)=zero 547 dwtermi(:,:)=zero 548 index=0 549 do iatom1=1,natom 550 do idir1=1,3 551 do iatom2=1,natom 552 ! Compute factor for SE term 553 if(phfreq(imod,iqpt)<tol6)then 554 factr = zero 555 else 556 factr=one/sqrt(amu(typat(iatom1))*amu(typat(iatom2)))/phfreq(imod,iqpt)/amu_emass 557 end if 558 559 do idir2=1,3 560 index = idir1 + 3*((iatom1 - 1) + natom * ((idir2-1)+3*(iatom2-1))) 561 562 ! Compute products of polarization vectors 563 vecr = eigvec(1,idir1,iatom1,imod)*eigvec(1,idir2,iatom2,imod)+& 564 & eigvec(2,idir1,iatom1,imod)*eigvec(2,idir2,iatom2,imod) 565 veci = eigvec(2,idir1,iatom1,imod)*eigvec(1,idir2,iatom2,imod)-& 566 & eigvec(1,idir1,iatom1,imod)*eigvec(2,idir2,iatom2,imod) 567 568 vec1r = eigvec(1,idir1,iatom1,imod)*eigvec(1,idir2,iatom1,imod)+& 569 & eigvec(2,idir1,iatom1,imod)*eigvec(2,idir2,iatom1,imod) 570 vec1i = eigvec(2,idir1,iatom1,imod)*eigvec(1,idir2,iatom1,imod)-& 571 & eigvec(1,idir1,iatom1,imod)*eigvec(2,idir2,iatom1,imod) 572 573 vec2r = eigvec(1,idir1,iatom2,imod)*eigvec(1,idir2,iatom2,imod)+& 574 & eigvec(2,idir1,iatom2,imod)*eigvec(2,idir2,iatom2,imod) 575 vec2i = eigvec(2,idir1,iatom2,imod)*eigvec(1,idir2,iatom2,imod)-& 576 & eigvec(1,idir1,iatom2,imod)*eigvec(2,idir2,iatom2,imod) 577 578 ! Compute factor for DW term 579 if(phfreq(imod,iqpt)<tol6)then 580 fact2r = zero 581 fact2i = zero 582 else 583 fact2r = -wghtq(iqpt)*(vec1r/amu(typat(iatom1)) + vec2r/amu(typat(iatom2)))/phfreq(imod,iqpt)/& 584 & amu_emass/2 !/norm(idir1)/norm(idir2) 585 fact2i = -wghtq(iqpt)*(vec1i/amu(typat(iatom1)) + vec2i/amu(typat(iatom2)))/phfreq(imod,iqpt)/& 586 & amu_emass/2 !/norm(idir1)/norm(idir2) 587 end if 588 589 multr(:,:) =(ddb_eig2%eig2dval(1,index,:,:)*vecr - ddb_eig2%eig2dval(2,index,:,:)*veci) !/(norm(idir1)*norm(idir2)) 590 multi(:,:) =(ddb_eig2%eig2dval(1,index,:,:)*veci + ddb_eig2%eig2dval(2,index,:,:)*vecr) !/(norm(idir1)*norm(idir2)) 591 592 593 ! Debye-Waller Term 594 if(thmflag==3 .or. thmflag==5 .or. thmflag==7) then 595 dwtermr(1:mband,1:nkpt)=dwtermr(1:mband,1:nkpt)+fact2r*eig2dGamma(1,index,:,:)-fact2i*eig2dGamma(2,index,:,:) 596 dwtermi(1:mband,1:nkpt)=dwtermi(1:mband,1:nkpt)+fact2r*eig2dGamma(2,index,:,:)+fact2i*eig2dGamma(1,index,:,:) 597 end if 598 599 ! Self-energy Term (Fan) 600 deigr(1:mband,1:nkpt) = deigr(1:mband,1:nkpt) + wghtq(iqpt)*factr*multr(1:mband,1:nkpt) 601 deigi(1:mband,1:nkpt) = deigi(1:mband,1:nkpt) + wghtq(iqpt)*factr*multi(1:mband,1:nkpt) 602 603 end do !idir2 604 end do !iatom2 605 end do !idir1 606 end do !iatom1 607 ! Eigenvalue derivative or broadening 608 if(thmflag==3 .or. thmflag==5 .or. thmflag==7) then 609 dednr(1:mband,1:nkpt,imod,iqpt) = deigr(1:mband,1:nkpt) + dwtermr(1:mband,1:nkpt) 610 dedni(1:mband,1:nkpt,imod,iqpt) = deigi(1:mband,1:nkpt) + dwtermi(1:mband,1:nkpt) 611 else if(thmflag==4 .or. thmflag==6 .or. thmflag==8) then 612 dednr(1:mband,1:nkpt,imod,iqpt) = pi*deigr(1:mband,1:nkpt) 613 dedni(1:mband,1:nkpt,imod,iqpt) = pi*deigi(1:mband,1:nkpt) 614 end if 615 616 end do ! imod 617 end do !iqpt 618 619 call ddb_hdr%close() 620 621 622 !============================================================================= 623 !3) Evaluation of the Eliashberg type spectral function 624 !and phonon DOS via gaussian broadning 625 !============================================================================= 626 627 if(telphint==1)then 628 ng2f = 500 ! number of frequencies 629 omega_min=zero 630 omega_max=zero 631 do iqpt=1,nqpt 632 do imod=1,3*natom 633 omega_min = min(omega_min,phfreq(imod,iqpt)) 634 omega_max = max(omega_max,phfreq(imod,iqpt)) 635 end do 636 end do 637 638 ABI_MALLOC(dos_phon,(ng2f)) 639 ABI_MALLOC(g2f,(mband,nkpt,ng2f)) 640 ABI_MALLOC(tmpg2f,(mband,nkpt,ng2f)) 641 ABI_MALLOC(tmpphondos,(ng2f)) 642 643 write(std_out,'(a,es13.6)') 'omega_min :', omega_min 644 write(std_out,'(a,es13.6)') 'omega_max :', omega_max 645 write(std_out,'(a,i8)') 'ng2f :', ng2f 646 647 omega_max = omega_max + 0.1 * omega_max 648 domega = (omega_max-omega_min)/(ng2f-one) 649 650 gaussprefactor = sqrt(piinv) / g2fsmear 651 gaussfactor = one / g2fsmear 652 653 g2f(:,:,:) = zero 654 dos_phon(:) = zero 655 656 do iqpt=1,nqpt 657 do imod=1,3*natom 658 omega = omega_min 659 tmpg2f(:,:,:) = zero 660 tmpphondos(:) = zero 661 do iomega=1,ng2f 662 xx = (omega-phfreq(imod,iqpt))*gaussfactor 663 gaussval = gaussprefactor*exp(-xx*xx) 664 tmpg2f(:,:,iomega) = tmpg2f(:,:,iomega) + gaussval*dednr(:,:,imod,iqpt) 665 tmpphondos(iomega) = tmpphondos(iomega) + gaussval 666 omega = omega+domega 667 end do 668 669 g2f(:,:,:) = g2f(:,:,:) + tmpg2f(:,:,:) 670 dos_phon(:) = dos_phon(:) + tmpphondos(:) 671 672 end do !imod 673 end do !iqpt 674 675 dos_phon(:) = dos_phon(:) / nqpt 676 677 ! output the g2f 678 kpnt(:,:,1) = ddb_eig2%kpt(:,:) 679 unit_g2f = 108 680 call outg2f(domega,omega_min,omega_max,elph_base_name,g2f,g2fsmear,kpnt,mband,ng2f,nkpt,nqpt,1,telphint,unit_g2f) 681 682 ! output the phonon DOS 683 unit_phdos = 108 684 call outphdos(domega,dos_phon,omega_min,omega_max,elph_base_name,g2fsmear,ng2f,nqpt,1,telphint,unit_g2f) 685 686 687 ABI_FREE(dos_phon) 688 ABI_FREE(g2f) 689 ABI_FREE(tmpg2f) 690 ABI_FREE(tmpphondos) 691 692 end if !telphint 693 694 !======================================================================= 695 !4) Evaluation of the Eliashberg type spectral function 696 !and phonon DOS via improved tetrahedron method 697 !======================================================================= 698 699 if(telphint==0)then 700 701 ! make dimension-ful rprimd and gprimd for transformation of derivatives to cartesian coordinates. 702 call mkrdim(acell,rprim,rprimd) 703 call matr3inv(rprimd,gprimd) 704 705 ! Q point Grid 706 qpt_full(:,:) = ddb%qpt(1:3,:) 707 708 ! Trivial Q point index 709 do iqpt=1,nqpt 710 indqpt(iqpt)=iqpt 711 qptnrm(iqpt)= qpt_full(1,iqpt)*qpt_full(1,iqpt)+qpt_full(2,iqpt)*qpt_full(2,iqpt)+qpt_full(3,iqpt)*qpt_full(3,iqpt) 712 end do 713 714 ! Build qlatt from scratch (for 5.7) 715 tol = 0.1_dp 716 ilatt = 0 717 call sort_dp(nqpt,qptnrm,indqpt,tol) 718 719 do iqpt1=1,nqpt-2 720 mesh(1:3,1) = qpt_full(1:3,indqpt(iqpt1)) 721 do iqpt2=iqpt1+1,nqpt-1 722 mesh(1:3,2)= qpt_full(1:3,indqpt(iqpt2)) 723 do iqpt3=iqpt2+1,nqpt 724 mesh(1:3,3)= qpt_full(1:3,indqpt(iqpt3)) 725 det = mesh(1,1)*mesh(2,2)*mesh(3,3) + mesh(1,2)*mesh(2,3)*mesh(3,1) + mesh(1,3)*mesh(2,1)*mesh(3,2) & 726 & -mesh(3,1)*mesh(2,2)*mesh(1,3) - mesh(3,2)*mesh(2,3)*mesh(1,1) - mesh(3,3)*mesh(2,1)*mesh(1,2) 727 invdet = one/det 728 if (abs(nint(invdet))==nqpt .and. abs(invdet)-nqpt < tol) then 729 ilatt = 1 730 qlatt(:,:) = mesh(:,:) 731 exit 732 end if 733 end do 734 if(ilatt==1) exit 735 end do 736 if(ilatt==1) exit 737 end do 738 739 ! error message if qlatt not found and stop 740 if(ilatt==0) then 741 write(message, '(a,a)' ) & 742 & ' Could not find homogeneous basis vectors for Q point grid ',ch10 743 call wrtout(std_out,message,'COLL') 744 call wrtout(ab_out,message,'COLL') 745 ABI_ERROR("Aborting now") 746 end if 747 748 ! test if qlatt is righthanded and possibly fixe it 749 if(invdet < 0) then 750 tempqlatt(:) = qlatt(:,2) 751 qlatt(:,2) = qlatt(:,1) 752 qlatt(:,1) = tempqlatt(:) 753 end if 754 755 write(std_out,*) 'qlatt',qlatt 756 757 ! test if qlatt generates all Q points TO DO 758 759 ! Get tetrahedra, ie indexes of the full kpoints at their summits 760 call init_tetra(indqpt,gprimd,qlatt,qpt_full,nqpt, tetrahedra, ierr, errstr, xmpi_comm_self) 761 !call htetra_init(tetra, indqpt, gprimd, qlatt, qpt_full, nqpt, kpt_ibz, nkpt_ibz, ierr, errstr, xmpi_comm_self 762 ABI_CHECK(ierr==0,errstr) 763 764 rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) & 765 & -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) & 766 & +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3))) 767 768 ! Calculate weights for phonon DOS 769 ! Special precautions must be taking for Gamma point 770 ! because of non-analytic term. 771 ! Non-analyticity must be taken out and treated separatly. 772 773 nene = 100 !nene=number of energies for DOS 774 enemin = minval(phfreq) 775 enemax = maxval(phfreq) 776 deltaene = (enemax-enemin)/dble(nene-1) 777 ! redefine enemin enemax to be at rounded multiples of deltaene 778 ! enemin = elph_ds%fermie - dble(ifermi)*deltaene 779 ! enemax = elph_ds%fermie + dble(nene-ifermi-1)*deltaene 780 781 ABI_MALLOC(tweight,(nqpt,nene)) 782 ABI_MALLOC(dtweightde,(nqpt,nene)) 783 ABI_MALLOC(intweight,(3*natom,nqpt,nene)) 784 ABI_MALLOC(indtweightde,(3*natom,nqpt,nene)) 785 786 do iband=1,3*natom 787 eigen_in(:) = phfreq(iband,:) 788 789 ! calculate general integration weights at each irred kpoint 790 ! as in Blochl et al PRB 49 16223 [[cite:Bloechl1994a]] 791 call get_tetra_weight(eigen_in,enemin,enemax,& 792 & one,nene,nqpt,tetrahedra,bcorr0,& 793 & tweight,dtweightde,xmpi_comm_self) 794 795 intweight(iband,:,:) = tweight(:,:) 796 indtweightde(iband,:,:) = dtweightde(:,:) 797 798 end do !iband 799 800 ! intdtweightse(nband,nqpt,nene) represents the weight in each energy bin for every kpt and every band 801 ! So phonon DOS is calculated (neglecting the non-analyticity contribution for now !!!) 802 803 ABI_MALLOC(total_dos,(nene)) 804 ABI_MALLOC(g2f,(mband,nkpt,nene)) 805 806 total_dos(:) = zero 807 do iband=1,3*natom 808 do iqpt=1,nqpt 809 total_dos(:) = total_dos + indtweightde(iband,iqpt,:) 810 end do 811 end do 812 813 ! For the g2f function 814 ! Right now for one electronic band and one K point: dednr(1:mband,1:nkpt,imod,iqpt) 815 ! Once again must pay close attention to the Gamma point 816 g2f(:,:,:) = zero 817 do ii=1,mband 818 do ikpt=1,nkpt 819 do iband=1,3*natom 820 do iqpt=1,nqpt 821 g2f(ii,ikpt,:) = g2f(ii,ikpt,:) + dednr(ii,ikpt,iband,iqpt) * indtweightde(iband,iqpt,:) 822 end do 823 end do 824 end do 825 end do 826 827 ! output the g2f 828 unit_g2f = 108 829 call outg2f(deltaene,enemin,enemax,elph_base_name,g2f,g2fsmear,kpnt,mband,nene,nkpt,nqpt,tetrahedra%ntetra,telphint,unit_g2f) 830 831 ! output the phonon DOS 832 unit_phdos = 108 833 call outphdos(deltaene,total_dos,enemin,enemax,elph_base_name,g2fsmear,nene,nqpt,tetrahedra%ntetra,telphint,unit_g2f) 834 835 ABI_FREE(tweight) 836 ABI_FREE(dtweightde) 837 ABI_FREE(intweight) 838 ABI_FREE(indtweightde) 839 ABI_FREE(total_dos) 840 ABI_FREE(g2f) 841 end if !telphint 842 843 !======================================================================= 844 !5) direct evaluation of thermal corrections 845 !======================================================================= 846 847 !open TBS file 848 outfile = trim(elph_base_name)//"_TBS" 849 if (open_file(outfile,message,newunit=unitout,form='formatted',status='unknown') /= 0) then 850 ABI_ERROR(message) 851 end if 852 write(unitout,'(a)')'thmeig: Thermal Eigenvalue corrections (eV)' 853 854 slope(:,:,:) = zero 855 zeropoint(:,:,:) = zero 856 !Loop on temperatures 857 do itemper= 1, ntemper 858 tmp=tempermin+temperinc*float(itemper-1) 859 thmeigen(:,:,:) = zero 860 861 ! Sum on all phonon wavevectors and modes 862 do iqpt=1,nqpt 863 do imod=1,3*natom 864 865 ! Bose-Einstein distribution 866 ! jmb overflow with exp(). So, select bosein to be still significant wrt half 867 if(phfreq(imod,iqpt)<tol6 .or. (phfreq(imod,iqpt)/(kb_HaK*tmp)) > -log(tol16))then 868 bosein = zero 869 else 870 bosein = one/(exp(phfreq(imod,iqpt)/(kb_HaK*tmp))-one) 871 end if 872 873 ! Calculate total 874 thmeigen(1,1:mband,1:nkpt) = thmeigen(1,1:mband,1:nkpt) + dednr(1:mband,1:nkpt,imod,iqpt)*(bosein+half) 875 thmeigen(2,1:mband,1:nkpt) = thmeigen(2,1:mband,1:nkpt) + dedni(1:mband,1:nkpt,imod,iqpt)*(bosein+half) 876 877 if(itemper==1)then 878 ! Calculate slope of linear regime 879 if(phfreq(imod,iqpt)<tol6)then 880 slope(1,1:mband,1:nkpt) = slope(1,1:mband,1:nkpt) 881 slope(2,1:mband,1:nkpt) = slope(2,1:mband,1:nkpt) 882 else 883 slope(1,1:mband,1:nkpt) = slope(1,1:mband,1:nkpt) + dednr(1:mband,1:nkpt,imod,iqpt)*(kb_HaK/phfreq(imod,iqpt)) 884 slope(2,1:mband,1:nkpt) = slope(2,1:mband,1:nkpt) + dedni(1:mband,1:nkpt,imod,iqpt)*(kb_HaK/phfreq(imod,iqpt)) 885 end if 886 ! Calculate zero-point renormalization 887 zeropoint(1,1:mband,1:nkpt) = zeropoint(1,1:mband,1:nkpt) + dednr(1:mband,1:nkpt,imod,iqpt)*half 888 zeropoint(2,1:mband,1:nkpt) = zeropoint(2,1:mband,1:nkpt) + dedni(1:mband,1:nkpt,imod,iqpt)*half 889 890 end if 891 end do ! imod 892 end do !iqpt 893 894 ! Write temperature independent results 895 if(itemper==1)then 896 write(unitout,'(a)')'Temperature independent results (zero-point renormalization and slope)' 897 do ikpt=1,nkpt 898 write(unitout,'(a,3es16.8)')' Kpt :', kpnt(:,ikpt,1) 899 do iband=1,mband 900 write(unitout,'(4d22.14)') Ha_eV*zeropoint(1,iband,ikpt),Ha_eV*zeropoint(2,iband,ikpt),& 901 & Ha_eV*slope(1,iband,ikpt),Ha_eV*slope(2,iband,ikpt) 902 end do 903 end do 904 write(unitout,'(a)')'Temperature dependent corrections' 905 end if 906 ! Write result in file for each temperature 907 write(unitout,'(a,es10.3,a)')'T :', tmp,' K' 908 do ikpt=1,nkpt 909 write(unitout,'(a,3es16.8)')' Kpt :', kpnt(:,ikpt,1) 910 do iband=1,mband 911 write(unitout,'(2d22.14)') Ha_eV*thmeigen(1,iband,ikpt), Ha_eV*thmeigen(2,iband,ikpt) 912 end do 913 end do 914 end do !itemper 915 916 close(unitout) 917 918 !Write temperature-independent results to the main output file 919 write(iout,'(a)')' ' 920 write(iout,'(80a)') ('-',ii=1,80) 921 write(iout,'(a)')' ' 922 write(iout,'(a)')' Electron-phonon change of electronic structure.' 923 write(iout,'(a)')' The temperature-dependent values are written in the _TBS file.' 924 write(iout,'(a)')' Here follows, for each electronic wavevector and band :' 925 write(iout,'(a)')' zero-point renormalisation (Ha) and linear slope (Ha/Kelvin)' 926 do ikpt=1,nkpt 927 write(iout,'(2a,i6,a,3es16.6)')ch10,' Kpt number ',ikpt,', with reduced coordinates :',kpnt(:,ikpt,1) 928 do iband=1,mband 929 write(iout,'(i6,2es20.6)') iband,zeropoint(1,iband,ikpt),slope(1,iband,ikpt) 930 end do 931 end do 932 933 ABI_FREE(typat) 934 ABI_FREE(atifc) 935 ABI_FREE(zion) 936 ABI_FREE(amu) 937 ABI_FREE(xcart) 938 ABI_FREE(xred) 939 ABI_FREE(symafm) 940 ABI_FREE(spinat) 941 ABI_FREE(symrel) 942 ABI_FREE(symrec) 943 ABI_FREE(indsym) 944 ABI_FREE(tnons) 945 ABI_FREE(deigi) 946 ABI_FREE(deigr) 947 ABI_FREE(dwtermi) 948 ABI_FREE(dwtermr) 949 ABI_FREE(multi) 950 ABI_FREE(multr) 951 ABI_FREE(slope) 952 ABI_FREE(thmeigen) 953 ABI_FREE(zeropoint) 954 955 ABI_FREE(dedni) 956 ABI_FREE(dednr) 957 if(allocated(indqpt)) then 958 ABI_FREE(indqpt) 959 end if 960 ABI_FREE(eigen_in) 961 ABI_FREE(qpt_full) 962 ABI_FREE(qptnrm) 963 ABI_FREE(wghtq) 964 ABI_FREE(spqpt) 965 ABI_FREE(eigvec) 966 ABI_FREE(phfreq) 967 968 ABI_FREE(eig2dGamma) 969 ABI_FREE(kpnt) 970 ABI_FREE(carflg_eig2) 971 972 call ddb_eig2%free() 973 call destroy_tetra(tetrahedra) 974 975 end subroutine thmeig