TABLE OF CONTENTS
ABINIT/outg2f [ 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
PARENTS
thmeig
CHILDREN
SOURCE
1127 subroutine outg2f(deltaene,enemin,enemax,filnam,g2f,g2fsmear,kpnt,mband,nene,nkpt,nqpt,ntetra,telphint,unit_g2f) 1128 1129 1130 !This section has been created automatically by the script Abilint (TD). 1131 !Do not modify the following lines by hand. 1132 #undef ABI_FUNC 1133 #define ABI_FUNC 'outg2f' 1134 !End of the abilint section 1135 1136 implicit none 1137 1138 !Arguments ------------------------------------ 1139 !scalars 1140 integer,intent(in) :: mband,nene,nkpt,nqpt,ntetra,telphint,unit_g2f 1141 character(len=fnlen),intent(in) :: filnam 1142 real(dp) :: deltaene,enemin,enemax,g2fsmear 1143 !arrays 1144 real(dp) :: g2f(mband,nkpt,nene),kpnt(3,nkpt,nqpt) 1145 1146 !Local variables------------------------------- 1147 !scalars 1148 integer :: iband,ikpt,iomega,iost 1149 real(dp) :: omega 1150 character(len=fnlen) :: outfile 1151 character(len=500) :: message 1152 !arrays 1153 1154 ! ************************************************************************* 1155 1156 !output the g2f 1157 outfile = trim(filnam) // '_G2F' 1158 open (unit=unit_g2f,file=outfile,status='unknown',iostat=iost) 1159 if (iost /= 0) then 1160 write (message,'(3a)')' thmeig : ERROR- opening file ',trim(outfile),' as new' 1161 MSG_ERROR(message) 1162 end if 1163 1164 write(std_out,*) ' g2f function' 1165 write (unit_g2f,'(a)') '#' 1166 write (unit_g2f,'(a)') '# ABINIT package : g2f file' 1167 write (unit_g2f,'(a)') '#' 1168 write (unit_g2f,'(a,I10)') '# number of qpoints integrated over : ', nqpt 1169 write (unit_g2f,'(a,I10)') '# number of energy points : ', nene 1170 write (unit_g2f,'(a,E16.6,a,E16.6,a)') '# between omega_min = ', enemin, & 1171 & ' Ha and omega_max = ', enemax, ' Ha' 1172 if(telphint==1)then 1173 write (unit_g2f,'(a,E16.6)') '# and the smearing width for gaussians is ', g2fsmear 1174 write (unit_g2f,'(a)') '#' 1175 end if 1176 if(telphint==0)then 1177 write (unit_g2f,'(a,I10)') '# number of tetrahedrons', ntetra 1178 write (unit_g2f,'(a)') '#' 1179 end if 1180 1181 !Write only the a2f function for the first K point 1182 !ikpt=1 1183 do ikpt=1,nkpt 1184 write(unit_g2f,'(a,3es16.8)')' Kpt :', kpnt(:,ikpt,1) 1185 do iband=1,mband 1186 write(unit_g2f,*) 'band :', iband 1187 omega = enemin 1188 do iomega=1,nene 1189 write (unit_g2f,*) omega*Ha_eV*1000, g2f(iband, ikpt,iomega) 1190 omega=omega+deltaene 1191 end do 1192 end do 1193 end do 1194 1195 close (unit=unit_g2f) 1196 1197 end subroutine outg2f
ABINIT/outphdos [ 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)
PARENTS
thmeig
CHILDREN
SOURCE
1025 subroutine outphdos(deltaene,dos_phon,enemin,enemax,filnam,g2fsmear,nene,nqpt,ntetra,telphint,unit_phdos) 1026 1027 1028 !This section has been created automatically by the script Abilint (TD). 1029 !Do not modify the following lines by hand. 1030 #undef ABI_FUNC 1031 #define ABI_FUNC 'outphdos' 1032 use interfaces_14_hidewrite 1033 !End of the abilint section 1034 1035 implicit none 1036 1037 !Arguments ------------------------------------ 1038 !scalars 1039 integer,intent(in) :: nene,nqpt,ntetra,telphint,unit_phdos 1040 character(len=fnlen),intent(in) :: filnam 1041 real(dp) :: deltaene,enemin,enemax,g2fsmear 1042 !arrays 1043 real(dp) :: dos_phon(nene) 1044 1045 !Local variables------------------------------- 1046 !scalars 1047 integer :: iomega,iost,step10 1048 real(dp) :: dos_effective,omega 1049 character(len=fnlen) :: outfile 1050 character(len=500) :: message 1051 !arrays 1052 1053 ! ************************************************************************* 1054 1055 outfile = trim(filnam) // '_PDS' 1056 write(message, '(3a)')ch10,& 1057 & ' Will write phonon DOS in file ',trim(outfile) 1058 call wrtout(ab_out,message,'COLL') 1059 call wrtout(std_out,message,'COLL') 1060 1061 write(message, '(4a)')ch10,& 1062 & ' For checking purposes, write ten values in the present file.',ch10,& 1063 & ' Index Energy (in Ha) DOS ' 1064 call wrtout(ab_out,message,'COLL') 1065 call wrtout(std_out,message,'COLL') 1066 1067 open (unit=unit_phdos,file=outfile,status='replace',iostat=iost) 1068 if (iost /= 0) then 1069 write (message,'(3a)')' Opening file ',trim(outfile),' as new' 1070 MSG_ERROR(message) 1071 end if 1072 1073 write (unit_phdos,'(a)') '#' 1074 write (unit_phdos,'(a)') '# ABINIT package : phonon DOS file' 1075 write (unit_phdos,'(a)') '#' 1076 write (unit_phdos,'(a,i10)') '# Number of Qpoints integrated over : ', nqpt 1077 write (unit_phdos,'(a,i10)') '# Number of energy points : ', nene 1078 write (unit_phdos,'(a,es16.6,a,es16.6,a)') '# between omega_min = ', enemin, & 1079 & ' Ha and omega_max = ', enemax, ' Ha' 1080 if(telphint==1)then 1081 write (unit_phdos,'(a,es16.6)') '# The smearing width for gaussians is ', g2fsmear 1082 end if 1083 if(telphint==0)then 1084 write (unit_phdos,'(a,i10)') '# Number of tetrahedrons', ntetra 1085 end if 1086 write (unit_phdos,'(a)') '#' 1087 write (unit_phdos,'(a)') '# Index Energy (in Ha) DOS ' 1088 1089 omega = enemin 1090 do iomega=1,nene 1091 dos_effective=dos_phon(iomega) 1092 if(abs(dos_effective)<tol16)then 1093 dos_effective=zero 1094 end if 1095 step10=nene/10 1096 if(mod(iomega,step10)==1)write (std_out,'(i10,es18.6,es18.6)')iomega, omega, dos_effective 1097 if(mod(iomega,step10)==1)write (ab_out,'(i10,es18.6,es18.6)')iomega, omega, dos_effective 1098 write (unit_phdos, '(i10,es18.6,es18.6)')iomega, omega, dos_effective 1099 omega=omega+deltaene 1100 end do 1101 1102 close (unit=unit_phdos) 1103 1104 end subroutine outphdos
ABINIT/thmeig [ Functions ]
NAME
thmeig
FUNCTION
This routine calculates the thermal corrections to the eigenvalues. The output is this quantity for the input k point.
COPYRIGHT
Copyright (C) 1999-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors .
INPUTS
elph_base_name = root filename for outputs eig2_filnam = name of the eig2 database file comm=MPI communicator
OUTPUT
PARENTS
anaddb
CHILDREN
SOURCE
31 #if defined HAVE_CONFIG_H 32 #include "config.h" 33 #endif 34 35 #include "abi_common.h" 36 37 subroutine thmeig(inp, ddb, crystal, & 38 & elph_base_name, eig2_filnam, ddbun, iout, & 39 & natom, mpert, msize, d2asr, & 40 & comm) 41 42 use defs_basis 43 use m_profiling_abi 44 use m_tetrahedron 45 use m_errors 46 use m_ddb 47 use m_ddb_hdr 48 use m_xmpi 49 use m_sort 50 51 use m_geometry, only : mkrdim, xred2xcart, metric 52 use m_crystal, only : crystal_t 53 use m_io_tools, only : open_file 54 use m_dynmat, only : asria_corr, dfpt_phfrq 55 use m_anaddb_dataset, only : anaddb_dataset_type 56 use m_pawtab, only : pawtab_type,pawtab_nullify,pawtab_free 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 'thmeig' 62 use interfaces_14_hidewrite 63 use interfaces_32_util 64 use interfaces_41_geometry 65 use interfaces_56_recipspace 66 !End of the abilint section 67 68 implicit none 69 70 !Arguments ------------------------------------ 71 !scalars 72 integer,intent(inout) :: natom 73 integer,intent(in) :: mpert,msize 74 integer,intent(in) :: comm 75 character(len=*),intent(in) :: elph_base_name, eig2_filnam 76 integer,intent(in) :: ddbun,iout 77 type(crystal_t), intent(inout) :: crystal 78 type(anaddb_dataset_type),intent(inout) :: inp 79 type(ddb_type),intent(inout) :: ddb 80 !arrays 81 real(dp),intent(inout) :: d2asr(2,3,natom,3,natom) 82 83 84 !Local variables------------------------------- 85 !scalars 86 integer,parameter :: msppol=2,master=0,bcorr0=0 87 integer :: msym 88 integer :: nkpt,mband,ntypat 89 integer :: usepaw,natifc 90 integer :: nsym,occopt,nblok2 91 integer :: ntemper,telphint,thmflag 92 integer :: brav,chksymbreak,found,gqpt,iatom1,iatom2,iband,iblok,iblok2,idir1,idir2,ii,ikpt,ilatt,imod,index 93 integer :: iomega,iqpt,iqpt1,iqpt2,iqpt2_previous,iqpt3,iscf_fake,itemper 94 integer :: mpert_eig2,msize2,nene,ng2f,nqshft,nsym_new,unit_g2f,nqpt,nqpt_computed,qptopt,rftyp 95 !integer :: mqpt,nqpt2,option 96 integer :: unit_phdos,unitout 97 integer :: isym 98 integer :: nptsym,use_inversion 99 integer :: ierr 100 real(dp) :: ucvol 101 real(dp) :: g2fsmear,temperinc,tempermin 102 real(dp) :: bosein,deltaene,det,domega,enemax,enemin,fact2i,fact2r,factr 103 real(dp) :: gaussfactor,gaussprefactor,gaussval,invdet,omega,omega_max,omega_min,qnrm,qptrlen 104 real(dp) :: rcvol,tmp,tol,vec1i,vec1r,vec2i,vec2r,veci,vecr,xx 105 real(dp) :: tolsym,tolsym8 !new 106 character(len=500) :: message 107 character(len=fnlen) :: outfile 108 type(ddb_type) :: ddb_eig2 109 type(ddb_hdr_type) :: ddb_hdr 110 !arrays 111 ! FIXME now these must be allocated 112 integer :: ngqpt(9),qptrlatt(3,3),rfelfd(4),rfphon(4),rfstrs(4),vacuum(3) 113 integer :: bravais(11) 114 integer,allocatable :: typat(:),atifc(:) 115 integer,allocatable :: symrel(:,:,:),symrec(:,:,:) 116 integer,allocatable :: indsym(:,:,:) 117 integer,allocatable :: indqpt(:) 118 integer,allocatable :: symafm(:),symafm_new(:) 119 integer,allocatable :: carflg_eig2(:,:,:,:) 120 integer,allocatable :: ptsymrel(:,:,:),symrel_new(:,:,:) 121 !integer,allocatable :: symrec_new(:,:,:) 122 real(dp) :: rprim(3,3),gprim(3,3),rmet(3,3),gmet(3,3) 123 real(dp) :: acell(3) 124 real(dp) :: diff_qpt(3) 125 real(dp) :: gprimd(3,3),mesh(3,3) 126 real(dp) :: qlatt(3,3),qphnrm(3),qpt_search(3,3) 127 real(dp) :: rprimd(3,3),shiftq(3,210),tempqlatt(3) 128 real(dp) :: dummy(0),dummy2(0,0) 129 real(dp),allocatable :: xcart(:,:),xred(:,:) 130 real(dp),allocatable :: amu(:),zion(:) 131 real(dp),allocatable :: tnons(:,:) 132 real(dp),allocatable :: deigi(:,:), deigr(:,:), multi(:,:), multr(:,:) 133 real(dp),allocatable :: dwtermi(:,:), dwtermr(:,:) 134 real(dp),allocatable :: slope(:,:,:),thmeigen(:,:,:),zeropoint(:,:,:) 135 real(dp),allocatable :: displ(:) 136 real(dp),allocatable :: dos_phon(:),dtweightde(:,:),d2cart(:,:) 137 real(dp),allocatable :: eigvec(:,:,:,:),eigval(:,:),g2f(:,:,:),intweight(:,:,:) 138 real(dp),allocatable :: indtweightde(:,:,:),tmpg2f(:,:,:),tmpphondos(:),total_dos(:),tweight(:,:) 139 real(dp),allocatable :: phfreq(:,:) 140 real(dp),allocatable :: blkval2(:,:,:,:),blkval2gqpt(:,:,:,:),kpnt(:,:,:) 141 real(dp),allocatable :: dedni(:,:,:,:),dednr(:,:,:,:) 142 real(dp),allocatable :: eigen_in(:) 143 real(dp),allocatable :: qpt_full(:,:),qptnrm(:) 144 real(dp),allocatable :: spqpt(:,:),tnons_new(:,:),spinat(:,:) 145 real(dp),allocatable :: wghtq(:) 146 147 type(t_tetrahedron) :: tetrahedra 148 character(len=80) :: errstr 149 150 ! ********************************************************************* 151 152 !DEBUG 153 ! write(std_out,*)'-thmeig : enter ' 154 !call flush(6) 155 !ENDDEBUG 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(ddb_hdr, eig2_filnam, ddbun, DDB_VERSION) 187 188 mband = ddb_hdr%mband 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_ALLOCATE(typat, (natom)) 199 ABI_ALLOCATE(atifc, (natom)) 200 ABI_ALLOCATE(zion, (ntypat)) 201 ABI_ALLOCATE(amu, (ntypat)) 202 203 ABI_ALLOCATE(xcart,(3,natom)) 204 ABI_ALLOCATE(xred,(3,natom)) 205 206 ABI_ALLOCATE(symafm, (nsym)) 207 ABI_ALLOCATE(spinat,(3,natom)) 208 209 ABI_ALLOCATE(symrel, (3,3,nsym)) 210 ABI_ALLOCATE(symrec, (3,3,nsym)) 211 ABI_ALLOCATE(tnons, (3,nsym)) 212 ABI_ALLOCATE(indsym, (4,nsym,natom)) 213 214 ABI_ALLOCATE(deigi, (mband,nkpt)) 215 ABI_ALLOCATE(deigr, (mband,nkpt)) 216 ABI_ALLOCATE(dwtermi, (mband,nkpt)) 217 ABI_ALLOCATE(dwtermr, (mband,nkpt)) 218 ABI_ALLOCATE(multi, (mband,nkpt)) 219 ABI_ALLOCATE(multr, (mband,nkpt)) 220 ABI_ALLOCATE(slope, (2,mband,nkpt)) 221 ABI_ALLOCATE(thmeigen, (2,mband,nkpt)) 222 ABI_ALLOCATE(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 call ddb_malloc(ddb_eig2,msize2,nblok2,natom,ntypat) 229 230 ABI_ALLOCATE(blkval2,(2,msize2,mband,nkpt)) 231 ABI_ALLOCATE(blkval2gqpt,(2,msize2,mband,nkpt)) 232 233 ABI_ALLOCATE(eigvec,(2,3,natom,3*natom)) 234 ABI_ALLOCATE(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 blkval2gqpt(:,:,:,:)=zero 286 287 ABI_ALLOCATE(carflg_eig2,(3,mpert_eig2,3,mpert_eig2)) 288 ABI_ALLOCATE(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 read_blok8(ddb_eig2,iblok2,mband,mpert_eig2,msize2,& 313 & nkpt,ddbun,blkval2(:,:,:,:),kpnt(:,:,1)) 314 315 316 qnrm = ddb_eig2%qpt(1,iblok2)*ddb_eig2%qpt(1,iblok2)+ & 317 & ddb_eig2%qpt(2,iblok2)*ddb_eig2%qpt(2,iblok2)+ & 318 & ddb_eig2%qpt(3,iblok2)*ddb_eig2%qpt(3,iblok2) 319 if(qnrm < DDB_QTOL) then 320 blkval2gqpt(:,:,:,:) = blkval2(:,:,:,:) 321 gqpt=iblok2 322 write(std_out,*)'-thmeig: found Gamma point in EIG2 DDB, blok number ',iblok2 323 found=1 324 exit 325 end if 326 end do 327 328 if(found==0)then 329 write(message,'(a,i3,2a)')& 330 & ' Was unable to find the blok for Gamma point in EIG2 DDB file, while thmflag= ',thmflag,ch10,& 331 & ' Action : compute the contribution from Gamma, and merge it in your EIG2 DDB file.' 332 MSG_ERROR(message) 333 end if 334 335 ! Put blkval2gqpt in cartesian coordinates 336 call carttransf(ddb_eig2%flg,blkval2gqpt,carflg_eig2,gprimd,gqpt,mband,& 337 & mpert_eig2,msize2,natom,nblok2,nkpt,rprimd) 338 339 end if 340 341 close(ddbun) 342 343 !========================================================================= 344 !2) Calculation of dE(n,k)/dn(Q,j) : consider all q and modes 345 !========================================================================= 346 347 if(thmflag==3 .or. thmflag==4)then 348 349 350 ! Use the first list of q wavevectors 351 nqpt=inp%nph1l 352 ABI_ALLOCATE(spqpt,(3,nqpt)) 353 do iqpt=1,inp%nph1l 354 spqpt(:,iqpt)=inp%qph1l(:,iqpt)/inp%qnrml1(iqpt) 355 end do 356 ABI_ALLOCATE(wghtq,(nqpt)) 357 wghtq(:)=one/nqpt 358 359 else if(thmflag>=5 .and. thmflag<=8)then 360 361 ! Generates the q point grid 362 ngqpt(1:3)=inp%ngqpt(1:3) 363 nqshft=inp%nqshft 364 qptrlatt(:,:)=0 365 qptrlatt(1,1)=ngqpt(1) 366 qptrlatt(2,2)=ngqpt(2) 367 qptrlatt(3,3)=ngqpt(3) 368 369 ABI_ALLOCATE(ptsymrel,(3,3,msym)) 370 ABI_ALLOCATE(symafm_new,(msym)) 371 ABI_ALLOCATE(symrel_new,(3,3,msym)) 372 ABI_ALLOCATE(tnons_new,(3,msym)) 373 if(thmflag==7 .or. thmflag==8) then 374 ! Re-generate symmetry operations from the lattice and atomic coordinates 375 tolsym=tol8 376 call symlatt(bravais,msym,nptsym,ptsymrel,rprimd,tolsym) 377 use_inversion=1 378 call symfind(0,(/zero,zero,zero/),gprimd,0,msym,natom,0,nptsym,nsym_new,0,& 379 & ptsymrel,spinat,symafm_new,symrel_new,tnons_new,tolsym,typat,use_inversion,xred) 380 write(std_out,*)' thmeig : found ',nsym_new,' symmetries ',ch10 381 qptopt=1 382 else 383 nsym_new=1 384 symrel_new(:,:,1)=0 ; symrel_new(1,1,1)=1 ; symrel_new(2,2,1)=1 ; symrel_new(3,3,1)=1 385 tnons_new(:,1)=zero 386 symafm_new(1)=1 387 qptopt=3 388 end if 389 390 brav=inp%brav 391 392 if(abs(brav)/=1)then 393 message = ' The possibility to have abs(brav)/=1 for thmeig was disabled.' 394 MSG_ERROR(message) 395 end if 396 397 ! Prepare to compute the q-point grid in the ZB or IZB 398 iscf_fake=5 ! Need the weights 399 chksymbreak=0 400 vacuum=0 401 shiftq(:,1:nqshft)=inp%q1shft(:,1:nqshft) 402 ! Compute the final number of q points 403 call getkgrid(chksymbreak,0,iscf_fake,dummy2,qptopt,qptrlatt,qptrlen,& 404 & nsym_new,0,nqpt,nqshft,nsym_new,rprimd,& 405 & shiftq,symafm_new,symrel_new,vacuum,dummy) 406 ABI_ALLOCATE(spqpt,(3,nqpt)) 407 ABI_ALLOCATE(wghtq,(nqpt)) 408 call getkgrid(chksymbreak,iout,iscf_fake,spqpt,qptopt,qptrlatt,qptrlen,& 409 & nsym_new,nqpt,nqpt_computed,nqshft,nsym_new,rprimd,& 410 & shiftq,symafm_new,symrel_new,vacuum,wghtq) 411 412 ABI_DEALLOCATE(ptsymrel) 413 ABI_DEALLOCATE(symafm_new) 414 ABI_DEALLOCATE(symrel_new) 415 ABI_DEALLOCATE(tnons_new) 416 417 end if 418 419 call ddb_hdr_free(ddb_hdr) 420 421 422 write(message,'(a,a)')ch10,' thmeig : list of q wavevectors, with integration weights ' 423 call wrtout(ab_out,message,'COLL') 424 call wrtout(std_out,message,'COLL') 425 do iqpt=1,nqpt 426 write(message,'(i6,3es16.6,es20.6)')iqpt,spqpt(:,iqpt),wghtq(iqpt) 427 call wrtout(ab_out,message,'COLL') 428 call wrtout(std_out,message,'COLL') 429 end do 430 431 if(.not.allocated(indqpt))then 432 ABI_ALLOCATE(indqpt,(nqpt)) 433 end if 434 ABI_ALLOCATE(dedni,(mband,nkpt,3*natom,nqpt)) 435 ABI_ALLOCATE(dednr,(mband,nkpt,3*natom,nqpt)) 436 ABI_ALLOCATE(eigen_in,(nqpt)) 437 ABI_ALLOCATE(qpt_full,(3,nqpt)) 438 ABI_ALLOCATE(qptnrm,(nqpt)) 439 440 dednr(:,:,:,:) = zero 441 dedni(:,:,:,:) = zero 442 443 !!Prepare the reading of the EIG2 files 444 call ddb_hdr_open_read(ddb_hdr, eig2_filnam, ddbun, DDB_VERSION, msym=msym) 445 446 call ddb_hdr_free(ddb_hdr) 447 448 !iqpt2 will be the index of the q point bloks inside the EIG2 file 449 iqpt2=0 450 451 !Sum on all phonon wavevectors and modes 452 do iqpt=1,nqpt 453 454 ! Finding the target wavevector in DDB file 455 qpt_search(:,:)=0.0d0 456 qpt_search(:,1)=spqpt(:,iqpt) 457 qphnrm(:)=one 458 rfphon(1:2)=1 459 ! NOTE : at present, no LO-TO splitting included !!! 460 rfelfd(1:2)=0 461 rfstrs(1:2)=0 462 rftyp=1 463 464 write(std_out,'(a,3es16.6)' )' Looking for spqpt=',qpt_search(:,1) 465 466 call gtblk9(ddb,iblok,qpt_search,qphnrm,rfphon,rfelfd,rfstrs,rftyp) 467 468 if(iblok==0) then 469 write(message,'(a,3es16.6,2a)')& 470 & ' Was unable to find in DDB file, the blok for point ',spqpt(:,iqpt),ch10,& 471 & ' Action : compute the contribution from this point, and merge it in your DDB file.' 472 MSG_ERROR(message) 473 end if 474 475 ABI_ALLOCATE(d2cart,(2,msize)) 476 ! Copy the dynamical matrix in d2cart 477 d2cart(:,1:msize)=ddb%val(:,:,iblok) 478 479 ! Eventually impose the acoustic sum rule based on previously calculated d2asr 480 !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart) 481 if (inp%asr==1 .or. inp%asr==2 .or. inp%asr==5) then 482 call asria_corr(inp%asr,d2asr,d2cart,mpert,natom) 483 end if 484 485 ! Calculation of the eigenvectors and eigenvalues 486 ! of the dynamical matrix 487 ABI_ALLOCATE(displ,(2*3*natom*3*natom)) 488 ABI_ALLOCATE(eigval,(3,natom)) 489 call dfpt_phfrq(amu,displ,d2cart,eigval,eigvec,indsym,& 490 & mpert,msym,natom,nsym,ntypat,phfreq(:,iqpt),qphnrm(1),spqpt(:,iqpt),rprimd,inp%symdynmat,& 491 & symrel,symafm,typat,ucvol) 492 ABI_DEALLOCATE(displ) 493 ABI_DEALLOCATE(eigval) 494 ABI_DEALLOCATE(d2cart) 495 496 497 ! Read the next bloks to find the next q point. 498 found=0 ; iqpt2_previous=iqpt2 499 do while (iqpt2<nblok2) 500 iqpt2=iqpt2+1 501 call read_blok8(ddb_eig2,iqpt2,mband,mpert_eig2,msize2,& 502 & nkpt,ddbun,blkval2(:,:,:,:),kpnt(:,:,1)) 503 !write (300,*) 'blkval2 _bis_ in thmeig' 504 !write (300,*) blkval2 505 diff_qpt(:)=ddb_eig2%qpt(1:3,iqpt2)/ddb_eig2%nrm(1,iqpt2)-spqpt(:,iqpt) 506 if(diff_qpt(1)**2+diff_qpt(2)**2+diff_qpt(3)**2 < DDB_QTOL )then 507 found=1 508 exit 509 end if 510 end do 511 512 ! Usually, the q points come in the right order. However, this is not always the case... 513 if(found==0)then 514 515 ! If the EIG2 database file has to be read again, close it, then search for the right q point, 516 ! from the beginning of the file 517 close(ddbun) 518 519 call ddb_hdr_open_read(ddb_hdr, eig2_filnam, ddbun, DDB_VERSION, msym=msym) 520 call ddb_hdr_free(ddb_hdr) 521 522 ! And examine again the EIG2 file. Still, not beyond the previously examined value. 523 found=0 524 do iqpt2=1,iqpt2_previous 525 call read_blok8(ddb_eig2,iqpt2,mband,mpert_eig2,msize2,& 526 & nkpt,ddbun,blkval2(:,:,:,:),kpnt(:,:,1)) 527 diff_qpt(:)=ddb_eig2%qpt(1:3,iqpt2)/ddb_eig2%nrm(1,iqpt2)-spqpt(:,iqpt) 528 if(diff_qpt(1)**2+diff_qpt(2)**2+diff_qpt(3)**2 < DDB_QTOL )then 529 found=1 530 exit 531 end if 532 end do 533 534 if(found==0)then 535 write(message,'(a,3es16.6,2a)')& 536 & ' Was unable to find in EIG2 DDB file, the blok for point ',spqpt(:,iqpt),ch10,& 537 & ' Action : compute the contribution from this point, and merge it in your EIG2 DDB file.' 538 MSG_ERROR(message) 539 end if 540 541 end if 542 543 ! Put blkval2 in cartesian coordinates 544 call carttransf(ddb_eig2%flg,blkval2,carflg_eig2,gprimd,iqpt,mband,& 545 & mpert_eig2,msize2,natom,nblok2,nkpt,rprimd) 546 547 do imod=1,3*natom 548 549 ! Calculate the derivative 550 deigr(:,:) = zero 551 deigi(:,:) = zero 552 dwtermr(:,:)=zero 553 dwtermi(:,:)=zero 554 index=0 555 do iatom1=1,natom 556 do idir1=1,3 557 do iatom2=1,natom 558 ! Compute factor for SE term 559 if(phfreq(imod,iqpt)<tol6)then 560 factr = zero 561 else 562 factr=one/sqrt(amu(typat(iatom1))*amu(typat(iatom2)))/phfreq(imod,iqpt)/amu_emass 563 end if 564 565 do idir2=1,3 566 index = idir1 + 3*((iatom1 - 1) + natom * ((idir2-1)+3*(iatom2-1))) 567 568 ! Compute products of polarization vectors 569 vecr = eigvec(1,idir1,iatom1,imod)*eigvec(1,idir2,iatom2,imod)+& 570 & eigvec(2,idir1,iatom1,imod)*eigvec(2,idir2,iatom2,imod) 571 veci = eigvec(2,idir1,iatom1,imod)*eigvec(1,idir2,iatom2,imod)-& 572 & eigvec(1,idir1,iatom1,imod)*eigvec(2,idir2,iatom2,imod) 573 574 vec1r = eigvec(1,idir1,iatom1,imod)*eigvec(1,idir2,iatom1,imod)+& 575 & eigvec(2,idir1,iatom1,imod)*eigvec(2,idir2,iatom1,imod) 576 vec1i = eigvec(2,idir1,iatom1,imod)*eigvec(1,idir2,iatom1,imod)-& 577 & eigvec(1,idir1,iatom1,imod)*eigvec(2,idir2,iatom1,imod) 578 579 vec2r = eigvec(1,idir1,iatom2,imod)*eigvec(1,idir2,iatom2,imod)+& 580 & eigvec(2,idir1,iatom2,imod)*eigvec(2,idir2,iatom2,imod) 581 vec2i = eigvec(2,idir1,iatom2,imod)*eigvec(1,idir2,iatom2,imod)-& 582 & eigvec(1,idir1,iatom2,imod)*eigvec(2,idir2,iatom2,imod) 583 584 ! Compute factor for DW term 585 if(phfreq(imod,iqpt)<tol6)then 586 fact2r = zero 587 fact2i = zero 588 else 589 fact2r = -wghtq(iqpt)*(vec1r/amu(typat(iatom1)) + vec2r/amu(typat(iatom2)))/phfreq(imod,iqpt)/& 590 & amu_emass/2 !/norm(idir1)/norm(idir2) 591 fact2i = -wghtq(iqpt)*(vec1i/amu(typat(iatom1)) + vec2i/amu(typat(iatom2)))/phfreq(imod,iqpt)/& 592 & amu_emass/2 !/norm(idir1)/norm(idir2) 593 end if 594 595 multr(:,:) =(blkval2(1,index,:,:)*vecr - blkval2(2,index,:,:)*veci) !/(norm(idir1)*norm(idir2)) 596 multi(:,:) =(blkval2(1,index,:,:)*veci + blkval2(2,index,:,:)*vecr) !/(norm(idir1)*norm(idir2)) 597 598 599 ! Debye-Waller Term 600 if(thmflag==3 .or. thmflag==5 .or. thmflag==7) then 601 dwtermr(1:mband,1:nkpt)=dwtermr(1:mband,1:nkpt)+fact2r*blkval2gqpt(1,index,:,:)-fact2i*blkval2gqpt(2,index,:,:) 602 dwtermi(1:mband,1:nkpt)=dwtermi(1:mband,1:nkpt)+fact2r*blkval2gqpt(2,index,:,:)+fact2i*blkval2gqpt(1,index,:,:) 603 end if 604 605 ! Self-energy Term (Fan) 606 deigr(1:mband,1:nkpt) = deigr(1:mband,1:nkpt) + wghtq(iqpt)*factr*multr(1:mband,1:nkpt) 607 deigi(1:mband,1:nkpt) = deigi(1:mband,1:nkpt) + wghtq(iqpt)*factr*multi(1:mband,1:nkpt) 608 609 end do !idir2 610 end do !iatom2 611 end do !idir1 612 end do !iatom1 613 ! Eigenvalue derivative or broadening 614 if(thmflag==3 .or. thmflag==5 .or. thmflag==7) then 615 dednr(1:mband,1:nkpt,imod,iqpt) = deigr(1:mband,1:nkpt) + dwtermr(1:mband,1:nkpt) 616 dedni(1:mband,1:nkpt,imod,iqpt) = deigi(1:mband,1:nkpt) + dwtermi(1:mband,1:nkpt) 617 else if(thmflag==4 .or. thmflag==6 .or. thmflag==8) then 618 dednr(1:mband,1:nkpt,imod,iqpt) = pi*deigr(1:mband,1:nkpt) 619 dedni(1:mband,1:nkpt,imod,iqpt) = pi*deigi(1:mband,1:nkpt) 620 end if 621 622 end do ! imod 623 end do !iqpt 624 625 close(ddbun) 626 627 628 !============================================================================= 629 !3) Evaluation of the Eliashberg type spectral function 630 !and phonon DOS via gaussian broadning 631 !============================================================================= 632 633 if(telphint==1)then 634 ng2f = 500 ! number of frequencies 635 omega_min=zero 636 omega_max=zero 637 do iqpt=1,nqpt 638 do imod=1,3*natom 639 omega_min = min(omega_min,phfreq(imod,iqpt)) 640 omega_max = max(omega_max,phfreq(imod,iqpt)) 641 end do 642 end do 643 644 ABI_ALLOCATE(dos_phon,(ng2f)) 645 ABI_ALLOCATE(g2f,(mband,nkpt,ng2f)) 646 ABI_ALLOCATE(tmpg2f,(mband,nkpt,ng2f)) 647 ABI_ALLOCATE(tmpphondos,(ng2f)) 648 649 write(std_out,'(a,es13.6)') 'omega_min :', omega_min 650 write(std_out,'(a,es13.6)') 'omega_max :', omega_max 651 write(std_out,'(a,i8)') 'ng2f :', ng2f 652 653 omega_max = omega_max + 0.1 * omega_max 654 domega = (omega_max-omega_min)/(ng2f-one) 655 656 gaussprefactor = sqrt(piinv) / g2fsmear 657 gaussfactor = one / g2fsmear 658 659 g2f(:,:,:) = zero 660 dos_phon(:) = zero 661 662 do iqpt=1,nqpt 663 do imod=1,3*natom 664 omega = omega_min 665 tmpg2f(:,:,:) = zero 666 tmpphondos(:) = zero 667 do iomega=1,ng2f 668 xx = (omega-phfreq(imod,iqpt))*gaussfactor 669 gaussval = gaussprefactor*exp(-xx*xx) 670 tmpg2f(:,:,iomega) = tmpg2f(:,:,iomega) + gaussval*dednr(:,:,imod,iqpt) 671 tmpphondos(iomega) = tmpphondos(iomega) + gaussval 672 omega = omega+domega 673 end do 674 675 g2f(:,:,:) = g2f(:,:,:) + tmpg2f(:,:,:) 676 dos_phon(:) = dos_phon(:) + tmpphondos(:) 677 678 end do !imod 679 end do !iqpt 680 681 dos_phon(:) = dos_phon(:) / nqpt 682 683 ! output the g2f 684 unit_g2f = 108 685 call outg2f(domega,omega_min,omega_max,elph_base_name,g2f,g2fsmear,kpnt,mband,ng2f,nkpt,nqpt,1,telphint,unit_g2f) 686 687 ! output the phonon DOS 688 unit_phdos = 108 689 call outphdos(domega,dos_phon,omega_min,omega_max,elph_base_name,g2fsmear,ng2f,nqpt,1,telphint,unit_g2f) 690 691 692 ABI_DEALLOCATE(dos_phon) 693 ABI_DEALLOCATE(g2f) 694 ABI_DEALLOCATE(tmpg2f) 695 ABI_DEALLOCATE(tmpphondos) 696 697 end if !telphint 698 699 !======================================================================= 700 !4) Evaluation of the Eliashberg type spectral function 701 !and phonon DOS via improved tetrahedron method 702 !======================================================================= 703 704 if(telphint==0)then 705 706 ! make dimension-ful rprimd and gprimd for transformation of derivatives to cartesian coordinates. 707 call mkrdim(acell,rprim,rprimd) 708 call matr3inv(rprimd,gprimd) 709 710 ! Q point Grid 711 qpt_full(:,:) = ddb%qpt(1:3,:) 712 713 ! Trivial Q point index 714 do iqpt=1,nqpt 715 indqpt(iqpt)=iqpt 716 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) 717 end do 718 719 ! Build qlatt from scratch (for 5.7) 720 tol = 0.1_dp 721 ilatt = 0 722 call sort_dp(nqpt,qptnrm,indqpt,tol) 723 724 do iqpt1=1,nqpt-2 725 mesh(1:3,1) = qpt_full(1:3,indqpt(iqpt1)) 726 do iqpt2=iqpt1+1,nqpt-1 727 mesh(1:3,2)= qpt_full(1:3,indqpt(iqpt2)) 728 do iqpt3=iqpt2+1,nqpt 729 mesh(1:3,3)= qpt_full(1:3,indqpt(iqpt3)) 730 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) & 731 & -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) 732 invdet = one/det 733 if (abs(nint(invdet))==nqpt .and. abs(invdet)-nqpt < tol) then 734 ilatt = 1 735 qlatt(:,:) = mesh(:,:) 736 exit 737 end if 738 end do 739 if(ilatt==1) exit 740 end do 741 if(ilatt==1) exit 742 end do 743 744 ! error message if qlatt not found and stop 745 if(ilatt==0) then 746 write(message, '(a,a)' ) & 747 & ' Could not find homogeneous basis vectors for Q point grid ',ch10 748 call wrtout(std_out,message,'COLL') 749 call wrtout(ab_out,message,'COLL') 750 MSG_ERROR("Aborting now") 751 end if 752 753 ! test if qlatt is righthanded and possibly fixe it 754 if(invdet < 0) then 755 tempqlatt(:) = qlatt(:,2) 756 qlatt(:,2) = qlatt(:,1) 757 qlatt(:,1) = tempqlatt(:) 758 end if 759 760 write(std_out,*) 'qlatt',qlatt 761 762 ! test if qlatt generates all Q points TO DO 763 764 765 766 ! Get tetrahedra, ie indexes of the full kpoints at their summits 767 call init_tetra(indqpt,gprimd,qlatt,qpt_full,nqpt,& 768 & tetrahedra, ierr, errstr) 769 ABI_CHECK(ierr==0,errstr) 770 771 rcvol = abs (gprimd(1,1)*(gprimd(2,2)*gprimd(3,3)-gprimd(3,2)*gprimd(2,3)) & 772 & -gprimd(2,1)*(gprimd(1,2)*gprimd(3,3)-gprimd(3,2)*gprimd(1,3)) & 773 & +gprimd(3,1)*(gprimd(1,2)*gprimd(2,3)-gprimd(2,2)*gprimd(1,3))) 774 775 ! Calculate weights for phonon DOS 776 ! Special precautions must be taking for Gamma point 777 ! because of non-analytic term. 778 ! Non-analyticity must be taken out and treated separatly. 779 780 nene = 100 !nene=number of energies for DOS 781 enemin = minval(phfreq) 782 enemax = maxval(phfreq) 783 deltaene = (enemax-enemin)/dble(nene-1) 784 ! redefine enemin enemax to be at rounded multiples of deltaene 785 ! enemin = elph_ds%fermie - dble(ifermi)*deltaene 786 ! enemax = elph_ds%fermie + dble(nene-ifermi-1)*deltaene 787 788 ABI_ALLOCATE(tweight,(nqpt,nene)) 789 ABI_ALLOCATE(dtweightde,(nqpt,nene)) 790 ABI_ALLOCATE(intweight,(3*natom,nqpt,nene)) 791 ABI_ALLOCATE(indtweightde,(3*natom,nqpt,nene)) 792 793 do iband=1,3*natom 794 eigen_in(:) = phfreq(iband,:) 795 796 ! calculate general integration weights at each irred kpoint as in Blochl et al PRB 49 16223 797 call get_tetra_weight(eigen_in,enemin,enemax,& 798 & one,nene,nqpt,tetrahedra,bcorr0,& 799 & tweight,dtweightde,xmpi_comm_self) 800 801 intweight(iband,:,:) = tweight(:,:) 802 indtweightde(iband,:,:) = dtweightde(:,:) 803 804 end do !iband 805 806 ! intdtweightse(nband,nqpt,nene) represents the weight in each energy bin for every kpt and every band 807 ! So phonon DOS is calculated (neglecting the non-analyticity contribution for now !!!) 808 809 ABI_ALLOCATE(total_dos,(nene)) 810 ABI_ALLOCATE(g2f,(mband,nkpt,nene)) 811 812 total_dos(:) = zero 813 do iband=1,3*natom 814 do iqpt=1,nqpt 815 total_dos(:) = total_dos + indtweightde(iband,iqpt,:) 816 end do 817 end do 818 819 ! For the g2f function 820 ! Right now for one electronic band and one K point: dednr(1:mband,1:nkpt,imod,iqpt) 821 ! Once again must pay close attention to the Gamma point 822 g2f(:,:,:) = zero 823 do ii=1,mband 824 do ikpt=1,nkpt 825 do iband=1,3*natom 826 do iqpt=1,nqpt 827 g2f(ii,ikpt,:) = g2f(ii,ikpt,:) + dednr(ii,ikpt,iband,iqpt) * indtweightde(iband,iqpt,:) 828 end do 829 end do 830 end do 831 end do 832 833 ! output the g2f 834 unit_g2f = 108 835 call outg2f(deltaene,enemin,enemax,elph_base_name,g2f,g2fsmear,kpnt,mband,nene,nkpt,nqpt,tetrahedra%ntetra,telphint,unit_g2f) 836 837 ! output the phonon DOS 838 unit_phdos = 108 839 call outphdos(deltaene,total_dos,enemin,enemax,elph_base_name,g2fsmear,nene,nqpt,tetrahedra%ntetra,telphint,unit_g2f) 840 841 ABI_DEALLOCATE(tweight) 842 ABI_DEALLOCATE(dtweightde) 843 ABI_DEALLOCATE(intweight) 844 ABI_DEALLOCATE(indtweightde) 845 ABI_DEALLOCATE(total_dos) 846 ABI_DEALLOCATE(g2f) 847 end if !telphint 848 849 !======================================================================= 850 !5) direct evaluation of thermal corrections 851 !======================================================================= 852 853 !open TBS file 854 outfile = trim(elph_base_name)//"_TBS" 855 if (open_file(outfile,message,newunit=unitout,form='formatted',status='unknown') /= 0) then 856 MSG_ERROR(message) 857 end if 858 write(unitout,'(a)')'thmeig: Thermal Eigenvalue corrections (eV)' 859 860 slope(:,:,:) = zero 861 zeropoint(:,:,:) = zero 862 !Loop on temperatures 863 do itemper= 1, ntemper 864 tmp=tempermin+temperinc*float(itemper-1) 865 thmeigen(:,:,:) = zero 866 867 ! Sum on all phonon wavevectors and modes 868 do iqpt=1,nqpt 869 do imod=1,3*natom 870 871 ! Bose-Einstein distribution 872 ! jmb overflow with exp(). So, select bosein to be still significant wrt half 873 if(phfreq(imod,iqpt)<tol6 .or. (phfreq(imod,iqpt)/(kb_HaK*tmp)) > -log(tol16))then 874 bosein = zero 875 else 876 bosein = one/(exp(phfreq(imod,iqpt)/(kb_HaK*tmp))-one) 877 end if 878 879 ! Calculate total 880 thmeigen(1,1:mband,1:nkpt) = thmeigen(1,1:mband,1:nkpt) + dednr(1:mband,1:nkpt,imod,iqpt)*(bosein+half) 881 thmeigen(2,1:mband,1:nkpt) = thmeigen(2,1:mband,1:nkpt) + dedni(1:mband,1:nkpt,imod,iqpt)*(bosein+half) 882 883 if(itemper==1)then 884 ! Calculate slope of linear regime 885 if(phfreq(imod,iqpt)<tol6)then 886 slope(1,1:mband,1:nkpt) = slope(1,1:mband,1:nkpt) 887 slope(2,1:mband,1:nkpt) = slope(2,1:mband,1:nkpt) 888 else 889 slope(1,1:mband,1:nkpt) = slope(1,1:mband,1:nkpt) + dednr(1:mband,1:nkpt,imod,iqpt)*(kb_HaK/phfreq(imod,iqpt)) 890 slope(2,1:mband,1:nkpt) = slope(2,1:mband,1:nkpt) + dedni(1:mband,1:nkpt,imod,iqpt)*(kb_HaK/phfreq(imod,iqpt)) 891 end if 892 ! Calculate zero-point renormalization 893 zeropoint(1,1:mband,1:nkpt) = zeropoint(1,1:mband,1:nkpt) + dednr(1:mband,1:nkpt,imod,iqpt)*half 894 zeropoint(2,1:mband,1:nkpt) = zeropoint(2,1:mband,1:nkpt) + dedni(1:mband,1:nkpt,imod,iqpt)*half 895 896 end if 897 end do ! imod 898 end do !iqpt 899 900 ! Write temperature independent results 901 if(itemper==1)then 902 write(unitout,'(a)')'Temperature independent results (zero-point renormalization and slope)' 903 do ikpt=1,nkpt 904 write(unitout,'(a,3es16.8)')' Kpt :', kpnt(:,ikpt,1) 905 do iband=1,mband 906 write(unitout,'(4d22.14)') Ha_eV*zeropoint(1,iband,ikpt),Ha_eV*zeropoint(2,iband,ikpt),& 907 & Ha_eV*slope(1,iband,ikpt),Ha_eV*slope(2,iband,ikpt) 908 end do 909 end do 910 write(unitout,'(a)')'Temperature dependent corrections' 911 end if 912 ! Write result in file for each temperature 913 write(unitout,'(a,es10.3,a)')'T :', tmp,' K' 914 do ikpt=1,nkpt 915 write(unitout,'(a,3es16.8)')' Kpt :', kpnt(:,ikpt,1) 916 do iband=1,mband 917 write(unitout,'(2d22.14)') Ha_eV*thmeigen(1,iband,ikpt), Ha_eV*thmeigen(2,iband,ikpt) 918 end do 919 end do 920 end do !itemper 921 922 close(unitout) 923 924 !Write temperature-independent results to the main output file 925 write(iout,'(a)')' ' 926 write(iout,'(80a)') ('-',ii=1,80) 927 write(iout,'(a)')' ' 928 write(iout,'(a)')' Electron-phonon change of electronic structure.' 929 write(iout,'(a)')' The temperature-dependent values are written in the _TBS file.' 930 write(iout,'(a)')' Here follows, for each electronic wavevector and band :' 931 write(iout,'(a)')' zero-point renormalisation (Ha) and linear slope (Ha/Kelvin)' 932 do ikpt=1,nkpt 933 write(iout,'(2a,i6,a,3es16.6)')ch10,' Kpt number ',ikpt,', with reduced coordinates :',kpnt(:,ikpt,1) 934 do iband=1,mband 935 write(iout,'(i6,2es20.6)') iband,zeropoint(1,iband,ikpt),slope(1,iband,ikpt) 936 end do 937 end do 938 939 ABI_DEALLOCATE(typat) 940 ABI_DEALLOCATE(atifc) 941 ABI_DEALLOCATE(zion) 942 ABI_DEALLOCATE(amu) 943 ABI_DEALLOCATE(xcart) 944 ABI_DEALLOCATE(xred) 945 ABI_DEALLOCATE(symafm) 946 ABI_DEALLOCATE(spinat) 947 ABI_DEALLOCATE(symrel) 948 ABI_DEALLOCATE(symrec) 949 ABI_DEALLOCATE(indsym) 950 ABI_DEALLOCATE(tnons) 951 ABI_DEALLOCATE(deigi) 952 ABI_DEALLOCATE(deigr) 953 ABI_DEALLOCATE(dwtermi) 954 ABI_DEALLOCATE(dwtermr) 955 ABI_DEALLOCATE(multi) 956 ABI_DEALLOCATE(multr) 957 ABI_DEALLOCATE(slope) 958 ABI_DEALLOCATE(thmeigen) 959 ABI_DEALLOCATE(zeropoint) 960 961 ABI_DEALLOCATE(dedni) 962 ABI_DEALLOCATE(dednr) 963 if(allocated(indqpt)) then 964 ABI_DEALLOCATE(indqpt) 965 end if 966 ABI_DEALLOCATE(eigen_in) 967 ABI_DEALLOCATE(qpt_full) 968 ABI_DEALLOCATE(qptnrm) 969 ABI_DEALLOCATE(wghtq) 970 ABI_DEALLOCATE(spqpt) 971 ABI_DEALLOCATE(eigvec) 972 ABI_DEALLOCATE(phfreq) 973 974 ABI_DEALLOCATE(blkval2) 975 ABI_DEALLOCATE(blkval2gqpt) 976 ABI_DEALLOCATE(kpnt) 977 ABI_DEALLOCATE(carflg_eig2) 978 979 980 981 call ddb_free(ddb_eig2) 982 983 call destroy_tetra(tetrahedra) 984 985 contains