TABLE OF CONTENTS
ABINIT/mka2f [ Functions ]
NAME
mka2f
FUNCTION
calculate the FS averaged alpha^2F function
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>=data type gathering info on the crystalline structure. Ifc<ifc_type>=Object containing the interatomic force constants. elph_ds elph_ds%gkk2 = gkk2 matrix elements on full FS grid for each phonon mode elph_ds%nbranch = number of phonon branches = 3*natom elph_ds%nFSband = number of bands included in the FS integration elph_ds%k_phon%nkpt = number of kpts included in the FS integration elph_ds%k_phon%kpt = coordinates of all FS kpoints elph_ds%k_phon%wtk = integration weights on the FS elph_ds%n0 = DOS at the Fermi level calculated from the k_phon integration weights (event. 2 spin pol) mustar = coulomb pseudopotential parameter natom = number of atoms
OUTPUT
a2f_1d = 1D alpha dos_phon = density of states for phonons elph_ds
PARENTS
elphon
CHILDREN
d2c_wtq,ep_ph_weights,ftgam,ftgam_init,gam_mult_displ,ifc_fourq phdispl_cart2red,simpson_int,wrtout,zgemm
NOTES
copied from ftiaf9.f
SOURCE
48 #if defined HAVE_CONFIG_H 49 #include "config.h" 50 #endif 51 52 #include "abi_common.h" 53 54 subroutine mka2f(Cryst,ifc,a2f_1d,dos_phon,elph_ds,kptrlatt,mustar) 55 56 use defs_basis 57 use defs_elphon 58 use m_errors 59 use m_profiling_abi 60 use m_special_funcs 61 62 use m_io_tools, only : open_file 63 use m_numeric_tools, only : simpson_int, simpson 64 use m_geometry, only : phdispl_cart2red 65 use m_crystal, only : crystal_t 66 use m_ifc, only : ifc_type, ifc_fourq 67 use m_dynmat, only : ftgam_init, ftgam 68 69 !This section has been created automatically by the script Abilint (TD). 70 !Do not modify the following lines by hand. 71 #undef ABI_FUNC 72 #define ABI_FUNC 'mka2f' 73 use interfaces_14_hidewrite 74 use interfaces_77_ddb, except_this_one => mka2f 75 !End of the abilint section 76 77 implicit none 78 79 !Arguments ------------------------------------ 80 !scalars 81 real(dp),intent(in) :: mustar 82 type(ifc_type),intent(in) :: ifc 83 type(crystal_t),intent(in) :: Cryst 84 type(elph_type),target,intent(inout) :: elph_ds 85 !arrays 86 integer, intent(in) :: kptrlatt(3,3) 87 real(dp),intent(out) :: a2f_1d(elph_ds%na2f),dos_phon(elph_ds%na2f) 88 89 !Local variables ------------------------- 90 !scalars 91 integer :: natom,iFSqpt,ibranch,iomega,nbranch,na2f,nsppol,nkpt,nrpt 92 integer :: isppol,jbranch,unit_a2f,unit_phdos,ep_scalprod 93 integer :: itemp, ntemp = 100 94 real(dp) :: temp 95 real(dp) :: a2fprefactor,avgelphg,avglambda,avgomlog,diagerr 96 real(dp) :: lambda_2,lambda_3,lambda_4,lambda_5 97 real(dp) :: spinfact 98 real(dp) :: lambda_iso(elph_ds%nsppol) 99 real(dp) :: lqn,omega 100 real(dp) :: omegalog(elph_ds%nsppol) 101 real(dp) :: omlog_qn 102 real(dp) :: tc_macmill,a2fsmear,domega,omega_min,omega_max 103 real(dp) :: gaussval, gaussprefactor, gaussfactor, gaussmaxval, xx 104 character(len=500) :: msg 105 character(len=fnlen) :: fname,base_name 106 !arrays 107 real(dp) :: displ_cart(2,elph_ds%nbranch,elph_ds%nbranch) 108 real(dp) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch) 109 real(dp) :: eigval(elph_ds%nbranch) 110 real(dp) :: gam_now(2,elph_ds%nbranch*elph_ds%nbranch) 111 real(dp) :: imeigval(elph_ds%nbranch) 112 ! real(dp) :: pheigvec(2*elph_ds%nbranch*elph_ds%nbranch),phfrq(elph_ds%nbranch) 113 real(dp) :: tmp_a2f(elph_ds%na2f) 114 real(dp) :: tmp_gam1(2,elph_ds%nbranch,elph_ds%nbranch) 115 real(dp) :: tmp_gam2(2,elph_ds%nbranch,elph_ds%nbranch) 116 real(dp) :: tmp_phondos(elph_ds%na2f),n0(elph_ds%nsppol) 117 real(dp),pointer :: kpt(:,:) 118 real(dp),allocatable :: phfrq(:,:) 119 real(dp),allocatable :: pheigvec(:,:) 120 real(dp),allocatable :: tmp_wtq(:,:,:) 121 real(dp),allocatable :: a2f1mom(:),a2f2mom(:),a2f3mom(:),a2f4mom(:) 122 real(dp),allocatable :: a2f_1mom(:),a2flogmom(:) 123 real(dp),allocatable :: a2flogmom_int(:) 124 real(dp),allocatable :: coskr(:,:) 125 real(dp),allocatable :: sinkr(:,:) 126 real(dp),allocatable :: linewidth_of_t(:) 127 real(dp),allocatable :: linewidth_integrand(:,:) 128 129 ! ********************************************************************* 130 !calculate a2f for frequencies between 0 and elph_ds%omega_max 131 132 DBG_ENTER("COLL") 133 134 !might need kptrlatt for finer interpolation later 135 ABI_UNUSED(kptrlatt(1,1)) 136 137 ! nrpt = number of real-space points for FT interpolation 138 nrpt = Ifc%nrpt 139 natom = Cryst%natom 140 141 nbranch = elph_ds%nbranch 142 na2f = elph_ds%na2f 143 nsppol = elph_ds%nsppol 144 base_name = elph_ds%elph_base_name 145 a2fsmear = elph_ds%a2fsmear 146 nkpt = elph_ds%k_phon%nkpt 147 kpt => elph_ds%k_phon%kpt 148 149 ep_scalprod = elph_ds%ep_scalprod 150 n0 = elph_ds%n0 151 152 !spinfact should be 1 for a normal non sppol calculation without spinorbit 153 !for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2 154 !for sppol 2 it should be 0.5 as we have 2 spin channels to sum 155 spinfact = one/elph_ds%nsppol !/elph_ds%nspinor 156 157 !maximum value of frequency (a grid has to be chosen for the representation of alpha^2 F) 158 !WARNING! supposes this value has been set in mkelph_linwid. 159 domega = (elph_ds%omega_max-elph_ds%omega_min)/(na2f-one) 160 elph_ds%domega = domega ! MG Why do we need to store domega in elph_ds? 161 omega_min = elph_ds%omega_min 162 omega_max = elph_ds%omega_max 163 164 gaussprefactor = sqrt(piinv) / a2fsmear 165 gaussfactor = one / a2fsmear 166 gaussmaxval = sqrt(-log(1.d-100)) 167 168 ! only open the file for the first sppol 169 fname = trim(base_name) // '_A2F' 170 if (open_file(fname,msg,newunit=unit_a2f,status="unknown") /= 0) then 171 MSG_ERROR(msg) 172 end if 173 174 !write (std_out,*) ' a2f function integrated over the FS' 175 176 !output the a2f_1d header 177 write (unit_a2f,'(a)') '#' 178 write (unit_a2f,'(a)') '# ABINIT package : a2f file' 179 write (unit_a2f,'(a)') '#' 180 write (unit_a2f,'(a)') '# a2f function integrated over the FS. omega in a.u.' 181 write (unit_a2f,'(a,I10)') '# number of kpoints integrated over : ',nkpt 182 write (unit_a2f,'(a,I10)') '# number of energy points : ',na2f 183 write (unit_a2f,'(a,E16.6,a,E16.6,a)') '# between omega_min = ',omega_min,' Ha and omega_max = ',omega_max,' Ha' 184 write (unit_a2f,'(a,E16.6)') '# and the smearing width for gaussians is ',a2fsmear 185 186 ! Open file for PH DOS 187 fname = trim(base_name) // '_PDS' 188 if (open_file(fname,msg,newunit=unit_phdos,status="replace") /= 0) then 189 MSG_ERROR(msg) 190 end if 191 192 ! output the phonon DOS header 193 write (unit_phdos,'(a)') '#' 194 write (unit_phdos,'(a)') '# ABINIT package : phonon DOS file' 195 write (unit_phdos,'(a)') '#' 196 write (unit_phdos,'(a)') '# Phonon DOS integrated over the FS. omega in a.u. EXPERIMENTAL!!!' 197 write (unit_phdos,'(a,I10)') '# number of kpoints integrated over : ',nkpt 198 write (unit_phdos,'(a,I10)') '# number of energy points : ',na2f 199 write (unit_phdos,'(a,E16.6,a,E16.6,a)')'# between omega_min = ',omega_min,' Ha and omega_max = ',omega_max,' Ha' 200 write (unit_phdos,'(a,i4,a,E16.6)') '# The DOS at Fermi level for spin ', 1, ' is ', n0(1) 201 if (nsppol==2) then 202 write (unit_phdos,'(a,i4,a,E16.6)') '# The DOS at Fermi level for spin ', 2, ' is ', n0(2) 203 end if 204 write (unit_phdos,'(a,E16.6)') '# and the smearing width for gaussians is ',a2fsmear 205 write (unit_phdos,'(a)') '#' 206 207 !Get the integration weights, using tetrahedron method or gaussian 208 ABI_ALLOCATE(tmp_wtq,(nbranch,elph_ds%k_fine%nkpt,na2f+1)) 209 ABI_ALLOCATE(elph_ds%k_fine%wtq,(nbranch,elph_ds%k_fine%nkpt,na2f)) 210 ABI_ALLOCATE(elph_ds%k_phon%wtq,(nbranch,nkpt,na2f)) 211 212 ABI_ALLOCATE(phfrq,(nbranch,elph_ds%k_fine%nkpt)) 213 ABI_ALLOCATE(pheigvec,(2*nbranch*nbranch,elph_ds%k_fine%nkpt)) 214 215 do iFSqpt=1,elph_ds%k_fine%nkpt 216 call ifc_fourq(ifc,cryst,elph_ds%k_fine%kpt(:,iFSqpt),phfrq(:,iFSqpt),displ_cart,out_eigvec=pheigvec(:,iFSqpt)) 217 end do 218 219 omega_min = omega_min - domega 220 221 call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,na2f+1,Cryst%gprimd,elph_ds%kptrlatt_fine, & 222 & elph_ds%nbranch,elph_ds%telphint,elph_ds%k_fine,tmp_wtq) 223 !call ep_ph_weights(phfrq,elph_ds%a2fsmear,omega_min,omega_max,na2f+1,Cryst%gprimd,elph_ds%kptrlatt_fine, & 224 !& elph_ds%nbranch,1,elph_ds%k_fine,tmp_wtq) 225 omega_min = omega_min + domega 226 227 do iomega = 1, na2f 228 elph_ds%k_fine%wtq(:,:,iomega) = tmp_wtq(:,:,iomega+1) 229 end do 230 ABI_DEALLOCATE(tmp_wtq) 231 232 if (elph_ds%use_k_fine == 1) then 233 call d2c_wtq(elph_ds) 234 end if 235 236 ABI_ALLOCATE(coskr, (nkpt,nrpt)) 237 ABI_ALLOCATE(sinkr, (nkpt,nrpt)) 238 call ftgam_init(Ifc%gprim, nkpt, nrpt, kpt, Ifc%rpt, coskr, sinkr) 239 240 ABI_DEALLOCATE(phfrq) 241 ABI_DEALLOCATE(pheigvec) 242 243 do isppol=1,nsppol 244 write (std_out,*) '##############################################' 245 write (std_out,*) 'mka2f : Treating spin polarization ', isppol 246 write (std_out,*) '##############################################' 247 248 ! Average of electron phonon coupling over the whole BZ 249 avgelphg = zero 250 ! MG20060607 Do the same for lambda and omega_log 251 avglambda = zero 252 avgomlog = zero 253 254 a2f_1d(:) = zero 255 dos_phon(:) = zero 256 257 ! reduce the dimenstion from fine to phon for phfrq and pheigvec 258 ABI_ALLOCATE(phfrq,(nbranch,elph_ds%k_phon%nkpt)) 259 ABI_ALLOCATE(pheigvec,(2*nbranch*nbranch,elph_ds%k_phon%nkpt)) 260 261 ! loop over qpoint in full kpt grid (presumably dense) 262 ! MG TODO : This loop can be performed using the IBZ and appropriated weights. 263 do iFSqpt=1,nkpt 264 ! 265 ! This reduced version of ftgkk supposes the kpoints have been integrated 266 ! in integrate_gamma. Do FT from real-space gamma grid to 1 qpt. 267 268 if (elph_ds%ep_int_gkk == 1) then 269 gam_now(:,:) = elph_ds%gamma_qpt(:,:,isppol,iFSqpt) 270 else 271 call ftgam(Ifc%wghatm,gam_now,elph_ds%gamma_rpt(:,:,isppol,:),natom,1,nrpt,0, & 272 & coskr(iFSqpt,:), sinkr(iFSqpt,:)) 273 end if 274 275 call ifc_fourq(ifc,cryst,kpt(:,iFSqpt),phfrq(:,iFSqpt),displ_cart,out_eigvec=pheigvec) 276 277 ! Diagonalize gamma matrix at qpoint (complex matrix). 278 279 ! if ep_scalprod==0 we have to dot in the displacement vectors here 280 if (ep_scalprod==0) then 281 282 call phdispl_cart2red(natom,Cryst%gprimd,displ_cart,displ_red) 283 284 tmp_gam2 = reshape (gam_now, (/2,nbranch,nbranch/)) 285 call gam_mult_displ(nbranch, displ_red, tmp_gam2, tmp_gam1) 286 287 do jbranch=1,nbranch 288 eigval(jbranch) = tmp_gam1(1, jbranch, jbranch) 289 imeigval(jbranch) = tmp_gam1(2, jbranch, jbranch) 290 291 if (abs(imeigval(jbranch)) > tol8) then 292 write (msg,'(a,i0,a,es16.8)')" imaginary values branch = ",jbranch,' imeigval = ',imeigval(jbranch) 293 MSG_WARNING(msg) 294 end if 295 296 end do 297 298 ! if ep_scalprod==1 we have to diagonalize the matrix we interpolated. 299 else if (ep_scalprod == 1) then 300 301 ! MJV NOTE : gam_now is being recast as a (3*natom)**2 matrix here 302 call ZGEMM ( 'N', 'N', 3*natom, 3*natom, 3*natom, cone, gam_now, 3*natom,& 303 & pheigvec, 3*natom, czero, tmp_gam1, 3*natom) 304 305 call ZGEMM ( 'C', 'N', 3*natom, 3*natom, 3*natom, cone, pheigvec, 3*natom,& 306 & tmp_gam1, 3*natom, czero, tmp_gam2, 3*natom) 307 308 diagerr = zero 309 do ibranch=1,nbranch 310 eigval(ibranch) = tmp_gam2(1,ibranch,ibranch) 311 do jbranch=1,ibranch-1 312 diagerr = diagerr + abs(tmp_gam2(1,jbranch,ibranch)) 313 end do 314 do jbranch=ibranch+1,nbranch 315 diagerr = diagerr + abs(tmp_gam2(1,jbranch,ibranch)) 316 end do 317 end do 318 319 if (diagerr > tol12) then 320 write(msg,'(a,es15.8)') 'mka2f: residual in diagonalization of gamma with phon eigenvectors: ', diagerr 321 MSG_WARNING(msg) 322 end if 323 324 else 325 write (msg,'(a,i0)')' Wrong value for ep_scalprod = ',ep_scalprod 326 MSG_BUG(msg) 327 end if 328 329 ! MG20060603MG 330 ! there was a bug in the calculation of the phonon DOS 331 ! since frequencies with small e-ph interaction were skipped inside the loop 332 ! In this new version all the frequencies (both positive and negative) are taken into account. 333 ! IDEA: it could be useful to calculate the PH-dos and the a2f 334 ! using several smearing values to perform a convergence study 335 ! Now the case ep_scalprod=1 is treated in the right way although it is not default anymore 336 ! FIXME to be checked 337 ! ENDMG 338 339 ! Add all contributions from the phonon modes at this qpoint to a2f and the phonon dos. 340 do ibranch=1,nbranch 341 342 ! if (abs(phfrq(ibranch,iFSqpt)) < tol10) then 343 if (abs(phfrq(ibranch,iFSqpt)) < tol7) then 344 a2fprefactor= zero 345 lqn = zero 346 omlog_qn = zero 347 else 348 a2fprefactor = eigval(ibranch)/(two_pi*abs(phfrq(ibranch,iFSqpt))*n0(isppol)) 349 lqn = eigval(ibranch)/(pi*phfrq(ibranch,iFSqpt)**2*n0(isppol)) 350 omlog_qn = lqn*log(abs(phfrq(ibranch,iFSqpt))) 351 end if 352 353 ! Add contribution to average elphon coupling 354 ! MANY ISSUES WITH FINITE T SUMS. THIS IS DEFINITELY 355 ! NOT A CORRECT FORMULATION YET. 356 357 ! Added avglambda and avgomglog to calculate lamda and omega_log using the sum over the kpt-grid. 358 ! If the k-grid is dense enough, these values should be better than the corresponding quantities 359 ! evaluated through the integration over omega that depends on the a2fsmear 360 361 avgelphg = avgelphg + eigval(ibranch) 362 avglambda = avglambda + lqn 363 avgomlog= avgomlog + omlog_qn 364 ! ENDMG 365 366 omega = omega_min 367 tmp_a2f(:) = zero 368 tmp_phondos(:) = zero 369 do iomega=1,na2f 370 xx = (omega-phfrq(ibranch,iFSqpt))*gaussfactor 371 omega = omega + domega 372 if (abs(xx) > gaussmaxval) cycle 373 374 gaussval = gaussprefactor*exp(-xx*xx) 375 tmp_a2f(iomega) = tmp_a2f(iomega) + gaussval*a2fprefactor 376 tmp_phondos(iomega) = tmp_phondos(iomega) + gaussval 377 end do 378 379 ! tmp_a2f(:) = zero 380 ! tmp_phondos(:) = zero 381 ! do iomega=1,na2f 382 ! tmp_a2f(iomega) = tmp_a2f(iomega) + a2fprefactor*elph_ds%k_phon%wtq(ibranch,iFSqpt,iomega) 383 ! tmp_phondos(iomega) = tmp_phondos(iomega) + elph_ds%k_phon%wtq(ibranch,iFSqpt,iomega) 384 ! end do 385 386 a2f_1d(:) = a2f_1d(:) + tmp_a2f(:) 387 dos_phon(:) = dos_phon(:) + tmp_phondos(:) 388 389 end do ! ibranch 390 end do ! iFSqpt do 391 392 393 ! second 1 / nkpt factor for the integration weights 394 a2f_1d(:) = a2f_1d(:) / nkpt 395 dos_phon(:) = dos_phon(:) / nkpt 396 397 ! MG 398 avglambda = avglambda/nkpt 399 avgomlog= avgomlog/nkpt 400 avgomlog = exp (avgomlog/avglambda) 401 write(std_out,*) ' from mka2f: for spin ', isppol 402 write(std_out,*) ' w/o interpolation lambda = ',avglambda,' omega_log= ',avgomlog 403 ! ENDMG 404 405 write (std_out,'(a,I4,a,E16.6)') '# The DOS at Fermi level for spin ',isppol,' is ',n0(isppol) 406 407 write (unit_a2f,'(a,I4,a,E16.6)') '# The DOS at Fermi level for spin ',isppol,' is ',n0(isppol) 408 write (unit_a2f,'(a)') '#' 409 410 omega = omega_min 411 do iomega=1,na2f 412 write (unit_a2f,*) omega, a2f_1d(iomega) 413 omega=omega + domega 414 end do 415 write (unit_a2f,*) 416 ! 417 ! output the phonon DOS, but only for the first sppol case 418 if (isppol == 1) then 419 omega = omega_min 420 do iomega=1,na2f 421 write (unit_phdos,*) omega, dos_phon(iomega) 422 omega=omega + domega 423 end do 424 end if 425 ! 426 ! Do isotropic calculation of lambda and output lambda, Tc(MacMillan) 427 ! 428 ABI_ALLOCATE(a2f_1mom,(na2f)) 429 ABI_ALLOCATE(a2f1mom,(na2f)) 430 ABI_ALLOCATE(a2f2mom,(na2f)) 431 ABI_ALLOCATE(a2f3mom,(na2f)) 432 ABI_ALLOCATE(a2f4mom,(na2f)) 433 ABI_ALLOCATE(linewidth_integrand,(na2f,ntemp)) 434 ABI_ALLOCATE(linewidth_of_t,(ntemp)) 435 436 a2f_1mom=zero 437 a2f1mom=zero; a2f2mom=zero 438 a2f3mom=zero; a2f4mom=zero 439 linewidth_integrand = zero 440 441 omega = omega_min 442 do iomega=1,na2f 443 if (abs(omega) > tol10) then 444 a2f_1mom(iomega) = two*spinfact*a2f_1d(iomega)/abs(omega) ! first inverse moment of alpha2F 445 a2f1mom(iomega) = two*spinfact*a2f_1d(iomega)*abs(omega) ! first positive moment of alpha2F 446 a2f2mom(iomega) = a2f1mom(iomega)*abs(omega) ! second positive moment of alpha2F 447 a2f3mom(iomega) = a2f2mom(iomega)*abs(omega) ! third positive moment of alpha2F 448 a2f4mom(iomega) = a2f3mom(iomega)*abs(omega) ! fourth positive moment of alpha2F 449 ! 450 ! electron lifetimes eq 4.48 in Grimvall electron phonon coupling in Metals (with T dependency). Also 5.69-5.72, 5.125, section 3.4 451 ! phonon lifetimes eq 19 in Savrasov PhysRevB.54.16487 (T=0) 452 ! a first T dependent expression in Allen PRB 6 2577 eq 10. Not sure about the units though 453 ! 454 do itemp = 1, ntemp 455 temp = (itemp-1)*10._dp*kb_HaK 456 linewidth_integrand(iomega, itemp) = a2f_1d(iomega) * (fermi_dirac(omega,zero,temp) + bose_einstein(omega,temp)) 457 end do 458 end if 459 omega=omega + domega 460 end do 461 ! 462 ! From Allen PRL 59 1460 463 ! \lambda <\omega^n> = 2 \int_0^{\infty} d\omega [\alpha^2F / \omega] \omega^n 464 ! 465 lambda_iso(isppol) = simpson(domega,a2f_1mom) 466 lambda_2 = simpson(domega,a2f1mom) 467 lambda_3 = simpson(domega,a2f2mom) 468 lambda_4 = simpson(domega,a2f3mom) 469 lambda_5 = simpson(domega,a2f4mom) 470 do itemp = 1, ntemp 471 linewidth_of_t(itemp) = simpson(domega,linewidth_integrand(:,itemp)) 472 ! print out gamma(T) here 473 temp = (itemp-1)*10._dp*kb_HaK 474 write (std_out,*) 'mka2f: T, average linewidth', temp, linewidth_of_t(itemp) 475 end do 476 477 478 ABI_DEALLOCATE(phfrq) 479 ABI_DEALLOCATE(pheigvec) 480 ABI_DEALLOCATE(a2f_1mom) 481 ABI_DEALLOCATE(a2f1mom) 482 ABI_DEALLOCATE(a2f2mom) 483 ABI_DEALLOCATE(a2f3mom) 484 ABI_DEALLOCATE(a2f4mom) 485 ABI_DEALLOCATE(linewidth_integrand) 486 ABI_DEALLOCATE(linewidth_of_t) 487 488 write (std_out,*) 'mka2f: elphon coupling lambdas for spin = ', isppol 489 write (std_out,*) 'mka2f: isotropic lambda', lambda_iso(isppol) 490 write (std_out,*) 'mka2f: positive moments of alpha2F:' 491 write (std_out,*) 'lambda <omega^2> = ', lambda_2 492 write (std_out,*) 'lambda <omega^3> = ', lambda_3 493 write (std_out,*) 'lambda <omega^4> = ', lambda_4 494 write (std_out,*) 'lambda <omega^5> = ', lambda_5 495 ! 496 ! Get log moment of alpha^2F 497 ABI_ALLOCATE(a2flogmom,(na2f)) 498 ABI_ALLOCATE(a2flogmom_int,(na2f)) 499 omega = omega_min 500 a2flogmom(:) = zero 501 do iomega=1,na2f 502 if (abs(omega) > tol10) then 503 a2flogmom(iomega) = a2f_1d(iomega)*log(abs(omega))/abs(omega) 504 end if 505 omega=omega + domega 506 end do 507 call simpson_int(na2f,domega,a2flogmom,a2flogmom_int) 508 509 ! NOTE: omegalog actually stores the log moment of a2F, which is the quantity to sum over spins, instead of 510 ! exp(moment/lambda) which is an actual frequency 511 omegalog(isppol) = two*spinfact*a2flogmom_int(na2f) 512 513 ABI_DEALLOCATE(a2flogmom) 514 ABI_DEALLOCATE(a2flogmom_int) 515 516 if (nsppol > 1) then 517 write (msg, '(3a)' ) ch10,& 518 & ' Warning : some of the following quantities should be integrated over spin', ch10 519 call wrtout(std_out,msg,'COLL') 520 call wrtout(ab_out,msg,'COLL') 521 end if 522 523 write (msg, '(3a)' ) ch10,& 524 & ' Superconductivity : isotropic evaluation of parameters from electron-phonon coupling.',ch10 525 call wrtout(std_out,msg,'COLL') 526 call wrtout(ab_out,msg,'COLL') 527 528 if (elph_ds%nsppol > 1) then 529 write (msg, '(a,i6,a,es16.6)' )' mka2f: isotropic lambda for spin ', isppol, ' = ', lambda_iso(isppol) 530 call wrtout(std_out,msg,'COLL') 531 call wrtout(ab_out,msg,'COLL') 532 end if 533 534 write (msg, '(a,es16.6)' )' mka2f: lambda <omega^2> = ', lambda_2 535 call wrtout(std_out,msg,'COLL') 536 call wrtout(ab_out,msg,'COLL') 537 538 write (msg, '(a,es16.6)' )' mka2f: lambda <omega^3> = ', lambda_3 539 call wrtout(std_out,msg,'COLL') 540 call wrtout(ab_out,msg,'COLL') 541 542 write (msg, '(a,es16.6)' )' mka2f: lambda <omega^4> = ', lambda_4 543 call wrtout(std_out,msg,'COLL') 544 call wrtout(ab_out,msg,'COLL') 545 546 write (msg, '(a,es16.6)' )' mka2f: lambda <omega^5> = ', lambda_5 547 call wrtout(std_out,msg,'COLL') 548 call wrtout(ab_out,msg,'COLL') 549 550 if (elph_ds%nsppol > 1) then 551 write (msg, '(a,i6,a,es16.6,a,es16.6,a)' )' mka2f: omegalog for spin ', isppol, ' = ',& 552 & exp(omegalog(isppol)/lambda_iso(isppol)), ' (Ha) ', exp(omegalog(isppol)/lambda_iso(isppol))/kb_HaK, ' (Kelvin) ' 553 call wrtout(std_out,msg,'COLL') 554 call wrtout(ab_out,msg,'COLL') 555 end if 556 557 end do ! isppol 558 559 560 561 !also print out spin-summed quantities 562 lambda_2 = sum(lambda_iso(1:elph_ds%nsppol)) 563 write (msg, '(a,es16.6)' )' mka2f: isotropic lambda = ', lambda_2 564 call wrtout(std_out,msg,'COLL') 565 call wrtout(ab_out,msg,'COLL') 566 567 omega = exp( sum(omegalog(1:elph_ds%nsppol))/lambda_2 ) 568 write (msg, '(a,es16.6,a,es16.6,a)' )' mka2f: omegalog = ', omega, ' (Ha) ', omega/kb_HaK, ' (Kelvin) ' 569 call wrtout(std_out,msg,'COLL') 570 call wrtout(ab_out,msg,'COLL') 571 572 write (msg, '(a,es16.6)' )' mka2f: input mustar = ', mustar 573 call wrtout(std_out,msg,'COLL') 574 call wrtout(ab_out,msg,'COLL') 575 576 tc_macmill = omega/1.2_dp * exp((-1.04_dp*(one+lambda_2)) / (lambda_2-mustar*(one+0.62_dp*lambda_2))) 577 write ( msg, '(a,es16.6,a,es16.6,a)')'-mka2f: MacMillan Tc = ', tc_macmill, ' (Ha) ', tc_macmill/kb_HaK, ' (Kelvin) ' 578 call wrtout(std_out,msg,'COLL') 579 call wrtout(ab_out,msg,'COLL') 580 581 close(unit=unit_a2f) 582 close(unit=unit_phdos) 583 584 ABI_DEALLOCATE(elph_ds%k_fine%wtq) 585 ABI_DEALLOCATE(elph_ds%k_phon%wtq) 586 587 ABI_DEALLOCATE(coskr) 588 ABI_DEALLOCATE(sinkr) 589 590 DBG_EXIT("COLL") 591 592 end subroutine mka2f