TABLE OF CONTENTS
- ABINIT/m_phgamma
- m_phgamma/a2fw_ee_write
- m_phgamma/a2fw_free
- m_phgamma/a2fw_get_moment
- m_phgamma/a2fw_init
- m_phgamma/a2fw_t
- m_phgamma/a2fw_tr_free
- m_phgamma/a2fw_tr_init
- m_phgamma/a2fw_tr_moment
- m_phgamma/a2fw_tr_t
- m_phgamma/a2fw_tr_write
- m_phgamma/a2fw_write
- m_phgamma/calc_dbldelta
- m_phgamma/eph_phgamma
- m_phgamma/find_ewin
- m_phgamma/phgamma_eval_qibz
- m_phgamma/phgamma_free
- m_phgamma/phgamma_init
- m_phgamma/phgamma_interp
- m_phgamma/phgamma_interp_setup
- m_phgamma/phgamma_linwid
- m_phgamma/phgamma_ncwrite
- m_phgamma/phgamma_setup_qpoint
- m_phgamma/phgamma_t
- m_phgamma/phgamma_vv_eval_qibz
- m_phgamma/phgamma_vv_interp
- m_phgamma/phgamma_vv_interp_setup
- m_phgamma/tgamma_symm
ABINIT/m_phgamma [ Modules ]
NAME
FUNCTION
Computation of phonon linewidths, isotropic superconducting properties and transport properties in metals within the LOVA approximation to the linearized Boltzmann equation.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
TODO
1) Implement restart capabilities (eph_restart) 2) Sum over the IBZ(q) on the FS instead of full BZ 3) Gaussian adaptive smearing for double delta (tetra version from libtetrabz really slow!) 4) Automatic detections of energy window, improve filtering techniques. 5) Interface with KERANGE trick 6) More examples and tutorials (using precomputed Netcd files) 7) SKW interpolation for ph linewidths and/or linear interpolation (I don't trust plain Fourier interpolation). 8) Perform more benchmarks with dense meshes to detect hotspots and memory bottlenecks 9) Test spin and SOC
SOURCE
27 #if defined HAVE_CONFIG_H 28 #include "config.h" 29 #endif 30 31 #include "abi_common.h" 32 33 module m_phgamma 34 35 use defs_basis 36 use m_abicore 37 use m_xmpi 38 use m_errors 39 use m_krank 40 use m_htetra 41 use m_tetrahedron 42 use libtetrabz 43 use m_ifc 44 use m_ebands 45 use m_fstab 46 use, intrinsic :: iso_c_binding 47 use m_nctk 48 use netcdf 49 use m_wfk 50 use m_ddb 51 use m_ddk 52 use m_dvdb 53 use m_crystal 54 use m_fft 55 use m_hamiltonian 56 use m_pawcprj 57 use m_dtset 58 use m_dtfil 59 use m_wfd 60 use m_ephtk 61 use m_mkffnl 62 63 use defs_abitypes, only : mpi_type 64 use m_time, only : cwtime, cwtime_report 65 use m_fstrings, only : toupper, itoa, sjoin, ktoa, ltoa, strcat 66 use m_numeric_tools, only : arth, wrap2_pmhalf, simpson_int, simpson, mkherm, get_diag, isdiagmat 67 use m_io_tools, only : open_file, iomode_from_fname 68 use m_symtk, only : littlegroup_q 69 use m_geometry, only : normv 70 use m_special_funcs, only : gaussian 71 use m_fftcore, only : ngfft_seq, get_kg 72 use m_cgtools, only : cg_zdotc 73 use m_kg, only : getph, mkkpg 74 use m_dynmat, only : symdyma, ftgam_init, ftgam, asrif9 75 use defs_datatypes, only : ebands_t, pseudopotential_type 76 use m_bz_mesh, only : kpath_t, kpath_new 77 use m_special_funcs, only : fermi_dirac 78 use m_kpts, only : kpts_ibz_from_kptrlatt, tetra_from_kptrlatt, listkk, kpts_timrev_from_kptopt, kpts_map 79 use defs_elphon, only : complete_gamma !, complete_gamma_tr 80 use m_getgh1c, only : getgh1c, rf_transgrid_and_pack, getgh1c_setup 81 use m_pawang, only : pawang_type 82 use m_pawrad, only : pawrad_type 83 use m_pawtab, only : pawtab_type 84 use m_pawfgr, only : pawfgr_type 85 86 implicit none 87 88 private
m_phgamma/a2fw_ee_write [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
a2fw_ee_write
FUNCTION
Write alpha^2F(e,e',w) to an external file in text form
INPUTS
a2f<a2fw_t>=Container storing the Eliashberg functions. basename=Filename for output.
OUTPUT
Output is written to file. This routine should be called by one MPI proc.
SOURCE
2425 subroutine a2fw_ee_write(a2f, basename) 2426 2427 !Arguments ------------------------------------ 2428 !scalars 2429 character(len=*),intent(in) :: basename 2430 type(a2fw_t),intent(in) :: a2f 2431 2432 !Local variables ------------------------- 2433 !scalars 2434 integer :: iw, spin, unt, ii, iene, jene 2435 real(dp) :: ene1, ene2 2436 character(len=500) :: msg 2437 character(len=fnlen) :: path 2438 2439 ! ********************************************************************* 2440 2441 ! Write spin-resolved a2F(e,e',w) 2442 path = strcat(basename, "_A2FEEW") 2443 if (open_file(path, msg, newunit=unt, form="formatted", action="write", status="unknown") /= 0) then 2444 ABI_ERROR(msg) 2445 end if 2446 2447 call write_a2fw_header() 2448 2449 write(unt,'(a)')"# en2, en1, Frequency, a2F_tot(w) (presently summed over spin)" 2450 do iw=1,a2f%nomega 2451 do iene=1,a2f%nene 2452 ene1 = a2f%enemin + (iene-1)*a2f%deltaene 2453 do jene=1,a2f%nene 2454 ene2 = a2f%enemin + (jene-1)*a2f%deltaene 2455 write(unt,'(3E20.10,2x,E20.10)') ene2, ene1, a2f%omega(iw), sum(a2f%vals_ee(jene,iene,iw,:)) ! TOT 2456 end do 2457 end do 2458 end do 2459 2460 close(unt) 2461 2462 ! Write spin-resolved a2F(ef, ef, w) 2463 path = strcat(basename, "_A2FW_reference") 2464 if (open_file(path, msg, newunit=unt, form="formatted", status="unknown") /= 0) then 2465 ABI_ERROR(msg) 2466 end if 2467 2468 call write_a2fw_header() 2469 2470 write(unt,'(a)')"# Frequency, a2F_tot(ef,ef,w) for comparison with normal a2F(w) (presently summed over spin)" 2471 iene = int(a2f%nene/2) 2472 jene = int(a2f%nene/2) 2473 do iw=1,a2f%nomega 2474 write(unt,'(E20.10,2x,E20.10)') a2f%omega(iw), sum(a2f%vals_ee(jene,iene,iw,:)) ! TOT 2475 end do 2476 2477 close(unt) 2478 2479 contains 2480 2481 subroutine write_a2fw_header() 2482 2483 ! Output the header. 2484 write(unt,'(a)') '#' 2485 write(unt,'(a)') '# ABINIT package: a2F(e,eprime,w) file' 2486 write(unt,'(a)') '#' 2487 write(unt,'(a)') '# a2F(e,eprime,w) function integrated over the FS. omega in a.u.' 2488 write(unt,'(2a)') '# ngqpt: ',trim(ltoa(a2f%ngqpt)) 2489 write(unt,'(a,i0,2(a,e16.6))')'# number of energies: ',a2f%nene," from enemin: ",a2f%enemin,' Ha with step ', a2f%deltaene 2490 write(unt,'(a,i0,3(a,e16.6))')'# number of frequencies: ',a2f%nomega," between omega_min: ",a2f%omega_min, & 2491 ' Ha and omega_max: ',a2f%omega_max,' Ha with step:',a2f%wstep 2492 write(unt,'(a,e16.6)') '# the smearing width for gaussians is ',a2f%smear 2493 write(unt,'(a,e16.6)')"# Total DOS at Fermi level ",sum(a2f%n0) 2494 do spin=1,a2f%nsppol 2495 write(unt,"(a,i0,a,e16.6)")"# The DOS at Fermi level for spin ",spin," is ",a2f%n0(spin) 2496 end do 2497 do ii=1,2 2498 write(unt,'(a)') "# " 2499 end do 2500 2501 end subroutine write_a2fw_header 2502 2503 end subroutine a2fw_ee_write
m_phgamma/a2fw_free [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
a2fw_free
FUNCTION
Free the memory allocated in a2f
SIDE EFFECTS
a2f<a2fw_t>=Structure storing the Eliashberg function a2F.
OUTPUT
SOURCE
1633 subroutine a2fw_free(a2f) 1634 1635 !Arguments ------------------------------------ 1636 class(a2fw_t),intent(inout) :: a2f 1637 1638 ! ********************************************************************* 1639 1640 ! integer 1641 ABI_SFREE(a2f%qshift) 1642 1643 ! real 1644 ABI_SFREE(a2f%n0) 1645 ABI_SFREE(a2f%omega) 1646 ABI_SFREE(a2f%vals) 1647 ABI_SFREE(a2f%vals_ee) 1648 ABI_SFREE(a2f%lambdaw) 1649 1650 end subroutine a2fw_free
m_phgamma/a2fw_get_moment [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
a2fw_get_moment
FUNCTION
Compute \int dw [a2F(w)/w] w^n From Allen PRL 59 1460 [[cite:Allen1987]]. See also [[cite:Grimvall1981]], Eq 6.72 page 175)
INPUTS
a2f<a2fw_t>=Structure storing the Eliashberg function. nn=Value of n spin=The spin component. 0 to sum over spins.
OUTPUT
a2fw_get_moment = \int dw [a2F(w)/w] w^n
[out_int(x)] = \int^{x} dw [a2F(w)/w] w^n
SOURCE
2139 real(dp) function a2fw_get_moment(a2f, nn, spin, out_int) 2140 2141 !Arguments ------------------------------------ 2142 !scalars 2143 integer,intent(in) :: spin,nn 2144 class(a2fw_t),intent(in) :: a2f 2145 !arrays 2146 real(dp),intent(out),optional :: out_int(a2f%nomega) 2147 2148 !Local variables ------------------------- 2149 !scalars 2150 integer :: iw 2151 real(dp) :: omg, omg_nm1 2152 !arrays 2153 real(dp) :: ff(a2f%nomega), int_ff(a2f%nomega), values(a2f%nomega) 2154 2155 ! ********************************************************************* 2156 2157 ! Construct the integrand function. [a2F(w)/w] w^n 2158 ff = zero; int_ff = zero 2159 2160 if (spin == 0) then 2161 if (a2f%nsppol == 1) then 2162 values = a2f%vals(:, 0, 1) 2163 !if (a2f%nspinor == 1) values = two * values 2164 else if (a2f%nsppol == 2) then 2165 values = sum(a2f%vals(:, 0, :), dim=2) 2166 end if 2167 2168 else 2169 values = a2f%vals(:, 0, spin) 2170 end if 2171 2172 if (nn - 1 >= 0) then 2173 do iw=1,a2f%nomega 2174 omg = a2f%omega(iw) 2175 omg_nm1 = omg ** (nn - 1) 2176 ff(iw) = values(iw) * omg_nm1 2177 end do 2178 else 2179 do iw=1,a2f%nomega 2180 omg = a2f%omega(iw) 2181 omg_nm1 = zero; if (abs(omg) > EPHTK_WTOL) omg_nm1 = omg**(nn-1) 2182 ff(iw) = values(iw) * omg_nm1 2183 end do 2184 end if 2185 2186 ! Integration with simpson rule on a linear mesh. 2187 call simpson_int(a2f%nomega, a2f%wstep, ff, int_ff) 2188 2189 a2fw_get_moment = int_ff(a2f%nomega) 2190 if (present(out_int)) out_int = int_ff 2191 2192 end function a2fw_get_moment
m_phgamma/a2fw_init [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
a2fw_init
FUNCTION
Calculates the FS averaged alpha^2F(w) function
INPUTS
cryst<crystal_t>=Info on the unit cell. ifc<ifc_type>=Interatomic force constants. gams<phgamma_t>=Structure storing the phonon linewidths. wstep=Step for linear frequency mesh in Ha. wminmax(2)=Minimum and maximum phonon frequency. Used to construct the linear mesh for A2F(w). ph_intmeth=Integration method for phonons: 1 for gaussian, 2 for tetrahedra smear=Gaussian broadening used to approximate the Dirac delta. ngqpt(3)=Divisions of the Q-mesh used for interpolating the phonon linewidths (see also nqshift and qshift). nqshift=Number of shifts used to generated the Q-mesh. qshift(3,nqshift)=The shifts. comm=MPI communicator [qintp]=If set to False, ngqgpt, nqshift, qshift and qptop are ignored and A2F(w) is computed from the IBZ values stored in gams. Default: True i.e use Fourier interpolation. [qptopt]=Controls the generation of the q-points. If not specified, the routine takes fully into account the symmetries of the system to generate the q points in the IBZone i.e. qptopt=1 Other values of qptopt can be used for debugging purpose.
OUTPUT
a2f<a2fw_t>=Structure storing the Eliashberg function a2F(w).
SOURCE
1685 subroutine a2fw_init(a2f, gams, cryst, ifc, ph_intmeth, wstep, wminmax, smear, ngqpt, nqshift, qshift, comm, & 1686 qintp, qptopt) ! optional 1687 1688 !Arguments ------------------------------------ 1689 !scalars 1690 integer,intent(in) :: ph_intmeth,nqshift,comm 1691 integer,intent(in),optional :: qptopt 1692 real(dp),intent(in) :: wstep,smear 1693 logical,optional,intent(in) :: qintp 1694 type(phgamma_t),intent(inout) :: gams 1695 type(ifc_type),intent(in) :: ifc 1696 type(crystal_t),intent(in) :: cryst 1697 type(a2fw_t),target,intent(out) :: a2f 1698 !arrays 1699 integer,intent(in) :: ngqpt(3) 1700 real(dp),intent(in) :: wminmax(2),qshift(3,nqshift) 1701 1702 !Local variables ------------------------- 1703 !scalars 1704 integer,parameter :: master = 0 1705 integer :: my_qptopt,iq_ibz,nqibz,ount,my_rank,nproc,cnt 1706 integer :: mu,iw,natom3,nsppol,spin,ierr,nomega,nqbz 1707 integer :: iene, jene, itemp, ntemp, jene_jump 1708 real(dp) :: cpu,wall,gflops 1709 real(dp) :: lambda_iso,omega,omega_log,xx,omega_min,omega_max,ww,mustar,tc_macmill 1710 real(dp) :: temp_el, min_temp, delta_temp, chempot, ene1, ene2, G0 1711 logical :: do_qintp 1712 character(len=500) :: msg 1713 type(htetra_t) :: qtetra 1714 !arrays 1715 integer :: qptrlatt(3,3),new_qptrlatt(3,3) 1716 real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom) 1717 real(dp) :: phfrq(gams%natom3),gamma_ph(gams%natom3),lambda_ph(gams%natom3) 1718 real(dp) :: invphfrq(gams%natom3) 1719 real(dp),allocatable :: my_qshift(:,:), gamma_ph_ee(:,:,:,:), tmp_a2f(:) 1720 real(dp), ABI_CONTIGUOUS pointer :: a2f_1d(:) 1721 real(dp),allocatable :: qibz(:,:),wtq(:),qbz(:,:) 1722 real(dp),allocatable :: a2f_1mom(:),a2flogmom(:),a2flogmom_int(:),wdt(:,:) 1723 real(dp),allocatable :: lambda_tetra(:,:,:),phfreq_tetra(:,:), tmp_gaussian(:,:) 1724 real(dp), allocatable :: a2feew_partial(:), a2feew_partial_int(:), a2feew_w(:), a2feew_w_int(:) 1725 1726 ! ********************************************************************* 1727 1728 call cwtime(cpu, wall, gflops, "start") 1729 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 1730 1731 my_qptopt = 1; if (present(qptopt)) my_qptopt = qptopt 1732 do_qintp = .True.; if (present(qintp)) do_qintp = qintp 1733 nsppol = gams%nsppol; natom3 = gams%natom3 1734 1735 if (gams%prteliash == 3) then 1736 ABI_MALLOC(gamma_ph_ee, (gams%nene, gams%nene, gams%natom3, gams%nsppol)) 1737 end if 1738 1739 if (do_qintp) then 1740 ! Generate fine q-mesh, find the IBZ and the corresponding weights. 1741 qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3) 1742 1743 call kpts_ibz_from_kptrlatt(cryst, qptrlatt, my_qptopt, nqshift, qshift, nqibz, qibz, wtq, nqbz, qbz, & 1744 new_kptrlatt=new_qptrlatt, new_shiftk=my_qshift) 1745 ABI_FREE(qbz) 1746 1747 ! Store quantities that cannot be easily (and safely) calculated if we only know the IBZ. 1748 a2f%ngqpt = ngqpt; a2f%nqshift = size(my_qshift, dim=2) 1749 ABI_MALLOC(a2f%qshift, (3, a2f%nqshift)) 1750 a2f%qshift = my_qshift 1751 ABI_FREE(my_qshift) 1752 1753 else 1754 ! No interpolation. Use q-mesh parameters from gams% 1755 a2f%ngqpt = gams%ngqpt; a2f%nqshift = 1 1756 nqibz = gams%nqibz 1757 ABI_MALLOC(qibz, (3, nqibz)) 1758 ABI_MALLOC(wtq, (nqibz)) 1759 ABI_MALLOC(a2f%qshift, (3, a2f%nqshift)) 1760 qibz = gams%qibz; wtq = gams%wtq 1761 a2f%qshift = zero ! Note: assuming q-mesh centered on Gamma. 1762 end if 1763 1764 call cwtime_report(" a2fw_init, q-setup:", cpu, wall, gflops) 1765 1766 ! Define Min and max frequency for the mesh (enlarge it a bit) 1767 omega_min = wminmax(1); omega_max = wminmax(2) 1768 omega_min = omega_min - 0.1*abs(omega_min) 1769 if (omega_min >= zero) omega_min = one/Ha_meV 1770 omega_max = omega_max + 0.1*abs(omega_max) 1771 1772 a2f%nsppol = nsppol; a2f%natom3 = gams%natom3; a2f%smear = smear 1773 a2f%omega_min = omega_min; a2f%omega_max = omega_max 1774 nomega = int((omega_max - omega_min) / wstep); a2f%nomega = nomega; a2f%wstep = wstep 1775 a2f%nene = gams%nene 1776 a2f%enemin = gams%enemin 1777 a2f%deltaene = gams%deltaene 1778 1779 ABI_MALLOC(a2f%n0, (nsppol)) 1780 a2f%n0 = gams%n0 1781 ! Build linear mesh. 1782 ABI_MALLOC(a2f%omega, (nomega)) 1783 a2f%omega = arth(omega_min, wstep, nomega) 1784 1785 ABI_CALLOC(a2f%vals, (nomega,0:natom3, nsppol)) 1786 ABI_CALLOC(a2f%lambdaw, (nomega,0:natom3, nsppol)) 1787 1788 if (gams%prteliash == 3) then 1789 ABI_CALLOC(a2f%vals_ee, (gams%nene, gams%nene, nomega, nsppol)) 1790 end if 1791 1792 ABI_MALLOC(tmp_a2f, (nomega)) 1793 1794 if (ph_intmeth == 2) then 1795 ! Prepare tetrahedron integration. 1796 call cwtime(cpu, wall, gflops, "start") 1797 1798 qptrlatt = 0; qptrlatt(1, 1) = a2f%ngqpt(1); qptrlatt(2, 2) = a2f%ngqpt(2); qptrlatt(3, 3) = a2f%ngqpt(3) 1799 qtetra = tetra_from_kptrlatt(cryst, my_qptopt, qptrlatt, a2f%nqshift, a2f%qshift, nqibz, qibz, comm, msg, ierr) 1800 if (ierr /= 0) ABI_ERROR(msg) 1801 1802 ABI_CALLOC(lambda_tetra, (nqibz, natom3, nsppol)) 1803 ABI_CALLOC(phfreq_tetra, (nqibz, natom3)) 1804 cnt = 0 1805 do iq_ibz = 1, nqibz 1806 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism. 1807 call ifc%fourq(cryst, qibz(:,iq_ibz), phfrq, displ_cart) 1808 ! save for tetrahedron interpolation 1809 phfreq_tetra(iq_ibz,:) = phfrq(:) 1810 end do 1811 call xmpi_sum(phfreq_tetra, comm, ierr) 1812 call cwtime_report(" a2fw_init%tetra", cpu, wall, gflops) 1813 end if 1814 1815 call cwtime(cpu, wall, gflops, "start") 1816 1817 ! DEV_MJV 1818 !open (unit=900, file="a2fvals_ee.dat") 1819 !write (900,*) '# do_qintp ', do_qintp 1820 1821 ! Loop over spins and q-points in the IBZ. For the moment parallelize over iq_ibz 1822 do spin=1,nsppol 1823 cnt = 0 1824 do iq_ibz=1,nqibz 1825 ! TODO: for the moment the memory is not distributed, only the calculation exception is vals_ee 1826 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism. 1827 1828 ! Interpolate or evaluate gamma directly. 1829 if (do_qintp) then 1830 if (gams%prteliash == 3) then 1831 call gams%interp(cryst,ifc,spin,qibz(:,iq_ibz),phfrq,gamma_ph,lambda_ph,displ_cart,gamma_ph_ee=gamma_ph_ee(:,:,:,spin)) 1832 else 1833 call gams%interp(cryst,ifc,spin,qibz(:,iq_ibz),phfrq,gamma_ph,lambda_ph,displ_cart) 1834 end if 1835 else 1836 if (gams%prteliash == 3) then 1837 call gams%eval_qibz(cryst,ifc,iq_ibz,spin,phfrq,gamma_ph,lambda_ph,displ_cart,gamma_ph_ee=gamma_ph_ee(:,:,:,spin)) 1838 else 1839 call gams%eval_qibz(cryst,ifc,iq_ibz,spin,phfrq,gamma_ph,lambda_ph, displ_cart) 1840 end if 1841 end if 1842 1843 select case (ph_intmeth) 1844 case (1) 1845 ! Gaussian: Add all contributions from the phonon modes at this qpoint to a2f 1846 ! (note that unstable modes are included). 1847 ABI_MALLOC(tmp_gaussian, (nomega, natom3)) 1848 do mu=1,natom3 1849 tmp_a2f = zero 1850 do iw=1,nomega 1851 xx = a2f%omega(iw) - phfrq(mu) 1852 tmp_gaussian(iw, mu) = gaussian(xx, smear) 1853 tmp_a2f(iw) = tmp_a2f(iw) + tmp_gaussian(iw,mu) * lambda_ph(mu) * abs(phfrq(mu)) 1854 end do 1855 a2f%vals(:,mu,spin) = a2f%vals(:,mu,spin) + tmp_a2f * wtq(iq_ibz) 1856 end do 1857 1858 if (gams%prteliash == 3) then 1859 ! reset phfrq for low freq modes 1860 invphfrq = zero 1861 do mu=1,natom3 1862 if (abs(phfrq(mu)) > EPHTK_WTOL) invphfrq(mu) = one / abs(phfrq(mu)) 1863 end do 1864 1865 do iene= 1, gams%nene 1866 do jene= 1, gams%nene 1867 tmp_a2f = zero 1868 ! TODO: following block is just a GEMM 1869 do mu=1,natom3 1870 do iw=1,nomega 1871 tmp_a2f(iw) = tmp_a2f(iw) + tmp_gaussian(iw,mu) * gamma_ph_ee(jene,iene,mu,spin) * invphfrq(mu) 1872 end do 1873 end do 1874 1875 a2f%vals_ee(jene, iene, :, spin) = a2f%vals_ee(jene, iene, :, spin) + tmp_a2f(:) * wtq(iq_ibz) 1876 !if (iene == gams%nene/2 .and. jene == gams%nene/2) then 1877 ! write (900, '(a,E20.10,2x,2x,I6,3E20.10)') '#', wtq(iq_ibz), iq_ibz, invphfrq(1:3) 1878 ! do iw=1,nomega 1879 ! write (900, '(i6,2x,E20.10,2x,3E20.10,2x,3E20.10)') iw, a2f%vals_ee(jene, iene, iw, spin), & 1880 ! gamma_ph_ee(jene,iene,:,spin), tmp_gaussian(iw,1:3) 1881 ! end do 1882 ! write (900,*) 1883 !end if 1884 end do 1885 end do 1886 end if 1887 1888 ABI_FREE(tmp_gaussian) 1889 1890 case (2) 1891 ! Tetra: store data. 1892 do mu=1,natom3 1893 lambda_tetra(iq_ibz, mu, spin) = lambda_ph(mu) * abs(phfrq(mu)) 1894 end do 1895 1896 if (gams%prteliash == 3) then 1897 ABI_ERROR("Eliashberg function with tetra not coded") 1898 end if 1899 1900 case default 1901 ABI_ERROR(sjoin("Wrong ph_intmeth:", itoa(ph_intmeth))) 1902 end select 1903 1904 end do ! iq_ibz 1905 end do ! spin 1906 1907 !DEV_MJV 1908 !close(900) 1909 1910 if (ph_intmeth == 2) then 1911 ! workspace for tetra. 1912 ABI_MALLOC(wdt, (nomega, 2)) 1913 1914 ! For each mode get its contribution 1915 do spin=1,nsppol 1916 do mu=1,natom3 1917 cnt = 0 1918 do iq_ibz=1,nqibz 1919 ! NB: if we are interpolating the gamma, nqibz > gams%nqibz 1920 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism. 1921 1922 call qtetra%get_onewk(iq_ibz, gams%bcorr, nomega, nqibz, phfreq_tetra(:,mu), omega_min, omega_max, one, wdt) 1923 wdt = wdt * wtq(iq_ibz) 1924 1925 ! Accumulate (Integral of a2F is computed afterwards) 1926 a2f%vals(:,mu,spin) = a2f%vals(:,mu,spin) + wdt(:,1) * lambda_tetra(iq_ibz,mu,spin) 1927 !a2f%lambdaw(:,mu,spin) = a2f%lambdaw(:,mu,spin) + wdt(:,2) * lambda_tetra(iq_ibz, mu, spin) 1928 end do 1929 end do 1930 end do 1931 1932 ! Free memory allocated for tetra. 1933 ABI_FREE(wdt) 1934 ABI_FREE(lambda_tetra) 1935 ABI_FREE(phfreq_tetra) 1936 call qtetra%free() 1937 end if 1938 1939 ! Collect final results on each node 1940 call xmpi_sum(a2f%vals, comm, ierr) 1941 do spin=1,nsppol 1942 a2f%vals(:,0,spin) = sum(a2f%vals(:,1:natom3,spin), dim=2) 1943 ! previously would divide by g(eF, spin) 1944 !a2f%vals(:,:,spin) = a2f%vals(:,:,spin) / (two_pi*a2f%n0(spin)) 1945 end do 1946 1947 if (gams%prteliash == 3) then 1948 ! For the moment vals_ee only works with gaussians 1949 call xmpi_sum(a2f%vals_ee, comm, ierr) 1950 do spin=1,nsppol 1951 a2f%vals_ee(:,:,:,spin) = a2f%vals_ee(:,:,:,spin) / (two * pi * gams%n0(spin)) 1952 end do 1953 end if 1954 1955 call cwtime_report(" a2fw_init, a2f_eval:", cpu, wall, gflops) 1956 1957 ABI_MALLOC(a2f_1mom, (nomega)) 1958 ABI_MALLOC(a2flogmom, (nomega)) 1959 ABI_MALLOC(a2flogmom_int, (nomega)) 1960 1961 !call a2fw_print_info() 1962 1963 ! Compute lambda(w) resolved in spin and phonon mode. 1964 ! lambda(w) = int a2F(w')/w' dw' 1965 do spin=1,nsppol 1966 do mu=0,natom3 1967 1968 do iw=1,nomega 1969 ww = a2f%omega(iw) 1970 if (abs(ww) > EPHTK_WTOL) then 1971 a2f_1mom(iw) = a2f%vals(iw, mu, spin) / abs(ww) 1972 else 1973 a2f_1mom(iw) = zero 1974 end if 1975 end do 1976 1977 call simpson_int(nomega, wstep, a2f_1mom, a2f%lambdaw(:, mu, spin)) 1978 end do 1979 end do 1980 1981 ! Logarithmic moment of alpha^2F: exp((2/\lambda) \int dw a2F(w) ln(w)/w) 1982 lambda_iso = a2f%get_moment(0, 0) 1983 1984 ! Get log moment of alpha^2F. 1985 a2flogmom = zero 1986 do spin=1,nsppol 1987 a2f_1d => a2f%vals(:,0,spin) 1988 do iw=1,nomega 1989 omega = a2f%omega(iw) 1990 if (abs(omega) > EPHTK_WTOL) then 1991 a2flogmom(iw) = a2flogmom(iw) + a2f_1d(iw) * log(abs(omega)) / abs(omega) 1992 end if 1993 end do 1994 end do 1995 call simpson_int(nomega, wstep, a2flogmom, a2flogmom_int) 1996 omega_log = exp((one / lambda_iso) * a2flogmom_int(nomega)) 1997 1998 mustar = 0.12 1999 tc_macmill = omega_log/1.2_dp * exp((-1.04_dp*(one+lambda_iso)) / (lambda_iso-mustar*(one+0.62_dp*lambda_iso))) 2000 2001 if (my_rank == master) then 2002 ount = std_out 2003 write(ount,'(a,es16.6)')' isotropic new_lambda = ',lambda_iso 2004 write(ount,'(a,es16.6,a,es16.6,a)' )' new_omegalog = ',omega_log,' (Ha) ', omega_log * Ha_K, ' (Kelvin) ' 2005 write(ount,'(a,es16.6,a,es16.6,a)')' MacMillan Tc = ',tc_macmill,' (Ha) ', tc_macmill * Ha_K, ' (Kelvin) ' 2006 end if 2007 2008 #if 1 2009 ! print log moments of the alpha^2 F functions. 2010 do spin=1,nsppol 2011 a2f_1d => a2f%vals(:,0,spin) 2012 2013 lambda_iso = a2f%get_moment(0, spin) 2014 2015 ! Get log moment of alpha^2F. 2016 a2flogmom = zero 2017 do iw=1,nomega 2018 omega = a2f%omega(iw) 2019 if (abs(omega) > EPHTK_WTOL) then 2020 a2flogmom(iw) = (two / lambda_iso) * a2f_1d(iw) * log(abs(omega)) / abs(omega) 2021 ! I think this is the correct expression. 2022 !a2flogmom(iw) = (one / lambda_iso) * a2f_1d(iw) * log(abs(omega)) / abs(omega) 2023 end if 2024 end do 2025 call simpson_int(nomega, wstep, a2flogmom, a2flogmom_int) 2026 omega_log = exp(a2flogmom_int(nomega)) 2027 2028 mustar = 0.12 2029 tc_macmill = omega_log/1.2_dp * exp((-1.04_dp*(one+lambda_iso)) / (lambda_iso-mustar*(one+0.62_dp*lambda_iso))) 2030 2031 if (my_rank == master) then 2032 ount = std_out 2033 if (nsppol > 1) then 2034 write(msg,'(3a)') ch10,'Warning: some of the following quantities should be integrated over spin', ch10 2035 call wrtout(ount, msg) 2036 end if 2037 2038 if (do_qintp) then 2039 write(ount,'(a)')' Superconductivity: isotropic evaluation of parameters from electron-phonon coupling (interpolated).' 2040 else 2041 write(ount,'(a)')' Superconductivity: isotropic evaluation of parameters from electron-phonon coupling (coarse grid).' 2042 endif 2043 write(ount,'(a,es16.6)')' isotropic lambda = ',lambda_iso 2044 write(ount,'(a,es16.6,a,es16.6,a)' )' omegalog = ',omega_log,' (Ha) ', omega_log * Ha_K, ' (Kelvin) ' 2045 write(ount,'(a,es16.6,a,es16.6,a)')' MacMillan Tc = ',tc_macmill,' (Ha) ', tc_macmill * Ha_K, ' (Kelvin) ' 2046 write(ount,"(a)")' positive moments of alpha2F:' 2047 write(ount,'(a,es16.6)' )' lambda <omega^2> = ',a2f%get_moment(2, spin) 2048 write(ount,'(a,es16.6)' )' lambda <omega^3> = ',a2f%get_moment(3, spin) 2049 write(ount,'(a,es16.6)' )' lambda <omega^4> = ',a2f%get_moment(4, spin) 2050 write(ount,'(a,es16.6)' )' lambda <omega^5> = ',a2f%get_moment(5, spin) 2051 end if 2052 end do 2053 #endif 2054 2055 ! Calculate the temperature dependence of the a2f(e,e',w) integrals (G_0(T_e) 2056 ! as in PRL 110 016405 (2013) [[cite:Arnaud2013]]) 2057 if (gams%prteliash == 3 .and. my_rank == master) then 2058 ntemp = 100 2059 min_temp = zero 2060 delta_temp = 40._dp ! Kelvin 2061 ABI_MALLOC(a2feew_partial, (a2f%nene)) 2062 ABI_MALLOC(a2feew_partial_int, (a2f%nene)) 2063 ABI_MALLOC(a2feew_w, (nomega)) 2064 ABI_MALLOC(a2feew_w_int, (nomega)) 2065 2066 if (open_file("EPC_strength_aafo_T.dat", msg, newunit=ount, form="formatted", action="write", status="unknown") /= 0) then 2067 ABI_ERROR(msg) 2068 end if 2069 2070 write(ount, "(a)")"# temp_el, G_0(T_e) in W/m^3/K, spin" 2071 do spin=1,nsppol 2072 do itemp = 1, ntemp 2073 temp_el = min_temp + (itemp-1)*delta_temp 2074 2075 ! TODO: need to evolve the chemical potential with T, but I do not have access to the full DOS here!! 2076 ! possible fix using local information on DOS variation near E_F... 2077 chempot = a2f%enemin + half * a2f%nene * a2f%deltaene 2078 2079 do iw=1,nomega 2080 omega = a2f%omega(iw) 2081 jene_jump = nint(omega/a2f%deltaene) 2082 a2feew_partial = zero 2083 do iene=1,a2f%nene 2084 ene1 = a2f%enemin + (iene-1)*a2f%deltaene 2085 ene2 = ene1 + omega 2086 a2feew_partial(iene) = a2f%vals_ee(min(a2f%nene,iene+jene_jump), iene, iw, spin) * & 2087 (fermi_dirac(ene1, chempot, temp_el/Ha_K) - fermi_dirac(ene2, chempot, temp_el/Ha_K)) 2088 end do 2089 call simpson_int(a2f%nene, a2f%deltaene, a2feew_partial, a2feew_partial_int) 2090 a2feew_w(iw) = a2feew_partial_int(a2f%nene) 2091 end do 2092 call simpson_int(nomega, wstep, a2feew_w, a2feew_w_int) 2093 G0 = a2feew_w_int(nomega) * two_pi * a2f%n0(spin) / cryst%ucvol 2094 ! conversion factor for G0 to SI units = Ha_J / Time_Sec / (Bohr_meter)**3 ~ 1.2163049915755545e+30 2095 write(ount, "(2(e20.10,2x),i5)") temp_el, G0 * kb_HaK / Time_Sec / (Bohr_meter)**3, spin !* Ha_J??? 2096 end do 2097 end do 2098 close(ount) 2099 2100 ABI_FREE(a2feew_partial) 2101 ABI_FREE(a2feew_partial_int) 2102 ABI_FREE(a2feew_w) 2103 ABI_FREE(a2feew_w_int) 2104 ABI_FREE(gamma_ph_ee) 2105 end if 2106 2107 ABI_FREE(tmp_a2f) 2108 ABI_FREE(a2f_1mom) 2109 ABI_FREE(a2flogmom) 2110 ABI_FREE(a2flogmom_int) 2111 ABI_FREE(qibz) 2112 ABI_FREE(wtq) 2113 2114 end subroutine a2fw_init
m_phgamma/a2fw_t [ Types ]
[ Top ] [ m_phgamma ] [ Types ]
NAME
a2fw_t
FUNCTION
Store the Eliashberg function a2F(w).
SOURCE
300 type,public :: a2fw_t 301 302 integer :: nomega 303 ! Number of frequency points in a2f(w). 304 305 integer :: nsppol 306 ! Number of independent spin polarizations. 307 308 integer :: natom3 309 ! Number of phonon modes. 310 311 integer :: nene 312 ! Number of chemical potential values used for inelastic integration 313 314 real(dp) :: enemin 315 ! Minimal chemical potential value used for inelastic integration Copied from fstab 316 317 real(dp) :: deltaene 318 ! Chemical potential increment for inelastic integration Copied from fstab 319 ! for simplicity could be made equal to phonon frequency step 320 321 real(dp) :: omega_min, omega_max 322 ! min and Max frequency (Ha) in the linear mesh. 323 324 real(dp) :: wstep 325 ! Step of the linear mesh 326 327 real(dp) :: smear 328 ! Gaussian broadening used to approximated the Dirac distribution. 329 330 integer :: nqshift 331 ! Number of shifts in the q-mesh 332 333 integer :: ngqpt(3) 334 ! The q-mesh used for calculating vals(w). 335 336 real(dp),allocatable :: qshift(:,:) 337 ! qshift(3,nqshift) 338 ! The shifts used to generate the q-mesh. 339 340 real(dp),allocatable :: n0(:) 341 ! n0(nsppol) 342 ! Electronic DOS at the Fermi level. 343 344 real(dp),allocatable :: omega(:) 345 ! omega(nomega) 346 ! Frequency mesh in Hartree (linear). 347 348 real(dp),allocatable :: vals(:,:,:) 349 ! vals(nomega, 0:natom3, nsppol) 350 ! Eliashberg function 351 ! vals(w, 1:natom3, 1:nsppol): a2f(w) decomposed per phonon branch and spin 352 ! vals(w, 0 , 1:nsppol): a2f(w) summed over phonons modes, decomposed in spin 353 354 real(dp),allocatable :: vals_ee(:,:,:,:) 355 ! vals_ee(nene,nene,nomega,nsppol) 356 ! Eliashberg function 357 ! vals(e,e',w,0,1:nsppol): a2f(e,e',w) summed over phonons modes, decomposed in spin 358 359 real(dp),allocatable :: lambdaw(:,:,:) 360 ! lambda(nomega,0:natom3,nsppol) 361 362 contains 363 364 procedure :: free => a2fw_free 365 ! Free the memory allocated in the structure. 366 367 procedure :: write => a2fw_write 368 ! Write alpha^2F(w) to an external file in text/netcdf format 369 370 procedure :: get_moment => a2fw_get_moment 371 ! Compute moments of alpha^2F(w)/w . 372 373 end type a2fw_t 374 375 public :: a2fw_init ! Calculates the FS averaged alpha^2F(w) function.
m_phgamma/a2fw_tr_free [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
a2fw_tr_free
FUNCTION
Free the memory allocated in a2f
SIDE EFFECTS
a2f<a2fw_tr_t>=Structure storing the Eliashberg function a2F.
OUTPUT
SOURCE
2522 subroutine a2fw_tr_free(a2f_tr) 2523 2524 !Arguments ------------------------------------ 2525 class(a2fw_tr_t),intent(inout) :: a2f_tr 2526 2527 ! ********************************************************************* 2528 2529 ! integer 2530 ABI_SFREE(a2f_tr%qshift) 2531 2532 ! real 2533 ABI_SFREE(a2f_tr%n0) 2534 ABI_SFREE(a2f_tr%omega) 2535 ABI_SFREE(a2f_tr%vals_in) 2536 ABI_SFREE(a2f_tr%vals_out) 2537 ABI_SFREE(a2f_tr%vals_tr) 2538 ABI_SFREE(a2f_tr%vals_tr_gen) 2539 ABI_SFREE(a2f_tr%lambdaw_tr) 2540 2541 end subroutine a2fw_tr_free
m_phgamma/a2fw_tr_init [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
a2fw_tr_init
FUNCTION
Calculates the FS averaged alpha^2F_tr,in,out(w) functions
INPUTS
cryst<crystal_t>=Info on the unit cell. ifc<ifc_type>=Interatomic force constants. gams<phgamma_t>=Structure storing the phonon linewidths. wstep=Step for linear frequency mesh in Ha. wminmax(2)=Minimum and maximum phonon frequency. Used to construct the linear mesh for A2F(w). ph_intmeth=Integration method for phonons: 1 for gaussian, 2 for tetrahedra smear=Gaussian broadening used to approximate the Dirac delta. ngqpt(3)=Divisions of the Q-mesh used for interpolating the phonon linewidths (see also nqshift and qshift). nqshift=Number of shifts used to generated the Q-mesh. qshift(3,nqshift)=The shifts. comm=MPI communicator [qintp]=If set to False, ngqgpt, nqshift, qshift and qptop are ignored and A2F(w) is computed from the IBZ values stored in gams. Default: True i.e use Fourier interpolation. [qptopt]=Controls the generation of the q-points. If not specified, the routine takes fully into account the symmetries of the system to generate the q points in the IBZone i.e. qptopt=1 Other values of qptopt can be used for debugging purpose.
OUTPUT
a2f_tr<a2fw_tr_t>=Structure storing the Eliashberg transport function a2F_tr(w).
SOURCE
2576 subroutine a2fw_tr_init(a2f_tr, gams, cryst, ifc, ph_intmeth, wstep, wminmax, smear, ngqpt, nqshift, qshift, comm,& 2577 qintp, qptopt) ! optional 2578 2579 !Arguments ------------------------------------ 2580 !scalars 2581 integer,intent(in) :: ph_intmeth,nqshift,comm 2582 integer,intent(in),optional :: qptopt 2583 real(dp),intent(in) :: wstep,smear 2584 logical,optional,intent(in) :: qintp 2585 type(phgamma_t),intent(inout) :: gams 2586 type(ifc_type),intent(in) :: ifc 2587 type(crystal_t),intent(in) :: cryst 2588 type(a2fw_tr_t),target,intent(out) :: a2f_tr 2589 !arrays 2590 integer,intent(in) :: ngqpt(3) 2591 real(dp),intent(in) :: wminmax(2),qshift(3,nqshift) 2592 2593 !Local variables ------------------------- 2594 !scalars 2595 integer,parameter :: master = 0 2596 integer :: my_qptopt,iq_ibz,nqibz,ount,my_rank,nproc,cnt 2597 integer :: mu,iw,natom3,nsppol,spin,ierr,nomega,nqbz, idir, jdir 2598 real(dp) :: omega,xx,omega_min,omega_max,ww, cpu, wall, gflops 2599 logical :: do_qintp 2600 character(len=500) :: msg 2601 type(htetra_t) :: qtetra 2602 !arrays 2603 integer :: qptrlatt(3,3),new_qptrlatt(3,3) 2604 real(dp),allocatable :: my_qshift(:,:) 2605 real(dp) :: lambda_iso(3,3), omega_log(3,3) 2606 real(dp) :: phfrq(gams%natom3) 2607 real(dp) :: gamma_in_ph(3,3,gams%natom3), gamma_out_ph(3,3,gams%natom3) 2608 real(dp) :: lambda_in_ph(3,3,gams%natom3), lambda_out_ph(3,3,gams%natom3) 2609 real(dp),allocatable :: tmp_a2f_in(:,:,:), tmp_a2f_out(:,:,:) 2610 real(dp), ABI_CONTIGUOUS pointer :: a2f_tr_1d(:) 2611 real(dp),allocatable :: qibz(:,:),wtq(:),qbz(:,:) 2612 real(dp),allocatable :: a2f_tr_1mom(:),a2f_tr_logmom(:),a2f_tr_logmom_int(:),wdt(:,:) 2613 real(dp),allocatable :: lambda_in_tetra(:,:,:,:,:),phfreq_tetra(:,:,:) 2614 real(dp),allocatable :: lambda_out_tetra(:,:,:,:,:) 2615 2616 ! ********************************************************************* 2617 2618 my_qptopt = 1; if (present(qptopt)) my_qptopt = qptopt 2619 do_qintp = .True.; if (present(qintp)) do_qintp = qintp 2620 2621 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 2622 nsppol = gams%nsppol; natom3 = gams%natom3 2623 2624 call cwtime(cpu, wall, gflops, "start") 2625 2626 if (do_qintp) then 2627 ! Generate the q-mesh by finding the IBZ and the corresponding weights. 2628 qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3) 2629 2630 call kpts_ibz_from_kptrlatt(cryst, qptrlatt, my_qptopt, nqshift, qshift, nqibz, qibz, wtq, nqbz, qbz, & 2631 new_kptrlatt=new_qptrlatt, new_shiftk=my_qshift) 2632 ABI_FREE(qbz) 2633 2634 ! Store quantities that cannot be easily (and safely) calculated if we only know the IBZ. 2635 a2f_tr%ngqpt = ngqpt; a2f_tr%nqshift = size(my_qshift, dim=2) 2636 ABI_MALLOC(a2f_tr%qshift, (3, a2f_tr%nqshift)) 2637 a2f_tr%qshift = my_qshift 2638 ABI_FREE(my_qshift) 2639 2640 else 2641 ! No interpolation. Use q-mesh parameters from gams% 2642 a2f_tr%ngqpt = gams%ngqpt; a2f_tr%nqshift = 1 2643 nqibz = gams%nqibz 2644 ABI_MALLOC(qibz, (3,nqibz)) 2645 ABI_MALLOC(wtq, (nqibz)) 2646 ABI_MALLOC(a2f_tr%qshift, (3, a2f_tr%nqshift)) 2647 qibz = gams%qibz; wtq = gams%wtq 2648 a2f_tr%qshift = zero ! Note: assuming q-mesh centered on Gamma. 2649 end if 2650 2651 call cwtime_report(" a2fw_tr_init, q-setup", cpu, wall, gflops) 2652 2653 ! Define Min and max frequency for the mesh (enlarge it a bit) 2654 omega_min = wminmax(1); omega_max = wminmax(2) 2655 omega_min = omega_min - 0.1*abs(omega_min) 2656 if (omega_min >= zero) omega_min = one/Ha_meV 2657 omega_max = omega_max + 0.1*abs(omega_max) 2658 2659 a2f_tr%nsppol = nsppol; a2f_tr%natom3 = gams%natom3; a2f_tr%smear = smear 2660 a2f_tr%omega_min = omega_min; a2f_tr%omega_max = omega_max 2661 nomega = int((omega_max - omega_min) / wstep); a2f_tr%nomega = nomega; a2f_tr%wstep = wstep 2662 2663 ABI_MALLOC(a2f_tr%n0, (nsppol)) 2664 a2f_tr%n0 = gams%n0 2665 ! Build linear mesh. 2666 ABI_MALLOC(a2f_tr%omega, (nomega)) 2667 a2f_tr%omega = arth(omega_min, wstep, nomega) 2668 2669 ABI_CALLOC(a2f_tr%vals_in, (nomega, 3, 3, 0:natom3, nsppol)) 2670 ABI_CALLOC(a2f_tr%vals_out, (nomega,3, 3, 0:natom3, nsppol)) 2671 ABI_CALLOC(a2f_tr%vals_tr, (nomega,3, 3, 0:natom3, nsppol)) 2672 ABI_CALLOC(a2f_tr%lambdaw_tr, (nomega, 3, 3, 0:natom3, nsppol)) 2673 ABI_MALLOC(tmp_a2f_in, (nomega, 3, 3)) 2674 ABI_MALLOC(tmp_a2f_out, (nomega, 3, 3)) 2675 2676 if (ph_intmeth == 2) then 2677 call cwtime(cpu, wall, gflops, "start") 2678 2679 ! Prepare tetrahedron integration. 2680 qptrlatt = 0; qptrlatt(1, 1) = a2f_tr%ngqpt(1); qptrlatt(2, 2) = a2f_tr%ngqpt(2); qptrlatt(3, 3) = a2f_tr%ngqpt(3) 2681 qtetra = tetra_from_kptrlatt(cryst, my_qptopt, qptrlatt, a2f_tr%nqshift, a2f_tr%qshift, nqibz, qibz, comm, msg, ierr) 2682 if (ierr/=0) ABI_ERROR(msg) 2683 2684 ABI_MALLOC_OR_DIE(lambda_in_tetra, (nqibz, 3, 3, natom3, nsppol), ierr) 2685 ABI_MALLOC_OR_DIE(lambda_out_tetra, (nqibz, 3, 3, natom3, nsppol), ierr) 2686 lambda_in_tetra = zero 2687 lambda_out_tetra = zero 2688 2689 ABI_MALLOC_OR_DIE(phfreq_tetra, (nqibz, natom3, nsppol), ierr) 2690 phfreq_tetra = zero 2691 2692 call cwtime_report(" a2fw_tr_init%tetra,", cpu, wall, gflops) 2693 end if 2694 2695 call cwtime(cpu, wall, gflops, "start") 2696 2697 ! Loop over spins and qpoints in the IBZ 2698 cnt = 0 2699 do spin=1,nsppol 2700 do iq_ibz=1,nqibz 2701 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle 2702 2703 ! Interpolate or evaluate lambda transport coefficients. 2704 if (do_qintp) then 2705 call phgamma_vv_interp(gams,cryst,ifc,spin,qibz(:,iq_ibz),phfrq,gamma_in_ph,gamma_out_ph,lambda_in_ph,lambda_out_ph) 2706 else 2707 call phgamma_vv_eval_qibz(gams,cryst,ifc,iq_ibz,spin,phfrq,gamma_in_ph,gamma_out_ph,lambda_in_ph,lambda_out_ph) 2708 end if 2709 2710 select case (ph_intmeth) 2711 case (1) 2712 ! Gaussian: Add all contributions from the phonon modes at this qpoint to a2f_tr 2713 ! (note that unstable modes are included). 2714 do mu=1,natom3 2715 tmp_a2f_in = zero 2716 tmp_a2f_out = zero 2717 do iw=1,nomega 2718 xx = a2f_tr%omega(iw) - phfrq(mu) 2719 tmp_a2f_in(iw,:,:) = tmp_a2f_in(iw,:,:) + gaussian(xx, smear) * lambda_in_ph(:,:,mu) * abs(phfrq(mu)) 2720 tmp_a2f_out(iw,:,:) = tmp_a2f_out(iw,:,:) + gaussian(xx, smear) * lambda_out_ph(:,:,mu) * abs(phfrq(mu)) 2721 end do 2722 a2f_tr%vals_in(:,:,:,mu,spin) = a2f_tr%vals_in(:,:,:,mu,spin) + tmp_a2f_in(:,:,:) * wtq(iq_ibz) 2723 a2f_tr%vals_out(:,:,:,mu,spin) = a2f_tr%vals_out(:,:,:,mu,spin) + tmp_a2f_out(:,:,:) * wtq(iq_ibz) 2724 end do 2725 2726 case (2) 2727 ! Tetra: store data. 2728 do mu=1,natom3 2729 lambda_in_tetra(iq_ibz, :,:, mu, spin) = lambda_in_ph(:,:,mu) * abs(phfrq(mu)) 2730 lambda_out_tetra(iq_ibz, :,:, mu, spin) = lambda_out_ph(:,:,mu) * abs(phfrq(mu)) 2731 phfreq_tetra(iq_ibz, mu, spin) = phfrq(mu) 2732 end do 2733 end select 2734 2735 end do ! iq_ibz 2736 end do ! spin 2737 2738 if (ph_intmeth == 2) then 2739 ! Collect results on each node. 2740 call xmpi_sum(lambda_in_tetra, comm, ierr) 2741 call xmpi_sum(lambda_out_tetra, comm, ierr) 2742 call xmpi_sum(phfreq_tetra, comm, ierr) 2743 2744 ! workspace for tetra. 2745 ABI_MALLOC(wdt, (nomega, 2)) 2746 2747 ! TODO: with the tetra_get_onewk call we can integrate this above 2748 ! and avoid allocating all of lambda_in_tetra and phfreq_tetra!!! 2749 ! For each mode get its contribution 2750 cnt = 0 2751 do spin=1,nsppol 2752 do mu=1,natom3 2753 do iq_ibz=1,nqibz 2754 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! mpi-parallelism 2755 2756 call qtetra%get_onewk(iq_ibz, gams%bcorr, nomega, nqibz, phfreq_tetra(:,mu,spin), omega_min, omega_max, one, wdt) 2757 wdt = wdt * wtq(iq_ibz) 2758 2759 ! Accumulate (Integral of a2F_tr is computed afterwards) 2760 do idir=1,3 2761 do jdir=1,3 2762 a2f_tr%vals_in(:,idir,jdir,mu,spin) = a2f_tr%vals_in(:,idir,jdir,mu,spin) & 2763 + wdt(:,1) * lambda_in_tetra(iq_ibz,idir,jdir,mu,spin) 2764 a2f_tr%vals_out(:,idir,jdir,mu,spin) = a2f_tr%vals_out(:,idir,jdir,mu,spin) & 2765 + wdt(:,1) * lambda_out_tetra(iq_ibz,idir,jdir,mu,spin) 2766 !a2f_tr%lambdaw(:,mu,spin) = a2f_tr%lambdaw(:,mu,spin) + wdt(:,2) * lambda_XX_tetra(iq_ibz, mu, spin) 2767 end do 2768 end do 2769 end do 2770 end do 2771 end do 2772 2773 ! Free memory allocated for qtetra. 2774 ABI_FREE(wdt) 2775 ABI_FREE(lambda_in_tetra) 2776 ABI_FREE(lambda_out_tetra) 2777 ABI_FREE(phfreq_tetra) 2778 call qtetra%free() 2779 end if 2780 2781 ! Collect final results on each node 2782 call xmpi_sum(a2f_tr%vals_in, comm, ierr) 2783 call xmpi_sum(a2f_tr%vals_out, comm, ierr) 2784 2785 ! TODO: divide by g(eF, spin) 2786 ! check normalization with N(0) and <v^2> from Savrasov eq 20 2787 do spin=1,nsppol 2788 a2f_tr%vals_in(:,:,:,0,spin) = sum(a2f_tr%vals_in(:,:,:,1:natom3,spin), dim=4) 2789 a2f_tr%vals_out(:,:,:,0,spin) = sum(a2f_tr%vals_out(:,:,:,1:natom3,spin), dim=4) 2790 !a2f_tr%vals(:,:,spin) = a2f_tr%vals(:,:,spin) / (two_pi * a2f_tr%n0(spin)) 2791 end do 2792 2793 a2f_tr%vals_tr = a2f_tr%vals_out - a2f_tr%vals_in 2794 2795 call cwtime_report(" a2fw_tr_init, a2f_eval", cpu, wall, gflops) 2796 2797 ! NOTE: for the moment calculate the log moms etc on just the trace of the a2f functions 2798 ABI_MALLOC(a2f_tr_1mom, (nomega)) 2799 ABI_MALLOC(a2f_tr_logmom, (nomega)) 2800 ABI_MALLOC(a2f_tr_logmom_int, (nomega)) 2801 2802 !call a2fw_tr_print_info() 2803 2804 ! Compute lambda_tr(w) (resolved in spin mode and directions) 2805 do spin=1,nsppol 2806 do mu=0,natom3 2807 do idir=1,3 2808 do jdir=1,3 2809 2810 do iw=1,nomega 2811 ww = a2f_tr%omega(iw) 2812 if (abs(ww) > EPHTK_WTOL) then 2813 a2f_tr_1mom(iw) = a2f_tr%vals_tr(iw,idir,jdir,mu,spin) / abs(ww) 2814 else 2815 a2f_tr_1mom(iw) = zero 2816 end if 2817 ! TODO: Strange that the first value of int(f) is not set to zero! 2818 ! FIXME: The frequency integration overestimates lambda if there are negative phonon frequencies! 2819 end do 2820 2821 call simpson_int(nomega, wstep, a2f_tr_1mom, a2f_tr%lambdaw_tr(:,idir,jdir,mu,spin)) 2822 end do ! jdir 2823 end do ! idir 2824 end do ! mu 2825 end do ! spin 2826 2827 if (my_rank == master) then 2828 ount = std_out 2829 if (nsppol > 1) then 2830 write(msg,'(3a)') ch10,'Comment: some of the following quantities should be integrated over spin', ch10 2831 call wrtout(ount, msg) 2832 end if 2833 end if 2834 2835 do spin=1,nsppol 2836 lambda_iso = a2fw_tr_moment(a2f_tr, 0, spin) 2837 2838 do idir=1,3 2839 do jdir=1,3 2840 a2f_tr_1d => a2f_tr%vals_tr(:,idir,jdir,0,spin) 2841 2842 ! Get log moment of alpha^2F. 2843 a2f_tr_logmom = zero 2844 do iw=1,nomega 2845 omega = a2f_tr%omega(iw) 2846 if (abs(omega) > EPHTK_WTOL .and. abs(lambda_iso(idir,jdir)) > EPHTK_WTOL) then 2847 a2f_tr_logmom(iw) = (two/lambda_iso(idir,jdir)) * a2f_tr_1d(iw)*log(abs(omega))/abs(omega) 2848 end if 2849 end do 2850 call simpson_int(nomega, wstep, a2f_tr_logmom, a2f_tr_logmom_int) 2851 !write(std_out,*)' iw,nomega,greatest_real,a2f_tr_logmom_int(nomega)=',& iw,nomega,greatest_real,a2f_tr_logmom_int(nomega) 2852 if(abs(a2f_tr_logmom_int(nomega)) < log(greatest_real*tol6)) then 2853 omega_log(idir,jdir) = exp(a2f_tr_logmom_int(nomega)) 2854 else 2855 omega_log(idir,jdir)=greatest_real*tol6 2856 endif 2857 end do 2858 end do 2859 2860 ! TODO: make output only for irred values xx yy zz and top half of matrix 2861 if (my_rank == master) then 2862 if (do_qintp) then 2863 write(ount,'(a)')' Evaluation of parameters analogous to electron-phonon coupling for 3x3 directions (interpolated) ' 2864 else 2865 write(ount,'(a)')' Evaluation of parameters analogous to electron-phonon coupling for 3x3 directions (coarse grid) ' 2866 endif 2867 write(ount,'(a,3(3es10.3,2x))') ' lambda = ',lambda_iso 2868 write(ount,'(a,3(3es10.3,2x),a)' )' omegalog = ',omega_log,' (Ha) ' 2869 write(ount,'(a,3(3es10.3,2x),a)' )' ',omega_log*Ha_K, ' (Kelvin) ' 2870 write(ount,"(a)")' positive moments of alpha2Ftr:' 2871 write(ount,'(a,9(es10.3,2x))' )' lambda <omega^2> = ',a2fw_tr_moment(a2f_tr, 2, spin) 2872 write(ount,'(a,9(es10.3,2x))' )' lambda <omega^3> = ',a2fw_tr_moment(a2f_tr, 3, spin) 2873 write(ount,'(a,9(es10.3,2x))' )' lambda <omega^4> = ',a2fw_tr_moment(a2f_tr, 4, spin) 2874 write(ount,'(a,9(es10.3,2x))' )' lambda <omega^5> = ',a2fw_tr_moment(a2f_tr, 5, spin) 2875 end if 2876 end do 2877 2878 ABI_FREE(tmp_a2f_in) 2879 ABI_FREE(tmp_a2f_out) 2880 ABI_FREE(a2f_tr_1mom) 2881 ABI_FREE(a2f_tr_logmom) 2882 ABI_FREE(a2f_tr_logmom_int) 2883 ABI_FREE(qibz) 2884 ABI_FREE(wtq) 2885 2886 end subroutine a2fw_tr_init
m_phgamma/a2fw_tr_moment [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
a2fw_tr_moment
FUNCTION
Compute \int dw [a2F_tr(w)/w] w^n From Allen PRL 59 1460 [[cite:Allen1987]] and later PRB papers. See also [[cite:Grimvall1981]] book.
INPUTS
a2f_tr<a2fw_tr_t>=Structure storing the Eliashberg function. nn=Value of n spin=The spin component
OUTPUT
a2fw_tr_moment = \int dw [a2F_tr(w)/w] w^n
[out_int(x)] = \int^{x} dw [a2F_tr(w)/w] w^n
SOURCE
2217 function a2fw_tr_moment(a2f_tr, nn, spin, out_int) 2218 2219 !Arguments ------------------------------------ 2220 !scalars 2221 integer,intent(in) :: spin,nn 2222 type(a2fw_tr_t),intent(in) :: a2f_tr 2223 !arrays 2224 real(dp),intent(out),optional :: out_int(a2f_tr%nomega,3,3) 2225 real(dp) :: a2fw_tr_moment(3,3) 2226 2227 !Local variables ------------------------- 2228 !scalars 2229 integer :: iw, idir, jdir 2230 real(dp) :: omg,omg_nm1 2231 !arrays 2232 real(dp) :: ff(a2f_tr%nomega),int_ff(a2f_tr%nomega) 2233 2234 ! ********************************************************************* 2235 2236 do jdir = 1, 3 2237 do idir = 1, 3 2238 ! Construct the integrand function. [a2F_tr(w)/w] w^n 2239 ff = zero; int_ff = zero 2240 2241 if (nn-1 >= 0) then 2242 do iw=1,a2f_tr%nomega 2243 omg = a2f_tr%omega(iw) 2244 omg_nm1 = omg ** (nn-1) 2245 ff(iw) = a2f_tr%vals_tr(iw,idir,jdir,0,spin) * omg_nm1 2246 end do 2247 else 2248 do iw=1,a2f_tr%nomega 2249 omg = a2f_tr%omega(iw) 2250 omg_nm1 = zero; if (abs(omg) > EPHTK_WTOL) omg_nm1 = omg**(nn-1) 2251 ff(iw) = a2f_tr%vals_tr(iw,idir,jdir,0,spin) * omg_nm1 2252 end do 2253 end if 2254 2255 ! Integration with simpson rule on a linear mesh. 2256 call simpson_int(a2f_tr%nomega, a2f_tr%wstep, ff, int_ff) 2257 2258 a2fw_tr_moment(idir, jdir) = int_ff(a2f_tr%nomega) 2259 if (present(out_int)) out_int(:,idir,jdir) = int_ff(:) 2260 end do 2261 end do 2262 2263 end function a2fw_tr_moment
m_phgamma/a2fw_tr_t [ Types ]
[ Top ] [ m_phgamma ] [ Types ]
NAME
a2fw_tr_t
FUNCTION
Store the Eliashberg transport spectral functions:
a2F_trin(w, x, x')
a2F_trout(w, x, x')
a2F_tr(w, x, x') = in - out
a2F_tr_gen(e, e', w, x, x')
SOURCE
391 type,public :: a2fw_tr_t 392 393 integer :: nomega 394 ! Number of frequency points in a2f_tr(w). 395 396 integer :: nene 397 ! Number of chemical potential values used for inelastic integration 398 ! Number of electron points in a2f_tr_gen(e,e',w). 399 400 integer :: nsppol 401 ! Number of independent spin polarizations. 402 403 integer :: natom3 404 ! Number of phonon modes. 405 406 real(dp) :: enemin 407 ! Minimal chemical potential value used for inelastic integration 408 409 real(dp) :: deltaene 410 ! Chemical potential increment for inelastic integration 411 412 real(dp) :: omega_min,omega_max 413 ! min and Max frequency (Ha) in the linear mesh. 414 415 real(dp) :: wstep 416 ! Step of the linear mesh 417 418 real(dp) :: smear 419 ! Gaussian broadening used to approximated the Dirac distribution. 420 421 integer :: nqshift 422 ! Number of shifts in the q-mesh 423 424 integer :: ngqpt(3) 425 ! The q-mesh used for calculating vals(w). 426 427 real(dp),allocatable :: qshift(:,:) 428 ! qshift(3,nqshift) 429 ! The shifts used to generate the q-mesh. 430 431 real(dp),allocatable :: n0(:) 432 ! n0(nsppol) 433 ! Electronic DOS at the Fermi level. 434 435 real(dp),allocatable :: omega(:) 436 ! omega(nomega) 437 ! Frequency mesh in Hartree (linear). 438 439 real(dp),allocatable :: vals_in(:,:,:,:,:) 440 real(dp),allocatable :: vals_out(:,:,:,:,:) 441 ! vals_in(nomega,3,3,0:natom3,nsppol) 442 ! Eliashberg transport functions for in and out scattering 443 ! vals_in(w,3,3,1:natom3,1:nsppol): a2f_tr(w) decomposed per phonon branch and spin 444 ! vals_in(w,3,3,0,1:nsppol): a2f_tr(w) summed over phonons modes, decomposed in spin 445 446 real(dp),allocatable :: vals_tr(:,:,:,:,:) 447 ! vals_tr(nomega,3,3,0:natom3,nsppol) 448 ! transport spectral function = in-out 449 450 real(dp),allocatable :: vals_tr_gen(:,:,:,:,:,:,:) 451 ! vals(nene,nene,nomega,3,3,nsppol) 452 ! generalized transport spectral function from PB Allen Phys. Rev. Lett. 59, 1460 (1987) [[cite:Allen1987]] 453 454 real(dp),allocatable :: lambdaw_tr(:,:,:,:,:) 455 ! lambda(nomega,3,3,0:natom3,nsppol) 456 457 contains 458 459 procedure :: free => a2fw_tr_free 460 ! Free the memory allocated in the structure. 461 462 procedure :: write => a2fw_tr_write 463 ! Write alpha^2F(w) to an external file in text/netcdf format 464 465 end type a2fw_tr_t 466 467 public :: a2fw_tr_init ! Calculates the FS averaged alpha^2F_tr,in,out(w, x, x') functions.
m_phgamma/a2fw_tr_write [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
a2fw_tr_write
FUNCTION
Write alpha^2F_tr(w) to an external file in text form
INPUTS
a2f_tr<a2fw_tr_t>=Container storing the Eliashberg transport functions. basename=Filename for output. post=String appended to netcdf variables e.g. _qcoarse, _qintp ncid=Netcdf file handler. Set it to nctk_noid to disable output.
OUTPUT
Output is written to file. This routine should be called by one MPI proc.
SOURCE
2909 subroutine a2fw_tr_write(a2f_tr, basename, post, ncid) 2910 2911 !Arguments ------------------------------------ 2912 !scalars 2913 integer,intent(in) :: ncid 2914 character(len=*),intent(in) :: basename, post 2915 class(a2fw_tr_t),intent(in) :: a2f_tr 2916 2917 !Local variables ------------------------- 2918 !scalars 2919 integer :: iw,spin,unt,ii,mu,idir, jdir 2920 integer :: ncerr 2921 character(len=500) :: dim1_name 2922 character(len=500) :: msg 2923 character(len=fnlen) :: path 2924 2925 ! ********************************************************************* 2926 2927 ! Write spin-resolved a2F_tr(w) 2928 path = strcat(basename, "_A2FW_tr") 2929 if (open_file(path, msg, newunit=unt, form="formatted", action="write", status="unknown") /= 0) then 2930 ABI_ERROR(msg) 2931 end if 2932 2933 call write_a2fw_tr_header() 2934 2935 ! Reminder: vals_tr(nomega,3,3,0:natom3,nsppol) 2936 2937 if (a2f_tr%nsppol == 1) then 2938 write(unt,'(a)')"# Frequency, a2F_tr(w,i,j)" 2939 do iw=1,a2f_tr%nomega 2940 write(unt,'(e16.6, 2x, 3(3e16.6,1x))') a2f_tr%omega(iw), a2f_tr%vals_tr(iw,:,:,0,1) 2941 end do 2942 write(unt,'(3a)')ch10, ch10, "# Frequency, lambda(w,i,j)" 2943 do iw=1,a2f_tr%nomega 2944 write(unt,'(e16.6, 2x, 3(3e16.6,1x))') a2f_tr%omega(iw), a2f_tr%lambdaw_tr(iw,:,:,0,1) 2945 end do 2946 2947 else 2948 write(unt,'(a)')"# Frequency, a2F_tr_tot(w,i,j), a2F_tr_spin1(w,i,j) ..." 2949 do iw=1,a2f_tr%nomega 2950 write(unt,'(e16.6, 2x, 3(3e16.6,1x), 2x, 3(3e16.6,1x), 2x, 3(3e16.6,1x))') a2f_tr%omega(iw), & 2951 ((sum(a2f_tr%vals_tr(iw,idir,jdir,0,:)), idir=1,3), jdir=1,3) , & 2952 a2f_tr%vals_tr(iw,:,:,0,1), & ! UP 2953 a2f_tr%vals_tr(iw,:,:,0,2) ! DOWN 2954 end do 2955 write(unt,'(3a)') ch10, ch10, "# Frequency, lambda_tr_tot(w,i,j)dw, lambda_tr_spin1(w,i,j) ..." 2956 do iw=1,a2f_tr%nomega 2957 write(unt,'(e16.6, 2x, 3(3e16.6,1x), 2x, 3(3e16.6,1x), 2x, 3(3e16.6,1x), 2x, 3(3e16.6,1x))') a2f_tr%omega(iw), & 2958 ((sum(a2f_tr%lambdaw_tr(iw,idir,jdir,0,:)), idir=1,3), jdir=1,3) , & ! TOT 2959 a2f_tr%lambdaw_tr(iw,:,:,0,1), & ! UP 2960 a2f_tr%lambdaw_tr(iw,:,:,0,2) ! DOWN 2961 end do 2962 end if 2963 2964 close(unt) 2965 2966 ! Write phonon mode contributions to a2F_tr(w,i,j) 2967 path = strcat(basename, "_PH_A2FW_tr") 2968 if (open_file(path, msg, newunit=unt, form="formatted", action="write", status="unknown") /= 0) then 2969 ABI_ERROR(msg) 2970 end if 2971 call write_a2fw_tr_header() 2972 2973 ! Reminder: vals_tr(nomega,3,3,0:natom3,nsppol) 2974 if (a2f_tr%nsppol == 1) then 2975 do mu=0,a2f_tr%natom3 2976 write(unt,'(a,i0)')"# Phonon mode ",mu 2977 write(unt,'(a)')"# Frequency, a2F_tr(w)" 2978 do iw=1,a2f_tr%nomega 2979 write(unt,'(e16.6,2x,3(3e16.6,1x))') a2f_tr%omega(iw), a2f_tr%vals_tr(iw,:,:,mu,1) 2980 end do 2981 do ii=1,2; write(unt,'(a)')""; end do 2982 write(unt,'(a)')"# Frequency, lambda_tr(w)" 2983 do iw=1,a2f_tr%nomega 2984 write(unt,'(e16.6,2x,3(3e16.6,1x))') a2f_tr%omega(iw), a2f_tr%lambdaw_tr(iw,:,:,mu,1) 2985 end do 2986 do ii=1,2; write(unt,'(a)')""; end do 2987 end do 2988 2989 else 2990 do mu=0,a2f_tr%natom3 2991 write(unt,'(a,i0)')"# Phonon mode ",mu 2992 write(unt,'(a)')"# Frequency, a2F_tot(w), a2F_spin1(w) ..." 2993 do iw=1,a2f_tr%nomega 2994 write(unt,*) a2f_tr%omega(iw), & 2995 sum(a2f_tr%vals_tr(iw,:,:,mu,:)), & ! TOT 2996 a2f_tr%vals_tr(iw,:,:,mu,1) , & ! UP 2997 a2f_tr%vals_tr(iw,:,:,mu,2) ! DOWN 2998 end do 2999 write(unt,'(3a)')ch10,ch10,"# Frequency, lambda_tot(w)dw, lambda_spin1(w) ..." 3000 do iw=1,a2f_tr%nomega 3001 write(unt,*) a2f_tr%omega(iw), & 3002 sum(a2f_tr%lambdaw_tr(iw,:,:,mu,:)), & ! TOT 3003 a2f_tr%lambdaw_tr(iw,:,:,mu,1), & ! UP 3004 a2f_tr%lambdaw_tr(iw,:,:,mu,2) ! DOWN 3005 end do 3006 end do 3007 end if 3008 3009 close(unt) 3010 3011 if (ncid /= nctk_noid) then 3012 ! Define dimensions. 3013 dim1_name = strcat("a2ftr_nomega", post) 3014 ncerr = nctk_def_dims(ncid, [ & 3015 nctkdim_t(dim1_name, a2f_tr%nomega), & 3016 nctkdim_t("natom3p1", a2f_tr%natom3 + 1), nctkdim_t("number_of_spins", a2f_tr%nsppol) & 3017 ], defmode=.True.) 3018 NCF_CHECK(ncerr) 3019 3020 ncerr = nctk_def_arrays(ncid, [ & 3021 nctkarr_t(strcat('a2ftr_mesh', post), "dp", dim1_name), & 3022 nctkarr_t(strcat('a2ftr_values', post), "dp", strcat(dim1_name, ", three, three, natom3p1, number_of_spins")), & 3023 nctkarr_t(strcat('a2ftr_lambdaw', post), "dp", strcat(dim1_name, ", three, three, natom3p1, number_of_spins")) & 3024 ]) 3025 NCF_CHECK(ncerr) 3026 3027 NCF_CHECK(nctk_set_datamode(ncid)) 3028 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2ftr_mesh", post)), a2f_tr%omega)) 3029 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2ftr_values", post)), a2f_tr%vals_tr)) 3030 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2ftr_lambdaw", post)), a2f_tr%lambdaw_tr)) 3031 end if 3032 3033 contains 3034 3035 subroutine write_a2fw_tr_header() 3036 3037 ! Output the header. 3038 write(unt,'(a)') '#' 3039 write(unt,'(a)') '# ABINIT package: a2F_tr(w) file' 3040 write(unt,'(a)') '#' 3041 write(unt,'(a)') '# a2F_tr(w) function integrated over the FS. omega in a.u.' 3042 write(unt,'(2a)') '# ngqpt: ',trim(ltoa(a2f_tr%ngqpt)) 3043 write(unt,'(a,i0,3(a,e16.6))')'# number of frequencies: ',a2f_tr%nomega," between omega_min: ",a2f_tr%omega_min,& 3044 ' Ha and omega_max: ',a2f_tr%omega_max,' Ha with step:',a2f_tr%wstep 3045 write(unt,'(a,e16.6)') '# the smearing width for gaussians is ',a2f_tr%smear 3046 if (a2f_tr%nsppol == 2) then 3047 write(unt,'(a,e16.6)')"# Total DOS at Fermi level ",sum(a2f_tr%n0) 3048 end if 3049 do spin=1,a2f_tr%nsppol 3050 write(unt,"(a,i0,a,e16.6)")"# The DOS at Fermi level for spin ",spin," is ",a2f_tr%n0(spin) 3051 end do 3052 do ii=1,2; write(unt,'(a)') "# "; end do 3053 3054 end subroutine write_a2fw_tr_header 3055 3056 end subroutine a2fw_tr_write
m_phgamma/a2fw_write [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
a2fw_write
FUNCTION
Write alpha^2F(w) to an external file in text form
INPUTS
a2f<a2fw_t>=Container storing the Eliashberg functions. basename=Filename for output. post=String appended to netcdf variables e.g. _qcoarse, _qintp ncid=Netcdf file handler. Set it to nctk_noid to disable output.
OUTPUT
Output is written to file. This routine should be called by one MPI proc.
SOURCE
2286 subroutine a2fw_write(a2f, basename, post, ncid) 2287 2288 !Arguments ------------------------------------ 2289 !scalars 2290 integer,intent(in) :: ncid 2291 character(len=*),intent(in) :: basename, post 2292 class(a2fw_t),intent(in) :: a2f 2293 2294 !Local variables ------------------------- 2295 !scalars 2296 integer :: iw,spin,unt,ii,mu 2297 integer :: ncerr 2298 character(len=500) :: dim1_name 2299 character(len=500) :: msg 2300 character(len=fnlen) :: path 2301 2302 ! ********************************************************************* 2303 2304 ! Write spin-resolved a2F(w) 2305 path = strcat(basename, "_A2FW") 2306 if (open_file(path, msg, newunit=unt, form="formatted", action="write", status="unknown") /= 0) then 2307 ABI_ERROR(msg) 2308 end if 2309 2310 call write_a2fw_header() 2311 2312 if (a2f%nsppol == 1) then 2313 write(unt,'(a)')"# Frequency, a2F(w), lambda(w)" 2314 do iw=1,a2f%nomega 2315 write(unt,*) a2f%omega(iw), a2f%vals(iw,0,1), a2f%lambdaw(iw,0,1) 2316 end do 2317 2318 else 2319 write(unt,'(a)')"# Frequency, a2F_tot(w), lambda_tot(w)dw, a2F_spin1(w), lambda_spin1(w) ..." 2320 do iw=1,a2f%nomega 2321 write(unt,'(E20.10,3x,3(2E20.10, 2x))') a2f%omega(iw), & 2322 sum(a2f%vals(iw,0,:)) , sum(a2f%lambdaw(iw,0,:)), & ! TOT 2323 a2f%vals(iw,0,1) , a2f%lambdaw(iw,0,1), & ! UP 2324 a2f%vals(iw,0,2) , a2f%lambdaw(iw,0,2) ! DOWN 2325 end do 2326 end if 2327 2328 close(unt) 2329 2330 ! Write phonon contributions to a2F(w) 2331 path = strcat(basename, "_PH_A2FW") 2332 if (open_file(path, msg, newunit=unt, form="formatted", action="write", status="unknown") /= 0) then 2333 ABI_ERROR(msg) 2334 end if 2335 call write_a2fw_header() 2336 2337 if (a2f%nsppol == 1) then 2338 do mu=0,a2f%natom3 2339 write(unt,'(a,i0)')"# Phonon mode ",mu 2340 write(unt,'(a)')"# Frequency, a2F(w), lambda(w)" 2341 do iw=1,a2f%nomega 2342 write(unt,*) a2f%omega(iw), a2f%vals(iw,mu,1), a2f%lambdaw(iw,mu,1) 2343 end do 2344 do ii=1,2; write(unt,'(a)')""; end do 2345 end do 2346 2347 else 2348 do mu=0,a2f%natom3 2349 write(unt,'(a,i0)')"# Phonon mode ",mu 2350 write(unt,'(a)')"# Frequency, a2F_tot(w), lambda_tot(w)dw, a2F_spin1(w), lambda_spin1(w) ..." 2351 do iw=1,a2f%nomega 2352 write(unt,'(E20.10,3x,3(2E20.10, 2x))') a2f%omega(iw), & 2353 sum(a2f%vals(iw,mu,:)), sum(a2f%lambdaw(iw,mu,:)), & ! TOT 2354 a2f%vals(iw,mu,1) , a2f%lambdaw(iw,mu,1), & ! UP 2355 a2f%vals(iw,mu,2) , a2f%lambdaw(iw,mu,2) ! DOWN 2356 end do 2357 end do 2358 end if 2359 2360 close(unt) 2361 2362 if (ncid /= nctk_noid) then 2363 ! Define dimensions. 2364 dim1_name = strcat("a2f_nomega", post) 2365 ncerr = nctk_def_dims(ncid, [ & 2366 nctkdim_t(dim1_name, a2f%nomega), & 2367 nctkdim_t("natom3p1", a2f%natom3 + 1), nctkdim_t("number_of_spins", a2f%nsppol) & 2368 ], defmode=.True.) 2369 NCF_CHECK(ncerr) 2370 2371 ncerr = nctk_def_arrays(ncid, [ & 2372 nctkarr_t(strcat('a2f_mesh', post), "dp", dim1_name), & 2373 nctkarr_t(strcat('a2f_values', post), "dp", strcat(dim1_name, ", natom3p1, number_of_spins")), & 2374 nctkarr_t(strcat('a2f_lambdaw', post), "dp", strcat(dim1_name, ", natom3p1, number_of_spins")) & 2375 ]) 2376 NCF_CHECK(ncerr) 2377 2378 NCF_CHECK(nctk_set_datamode(ncid)) 2379 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2f_mesh", post)), a2f%omega)) 2380 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2f_values", post)), a2f%vals)) 2381 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, strcat("a2f_lambdaw", post)), a2f%lambdaw)) 2382 end if 2383 2384 contains 2385 2386 subroutine write_a2fw_header() 2387 2388 ! Output the header. 2389 write(unt,'(a)') '#' 2390 write(unt,'(a)') '# ABINIT package: a2F(w) file' 2391 write(unt,'(a)') '#' 2392 write(unt,'(a)') '# a2F(w) function integrated over the FS. omega in a.u.' 2393 write(unt,'(2a)') '# ngqpt: ',trim(ltoa(a2f%ngqpt)) 2394 write(unt,'(a,i0,3(a,e16.6))')'# number of frequencies: ',a2f%nomega," between omega_min: ",a2f%omega_min,& 2395 ' Ha and omega_max: ',a2f%omega_max,' Ha with step:',a2f%wstep 2396 write(unt,'(a,e16.6)') '# the smearing width for gaussians is ',a2f%smear 2397 write(unt,'(a,e16.6)')"# Total DOS at Fermi level ",sum(a2f%n0) 2398 do spin=1,a2f%nsppol 2399 write(unt,"(a,i0,a,e16.6)")"# The DOS at Fermi level for spin ",spin," is ",a2f%n0(spin) 2400 end do 2401 do ii=1,2; write(unt,'(a)') "# "; end do 2402 2403 end subroutine write_a2fw_header 2404 2405 end subroutine a2fw_write
m_phgamma/calc_dbldelta [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
calc_dbldelta
FUNCTION
Compute Tetrahedron weights for the double delta.
SOURCE
4721 subroutine calc_dbldelta(cryst, ebands, ltetra, bstart, bstop, nqibz, qibz, wtqs, comm) 4722 4723 !Arguments ------------------------------------ 4724 type(crystal_t),intent(in) :: cryst 4725 type(ebands_t),intent(in) :: ebands 4726 integer,intent(in) :: ltetra, bstart, bstop, nqibz, comm 4727 !arrays 4728 real(dp),intent(in) :: qibz(3,nqibz) 4729 real(dp),intent(out) :: wtqs(nqibz, ebands%nsppol) 4730 4731 !Local variables------------------------------- 4732 !scalars 4733 integer,parameter :: enough = 50 4734 integer :: nkbz, ierr, nb, ik_bz, ik_ibz, ikq_ibz, i1, i2, i3, nene, spin, iq_ibz 4735 integer :: ib1, ib2, my_rank, nproc, cnt 4736 real(dp),parameter :: max_occ1 = one 4737 real(dp) :: enemin, enemax !cpu, wall, gflops, 4738 !character(len=500) :: msg 4739 type(krank_t) :: ibz_krank 4740 type(t_tetrahedron) :: tetra 4741 !character(len=80) :: errorstring 4742 !arrays 4743 integer :: nge(3), ngw(3) 4744 integer,allocatable :: indkpt(:) !, symrecfm(:,:,:) 4745 real(dp) :: qpt(3), kk(3), kq(3) 4746 real(dp),allocatable :: eig_k(:,:), eig_kq(:,:), wght_bz(:,:,:), kbz(:,:) 4747 real(dp),allocatable :: work_k(:), work_kq(:), dtweightde(:,:,:), tweight(:,:,:) 4748 4749 ! ************************************************************************* 4750 4751 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 4752 !call cwtime(cpu, wall, gflops, "start") 4753 4754 ABI_CHECK(nqibz > 1, "Need more that 1 q-point") 4755 ABI_CHECK(isdiagmat(ebands%kptrlatt), "kptrlatt must be diagonal when tetra is used.") 4756 ABI_CHECK(ebands%nshiftk == 1, "nshiftk must be 1 when tetra is used") 4757 nge = get_diag(ebands%kptrlatt); ngw = nge 4758 nkbz = product(nge(1:3)) 4759 4760 ! TODO: Handle symmetries in a cleaner way. Change API of krank_new to pass symafm and kptopt 4761 ibz_krank = krank_new(ebands%nkpt, ebands%kptns, nsym=cryst%nsym, symrec=cryst%symrec, & 4762 time_reversal=kpts_timrev_from_kptopt(ebands%kptopt) == 1) 4763 4764 nb = bstop - bstart + 1 4765 ABI_MALLOC(kbz, (3, nkbz)) 4766 ABI_MALLOC(indkpt, (nkbz)) 4767 ABI_MALLOC(eig_k, (nb, nkbz)) 4768 ABI_MALLOC(eig_kq, (nb, nkbz)) 4769 4770 wtqs = zero 4771 cnt = 0 4772 do spin=1,ebands%nsppol 4773 do iq_ibz=1,nqibz 4774 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism. 4775 qpt = qibz(:, iq_ibz) 4776 4777 ! The double delta is ill-defined for q == 0 4778 if (all(abs(qpt) < tol12)) cycle 4779 4780 ierr = 0; ik_bz = 0 4781 do i3=0,nge(3) - 1 4782 do i2=0,nge(2) - 1 4783 do i1=0,nge(1) - 1 4784 ik_bz = ik_bz + 1 4785 indkpt(ik_bz) = ik_bz 4786 4787 ! Find correspondence between the grid and the IBZ 4788 kk = ([i1, i2, i3] + ebands%shiftk(:, 1)) / nge(:) 4789 kbz(:, ik_bz) = kk 4790 4791 ik_ibz = ibz_krank%get_index(kk) 4792 if (ik_ibz < 1) then 4793 if (ierr <= enough) then 4794 ABI_WARNING(sjoin('kpt:', trim(ktoa(kk)), 'has no symmetric among the k-points!')) 4795 end if 4796 ierr = ierr + 1; cycle 4797 end if 4798 4799 eig_k(:, ik_bz) = ebands%eig(bstart:bstop, ik_ibz, spin) 4800 4801 ! Find correspondence between the BZ grid and the IBZ. 4802 kq = kk + qpt 4803 ikq_ibz = ibz_krank%get_index(kq) 4804 4805 if (ikq_ibz < 1) then 4806 if (ierr <= enough) then 4807 ABI_WARNING(sjoin('kpt + qpt:', trim(ktoa(kq)), 'has no symmetric among the k-points!')) 4808 end if 4809 ierr = ierr + 1; cycle 4810 end if 4811 4812 eig_kq(:, ik_bz) = ebands%eig(bstart:bstop, ikq_ibz, spin) 4813 end do 4814 end do 4815 end do 4816 4817 ABI_CHECK(ierr == 0, "See above warnings") 4818 4819 if (any(ltetra == [1, 2])) then 4820 ! Compute weights for double delta integration. Note that libtetra assumes Ef set to zero. 4821 ! TODO: Average weights over degenerate states? 4822 !write(std_out,"(a,i0,2a)")" Calling libtetrabz_dbldelta with ltetra: ", ltetra, " for q-point:", trim(ktoa(qpt)) 4823 eig_k = eig_k - ebands%fermie; eig_kq = eig_kq - ebands%fermie 4824 ABI_MALLOC(wght_bz, (nb, nb, nkbz)) 4825 call libtetrabz_dbldelta(ltetra, cryst%gprimd, nb, nge, eig_k, eig_kq, ngw, wght_bz, comm=comm) 4826 eig_k = eig_k + ebands%fermie !; eig_kq = eig_kq + ebands%fermie 4827 4828 wtqs(iq_ibz, spin) = sum(abs(wght_bz)) 4829 ABI_FREE(wght_bz) 4830 4831 else if (ltetra == 3) then 4832 ! Tetrahedron method with Allen's approach for double delta. 4833 !write(std_out,"(2a)")" Calling Allen's version for q-point: ", trim(ktoa(qpt)) 4834 nene = 3 4835 enemin = ebands%fermie - tol6 4836 enemax = ebands%fermie + tol6 4837 4838 ABI_MALLOC(dtweightde, (nkbz, nene, nene)) 4839 ABI_MALLOC(tweight, (nkbz, nene, nene)) 4840 ABI_MALLOC(work_k, (nkbz)) 4841 ABI_MALLOC(work_kq, (nkbz)) 4842 4843 ! TODO 4844 ABI_ERROR("Not Implemented") 4845 !call init_tetra(indkpt, cryst%gprimd, fs%klatt, kbz, nkbz, tetra, ierr, errorstring, comm) 4846 ABI_CHECK(ierr == 0, "init_tetra returned ierr /= 0") 4847 4848 !fs%dbldelta_tetra_weights_kfs = zero 4849 do ib2=1,nb 4850 work_k = eig_k(ib2, :) 4851 do ib1=1,nb 4852 work_kq = eig_kq(ib1, :) 4853 ! TODO: Average weights over degenerate states? 4854 call get_dbl_tetra_weight(work_k, work_kq, enemin, enemax, enemin, enemax, & 4855 max_occ1, nene, nene, nkbz, tetra, tweight, dtweightde, ierr) 4856 ABI_CHECK(ierr == 0, "get_dbldelta_weights returned ierr /= 0") 4857 4858 wtqs(iq_ibz, spin) = wtqs(iq_ibz, spin) + sum(abs(dtweightde(:, 2, 2))) 4859 end do 4860 end do 4861 4862 call destroy_tetra(tetra) 4863 ABI_FREE(work_k) 4864 ABI_FREE(work_kq) 4865 ABI_FREE(dtweightde) 4866 ABI_FREE(tweight) 4867 else 4868 ABI_ERROR(sjoin("Invalid value of ltetra:", itoa(ltetra))) 4869 end if 4870 end do ! iq_ibz 4871 end do ! spin 4872 4873 ABI_FREE(indkpt) 4874 ABI_FREE(kbz) 4875 ABI_FREE(eig_k) 4876 ABI_FREE(eig_kq) 4877 4878 call ibz_krank%free() 4879 4880 call xmpi_sum(wtqs, comm, ierr) 4881 !call cwtime_report(" calc_dbldelta", cpu, wall, gflops) 4882 4883 end subroutine calc_dbldelta
m_phgamma/eph_phgamma [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
eph_phgamma
FUNCTION
Compute phonon linewidths in metals.
INPUTS
wk0_path=String with the path to the GS unperturbed WFK file. ngfft(18),ngfftf(18)=Coarse and Fine FFT meshes. dtset<dataset_type>=All input variables for this dataset. ebands<ebands_t>=The GS KS band structure (energies, occupancies, k-weights...) dvdb<dbdb_type>=Database with the DFPT SCF potentials. ifc<ifc_type>=interatomic force constants and corresponding real space grid info. pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawang<pawang_type)>=PAW angular mesh and related data. pawrad(ntypat*usepaw)<pawrad_type>=Paw radial mesh and related data. pawtab(ntypat*usepaw)<pawtab_type>=Paw tabulated starting data. psps<pseudopotential_type>=Variables related to pseudopotentials. comm=MPI communicator.
OUTPUT
SOURCE
3086 subroutine eph_phgamma(wfk0_path, dtfil, ngfft, ngfftf, dtset, cryst, ebands, dvdb, ifc, & 3087 pawfgr, pawang, pawrad, pawtab, psps, mpi_enreg, comm) 3088 3089 !Arguments ------------------------------------ 3090 !scalars 3091 character(len=*),intent(in) :: wfk0_path 3092 integer,intent(in) :: comm 3093 type(datafiles_type),intent(in) :: dtfil 3094 type(dataset_type),intent(in) :: dtset 3095 type(crystal_t),intent(in) :: cryst 3096 type(ebands_t),intent(in) :: ebands 3097 type(dvdb_t),intent(inout) :: dvdb 3098 type(pawang_type),intent(in) :: pawang 3099 type(pseudopotential_type),intent(in) :: psps 3100 type(pawfgr_type),intent(in) :: pawfgr 3101 type(ifc_type),intent(in) :: ifc 3102 type(mpi_type),intent(in) :: mpi_enreg 3103 !arrays 3104 integer,intent(in) :: ngfft(18),ngfftf(18) 3105 type(pawrad_type),intent(in) :: pawrad(psps%ntypat*psps%usepaw) 3106 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 3107 3108 !Local variables ------------------------------ 3109 !scalars 3110 integer,parameter :: tim_getgh1c = 1, berryopt0 = 0, ider0 = 0, idir0 = 0, qptopt1 = 1 3111 integer,parameter :: useylmgr = 0, useylmgr1 = 0, master = 0, ndat1 = 1, eph_scalprod0 = 0 3112 integer :: my_rank,nproc,mband,nsppol,nkibz,idir,ipert,iq_ibz, ebands_timrev 3113 integer :: cplex,db_iqpt,natom,natom3,ipc,ipc1,ipc2,nspinor,onpw 3114 integer :: bstart_k,bstart_kq,nband_k,nband_kq,band_k, band_kq, ib_k, ib_kq !ib1,ib2, 3115 integer :: ik_ibz,ik_bz,ikq_bz,ikq_ibz,isym_k,isym_kq,trev_k,trev_kq,timerev_q 3116 integer :: ik_fs, myik, mys !, myiq 3117 integer :: spin,istwf_k,istwf_kq,npw_k,npw_kq 3118 integer :: ii,jj,ipw,mpw,my_mpw,mnb,ierr,cnt,ncid 3119 integer :: n1,n2,n3,n4,n5,n6,nspden,do_ftv1q, ltetra 3120 integer :: sij_opt,usecprj,usevnl,optlocal,optnl,opt_gvnlx1 3121 integer :: nfft,nfftf,mgfft,mgfftf,kq_count,nkpg,nkpg1,edos_intmeth 3122 integer :: jene, iene, comm_rpt, nesting, my_npert, imyp, imyq 3123 integer :: ncerr 3124 real(dp) :: cpu, wall, gflops, cpu_q, wall_q, gflops_q, cpu_k, wall_k, gflops_k, cpu_all, wall_all, gflops_all 3125 real(dp) :: edos_step, edos_broad, sigma, ecut, eshift, eig0nk 3126 logical :: gen_eigenpb, need_velocities, isirr_k, isirr_kq 3127 type(wfd_t) :: wfd 3128 type(fstab_t),pointer :: fs 3129 type(gs_hamiltonian_type) :: gs_hamkq 3130 type(rf_hamiltonian_type) :: rf_hamkq 3131 type(edos_t) :: edos 3132 type(phgamma_t) :: gams 3133 type(a2fw_t) :: a2fw 3134 type(a2fw_tr_t) :: a2fw_tr 3135 type(ddkop_t) :: ddkop 3136 type(xcomm_t) :: pert_comm, qs_comm, qpt_comm, bsum_comm, kpt_comm, spin_comm, pkb_comm !, ncwrite_comm 3137 type(krank_t) :: krank 3138 character(len=500) :: msg 3139 character(len=fnlen) :: path 3140 !arrays 3141 integer :: g0_k(3),g0bz_kq(3),g0_kq(3),symq(4,2,cryst%nsym) 3142 integer :: work_ngfft(18),gmax(3),my_gmax(3),gamma_ngqpt(3) !g0ibz_kq(3), 3143 integer :: indkk_kq(6,1) 3144 integer,allocatable :: kg_k(:,:),kg_kq(:,:),gtmp(:,:),nband(:,:),wfd_istwfk(:) 3145 integer,allocatable :: my_pinfo(:,:), pert_table(:,:) !, qibz_done(:) 3146 real(dp) :: kk(3),kq(3),kk_ibz(3),kq_ibz(3),qpt(3), lf(2),rg(2),res(2), vk(3), vkq(3) 3147 real(dp) :: wminmax(2), n0(ebands%nsppol) 3148 real(dp) :: resvv_in(2,9), resvv_out(2,9), phfrq(3*cryst%natom) 3149 real(dp) :: ylmgr_dum(1,1,1) 3150 real(dp),allocatable :: displ_cart(:,:,:,:), displ_red(:,:,:,:) 3151 real(dp),allocatable :: grad_berry(:,:), kinpw1(:), kpg1_k(:,:), kpg_k(:,:), dkinpw(:) 3152 real(dp),allocatable :: ffnlk(:,:,:,:), ffnl1(:,:,:,:), ph3d(:,:,:), ph3d1(:,:,:) 3153 real(dp),allocatable :: v1scf(:,:,:,:), tgam(:,:,:), gkk_atm(:,:,:,:) !,gkq_nu(:,:,:,:) 3154 real(dp),allocatable :: bras_kq(:,:,:), kets_k(:,:,:), h1kets_kq(:,:,:), cgwork(:,:) 3155 real(dp),allocatable :: ph1d(:,:), vlocal(:,:,:,:), vlocal1(:,:,:,:,:) 3156 real(dp),allocatable :: ylm_kq(:,:), ylm_k(:,:), ylmgr_kq(:,:,:) 3157 real(dp),allocatable :: dummy_vtrial(:,:), gvnlx1(:,:), work(:,:,:,:) 3158 real(dp),allocatable :: gs1c(:,:), v1_work(:,:,:,:), vcar_ibz(:,:,:,:) 3159 real(dp),allocatable :: wt_ek(:,:), wt_ekq(:,:), dbldelta_wts(:,:) 3160 real(dp), allocatable :: tgamvv_in(:,:,:,:), vv_kk(:,:,:), tgamvv_out(:,:,:,:), vv_kkq(:,:,:) 3161 real(dp), allocatable :: tmp_vals_ee(:,:,:,:,:), emesh(:) 3162 logical,allocatable :: bks_mask(:,:,:),keep_ur(:,:,:) 3163 type(fstab_t),target,allocatable :: fstab(:) 3164 type(pawcprj_type),allocatable :: cwaveprj0(:,:) 3165 real(dp) :: abc(3) 3166 #ifdef HAVE_MPI 3167 integer :: ndims, comm_cart, me_cart 3168 logical :: reorder 3169 integer :: coords(5) 3170 integer,allocatable :: dims(:) 3171 logical,allocatable :: periods(:), keepdim(:) 3172 #endif 3173 3174 !************************************************************************ 3175 3176 if (psps%usepaw == 1) then 3177 ABI_ERROR("PAW not implemented") 3178 ABI_UNUSED((/pawang%nsym, pawrad(1)%mesh_size/)) 3179 end if 3180 3181 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 3182 call cwtime(cpu_all, wall_all, gflops_all, "start") 3183 3184 ! Copy important dimensions 3185 natom = cryst%natom; natom3 = 3 * natom; nsppol = ebands%nsppol; nspinor = ebands%nspinor; nspden = dtset%nspden 3186 nkibz = ebands%nkpt; mband = ebands%mband 3187 ebands_timrev = kpts_timrev_from_kptopt(ebands%kptopt) 3188 3189 ! FFT meshes 3190 nfftf = product(ngfftf(1:3)); mgfftf = maxval(ngfftf(1:3)) 3191 nfft = product(ngfft(1:3)) ; mgfft = maxval(ngfft(1:3)) 3192 n1 = ngfft(1); n2 = ngfft(2); n3 = ngfft(3) 3193 n4 = ngfft(4); n5 = ngfft(5); n6 = ngfft(6) 3194 3195 ! Compute electron DOS. 3196 edos_intmeth = 2; if (dtset%prtdos /= 0) edos_intmeth = dtset%prtdos 3197 edos_step = dtset%dosdeltae; edos_broad = dtset%tsmear 3198 edos_step = 0.01 * eV_Ha; edos_broad = 0.3 * eV_Ha 3199 edos = ebands_get_edos(ebands, cryst, edos_intmeth, edos_step, edos_broad, comm) 3200 3201 ! Store DOS per spin channel 3202 n0(:) = edos%gef(1:edos%nsppol) 3203 if (my_rank == master) then 3204 call edos%print(unit=ab_out) 3205 path = strcat(dtfil%filnam_ds(4), "_EDOS") 3206 call wrtout(ab_out, sjoin("- Writing electron DOS to file:", path, ch10)) 3207 call edos%write(path) 3208 end if 3209 3210 ! Find Fermi surface k-points 3211 ! TODO: support kptopt, change setup of k-points if tetra: fist tetra weights then k-points on the Fermi surface! 3212 ABI_MALLOC(fstab, (nsppol)) 3213 call fstab_init(fstab, ebands, cryst, dtset, comm) 3214 if (my_rank == master) then 3215 call fstab_print(fstab, unit=std_out) 3216 call fstab_print(fstab, unit=ab_out) 3217 end if 3218 3219 ! Define q-mesh. eph_ngqpt_fine activates the Fourier interpolation of the DFPT potentials. 3220 gamma_ngqpt = ifc%ngqpt; if (all(dtset%eph_ngqpt_fine /= 0)) gamma_ngqpt = dtset%eph_ngqpt_fine 3221 3222 call phgamma_init(gams, cryst, ifc, ebands, fstab(1), dtset, eph_scalprod0, gamma_ngqpt, n0, comm) 3223 3224 call wrtout(std_out, sjoin("q-mesh for the phonon linewidths:", ltoa(gamma_ngqpt))) 3225 call wrtout(std_out, sjoin("Will compute", itoa(gams%nqibz), "q-points in the IBZ")) 3226 3227 ! Select option for double delta with tetra. 3228 ! 2 for the optimized tetrahedron method. 3229 ! -2 for the linear tetrahedron method. 3230 ltetra = 0 3231 if (dtset%eph_intmeth == 2) ltetra = 2 3232 if (dtset%eph_intmeth == -2) ltetra = 1 3233 3234 ! Compute optimal energy window for tetra. 3235 !call find_ewin(gams%nqibz, gams%qibz, cryst, ebands, ltetra, sigma, comm) 3236 3237 ! ======================== 3238 ! === MPI DISTRIBUTION === 3239 ! ======================== 3240 ! Init for sequential execution. 3241 my_npert = natom3 3242 3243 if (any(dtset%eph_np_pqbks /= 0)) then 3244 ! Use parameters from input file. 3245 pert_comm%nproc = dtset%eph_np_pqbks(1) 3246 qpt_comm%nproc = dtset%eph_np_pqbks(2) 3247 bsum_comm%nproc = dtset%eph_np_pqbks(3) 3248 kpt_comm%nproc = dtset%eph_np_pqbks(4) 3249 spin_comm%nproc = dtset%eph_np_pqbks(5) 3250 my_npert = natom3 / pert_comm%nproc 3251 ABI_CHECK(my_npert > 0, "pert_comm_nproc cannot be greater than 3 * natom.") 3252 ABI_CHECK(mod(natom3, pert_comm%nproc) == 0, "pert_comm_nproc must divide 3 * natom.") 3253 else 3254 ! Automatic grid generation 3255 ! Not the most efficient distribution if large number of MPI procs. 3256 ! TODO: Add parallelism over perturbations although it's less efficient than the parallelism over k 3257 ! It starts to be interesting if we implement symmetries in the k-integration though. 3258 3259 ! TODO: Spin 3260 ! Automatic grid generation over q-points and spins. 3261 !if (new%nsppol == 2 .and. mod(nprocs, 2) == 0) then 3262 ! spin_comm%nproc = 2 3263 ! new%qpt_comm%nproc = nprocs / 2 3264 !else 3265 ! new%qpt_comm%nproc = nprocs 3266 !end if 3267 3268 ! Handle parallelism over perturbations first. 3269 ! Use MPI communicator to distribute the 3 * natom perturbations to reduce memory requirements for DFPT potentials. 3270 ! Ideally, perturbations are equally distributed --> total number of CPUs should be divisible by 3 * natom. 3271 ! or at least, divisible by one integer i for i in [2, 3 * natom - 1]. 3272 3273 if (nsppol == 2 .and. mod(nproc, 2) == 0) then 3274 spin_comm%nproc = 2 3275 kpt_comm%nproc = nproc / 2 3276 else 3277 3278 ! Try to have 3 perts per proc first because the q-point parallelism is more efficient. 3279 ! The memory for W(R,r,ipert) will increase though. 3280 !do cnt=natom,2,-1 3281 ! if (mod(nprocs, cnt) == 0 .and. mod(natom3, cnt) == 0) then 3282 ! pert_comm%nproc = cnt; new%my_npert = natom3 / cnt; exit 3283 ! end if 3284 !end do 3285 3286 !if (pert_comm%nproc == 1) then 3287 ! ! Try again with more procs. 3288 ! do cnt=natom3,2,-1 3289 ! if (mod(nprocs, cnt) == 0 .and. mod(natom3, cnt) == 0) then 3290 ! pert_comm%nproc = cnt; new%my_npert = natom3 / cnt; exit 3291 ! end if 3292 ! end do 3293 !end if 3294 3295 !if (new%my_npert == natom3 .and. nprocs > 1) then 3296 ! ABI_WARNING("The number of MPI procs should be divisible by 3*natom to reduce memory requirements!") 3297 !end if 3298 3299 kpt_comm%nproc = nproc 3300 end if 3301 end if 3302 3303 ! Consistency check. 3304 if (pert_comm%nproc * qpt_comm%nproc * bsum_comm%nproc * kpt_comm%nproc * spin_comm%nproc /= nproc) then 3305 write(msg, "(a,i0,3a, 6(a,1x,i0))") & 3306 "Cannot create 5d Cartesian grid with total nproc: ", nproc, ch10, & 3307 "Idle processes are not supported. The product of the `nproc_*` vars should be equal to nproc.", ch10, & 3308 "pert_nproc (", pert_comm%nproc, ") x qpt_nproc (", qpt_comm%nproc, ") x bsum_nproc (", bsum_comm%nproc, & 3309 ") x kcalc_nproc (", kpt_comm%nproc, ") x spin_nproc (", spin_comm%nproc, ") != ", & 3310 pert_comm%nproc * qpt_comm%nproc * bsum_comm%nproc * kpt_comm%nproc * spin_comm%nproc 3311 ABI_ERROR(msg) 3312 end if 3313 3314 ABI_CHECK(bsum_comm%nproc == 1, "Band parallelism not implemented in m_phgamma") 3315 3316 #ifdef HAVE_MPI 3317 ! Create 5d cartesian communicator: 3*natom perturbations, q-points in IBZ, bands in FS, kpoints in FS, spins 3318 ndims = 5 3319 ABI_MALLOC(dims, (ndims)) 3320 ABI_MALLOC(periods, (ndims)) 3321 ABI_MALLOC(keepdim, (ndims)) 3322 periods(:) = .False.; reorder = .False. 3323 dims = [pert_comm%nproc, qpt_comm%nproc, bsum_comm%nproc, kpt_comm%nproc, spin_comm%nproc] 3324 3325 call MPI_CART_CREATE(comm, ndims, dims, periods, reorder, comm_cart, ierr) 3326 ! Find the index and coordinates of the current processor 3327 call MPI_COMM_RANK(comm_cart, me_cart, ierr) 3328 call MPI_CART_COORDS(comm_cart, me_cart, ndims, coords, ierr) 3329 3330 ! Create communicator to distribute natom3 perturbations. 3331 keepdim = .False.; keepdim(1) = .True. 3332 call MPI_CART_SUB(comm_cart, keepdim, pert_comm%value, ierr); pert_comm%me = xmpi_comm_rank(pert_comm%value) 3333 3334 ! Create communicator for qpoints in self-energy integration. 3335 keepdim = .False.; keepdim(2) = .True. 3336 call MPI_CART_SUB(comm_cart, keepdim, qpt_comm%value, ierr); qpt_comm%me = xmpi_comm_rank(qpt_comm%value) 3337 3338 ! Create communicator for bands for band summation 3339 keepdim = .False.; keepdim(3) = .True. 3340 call MPI_CART_SUB(comm_cart, keepdim, bsum_comm%value, ierr); bsum_comm%me = xmpi_comm_rank(bsum_comm%value) 3341 3342 ! Create communicator for kpoints. 3343 keepdim = .False.; keepdim(4) = .True. 3344 call MPI_CART_SUB(comm_cart, keepdim, kpt_comm%value, ierr); kpt_comm%me = xmpi_comm_rank(kpt_comm%value) 3345 3346 ! Create communicator for spins. 3347 keepdim = .False.; keepdim(5) = .True. 3348 call MPI_CART_SUB(comm_cart, keepdim, spin_comm%value, ierr); spin_comm%me = xmpi_comm_rank(spin_comm%value) 3349 3350 ! Create communicator for the (qpoint, spin) loops 3351 keepdim = .False.; keepdim(2) = .True.; keepdim(5) = .True. 3352 call MPI_CART_SUB(comm_cart, keepdim, qs_comm%value, ierr); qs_comm%me = xmpi_comm_rank(qs_comm%value) 3353 3354 ! Create communicator for the (perturbation, k-point, band_sum) 3355 keepdim = .False.; keepdim(1) = .True.; keepdim(3:4) = .True. 3356 call MPI_CART_SUB(comm_cart, keepdim, pkb_comm%value, ierr); pkb_comm%me = xmpi_comm_rank(pkb_comm%value) 3357 3358 ABI_FREE(dims) 3359 ABI_FREE(periods) 3360 ABI_FREE(keepdim) 3361 call xmpi_comm_free(comm_cart) 3362 #endif 3363 3364 ! Distribute spins and q-points in IBZ. 3365 call xmpi_split_cyclic(nsppol, spin_comm%value, gams%my_nspins, gams%my_spins) 3366 ABI_CHECK(gams%my_nspins > 0, sjoin("nsppol (", itoa(nsppol), ") < spin_comm_nproc (", itoa(spin_comm%nproc), ")")) 3367 3368 call xmpi_split_cyclic(gams%nqibz, qpt_comm%value, gams%my_nqibz, gams%my_iqibz) 3369 ABI_CHECK(gams%my_nqibz > 0, sjoin("nqibz (", itoa(gams%nqibz), ") < qpt_comm_nproc (", itoa(qpt_comm%nproc), ")")) 3370 3371 !ABI_CALLOC(qibz_done, (gams%nqibz)) 3372 3373 path = strcat(dtfil%filnam_ds(4), "_A2F.nc") 3374 ncid = nctk_noid 3375 if (my_rank == master) then 3376 3377 write(std_out, "(/,a)")" === MPI parallelism ===" 3378 !write(std_out, "(2(a,i0))")"P Allocating and summing bands from my_bsum_start: ", self%my_bsum_start, & 3379 ! " up to my_bsum_stop: ", self%my_bsum_stop 3380 write(std_out, "(a,i0)")"P Number of CPUs for parallelism over perturbations: ", pert_comm%nproc 3381 write(std_out, "(a,i0)")"P Number of perturbations treated by this CPU: ", my_npert 3382 write(std_out, "(a,i0)")"P Number of CPUs for parallelism over q-points: ", qpt_comm%nproc 3383 !write(std_out, "(2(a,i0))")"P Number of q-points in the IBZ treated by this proc: " , & 3384 ! count(self%itreat_qibz == 1), " of ", self%nqibz 3385 write(std_out, "(a,i0)")"P Number of CPUs for parallelism over bands: ", bsum_comm%nproc 3386 write(std_out, "(a,i0)")"P Number of CPUs for parallelism over spins: ", spin_comm%nproc 3387 write(std_out, "(a,i0)")"P Number of CPUs for parallelism over k-points: ", kpt_comm%nproc 3388 !write(std_out, "(2(a,i0),/)")"P Number of k-point in Sigma_nk treated by this proc: ", self%my_nkcalc, " of ", self%nkcalc 3389 3390 ! Master creates the netcdf file used to store the results of the calculation. 3391 NCF_CHECK(nctk_open_create(ncid, path, xmpi_comm_self)) 3392 NCF_CHECK(cryst%ncwrite(ncid)) 3393 NCF_CHECK(ebands_ncwrite(ebands, ncid)) 3394 NCF_CHECK(edos%ncwrite(ncid)) 3395 3396 ! Add dimensions. 3397 ncerr = nctk_def_dims(ncid, [nctkdim_t("nqibz", gams%nqibz), nctkdim_t("natom3", natom3)], defmode=.True.) 3398 NCF_CHECK(ncerr) 3399 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "eph_intmeth", "eph_transport", "symdynmat"]) 3400 NCF_CHECK(ncerr) 3401 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "ph_intmeth"]) 3402 NCF_CHECK(ncerr) 3403 ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "eph_fsewin", "eph_fsmear", "eph_extrael", "eph_fermie"]) 3404 NCF_CHECK(ncerr) 3405 ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "ph_wstep", "ph_smear"]) 3406 NCF_CHECK(ncerr) 3407 3408 ! Define arrays for results. 3409 ncerr = nctk_def_arrays(ncid, [ & 3410 nctkarr_t("ngqpt", "int", "three"), & 3411 nctkarr_t("eph_ngqpt_fine", "int", "three"), & 3412 nctkarr_t("ddb_ngqpt", "int", "three"), & 3413 nctkarr_t("ph_ngqpt", "int", "three"), & 3414 nctkarr_t('qibz', "dp", "number_of_reduced_dimensions, nqibz"), & 3415 nctkarr_t('wtq', "dp", "nqibz"), & 3416 !nctkarr_t('qibz_done', "int", "nqibz"), & 3417 nctkarr_t('phfreq_qibz', "dp", "natom3, nqibz"), & 3418 nctkarr_t('phdispl_cart_qibz', "dp", "two, natom3, natom3, nqibz"), & 3419 nctkarr_t('phgamma_qibz', "dp", "natom3, nqibz, number_of_spins"), & 3420 nctkarr_t('phlambda_qibz', "dp", "natom3, nqibz, number_of_spins") & 3421 ]) 3422 NCF_CHECK(ncerr) 3423 3424 ! ====================================================== 3425 ! Write data that do not depend on the (kpt, spin) loop. 3426 ! ====================================================== 3427 NCF_CHECK(nctk_set_datamode(ncid)) 3428 3429 ncerr = nctk_write_iscalars(ncid, & 3430 [character(len=nctk_slen) :: "eph_intmeth", "eph_transport", "symdynmat", "ph_intmeth"], & 3431 [dtset%eph_intmeth, dtset%eph_transport, dtset%symdynmat, dtset%ph_intmeth]) 3432 NCF_CHECK(ncerr) 3433 ncerr = nctk_write_dpscalars(ncid, & 3434 [character(len=nctk_slen) :: "eph_fsewin", "eph_fsmear", "eph_extrael", "eph_fermie", "ph_wstep", "ph_smear"], & 3435 [dtset%eph_fsewin, dtset%eph_fsmear, dtset%eph_extrael, dtset%eph_fermie, dtset%ph_wstep, dtset%ph_smear]) 3436 NCF_CHECK(ncerr) 3437 3438 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'qibz'), gams%qibz)) 3439 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'wtq'), gams%wtq)) 3440 !NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'qibz_done'), qibz_done)) 3441 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ngqpt"), gamma_ngqpt)) 3442 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "eph_ngqpt_fine"), dtset%eph_ngqpt_fine)) 3443 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ph_ngqpt"), dtset%ph_ngqpt)) 3444 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "ddb_ngqpt"), dtset%ddb_ngqpt)) 3445 !NCF_CHECK(nf90_close(ncid)) 3446 end if 3447 3448 !call xmpi_barrier(comm) 3449 3450 ! Now reopen the file inside ncwrite_comm to perform parallel-IO (required for q-point parallelism). 3451 !if (self%ncwrite_comm%value /= xmpi_comm_null) then 3452 ! NCF_CHECK(nctk_open_modify(self%ncid, path, self%ncwrite_comm%value)) 3453 ! NCF_CHECK(nctk_set_datamode(self%ncid)) 3454 !end if 3455 call edos%free() 3456 3457 ! Open the DVDB file 3458 call dvdb%open_read(ngfftf, xmpi_comm_self) 3459 3460 if (pert_comm%nproc > 1) then 3461 ! Activate parallelism over perturbations 3462 ! Build table with list of perturbations treated by this CPU inside pert_comm 3463 call ephtk_set_pertables(cryst%natom, my_npert, pert_table, my_pinfo, pert_comm%value) 3464 call dvdb%set_pert_distrib(my_npert, natom3, my_pinfo, pert_table, pert_comm%value) 3465 ABI_FREE(my_pinfo) 3466 ABI_FREE(pert_table) 3467 end if 3468 3469 ! Check whether all q-points are available in the input DVDB file. 3470 do_ftv1q = 0 3471 do iq_ibz=1,gams%nqibz 3472 qpt = gams%qibz(:,iq_ibz) 3473 if (dvdb%findq(qpt) == -1) do_ftv1q = do_ftv1q + 1 3474 do spin=1,nsppol 3475 fs => fstab(spin) 3476 kq_count = 0 3477 do ik_bz=1,fs%nkfs 3478 kk = fs%kpts(:, ik_bz); kq = kk + qpt 3479 if (fs%findkg0(kq, g0bz_kq) == -1) cycle 3480 kq_count = kq_count + 1 3481 end do 3482 !write(std_out,"((a,i0,2a,a,i0))")" For spin: ",spin,", qpt: ",trim(ktoa(qpt)),", number of (k,q) pairs: ",kq_count 3483 end do 3484 end do 3485 call wrtout(std_out, " ", do_flush=.True.) 3486 3487 if (do_ftv1q /= 0) then 3488 call wrtout([std_out, ab_out], " Cannot find eph_ngqpt_fine q-points in DVDB --> Activating Fourier interpolation.") 3489 ! Prepare Fourier interpolation of DFPT potentials. 3490 comm_rpt = xmpi_comm_self 3491 !comm_rpt = bqs_comm%value 3492 call dvdb%ftinterp_setup(dtset%ddb_ngqpt, qptopt1, 1, dtset%ddb_shiftq, nfftf, ngfftf, comm_rpt) 3493 end if 3494 3495 ! Initialize the wave function descriptor. 3496 ! Only wavefunctions on the FS are stored in wfd. 3497 ! Need all k-points on the FS because of k+q, spin is not distributed for the time being. 3498 ! It would be possible to reduce the memory allocated per MPI-rank via OpenMP. 3499 3500 ABI_MALLOC(nband, (nkibz, nsppol)) 3501 ABI_MALLOC(bks_mask, (mband, nkibz, nsppol)) 3502 ABI_MALLOC(keep_ur, (mband, nkibz ,nsppol)) 3503 nband = mband; bks_mask = .False.; keep_ur = .False. 3504 3505 do mys=1,gams%my_nspins 3506 spin = gams%my_spins(mys) 3507 fs => fstab(spin) 3508 do ik_bz=1,fs%nkfs 3509 ik_ibz = fs%indkk_fs(1, ik_bz) 3510 bstart_k = fs%bstart_cnt_ibz(1, ik_ibz); nband_k = fs%bstart_cnt_ibz(2, ik_ibz) 3511 bks_mask(bstart_k:bstart_k+nband_k-1, ik_ibz, spin) = .True. 3512 end do 3513 end do 3514 !bks_mask(1:mband,:,:) = .True. ! no memory distribution, each node has the full set of states. 3515 3516 ! Impose istwfk = 1 for all k-points. This is also done in respfn (see inkpts) 3517 ! wfd_read_wfk will handle a possible conversion if WFK contains istwfk /= 1. 3518 ABI_MALLOC(wfd_istwfk, (nkibz)) 3519 wfd_istwfk = 1 3520 3521 ecut = dtset%ecut ! dtset%dilatmx 3522 call wfd_init(wfd, cryst, pawtab, psps, keep_ur, mband, nband, nkibz, nsppol, bks_mask,& 3523 nspden, nspinor, ecut, dtset%ecutsm, dtset%dilatmx, wfd_istwfk, ebands%kptns, ngfft,& 3524 dtset%nloalg, dtset%prtvol, dtset%pawprtvol, comm) 3525 3526 call wfd%print(header="Wavefunctions on the Fermi surface") 3527 3528 ABI_FREE(nband) 3529 ABI_FREE(keep_ur) 3530 ABI_FREE(wfd_istwfk) 3531 3532 ! Read wavefunctions. 3533 call wfd%read_wfk(wfk0_path, iomode_from_fname(wfk0_path)) 3534 3535 ! one-dimensional structure factor information on the coarse grid. 3536 ABI_MALLOC(ph1d, (2, 3*(2*mgfft+1)*natom)) 3537 call getph(cryst%atindx, natom, n1, n2, n3, ph1d, cryst%xred) 3538 3539 ! mpw is the maximum number of plane-waves over k and k+q where k and k+q are in the BZ. 3540 ! we also need the max components of the G-spheres (k, k+q) in order to allocate the workspace array work 3541 ! that will be used to symmetrize the wavefunctions in G-space. 3542 call cwtime(cpu, wall, gflops, "start") 3543 call wrtout(std_out, " Computing mpw. This may take some time for dense k/q meshes...") 3544 !call fstab_get_mpw_gmax(nsppol, fstab, mpw, gmax, comm) 3545 mpw = 0; gmax = 0; cnt = 0 3546 3547 ! FIXME: This is an hotspot due to the loop over nkfs. 3548 ! Should use a geometric approach to compute mpw gmax. 3549 do spin=1,nsppol 3550 fs => fstab(spin) 3551 do ik_bz=1,fs%nkfs 3552 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism inside comm 3553 kk = fs%kpts(:, ik_bz) 3554 ! Compute G sphere, returning npw. Note istwfk == 1. 3555 call get_kg(kk, 1, ecut, cryst%gmet, onpw, gtmp) 3556 mpw = max(mpw, onpw) 3557 do ipw=1,onpw 3558 do ii=1,3 3559 gmax(ii) = max(gmax(ii), abs(gtmp(ii,ipw))) 3560 end do 3561 end do 3562 ABI_FREE(gtmp) 3563 3564 do iq_ibz=1,gams%nqibz 3565 qpt = gams%qibz(:,iq_ibz) 3566 ! TODO: g0 umklapp here can enter into play! 3567 ! fstab should contains the max of the umlapp G-vectors. 3568 ! gmax could not be large enough! 3569 kq = kk + qpt 3570 call get_kg(kq, 1, ecut, cryst%gmet, onpw, gtmp) 3571 mpw = max(mpw, onpw) 3572 do ipw=1,onpw 3573 do ii=1,3 3574 gmax(ii) = max(gmax(ii), abs(gtmp(ii,ipw))) 3575 end do 3576 end do 3577 ABI_FREE(gtmp) 3578 end do 3579 3580 end do ! ik_bz 3581 end do ! spin 3582 3583 my_mpw = mpw; call xmpi_max(my_mpw, mpw, comm, ierr) 3584 my_gmax = gmax; call xmpi_max(my_gmax, gmax, comm, ierr) 3585 call wrtout(std_out, sjoin(' Optimal value of mpw: ', itoa(mpw))) 3586 call cwtime_report(" gmax and mpw", cpu, wall, gflops) 3587 3588 ! Init work_ngfft 3589 gmax = gmax + 4 ! FIXME: this is to account for umklapp 3590 gmax = 2*gmax + 1 3591 call ngfft_seq(work_ngfft, gmax) 3592 !write(std_out,*)"work_ngfft(1:3): ",work_ngfft(1:3) 3593 ABI_MALLOC(work, (2, work_ngfft(4), work_ngfft(5), work_ngfft(6))) 3594 3595 ! Allow PW-arrays dimensioned with mpw 3596 ABI_MALLOC(kg_k, (3, mpw)) 3597 ABI_MALLOC(kg_kq, (3, mpw)) 3598 3599 ! Spherical Harmonics for useylm == 1. 3600 ABI_MALLOC(ylm_k, (mpw, psps%mpsang*psps%mpsang*psps%useylm)) 3601 ABI_MALLOC(ylm_kq, (mpw, psps%mpsang*psps%mpsang*psps%useylm)) 3602 ABI_MALLOC(ylmgr_kq, (mpw, 3, psps%mpsang*psps%mpsang*psps%useylm*useylmgr1)) 3603 3604 ! TODO FOR PAW 3605 usecprj = 0 3606 ABI_MALLOC(cwaveprj0, (natom, nspinor*usecprj)) 3607 3608 ! Prepare call to getgh1c 3609 usevnl = 0 3610 optlocal = 1 ! local part of H^(1) is computed in gh1c=<G|H^(1)|C> 3611 optnl = 2 ! non-local part of H^(1) is totally computed in gh1c=<G|H^(1)|C> 3612 opt_gvnlx1 = 0 ! gvnlx1 is output 3613 ABI_MALLOC(gvnlx1, (2, usevnl)) 3614 ABI_MALLOC(grad_berry, (2, nspinor*(berryopt0/4))) 3615 3616 ! This part is taken from dfpt_vtorho 3617 !==== Initialize most of the Hamiltonian (and derivative) ==== 3618 !1) Allocate all arrays and initialize quantities that do not depend on k and spin. 3619 !2) Perform the setup needed for the non-local factors: 3620 ! Norm-conserving: Constant kleimann-Bylander energies are copied from psps to gs_hamk. 3621 ! PAW: Initialize the overlap coefficients and allocate the Dij coefficients. 3622 3623 call init_hamiltonian(gs_hamkq, psps, pawtab, nspinor, nsppol, nspden, natom, & 3624 dtset%typat, cryst%xred, nfft, mgfft, ngfft, cryst%rprimd, dtset%nloalg, & 3625 comm_atom=mpi_enreg%comm_atom, mpi_atmtab=mpi_enreg%my_atmtab, mpi_spintab=mpi_enreg%my_isppoltab, & 3626 usecprj=usecprj, ph1d=ph1d, nucdipmom=dtset%nucdipmom, gpu_option=dtset%gpu_option) 3627 3628 ! Allocate vlocal. Note nvloc 3629 ! I set vlocal to huge to trigger possible bugs (DFPT routines should not access the data) 3630 ABI_MALLOC(vlocal, (n4, n5, n6, gs_hamkq%nvloc)) 3631 vlocal = huge(one) 3632 3633 ! Allocate work space arrays. 3634 ABI_MALLOC(tgam, (2, natom3, natom3)) 3635 ABI_MALLOC(displ_cart, (2, 3, cryst%natom, natom3)) 3636 ABI_MALLOC(displ_red, (2, 3, cryst%natom, natom3)) 3637 3638 ABI_CALLOC(dummy_vtrial, (nfftf, nspden)) 3639 ! TODO: Save data to netcdf file for each q in IBZ and then read data to build a2Fw once all big 3640 ! datastructures (wfd, dvdb) have been deallocated. 3641 ! As a side effect, one can also implement restart over q-points 3642 3643 ! Create ddkop object to compute group velocities if needed. 3644 ! 3645 ! 1) precompute group velocities in the IBZ and the ihave_ikibz_spin file (common to all procs) 3646 ! 2) Use symmetries to reconstruct v_kq from vcar_ibz 3647 ! 3648 ! NB: All procs store in memory the same set of Bloch states inside the energy window. 3649 3650 ddkop = ddkop_new(dtset, cryst, pawtab, psps, wfd%mpi_enreg, mpw, wfd%ngfft) 3651 3652 call cwtime(cpu, wall, gflops, "start", msg=" Computing v_nk matrix elements for all states on the FS...") 3653 ii = huge(1); jj = -1 3654 do spin=1,nsppol 3655 ii = min(ii, fstab(spin)%bmin) 3656 jj = max(jj, fstab(spin)%bmax) 3657 end do 3658 ABI_CALLOC(vcar_ibz, (3, ii:jj, nkibz, nsppol)) 3659 ABI_MALLOC(cgwork, (2, mpw * wfd%nspinor)) 3660 3661 cnt = 0 3662 do spin=1,nsppol 3663 fs => fstab(spin) 3664 do ik_ibz=1,ebands%nkpt 3665 kk = ebands%kptns(:, ik_ibz) 3666 npw_k = wfd%npwarr(ik_ibz); istwf_k = wfd%istwfk(ik_ibz) 3667 ! NB: The two checks below are global --> all procs will cycle. 3668 if (all(bks_mask(:, ik_ibz, spin) .eqv. .False.)) cycle 3669 if (npw_k == 1) cycle 3670 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle ! MPI parallelism. 3671 3672 call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, kk, istwf_k, npw_k, wfd%kdata(ik_ibz)%kg_k) 3673 3674 do band_k=fs%bmin,fs%bmax 3675 if (.not. bks_mask(band_k, ik_ibz, spin)) cycle 3676 !if (.not. wfd%ihave_ug(band_k, ik_ibz, spin)) cycle 3677 call wfd%copy_cg(band_k, ik_ibz, spin, cgwork) 3678 eig0nk = ebands%eig(band_k, ik_ibz, spin) 3679 vk = ddkop%get_vdiag(eig0nk, istwf_k, npw_k, wfd%nspinor, cgwork, cwaveprj0) 3680 vcar_ibz(:, band_k, ik_ibz, spin) = vk 3681 3682 ! TODO: Use ebands_get_edos_matrix_elements 3683 ! reald(dp) :: vv_fs(3,3,nsppol) 3684 !vv_fs = zero 3685 !sigma = fs%eph_fsmear 3686 !if (fs%eph_fsmear < zero) then 3687 ! sigma = max(maxval([(abs(dot_product(fs%vk(:, ib2), fs%kmesh_cartvec(:,ii))), ii=1,3)]), fs%min_smear) 3688 !end if 3689 !fs_wtk = gaussian(emesh - ebands%eig(band_k, ik_ibz, spin), sigma) / (one * fs%nktot) 3690 !do ii=1,3 3691 ! do jj=1,3 3692 ! vv_fs(ii,jj,spin) = vv_fs(ii,jj,spin) + vk(ii) * vk(jj) * fs_wtk 3693 ! end do 3694 !end do 3695 !call xmpi_sum(vv_fs, comm, ierr) 3696 3697 end do 3698 end do 3699 end do ! spin 3700 3701 call xmpi_sum(vcar_ibz, comm, ierr) 3702 ABI_FREE(cgwork) 3703 call cwtime_report(" Velocities", cpu, wall, gflops) 3704 3705 ABI_FREE(bks_mask) 3706 3707 ! Write v_nk to disk. 3708 ! if (my_rank == master) then 3709 ! NCF_CHECK(nf90_put_var(sigma%ncid, nctk_idname(sigma%ncid, "vcar_ibz"), vcar_ibz)) 3710 ! end if 3711 3712 if (dtset%eph_transport > 0) then 3713 ABI_MALLOC(tgamvv_in, (2, 9, natom3, natom3)) 3714 ABI_MALLOC(tgamvv_out, (2, 9, natom3, natom3)) 3715 end if 3716 3717 ! Build krank object to find k-points 3718 krank = krank_from_kptrlatt(ebands%nkpt, ebands%kptns, ebands%kptrlatt, compute_invrank=.False.) 3719 3720 ! Loop over my q-points in the IBZ. 3721 do imyq=1,gams%my_nqibz 3722 iq_ibz = gams%my_iqibz(imyq) 3723 3724 ! Check if this (kpoint, spin) was already calculated 3725 !if (all(sigma%qp_done(ikcalc, :) == 1)) cycle 3726 3727 call cwtime(cpu_q, wall_q, gflops_q, "start") 3728 3729 qpt = gams%qibz(:, iq_ibz) 3730 msg = sjoin("[", itoa(iq_ibz), "/", itoa(gams%nqibz), "]") 3731 call wrtout(std_out, sjoin(" Computing phonon linewidths for IBZ q-point:", ktoa(qpt), msg)) 3732 3733 tgam = zero 3734 if (dtset%eph_transport > 0) then 3735 tgamvv_in = zero; tgamvv_out = zero 3736 end if 3737 3738 if (do_ftv1q == 0) then 3739 ! No interpolation --> find the index of the q-point in the DVDB. 3740 db_iqpt = dvdb%findq(qpt) 3741 3742 if (db_iqpt /= -1) then 3743 ! Read and reconstruct the dvscf potentials for all 3*natom perturbations. 3744 ! This call allocates v1scf(cplex, nfftf, nspden, 3*natom)) 3745 call dvdb%readsym_allv1(db_iqpt, cplex, nfftf, ngfftf, v1scf, pkb_comm%value) 3746 3747 if (dvdb%my_npert /= natom3) then 3748 ! Extract my npert from v1scf 3749 ABI_MALLOC(v1_work, (cplex, nfftf, nspden, dvdb%my_npert)) 3750 do imyp=1,dvdb%my_npert 3751 v1_work(:,:,:,imyp) = v1scf(:,:,:,dvdb%my_pinfo(3, imyp)) 3752 end do 3753 ABI_MOVE_ALLOC(v1_work, v1scf) 3754 end if 3755 3756 else 3757 ABI_ERROR(sjoin("Cannot find q-point:", ktoa(qpt), "in the DVDB file.")) 3758 end if 3759 3760 else 3761 ! Use Fourier interpolation of DFPT potentials to get my_npert potentials. 3762 cplex = 2 3763 ABI_MALLOC(v1scf, (cplex, nfft, nspden, dvdb%my_npert)) 3764 call dvdb%ftinterp_qpt(qpt, nfftf, ngfftf, v1scf, dvdb%comm_rpt) 3765 end if 3766 3767 ! Examine the symmetries of the q wavevector. 3768 call littlegroup_q(cryst%nsym, qpt, symq, cryst%symrec, cryst%symafm, timerev_q, prtvol=dtset%prtvol) 3769 3770 ! Get phonon frequencies and eigenvectors for this q-point. 3771 call ifc%fourq(cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red) 3772 wminmax(1) = min(wminmax(1), phfrq(1)) 3773 wminmax(2) = max(wminmax(2), phfrq(3*cryst%natom)) 3774 3775 ! Allocate vlocal1 with correct cplex. Note nvloc and my_npert. 3776 ABI_MALLOC(vlocal1, (cplex*n4, n5, n6, gs_hamkq%nvloc, my_npert)) 3777 3778 ! Loop over my spins. 3779 do mys=1,gams%my_nspins 3780 spin = gams%my_spins(mys) 3781 fs => fstab(spin) 3782 3783 if (dtset%prteliash == 3) then 3784 ABI_MALLOC_OR_DIE(tmp_vals_ee, (2, gams%nene, gams%nene, natom3, natom3), ierr) 3785 tmp_vals_ee = zero 3786 ! Energy mesh for electrons in a2F(e,e',w) NB: It depends on the spin through enemin 3787 ABI_MALLOC(emesh, (gams%nene)) 3788 emesh = arth(fs%enemin, fs%deltaene, gams%nene) 3789 end if 3790 3791 ! Set up local potential vlocal1 with proper dimensioning from vtrial1 taking into account the spin. 3792 do imyp=1,my_npert 3793 call rf_transgrid_and_pack(spin, nspden, psps%usepaw, cplex, nfftf, nfft, ngfft, gs_hamkq%nvloc,& 3794 pawfgr, mpi_enreg, dummy_vtrial, v1scf(:,:,:,imyp), vlocal, vlocal1(:,:,:,:,imyp)) 3795 end do 3796 3797 ! Continue to initialize the GS Hamiltonian 3798 call gs_hamkq%load_spin(spin, vlocal=vlocal, with_nonlocal=.true.) 3799 3800 ! Allocate workspace for wavefunctions. Make npw larger than expected. 3801 ! maxnb is the maximum number of bands crossing the FS, used to dimension arrays. 3802 mnb = fs%maxnb 3803 ABI_MALLOC(bras_kq, (2, mpw*nspinor, mnb)) 3804 ABI_MALLOC(kets_k, (2, mpw*nspinor, mnb)) 3805 ABI_MALLOC(h1kets_kq, (2, mpw*nspinor, mnb)) 3806 ABI_MALLOC(gkk_atm, (2, mnb, mnb, natom3)) 3807 !ABI_MALLOC(gkq_nu, (2, mnb, natom3)) 3808 3809 ! The weights for the integration of the double-delta. 3810 ABI_MALLOC(dbldelta_wts, (mnb, mnb)) 3811 3812 if (dtset%eph_transport > 0) then 3813 ABI_MALLOC(vv_kk, (9, mnb, mnb)) 3814 ABI_MALLOC(vv_kkq, (9, mnb, mnb)) 3815 end if 3816 3817 if (dtset%prteliash == 3) then 3818 ABI_MALLOC(wt_ek, (gams%nene, mnb)) 3819 ABI_MALLOC(wt_ekq, (gams%nene, mnb)) 3820 end if 3821 3822 ! ===================================== 3823 ! Integration over the FS for this spin 3824 ! ===================================== 3825 ! Compute integration weights and distribute k-points (gams%my_nfsk_q) 3826 call phgamma_setup_qpoint(gams, fs, cryst, ebands, spin, ltetra, qpt, nesting, kpt_comm%value) 3827 3828 do myik=1,gams%my_nfsk_q 3829 call cwtime(cpu_k, wall_k, gflops_k, "start") 3830 3831 ! The k-point and the symmetries relating the BZ k-point to the IBZ. 3832 ik_fs = gams%my_ifsk_q(myik) 3833 kk = fs%kpts(:, ik_fs) 3834 ik_ibz = fs%indkk_fs(1, ik_fs); isym_k = fs%indkk_fs(2, ik_fs) 3835 trev_k = fs%indkk_fs(6, ik_fs); g0_k = fs%indkk_fs(3:5,ik_fs) 3836 isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0)) 3837 kk_ibz = ebands%kptns(:,ik_ibz) 3838 3839 ! Number of bands crossing the Fermi level at k 3840 bstart_k = fs%bstart_cnt_ibz(1, ik_ibz); nband_k = fs%bstart_cnt_ibz(2, ik_ibz) 3841 3842 ! Find k+q in the extended zone and extract symmetry info. cycle if k+q not in FS. 3843 ! Be careful here because there are two umklapp vectors to be considered: 3844 ! 3845 ! k + q = k_bz + g0_bz = IS(k_ibz) + g0_ibz + g0_bz 3846 ! 3847 kq = kk + qpt; ikq_bz = fs%findkg0(kq, g0bz_kq) 3848 3849 ! Skip this point if kq does not belong to the FS window. 3850 if (ikq_bz == -1) cycle 3851 3852 if (kpts_map("symrel", ebands_timrev, cryst, krank, 1, kq, indkk_kq) /= 0) then 3853 write(msg, '(9a)' ) & 3854 "The WFK file cannot be used to compute phonon linewidths.",ch10, & 3855 "At least one of the k-points on the FS could not be generated from a symmetrical one.", ch10, & 3856 "q-mesh: ", trim(ltoa(gamma_ngqpt)), ", k-mesh (from kptrlatt) ", trim(ltoa(get_diag(ebands%kptrlatt))), & 3857 'Action: check your WFK file and the (k, q) point input variables.' 3858 ABI_ERROR(msg) 3859 end if 3860 3861 ikq_ibz = indkk_kq(1, 1); isym_kq = indkk_kq(2, 1) 3862 trev_kq = indkk_kq(6, 1); g0_kq = indkk_kq(3:5, 1) 3863 isirr_kq = (isym_kq == 1 .and. trev_kq == 0 .and. all(g0_kq == 0)) 3864 kq_ibz = ebands%kptns(:, ikq_ibz) 3865 3866 ! If we have used the KERANGE trick, we may have k or k+q points with just one G component set to zero 3867 ! so we skip this transition immediately. This should happen only if fsewin > sigma_erange. 3868 if (wfd%npwarr(ik_ibz) == 1 .or. wfd%npwarr(ikq_ibz) == 1) cycle 3869 3870 ! Number of bands crossing the Fermi level at k+q 3871 bstart_kq = fs%bstart_cnt_ibz(1, ikq_ibz); nband_kq = fs%bstart_cnt_ibz(2, ikq_ibz) 3872 ABI_CHECK(nband_k <= mnb .and. nband_kq <= mnb, "wrong nband") 3873 3874 ! Get npw_k, kg_k and symmetrize wavefunctions from IBZ (if needed). 3875 call wfd%sym_ug_kg(ecut, kk, kk_ibz, bstart_k, nband_k, spin, mpw, fs%indkk_fs(:,ik_fs), cryst, & 3876 work_ngfft, work, istwf_k, npw_k, kg_k, kets_k) 3877 3878 ! Get npw_kq, kg_kq and symmetrize wavefunctions from IBZ (if needed). 3879 call wfd%sym_ug_kg(ecut, kq, kq_ibz, bstart_kq, nband_kq, spin, mpw, indkk_kq(:,1), cryst, & 3880 work_ngfft, work, istwf_kq, npw_kq, kg_kq, bras_kq) 3881 3882 ! if PAW, one has to solve a generalized eigenproblem 3883 ! Be careful here because I will need sij_opt==-1 3884 gen_eigenpb = psps%usepaw == 1; sij_opt = 0; if (gen_eigenpb) sij_opt = 1 3885 ABI_MALLOC(gs1c, (2, npw_kq*nspinor*((sij_opt+1)/2))) 3886 3887 ! Set up the spherical harmonics (Ylm) at k and k+q. See also dfpt_looppert 3888 !if (psps%useylm == 1) then 3889 ! optder = 0; if (useylmgr == 1) optder = 1 3890 ! call initylmg(cryst%gprimd, kg_k, kk, mkmem1, mpi_enreg, psps%mpsang, mpw, nband, mkmem1,& 3891 ! [npw_k], dtset%nsppol, optder, cryst%rprimd, ylm_k, ylmgr) 3892 ! call initylmg(cryst%gprimd, kg_kq, kq, mkmem1, mpi_enreg, psps%mpsang, mpw, nband, mkmem1,& 3893 ! [npw_kq], dtset%nsppol, optder, cryst%rprimd, ylm_kq, ylmgr_kq) 3894 !end if 3895 3896 ! Compute k+G vectors 3897 nkpg = 3 * dtset%nloalg(3) 3898 ABI_MALLOC(kpg_k, (npw_k, nkpg)) 3899 if (nkpg > 0) call mkkpg(kg_k, kpg_k, kk, nkpg, npw_k) 3900 3901 ! Compute nonlocal form factors ffnlk at (k+G) 3902 ABI_MALLOC(ffnlk, (npw_k, 1, psps%lmnmax, psps%ntypat)) 3903 3904 call mkffnl_objs(cryst, psps, 1, ffnlk, ider0, idir0, kg_k, kpg_k, kk, nkpg, npw_k, ylm_k, ylmgr_dum, & 3905 comm=pert_comm%value) 3906 3907 ! Compute k+q+G vectors 3908 nkpg1 = 3 * dtset%nloalg(3) 3909 ABI_MALLOC(kpg1_k, (npw_kq, nkpg1)) 3910 if (nkpg1 > 0) call mkkpg(kg_kq, kpg1_k, kq, nkpg1, npw_kq) 3911 3912 ! Compute nonlocal form factors ffnl1 at (k+q+G) 3913 ABI_MALLOC(ffnl1, (npw_kq, 1, psps%lmnmax, psps%ntypat)) 3914 3915 call mkffnl_objs(cryst, psps, 1, ffnl1, ider0, idir0, kg_kq, kpg1_k, kq, nkpg1, npw_kq, ylm_kq, ylmgr_kq, & 3916 comm=pert_comm%value) 3917 3918 ! Loop over all my atomic perturbations and compute gkk_atm. 3919 gkk_atm = zero 3920 do imyp=1,my_npert 3921 idir = dvdb%my_pinfo(1, imyp); ipert = dvdb%my_pinfo(2, imyp); ipc = dvdb%my_pinfo(3, imyp) 3922 3923 ! Prepare application of the NL part. 3924 call init_rf_hamiltonian(cplex, gs_hamkq, ipert, rf_hamkq, has_e1kbsc=.true.) 3925 3926 call rf_hamkq%load_spin(spin, vlocal1=vlocal1(:,:,:,:,imyp), with_nonlocal=.true.) 3927 3928 ! This call is not optimal because there are quantities in out that do not depend on idir,ipert 3929 call getgh1c_setup(gs_hamkq, rf_hamkq, dtset, psps, kk, kq, idir, ipert, & ! In 3930 cryst%natom, cryst%rmet, cryst%gprimd, cryst%gmet, istwf_k, & ! In 3931 npw_k, npw_kq, useylmgr1, kg_k, ylm_k, kg_kq, ylm_kq, ylmgr_kq, & ! In 3932 dkinpw, nkpg, nkpg1, kpg_k, kpg1_k, kinpw1, ffnlk, ffnl1, ph3d, ph3d1, & ! Out 3933 reuse_kpg_k=1, reuse_kpg1_k=1, reuse_ffnlk=1, reuse_ffnl1=1) ! Reuse some arrays 3934 3935 ! Calculate dvscf * psi_k, results stored in h1kets_kq on the k+q sphere. 3936 ! Compute H(1) applied to GS wavefunction Psi(0) 3937 do ib_k=1,nband_k 3938 band_k = ib_k + bstart_k - 1 3939 eig0nk = ebands%eig(band_k, ik_ibz, spin) 3940 ! Use scissor shift on 0-order eigenvalue 3941 eshift = eig0nk - dtset%dfpt_sciss 3942 3943 call getgh1c(berryopt0, kets_k(:,:,ib_k), cwaveprj0, h1kets_kq(:,:,ib_k), & 3944 grad_berry, gs1c, gs_hamkq, gvnlx1, idir, ipert, eshift, mpi_enreg, optlocal, & 3945 optnl, opt_gvnlx1, rf_hamkq, sij_opt, tim_getgh1c, usevnl) 3946 end do 3947 3948 call rf_hamkq%free() 3949 3950 ABI_FREE(kinpw1) 3951 ABI_FREE(dkinpw) 3952 ABI_FREE(ph3d) 3953 ABI_SFREE(ph3d1) 3954 3955 ! Calculate elphmat(j,i) = <psi_{k+q,j}|dvscf_q*psi_{k,i}> for this perturbation. 3956 ! No need to handle istwf_kq because it's always 1. 3957 ! The array eig1_k contains: 3958 ! 3959 ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)|u_(band,k)^(0)> (NC psps) 3960 ! <u_(band,k+q)^(0)|H_(k+q,k)^(1)-(eig0_k+eig0_k+q)/2.S^(1)|u_(band,k)^(0)> (PAW) 3961 do ib_k=1,nband_k 3962 do ib_kq=1,nband_kq 3963 gkk_atm(:, ib_kq, ib_k, ipc) = cg_zdotc(npw_kq*nspinor, bras_kq(1,1,ib_kq), h1kets_kq(1,1,ib_k)) 3964 end do 3965 end do 3966 3967 end do ! imyp (loop over my_npert atomic perturbations) 3968 3969 ABI_FREE(gs1c) 3970 ABI_FREE(ffnlk) 3971 ABI_FREE(ffnl1) 3972 ABI_FREE(kpg1_k) 3973 ABI_FREE(kpg_k) 3974 3975 ! Collect gkk_atm inside pert_comm so that all procs can operate on the data. 3976 if (pert_comm%nproc > 1) call xmpi_sum(gkk_atm, pert_comm%value, ierr) 3977 3978 ! Get gkq in the phonon representation. 3979 !call ephtk_gkknu_from_atm(mnb, mnb, 1, natom, gkq_atm, phfrq, displ_red, gkq_nu) 3980 3981 ! Compute group velocities if we are in transport mode or adaptive gaussian or 3982 ! tetrahedron with libtetrabz returning nesting condition. 3983 need_velocities = dtset%eph_transport > 0 .or. fs%eph_fsmear < zero .or. nesting /= 0 3984 3985 if (need_velocities) then 3986 ! Compute diagonal matrix elements of velocity operator with DFPT routines 3987 ! Velocities are in Cartesian coordinates. 3988 ! 3989 ! If k+q is not in the IBZ, we need to recostruct the value by symmetry using v(Sq) = S v(q). 3990 ! Use transpose(R) because we are using the tables for the wavefunctions 3991 ! In this case listkk has been called with symrec and use_symrec=False 3992 ! so q_bz = S^T q_ibz where S is the isym_kq symmetry 3993 3994 !call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, kk, istwf_k, npw_k, kg_k) 3995 do ib_k=1,nband_k 3996 band_k = ib_k + bstart_k - 1 3997 vk = vcar_ibz(:, band_k, ik_ibz, spin) 3998 if (.not. isirr_k) then 3999 vk = matmul(transpose(cryst%symrel_cart(:,:,isym_k)), vk) 4000 if (trev_k /= 0) vk = -vk 4001 end if 4002 !vk = ddkop%get_vdiag(ebands%eig(band_k, ik_ibz, spin), & 4003 ! istwf_k, npw_k, wfd%nspinor, kets_k(:,:,ib_k), cwaveprj0) 4004 fs%vk(:,ib_k) = vk 4005 end do 4006 4007 !call ddkop%setup_spin_kpoint(dtset, cryst, psps, spin, kq, istwf_kq, npw_kq, kg_kq) 4008 do ib_kq=1,nband_kq 4009 band_kq = ib_kq + bstart_kq - 1 4010 vkq = vcar_ibz(:, band_kq, ikq_ibz, spin) 4011 if (.not. isirr_kq) then 4012 vkq = matmul(transpose(cryst%symrel_cart(:,:,isym_kq)), vkq) 4013 if (trev_kq /= 0) vkq = -vkq 4014 end if 4015 !vkq = ddkop%get_vdiag(ebands%eig(band_kq, ikq_ibz, spin), & 4016 ! istwf_kq, npw_kq, wfd%nspinor, bras_kq(:,:,ib_kq), cwaveprj0) 4017 fs%vkq(:,ib_kq) = vkq 4018 end do 4019 end if 4020 4021 ! Compute weights for double delta integration at the Fermi level. 4022 ! Note that we have to call the routine here after the initialization of fs%vk and fs%vkq 4023 ! TODO: Could precompute weights before entering the loop, apply filter and compute MPI-distribution 4024 call fs%get_dbldelta_weights(ebands, ik_fs, ik_ibz, ikq_ibz, spin, nesting, dbldelta_wts) 4025 ! Multiply by the weight of the q-point if we are summing over the IBZ(q). 4026 !dbldelta_wts = dbldelta_wts * wtk_lgq 4027 4028 ! Accumulate results in tgam (sum over FS k-points and bands for this spin). 4029 do ipc2=1,natom3 4030 do ipc1=1,natom3 4031 do ib_k=1,nband_k 4032 do ib_kq=1,nband_kq 4033 lf = gkk_atm(:, ib_kq, ib_k, ipc1) 4034 rg = gkk_atm(:, ib_kq, ib_k, ipc2) 4035 res(1) = lf(1) * rg(1) + lf(2) * rg(2) 4036 res(2) = lf(1) * rg(2) - lf(2) * rg(1) 4037 tgam(:,ipc1,ipc2) = tgam(:,ipc1,ipc2) + res(:) * dbldelta_wts(ib_kq, ib_k) 4038 end do 4039 end do 4040 end do 4041 end do 4042 4043 if (dtset%eph_transport > 0) then 4044 ! TODO: could almost make this a BLAS call plus a reshape... 4045 do ib_k = 1,nband_k 4046 do ib_kq = 1,nband_kq 4047 do ipc2 = 1,3 4048 do ipc1 = 1,3 4049 ! vk x vk 4050 vv_kk(ipc1 + (ipc2-1)*3, ib_kq, ib_k) = fs%vk(ipc1, ib_kq) * fs%vk(ipc2, ib_k) 4051 ! vk x vk+q 4052 vv_kkq(ipc1 + (ipc2-1)*3, ib_kq, ib_k) = fs%vkq(ipc1, ib_kq) * fs%vk(ipc2, ib_k) 4053 end do 4054 end do 4055 end do 4056 end do 4057 4058 ! Accumulate results (sum over FS and bands). 4059 do ipc2=1,natom3 4060 do ipc1=1,natom3 4061 do ib_k=1,nband_k 4062 do ib_kq=1,nband_kq 4063 lf = gkk_atm(:, ib_kq, ib_k, ipc1) 4064 rg = gkk_atm(:, ib_kq, ib_k, ipc2) 4065 ! res was missing in MJV version! 4066 res(1) = lf(1) * rg(1) + lf(2) * rg(2) 4067 res(2) = lf(1) * rg(2) - lf(2) * rg(1) 4068 resvv_in(1,:) = res(1) * vv_kkq(:,ib_kq, ib_k) 4069 resvv_in(2,:) = res(2) * vv_kkq(:,ib_kq, ib_k) 4070 resvv_out(1,:) = res(1) * vv_kk(:,ib_kq, ib_k) 4071 resvv_out(2,:) = res(2) * vv_kk(:,ib_kq, ib_k) 4072 tgamvv_in(:,:,ipc1, ipc2) = tgamvv_in(:,:,ipc1,ipc2) + resvv_in * dbldelta_wts(ib_kq, ib_k) 4073 tgamvv_out(:,:,ipc1, ipc2) = tgamvv_out(:,:,ipc1,ipc2) + resvv_out * dbldelta_wts(ib_kq, ib_k) 4074 end do 4075 end do 4076 end do 4077 end do 4078 4079 end if ! add transport things 4080 4081 if (dtset%prteliash == 3) then 4082 ! 4083 ! Precompute deltas with gaussian (tetra is not supported here 4084 ! also because one should allocate (nene, nene, mband, mband) 4085 !call get_dbldelta_weights_emesh() 4086 do ib_k=1,nband_k 4087 band_k = ib_k + bstart_k - 1 4088 sigma = fs%eph_fsmear 4089 if (fs%eph_fsmear < zero) then 4090 !ori sigma = max(maxval([(abs(dot_product(fs%vk(:, ib_k), fs%kmesh_cartvec(:,ii))), ii=1,3)]), fs%min_smear) 4091 !replace the implicit loop by an explicit one 4092 !workaround works with both ifort and ifx on oneapi 2024 4093 do ii=1,3 4094 abc(ii) = abs(dot_product(fs%vk(:, ib_k), fs%kmesh_cartvec(:,ii))) 4095 end do 4096 sigma = max(maxval(abc), fs%min_smear) 4097 end if 4098 wt_ek(:, ib_k) = gaussian(emesh - ebands%eig(band_k, ik_ibz, spin), sigma) / sqrt(one * fs%nktot) 4099 end do 4100 4101 do ib_kq=1,nband_kq 4102 band_kq = ib_kq + bstart_kq - 1 4103 sigma = fs%eph_fsmear 4104 if (fs%eph_fsmear < zero) then 4105 !ori sigma = max(maxval([(abs(dot_product(fs%vkq(:, ib_kq), fs%kmesh_cartvec(:,ii))), ii=1,3)]), fs%min_smear) 4106 !replace the implicit loop by an explicit one 4107 !workaround works with both ifort and ifx on oneapi 2024 4108 do ii=1,3 4109 abc(ii) = abs(dot_product(fs%vkq(:, ib_kq), fs%kmesh_cartvec(:,ii))) 4110 end do 4111 sigma = max(maxval(abc), fs%min_smear) 4112 end if 4113 wt_ekq(:, ib_kq) = gaussian(emesh - ebands%eig(band_kq, ikq_ibz, spin), sigma) / sqrt(one * fs%nktot) 4114 end do 4115 4116 do ib_k=1,nband_k 4117 do ib_kq=1,nband_kq 4118 do ipc2=1,natom3 4119 do ipc1=1,natom3 4120 lf = gkk_atm(:, ib_kq, ib_k, ipc1) 4121 rg = gkk_atm(:, ib_kq, ib_k, ipc2) 4122 res(1) = lf(1) * rg(1) + lf(2) * rg(2) 4123 res(2) = lf(1) * rg(2) - lf(2) * rg(1) 4124 do jene = 1, gams%nene 4125 do iene = 1, gams%nene 4126 ! TODO: distribute this in procs over q. Make a temp array here for 1 q 4127 ! then mpi sync it and save it only on 1 processor below after mpisum over k 4128 tmp_vals_ee(:,iene,jene,ipc1,ipc2) = tmp_vals_ee(:,iene,jene,ipc1,ipc2) + & 4129 res(:) * wt_ekq(iene, ib_kq) * wt_ek(jene, ib_k) 4130 end do 4131 end do 4132 end do 4133 end do 4134 end do 4135 end do 4136 end if 4137 4138 if (myik < 20 .or. (fs%nkfs > 100 .and. mod(myik, 200) == 0)) then 4139 write(msg,'(4(a,i0),a,f8.2)')" q-point [", iq_ibz, "/", gams%nqibz, "] k-point [", myik, "/", gams%my_nfsk_q, "]" 4140 call cwtime_report(msg, cpu_k, wall_k, gflops_k) 4141 end if 4142 4143 end do ! ik_fs: sum over k-points on the BZ FS for this spin. 4144 4145 ABI_FREE(dbldelta_wts) 4146 ABI_FREE(bras_kq) 4147 ABI_FREE(kets_k) 4148 ABI_FREE(h1kets_kq) 4149 ABI_FREE(gkk_atm) 4150 !ABI_FREE(gkq_nu) 4151 ABI_SFREE(wt_ek) 4152 ABI_SFREE(wt_ekq) 4153 4154 ! Collect tgam values inside kpt comm and divide by the total number of k-points in the full mesh. 4155 call xmpi_sum(tgam, kpt_comm%value, ierr) 4156 4157 ! Symmetrize gamma matrices (default). Note that this call is not executed in elphon! 4158 if (gams%symgamma == 1) call tgamma_symm(cryst, gams%qibz(:,iq_ibz), tgam) 4159 4160 ! Save results in gams. 4161 gams%vals_qibz(:,:,:,iq_ibz, spin) = tgam 4162 4163 if (dtset%eph_transport > 0) then 4164 ABI_FREE(vv_kk) 4165 ABI_FREE(vv_kkq) 4166 call xmpi_sum(tgamvv_in, kpt_comm%value, ierr) 4167 call xmpi_sum(tgamvv_out, kpt_comm%value, ierr) 4168 gams%vals_in_qibz(:,:,:,:,iq_ibz,spin) = tgamvv_in 4169 gams%vals_out_qibz(:,:,:,:,iq_ibz,spin) = tgamvv_out 4170 end if 4171 4172 if (dtset%prteliash == 3) then 4173 call xmpi_sum(tmp_vals_ee, kpt_comm%value, ierr) 4174 !TODO: write arrays to ncfile. 4175 gams%vals_ee(:,:,:,:,:, iq_ibz, spin) = tmp_vals_ee 4176 ABI_FREE(tmp_vals_ee) 4177 ABI_FREE(emesh) 4178 end if 4179 4180 end do ! spin 4181 4182 ABI_FREE(v1scf) 4183 ABI_FREE(vlocal1) 4184 4185 write(msg,'(2(a,i0),a)')" Computation of q-point [", iq_ibz, "/", gams%nqibz, "]" 4186 call cwtime_report(msg, cpu_q, wall_q, gflops_q, end_str=ch10) 4187 end do ! iq_ibz 4188 4189 call cwtime_report(" phonon linewidths k-loop", cpu_all, wall_all, gflops_all, pre_str=ch10, end_str=ch10) 4190 4191 ! Free memory 4192 ABI_FREE(gvnlx1) 4193 ABI_FREE(grad_berry) 4194 ABI_FREE(dummy_vtrial) 4195 ABI_FREE(work) 4196 ABI_FREE(ph1d) 4197 ABI_FREE(vlocal) 4198 ABI_FREE(kg_k) 4199 ABI_FREE(kg_kq) 4200 ABI_FREE(ylm_k) 4201 ABI_FREE(ylm_kq) 4202 ABI_FREE(ylmgr_kq) 4203 ABI_FREE(tgam) 4204 ABI_FREE(displ_cart) 4205 ABI_FREE(displ_red) 4206 !ABI_FREE(qibz_done) 4207 ABI_SFREE(vcar_ibz) 4208 call krank%free() 4209 4210 if (dtset%eph_transport > 0) then 4211 ABI_FREE(tgamvv_in) 4212 ABI_FREE(tgamvv_out) 4213 end if 4214 4215 call pawcprj_free(cwaveprj0) 4216 ABI_FREE(cwaveprj0) 4217 call ddkop%free() 4218 call gs_hamkq%free() 4219 call wfd%free() 4220 do spin=1,ebands%nsppol 4221 call fstab(spin)%free() 4222 end do 4223 ABI_FREE(fstab) 4224 4225 ! Collect results on each node 4226 call xmpi_sum(gams%vals_qibz, qs_comm%value, ierr) 4227 if (dtset%eph_transport > 0) then 4228 call xmpi_sum(gams%vals_in_qibz, qs_comm%value, ierr) 4229 call xmpi_sum(gams%vals_out_qibz, qs_comm%value, ierr) 4230 end if 4231 if (dtset%prteliash == 3) call xmpi_sum(gams%vals_ee, qs_comm%value, ierr) 4232 4233 ! Close the netcdf file then master reopens it 4234 !if (ncwrite_comm%value /= xmpi_comm_null) then 4235 ! NCF_CHECK(nf90_close(ncid)) 4236 !end if 4237 !call xmpi_barrier(comm) 4238 4239 !ncid = nctk_noid 4240 !if (my_rank == master) then 4241 ! NCF_CHECK(nctk_open_modify(ncid, path, xmpi_comm_self)) 4242 ! NCF_CHECK(nctk_set_datamode(ncid)) 4243 !end if 4244 4245 ! Deallocate MPI communicators. 4246 call pert_comm%free() 4247 call qpt_comm%free() 4248 call bsum_comm%free() 4249 call qs_comm%free() 4250 call kpt_comm%free() 4251 call spin_comm%free() 4252 call pkb_comm%free() 4253 !call ncwrite_comm%free() 4254 4255 ! Print gamma(IBZ) to ab_out and ncid 4256 if (my_rank == master) call phgamma_ncwrite(gams, cryst, ifc, ncid) 4257 4258 ! Interpolate linewidths along the q-path. 4259 if (dtset%ph_nqpath <= 0) then 4260 write(msg, '(7a,es16.6,4a)' ) & 4261 'You have not specified a path for the linewidth calculation - no interpolation or output will be done ',ch10,& 4262 'Action: check your input variables ph_nqpath and ph_qpath' 4263 ABI_WARNING(msg) 4264 else 4265 call gams%linwid(cryst, ifc, dtset%ph_ndivsm, dtset%ph_nqpath, dtset%ph_qpath, dtfil%filnam_ds(4), ncid, wminmax, comm) 4266 end if 4267 4268 ! Compute a2Fw using the ab-initio q-points (no interpolation here) 4269 call a2fw_init(a2fw, gams, cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, wminmax, dtset%ph_smear, & 4270 dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, comm, qintp=.False., qptopt=1) 4271 4272 if (my_rank == master) then 4273 call a2fw%write(strcat(dtfil%filnam_ds(4), "_NOINTP"), "_qcoarse", ncid) 4274 if (dtset%prteliash == 3) call a2fw_ee_write(a2fw, strcat(dtfil%filnam_ds(4), "_NOINTP")) 4275 end if 4276 call a2fw%free() 4277 4278 ! Compute a2Fw using Fourier interpolation (R -> q) and ph_ngqpt grid. 4279 call a2fw_init(a2fw, gams, cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, wminmax, dtset%ph_smear, & 4280 dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, comm, qptopt=1) 4281 4282 if (my_rank == master) then 4283 call a2fw%write(dtfil%filnam_ds(4), "_qintp", ncid) 4284 if (dtset%prteliash == 3) call a2fw_ee_write(a2fw, dtfil%filnam_ds(4)) 4285 end if 4286 4287 call a2fw%free() 4288 4289 ! Compute A2fw using Fourier interpolation and full BZ for debugging purposes. 4290 !call a2fw_init(a2fw, gams, cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, wminmax, dtset%ph_smear,& 4291 ! dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, comm, qptopt=3) 4292 !if (my_rank == master) call a2fw%write(strcat(dtfil%filnam_ds(4), "_A2FW_QPTOPT3"), "fake", nctk_noid) 4293 !call a2fw%free() 4294 4295 if (dtset%eph_transport == 1) then 4296 ! Calculate and output transport quantities 4297 4298 ! Compute a2Fw_tr using ab-initio q-points (no interpolation) 4299 call a2fw_tr_init(a2fw_tr, gams, cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, wminmax, dtset%ph_smear,& 4300 dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, comm, qintp=.False., qptopt=1) 4301 4302 if (my_rank == master) call a2fw_tr%write(strcat(dtfil%filnam_ds(4), "_NOINTP"), "_qcoarse", ncid) 4303 call a2fw_tr%free() 4304 4305 ! Compute a2Fw_tr using Fourier interpolation (R --> q) and ph_ngqpt grid 4306 call a2fw_tr_init(a2fw_tr, gams, cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, wminmax, dtset%ph_smear,& 4307 dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, comm, qptopt=1) 4308 4309 if (my_rank == master) call a2fw_tr%write(dtfil%filnam_ds(4), "_qintp", ncid) 4310 call a2fw_tr%free() 4311 4312 ! Compute A2fw_tr using Fourier interpolation and full BZ for debugging purposes. 4313 !call a2fw_tr_init(a2fw_tr, gams, cryst, ifc, dtset%ph_intmeth, dtset%ph_wstep, wminmax, dtset%ph_smear,& 4314 ! dtset%ph_ngqpt, dtset%ph_nqshift, dtset%ph_qshift, comm, qptopt=3) 4315 !if (my_rank == master) call a2fw_tr%write(strcat(dtfil%filnam_ds(4), "_A2FWTR_QPTOPT3"), "fake", nctk_noid) 4316 !call a2fw_tr%free() 4317 end if 4318 4319 if (my_rank == master) then 4320 NCF_CHECK(nf90_close(ncid)) 4321 end if 4322 4323 call gams%free() 4324 4325 end subroutine eph_phgamma
m_phgamma/find_ewin [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
find_ewin
FUNCTION
Use bisection to find the optimal energy window around the Fermi level
SOURCE
4599 subroutine find_ewin(nqibz, qibz, cryst, ebands, ltetra, fs_ewin, comm) 4600 4601 !Arguments ------------------------------------ 4602 integer,intent(in) :: nqibz 4603 type(crystal_t),intent(in) :: cryst 4604 type(ebands_t),intent(in) :: ebands 4605 integer,intent(in) :: ltetra, comm 4606 real(dp),intent(out) :: fs_ewin 4607 !arrays 4608 real(dp),intent(in) :: qibz(3, nqibz) 4609 4610 !Local variables------------------------------- 4611 !scalars 4612 integer,parameter :: master = 0 4613 integer :: my_rank, iq_ibz, spin, ii, unt 4614 real(dp) :: cpu, wall, gflops 4615 character(len=500) :: msg 4616 !arrays 4617 integer :: bstarts(3), bstops(3) 4618 real(dp) :: elows(3), ehighs(3), ewins(3), qsums(3) 4619 real(dp), allocatable :: wtqs(:,:,:) 4620 4621 ! ************************************************************************* 4622 4623 call cwtime(cpu, wall, gflops, "start") 4624 my_rank = xmpi_comm_rank(comm) 4625 4626 unt = -1 4627 if (xmpi_comm_rank(comm) == master) then 4628 if (open_file("tetra.dat", msg, newunit=unt, action="write", form="formatted", status="unknown") /= 0) then 4629 ABI_ERROR(msg) 4630 end if 4631 end if 4632 4633 ! 1 is low, 2 is high, 3 is for the workspace 4634 ABI_MALLOC(wtqs, (nqibz, ebands%nsppol, 3)) 4635 ewins(1) = half * eV_Ha; elows(1) = ebands%fermie - ewins(1); ehighs(1) = ebands%fermie + ewins(1) 4636 call ebands_get_bands_from_erange(ebands, elows(1), ehighs(1), bstarts(1), bstops(1)) 4637 call calc_dbldelta(cryst, ebands, ltetra, bstarts(1), bstops(1), nqibz, qibz, wtqs(:,:,1), comm) 4638 if (unt /= -1) write(unt, *) wtqs(:,:,1) 4639 4640 !ewins(2) = five * eV_Ha; elows(2) = ebands%fermie - ewins(2); ehighs(2) = ebands%fermie + ewins(2) 4641 ewins(2) = two * eV_Ha; elows(2) = ebands%fermie - ewins(2); ehighs(2) = ebands%fermie + ewins(2) 4642 ewins(2) = one * eV_Ha; elows(2) = ebands%fermie - ewins(2); ehighs(2) = ebands%fermie + ewins(2) 4643 call ebands_get_bands_from_erange(ebands, elows(2), ehighs(2), bstarts(2), bstops(2)) 4644 call calc_dbldelta(cryst, ebands, ltetra, bstarts(2), bstops(2), nqibz, qibz, wtqs(:,:,2), comm) 4645 if (unt /= -1) write(unt, *) wtqs(:,:,2) 4646 4647 if (abs(sum(abs(wtqs(:,:,1)) - sum(abs(wtqs(:,:,2))))) / sum(abs(wtqs(:,:,2))) < tol2) then 4648 fs_ewin = ewins(1) 4649 call print_weights_index(1) 4650 goto 100 4651 end if 4652 4653 ! Bisection part. 4654 do 4655 ewins(3) = (ewins(1) + ewins(2)) / two; elows(3) = ebands%fermie - ewins(3); ehighs(3) = ebands%fermie + ewins(3) 4656 call ebands_get_bands_from_erange(ebands, elows(3), ehighs(3), bstarts(3), bstops(3)) 4657 call calc_dbldelta(cryst, ebands, ltetra, bstarts(3), bstops(3), nqibz, qibz, wtqs(:,:,3), comm) 4658 if (unt /= -1) write(unt, *) wtqs(:,:,3) 4659 4660 do ii=1,3 4661 qsums(ii) = sum(abs(wtqs(:,:,ii))) / nqibz 4662 end do 4663 4664 write(std_out,*)qsums(1), ", for ewins(1)_eV", ewins(1) * Ha_eV, " bstart(1)", bstarts(1), " bstopt(1)", bstops(1) 4665 write(std_out,*)qsums(2), ", for ewins(2)_eV", ewins(2) * Ha_eV, " bstart(1)", bstarts(2), " bstopt(1)", bstops(2) 4666 write(std_out,*)qsums(3), ", for ewins(3)_eV", ewins(3) * Ha_eV, " bstart(1)", bstarts(3), " bstopt(1)", bstops(3) 4667 4668 !if (abs(sum(abs(wtqs(:,:,3)) - sum(abs(wtqs(:,:,2))))) / sum(abs(wtqs(:,:,2))) < tol3) then 4669 ! ! If 3 is close to 2, try to reduce the window by moving towards 1 4670 ! ewins(2) = ewins(3); elows(2) = elows(3); ehighs(2) = ehighs(3) 4671 4672 if (abs(qsums(3) - qsums(2)) / qsums(2) < tol2) then 4673 fs_ewin = ewins(3) 4674 call print_weights_index(3) 4675 exit 4676 4677 else 4678 ewins(1) = ewins(3); elows(1) = elows(3); ehighs(1) = ehighs(3) 4679 end if 4680 4681 end do 4682 4683 100 continue 4684 ABI_FREE(wtqs) 4685 if (unt /= -1) close(unt) 4686 4687 call cwtime_report(" find_ewin", cpu, wall, gflops) 4688 4689 contains 4690 4691 subroutine print_weights_index(ind) 4692 4693 integer,intent(in) :: ind 4694 4695 if (my_rank == master) then 4696 write(std_out, "(a, f8.3, a)")" Optimal FS energy window: ", fs_ewin * Ha_eV, " (eV)" 4697 do iq_ibz=1,nqibz 4698 write(std_out, "(2a)")" For q-point:", trim(ktoa(qibz(:, iq_ibz))) 4699 do spin=1,ebands%nsppol 4700 write(std_out, *) wtqs(iq_ibz, spin, ind) 4701 end do 4702 end do 4703 end if 4704 4705 end subroutine print_weights_index 4706 4707 end subroutine find_ewin
m_phgamma/phgamma_eval_qibz [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
phgamma_eval_qibz
FUNCTION
Compute the phonon linewidths for q-points in the IBZ without performing the interpolation.
INPUTS
gams<phgamma_t> cryst<crystal_t>=Crystal structure. ifc<ifc_type>=Interatomic force constants. iq_ibz=Index of the q-point in the IBZ array. spin=Spin index. See notes below.
OUTPUT
phfrq(gams%natom3)=Phonon frequencies gamma_ph(gams%natom3)=Phonon linewidths. lambda_ph(gams%natom3)=coupling strength coefficients displ_cart(2,3,cry%natom,3*cryst%natom)=Phonon displacement in cartesian coordinates. [gamma_ph_ee]
NOTES
If nsppol == 1 and nspinor == 1, lambda and gamma are already summed over the two equivalent spin channels. If nsppol == 2, lambda and gamma are the particual contributions given by the input spin index. Client code is responsible for assembling the final observables by summing over spins.
SOURCE
801 subroutine phgamma_eval_qibz(gams, cryst, ifc, iq_ibz, spin, phfrq, gamma_ph, lambda_ph, displ_cart, gamma_ph_ee) 802 803 !Arguments ------------------------------------ 804 !scalars 805 integer,intent(in) :: iq_ibz,spin 806 class(phgamma_t),intent(inout) :: gams 807 type(crystal_t),intent(in) :: cryst 808 type(ifc_type),intent(in) :: ifc 809 !arrays 810 real(dp),intent(out) :: phfrq(gams%natom3),gamma_ph(gams%natom3),lambda_ph(gams%natom3) 811 real(dp),intent(out) :: displ_cart(2,3,cryst%natom,3*cryst%natom) 812 real(dp),intent(out),optional :: gamma_ph_ee(gams%nene,gams%nene,gams%natom3) 813 814 !Local variables------------------------------- 815 !scalars 816 integer :: natom3, nu1, iene, jene 817 real(dp) :: spinfact 818 !character(len=500) :: msg 819 !arrays 820 real(dp) :: displ_red(2,gams%natom3,gams%natom3), work_qnu(gams%natom3), gam_atm(2,gams%natom3,gams%natom3) 821 822 ! ************************************************************************* 823 824 natom3 = gams%natom3 825 826 ! Get phonon frequencies and eigenvectors. 827 call ifc%fourq(cryst, gams%qibz(:,iq_ibz), phfrq, displ_cart, out_displ_red=displ_red) 828 829 ! Scalar product with the displ_red vectors. 830 ! Note that the factor 1 / (2 * omega_qnu) coming from |g|^2 is absorbed in the expressions below. 831 gam_atm = reshape(gams%vals_qibz(:,:,:,iq_ibz,spin), [2, natom3, natom3]) 832 call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, gamma_ph) 833 834 if (present(gamma_ph_ee) .and. gams%my_iqibz(iq_ibz) /= -1) then 835 do iene = 1, gams%nene 836 do jene = 1, gams%nene 837 gam_atm = reshape(gams%vals_ee(:,jene,iene,:,:,iq_ibz,spin), [2, natom3, natom3]) 838 call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, work_qnu) 839 gamma_ph_ee(jene, iene, :) = work_qnu 840 end do 841 end do 842 end if 843 844 ! Compute lambda 845 ! TODO: check this - looks like a factor of 2 wrt the inline documentation! 846 ! NB: one factor of 2 comes from the phonon propagator and BE factor, 847 ! then you have to be careful with the convention for the Fermi level DOS 848 ! 849 ! spinfact should be 1 for a normal non sppol calculation without spinorbit 850 ! for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2 851 ! for nsppol 2 it should be 0.5 as we have 2 spin channels to sum 852 spinfact = two / (gams%nsppol * gams%nspinor) 853 854 do nu1=1,gams%natom3 855 gamma_ph(nu1) = gamma_ph(nu1) * pi * spinfact 856 lambda_ph(nu1) = zero 857 if (abs(phfrq(nu1)) > EPHTK_WTOL) lambda_ph(nu1) = gamma_ph(nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2) 858 if (present(gamma_ph_ee)) gamma_ph_ee(:,:,nu1) = gamma_ph_ee(:,:,nu1) * pi * spinfact 859 end do 860 861 ! This to avoid spurious results for the acoustic modes as gamma(q) --> 0 and lambda(q) --> 0 for q --> 0. 862 ! Here we set everything to zero when we are inside a sphere of radius EPH_Q0TOL as this acoustic rule 863 ! is not fulfilled due to numerical inaccuracies. 864 if (normv(gams%qibz(:,iq_ibz), cryst%gmet, "G") < EPH_Q0TOL) then 865 gamma_ph(1:3) = zero 866 lambda_ph(1:3) = zero 867 end if 868 869 end subroutine phgamma_eval_qibz
m_phgamma/phgamma_free [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
phgamma_free
FUNCTION
Free the dynamic memory in a <phgamma_t> datatype
SOURCE
486 subroutine phgamma_free(gams) 487 488 !Arguments ------------------------------------ 489 class(phgamma_t),intent(inout) :: gams 490 491 ! ************************************************************************* 492 493 !real 494 ABI_SFREE(gams%n0) 495 ABI_SFREE(gams%qibz) 496 ABI_SFREE(gams%wtq) 497 ABI_SFREE(gams%qbz) 498 ABI_SFREE(gams%rpt) 499 ABI_SFREE(gams%wghatm) 500 ABI_SFREE(gams%vals_qibz) 501 ABI_SFREE(gams%vals_rpt) 502 ABI_SFREE(gams%vals_in_qibz) 503 ABI_SFREE(gams%vals_in_rpt) 504 ABI_SFREE(gams%vals_out_qibz) 505 ABI_SFREE(gams%vals_out_rpt) 506 ABI_SFREE(gams%vals_ee) 507 ABI_SFREE(gams%my_iqibz) 508 ABI_SFREE(gams%my_ifsk_q) 509 ABI_SFREE(gams%my_spins) 510 511 end subroutine phgamma_free
m_phgamma/phgamma_init [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
phgamma_init
FUNCTION
Creation method for the phgamma_t datatype.
INPUTS
cryst<crystal_t> ifc<ifc_type>=Interatomic force constants. symdynmat=1 to activa symmetrization of gamma matrices. ngqpt(3)=Q-mesh divisions n0(dtset%nsppol)=Density of states at the Fermi level per spin.
OUTPUT
gams<phgamma_t>
SOURCE
535 subroutine phgamma_init(gams, cryst, ifc, ebands, fstab, dtset, eph_scalprod, ngqpt, n0, comm) 536 537 !Arguments ------------------------------------ 538 !scalars 539 integer,intent(in) :: eph_scalprod 540 integer,intent(in) :: comm 541 type(crystal_t),intent(in) :: cryst 542 type(ifc_type),intent(in) :: ifc 543 type(phgamma_t),intent(out) :: gams 544 type(ebands_t),intent(in) :: ebands 545 type(fstab_t), intent(in) :: fstab 546 type(dataset_type),intent(in) :: dtset 547 !arrays 548 integer,intent(in) :: ngqpt(3) 549 real(dp),intent(in) :: n0(dtset%nsppol) 550 551 !Local variables------------------------------- 552 !scalars 553 integer :: my_rank, nproc, ierr, nsppol, natom3, qptopt, qtimrev 554 !arrays 555 integer :: qptrlatt(3,3) 556 557 ! ************************************************************************* 558 559 my_rank = xmpi_comm_rank(comm); nproc = xmpi_comm_size(comm) 560 561 ! Set basic dimensions. 562 nsppol = dtset%nsppol; gams%nsppol = nsppol; gams%nspinor = dtset%nspinor 563 564 gams%natom = cryst%natom; gams%natom3 = 3*cryst%natom; natom3 = gams%natom3 565 566 gams%symgamma = dtset%symdynmat; gams%eph_scalprod = eph_scalprod 567 !gams%asr = ifc%asr 568 gams%asr = 0 569 gams%prteliash = dtset%prteliash 570 571 gams%ndir_transp = 0; if (dtset%eph_transport > 0) gams%ndir_transp = 3 572 573 ! Copy DOS. 574 ABI_MALLOC(gams%n0, (nsppol)) 575 gams%n0 = n0 576 577 gams%nene = fstab%nene 578 gams%enemin = fstab%enemin 579 gams%deltaene = fstab%deltaene 580 581 ! Setup IBZ, weights and BZ. Always use q --> -q symmetry for phonons even in systems wo inversion 582 gams%ngqpt = ngqpt 583 qptrlatt = 0; qptrlatt(1,1) = ngqpt(1); qptrlatt(2,2) = ngqpt(2); qptrlatt(3,3) = ngqpt(3) 584 qptopt = ebands%kptopt; if (dtset%qptopt /= 0) qptopt = dtset%qptopt 585 qtimrev = kpts_timrev_from_kptopt(qptopt) 586 587 call kpts_ibz_from_kptrlatt(cryst, qptrlatt, qptopt, 1, [zero, zero, zero], & 588 gams%nqibz, gams%qibz, gams%wtq, gams%nqbz, gams%qbz) 589 590 ! Allocate matrices in the IBZ. 591 ABI_MALLOC_OR_DIE(gams%vals_qibz, (2, natom3, natom3, gams%nqibz, nsppol), ierr) 592 gams%vals_qibz = zero 593 594 if (dtset%eph_transport > 0) then 595 ABI_MALLOC_OR_DIE(gams%vals_in_qibz, (2, 9, natom3, natom3, gams%nqibz, nsppol), ierr) 596 gams%vals_in_qibz = zero 597 ABI_MALLOC_OR_DIE(gams%vals_out_qibz, (2, 9, natom3, natom3, gams%nqibz, nsppol), ierr) 598 gams%vals_out_qibz = zero 599 end if 600 601 if (gams%prteliash == 3) then 602 ABI_MALLOC_OR_DIE(gams%vals_ee, (2, gams%nene, gams%nene, natom3, natom3, gams%nqibz, nsppol), ierr) 603 gams%vals_ee = zero 604 end if 605 606 ! Prepare Fourier interpolation. 607 gams%gprim = ifc%gprim 608 gams%nrpt = ifc%nrpt 609 ABI_MALLOC(gams%rpt, (3, gams%nrpt)) 610 gams%rpt = ifc%rpt 611 ABI_MALLOC(gams%wghatm, (gams%natom, gams%natom, gams%nrpt)) 612 gams%wghatm = ifc%wghatm 613 614 end subroutine phgamma_init
m_phgamma/phgamma_interp [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
phgamma_interp
FUNCTION
Interpolate the phonon linewidths at a given q-point.
INPUTS
gams<phgamma_t> cryst<crystal_t>=crystalline structure. ifc<ifc_type>=Interatomic force constants. spin=Spin index qpt(3)=q-point in reduced coordinates gamma_ph(3*natom)=Phonon linewidths
OUTPUT
gamma_ph(gams%natom3)=Interpolated Phonon linewidths. lamda_ph(3*natom)=Lambda coefficients for the different phonon modes. phfrq(3*natom)=phonon frequencies at current q displ_cart(2,3,natom,3*natom) = Phonon displacement in Cartesian coordinates
SOURCE
897 subroutine phgamma_interp(gams, cryst, ifc, spin, qpt, phfrq, gamma_ph, lambda_ph, displ_cart, gamma_ph_ee) 898 899 !Arguments ------------------------------------ 900 !scalars 901 integer,intent(in) :: spin 902 class(phgamma_t),intent(inout) :: gams 903 type(crystal_t),intent(in) :: cryst 904 type(ifc_type),intent(in) :: ifc 905 !arrays 906 real(dp),intent(in) :: qpt(3) 907 real(dp),intent(out) :: phfrq(gams%natom3),gamma_ph(gams%natom3),lambda_ph(gams%natom3) 908 real(dp),intent(out) :: displ_cart(2,3,cryst%natom,3*cryst%natom) 909 real(dp),intent(out),optional :: gamma_ph_ee(gams%nene,gams%nene,gams%natom3) 910 911 !Local variables------------------------------- 912 !scalars 913 integer,parameter :: qtor0 = 0 914 integer, save :: icall=0 915 integer :: natom3,nu1 916 real(dp) :: spinfact 917 character(len=500) :: msg 918 !arrays 919 real(dp) :: displ_red(2,gams%natom3,gams%natom3) 920 real(dp) :: gam_now(2,gams%natom3**2), gam_atm(2,gams%natom3,gams%natom3) !,work_qnu(gams%natom3) 921 real(dp),allocatable :: coskr(:,:),sinkr(:,:) 922 923 ! ************************************************************************* 924 925 ! Compute internal tables used for Fourier interpolation. 926 if (.not. allocated(gams%vals_rpt)) call gams%interp_setup(cryst) 927 928 if (present(gamma_ph_ee) .and. icall == 0) then 929 gamma_ph_ee = zero 930 write (msg,'(3a)')& 931 " For the moment gams_ee matrix elements are not FT interpolated wrt q,",ch10,& 932 " only evaluated on the electron k grid. The resulting a2feew will be 0" 933 ABI_WARNING(msg) 934 icall = 1 935 end if 936 937 natom3 = gams%natom3 938 939 ! Taken from mkph_linwid 940 ! This reduced version of ftgkk supposes the kpoints have been integrated 941 ! in integrate_gamma. Do FT from real-space gamma grid to 1 qpt. 942 ABI_MALLOC(coskr, (1, gams%nrpt)) 943 ABI_MALLOC(sinkr, (1, gams%nrpt)) 944 ! TODO: This is not optimal 945 call ftgam_init(gams%gprim, 1, gams%nrpt, qpt, gams%rpt, coskr, sinkr) 946 947 call ftgam(gams%wghatm, gam_now, gams%vals_rpt(:,:,:,spin), gams%natom, 1, gams%nrpt, qtor0, coskr, sinkr) 948 949 ! This call is not executed in elphon! 950 if (gams%symgamma == 1) call tgamma_symm(cryst, qpt, gam_now) 951 952 ABI_FREE(coskr) 953 ABI_FREE(sinkr) 954 955 ! Get phonon frequencies and eigenvectors. 956 call ifc%fourq(cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red) 957 958 ! Scalar product with the displ_red 959 gam_atm = reshape(gam_now, [2, natom3, natom3]) 960 call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, gamma_ph) 961 962 ! Compute lambda 963 ! spinfact should be 1 for a normal non sppol calculation without spinorbit 964 ! for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2 965 ! for nsppol 2 it should be 0.5 as we have 2 spin channels to sum 966 spinfact = two / (gams%nsppol * gams%nspinor) 967 968 ! Compute lambda 969 do nu1=1,gams%natom3 970 gamma_ph(nu1) = gamma_ph(nu1) * pi * spinfact 971 lambda_ph(nu1) = zero 972 if (abs(phfrq(nu1)) > EPHTK_WTOL) lambda_ph(nu1) = gamma_ph(nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2) 973 end do 974 975 ! This to avoid spurious results for the acoustic modes. 976 ! In principle, gamma(q) --> 0 and lambda(q) --> 0 for q --> 0 977 ! but the Fourier interpolated gammas do not fulfill this property so we set everything 978 ! to zero when we are inside a sphere or radius 979 if (normv(qpt, cryst%gmet, "G") < EPH_Q0TOL) then 980 !write(std_out,*)"Setting values to zero." 981 gamma_ph(1:3) = zero 982 lambda_ph(1:3) = zero 983 end if 984 985 end subroutine phgamma_interp
m_phgamma/phgamma_interp_setup [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
phgamma_interp_setup
FUNCTION
This routines prepares the internal tables used to interpolate the linewidths in q-space
INPUTS
SOURCE
1001 subroutine phgamma_interp_setup(gams, cryst) 1002 1003 !Arguments ------------------------------------ 1004 !scalars 1005 class(phgamma_t),intent(inout) :: gams 1006 type(crystal_t),intent(in) :: cryst 1007 1008 !Local variables------------------------------- 1009 !scalars 1010 integer,parameter :: qtor1 = 1 1011 integer :: iq_bz,iq_ibz,spin,ierr,ii 1012 !character(len=500) :: msg 1013 !arrays 1014 integer,allocatable :: qirredtofull(:),qpttoqpt(:,:,:) 1015 real(dp),allocatable :: coskr(:,:),sinkr(:,:) 1016 real(dp),allocatable :: gamma_qpt(:,:,:,:),atmfrc(:,:) 1017 real(dp),allocatable :: vals_bz(:,:,:,:) 1018 1019 ! ************************************************************************* 1020 1021 ABI_MALLOC_OR_DIE(vals_bz, (2, gams%natom3**2, gams%nqbz, gams%nsppol), ierr) 1022 vals_bz = zero 1023 1024 ! Build tables needed by complete_gamma. 1025 call ephtk_mkqtabs(cryst, gams%nqibz, gams%qibz, gams%nqbz, gams%qbz, qirredtofull, qpttoqpt) 1026 1027 ! Fill BZ array with IBZ data. 1028 do spin=1,gams%nsppol 1029 do iq_ibz=1,gams%nqibz 1030 iq_bz = qirredtofull(iq_ibz) 1031 vals_bz(:,:,iq_bz,spin) = reshape(gams%vals_qibz(:,:,:,iq_ibz,spin), [2, gams%natom3**2]) 1032 end do 1033 end do 1034 1035 ! Complete vals_bz in the full BZ. 1036 ! FIXME: Change complete_gamma API to pass (..., nsppol) 1037 ABI_MALLOC(gamma_qpt, (2, gams%natom3**2, gams%nsppol, gams%nqbz)) 1038 do spin=1,gams%nsppol 1039 gamma_qpt(:, :, spin, :) = vals_bz(:, :, :, spin) 1040 end do 1041 1042 call complete_gamma(cryst, gams%natom3, gams%nsppol, gams%nqibz, gams%nqbz, & 1043 gams%eph_scalprod, qirredtofull, qpttoqpt, gamma_qpt) 1044 1045 do spin=1,gams%nsppol 1046 vals_bz(:, :, :, spin) = gamma_qpt(:, :, spin, :) 1047 end do 1048 1049 ABI_FREE(gamma_qpt) 1050 ABI_FREE(qirredtofull) 1051 ABI_FREE(qpttoqpt) 1052 1053 ! This call is not executed in elphon! 1054 if (gams%symgamma == 1) then 1055 do spin=1,gams%nsppol 1056 do iq_bz=1,gams%nqbz 1057 call tgamma_symm(cryst, gams%qbz(:,iq_bz), vals_bz(:,:,iq_bz,spin)) 1058 end do 1059 end do 1060 end if 1061 1062 ! Now FT to real space too 1063 ! NOTE: gprim (not gprimd) is used for all FT interpolations, 1064 ! to be consistent with the dimensions of the rpt, which come from anaddb. 1065 ABI_MALLOC_OR_DIE(gams%vals_rpt, (2, gams%natom3**2, gams%nrpt, gams%nsppol), ierr) 1066 gams%vals_rpt = zero 1067 1068 ! q --> r 1069 ABI_MALLOC(coskr, (gams%nqbz, gams%nrpt)) 1070 ABI_MALLOC(sinkr, (gams%nqbz, gams%nrpt)) 1071 call ftgam_init(gams%gprim, gams%nqbz, gams%nrpt, gams%qbz, gams%rpt, coskr, sinkr) 1072 1073 do spin=1,gams%nsppol 1074 call ftgam(gams%wghatm, vals_bz(:,:,:,spin), gams%vals_rpt(:,:,:,spin), gams%natom, gams%nqbz,& 1075 gams%nrpt, qtor1, coskr, sinkr) 1076 1077 ! Enforce "acoustic" rule on vals_rpt 1078 ! This call is not executed in elphon! 1079 if (gams%asr /= 0) then 1080 ABI_MALLOC(atmfrc, (3*gams%natom*3*gams%natom, gams%nrpt)) 1081 do ii=1,2 1082 atmfrc = gams%vals_rpt(ii,:,:,spin) 1083 !gals%vals_rpt(2,:,,spin) = zero 1084 call asrif9(gams%asr, atmfrc, gams%natom, gams%nrpt, gams%rpt, gams%wghatm) 1085 gams%vals_rpt(ii,:,:,spin) = atmfrc 1086 end do 1087 ABI_FREE(atmfrc) 1088 end if 1089 end do 1090 1091 ABI_FREE(vals_bz) 1092 ABI_FREE(coskr) 1093 ABI_FREE(sinkr) 1094 1095 end subroutine phgamma_interp_setup
m_phgamma/phgamma_linwid [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
phgamma_linwid
FUNCTION
Interpolate the phonon linewidths along an arbitrary q-path (use Fourier interpolation).
INPUTS
cryst<crystal_t>=Info on the unit cell and symmetries. ifc<ifc_type>=Interatomic force constants. ndivsm=Number of points used to sample the smallest segment nvert = Number of extrema in qverts qverts(3,nvert) = vertices of reciprocal space trajectory basename=name used to create the different output files (text format). ncid=Netcdf file handler (already open in the caller). comm=MPI communicator
OUTPUT
wminmax=Minimum and max phonon frequency obtained on the path (Hartree units)
SOURCE
1457 subroutine phgamma_linwid(gams, cryst, ifc, ndivsm, nvert, qverts, basename, ncid, wminmax, comm) 1458 1459 !Arguments ------------------------------------ 1460 !scalars 1461 integer,intent(in) :: nvert,ndivsm,comm,ncid 1462 type(crystal_t),intent(in) :: cryst 1463 type(ifc_type),intent(in) :: ifc 1464 class(phgamma_t),intent(inout) :: gams 1465 character(len=*),intent(in) :: basename 1466 !arrays 1467 real(dp),intent(in) :: qverts(3,nvert) 1468 real(dp),intent(out) :: wminmax(2) 1469 1470 !Local variables------------------------------- 1471 !scalars 1472 integer,parameter :: master = 0 1473 integer :: natom,ii,mu,iqpt,natom3,nsppol,ierr 1474 integer :: spin,unt,nqpt,nrpt,cnt,nproc,my_rank 1475 integer :: ncerr 1476 real(dp) :: omega_min,omega_max,wtmp,omega 1477 character(len=500) :: msg 1478 type(kpath_t) :: qpath 1479 !arrays 1480 real(dp) :: gamma_spin(gams%nsppol),lambda_spin(gams%nsppol) 1481 real(dp) :: displ_cart(2,3*cryst%natom,3*cryst%natom) 1482 real(dp) :: phfrq(3*cryst%natom),gamma_ph(3*cryst%natom),lambda_ph(3*cryst%natom) 1483 real(dp) :: qpt(3),shift(3) 1484 real(dp),allocatable :: all_phfreq(:,:),all_gammaq(:,:,:),all_lambdaq(:,:,:),all_displ_cart(:,:,:,:) 1485 1486 ! ********************************************************************* 1487 1488 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 1489 natom = cryst%natom; natom3 = gams%natom3; nsppol = gams%nsppol; nrpt = gams%nrpt 1490 1491 ! Define the q-path along which phonon linwid will be interpolated. 1492 qpath = kpath_new(qverts, cryst%gprimd, ndivsm) 1493 nqpt = qpath%npts 1494 1495 ! Allocate workspace arrays for MPI. 1496 ABI_CALLOC(all_phfreq, (natom3, nqpt)) 1497 ABI_CALLOC(all_gammaq, (natom3, nqpt, nsppol)) 1498 ABI_CALLOC(all_lambdaq, (natom3, nqpt, nsppol)) 1499 ABI_CALLOC(all_displ_cart, (2, natom3, natom3, nqpt)) 1500 1501 ! initialize the minimum and maximum phonon frequency 1502 omega_min = huge(one); omega_max = -huge(one) 1503 1504 ! Interpolation along specified path in q space (keep spin dep.) 1505 cnt = 0 1506 do spin=1,nsppol 1507 do iqpt=1,nqpt 1508 cnt = cnt + 1; if (mod(cnt, nproc) /= my_rank) cycle 1509 call wrap2_pmhalf(qpath%points(:,iqpt), qpt, shift) 1510 1511 ! Get phgamma 1512 call gams%interp(cryst, ifc, spin, qpt, phfrq, gamma_ph, lambda_ph, displ_cart) 1513 all_gammaq(:, iqpt, spin) = gamma_ph 1514 all_lambdaq(:, iqpt, spin) = lambda_ph 1515 if (spin == 1) then 1516 all_phfreq(:, iqpt) = phfrq 1517 all_displ_cart(:, :, :, iqpt) = displ_cart 1518 end if 1519 1520 ! Find max/min phonon frequency along path chosen 1521 ! presumed to be representative of full BZ to within 10 percent 1522 omega_min = min(omega_min, phfrq(1)) 1523 omega_max = max(omega_max, phfrq(natom3)) 1524 end do ! end iqpt do 1525 end do ! spin 1526 1527 ! Collect results 1528 if (omega_min > tol12) omega_min = zero 1529 wtmp = omega_min; call xmpi_min(wtmp, omega_min, comm, ierr) 1530 wtmp = omega_max; call xmpi_max(wtmp, omega_max, comm, ierr) 1531 wminmax = [omega_min, omega_max] 1532 1533 call xmpi_sum_master(all_gammaq, master, comm, ierr) 1534 call xmpi_sum_master(all_lambdaq, master, comm, ierr) 1535 call xmpi_sum_master(all_phfreq, master, comm, ierr) 1536 call xmpi_sum_master(all_displ_cart, master, comm, ierr) 1537 1538 ! Master writes text file with final results. 1539 ! 3 * natom blocks, one block for each phonon mode. 1540 ! Each block contains: 1541 ! iqpt omega_(q) gamma(q) lambda(q) nesting(q) .... 1542 1543 if (xmpi_comm_rank(comm) == master) then 1544 if (open_file(strcat(basename, '_PHGAMMA'), msg, newunit=unt, action="write", form="formatted", status="unknown") /= 0) then 1545 ABI_ERROR(msg) 1546 end if 1547 1548 write(unt,'(a)') '#' 1549 write(unt,'(a)') '# ABINIT package: E-PH band structure file. Hartree units' 1550 write(unt,'(a)') '#' 1551 write(unt,'(a,i0,a)')'# Phonon frequencies, ph linewidths and lambda calculated on ',nqpt,' q-points' 1552 call qpath%print(header="Description of the q-path:", unit=unt, pre="#") 1553 do ii=1,2; write(unt,'(a)') "# "; end do 1554 1555 write(unt,'(a,e16.6)')"# Total DOS at Fermi level ",sum(gams%n0) 1556 do spin=1,nsppol 1557 write(unt,"(a,i0,a,e16.6)")"# The DOS at Fermi level for spin ",spin," is ",gams%n0(spin) 1558 end do 1559 1560 do mu=1,natom3 1561 write(unt,'(a)')"#" 1562 if (nsppol == 1) write(unt,'(a,i0,a)')"# phonon mode [",mu,"] q-index omega gamma lambda" 1563 if (nsppol == 2) write(unt,'(a,i0,a)')& 1564 "# phonon mode [",mu,"] q-index omega gamma_tot lambda_tot gamma[spin=1] lambda[spin=1] gamma[2] lambda[2]" 1565 write(unt,'(a)')"#" 1566 do iqpt=1,nqpt 1567 omega = all_phfreq(mu, iqpt) 1568 gamma_spin = all_gammaq(mu, iqpt, :) 1569 lambda_spin = all_lambdaq(mu, iqpt, :) 1570 if (nsppol == 1) then 1571 write(unt,'(i8,3es16.6)' )iqpt,omega,gamma_spin(1),lambda_spin(1) 1572 else 1573 write(unt,'(i8,es20.6,6es16.6)' )iqpt,omega,& 1574 sum(gamma_spin),sum(lambda_spin),& 1575 gamma_spin(1),lambda_spin(1),& 1576 gamma_spin(2),lambda_spin(2) 1577 end if 1578 end do 1579 end do 1580 1581 close(unt) 1582 1583 ! Write data to netcdf file 1584 if (ncid /= nctk_noid) then 1585 ncerr = nctk_def_dims(ncid, [& 1586 nctkdim_t("natom3", 3*natom), nctkdim_t("nqpath", nqpt), nctkdim_t("number_of_spins", nsppol) & 1587 ], defmode=.True.) 1588 NCF_CHECK(ncerr) 1589 1590 ncerr = nctk_def_arrays(ncid, [& 1591 nctkarr_t('qpath', "dp", "number_of_reduced_dimensions, nqpath"), & 1592 nctkarr_t('phfreq_qpath', "dp", "natom3, nqpath"), & 1593 nctkarr_t('phdispl_cart_qpath', "dp", "two, natom3, natom3, nqpath"), & 1594 nctkarr_t('phgamma_qpath', "dp", "natom3, nqpath, number_of_spins"), & 1595 nctkarr_t('phlambda_qpath', "dp", "natom3, nqpath, number_of_spins")]) 1596 NCF_CHECK(ncerr) 1597 1598 NCF_CHECK(nctk_set_datamode(ncid)) 1599 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "qpath"), qpath%points)) 1600 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phfreq_qpath"), all_phfreq)) 1601 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phdispl_cart_qpath"), all_displ_cart)) 1602 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phgamma_qpath"), all_gammaq)) 1603 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phlambda_qpath"), all_lambdaq)) 1604 end if 1605 end if ! master 1606 1607 ABI_FREE(all_phfreq) 1608 ABI_FREE(all_gammaq) 1609 ABI_FREE(all_lambdaq) 1610 ABI_FREE(all_displ_cart) 1611 1612 call qpath%free() 1613 1614 end subroutine phgamma_linwid
m_phgamma/phgamma_ncwrite [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
phgamma_ncwrite
FUNCTION
Write the results stored in the phgamma_t datatype to netcdf file.
INPUTS
cryst<crystal_t>=Crystalline structure. ifc<ifc_type>=Interatomic force constants. ncid=Netcdf file handler (already open in the caller).
OUTPUT
SOURCE
635 subroutine phgamma_ncwrite(gams, cryst, ifc, ncid) 636 637 !Arguments ------------------------------------ 638 !scalars 639 type(phgamma_t),intent(inout) :: gams 640 type(crystal_t),intent(in) :: cryst 641 type(ifc_type),intent(in) :: ifc 642 integer,intent(in) :: ncid 643 644 !Local variables------------------------------- 645 !scalars 646 integer :: iq_ibz,spin,mu 647 real(dp) :: lambda_tot 648 character(len=500) :: msg 649 !arrays 650 real(dp) :: phfrq(3*cryst%natom), gamma_ph(3*cryst%natom), lambda_ph(3*cryst%natom) 651 real(dp) :: displ_cart(2,3*cryst%natom,3*cryst%natom) 652 653 ! ************************************************************************* 654 655 ! Write data to files for each q point, also compute total lambda. 656 lambda_tot = zero 657 do spin=1,gams%nsppol 658 do iq_ibz=1,gams%nqibz 659 660 ! Get phonon frequencies, gamma(q,nu) and lambda(q,nu) 661 ! Quantities are already summed over (collinear) spin channels if nsppol == 1 662 call gams%eval_qibz(cryst, ifc, iq_ibz, spin, phfrq, gamma_ph, lambda_ph, displ_cart) 663 664 do mu=1,gams%natom3 665 lambda_tot = lambda_tot + lambda_ph(mu) * gams%wtq(iq_ibz) 666 end do 667 668 ! Write data to netcdf file 669 if (ncid /= nctk_noid) then 670 if (spin == 1) then 671 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "phfreq_qibz"), phfrq, start=[1, iq_ibz])) 672 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phdispl_cart_qibz'), displ_cart, start=[1, 1, 1, iq_ibz])) 673 end if 674 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phgamma_qibz'), gamma_ph, start=[1, iq_ibz, spin])) 675 NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, 'phlambda_qibz'), lambda_ph, start=[1, iq_ibz, spin])) 676 end if 677 678 ! Output to the main output file 679 if (gams%nsppol == 2) then 680 write(msg,'(2a,3es16.6,a,i1,a,a)')ch10,& 681 ' q-point =',gams%qibz(:, iq_ibz),' spin = ',spin,ch10,& 682 ' Mode number Frequency (Ha) Linewidth (Ha) Lambda(q,n)' 683 else 684 write(msg,'(2a,3es16.6,a,a)')ch10,& 685 ' q-point =',gams%qibz(:, iq_ibz),ch10,& 686 ' Mode number Frequency (Ha) Linewidth (Ha) Lambda(q,n)' 687 end if 688 call wrtout([std_out, ab_out], msg) 689 690 do mu=1,gams%natom3 691 write(msg,'(i5,es20.6,2es16.6)')mu, phfrq(mu), gamma_ph(mu), lambda_ph(mu) 692 call wrtout([std_out, ab_out], msg) 693 end do 694 695 end do 696 ! Add blank lines to output files between spins 697 call wrtout([std_out, ab_out], "") 698 end do 699 700 write(ab_out,"(a,f8.4)")" lambda= ",lambda_tot 701 !write(ab_out,"(a,f8.4)")" omega_log= ",omega_log 702 703 end subroutine phgamma_ncwrite
m_phgamma/phgamma_setup_qpoint [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
phgamma_setup_qpoint
FUNCTION
SOURCE
4338 subroutine phgamma_setup_qpoint(gams, fs, cryst, ebands, spin, ltetra, qpt, nesting, comm) 4339 4340 !Arguments ------------------------------------ 4341 type(phgamma_t),intent(inout) :: gams 4342 type(fstab_t),intent(inout) :: fs 4343 type(crystal_t),intent(in) :: cryst 4344 type(ebands_t),intent(in) :: ebands 4345 integer,intent(in) :: spin, ltetra, comm 4346 integer,intent(out) :: nesting 4347 !arrays 4348 real(dp),intent(in) :: qpt(3) 4349 4350 !Local variables------------------------------- 4351 !scalars 4352 integer,parameter :: enough = 5 4353 integer :: nkbz, ierr, nb, ik_bz, ik_ibz, ikq_ibz, ikq_fs, ik_fs, i1, i2, i3, nkfs_q, nene 4354 integer :: ib1, ib2 ! band_k, band_kq 4355 real(dp),parameter :: max_occ1 = one 4356 real(dp) :: cpu, wall, gflops, enemin, enemax 4357 character(len=500) :: msg 4358 type(krank_t) :: ibz_krank 4359 type(t_tetrahedron) :: tetra 4360 !type(lgroup_t) :: lgq 4361 character(len=80) :: errorstring 4362 !arrays 4363 integer :: nge(3), ngw(3), g0bz_kq(3) 4364 integer,allocatable :: select_ikfs(:), indkpt(:), kbz2fs(:) !, symrecfm(:,:,:) 4365 real(dp) :: kk(3), kq(3) 4366 real(dp),allocatable :: eig_k(:,:), eig_kq(:,:), wght_bz(:,:,:), kbz(:,:) 4367 real(dp),allocatable :: work_k(:), work_kq(:), dtweightde(:,:,:), tweight(:,:,:) 4368 4369 ! ************************************************************************* 4370 4371 call cwtime(cpu, wall, gflops, "start") 4372 4373 ! The double delta with tetra is ill-defined for q == 0. 4374 ! Set nesting to 1 and return 4375 nesting = 0 4376 if (abs(fs%eph_intmeth) == 2 .and. all(abs(qpt) < tol12)) then 4377 ABI_COMMENT("Tetrahedron for double grid with q = 0 is ill-defined. Using adaptive gaussian.") 4378 nesting = 1 4379 end if 4380 4381 ! Compute little group of the q-point. Map fs%kpts to ebands%kptns (IBZ) 4382 !lg_q = lgroup_new(cryst, qpt, self%timrev, fs%nkfs, fs%kpts, ebands%nkpt, ebands%kpnts, comm) 4383 !do ik_fs=1,fs%nkfs 4384 ! ik_lgibz = lg_q%bz2ibz_smap(1, ik_fs) 4385 ! lg_q%weights(iklg_ibz) 4386 !end do 4387 !call lg_q%free() 4388 4389 if (fs%eph_intmeth == 1 .or. nesting == 1) then 4390 ! Gaussian method: 4391 ! Distribute k-points within the FS window inside comm. 4392 ! 1) Select k-points such that k+q is stil inside the FS window 4393 ! 2) Distribute effective k-points assuming all procs in comm have all FS k-points (no filtering) 4394 ABI_MALLOC(select_ikfs, (fs%nkfs)) 4395 nkfs_q = 0 4396 do ik_fs=1,fs%nkfs 4397 kq = fs%kpts(:, ik_fs) + qpt 4398 ikq_fs = fs%findkg0(kq, g0bz_kq); if (ikq_fs == -1) cycle 4399 nkfs_q = nkfs_q + 1 4400 select_ikfs(nkfs_q) = ik_fs 4401 end do 4402 ABI_SFREE(gams%my_ifsk_q) 4403 call xmpi_split_list(nkfs_q, select_ikfs, comm, gams%my_nfsk_q, gams%my_ifsk_q) 4404 ABI_FREE(select_ikfs) 4405 call cwtime_report(" phgamma_setup_qpoint", cpu, wall, gflops) 4406 return 4407 end if 4408 4409 ! Tetrahedron method: 4410 ! 4411 ! 1) Compute weights for double delta integration. 4412 ! 2) Filter k-points according to the weights (k must be in FS window with non-zero weight) 4413 ! 3) Distribute effective k-points assuming all procs have all FS k-points. 4414 4415 fs%dbldelta_tetra_weights_kfs = zero 4416 4417 ABI_CHECK(isdiagmat(ebands%kptrlatt), "kptrlatt must be diagonal when tetra is used.") 4418 ABI_CHECK(ebands%nshiftk == 1, "nshiftk must be 1 when tetra is used") 4419 nge = get_diag(ebands%kptrlatt); ngw = nge 4420 nkbz = product(nge(1:3)) 4421 4422 ! TODO: Handle symmetries in a cleaner way. Change API of krank_new to pass symafm and kptopt 4423 ibz_krank = krank_new(ebands%nkpt, ebands%kptns, nsym=cryst%nsym, symrec=cryst%symrec, & 4424 time_reversal=kpts_timrev_from_kptopt(ebands%kptopt) == 1) 4425 4426 ! Compute eig_k and eig_kq in full BZ for the relevant bands around Ef. 4427 nb = fs%maxnb 4428 ABI_MALLOC(kbz, (3, nkbz)) 4429 ABI_MALLOC(indkpt, (nkbz)) 4430 ABI_MALLOC(eig_k, (nb, nkbz)) 4431 ABI_MALLOC(eig_kq, (nb, nkbz)) 4432 ABI_MALLOC(kbz2fs, (nkbz)) 4433 4434 kbz2fs = -1; ierr = 0; ik_bz = 0 4435 do i3=0,nge(3) - 1 4436 do i2=0,nge(2) - 1 4437 do i1=0,nge(1) - 1 4438 ik_bz = ik_bz + 1 4439 indkpt(ik_bz) = ik_bz 4440 4441 ! Find correspondence between libtetra mesh and the IBZ. 4442 kk = ([i1, i2, i3] + ebands%shiftk(:, 1)) / nge(:) 4443 kbz(:, ik_bz) = kk 4444 4445 ik_fs = fs%krank%get_index(kk) 4446 if (ik_fs /= -1) then 4447 kbz2fs(ik_bz) = ik_fs 4448 else 4449 !ABI_ERROR(sjoin('kpt:', trim(ktoa(kk)), 'is not in FS!!')) 4450 end if 4451 4452 ik_ibz = ibz_krank%get_index(kk) 4453 if (ik_ibz < 1) then 4454 if (ierr <= enough) then 4455 ABI_WARNING(sjoin('kpt:', trim(ktoa(kk)), 'has no symmetric among the k-points!')) 4456 end if 4457 ierr = ierr + 1; cycle 4458 end if 4459 4460 eig_k(:, ik_bz) = ebands%eig(fs%bmin:fs%bmax, ik_ibz, spin) 4461 4462 ! Find correspondence between the k+q in the BZ grid and the IBZ. 4463 kq = kk + qpt 4464 ikq_ibz = ibz_krank%get_index(kq) 4465 4466 if (ikq_ibz < 1) then 4467 if (ierr <= enough) then 4468 ABI_WARNING(sjoin('kpt + qpt:', trim(ktoa(kq)), 'has no symmetric among the k-points!')) 4469 end if 4470 ierr = ierr + 1; cycle 4471 end if 4472 4473 eig_kq(:, ik_bz) = ebands%eig(fs%bmin:fs%bmax, ikq_ibz, spin) 4474 end do 4475 end do 4476 end do 4477 4478 ABI_CHECK(ierr == 0, "See above warnings") 4479 call ibz_krank%free() 4480 4481 select case (ltetra) 4482 case (1, 2) 4483 ! Use libtetra routines. 4484 ! Compute weights for double delta integration. Note that libtetra assumes Ef set to zero. 4485 ! TODO: Average weights over degenerate states? 4486 if (ltetra == 1) call wrtout(std_out, " Using linear tetrahedron method from libtetrabz (ltetra 1)") 4487 if (ltetra == 2) call wrtout(std_out, " Using optimized tetrahedron method from libtetrabz (ltetra 2)") 4488 eig_k = eig_k - ebands%fermie; eig_kq = eig_kq - ebands%fermie 4489 ABI_MALLOC(wght_bz, (nb, nb, nkbz)) 4490 call libtetrabz_dbldelta(ltetra, cryst%gprimd, nb, nge, eig_k, eig_kq, ngw, wght_bz, comm=comm) 4491 ! Revert changes to eig arrays 4492 eig_k = eig_k + ebands%fermie !; eig_kq = eig_kq + ebands%fermie 4493 4494 ! Reindex from full BZ to fs% kpoints. 4495 do ik_bz=1,nkbz 4496 ik_fs = kbz2fs(ik_bz) 4497 if (ik_fs /= -1) then 4498 fs%dbldelta_tetra_weights_kfs(:,:,ik_fs) = wght_bz(:,:,ik_bz) 4499 else 4500 !write(std_out,*)"should be zero :", wght_bz(:,:,ik_bz) 4501 end if 4502 end do 4503 ABI_FREE(wght_bz) 4504 4505 case (3) 4506 ! Tetrahedron method with Allen's approach for double delta. 4507 ! WARNING: wrong results maybe bug somewhere. 4508 call wrtout(std_out, " Calling Allen's version for q-point") 4509 nene = 3 4510 enemin = ebands%fermie - tol6 4511 enemax = ebands%fermie + tol6 4512 4513 ABI_MALLOC(dtweightde, (nkbz, nene, nene)) 4514 ABI_MALLOC(tweight, (nkbz, nene, nene)) 4515 ABI_MALLOC(work_k, (nkbz)) 4516 ABI_MALLOC(work_kq, (nkbz)) 4517 4518 call init_tetra(indkpt, cryst%gprimd, fs%klatt, kbz, nkbz, tetra, ierr, errorstring, comm) 4519 ABI_CHECK(ierr == 0, "init_tetra returned ierr /= 0") 4520 4521 fs%dbldelta_tetra_weights_kfs = zero 4522 do ib2=1,nb 4523 work_k = eig_k(ib2, :) 4524 do ib1=1,nb 4525 work_kq = eig_kq(ib1, :) 4526 ! TODO: Average weights over degenerate states? 4527 call get_dbl_tetra_weight(work_k, work_kq, enemin, enemax, enemin, enemax, & 4528 max_occ1, nene, nene, nkbz, tetra, tweight, dtweightde, ierr) 4529 ABI_CHECK(ierr == 0, "get_dbldelta_weights returned ierr /= 0") 4530 4531 ! Reindex from full BZ to fs% kpoints. 4532 do ik_bz=1,nkbz 4533 ik_fs = kbz2fs(ik_bz) 4534 if (ik_fs /= -1) then 4535 fs%dbldelta_tetra_weights_kfs(ib1,ib2,ik_fs) = dtweightde(ik_bz, 2, 2) ! / nkbz 4536 else 4537 !write(std_out,*)"should be zero :", wght_bz(:,:,ik_bz) 4538 end if 4539 end do 4540 end do 4541 end do 4542 4543 call destroy_tetra(tetra) 4544 ABI_FREE(work_k) 4545 ABI_FREE(work_kq) 4546 ABI_FREE(dtweightde) 4547 ABI_FREE(tweight) 4548 4549 case default 4550 ABI_ERROR(sjoin("Invalid value of ltetra:", itoa(ltetra))) 4551 end select 4552 4553 ! Now we can filter the k-points according to the tetra weights and distribute them inside comm. 4554 ! Assuming all procs in comm have all k-points in the IBZ. 4555 ! nkfs_q is the total number of BZ k-points on the FS contributing for this q-point. 4556 ABI_MALLOC(select_ikfs, (fs%nkfs)) 4557 nkfs_q = 0 4558 do ik_fs=1,fs%nkfs 4559 if (any(abs(fs%dbldelta_tetra_weights_kfs(:,:,ik_fs)) > zero)) then 4560 nkfs_q = nkfs_q + 1 4561 select_ikfs(nkfs_q) = ik_fs 4562 end if 4563 end do 4564 4565 ! Here we compute my_ifsk_q. 4566 ABI_SFREE(gams%my_ifsk_q) 4567 call xmpi_split_list(nkfs_q, select_ikfs, comm, gams%my_nfsk_q, gams%my_ifsk_q) 4568 ABI_FREE(select_ikfs) 4569 4570 !write(std_out,"(2(a,i0),/)")" Treating ", gams%my_nfsk_q, " my k-points in the FS window over total nkfs: ", fs%nkfs 4571 4572 write(msg, "(2(a,i0),a)") & 4573 " Number of k-points in the FS window treated by this MPI proc: ", gams%my_nfsk_q, " over: ", nkfs_q, ch10 4574 !" Number of MPI procs in kpt_comm: ", gams%kpt_comm%nproc 4575 call wrtout(std_out, msg) 4576 4577 ABI_FREE(kbz2fs) 4578 ABI_FREE(indkpt) 4579 ABI_FREE(kbz) 4580 ABI_FREE(eig_k) 4581 ABI_FREE(eig_kq) 4582 4583 call cwtime_report(" phgamma_setup_qpoint", cpu, wall, gflops) 4584 4585 end subroutine phgamma_setup_qpoint
m_phgamma/phgamma_t [ Types ]
[ Top ] [ m_phgamma ] [ Types ]
NAME
phgamma_t
FUNCTION
Provides methods for computing phonon linewidths, interpolating the results in q-space and evaluate superconducting properties withing the isotropic formalism.
SOURCE
106 type,public :: phgamma_t 107 108 integer :: natom 109 ! Number of atoms per unit cell. 110 111 integer :: natom3 112 ! Number of phonon branches i.e. 3*natom. 113 114 integer :: nsppol 115 ! Number of independent spin polarizations. 116 117 integer :: nspinor 118 ! Number of spinorial components. 119 120 integer :: nqibz 121 ! Number of q-points in the IBZ. 122 123 integer :: my_nqibz 124 ! Number of q-points in the IBZ treated by the current MPI processor 125 126 integer :: nqbz 127 ! Number of q-points in the BZ. 128 129 integer :: eph_scalprod = 0 130 ! This to call anaddb routines. Note that eph_scalprod 1 is not supported in eph. 131 132 integer :: bcorr = 0 133 ! 1 to include Blochl correction in the tetrahedron method else 0. 134 135 integer :: prteliash = 0 136 ! This flag activates the computation of the Eliashberg function. 137 138 integer :: nrpt 139 ! Number of R-points in the real space representation of the gamma matrices. 140 141 integer :: symgamma 142 ! 1 if gamma matrices should be symmetrized by symdyma when using Fourier interpolation 143 144 integer :: asr 145 ! If the "Acoustic rule" at Gamma should be enforced. 146 147 integer :: ndir_transp 148 ! 0 if no transport, otherwise 3 149 150 integer :: ngqpt(3) 151 ! Number of divisions in the q-mesh. 152 153 integer :: nene 154 ! Number of chemical potential values used for inelastic integration 155 156 !integer :: my_nqpt 157 !integer,allocatable :: my_iqpt(:) 158 159 integer :: my_nfsk_q 160 ! Number of k-points in the FS treated by this MPI processor for a given q. 161 ! Computed in phgamma_setup_qpoint 162 163 integer,allocatable :: my_ifsk_q(:) 164 ! Index of the FS k-points treated by this processor for a given q 165 ! Computed in phgamma_setup_qpoint 166 167 integer :: my_nspins 168 ! Number of spins treated by the MPI rank 169 170 integer,allocatable :: my_spins(:) 171 ! my_spins(my_nspins) 172 ! Indirect table giving the spin indices treated by this rank. 173 ! Used only the collinear case with nspinor == 1 174 175 !integer :: my_npert 176 ! Number of atomic perturbations or phonon modes treated by this MPI rank. 177 178 !integer(i1b),allocatable :: itreat_qibz(:) 179 ! itreat_qibz(nqibz) 180 ! Table used to distribute potentials over q-points in the IBZ. 181 ! The loop over qpts in the IBZ(k) is MPI distributed inside qpt_comm accordinging to this table. 182 ! 0 if this IBZ point is not treated by this proc. 183 ! 1 if this IBZ is treated. 184 185 !integer,allocatable :: my_pinfo(:,:) 186 ! my_pinfo(3, my_npert) 187 ! my_pinfo(1, ip) gives the `idir` index of the ip-th perturbation. 188 ! my_pinfo(2, ip) gives the `ipert` index of the ip-th perturbation. 189 ! my_pinfo(3, ip) gives `pertcase`=idir + (ipert-1)*3 190 191 !integer,allocatable :: pert_table(:,:) 192 ! pert_table(2, natom3) 193 ! pert_table(1, npert): rank of the processor treating this atomic perturbation. 194 ! pert_table(2, npert): imyp index in my_pinfo table, -1 if this rank is not treating ipert. 195 196 integer, allocatable :: my_iqibz(:) 197 ! indices of ibz iq in local array. -1 if iq does not belong to current proc 198 199 integer, allocatable :: my_iqbz(:) 200 ! List of iq_ibz indices treated by current proc 201 202 real(dp) :: enemin 203 ! Minimal chemical potential value used for inelastic integration Copied from fstab 204 205 real(dp) :: deltaene 206 ! Chemical potential increment for inelastic integration Copied from fstab 207 ! for simplicity could be made equal to phonon frequency step 208 209 real(dp) :: gprim(3,3) 210 ! Needed for Fourier interpolation. 211 ! NOTE: gprim (not gprimd) is used for all FT interpolations, 212 ! to be consistent with the dimensions of the rpt, which come from anaddb. 213 214 real(dp),allocatable :: n0(:) 215 ! (%nsppol) 216 ! Density of states at the Fermi level per spin in a.u. 217 218 real(dp),allocatable :: qibz(:,:) 219 ! qibz(3,nqibz) 220 ! Reduced coordinates of the q-points in the IBZ. 221 222 real(dp),allocatable :: wtq(:) 223 ! wtq(nqibz) 224 ! Weights of the q-points in the IBZ (normalized to one) 225 226 real(dp),allocatable :: qbz(:,:) 227 ! qbz(3,nqbz) 228 ! Reduced coordinates of the q-points in the BZ. 229 230 real(dp),allocatable :: rpt(:,:) 231 ! rpt(3,nrpt) 232 ! Reduced coordinates ***in terms of rprim*** of the lattice points used 233 ! for the Fourier transform of the phonon linewidths. 234 235 real(dp),allocatable :: wghatm(:,:,:) 236 ! wghatm(natom,natom,nrpt) 237 ! Weights used in the FT of the phonon linewidths. 238 239 real(dp),allocatable :: vals_qibz(:,:,:,:,:) 240 ! vals_qibz(2,natom3,natom3,nqibz,nsppol)) in reduced coordinates for each q-point in the IBZ. 241 ! vals_qibz {\tau'\alpha',\tau\alpha} = sum over k, ib_kq, ib_k 242 ! <psi_{k+q,ib_kq} | H(1)_{\tau'\alpha'} | psi_{k,ib_k}>* \cdot 243 ! <psi_{k+q,ib_kq} | H(1)_{\tau \alpha } | psi_{k,ib_k}> 244 245 !NOTE: choice to put nsppol before or after nqbz is a bit arbitrary 246 ! abinit uses nband,nkpt,nsppol, but here for convenience nkpt_phon,nsppol,nqbz interpolation is on qpt 247 !MG: I think that nsppol should be the last dimension set call to ftgam. 248 249 real(dp),allocatable :: vals_rpt(:,:,:,:) 250 ! vals_rpt(2,natom3**2,nrpt,nsppol) 251 ! tgamma matrices in real space in reduced coordinates. Used for the Fourier interpolation. 252 253 ! transport stuff with velocity factors 254 real(dp),allocatable :: vals_in_qibz(:,:,:,:,:,:) 255 real(dp),allocatable :: vals_out_qibz(:,:,:,:,:,:) 256 ! vals_XX_qibz(2,9,natom3,natom3,nqibz,nsppol)) in reduced coordinates for each q-point in the IBZ. 257 258 real(dp),allocatable :: vals_in_rpt(:,:,:,:,:) 259 real(dp),allocatable :: vals_out_rpt(:,:,:,:,:) 260 ! vals_XX_rpt(2,9,natom3**2,nrpt,nsppol) 261 ! tgamma matrices in real space in reduced coordinates. 262 263 ! gamma matrices keeping full electron energy dependency 264 real(dp),allocatable :: vals_ee(:,:,:,:,:,:,:) 265 ! vals_eew(2, nene, nene, natom3, natom3, nqibz, nsppol) 266 267 contains 268 269 procedure :: free => phgamma_free 270 ! Free memory. 271 272 procedure :: interp => phgamma_interp 273 ! Interpolates the phonon linewidths. 274 275 procedure :: eval_qibz => phgamma_eval_qibz 276 ! Evaluate phonon linewidths without Fourier interpolation. 277 278 procedure :: interp_setup => phgamma_interp_setup 279 ! Compute the tables used for the interpolation in q-space. 280 281 procedure :: linwid => phgamma_linwid 282 ! Interpolate linewidths along an arbitrary q-path. 283 284 end type phgamma_t 285 286 public :: phgamma_init ! Creation method.
m_phgamma/phgamma_vv_eval_qibz [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
phgamma_vv_eval_qibz
FUNCTION
Compute the phonon linewidths times velocity squared, for q-points in the IBZ without performing interpolation.
INPUTS
gams<phgamma_t> cryst<crystal_t>=Crystal structure. ifc<ifc_type>=Interatomic force constants. iq_ibz=Index of the q-point in the IBZ array. spin=Spin index
OUTPUT
phfrq(gams%natom3)=Phonon frequencies gamma_in_ph(gams%natom3)=Phonon linewidths. gamma_out_ph(gams%natom3)=Phonon linewidths. lambda_in_ph(gams%natom3)=Phonon linewidths. lambda_out_ph(gams%natom3)=Phonon linewidths.
SOURCE
1123 subroutine phgamma_vv_eval_qibz(gams, cryst, ifc, iq_ibz, spin, phfrq, gamma_in_ph, gamma_out_ph, lambda_in_ph, lambda_out_ph) 1124 1125 !Arguments ------------------------------------ 1126 !scalars 1127 integer,intent(in) :: iq_ibz,spin 1128 type(phgamma_t),intent(inout) :: gams 1129 type(crystal_t),intent(in) :: cryst 1130 type(ifc_type),intent(in) :: ifc 1131 !arrays 1132 real(dp),intent(out) :: phfrq(gams%natom3) 1133 real(dp),intent(out) :: gamma_in_ph(3,3,gams%natom3),lambda_in_ph(3,3,gams%natom3) 1134 real(dp),intent(out) :: gamma_out_ph(3,3,gams%natom3),lambda_out_ph(3,3,gams%natom3) 1135 1136 !Local variables------------------------------- 1137 !scalars 1138 integer,parameter :: qtor0 = 0 1139 integer :: natom3, nu1, idir, jdir, ii 1140 real(dp) :: spinfact 1141 !character(len=500) :: msg 1142 !arrays 1143 real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom), displ_red(2,gams%natom3,gams%natom3) 1144 real(dp) :: work_qnu(gams%natom3), gam_atm(2,gams%natom3,gams%natom3) 1145 1146 ! ************************************************************************* 1147 1148 natom3 = gams%natom3 1149 1150 ! Get phonon frequencies and eigenvectors. 1151 call ifc%fourq(cryst, gams%qibz(:,iq_ibz), phfrq, displ_cart, out_displ_red=displ_red) 1152 1153 do jdir=1,3 1154 do idir=1,3 1155 ii = idir + gams%ndir_transp * (jdir - 1) 1156 1157 ! Scalar product with the displ_red vectors. 1158 gam_atm = reshape(gams%vals_in_qibz(:,ii,:,:,iq_ibz,spin), [2, natom3, natom3]) 1159 call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, work_qnu) 1160 gamma_in_ph(idir,jdir,:) = work_qnu 1161 1162 gam_atm = reshape(gams%vals_out_qibz(:,ii,:,:,iq_ibz,spin), [2, natom3, natom3]) 1163 call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, work_qnu) 1164 gamma_out_ph(idir,jdir,:) = work_qnu 1165 end do ! idir 1166 end do ! jdir 1167 1168 ! TODO : check this - looks like a factor of 2 wrt the inline documentation! 1169 !spinfact should be 1 for a normal non sppol calculation without spinorbit 1170 !for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2 1171 !for sppol 2 it should be 0.5 as we have 2 spin channels to sum 1172 spinfact = two / (gams%nsppol*gams%nspinor) 1173 1174 ! Compute lambda transport. 1175 do jdir =1,3 1176 do idir =1,3 1177 do nu1=1,gams%natom3 1178 gamma_in_ph(idir,jdir,nu1) = gamma_in_ph(idir,jdir,nu1) * pi * spinfact 1179 gamma_out_ph(idir,jdir,nu1) = gamma_out_ph(idir,jdir,nu1) * pi * spinfact 1180 lambda_in_ph(idir,jdir,nu1) = zero 1181 lambda_out_ph(idir,jdir,nu1) = zero 1182 if (abs(phfrq(nu1)) > EPHTK_WTOL) then 1183 lambda_in_ph(idir,jdir,nu1) = gamma_in_ph(idir,jdir,nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2) 1184 lambda_out_ph(idir,jdir,nu1) = gamma_out_ph(idir,jdir,nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2) 1185 end if 1186 end do 1187 end do ! idir 1188 end do ! jdir 1189 1190 end subroutine phgamma_vv_eval_qibz
m_phgamma/phgamma_vv_interp [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
phgamma_vv_interp
FUNCTION
Interpolate the linewidths at a given q-point.
INPUTS
gams<phgamma_t> cryst<crystal_t>=crystalline structure. ifc<ifc_type>=Interatomic force constants. spin=Spin index qpt(3)=q-point in reduced coordinates gamma_ph(3*natom)=Phonon linewidths displ_cart(2,3,cryst%natom,3*cryst%natom)=Phonon displacement in cartesian coordinates.
OUTPUT
gamma_ph(gams%natom3)=Interpolated Phonon linewidths. lamda_ph(3*natom)=Lambda coefficients for the different phonon modes. phfrq(3*natom)=phonon frequencies at current q
SOURCE
1218 subroutine phgamma_vv_interp(gams, cryst, ifc, spin, qpt, phfrq, gamma_in_ph, gamma_out_ph, lambda_in_ph, lambda_out_ph) 1219 1220 !Arguments ------------------------------------ 1221 !scalars 1222 integer,intent(in) :: spin 1223 type(phgamma_t),intent(inout) :: gams 1224 type(crystal_t),intent(in) :: cryst 1225 type(ifc_type),intent(in) :: ifc 1226 !arrays 1227 real(dp),intent(in) :: qpt(3) 1228 real(dp),intent(out) :: phfrq(gams%natom3) 1229 real(dp),intent(out) :: gamma_in_ph(3,3,gams%natom3),lambda_in_ph(3,3,gams%natom3) 1230 real(dp),intent(out) :: gamma_out_ph(3,3,gams%natom3),lambda_out_ph(3,3,gams%natom3) 1231 1232 !Local variables------------------------------- 1233 !scalars 1234 integer,parameter :: qtor0 = 0 1235 integer :: natom3, nu1, idir,jdir,ii 1236 real(dp) :: spinfact 1237 !character(len=500) :: msg 1238 !arrays 1239 real(dp) :: displ_cart(2,3,cryst%natom,3*cryst%natom) 1240 real(dp) :: displ_red(2,gams%natom3,gams%natom3),work_qnu(gams%natom3) 1241 real(dp) :: gam_in_now(2,3,3,gams%natom3**2) 1242 real(dp) :: gam_out_now(2,3,3,gams%natom3**2) 1243 real(dp) :: gam_atm(2,gams%natom3,gams%natom3) 1244 real(dp),allocatable :: coskr(:,:),sinkr(:,:) 1245 1246 ! ************************************************************************* 1247 1248 ! Compute internal tables used for Fourier interpolation. 1249 if (.not.allocated(gams%vals_in_rpt)) call phgamma_vv_interp_setup(gams, cryst) 1250 1251 natom3 = gams%natom3 1252 1253 ! Get phonon frequencies and eigenvectors. 1254 call ifc%fourq(cryst, qpt, phfrq, displ_cart, out_displ_red=displ_red) 1255 1256 ! Taken from mkph_linwid 1257 ! This reduced version of ftgkk supposes the kpoints have been integrated 1258 ! in integrate_gamma. Do FT from real-space gamma grid to 1 qpt. 1259 ABI_MALLOC(coskr, (1,gams%nrpt)) 1260 ABI_MALLOC(sinkr, (1,gams%nrpt)) 1261 ! TODO: This is not optimal 1262 call ftgam_init(gams%gprim, 1, gams%nrpt, qpt, gams%rpt, coskr, sinkr) 1263 1264 do idir=1,3 1265 do jdir=1,3 1266 ii = idir+gams%ndir_transp*(jdir-1) 1267 1268 call ftgam(gams%wghatm, gam_in_now, gams%vals_in_rpt(:,ii,:,:,spin), gams%natom, 1, gams%nrpt, qtor0, coskr, sinkr) 1269 call ftgam(gams%wghatm, gam_out_now, gams%vals_out_rpt(:,ii,:,:,spin), gams%natom, 1, gams%nrpt, qtor0, coskr, sinkr) 1270 1271 ! This call is not executed in elphon! 1272 ! TODO: needs to take into account the matrix nature of _in_ and _out_ gammas 1273 if (gams%symgamma == 1) then 1274 call tgamma_symm(cryst, qpt, gam_in_now) 1275 call tgamma_symm(cryst, qpt, gam_out_now) 1276 end if 1277 1278 ! Scalar product with the displ_red 1279 gam_atm = reshape(gam_in_now, [2, natom3, natom3]) 1280 call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, work_qnu) 1281 gamma_in_ph(idir,jdir,:) = work_qnu 1282 1283 gam_atm = reshape(gam_out_now, [2, natom3, natom3]) 1284 call ephtk_gam_atm2qnu(natom3, displ_red, gam_atm, work_qnu) 1285 gamma_out_ph(idir,jdir,:) = work_qnu 1286 end do 1287 end do 1288 1289 ABI_FREE(coskr) 1290 ABI_FREE(sinkr) 1291 1292 ! Compute lambda 1293 ! spinfact should be 1 for a normal non sppol calculation without spinorbit 1294 ! for spinors it should also be 1 as bands are twice as numerous but n0 has been divided by 2 1295 ! for sppol 2 it should be 0.5 as we have 2 spin channels to sum 1296 spinfact = two / (gams%nsppol * gams%nspinor) 1297 1298 ! Compute lambda 1299 do nu1=1,gams%natom3 1300 do idir=1,3 1301 do jdir=1,3 1302 gamma_in_ph(idir,jdir,nu1) = gamma_in_ph(idir,jdir,nu1) * pi * spinfact 1303 gamma_out_ph(idir,jdir,nu1) = gamma_out_ph(idir,jdir,nu1) * pi * spinfact 1304 lambda_in_ph(idir,jdir,nu1) = zero 1305 lambda_out_ph(idir,jdir,nu1) = zero 1306 if (abs(phfrq(nu1)) > EPHTK_WTOL) then 1307 lambda_in_ph(idir,jdir,nu1) = gamma_in_ph(idir,jdir,nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2) 1308 lambda_out_ph(idir,jdir,nu1) = gamma_out_ph(idir,jdir,nu1) / (two * pi * gams%n0(spin) * phfrq(nu1)**2) 1309 end if 1310 end do 1311 end do 1312 end do 1313 1314 end subroutine phgamma_vv_interp
m_phgamma/phgamma_vv_interp_setup [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
phgamma_vv_interp_setup
FUNCTION
This routines prepare the internal tables used to interpolate the vv_linewidths in q-space
INPUTS
action =
"INIT" to allocate and compute the internal tables (default)
"FREE" to deallocate the internal tables.
SIDE EFFECTS
gams<phgamma_t>= gams%vals_in_rpt, etc... depending on action.
SOURCE
1336 subroutine phgamma_vv_interp_setup(gams, cryst) 1337 1338 !Arguments ------------------------------------ 1339 !scalars 1340 type(phgamma_t),intent(inout) :: gams 1341 type(crystal_t),intent(in) :: cryst 1342 1343 !Local variables------------------------------- 1344 !scalars 1345 integer,parameter :: qtor1 = 1 1346 integer :: iq_bz,iq_ibz,spin,ierr, ii, idir, jdir 1347 !character(len=500) :: msg 1348 !arrays 1349 integer,allocatable :: qirredtofull(:),qpttoqpt(:,:,:) 1350 real(dp),allocatable :: coskr(:,:),sinkr(:,:) 1351 real(dp),allocatable :: vals_in_bz(:,:,:,:,:), vals_out_bz(:,:,:,:,:) 1352 1353 ! ************************************************************************* 1354 1355 ABI_MALLOC_OR_DIE(vals_in_bz,(2, 9, gams%natom3**2, gams%nqbz, gams%nsppol), ierr) 1356 ABI_MALLOC_OR_DIE(vals_out_bz,(2, 9, gams%natom3**2, gams%nqbz, gams%nsppol), ierr) 1357 vals_in_bz = zero 1358 vals_out_bz = zero 1359 1360 ! Build tables needed by complete_gamma. 1361 call ephtk_mkqtabs(cryst, gams%nqibz, gams%qibz, gams%nqbz, gams%qbz, qirredtofull, qpttoqpt) 1362 1363 ! Fill BZ array with IBZ data. 1364 do spin=1,gams%nsppol 1365 do iq_ibz=1,gams%nqibz 1366 iq_bz = qirredtofull(iq_ibz) 1367 vals_in_bz(:,:,:,iq_bz,spin) = reshape(gams%vals_in_qibz(:,:,:,:,iq_ibz,spin), [2, 9, gams%natom3**2]) 1368 vals_out_bz(:,:,:,iq_bz,spin) = reshape(gams%vals_out_qibz(:,:,:,:,iq_ibz,spin), [2, 9, gams%natom3**2]) 1369 end do 1370 end do 1371 1372 ! Complete vals_bz in the full BZ. 1373 !TODO!!! rotate the vv in and out matrices, according to the symmetry operation, instead of just copying them 1374 !call complete_gamma_vv(cryst, gams%natom3, gams%nsppol, gams%nqibz, gams%nqbz,& 1375 ! gams%eph_scalprod, qirredtofull, qpttoqpt, vals_in_bz(:,:,:,:,spin)) 1376 !call complete_gamma_vv(cryst, gams%natom3, gams%nsppol, gams%nqibz, gams%nqbz,& 1377 ! gams%eph_scalprod, qirredtofull, qpttoqpt, vals_out_bz(:,:,:,:,spin)) 1378 1379 ! TODO: replace the above with these calls from anaddb 1380 !call complete_gamma_tr(cryst,elph_ds%ep_scalprod,elph_ds%nbranch,elph_ds%nqptirred,& 1381 !& elph_ds%nqpt_full,elph_ds%nsppol,elph_tr_ds%gamma_qpt_trout,elph_ds%qirredtofull,qpttoqpt) 1382 1383 ! TODO: idem for vv_vals 3x3 matrices 1384 1385 ABI_FREE(qirredtofull) 1386 ABI_FREE(qpttoqpt) 1387 1388 !! This call is not executed in elphon! 1389 !if (gams%symgamma == 1) then 1390 ! do spin=1,gams%nsppol 1391 ! do iq_bz=1,gams%nqbz 1392 ! call tgamma_symm_vv(cryst, gams%qbz(:,iq_bz), vals_in_bz(:,:,:,iq_bz,spin)) 1393 ! call tgamma_symm_vv(cryst, gams%qbz(:,iq_bz), vals_out_bz(:,:,:,iq_bz,spin)) 1394 ! end do 1395 ! end do 1396 !end if 1397 1398 ! Now FT to real space 1399 ! NOTE: gprim (not gprimd) is used for all FT interpolations, 1400 ! to be consistent with the dimensions of the rpt, which come from anaddb. 1401 ! TODO: this is needed only if FT is used, not when the linear interpolation is employed. 1402 !if (.not. allocated(gams%vals_in_rpt)) then 1403 ABI_MALLOC_OR_DIE(gams%vals_in_rpt, (2, 9, gams%natom3**2, gams%nrpt, gams%nsppol), ierr) 1404 ABI_MALLOC_OR_DIE(gams%vals_out_rpt, (2, 9, gams%natom3**2, gams%nrpt, gams%nsppol), ierr) 1405 gams%vals_in_rpt = zero 1406 gams%vals_out_rpt = zero 1407 1408 ! q --> r 1409 ABI_MALLOC(coskr, (gams%nqbz,gams%nrpt)) 1410 ABI_MALLOC(sinkr, (gams%nqbz,gams%nrpt)) 1411 call ftgam_init(gams%gprim, gams%nqbz, gams%nrpt, gams%qbz, gams%rpt, coskr, sinkr) 1412 1413 do spin=1,gams%nsppol 1414 do idir=1,3 1415 do jdir=1,3 1416 ii = idir+gams%ndir_transp*(jdir-1) 1417 ! TODO: this is no contiguous in memory and will be slow. Make adapted ftgam? 1418 call ftgam(gams%wghatm, vals_in_bz(:,ii,:,:,spin), gams%vals_in_rpt(:,ii,:,:,spin), gams%natom, gams%nqbz,& 1419 gams%nrpt, qtor1, coskr, sinkr) 1420 call ftgam(gams%wghatm, vals_out_bz(:,ii,:,:,spin), gams%vals_out_rpt(:,ii,:,:,spin), gams%natom, gams%nqbz,& 1421 gams%nrpt, qtor1, coskr, sinkr) 1422 end do 1423 end do 1424 end do 1425 1426 ABI_FREE(vals_in_bz) 1427 ABI_FREE(vals_out_bz) 1428 ABI_FREE(coskr) 1429 ABI_FREE(sinkr) 1430 1431 end subroutine phgamma_vv_interp_setup
m_phgamma/tgamma_symm [ Functions ]
[ Top ] [ m_phgamma ] [ Functions ]
NAME
tgamma_symm
FUNCTION
Symmetrize the tgamma matrix
INPUTS
qpt(3)=phonon wavevector in reduced coordinates. cryst<crystal_t>=Crystalline structure.
SOURCE
721 subroutine tgamma_symm(cryst, qpt, tgamma) 722 723 !Arguments ------------------------------------ 724 !scalars 725 type(crystal_t),intent(in) :: cryst 726 !arrays 727 real(dp),intent(in) :: qpt(3) 728 real(dp),intent(inout) :: tgamma(2,3*cryst%natom,3*cryst%natom) 729 730 !Local variables------------------------------- 731 !scalars 732 integer :: ii,natom3,k 733 !arrays 734 real(dp) :: tgcart(2,3*cryst%natom,3*cryst%natom) 735 real(dp) :: umat(2,3*cryst%natom,3*cryst%natom),tmp_mat(2,3*cryst%natom,3*cryst%natom) 736 737 ! ********************************************************************* 738 739 ! Build U matrix. 740 umat = zero; k = 1 741 do ii=1,cryst%natom 742 umat(1,k:k+2, k:k+2) = cryst%gprimd 743 k = k + 3 744 end do 745 746 natom3 = 3 * cryst%natom 747 748 ! Reduced --> Cartesian 749 call zgemm('N', 'N', natom3, natom3, natom3, cone, tgamma, natom3, umat, natom3, czero, tmp_mat, natom3) 750 call zgemm('T', 'N', natom3, natom3, natom3, cone, umat, natom3, tmp_mat, natom3, czero, tgcart, natom3) 751 752 ! Make the matrix hermitian 753 call mkherm(tgcart, 3*cryst%natom) 754 755 ! Symmetrize tgamma matrix. 756 call symdyma(tgcart, cryst%indsym, cryst%natom, cryst%nsym, qpt, cryst%rprimd, cryst%symrel, cryst%symafm) 757 758 umat = zero; k = 1 759 do ii=0,cryst%natom-1 760 umat(1,k:k+2, k:k+2) = cryst%rprimd 761 k = k + 3 762 end do 763 764 ! Cartesian --> Reduced 765 call zgemm('N', 'N', natom3, natom3, natom3, cone, tgcart, natom3, umat, natom3, czero, tmp_mat, natom3) 766 call zgemm('T', 'N', natom3, natom3, natom3, cone, umat, natom3, tmp_mat, natom3, czero, tgamma, natom3) 767 768 end subroutine tgamma_symm