TABLE OF CONTENTS
- ABINIT/m_iogkk
- ABINIT/nmsq_gam
- ABINIT/nmsq_gam_sumfs
- ABINIT/nmsq_pure_gkk
- ABINIT/nmsq_pure_gkk_sumfs
- ABINIT/normsq_gkq
- m_iogkk/completeperts
- m_iogkk/inpgkk
- m_iogkk/outgkk
- m_iogkk/prt_gkk_yambo
- m_iogkk/read_el_veloc
- m_iogkk/read_gkk
ABINIT/m_iogkk [ Modules ]
NAME
m_iogkk
FUNCTION
IO routines for GKK files
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MVer) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_iogkk 23 24 use defs_basis 25 use defs_elphon 26 use m_errors 27 use m_abicore 28 use m_xmpi 29 use m_krank 30 use m_hdr 31 32 use defs_datatypes, only : ebands_t 33 use defs_abitypes, only : MPI_type 34 use m_numeric_tools, only : wrap2_pmhalf 35 use m_io_tools, only : open_file, get_unit 36 use m_symtk, only : mati3inv, littlegroup_q 37 use m_geometry, only : phdispl_cart2red, littlegroup_pert 38 use m_crystal, only : crystal_t 39 use m_ifc, only : ifc_type 40 use m_dynmat, only : d2sym3 41 42 implicit none 43 44 private
ABINIT/nmsq_gam [ Functions ]
NAME
nmsq_gam
FUNCTION
Calculate gamma matrices keeping full dependence on bands from original h1_mat_el_sq matrix elements (no averaging over bands near the Fermi surface)
INPUTS
displ_red = phonon mode displacement vectors, post-multiplied by gprim matrix
(ie. turned to reduced coordinates)
eigvec = phonon eigenvectors
elph_ds = datastructure with gkk matrix elements
FSfullpqtofull = mapping of k+q to k
kpt_phon = coordinates of kpoints near to FS
h1_mat_el_sq = matrix elements $<psi_{k+q,m} | H^{1} | psi_{k,n}>$ squared
iqptirred = index of present qpoint
OUTPUT
accum_mat = matrix for accumulating FS average of gkk (gamma matrix -> linewidths) accum_mat2 = matrix for accumulating FS average of gamma matrix with good prefactors
SOURCE
1596 subroutine nmsq_gam (accum_mat,accum_mat2,displ_red,eigvec,elph_ds,FSfullpqtofull,& 1597 & h1_mat_el_sq,iqptirred) 1598 1599 !Arguments ------------------------------------ 1600 !scalars 1601 integer,intent(in) :: iqptirred 1602 type(elph_type),intent(inout) :: elph_ds 1603 !arrays 1604 integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full) 1605 real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch) 1606 real(dp),intent(in) :: eigvec(2,elph_ds%nbranch,elph_ds%nbranch) 1607 real(dp),intent(inout) :: & 1608 & h1_mat_el_sq(2,elph_ds%nFSband*elph_ds%nFSband,elph_ds%nbranch*elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol) 1609 real(dp),intent(inout) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol) 1610 real(dp),intent(inout) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol) 1611 1612 !Local variables------------------------------- 1613 ! tmp variables for diagonalization 1614 !scalars 1615 integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ibranch,isppol,ipert1 1616 integer :: jbranch 1617 integer :: iqpt_fullbz 1618 integer :: ik_this_proc 1619 real(dp) :: sd1,sd2 1620 character(len=500) :: message 1621 !arrays 1622 real(dp) :: gkq_1band(2,elph_ds%nbranch,elph_ds%nbranch) 1623 real(dp) :: tmp_mat2(2,elph_ds%nbranch,elph_ds%nbranch) 1624 real(dp) :: zgemm_tmp_mat(2,elph_ds%nbranch,elph_ds%nbranch) 1625 1626 ! ************************************************************************* 1627 1628 if (elph_ds%ep_keepbands == 0) then 1629 write (message,'(a,i0)')' elph_ds%ep_keepbands should be 1 while is ',elph_ds%ep_keepbands 1630 ABI_ERROR(message) 1631 end if 1632 1633 !MG20060603 NOTE: 1634 !accum_mat and accum_mat2 are real, the imaginary part is used for debugging purpose 1635 !accum_mat2 is used to store the phonon-linewidhts before interpolation 1636 1637 iqpt_fullbz = elph_ds%qirredtofull(iqptirred) 1638 write(std_out,*) 'nmsq_gam : iqptirred = ', iqptirred 1639 1640 do isppol=1,elph_ds%nsppol 1641 do ik_this_proc =1, elph_ds%k_phon%my_nkpt 1642 ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc) 1643 1644 ikpt_phonq = FSfullpqtofull(ikpt_phon,iqpt_fullbz) 1645 1646 do ib1=1,elph_ds%nFSband 1647 sd1 = elph_ds%k_phon%wtk(ib1,ikpt_phon,isppol) !weights for distance from the fermi surface 1648 1649 do ib2=1,elph_ds%nFSband 1650 sd2 = elph_ds%k_phon%wtk(ib2,ikpt_phonq,isppol) !weights for distance from the fermi surface 1651 ibeff = ib2+elph_ds%nFSband*(ib1-1) 1652 1653 gkq_1band(:,:,:) = zero 1654 1655 zgemm_tmp_mat= reshape (h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol),(/2,elph_ds%nbranch,elph_ds%nbranch/)) 1656 1657 call gam_mult_displ(elph_ds%nbranch, displ_red, zgemm_tmp_mat, tmp_mat2) 1658 1659 ! sum over bands 1660 do ipert1=1,elph_ds%nbranch 1661 gkq_1band(1,ipert1,ipert1) = gkq_1band(1,ipert1,ipert1) + tmp_mat2(1,ipert1,ipert1) 1662 end do 1663 1664 ! summing over k points and bands, still diagonal in jbranch 1665 accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_1band(:,:,:)*sd1*sd2 1666 1667 ! MG20060603 : summing over bands and kpoints with weights to calculate the phonon linewidth 1668 do jbranch=1,elph_ds%nbranch 1669 accum_mat2(:,jbranch,jbranch,isppol) = accum_mat2(:,jbranch,jbranch,isppol) + gkq_1band(:,jbranch,jbranch)*sd1*sd2 1670 end do 1671 ! END MG 1672 1673 1674 ! now turn to cartesian coordinates 1675 1676 ! Final Gamma matrix (hermitian) = E * D_g * E^{+} 1677 ! Where E^{+} is the hermitian conjugate of the eigenvector matrix E 1678 ! And D_g is the diagonal matrix of values of gamma for this qpoint 1679 1680 ! Here gkq_1band is indexed with real phonon modes (not atom+idir) 1681 ! turn gkq_1band to atom+cartesian coordinates (instead of normal coordinates for qpoint) 1682 tmp_mat2(:,:,:) = zero 1683 do ibranch =1,elph_ds%nbranch 1684 do jbranch =1,elph_ds%nbranch 1685 tmp_mat2(1,ibranch,jbranch) = tmp_mat2(1,ibranch,jbranch) + & 1686 & eigvec(1,ibranch,jbranch) * gkq_1band(1,jbranch,jbranch) 1687 tmp_mat2(2,ibranch,jbranch) = tmp_mat2(2,ibranch,jbranch) + & 1688 & eigvec(2,ibranch,jbranch) * gkq_1band(1,jbranch,jbranch) 1689 end do 1690 end do 1691 gkq_1band(:,:,:) = zero 1692 1693 ! here eigvec is transposed and complexconjugated. 1694 zgemm_tmp_mat=zero 1695 call zgemm('n','c',elph_ds%nbranch,elph_ds%nbranch,elph_ds%nbranch,cone,& 1696 & tmp_mat2,elph_ds%nbranch,eigvec,elph_ds%nbranch,czero,zgemm_tmp_mat,elph_ds%nbranch) 1697 1698 gkq_1band = zgemm_tmp_mat 1699 1700 ! gamma matrix contribution in cartesian coordinates (ie interpolatable form) 1701 h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol) = reshape(gkq_1band,(/2,elph_ds%nbranch*elph_ds%nbranch/)) 1702 1703 end do 1704 end do 1705 ! END loop over bands ib1 ib2 1706 1707 end do 1708 ! END loop over kpt_phon 1709 end do 1710 !END loop over nsppol 1711 1712 1713 end subroutine nmsq_gam
ABINIT/nmsq_gam_sumfs [ Functions ]
NAME
nmsq_gam_sumfs
FUNCTION
Calculate gamma matrices from original h1_mat_el_sq matrix elements averaging over bands near the Fermi surface
INPUTS
displ_red = phonon mode displacement vectors, post-multiplied by gprim matrix
(ie. turned to reduced coordinates)
eigvec = eigenvectors of phonons (to turn to cartesian coord frame)
elph_ds = datastructure with gkk matrix elements
FSfullpqtofull = mapping of k+q to k
kpt_phon = coordinates of kpoints near to FS
h1_mat_el_sq = matrix elements $<psi_{k+q,m} | H^{1} | psi_{k,n}>$ squared
iqptirred = index of present qpoint
OUTPUT
accum_mat = matrix for accumulating FS average of gkk (gamma matrix -> linewidths) accum_mat2 = matrix for accumulating FS average of gamma matrix with good prefactors
SOURCE
1741 subroutine nmsq_gam_sumFS(accum_mat,accum_mat2,displ_red,eigvec,elph_ds,FSfullpqtofull,& 1742 & h1_mat_el_sq,iqptirred) 1743 1744 !Arguments ------------------------------------ 1745 !scalars 1746 integer,intent(in) :: iqptirred 1747 type(elph_type),intent(inout) :: elph_ds 1748 !arrays 1749 integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full) 1750 real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch) 1751 real(dp),intent(in) :: eigvec(2,elph_ds%nbranch,elph_ds%nbranch) 1752 real(dp),intent(inout) :: & 1753 & h1_mat_el_sq(2,elph_ds%nFSband*elph_ds%nFSband,elph_ds%nbranch*elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol) 1754 real(dp),intent(inout) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol) 1755 real(dp),intent(inout) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol) 1756 1757 !Local variables------------------------------- 1758 !scalars 1759 integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ibranch,ipert1,isppol,jbranch,iqpt_fullbz 1760 integer :: ik_this_proc 1761 real(dp) :: sd1,sd2 1762 character(len=500) :: message 1763 !arrays 1764 real(dp) :: gkq_sum_bands(2,elph_ds%nbranch,elph_ds%nbranch) 1765 real(dp) :: tmp_gkq_sum_bands(2,elph_ds%nbranch,elph_ds%nbranch) 1766 real(dp) :: tmp_mat2(2,elph_ds%nbranch,elph_ds%nbranch) 1767 real(dp),allocatable :: zgemm_tmp_mat(:,:,:) 1768 1769 ! ************************************************************************* 1770 1771 if (elph_ds%ep_keepbands /= 0) then 1772 write (message,'(a,i0)')' elph_ds%ep_keepbands should be 0 in order to average over bands!',elph_ds%ep_keepbands 1773 ABI_ERROR(message) 1774 end if 1775 1776 iqpt_fullbz = elph_ds%qirredtofull(iqptirred) 1777 1778 1779 !MG20060603 NOTE: 1780 !accum_mat and accum_mat2 are real, the imaginary part is used for debugging purpose 1781 !accum_mat2 is used to store the phonon-linewidhts before interpolation 1782 1783 ABI_MALLOC(zgemm_tmp_mat ,(2,elph_ds%nbranch,elph_ds%nbranch)) 1784 1785 do isppol=1,elph_ds%nsppol 1786 do ik_this_proc =1, elph_ds%k_phon%my_nkpt 1787 ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc) 1788 1789 ikpt_phonq = FSfullpqtofull(ikpt_phon,iqpt_fullbz) 1790 1791 gkq_sum_bands = zero 1792 tmp_gkq_sum_bands = zero 1793 1794 1795 do ib1=1,elph_ds%nFSband 1796 ! weights for distance from the fermi surface 1797 sd1 = elph_ds%k_phon%wtk(ib1,ikpt_phon,isppol) 1798 1799 do ib2=1,elph_ds%nFSband 1800 ! weights for distance from the fermi surface 1801 sd2 = elph_ds%k_phon%wtk(ib2,ikpt_phonq,isppol) 1802 ibeff=ib2+(ib1-1)*elph_ds%nFSband 1803 1804 zgemm_tmp_mat = reshape(h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol),(/2,elph_ds%nbranch,elph_ds%nbranch/)) 1805 1806 call gam_mult_displ(elph_ds%nbranch, displ_red, zgemm_tmp_mat, tmp_mat2) 1807 1808 ! sum over bands in gkq_sum_bands 1809 do ipert1=1,elph_ds%nbranch 1810 gkq_sum_bands(1,ipert1,ipert1) = gkq_sum_bands(1,ipert1,ipert1) + sd1*sd2*tmp_mat2(1,ipert1,ipert1) 1811 end do 1812 1813 1814 1815 end do 1816 end do 1817 ! END loop over bands 1818 1819 ! summing over k points, still diagonal in jbranch 1820 accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_sum_bands(:,:,:) 1821 accum_mat2(:,:,:,isppol) = accum_mat2(:,:,:,isppol) + gkq_sum_bands(:,:,:) 1822 1823 ! summed over bands, now turn to cartesian coordinates 1824 1825 ! Final Gamma matrix (hermitian) = E * D_g * E^{+} 1826 ! Where E^{+} is the hermitian conjugate of the eigenvector matrix E 1827 ! And D_g is the diagonal matrix of values of gamma for this qpoint 1828 1829 ! Here gkq_sum_bands is indexed with real phonon modes (not atom+idir) 1830 ! turn gkq_sum_bands to atom+cartesian coordinates (instead of normal coordinates for qpoint) 1831 ! This is not a full matrix multiplication, just vector one, by 1832 ! gkq_sum_bands(1,jbranch,jbranch) 1833 tmp_mat2(:,:,:) = zero 1834 do ibranch =1,elph_ds%nbranch 1835 do jbranch =1,elph_ds%nbranch 1836 tmp_mat2(1,ibranch,jbranch) = tmp_mat2(1,ibranch,jbranch) + & 1837 & eigvec(1,ibranch,jbranch) * & 1838 & gkq_sum_bands(1,jbranch,jbranch) 1839 tmp_mat2(2,ibranch,jbranch) = tmp_mat2(2,ibranch,jbranch) + & 1840 & eigvec(2,ibranch,jbranch) * & 1841 & gkq_sum_bands(1,jbranch,jbranch) 1842 end do 1843 end do 1844 1845 ! here eigvec is transposed and complex conjugated. 1846 zgemm_tmp_mat=zero 1847 call zgemm('n','c',elph_ds%nbranch,elph_ds%nbranch,elph_ds%nbranch,cone,& 1848 & tmp_mat2,elph_ds%nbranch,eigvec,elph_ds%nbranch,czero,zgemm_tmp_mat,elph_ds%nbranch) 1849 1850 gkq_sum_bands = zgemm_tmp_mat 1851 1852 ! ! gamma matrix contribution in cartesian coordinates (ie interpolatable form) 1853 ! gamma matrix contribution in reduced coordinates (ie interpolatable form) 1854 h1_mat_el_sq(:,1,:,ik_this_proc,isppol) = reshape(gkq_sum_bands(:,:,:),(/2,elph_ds%nbranch*elph_ds%nbranch/)) 1855 1856 ! accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_sum_bands(:,:,:) 1857 end do 1858 ! END loop over kpt_phon 1859 end do 1860 !END loop over sppol 1861 1862 ABI_FREE(zgemm_tmp_mat) 1863 1864 end subroutine nmsq_gam_sumFS
ABINIT/nmsq_pure_gkk [ Functions ]
NAME
nmsq_pure_gkk
FUNCTION
Calculate gamma matrices for pure gkk case, ie when the scalar product with the displacement vector is done later Sum over bands is carried out later.
INPUTS
displ_red = phonon displacement in reduced coordinates (used to calculate the ph linewidth)
elph_ds = datastructure with gkk matrix elements
FSfullpqtofull = mapping of k+q to k
kpt_phon = coordinates of kpoints near to FS
h1_mat_el_sq = matrix elements $<psi_{k+q,m} | H^{1} | psi_{k,n}>$ squared
iqptirred = index of present qpoint
OUTPUT
elph_ds%gkq filled accum_mat = matrix for accumulating FS average of gkk (gamma matrix -> linewidths) accum_mat2 = complex array whose real part contains the phonon linewidth
SOURCE
1893 subroutine nmsq_pure_gkk(accum_mat,accum_mat2,displ_red,elph_ds,FSfullpqtofull,& 1894 & h1_mat_el_sq,iqptirred) 1895 1896 !Arguments ------------------------------------ 1897 !scalars 1898 integer,intent(in) :: iqptirred 1899 type(elph_type),intent(inout) :: elph_ds 1900 !arrays 1901 integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full) 1902 real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch) 1903 real(dp),intent(inout) :: & 1904 & h1_mat_el_sq(2,elph_ds%nFSband*elph_ds%nFSband,elph_ds%nbranch*elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol) 1905 real(dp),intent(inout) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol) 1906 real(dp),intent(inout) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol) 1907 1908 !Local variables------------------------------- 1909 !scalars 1910 integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ipert1,isppol 1911 integer :: iqpt_fullbz 1912 integer :: ik_this_proc 1913 real(dp) :: sd1,sd2 1914 character(len=500) :: message 1915 !arrays 1916 real(dp) :: gkq_sum_bands(2,elph_ds%nbranch,elph_ds%nbranch) 1917 real(dp) :: tmp_mat2(2,elph_ds%nbranch,elph_ds%nbranch) 1918 real(dp) :: zgemm_tmp_mat(2,elph_ds%nbranch,elph_ds%nbranch) 1919 1920 ! ************************************************************************* 1921 1922 if (elph_ds%ep_keepbands /= 1) then 1923 message = ' elph_ds%ep_keepbands should be 1 to keep bands!' 1924 ABI_ERROR(message) 1925 end if 1926 1927 iqpt_fullbz = elph_ds%qirredtofull(iqptirred) 1928 1929 !h1_mat_el_sq is already fine here - nothing to do 1930 1931 1932 !MG20060603 NOTE: 1933 !accum_mat and accum_mat2 are real, the imaginary part is used for debugging purpose 1934 !accum_mat2 is used to store the phonon-linewidhts before interpolation 1935 1936 !MJV 20070525 NOTE: 1937 !in some of the nmsq routines, in particular this one, the work done to 1938 !calculate accum_mat,accum_mat2 is completely superfluous and will be re-done 1939 !on the interpolated values. 1940 !MG uses them for the QPT output, however, so keep it for consistency for the 1941 !moment. 1942 1943 do isppol=1,elph_ds%nsppol 1944 do ik_this_proc =1, elph_ds%k_phon%my_nkpt 1945 ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc) 1946 1947 ikpt_phonq = FSfullpqtofull(ikpt_phon,iqpt_fullbz) 1948 1949 gkq_sum_bands(:,:,:) = zero 1950 1951 ! gkq_sum_bands = \sum_{ib1,ib2} \langle k+q \mid H^{(1)}_{q,\tau_i,\alpha_i} \mid k \rangle 1952 ! \cdot \langle k \mid H^{(1)}_{q,\tau_j,\alpha_j} \mid k+q \rangle 1953 ! where ibranch -> \tau_i,\alpha_i and jbranch -> \tau_j,\alpha_j 1954 1955 do ib1=1,elph_ds%nFSband 1956 1957 sd1 = elph_ds%k_phon%wtk(ib1,ikpt_phon,isppol) ! weights for distance from the fermi surface 1958 1959 do ib2=1,elph_ds%nFSband 1960 1961 sd2 = elph_ds%k_phon%wtk(ib2,ikpt_phonq,isppol) ! weights for distance from the fermi surface 1962 ibeff = ib2+(ib1-1)*elph_ds%nFSband 1963 1964 gkq_sum_bands = gkq_sum_bands + & 1965 & sd1*sd2*reshape(h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol),(/2,elph_ds%nbranch,elph_ds%nbranch/)) 1966 1967 end do !ib2 1968 end do !ib1 1969 ! END loops over bands 1970 1971 1972 accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_sum_bands(:,:,:) 1973 end do 1974 ! END loop over kpt_phon 1975 1976 ! MG20060603 1977 ! do scalar product with the displ_red to calculate the ph lwdth before interpolation (stored in accum_mat2) 1978 1979 zgemm_tmp_mat = accum_mat(:,:,:,isppol) 1980 1981 call gam_mult_displ(elph_ds%nbranch, displ_red, zgemm_tmp_mat, tmp_mat2) 1982 1983 do ipert1=1,elph_ds%nbranch 1984 accum_mat2(1,ipert1,ipert1,isppol) = accum_mat2(1,ipert1,ipert1,isppol) + tmp_mat2(1,ipert1,ipert1) 1985 end do 1986 1987 ! ENDMG 1988 1989 end do ! isppol 1990 1991 end subroutine nmsq_pure_gkk
ABINIT/nmsq_pure_gkk_sumfs [ Functions ]
NAME
nmsq_pure_gkk_sumfs
FUNCTION
Calculate gamma matrices for pure gkk case, i.e, when the scalar product with the displacement vector is done later Sum over bands is carried out now.
INPUTS
displ_red = phonon displacement in reduced coordinates (used to calculate the ph linewidth)
elph_ds = datastructure with gkk matrix elements
FSfullpqtofull = mapping of k+q to k
kpt_phon = coordinates of kpoints near to FS
h1_mat_el_sq = matrix elements $<psi_{k+q,m} | H^{1} | psi_{k,n}>$ squared
iqptirred = index of present qpoint
OUTPUT
accum_mat = matrix for accumulating FS average of gkk (gamma matrix -> linewidths) accum_mat2 = complex array whose real part contains the phonon linewidth
SOURCE
2018 subroutine nmsq_pure_gkk_sumfs(accum_mat,accum_mat2,displ_red,elph_ds,FSfullpqtofull,h1_mat_el_sq,iqptirred) 2019 2020 !Arguments ------------------------------------ 2021 !scalars 2022 integer,intent(in) :: iqptirred 2023 type(elph_type),intent(in) :: elph_ds 2024 !arrays 2025 integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full) 2026 real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch) 2027 real(dp),intent(inout) :: & 2028 & h1_mat_el_sq(2,elph_ds%nFSband*elph_ds%nFSband,elph_ds%nbranch*elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol) 2029 real(dp),intent(inout) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol) 2030 real(dp),intent(inout) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol) 2031 2032 !Local variables------------------------------- 2033 !scalars 2034 integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ipert1,isppol,iqpt_fullbz 2035 integer :: nbranch,nsppol,nFSband,nkpt_phon 2036 integer :: ik_this_proc 2037 real(dp) :: sd1,sd2 2038 !character(len=500) :: message 2039 !arrays 2040 real(dp) :: gkq_sum_bands(2,elph_ds%nbranch,elph_ds%nbranch) 2041 real(dp) :: tmp_mat2(2,elph_ds%nbranch,elph_ds%nbranch) 2042 real(dp) :: zgemm_tmp_mat(2,elph_ds%nbranch,elph_ds%nbranch) 2043 2044 ! ************************************************************************* 2045 2046 if (elph_ds%ep_keepbands /= 0) then 2047 ABI_BUG('ep_keepbands should be 0 to average over bands!') 2048 end if 2049 2050 nbranch = elph_ds%nbranch 2051 nsppol = elph_ds%nsppol 2052 nFSband = elph_ds%nFSband 2053 nkpt_phon = elph_ds%k_phon%nkpt 2054 2055 iqpt_fullbz = elph_ds%qirredtofull(iqptirred) 2056 2057 !MG20060603 NOTE: 2058 !accum_mat and accum_mat2 are real, the imaginary part is used for debugging purpose 2059 !accum_mat2 is used to store the phonon-linewidhts before interpolation 2060 2061 do isppol=1,nsppol 2062 do ik_this_proc =1, elph_ds%k_phon%my_nkpt 2063 ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc) 2064 2065 ! 2066 ! The index of k+q in the BZ. 2067 ikpt_phonq = FSfullpqtofull(ikpt_phon,iqpt_fullbz) 2068 ! 2069 ! gkq_sum_bands = 2070 ! \sum_{ib1,ib2} <k+q| H^{(1)}_{q,\tau_i,\alpha_i} |k> \cdot <k| H^{(1)}_{q,\tau_j,\alpha_j}|k+q> 2071 ! 2072 ! where ibranch = (\tau_i,\alpha_i) and jbranch = (\tau_j,\alpha_j). 2073 gkq_sum_bands(:,:,:) = zero 2074 2075 do ib1=1,nFSband 2076 sd1 = elph_ds%k_phon%wtk(ib1,ikpt_phon,isppol) ! weights for distance from the fermi surface 2077 2078 do ib2=1,nFSband 2079 sd2 = elph_ds%k_phon%wtk(ib2,ikpt_phonq,isppol) ! weights for distance from the fermi surface 2080 ibeff=ib2+(ib1-1)*nFSband 2081 2082 gkq_sum_bands = gkq_sum_bands + & 2083 & sd1*sd2* reshape(h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol),(/2,nbranch,nbranch/)) 2084 end do !ib2 2085 end do !ib1 2086 ! 2087 ! gamma matrix contribution in reduced coordinates (ie interpolatable form) 2088 ! The sum over Fermi surface bands is done here, and fed into (ib1,ib2)=(1,1) 2089 h1_mat_el_sq(:,1,:,ik_this_proc,isppol) = reshape(gkq_sum_bands,(/2,nbranch**2/)) 2090 2091 accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_sum_bands(:,:,:) 2092 end do ! kpt_phon 2093 end do ! isppol 2094 ! 2095 !MG20060603 2096 !do scalar product wit displ_red to calculate the ph lwdth before interpolation (stored in accum_mat2) 2097 do isppol=1,nsppol 2098 zgemm_tmp_mat = accum_mat(:,:,:,isppol) 2099 ! 2100 call gam_mult_displ(nbranch, displ_red, zgemm_tmp_mat, tmp_mat2) 2101 2102 do ipert1=1,nbranch 2103 accum_mat2(1,ipert1,ipert1,isppol) = accum_mat2(1,ipert1,ipert1,isppol) + tmp_mat2(1,ipert1,ipert1) 2104 end do 2105 ! 2106 end do 2107 2108 end subroutine nmsq_pure_gkk_sumfs
ABINIT/normsq_gkq [ Functions ]
NAME
normsq_gkq
FUNCTION
This routine takes the gkq matrix elements for a given qpoint, does the scalar product with the phonon displacement vector, squares the gkq matrix elements multiplies by the appropriate weights and puts them in a uniform (atom,icart) basis
INPUTS
displ_red = phonon mode displacement vectors in reduced coordinated.
eigvec = eigenvectors of phonons (to turn to cartesian coord frame)
elph_ds = datastructure with gkk matrix elements
FSfullpqtofull = mapping of k + q to k
h1_mat_el_sq = matrix elements $<psi_{k+q,m} | H^{1} | psi_{k,n}>$ matrix-squared
iqptirred = index of present qpoint
phfrq_tmp = phonon frequencies
qpt_irred = array of qpoint coordinates
OUTPUT
elph_ds%gkq filled
qdata(elph_ds%nbranch,elph_ds%nsppol,3) = array containing the phonon frequency, the linewidth and $\lambda_{q,\nu}$.
SOURCE
1393 subroutine normsq_gkq(displ_red,eigvec,elph_ds,FSfullpqtofull,& 1394 & h1_mat_el_sq,iqptirred,phfrq_tmp,qpt_irred,qdata) 1395 1396 !Arguments ------------------------------------ 1397 !scalars 1398 integer,intent(in) :: iqptirred 1399 type(elph_type),intent(inout) :: elph_ds 1400 !arrays 1401 integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full) 1402 real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch) 1403 real(dp),intent(in) :: eigvec(2,elph_ds%nbranch,elph_ds%nbranch) 1404 real(dp),intent(inout) :: & 1405 & h1_mat_el_sq(2,elph_ds%nFSband*elph_ds%nFSband,elph_ds%nbranch*elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol) 1406 real(dp),intent(in) :: phfrq_tmp(elph_ds%nbranch),qpt_irred(3,elph_ds%nqptirred) 1407 real(dp),intent(out) :: qdata(elph_ds%nbranch,elph_ds%nsppol,3) 1408 1409 !Local variables------------------------------- 1410 !scalars 1411 integer :: i1,i2,ier,ii,isppol,jbranch,comm 1412 real(dp) :: lambda_tot 1413 character(len=500) :: message 1414 !arrays 1415 real(dp) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol) 1416 real(dp) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol) 1417 real(dp) :: gam_now2(2,elph_ds%nbranch,elph_ds%nbranch) 1418 real(dp) :: lambda(elph_ds%nsppol) 1419 real(dp),allocatable :: matrx(:,:),val(:),vec(:,:,:) 1420 real(dp),allocatable :: zhpev1(:,:),zhpev2(:) 1421 1422 ! ************************************************************************* 1423 1424 DBG_ENTER("COLL") 1425 1426 accum_mat = zero 1427 accum_mat2 = zero 1428 comm = xmpi_world 1429 1430 if (elph_ds%ep_scalprod == 1) then 1431 ! 1432 if (elph_ds%ep_keepbands == 0) then 1433 call wrtout(std_out,' normsq_gkq : calling nmsq_gam_sumFS',"COLL") 1434 call nmsq_gam_sumFS (accum_mat,accum_mat2,displ_red,eigvec,elph_ds,FSfullpqtofull,& 1435 & h1_mat_el_sq,iqptirred) 1436 1437 else if (elph_ds%ep_keepbands == 1) then 1438 call wrtout(std_out,' normsq_gkq : calling nmsq_gam',"COLL") 1439 call nmsq_gam (accum_mat,accum_mat2,displ_red,eigvec,elph_ds,FSfullpqtofull,& 1440 & h1_mat_el_sq,iqptirred) 1441 1442 else 1443 write (message,'(a,i0)')' Wrong value for elph_ds%ep_keepbands = ',elph_ds%ep_keepbands 1444 ABI_BUG(message) 1445 end if 1446 ! 1447 else if (elph_ds%ep_scalprod == 0) then ! Interpolate on the pure "matrix of matrix elements" and do the scalar products later. 1448 ! 1449 if (elph_ds%ep_keepbands == 0) then 1450 call wrtout(std_out,' normsq_gkq : calling nmsq_pure_gkk_sumFS',"COLL") 1451 call nmsq_pure_gkk_sumFS (accum_mat,accum_mat2,displ_red,elph_ds,FSfullpqtofull,& 1452 & h1_mat_el_sq,iqptirred) 1453 1454 else if (elph_ds%ep_keepbands == 1) then 1455 call wrtout(std_out,' normsq_gkq : calling nmsq_pure_gkk',"COLL") 1456 1457 call nmsq_pure_gkk (accum_mat,accum_mat2,displ_red,elph_ds,FSfullpqtofull,& 1458 & h1_mat_el_sq,iqptirred) 1459 else 1460 write (message,'(a,i0)')' Wrong value for elph_ds%ep_keepbands = ',elph_ds%ep_keepbands 1461 ABI_BUG(message) 1462 end if 1463 ! 1464 else 1465 write (message,'(a,i0)')' Wrong value for elph_ds%ep_scalprod = ',elph_ds%ep_scalprod 1466 ABI_BUG(message) 1467 end if 1468 !end if flag for doing scalar product now. 1469 1470 1471 !MG: values without the good prefactor 1472 accum_mat = accum_mat * elph_ds%occ_factor/elph_ds%k_phon%nkpt 1473 1474 !MG: accum_mat2 contains the line-widhts before the Fourier interpolation 1475 accum_mat2 = accum_mat2 * elph_ds%occ_factor/elph_ds%k_phon%nkpt 1476 1477 !mpi sum over procs for accum_mat2 1478 call xmpi_sum (accum_mat, comm, ier) 1479 call xmpi_sum (accum_mat2, comm, ier) 1480 1481 !MG20060531i 1482 !write e-ph quantities before Fourier interpolation 1483 !save e-ph values in the temporary array qdata that will be copied into elph_ds%qgrid_data 1484 1485 write (message,'(4a,3es16.6,63a)')ch10, & 1486 & ' Phonon linewidths before interpolation ',ch10, & 1487 & ' Q point = ',qpt_irred(:,iqptirred),ch10,('=',ii=1,60),ch10,& 1488 & ' Mode Frequency (Ha) Linewidth (Ha) Lambda ' 1489 call wrtout(std_out,message,'COLL') 1490 1491 lambda_tot = zero 1492 do isppol=1,elph_ds%nsppol 1493 do ii=1,elph_ds%nbranch 1494 lambda(isppol)=zero 1495 ! MG: the tolerance factor is somehow arbitrary 1496 if (abs(phfrq_tmp(ii)) > tol10) lambda(isppol)=accum_mat2(1,ii,ii,isppol)/& 1497 & (pi*elph_ds%n0(isppol)*phfrq_tmp(ii)**2) 1498 lambda_tot=lambda_tot+lambda(isppol) 1499 write(message,'(i8,es20.6,2es16.6)' )ii,phfrq_tmp(ii),accum_mat2(1,ii,ii,isppol),lambda(isppol) 1500 call wrtout(std_out,message,'COLL') 1501 ! save values 1502 qdata(ii,isppol,1)=phfrq_tmp(ii) 1503 qdata(ii,isppol,2)=accum_mat2(1,ii,ii,isppol) 1504 qdata(ii,isppol,3)=lambda(isppol) 1505 end do !loop over branch 1506 end do !loop over sppol 1507 1508 !normalize for number of spins 1509 lambda_tot = lambda_tot / elph_ds%nsppol 1510 1511 write(message,'(61a,44x,es16.6,62a)' )('=',ii=1,60),ch10,lambda_tot,ch10,('=',ii=1,60),ch10 1512 call wrtout(std_out,message,'COLL') 1513 !ENDMG20060531 1514 1515 !immediately calculate linewidths: 1516 write(std_out,*) 'summed accum_mat = ' 1517 write(std_out,'(3(2E18.6,1x))') accum_mat(:,:,:,1) 1518 write(std_out,*) 'summed accum_mat2 = ' 1519 write(std_out,'(3(2E18.6,1x))') (accum_mat2(:,ii,ii,1),ii=1,elph_ds%nbranch) 1520 write(std_out,*) 'displ_red = ' 1521 write(std_out,'(3(2E18.6,1x))') displ_red 1522 1523 if (elph_ds%ep_scalprod == 1) then 1524 do isppol=1,elph_ds%nsppol 1525 ! Diagonalize gamma matrix at qpoint (complex matrix). Copied from dfpt_phfrq 1526 ier=0 1527 ii=1 1528 ABI_MALLOC(matrx,(2,(elph_ds%nbranch*(elph_ds%nbranch+1))/2)) 1529 do i2=1,elph_ds%nbranch 1530 do i1=1,i2 1531 matrx(1,ii)=accum_mat2(1,i1,i2,isppol) 1532 matrx(2,ii)=accum_mat2(2,i1,i2,isppol) 1533 ii=ii+1 1534 end do 1535 end do 1536 ABI_MALLOC(zhpev1,(2,2*elph_ds%nbranch-1)) 1537 ABI_MALLOC(zhpev2,(3*elph_ds%nbranch-2)) 1538 ABI_MALLOC(val,(elph_ds%nbranch)) 1539 ABI_MALLOC(vec,(2,elph_ds%nbranch,elph_ds%nbranch)) 1540 call ZHPEV ('V','U',elph_ds%nbranch,matrx,val,vec,elph_ds%nbranch,zhpev1,zhpev2,ier) 1541 1542 write (std_out,*) ' normsq_gkq : accumulated eigenvalues isppol ',isppol, ' = ' 1543 write (std_out,'(3E18.6)') val 1544 ABI_FREE(matrx) 1545 ABI_FREE(zhpev1) 1546 ABI_FREE(zhpev2) 1547 ABI_FREE(vec) 1548 ABI_FREE(val) 1549 end do ! isppol 1550 1551 else if (elph_ds%ep_scalprod == 0) then 1552 1553 1554 do isppol=1,elph_ds%nsppol 1555 call gam_mult_displ(elph_ds%nbranch, displ_red, accum_mat(:,:,:,isppol), gam_now2) 1556 1557 write (std_out,*) ' normsq_gkq : accumulated eigenvalues isppol ', isppol, ' = ' 1558 write (std_out,'(3(E14.6,1x))') (gam_now2(1,jbranch,jbranch), jbranch=1,elph_ds%nbranch) 1559 write (std_out,*) ' normsq_gkq : imag part = ' 1560 write (std_out,'(3(E14.6,1x))') (gam_now2(2,jbranch,jbranch), jbranch=1,elph_ds%nbranch) 1561 end do ! isppol 1562 1563 end if 1564 1565 DBG_EXIT("COLL") 1566 1567 end subroutine normsq_gkq
m_iogkk/completeperts [ Functions ]
[ Top ] [ m_iogkk ] [ Functions ]
NAME
completeperts
FUNCTION
Complete perturbations wrt atoms and reduced directions for a fixed qpoint. Normally there is a test in read_gkk which guarantees that enough irreducible perturbations are present to generate everything. h1_mat_el is first squared, making a (ipert,jpert) matrix which has the same symmetry properties as the dynamical matrix.
INPUTS
Cryst<crystal_t>=Info on the unit cell and symmetries. nbranch=Number of phonon branches. nFSband=Number of bands in H1 matrix elements. nkpt=Number of k-points in matrix elements. nsppol=Number of independent spin polarizations. gkk_flag = flags for presence of gkk matrix elements h1_mat_el = irreducible matrix elements to be completed and squared qpt = qpoint symq = flags for symmetry elements conserving the present qpoint tnons = translation vectors associated with symops
OUTPUT
h1_mat_el_sq = irreducible matrix elements squared and completed gkk_flag = changed on output
SOURCE
1251 subroutine completeperts(Cryst,nbranch,nFSband,nkpt,nsppol,gkk_flag,h1_mat_el,h1_mat_el_sq,& 1252 & qpt,symq,qtimrev) 1253 1254 !Arguments ------------------------------------ 1255 !scalars 1256 integer,intent(in) :: qtimrev,nbranch,nFSband,nkpt,nsppol 1257 type(crystal_t),intent(in) :: Cryst 1258 !arrays 1259 integer,intent(in) :: symq(4,2,Cryst%nsym) 1260 integer,intent(inout) :: gkk_flag(nbranch,nbranch,nkpt,nsppol) 1261 real(dp),intent(in) :: qpt(3) 1262 real(dp),intent(in) :: h1_mat_el(2,nFSband**2,nbranch,nkpt,nsppol) 1263 real(dp),intent(out) :: h1_mat_el_sq(2,nFSband**2,nbranch**2,nkpt,nsppol) 1264 1265 !Local variables------------------------------- 1266 !scalars 1267 integer :: ikpt_phon,iatom1,iatom2,ibb,idir1,idir2,ipert1,ipert2,isppol,mpert,natom 1268 real(dp) :: im1,im2,re1,re2,res 1269 character(len=500) :: msg 1270 !arrays 1271 integer,allocatable :: tmpflg(:,:,:,:) 1272 real(dp),allocatable :: tmpval(:,:,:,:,:) 1273 1274 ! ************************************************************************* 1275 1276 !WARNING! Stupid patch in d2sym3 imposes these matrices to have size natom+2 1277 natom = Cryst%natom 1278 mpert = natom+2 1279 1280 ABI_MALLOC(tmpflg,(3,mpert,3,mpert)) 1281 ABI_MALLOC(tmpval,(2,3,mpert,3,mpert)) 1282 1283 h1_mat_el_sq = zero 1284 1285 write (std_out,*) ' completeperts: shape(h1_mat_el_sq) = ', shape(h1_mat_el_sq) 1286 1287 do isppol=1,nsppol 1288 write(std_out,*)'completeperts: isppol = ', isppol 1289 ! 1290 do ikpt_phon=1,nkpt 1291 do ibb=1,nFSband**2 1292 ! 1293 tmpval = zero 1294 tmpflg = 0 1295 ! for a fixed k (q) band and sppol construct the gamma matrix for (3 natom)^2 perturbation pairs 1296 do iatom1=1,natom 1297 do idir1=1,3 1298 ipert1 = (iatom1-1)*3+idir1 1299 if (gkk_flag(ipert1,ipert1,ikpt_phon,isppol) < 0) cycle 1300 re1 = h1_mat_el(1,ibb,ipert1,ikpt_phon,isppol) 1301 im1 = h1_mat_el(2,ibb,ipert1,ikpt_phon,isppol) 1302 1303 do iatom2=1,natom 1304 do idir2=1,3 1305 ipert2 = (iatom2-1)*3+idir2 1306 if (gkk_flag(ipert2,ipert2,ikpt_phon,isppol) < 0) cycle 1307 tmpflg(idir1,iatom1,idir2,iatom2) = 1 1308 re2 = h1_mat_el(1,ibb,ipert2,ikpt_phon,isppol) 1309 im2 = h1_mat_el(2,ibb,ipert2,ikpt_phon,isppol) 1310 ! 1311 ! conjg(h1_mat_el_2) * h1_mat_el_1 1312 res = re1*re2 + im1*im2 1313 tmpval(1,idir1,iatom1,idir2,iatom2) = res 1314 res = re1*im2 - im1*re2 1315 tmpval(2,idir1,iatom1,idir2,iatom2) = res 1316 1317 end do !idir2 1318 end do !iatom2 1319 end do !idir1 1320 end do !iatom1 1321 1322 ! matrix is symmetrized like a dynamical matrix. No change of band or k 1323 ! in here. This should be checked (if we have to restrict further the symmetry operations) 1324 call d2sym3(tmpflg,tmpval,Cryst%indsym,mpert,natom,Cryst%nsym,qpt,symq,Cryst%symrec,Cryst%symrel,qtimrev,1) 1325 if (sum(tmpflg(:,1:natom,:,1:natom)) /= 3*natom*3*natom) then 1326 write(msg,'(3a,4i0)')& 1327 & 'A perturbation is missing after completion with d2sym3',ch10,& 1328 & 'tmpflg, ikpt_phon, isppol: ',tmpflg,ikpt_phon,isppol 1329 ABI_ERROR(msg) 1330 end if 1331 ! 1332 ! Save values for calculation of |gkk|^2 1333 do iatom1=1,natom 1334 do idir1=1,3 1335 ipert1 = (iatom1-1)*3+idir1 1336 do iatom2=1,natom 1337 do idir2=1,3 1338 ! 1339 ! mjv 29/10/2007 ipert2 now contains the composite index ip1*nperts+ip2 1340 ipert2 = (iatom2-1)*3 + idir2 + (ipert1-1)*3*natom 1341 h1_mat_el_sq(1,ibb,ipert2,ikpt_phon,isppol) = pi*tmpval(1,idir2,iatom2,idir1,iatom1) 1342 h1_mat_el_sq(2,ibb,ipert2,ikpt_phon,isppol) = pi*tmpval(2,idir2,iatom2,idir1,iatom1) 1343 end do 1344 end do 1345 end do 1346 end do 1347 ! 1348 end do !end ibb band dos 1349 ! 1350 ! Set flags. 1351 do ipert1=1,3*natom 1352 do ipert2=1,3*natom 1353 if (gkk_flag(ipert2,ipert1,ikpt_phon,isppol) < 0) gkk_flag(ipert2,ipert1,ikpt_phon,isppol) = 1 1354 end do 1355 end do 1356 1357 end do !end kpt_phon do 1358 end do !end sppol do 1359 1360 ABI_FREE(tmpflg) 1361 ABI_FREE(tmpval) 1362 1363 end subroutine completeperts
m_iogkk/inpgkk [ Functions ]
[ Top ] [ m_iogkk ] [ Functions ]
NAME
inpgkk
FUNCTION
read in gkk file and return eigenvalue matrix Only works for a single gkk matrix (1 perturbation and qpoint) in the file like the files produced by outgkk
INPUTS
filegkk= filename
OUTPUT
eigen1 = response function 1st order eigenvalue matrix
SOURCE
1139 subroutine inpgkk(eigen1,filegkk,hdr1) 1140 1141 !Arguments ------------------------------------ 1142 !scalars 1143 character(len=fnlen),intent(in) :: filegkk 1144 type(hdr_type), intent(out) :: hdr1 1145 !arrays 1146 real(dp),allocatable,intent(out) :: eigen1(:) 1147 1148 !Local variables------------------------------- 1149 !scalars 1150 integer :: bantot1 1151 integer :: isppol, ikpt, mband, ikb 1152 integer :: unitgkk, fform, ierr, n1wf, i1wf 1153 type(hdr_type) :: hdr0 1154 real(dp), allocatable :: eigen(:) 1155 character(len=500) :: message 1156 1157 ! ************************************************************************* 1158 1159 if (open_file(filegkk,message,newunit=unitgkk,form='unformatted',status='old') /= 0) then 1160 ABI_ERROR(message) 1161 end if 1162 1163 1164 !read in header of GS file and eigenvalues 1165 call hdr_fort_read(hdr0, unitgkk, fform) 1166 ABI_CHECK(fform /= 0, "hdr_fort_read returned fform == 0") 1167 1168 mband = maxval(hdr0%nband(:)) 1169 ABI_MALLOC(eigen,(mband)) 1170 call wrtout(std_out,'inpgkk : try to reread GS eigenvalues','COLL') 1171 1172 do isppol=1,hdr0%nsppol 1173 do ikpt=1,hdr0%nkpt 1174 read (unitgkk,IOSTAT=ierr) eigen(1:hdr0%nband(ikpt)) 1175 ABI_CHECK(ierr==0,'reading eigen from gkk file') 1176 end do 1177 end do 1178 1179 read(unitgkk,IOSTAT=ierr) n1wf 1180 ABI_CHECK(ierr==0,"reading n1wf from gkk file") 1181 1182 ABI_FREE(eigen) 1183 call hdr0%free() 1184 1185 if (n1wf > 1) then 1186 write(message,'(3a)')& 1187 & 'several 1wf records were found in the file,',ch10, & 1188 & 'which is not allowed for reading with this routine' 1189 ABI_ERROR(message) 1190 end if 1191 1192 !read in header of 1WF file 1193 call hdr_fort_read(hdr1, unitgkk, fform) 1194 if (fform == 0) then 1195 write(message,'(a,i0,a)')' 1WF header number ',i1wf,' was mis-read. fform == 0' 1196 ABI_ERROR(message) 1197 end if 1198 1199 bantot1 = 2*hdr1%nsppol*hdr1%nkpt*mband**2 1200 ABI_MALLOC(eigen1, (bantot1)) 1201 1202 1203 !retrieve 1WF <psi_k+q | H | psi_k> from gkk file and echo to output 1204 ikb = 0 1205 do isppol=1,hdr1%nsppol 1206 do ikpt=1,hdr1%nkpt 1207 read (unitgkk,IOSTAT=ierr) eigen1(ikb+1:ikb+2*hdr1%nband(ikpt)**2) 1208 ikb = ikb + 2*hdr1%nband(ikpt)**2 1209 if (ierr /= 0) then 1210 write(message,'(a,2i0)')'reading eigen1 from gkk file, spin, kpt_idx',isppol,ikpt 1211 ABI_ERROR(message) 1212 end if 1213 end do 1214 end do 1215 1216 close(unitgkk) 1217 1218 end subroutine inpgkk
m_iogkk/outgkk [ Functions ]
[ Top ] [ m_iogkk ] [ Functions ]
NAME
outgkk
FUNCTION
output gkk file for one perturbation (used for elphon calculations in anaddb)
INPUTS
bantot0 = total number of bands for all kpoints bantot1 = total number of matrix elements for 1st order eigenvalues eigen0 = GS eigenvalues eigen1 = response function 1st order eigenvalue matrix hdr0 = GS header hdr1 = RF header mpi_enreg=information about MPI parallelization
SOURCE
700 subroutine outgkk(bantot0,bantot1,outfile,eigen0,eigen1,hdr0,hdr1,mpi_enreg,phasecg) 701 702 !Arguments ------------------------------------ 703 !scalars 704 integer,intent(in) :: bantot0,bantot1 705 character(len=fnlen),intent(in) :: outfile 706 type(MPI_type),intent(in) :: mpi_enreg 707 type(hdr_type),intent(inout) :: hdr0,hdr1 708 !arrays 709 real(dp),intent(in) :: eigen0(bantot0),eigen1(2*bantot1) 710 real(dp),intent(in) :: phasecg(2,bantot1) 711 712 !Local variables------------------------------- 713 !scalars 714 integer :: fform,iband,ikpt,isppol,me,ntot,unitout 715 integer :: iband_off, mband, ierr 716 character(len=500) :: msg 717 real(dp), allocatable :: tmpeig(:) 718 719 ! ************************************************************************* 720 721 !only master should be writing to disk 722 !Init me 723 me=mpi_enreg%me_kpt 724 if (me /= 0) return 725 726 call wrtout(std_out,' writing gkk file: '//outfile,"COLL") 727 728 !initializations 729 fform = 42 730 ntot = 1 731 732 !open gkk file 733 if (open_file(outfile, msg, newunit=unitout, form='unformatted', status='unknown', action="write") /= 0) then 734 ABI_ERROR(msg) 735 end if 736 737 !output GS header 738 call hdr0%fort_write(unitout, fform, ierr) 739 ABI_CHECK(ierr == 0 , "hdr_fort_write returned ierr != 0") 740 741 !output GS eigenvalues 742 iband=0 743 do isppol=1,hdr0%nsppol 744 do ikpt=1,hdr0%nkpt 745 write (unitout) eigen0(iband+1:iband+hdr0%nband(ikpt)) 746 iband=iband+hdr0%nband(ikpt) 747 end do 748 end do 749 750 !output number of gkk in this file (1) 751 write (unitout) ntot 752 753 !output RF header 754 call hdr1%fort_write(unitout, fform, ierr) 755 ABI_CHECK(ierr == 0 , "hdr_fort_write returned ierr != 0") 756 757 !output RF eigenvalues 758 mband = maxval(hdr1%nband(:)) 759 ABI_MALLOC(tmpeig,(2*mband**2)) 760 iband_off = 0 761 tmpeig(1) = phasecg(1, 1) 762 do isppol = 1, hdr1%nsppol 763 do ikpt = 1, hdr1%nkpt 764 tmpeig = zero 765 do iband = 1, hdr1%nband(ikpt)**2 766 tmpeig (2*(iband-1)+1) = eigen1(2*(iband_off+iband-1)+1) 767 tmpeig (2*(iband-1)+2) = eigen1(2*(iband_off+iband-1)+2) 768 end do 769 write (unitout) tmpeig(1:2*hdr1%nband(ikpt)**2) 770 iband_off = iband_off + hdr1%nband(ikpt)**2 771 end do 772 end do 773 ABI_FREE(tmpeig) 774 775 !close gkk file 776 close (unitout) 777 778 end subroutine outgkk
m_iogkk/prt_gkk_yambo [ Functions ]
[ Top ] [ m_iogkk ] [ Functions ]
NAME
prt_gkk_yambo
FUNCTION
This routine outputs el-phon related quantities for the yambo code at 1 q-point
INPUTS
displ_cart = phonon displacement vectors for this q-point in Cartesian coordinates. displ_red = phonon displacement vectors for this q-point, in reduced coordinates elph_ds = datastructure containing elphon matrix elements h1_mat_el = matrix elements of first order hamiltonian for present q-point, all perturbations iqptfull = index of present q-point in full array of q-points irredpert = index of irreducible perturbation (atom displaced) natom = number of atoms phfrq = phonon frequencies at present q-point qptn = q-point we will print for
OUTPUT
only writes to a file
NOTES
SOURCE
809 subroutine prt_gkk_yambo(displ_cart,displ_red,kpt_phon,h1_mat_el,iqpt,& 810 & natom,nFSband,nkpt_phon,phfrq,qptn) 811 812 !Arguments ------------------------------------ 813 !scalars 814 integer,intent(in) :: natom,iqpt 815 integer,intent(in) :: nFSband,nkpt_phon 816 !arrays 817 real(dp),intent(in) :: kpt_phon(3,nkpt_phon) 818 real(dp),intent(in) :: h1_mat_el(2,nFSband*nFSband,3*natom,nkpt_phon,1) 819 real(dp),intent(in) :: phfrq(3*natom) 820 real(dp),intent(in) :: displ_cart(2,3*natom,3*natom) 821 real(dp),intent(in) :: displ_red(2,3*natom,3*natom) 822 real(dp),intent(in) :: qptn(3) 823 824 !Local variables------------------------------- 825 !scalars 826 integer, save :: firsttime=1 827 integer :: outunit,ikpt,imode,iband,ibandp,iatom,idir,ibandindex 828 integer :: jmode, outunit2, outunit3 829 !arrays 830 real(dp) :: gkk_mode_dep(2) 831 ! ************************************************************************* 832 833 !if first time round: 834 if (firsttime==1) then 835 firsttime=0 836 ! squash file 837 outunit=get_unit() 838 open (unit=outunit,file="yambo_elphon_data",status="REPLACE") 839 outunit2=get_unit() 840 open (unit=outunit2,file="yambo_elphon_gkk_bymode",status="replace") 841 outunit3=get_unit() 842 open (unit=outunit3,file="yambo_elphon_gkksqtw_bymode",status="replace") 843 844 ! write dimensions 845 write (outunit,'(a,I6)') 'number of el atoms ', natom 846 write (outunit2,'(a,I6)') 'number of el atoms ', natom 847 write (outunit3,'(a,I6)') 'number of el atoms ', natom 848 write (outunit,'(a,I6)') 'number of ph modes ', 3*natom 849 write (outunit2,'(a,I6)') 'number of ph modes ', 3*natom 850 write (outunit3,'(a,I6)') 'number of ph modes ', 3*natom 851 write (outunit,'(a,I6)') 'number of el bands ', nFSband 852 write (outunit2,'(a,I6)') 'number of el bands ', nFSband 853 write (outunit3,'(a,I6)') 'number of el bands ', nFSband 854 855 ! write k-points 856 write (outunit,'(a,I6)') 'number of k-points ', nkpt_phon 857 write (outunit2,'(a,I6)') 'number of k-points ', nkpt_phon 858 write (outunit3,'(a,I6)') 'number of k-points ', nkpt_phon 859 do ikpt=1,nkpt_phon 860 write (outunit,'(a,I6,3E20.10)') 'reduced coord kpoint no ', ikpt, kpt_phon(:,ikpt) 861 write (outunit2,'(a,I6,3E20.10)') 'reduced coord kpoint no ', ikpt, kpt_phon(:,ikpt) 862 write (outunit3,'(a,I6,3E20.10)') 'reduced coord kpoint no ', ikpt, kpt_phon(:,ikpt) 863 end do 864 865 ! band energies are not accessible this deep in the code: simpler to get them 866 ! from elsewhere 867 868 close (outunit) 869 close (outunit2) 870 close (outunit3) 871 end if ! first time round 872 873 !open file 874 outunit=get_unit() 875 open (unit=outunit,file="yambo_elphon_data",status="unknown",position="append") 876 877 !qpoint 878 write (outunit,'(a,I6,3E20.10)') 'reduced coord qpoint no ', iqpt, qptn(:) 879 880 !frequencies 881 do imode=1,3*natom 882 write (outunit,'(a,I6,3E20.10)') 'phonon freq no ', imode, phfrq(imode) 883 end do 884 885 !displacement vector 886 do imode=1,3*natom 887 write (outunit,'(a,I6,3E20.10)') 'phonon displ vec no ', imode 888 do iatom=1,natom 889 write (outunit,'(3(2E20.10,2x))') displ_cart(:,(iatom-1)*3+1:iatom*3,imode) 890 end do 891 end do 892 893 !the beef: matrix elements of the first order hamiltonian for displacement of 894 !all atoms along all reduced directions 895 write (outunit,'(a)') ' matrix elements of all perturbations for this q-point' 896 do ikpt=1,nkpt_phon 897 write (outunit,'(a,I6)') ' kpoint ', ikpt 898 imode=0 899 do iatom=1,natom 900 do idir=1,3 901 imode=imode+1 902 write (outunit,'(a,I6,I6)') ' atom, direction = ', iatom,idir 903 ibandindex=0 904 do iband=1,nFSband 905 do ibandp=1,nFSband 906 ibandindex=ibandindex+1 907 write (outunit,'(a,I6,I6,2E20.10)') ' mat el for n,np ', iband,ibandp,& 908 & h1_mat_el(:,ibandindex,imode,ikpt,1) 909 end do !bandp 910 end do !band 911 end do !dir 912 end do !atom 913 end do 914 915 !blank line 916 write (outunit,*) 917 close (outunit) 918 919 outunit2=get_unit() 920 open (unit=outunit2,file="yambo_elphon_gkk_bymode",status="unknown",position="append") 921 outunit3=get_unit() 922 open (unit=outunit3,file="yambo_elphon_gkksqtw_bymode",status="unknown",position="append") 923 924 !qpoint 925 write (outunit2,'(a,I6,3E20.10)') 'reduced coord qpoint no ', iqpt, qptn(:) 926 write (outunit3,'(a,I6,3E20.10)') 'reduced coord qpoint no ', iqpt, qptn(:) 927 928 !print out mode-dependent matrix elements 929 write (outunit2,'(a)') ' matrix elements of all phonon modes for this q-point' 930 write (outunit3,'(a)') ' 1/w**1/2 times matrix elements of all phonon modes for this q-point' 931 do ikpt=1,nkpt_phon 932 write (outunit2,'(a,I6)') ' kpoint ', ikpt 933 write (outunit3,'(a,I6)') ' kpoint ', ikpt 934 ibandindex=0 935 do iband=1,nFSband 936 do ibandp=1,nFSband 937 ibandindex=ibandindex+1 938 write (outunit2,'(a,I6,I6)') ' el bands n,np ', iband,ibandp 939 write (outunit3,'(a,I6,I6)') ' el bands n,np ', iband,ibandp 940 do imode=1,3*natom 941 ! gkk_mode_dep = cg_zdotc(3*natom,displ_red(:,:,imode),h1_mat_el(:,ibandindex,:,ikpt,1)) 942 gkk_mode_dep = zero 943 do jmode=1,3*natom 944 gkk_mode_dep(1) = gkk_mode_dep(1) & 945 & + displ_red(1,jmode,imode)*h1_mat_el(1,ibandindex,jmode,ikpt,1) & 946 & + displ_red(2,jmode,imode)*h1_mat_el(2,ibandindex,jmode,ikpt,1) 947 gkk_mode_dep(2) = gkk_mode_dep(2) & 948 & + displ_red(1,jmode,imode)*h1_mat_el(2,ibandindex,jmode,ikpt,1) & 949 & - displ_red(2,jmode,imode)*h1_mat_el(1,ibandindex,jmode,ikpt,1) 950 end do 951 write (outunit2,'(a,I6,2E20.10)') ' mat el for phonon mode num = ', imode, gkk_mode_dep 952 write (outunit3,'(a,I6,2E20.10)') ' 1/w**1/2 * mat el for phonon mode num = ', & 953 & imode, gkk_mode_dep/sqrt(two*abs(phfrq(imode))+tol10) 954 end do !imode 955 end do !bandp 956 end do !band 957 end do 958 !blank line 959 write (outunit2,*) 960 write (outunit3,*) 961 962 close (outunit2) 963 close (outunit3) 964 965 end subroutine prt_gkk_yambo
m_iogkk/read_el_veloc [ Functions ]
[ Top ] [ m_iogkk ] [ Functions ]
NAME
read_el_veloc
FUNCTION
This routine reads the velocities of the electronic GS for all kpts and bands then maps them into the FS kpt states
COPYRIGHT
Copyright (C) 2002-2024 ABINIT group (JPCroc) based on conducti This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
nkpt_in = number of kpoints according to parent routine nband_in = number of bands according to parent routine nsppol_in = number of spin polarizations
OUTPUT
el_veloc(nkpt_in,nband_in,3)
SOURCE
995 subroutine read_el_veloc(nband_in,nkpt_in,kpt_in,nsppol_in,elph_tr_ds) 996 997 !Arguments ----------------------------------- 998 !scalars 999 integer, intent(in) :: nband_in,nkpt_in,nsppol_in 1000 type(elph_tr_type), intent(inout) :: elph_tr_ds 1001 real(dp), intent(in) :: kpt_in(3,nkpt_in) 1002 1003 !Local variables------------------------------- 1004 !scalars 1005 integer :: bd2tot_index 1006 integer :: iband,ii,ikpt, ikpt_ddk 1007 integer :: isppol,l1,mband 1008 integer :: bantot1 1009 integer :: unit_ddk 1010 integer :: symrankkpt 1011 character(len=fnlen) :: filnam1,filnam2,filnam3 1012 character(len=500) :: msg 1013 type(hdr_type) :: hdr1 1014 type(krank_t) :: krank 1015 !arrays 1016 real(dp) :: im_el_veloc(3) 1017 real(dp),allocatable :: eig1_k(:,:) 1018 real(dp),allocatable :: eigen11(:),eigen12(:),eigen13(:) 1019 1020 ! ********************************************************************************* 1021 1022 !Read data file name 1023 !TODO: this should be standardized and read in anaddb always, not 1024 !conditionally. Otherwise when new files are added to the anaddb files 1025 !file... Catastrophe! 1026 1027 write(std_out,*)'enter read_el_veloc ' 1028 1029 !Read data file 1030 if (open_file(elph_tr_ds%ddkfilename,msg,newunit=unit_ddk,form='formatted') /= 0) then 1031 ABI_ERROR(msg) 1032 end if 1033 1034 rewind(unit_ddk) 1035 read(unit_ddk,'(a)')filnam1 ! first ddk file 1036 read(unit_ddk,'(a)')filnam2 ! second ddk file 1037 read(unit_ddk,'(a)')filnam3 ! third ddk file 1038 close (unit_ddk) 1039 1040 bantot1 = 2*nband_in**2*nkpt_in*nsppol_in 1041 1042 call inpgkk(eigen11,filnam1,hdr1) 1043 call hdr1%free() 1044 1045 call inpgkk(eigen12,filnam2,hdr1) 1046 call hdr1%free() 1047 1048 !we use the hdr1 from the last call - should add some consistency 1049 !testing here, we are trusting users not to mix different ddk files... 1050 call inpgkk(eigen13,filnam3,hdr1) 1051 1052 !Extract info from the header 1053 if(hdr1%nsppol /= nsppol_in) then 1054 ABI_ERROR('nsspol /= input nsppol') 1055 end if 1056 1057 !Get mband, as the maximum value of nband(nkpt) 1058 mband=maxval(hdr1%nband(1:hdr1%nkpt)) 1059 if (mband /= nband_in) then 1060 ABI_ERROR('nband_in input to read_el_veloc is inconsistent with mband') 1061 end if 1062 1063 write(std_out,*) 1064 write(std_out,*) 'readings from read_el_veloc header' 1065 write(std_out,'(a,i8)') ' natom =',hdr1%natom 1066 write(std_out,'(a,3i8)') ' nkpt,nband_in,mband =',hdr1%nkpt,nband_in,mband 1067 write(std_out,'(a, f10.5,a)' ) ' ecut =',hdr1%ecut,' Ha' 1068 write(std_out,'(a,e15.5,a,e15.5,a)' )' fermie =',hdr1%fermie,' Ha ',hdr1%fermie*Ha_eV,' eV' 1069 1070 ABI_MALLOC(eig1_k,(2*nband_in**2,3)) 1071 bd2tot_index = 0 1072 elph_tr_ds%el_veloc=zero 1073 1074 !need correspondence between the DDK kpoints and the kpt_phon 1075 krank = krank_new(hdr1%nkpt, hdr1%kptns) 1076 1077 do isppol=1,nsppol_in 1078 im_el_veloc(:)=zero 1079 do ikpt=1,nkpt_in 1080 symrankkpt = krank%get_rank (kpt_in(:,ikpt)) 1081 ikpt_ddk = krank%invrank(symrankkpt) 1082 if (ikpt_ddk == -1) then 1083 write(std_out,*)'read_el_veloc ******** error in correspondence between ddk and gkk kpoint sets' 1084 write(std_out,*)' kpt sets in gkk and ddk files must agree.' 1085 ABI_ERROR("Aborting now") 1086 end if 1087 bd2tot_index=2*nband_in**2*(ikpt_ddk-1) 1088 1089 ! first derivative eigenvalues for k-point 1090 eig1_k(:,1)=eigen11(1+bd2tot_index:2*nband_in**2+bd2tot_index) 1091 eig1_k(:,2)=eigen12(1+bd2tot_index:2*nband_in**2+bd2tot_index) 1092 eig1_k(:,3)=eigen13(1+bd2tot_index:2*nband_in**2+bd2tot_index) 1093 1094 ! turn el_veloc to cartesian coordinates 1095 do iband=1,nband_in 1096 do l1=1,3 1097 do ii=1,3 1098 elph_tr_ds%el_veloc(ikpt,iband,l1,isppol)=elph_tr_ds%el_veloc(ikpt,iband,l1,isppol)+& 1099 & hdr1%rprimd(l1,ii)*eig1_k(2*iband-1+(iband-1)*2*nband_in,ii)/two_pi 1100 im_el_veloc(l1)=im_el_veloc(l1)+& 1101 & hdr1%rprimd(l1,ii)*eig1_k(2*iband+(iband-1)*2*nband_in,ii)/two_pi 1102 end do 1103 end do ! l1 1104 end do 1105 end do 1106 end do ! end isppol 1107 1108 call krank%free() 1109 ABI_FREE(eig1_k) 1110 ABI_FREE(eigen11) 1111 ABI_FREE(eigen12) 1112 ABI_FREE(eigen13) 1113 1114 call hdr1%free() 1115 1116 write(std_out,*)'out of read_el_veloc ' 1117 1118 end subroutine read_el_veloc
m_iogkk/read_gkk [ Functions ]
[ Top ] [ m_iogkk ] [ Functions ]
NAME
read_gkk
FUNCTION
This routine reads in elphon matrix elements and completes them using the appropriate symmetries
INPUTS
elph_ds = datastructure containing elphon matrix elements Cryst<crystal_t>=Info on the crystal unit cell. Ifc<ifc_type>=Object containing the interatomic force constants. FSfullpqtofull = mapping of k+q to k n1wf = number of 1WF files to be read and analyzed nband = number of bands per kpoint unitgkk = unit of GKK file for reading
OUTPUT
elph_ds = modified gkq
gkk_qpt = el-ph matrix elements for irreducible qpoints and
kpoints (as a function of the reduced symmetry for the qpoint)
gkk_flag = flag array:
-1 -> element is missing
0 -> element is from symmetric qpt (Now done in complete_gkk)
1 -> element is from symmetric pert
2 -> element is kptsym of gkk file
3 -> element was read from gkk file
SOURCE
86 subroutine read_gkk(elph_ds,Cryst,ifc,Bst,FSfullpqtofull,gkk_flag,n1wf,nband,ep_prt_yambo,unitgkk) 87 88 !Arguments ------------------------------------ 89 !scalars 90 integer,intent(in) :: n1wf,nband,unitgkk,ep_prt_yambo 91 type(crystal_t),intent(in) :: Cryst 92 type(ifc_type),intent(in) :: ifc 93 type(ebands_t),intent(in) :: Bst 94 type(elph_type),intent(inout) :: elph_ds 95 !arrays 96 integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full) 97 integer,intent(out) :: gkk_flag(elph_ds%nbranch,elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol,elph_ds%nqpt_full) 98 99 !Local variables------------------------------- 100 !scalars 101 integer :: nsppol,nbranch,nFSband,minFSband,comm,use_sym 102 integer :: fform,i1wf,ikpt_phon,iatom1,iatom2 103 integer :: ib,ib1,ib2,ibb,ibranch,idir,idir1,idir2,ierr,ii,ikpt1 104 integer :: ipert,ipert1,ipert2,iqptirred,iqptfull,isppol,isym1 105 integer :: itim1,jkpt_phon,new 106 integer :: nsym1,qtimrev,syuse 107 integer :: tdonecompl,test_flag,verify 108 integer :: nqptirred_local 109 integer :: master, me 110 integer :: symrankkpt, ikpt1_phon, ik_this_proc 111 real(dp) :: res,ss,timsign 112 character(len=500) :: msg 113 type(hdr_type) :: hdr1 114 !arrays 115 integer :: FSirrtok(3,elph_ds%k_phon%nkpt) 116 integer :: symaf1(Cryst%nsym),symq(4,2,Cryst%nsym) 117 integer :: symrc1(3,3,Cryst%nsym),symrl1(3,3,Cryst%nsym) 118 integer :: tmpflg(3,Cryst%natom+2,3,Cryst%natom+2) 119 real(dp) :: displ_cart(2,3*Cryst%natom,3*Cryst%natom) 120 real(dp) :: displ_red(2,3*Cryst%natom,3*Cryst%natom) 121 real(dp) :: eigvec(2,3*Cryst%natom,3*Cryst%natom),kpt(3),phfrq_tmp(3*Cryst%natom),redkpt(3) 122 real(dp) :: qptirred_local(3,n1wf) 123 real(dp) :: tnons1(3,Cryst%nsym) 124 real(dp),allocatable :: eigen1(:,:,:),gkk_qpt_tmp(:,:,:,:) 125 real(dp),allocatable :: h1_mat_el(:,:,:,:,:),h1_mat_el_sq(:,:,:,:,:) 126 real(dp),allocatable :: qdata(:,:,:),qdata_tmp(:,:,:,:) 127 128 ! ************************************************************************* 129 130 ABI_UNUSED(Bst%bantot) 131 132 use_sym = 1 133 nsppol = elph_ds%nsppol 134 nbranch = elph_ds%nbranch 135 nFSband = elph_ds%nFSband 136 minFSband = elph_ds%minFSband 137 138 !init values for parallelization 139 comm = xmpi_world 140 me = xmpi_comm_rank(comm) 141 master = 0 142 143 ABI_MALLOC_OR_DIE(h1_mat_el,(2, nFSband**2, nbranch, elph_ds%k_phon%my_nkpt, nsppol), ierr) 144 h1_mat_el= zero 145 146 ABI_MALLOC_OR_DIE(h1_mat_el_sq,(2, nFSband**2, nbranch**2,elph_ds%k_phon%my_nkpt, nsppol), ierr) 147 h1_mat_el_sq = zero 148 149 ABI_MALLOC(elph_ds%qirredtofull,(elph_ds%nqptirred)) 150 151 !MG array to store the e-ph quantities calculated over the input Q-grid 152 ABI_MALLOC(qdata_tmp,(elph_ds%nqptirred,nbranch,nsppol,3)) 153 qdata_tmp=zero 154 155 nqptirred_local=0 !zero number of irred q-points found 156 qptirred_local(:,:)=zero 157 158 gkk_flag = -1 159 160 if (elph_ds%gkqwrite ==0) then 161 elph_ds%gkk_qpt = zero 162 163 else if (elph_ds%gkqwrite == 1) then 164 ABI_MALLOC_OR_DIE(gkk_qpt_tmp,(2,elph_ds%ngkkband**2,nbranch**2,nsppol), ierr) 165 gkk_qpt_tmp = zero 166 do iqptirred=1,elph_ds%nqptirred*elph_ds%k_phon%nkpt 167 write (elph_ds%unitgkq,REC=iqptirred) gkk_qpt_tmp 168 end do 169 ABI_FREE(gkk_qpt_tmp) 170 171 else 172 write (msg,'(a,i0)')' Wrong values for gkqwrite = ',elph_ds%gkqwrite 173 ABI_BUG(msg) 174 end if !gkqwrite 175 176 !=========================================================== 177 !Loop over all files we have 178 !read in header for perturbation 179 !should check that all files are complete, have same header 180 !(taking into account the symmetries for the qpoint), 181 !represent the correct qpoints ... 182 !MG: this task should be performed in mrggkk 183 !=========================================================== 184 185 ABI_MALLOC(eigen1,(2,nband,nband)) 186 do i1wf=1,n1wf 187 188 if (master == me) then 189 write (msg,'(2a,i4,a,i4)')ch10,' read_gkk : reading 1WF header # ',i1wf,' /',n1wf 190 call wrtout(std_out,msg,'COLL') 191 192 ! Could check for compatibility of natom, kpt grids, ecut, qpt with DDB grid... 193 ! MG: Also this task should be done in mrggkk 194 195 call hdr_fort_read(hdr1, unitgkk, fform) 196 if (fform == 0) then 197 write (msg,'(a,i0,a)')' 1WF header number ',i1wf,' was mis-read. fform == 0' 198 ABI_ERROR(msg) 199 end if 200 201 write(msg,'(a,i4)')' read_gkk : have read 1WF header #',i1wf 202 call wrtout(std_out,msg,'COLL') 203 write (msg,'(2a,i4,a)')ch10,' read_gkk : # of kpt for this perturbation: ',hdr1%nkpt,ch10 204 call wrtout(std_out,msg,'COLL') 205 206 end if 207 208 ! broadcast data to all nodes: 209 call hdr1%bcast(master, me, comm) 210 211 ! Find qpoint in full grid 212 new=1 213 do iqptfull=1,elph_ds%nqpt_full 214 kpt(:) = hdr1%qptn(:) - elph_ds%qpt_full(:,iqptfull) 215 call wrap2_pmhalf(kpt(1),redkpt(1),res) 216 call wrap2_pmhalf(kpt(2),redkpt(2),res) 217 call wrap2_pmhalf(kpt(3),redkpt(3),res) 218 ss=redkpt(1)**2+redkpt(2)**2+redkpt(3)**2 219 if(ss < tol6) then 220 new = 0 221 exit !exit with iqptfull 222 end if 223 end do !iqptfull 224 225 if (new == 1) then 226 ! Test should be at the end: dont care if there are additional 227 ! qpts in gkk file which are not on the main grid. Ignore them. 228 write (msg,'(4a,3es16.6,2a)')ch10,& 229 & ' read_gkk : WARNING- ',ch10,& 230 & ' qpoint = ',hdr1%qptn(:),ch10,& 231 & ' not found in the input q-grid. Ignoring this point ' 232 call wrtout(ab_out,msg,'COLL') 233 call wrtout(std_out,msg,'COLL') 234 if (me == master) then 235 do isppol=1,hdr1%nsppol 236 do ikpt1=1,hdr1%nkpt 237 read(unitgkk) ((eigen1(:,ii,ib),ii=1,nband),ib=1,nband) 238 end do 239 end do 240 end if 241 242 cycle !cycle the loop on i1wf 243 end if !end if (new ==1) 244 245 246 ! Check whether other pieces of the DDB have used this qpt already 247 new=1 248 do iqptirred=1,nqptirred_local 249 kpt(:) = qptirred_local(:,iqptirred) - hdr1%qptn(:) 250 call wrap2_pmhalf(kpt(1),redkpt(1),res) 251 call wrap2_pmhalf(kpt(2),redkpt(2),res) 252 call wrap2_pmhalf(kpt(3),redkpt(3),res) 253 ss=redkpt(1)**2+redkpt(2)**2+redkpt(3)**2 254 if(ss < tol6) then 255 new=0 256 exit !MG We can use this information to avoid recalculating the dynamical matrix 257 end if !but we need to use a fixed format in GKK! 258 end do !iqptirred 259 260 if (new==1) then !we have a new valid irreducible qpoint, add it! 261 nqptirred_local = nqptirred_local+1 262 if (nqptirred_local > elph_ds%nqptirred) then 263 write (msg, '(a,a,a,i6,i6)') & 264 & 'found too many qpoints in GKK file wrt anaddb input ', ch10, & 265 & 'nqpt_anaddb nqpt_gkk = ', elph_ds%nqptirred, nqptirred_local 266 ABI_ERROR(msg) 267 end if 268 qptirred_local(:,nqptirred_local) = hdr1%qptn(:) 269 iqptirred = nqptirred_local 270 tdonecompl = 0 271 h1_mat_el = zero 272 end if 273 274 ! now iqptirred is the index of the present qpoint in the array qptirred_local 275 ! and iqptfull is the index in the full qpt_full array for future reference 276 elph_ds%qirredtofull(iqptirred) = iqptfull 277 278 write (msg,'(a,i5,a,3es16.8)')& 279 & ' read_gkk : full zone qpt number ',iqptfull,' is ',elph_ds%qpt_full(:,iqptfull) 280 call wrtout(std_out,msg,'COLL') 281 282 ! if this perturbation has already been filled (overcomplete gkk) 283 ! check only 1st kpoint and spinpol, then check others 284 verify = 0 285 if (gkk_flag(hdr1%pertcase,hdr1%pertcase,1,1,elph_ds%qirredtofull(iqptirred)) /= -1) then 286 ! 287 do isppol=1,nsppol 288 do ik_this_proc=1,elph_ds%k_phon%my_nkpt 289 if (gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) == -1) then 290 write (std_out,*)" hdr1%pertcase,ik_this_proc,iqptirred",hdr1%pertcase,ik_this_proc,iqptirred 291 ABI_ERROR('Partially filled perturbation ') 292 end if 293 end do ! ikpt_phon 294 end do ! isppol 295 ! 296 ABI_WARNING(' gkk perturbation is already filled') 297 write(std_out,*)' hdr1%pertcase,iqptirred,iqptfull = ',hdr1%pertcase,iqptirred,iqptfull,& 298 & gkk_flag(hdr1%pertcase,hdr1%pertcase,1,1,elph_ds%qirredtofull(iqptirred)) 299 verify = 1 300 write (125,*) '# matrix elements for symmetric perturbation' 301 ! Instead of reading eigen1 into void, verify == 1 checks them later on wrt values in memory 302 end if !gkk_flag 303 304 ! Examine the symmetries of the q wavevector 305 ! these will be used to complete the perturbations for other atoms and idir 306 if (ep_prt_yambo==1) then 307 ! If one wants to print GKKs along phonon modes, it mean mixing of 308 ! perturbations with differnt jauge. Symmetries must then be disable. 309 call littlegroup_q(Cryst%nsym,qptirred_local(:,iqptirred),symq,Cryst%symrec,Cryst%symafm,qtimrev,prtvol=0,use_sym=0) 310 else 311 call littlegroup_q(Cryst%nsym,qptirred_local(:,iqptirred),symq,Cryst%symrec,Cryst%symafm,qtimrev,prtvol=0) 312 end if 313 314 ! Determine dynamical matrix, phonon frequencies and displacement vector for qpoint 315 !call wrtout(std_out,' read_gkk: calling inpphon to calculate the dynamical matrix','COLL') 316 317 call ifc%fourq(cryst,qptirred_local(:,iqptirred),phfrq_tmp,displ_cart,out_eigvec=eigvec) 318 319 ! Get displacement vectors for all branches in reduced coordinates 320 ! used in scalar product with H(1)_atom,idir matrix elements 321 ! Calculate $displ_red = displ_cart \cdot gprimd$ for each phonon branch 322 323 call phdispl_cart2red(Cryst%natom,Cryst%gprimd,displ_cart,displ_red) 324 325 ! prefactors for gk+q,n\prime;k,n matrix element 326 ! COMMENT : in decaft there is a weird term in the mass factor, of M-zval(species) 327 ! dont know why. Not needed to reproduce decaft results, though... 328 ! weight is squared in evaluation of 329 ! gamma_{k,q,j} = 2 \pi omega_{q,j} sum_{nu,nu\prime} |g^{q,j}_{k+q,nu\prime; k,nu}|^2 330 ! normally cancels with the 2 \pi omega_{q,j} factor in front of the sum... 331 332 ! hdr1%pertcase = idir + (ipert-1)*3 where ipert=iatom in the interesting cases 333 idir = mod (hdr1%pertcase-1,3)+1 334 ipert = int(dble(hdr1%pertcase-idir)/three)+1 335 336 write (msg,'(4a,i3,a,i3,a,i4,a)')ch10,& 337 & ' read_gkk : calling littlegroup_pert to examine the symmetries of the full perturbation ',ch10,& 338 & ' idir = ',idir,' ipert = ',ipert,' and Q point = ',iqptirred,ch10 339 call wrtout(std_out,msg,'COLL') 340 341 ! Examine the symmetries of the full perturbation these will be used to complete the kpoints 342 ! DOESNT USE TIME REVERSAL IN littlegroup_pert except for gamma 343 344 syuse=0 345 346 call littlegroup_pert(Cryst%gprimd,idir,Cryst%indsym,ab_out,ipert,Cryst%natom,Cryst%nsym,nsym1,2,Cryst%symafm,symaf1,& 347 & symq,Cryst%symrec,Cryst%symrel,symrl1,syuse,Cryst%tnons,tnons1) 348 349 do isym1=1,nsym1 350 call mati3inv(symrl1(:,:,isym1),symrc1(:,:,isym1)) 351 end do 352 FSirrtok = 0 353 354 ! ======================================================== 355 ! Loop over irred kpts in file, and fill the default gkk 356 ! ======================================================== 357 358 ! MG NOTE : in the present implementation, if nsppol /=1 the code stops in rchkGSheader! 359 do isppol=1,hdr1%nsppol !Loop over spins is trivial? Not tested. 360 write (std_out,*) ' read_gkk : isppol = ', isppol 361 362 do ikpt1=1,hdr1%nkpt !Loop over irred kpoints, WARNING nkpt depends on qpoint and symmetry! 363 ! 364 ! this is the main read of the gkk matrix elements from the file (eigen1 arrays) 365 ! it has to be done exactly nsppol*nkpt times, and the kpt_phon are completed 366 ! where appropriate in the loop below (normally succeeding only once for each kpt) 367 ! 368 if (master == me) then 369 read(unitgkk) ((eigen1(:,ii,ib),ii=1,nband),ib=1,nband) 370 end if 371 372 ! MPI broadcast data to all nodes: 373 call xmpi_bcast(eigen1, master, comm, ierr) 374 375 ! find place of irred k in k_phon 376 ! the kpoints in the file (kptns) could be ordered arbitrarily 377 symrankkpt = elph_ds%k_phon%krank%get_rank (hdr1%kptns(:,ikpt1)-qptirred_local(:,iqptirred)) 378 ikpt1_phon = elph_ds%k_phon%krank%invrank(symrankkpt) 379 if (ikpt1_phon < 0) then 380 write (msg,'(a,3es16.6,a)')' irred k ',hdr1%kptns(:,ikpt1),' was not found in full grid' 381 ABI_ERROR(msg) 382 end if 383 ! find correspondence between this kpt_phon and the others 384 ! symrc1 conserves perturbation as well as qpoint 385 ! add to FSirrtok list 386 do isym1=1,nsym1 387 do itim1=0,qtimrev 388 timsign=one-two*itim1 389 kpt(:) = timsign*matmul(symrc1(:,:,isym1), elph_ds%k_phon%kpt(:,ikpt1_phon)) 390 391 symrankkpt = elph_ds%k_phon%krank%get_rank (kpt) 392 jkpt_phon = elph_ds%k_phon%krank%invrank(symrankkpt) 393 394 if (jkpt_phon > 0) then 395 FSirrtok(1,jkpt_phon) = ikpt1_phon 396 FSirrtok(2,jkpt_phon) = isym1 397 FSirrtok(3,jkpt_phon) = itim1 398 else 399 write (msg,'(a,3es16.6,a,i5,a,i4,a)')& 400 & ' sym equivalent of kpt ',hdr1%kptns(:,ikpt1),' by sym ',& 401 & isym1,' and itime ',itim1,' was not found' 402 ABI_ERROR(msg) 403 end if 404 end do !itim1 405 end do !isim1 406 407 408 ! 409 ! Here check if the symmetry-copied gkk at new k point is equal to the one found in the file for non-irreducible point 410 ! NB This is DEBUG code 411 ! 412 if (verify == 1 .and. elph_ds%k_phon%my_kpt(ikpt1_phon) == me) then 413 do ik_this_proc = 1, elph_ds%k_phon%my_nkpt 414 if (elph_ds%k_phon%my_ikpt(ik_this_proc) == ikpt1_phon) exit 415 end do 416 do ib1=1,nFSband 417 do ib2=1,nFSband 418 ibb = (ib1-1)*nFSband+ib2 419 write (125,'(2(2E16.6,2x))') h1_mat_el(:,ibb,hdr1%pertcase,ik_this_proc,isppol),& 420 & eigen1(:,minFSband-1+ib2,minFSband-1+ib1) 421 end do 422 end do 423 end if !verify end DEBUG code 424 425 426 do ik_this_proc = 1, elph_ds%k_phon%my_nkpt 427 ! should I be dealing with this k-point? 428 jkpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc) 429 430 ! does present ikpt1 contribute to this k-point? 431 if (FSirrtok(1,jkpt_phon) /= ikpt1_phon) cycle 432 433 ! if this kpoint has already been filled (overcomplete gkk) 434 if (gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) /= -1) then 435 ABI_WARNING("gkk element is already filled") 436 write(std_out,*)' hdr1%pertcase,ik_this_proc,isppol,iqptirred = ',& 437 & hdr1%pertcase,ik_this_proc,isppol,iqptirred,& 438 & gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) 439 ! exit 440 end if !gkk_flag 441 442 ! =============================================================== 443 ! TODO: if there is a phase factor in swapping k-points, insert it here in copy to h1_mat_el 444 ! as a function of symops in FSirrtok 445 ! complete gkk for symmetric ikpt_phon with sym1 which conserve 446 ! the full perturbation+qpoint 447 ! Not tested explicitly, but the results for Pb using reduced kpts look good 448 ! should do same RF calculation with nsym=1 and check 449 ! =============================================================== 450 451 ! save this kpoint 452 do ib1=1,nFSband 453 do ib2=1,nFSband 454 ibb = (ib1-1)*nFSband+ib2 455 456 ! real 457 res=eigen1(1,minFSband-1+ib2,minFSband-1+ib1) 458 h1_mat_el(1,ibb,hdr1%pertcase,ik_this_proc,isppol) = res 459 460 ! imag 461 res=eigen1(2,minFSband-1+ib2,minFSband-1+ib1) 462 h1_mat_el(2,ibb,hdr1%pertcase,ik_this_proc,isppol) = res 463 end do !ib2 464 end do !ib1 465 ! if jkpt is equal to ikpt1_phon (if clause above) flag == 3 466 if (FSirrtok(2,jkpt_phon) == 1) then 467 gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) = 3 468 ! if jkpt_phon comes from ikpt1_phon flag == 2 with some symop 469 else 470 gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) = 2 471 end if 472 473 end do !jkpt_phon 474 475 ! =============================================================== 476 ! we now have contribution to g(k+q,k; \kappa,\alpha) from one 477 ! kpoint,and one perturbation, 478 ! NB: each perturbation will contribute to all the modes later! 479 ! 480 ! SHOULD ONLY DO THIS FOR THE SYMS NEEDED 481 ! TO COMPLETE THE PERTURBATIONS!!! 482 ! ================================================================ 483 484 end do !ikpt1 485 end do !isppol 486 487 ! 14 Jan 2014 removed test on verify - in new scheme full BZ is read in and should be used to avoid phase errors 488 ! if (verify == 1) cycle 489 490 ! Checks on irred grid provided and on gkk_flag accumulated up to now 491 if (elph_ds%tuniformgrid == 1) then ! check if irred kpoints found reconstitute the FS kpts 492 do ikpt_phon=1,elph_ds%k_phon%nkpt 493 if (FSirrtok(1,ikpt_phon) == 0) then 494 write(msg,'(a,3es16.6,2a)')& 495 & ' kpt = ',elph_ds%k_phon%kpt(:,ikpt_phon),ch10,& 496 & ' is not the symmetric of one of those found in the GKK file' 497 ABI_ERROR(msg) 498 end if 499 end do !ikpt_phon 500 501 ! normally at this point we have used all the gkk for all kpoints on the FS 502 ! for the given irred perturbation: check 503 do ik_this_proc = 1, elph_ds%k_phon%my_nkpt 504 ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc) 505 506 if (gkk_flag(hdr1%pertcase, hdr1%pertcase, ik_this_proc, 1, elph_ds%qirredtofull(iqptirred)) == -1) then 507 write (msg,'(a,i3,a,3es18.6,2a,i3,a,i3,a,3es18.6,a,a,i4,a,a)')& 508 & ' For irreducible qpt ', iqptirred,' = ',qptirred_local(:,iqptirred),ch10, & 509 & ' the gkk element : pertcase = ',hdr1%pertcase,' ik_this_proc = ',ik_this_proc, & 510 & ' kpt = ',elph_ds%k_phon%kpt(:,ikpt_phon),ch10,& 511 & ' and isppol ',1,ch10,& 512 & ' was not found by symmetry operations on the irreducible kpoints given' 513 ABI_ERROR(msg) 514 end if 515 end do !ikpt_phon 516 end if ! end elph_ds%tuniformgrid == 1 checks 517 518 write(msg,'(a,i0)')' read_gkk : Done completing the kpoints for pertcase ',hdr1%pertcase 519 call wrtout(std_out,msg,'COLL') 520 521 tmpflg(:,:,:,:) = 0 522 523 do idir1=1,3 524 do iatom1=1,Cryst%natom 525 ipert1 = (iatom1-1)*3+idir1 526 do idir2=1,3 527 do iatom2=1,Cryst%natom 528 ipert2 = (iatom2-1)*3+idir2 529 if (gkk_flag(ipert1,ipert1,1,1,elph_ds%qirredtofull(iqptirred)) >= 0 .and. & 530 & gkk_flag(ipert2,ipert2,1,1,elph_ds%qirredtofull(iqptirred)) >= 0) then 531 tmpflg(idir1,iatom1,idir2,iatom2) = 1 532 end if 533 end do 534 end do 535 end do 536 end do 537 538 539 ! =============================================== 540 ! Full test: need all perturbations explicitly 541 ! =============================================== 542 543 test_flag = 0 544 if (sum(tmpflg(:,1:Cryst%natom,:,1:Cryst%natom)) == (3*Cryst%natom)**2 .and. tdonecompl == 0) test_flag = 1 545 546 write(std_out,*)'read_gkk: tdonecompl = ', tdonecompl 547 548 ! de-activate completion of perts by symmetry for now. 549 ! Must be called when all irreducible perturbations are in memory!!!! 550 if (test_flag == 1 .and. tdonecompl == 0) then 551 552 ! write(std_out,*) ' read_gkk : enter fxgkkphase before completeperts' 553 ! call fxgkkphase(elph_ds,gkk_flag,h1_mat_el,iqptirred) 554 555 if (ep_prt_yambo==1) then 556 if (elph_ds%k_phon%my_nkpt /= elph_ds%k_phon%nkpt) then 557 write (msg, '(a)') 'prt_gkk_yambo can not handle parallel anaddb yet' 558 ABI_ERROR(msg) 559 end if 560 call prt_gkk_yambo(displ_cart,displ_red,elph_ds%k_phon%kpt,h1_mat_el,iqptirred,& 561 & Cryst%natom,nFSband,elph_ds%k_phon%my_nkpt,phfrq_tmp,hdr1%qptn) 562 end if 563 564 ! ======================================================================== 565 ! Now use more general symops to complete the other equivalent 566 ! perturbations: the kpoints are also shuffled by these symops 567 ! afterwards h1_mat_el_sq contains gamma_\tau\alpha,\tau'\alpha' in reduced coordinates 568 ! 569 ! \gamma_{\tau'\alpha',\tau\alpha} = 570 ! <psi_{k+q,ib2}| H(1)_{\tau'\alpha'}| psi_{k,ib1}>* \cdot 571 ! <psi_{k+q,ib2}| H(1)_{\tau \alpha }| psi_{k,ib1}> 572 ! 573 ! ======================================================================== 574 575 call completeperts(Cryst,nbranch,nFSband,elph_ds%k_phon%my_nkpt,nsppol,& 576 & gkk_flag(:,:,:,:,elph_ds%qirredtofull(iqptirred)),h1_mat_el,h1_mat_el_sq,qptirred_local(:,iqptirred),symq,qtimrev) 577 578 tdonecompl = 1 579 end if 580 581 ! ============================================================== 582 ! if we have all the perturbations for this qpoint, proceed 583 ! with scalar product, norm squared, and add weight factors 584 ! 585 ! SHOULD HAVE A TEST SO h1_mat_el IS NOT OVERWRITTEN 586 ! BEFORE PREVIOUS QPOINT IS FINISHED!!!!! 587 ! ============================================================== 588 589 test_flag = 1 590 do isppol=1,nsppol 591 do ik_this_proc = 1, elph_ds%k_phon%my_nkpt 592 do ibranch=1,nbranch 593 if (gkk_flag (ibranch,ibranch,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) == -1) then 594 test_flag = 0 595 exit 596 end if 597 end do 598 end do 599 end do 600 601 if (test_flag /= 0) then 602 call wrtout(std_out,' read_gkk : enter normsq_gkq',"COLL") 603 604 ! MG temporary array to save ph-linewidths before Fourier interpolation 605 ABI_MALLOC(qdata,(nbranch,nsppol,3)) 606 qdata(:,:,:)=zero 607 608 call normsq_gkq(displ_red,eigvec,elph_ds,FSfullpqtofull,& 609 & h1_mat_el_sq,iqptirred,phfrq_tmp,qptirred_local,qdata) 610 611 ! save gkk_qpt, eventually to disk, for bands up to ngkkband, 612 ! NB: if the sum over bands has been performed ngkkband is 1 instead of nFSband 613 if (elph_ds%gkqwrite == 0) then 614 elph_ds%gkk_qpt(:,:,:,:,:,iqptirred) = h1_mat_el_sq(:,1:elph_ds%ngkkband*elph_ds%ngkkband,:,:,:) 615 else 616 ! write all kpoints to disk 617 write (std_out,*) 'size of record to be written: ', 8 * 2*elph_ds%ngkkband*elph_ds%ngkkband*& 618 & elph_ds%nbranch*elph_ds%nbranch*elph_ds%k_phon%my_nkpt*elph_ds%nsppol 619 inquire(unit=elph_ds%unitgkq, recl=isppol) 620 write (std_out,*) 'recl =', isppol 621 write (std_out,*) 'iqptirred ', iqptirred 622 do ik_this_proc = 1, elph_ds%k_phon%my_nkpt 623 write (elph_ds%unitgkq,REC=((iqptirred-1)*elph_ds%k_phon%my_nkpt+ik_this_proc)) & 624 & h1_mat_el_sq(:,1:elph_ds%ngkkband*elph_ds%ngkkband,:,ik_this_proc,:) 625 end do 626 end if 627 628 qdata_tmp(iqptirred,:,:,:)=qdata(:,:,:) 629 ABI_FREE(qdata) 630 end if 631 632 call hdr1%free() 633 634 end do !of i1wf 635 636 !got all the gkk perturbations 637 638 ABI_FREE(eigen1) 639 ABI_FREE(h1_mat_el) 640 ABI_FREE(h1_mat_el_sq) 641 642 if (nqptirred_local /= elph_ds%nqptirred) then 643 write (msg, '(3a,i0,i0)') & 644 & ' Found wrong number of qpoints in GKK file wrt anaddb input ', ch10, & 645 & ' nqpt_anaddb nqpt_gkk = ', elph_ds%nqptirred, nqptirred_local 646 ABI_ERROR(msg) 647 end if 648 649 !normally at this point we have the gkk for all kpoints on the FS 650 !for all the perturbations. Otherwise a 1WF file is missing. 651 !NOTE: still havent checked the qpoint grid completeness 652 do iqptirred=1,elph_ds%nqptirred 653 do isppol=1,nsppol 654 do ik_this_proc = 1, elph_ds%k_phon%my_nkpt 655 ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc) 656 do ipert=1,nbranch 657 if (gkk_flag(ipert,ipert,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) == -1) then 658 write (msg,'(a,i5,1x,i5,1x,i5,1x,i5,a,a)')& 659 & ' gkk element',ipert,ikpt_phon,isppol,iqptirred,' was not found by symmetry operations ',& 660 & ' on the irreducible perturbations and qpoints given' 661 ABI_ERROR(msg) 662 end if 663 end do !ipert 664 end do !ik_this_proc 665 end do !isppol 666 end do !iqptirred 667 668 call wrtout(std_out,'read_gkk : done completing the perturbations (and checked!)','COLL') 669 670 !MG save phonon frequencies, ph-linewidths and lambda(q,n) values before Fourier interpolation 671 ABI_MALLOC(elph_ds%qgrid_data,(elph_ds%nqptirred,nbranch,nsppol,3)) 672 673 do iqptirred=1,elph_ds%nqptirred 674 elph_ds%qgrid_data(iqptirred,:,:,:)=qdata_tmp(iqptirred,:,:,:) 675 end do 676 677 ABI_FREE(qdata_tmp) 678 679 end subroutine read_gkk