TABLE OF CONTENTS


ABINIT/m_iogkk [ Modules ]

[ Top ] [ Modules ]

NAME

  m_iogkk

FUNCTION

  IO routines for GKK files

COPYRIGHT

  Copyright (C) 2008-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_iogkk
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use defs_elphon
33  use m_errors
34  use m_abicore
35  use m_xmpi
36  use m_kptrank
37  use m_hdr
38 
39  use m_numeric_tools,   only : wrap2_pmhalf
40  use m_io_tools,        only : open_file, get_unit
41  use m_symtk,           only : mati3inv, littlegroup_q
42  use m_geometry,        only : phdispl_cart2red, littlegroup_pert
43  use m_crystal,         only : crystal_t
44  use m_ifc,             only : ifc_type, ifc_fourq
45  use m_dynmat,          only : d2sym3
46 
47  implicit none
48 
49  private

ABINIT/nmsq_gam [ Functions ]

[ Top ] [ 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

PARENTS

      normsq_gkq

CHILDREN

      gam_mult_displ,zgemm

SOURCE

1725 subroutine nmsq_gam (accum_mat,accum_mat2,displ_red,eigvec,elph_ds,FSfullpqtofull,&
1726 &  h1_mat_el_sq,iqptirred)
1727 
1728 
1729 !This section has been created automatically by the script Abilint (TD).
1730 !Do not modify the following lines by hand.
1731 #undef ABI_FUNC
1732 #define ABI_FUNC 'nmsq_gam'
1733 !End of the abilint section
1734 
1735  implicit none
1736 
1737 !Arguments ------------------------------------
1738 !scalars
1739  integer,intent(in) :: iqptirred
1740  type(elph_type),intent(inout) :: elph_ds
1741 !arrays
1742  integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)
1743  real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
1744  real(dp),intent(in) :: eigvec(2,elph_ds%nbranch,elph_ds%nbranch)
1745  real(dp),intent(inout) :: &
1746 & 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)
1747  real(dp),intent(inout) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1748  real(dp),intent(inout) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1749 
1750 !Local variables-------------------------------
1751 ! tmp variables for diagonalization
1752 !scalars
1753  integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ibranch,isppol,ipert1
1754  integer :: jbranch
1755  integer :: iqpt_fullbz
1756  integer :: ik_this_proc
1757  real(dp) :: sd1,sd2
1758  character(len=500) :: message
1759 !arrays
1760  real(dp) :: gkq_1band(2,elph_ds%nbranch,elph_ds%nbranch)
1761  real(dp) :: tmp_mat2(2,elph_ds%nbranch,elph_ds%nbranch)
1762  real(dp) :: zgemm_tmp_mat(2,elph_ds%nbranch,elph_ds%nbranch)
1763 
1764 ! *************************************************************************
1765 
1766  if (elph_ds%ep_keepbands == 0) then
1767    write (message,'(a,i0)')' elph_ds%ep_keepbands should be 1 while is ',elph_ds%ep_keepbands
1768    MSG_ERROR(message)
1769  end if
1770 
1771 !MG20060603 NOTE:
1772 !accum_mat and accum_mat2 are real, the imaginary part is used for debugging purpose
1773 !accum_mat2 is used to store the phonon-linewidhts before interpolation
1774 
1775  iqpt_fullbz = elph_ds%qirredtofull(iqptirred)
1776  write(std_out,*) 'nmsq_gam : iqptirred = ', iqptirred
1777 
1778  do isppol=1,elph_ds%nsppol
1779    do ik_this_proc =1, elph_ds%k_phon%my_nkpt
1780      ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
1781 
1782      ikpt_phonq = FSfullpqtofull(ikpt_phon,iqpt_fullbz)
1783 
1784      do ib1=1,elph_ds%nFSband
1785        sd1 = elph_ds%k_phon%wtk(ib1,ikpt_phon,isppol) !weights for distance from the fermi surface
1786 
1787        do ib2=1,elph_ds%nFSband
1788          sd2 = elph_ds%k_phon%wtk(ib2,ikpt_phonq,isppol) !weights for distance from the fermi surface
1789          ibeff = ib2+elph_ds%nFSband*(ib1-1)
1790 
1791          gkq_1band(:,:,:) = zero
1792 
1793          zgemm_tmp_mat= reshape (h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol),(/2,elph_ds%nbranch,elph_ds%nbranch/))
1794 
1795          call gam_mult_displ(elph_ds%nbranch, displ_red, zgemm_tmp_mat, tmp_mat2)
1796 
1797 !        sum over bands
1798          do ipert1=1,elph_ds%nbranch
1799            gkq_1band(1,ipert1,ipert1) = gkq_1band(1,ipert1,ipert1) + tmp_mat2(1,ipert1,ipert1)
1800          end do
1801 
1802 !        summing over k points and bands, still diagonal in jbranch
1803          accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_1band(:,:,:)*sd1*sd2
1804 
1805 !        MG20060603 : summing over bands and kpoints with weights to calculate the phonon linewidth
1806          do jbranch=1,elph_ds%nbranch
1807            accum_mat2(:,jbranch,jbranch,isppol) = accum_mat2(:,jbranch,jbranch,isppol) + gkq_1band(:,jbranch,jbranch)*sd1*sd2
1808          end do
1809 !        END MG
1810 
1811 
1812 !        now turn to cartesian coordinates
1813 
1814 !        Final Gamma matrix (hermitian) = E * D_g * E^{+}
1815 !        Where E^{+} is the hermitian conjugate of the eigenvector matrix E
1816 !        And D_g is the diagonal matrix of values of gamma for this qpoint
1817 
1818 !        Here gkq_1band is indexed with real phonon modes (not atom+idir)
1819 !        turn gkq_1band to atom+cartesian coordinates (instead of normal coordinates for qpoint)
1820          tmp_mat2(:,:,:) = zero
1821          do ibranch =1,elph_ds%nbranch
1822            do jbranch =1,elph_ds%nbranch
1823              tmp_mat2(1,ibranch,jbranch) = tmp_mat2(1,ibranch,jbranch) + &
1824 &             eigvec(1,ibranch,jbranch) * gkq_1band(1,jbranch,jbranch)
1825              tmp_mat2(2,ibranch,jbranch) = tmp_mat2(2,ibranch,jbranch) + &
1826 &             eigvec(2,ibranch,jbranch) * gkq_1band(1,jbranch,jbranch)
1827            end do
1828          end do
1829          gkq_1band(:,:,:) = zero
1830 
1831 !        here eigvec is transposed and complexconjugated.
1832          zgemm_tmp_mat=zero
1833          call zgemm('n','c',elph_ds%nbranch,elph_ds%nbranch,elph_ds%nbranch,cone,&
1834 &         tmp_mat2,elph_ds%nbranch,eigvec,elph_ds%nbranch,czero,zgemm_tmp_mat,elph_ds%nbranch)
1835 
1836          gkq_1band = zgemm_tmp_mat
1837 
1838 !        gamma matrix contribution in cartesian coordinates (ie interpolatable form)
1839          h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol) = reshape(gkq_1band,(/2,elph_ds%nbranch*elph_ds%nbranch/))
1840 
1841        end do
1842      end do
1843 !    END loop over bands ib1 ib2
1844 
1845    end do
1846 !  END loop over kpt_phon
1847  end do
1848 !END loop over nsppol
1849 
1850 
1851 end subroutine nmsq_gam

ABINIT/nmsq_gam_sumfs [ Functions ]

[ Top ] [ 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

PARENTS

      normsq_gkq

CHILDREN

      gam_mult_displ,zgemm

SOURCE

1885 subroutine nmsq_gam_sumFS(accum_mat,accum_mat2,displ_red,eigvec,elph_ds,FSfullpqtofull,&
1886 &   h1_mat_el_sq,iqptirred)
1887 
1888 
1889 !This section has been created automatically by the script Abilint (TD).
1890 !Do not modify the following lines by hand.
1891 #undef ABI_FUNC
1892 #define ABI_FUNC 'nmsq_gam_sumFS'
1893 !End of the abilint section
1894 
1895  implicit none
1896 
1897 !Arguments ------------------------------------
1898 !scalars
1899  integer,intent(in) :: iqptirred
1900  type(elph_type),intent(inout) :: elph_ds
1901 !arrays
1902  integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)
1903  real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
1904  real(dp),intent(in) :: eigvec(2,elph_ds%nbranch,elph_ds%nbranch)
1905  real(dp),intent(inout) :: &
1906 & 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)
1907  real(dp),intent(inout) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1908  real(dp),intent(inout) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1909 
1910 !Local variables-------------------------------
1911 !scalars
1912  integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ibranch,ipert1,isppol,jbranch,iqpt_fullbz
1913  integer :: ik_this_proc
1914  real(dp) :: sd1,sd2
1915  character(len=500) :: message
1916 !arrays
1917  real(dp) :: gkq_sum_bands(2,elph_ds%nbranch,elph_ds%nbranch)
1918  real(dp) :: tmp_gkq_sum_bands(2,elph_ds%nbranch,elph_ds%nbranch)
1919  real(dp) :: tmp_mat2(2,elph_ds%nbranch,elph_ds%nbranch)
1920  real(dp),allocatable :: zgemm_tmp_mat(:,:,:)
1921 
1922 ! *************************************************************************
1923 
1924  if (elph_ds%ep_keepbands /= 0) then
1925    write (message,'(a,i0)')' elph_ds%ep_keepbands should be 0 in order to average over bands!',elph_ds%ep_keepbands
1926    MSG_ERROR(message)
1927  end if
1928 
1929  iqpt_fullbz = elph_ds%qirredtofull(iqptirred)
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  ABI_ALLOCATE(zgemm_tmp_mat ,(2,elph_ds%nbranch,elph_ds%nbranch))
1937 
1938  do isppol=1,elph_ds%nsppol
1939    do ik_this_proc =1, elph_ds%k_phon%my_nkpt
1940      ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
1941 
1942      ikpt_phonq = FSfullpqtofull(ikpt_phon,iqpt_fullbz)
1943 
1944      gkq_sum_bands = zero
1945      tmp_gkq_sum_bands = zero
1946 
1947 
1948      do ib1=1,elph_ds%nFSband
1949 !      weights for distance from the fermi surface
1950        sd1 = elph_ds%k_phon%wtk(ib1,ikpt_phon,isppol)
1951 
1952        do ib2=1,elph_ds%nFSband
1953 !        weights for distance from the fermi surface
1954          sd2 = elph_ds%k_phon%wtk(ib2,ikpt_phonq,isppol)
1955          ibeff=ib2+(ib1-1)*elph_ds%nFSband
1956 
1957          zgemm_tmp_mat = reshape(h1_mat_el_sq(:,ibeff,:,isppol,ik_this_proc),(/2,elph_ds%nbranch,elph_ds%nbranch/))
1958 
1959          call gam_mult_displ(elph_ds%nbranch, displ_red, zgemm_tmp_mat, tmp_mat2)
1960 
1961 !        sum over bands in gkq_sum_bands
1962          do ipert1=1,elph_ds%nbranch
1963            gkq_sum_bands(1,ipert1,ipert1) = gkq_sum_bands(1,ipert1,ipert1) + sd1*sd2*tmp_mat2(1,ipert1,ipert1)
1964          end do
1965 
1966 
1967 
1968        end do
1969      end do
1970 !    END loop over bands
1971 
1972 !    summing over k points, still diagonal in jbranch
1973      accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_sum_bands(:,:,:)
1974      accum_mat2(:,:,:,isppol) = accum_mat2(:,:,:,isppol) + gkq_sum_bands(:,:,:)
1975 
1976 !    summed over bands, now turn to cartesian coordinates
1977 
1978 !    Final Gamma matrix (hermitian) = E * D_g * E^{+}
1979 !    Where E^{+} is the hermitian conjugate of the eigenvector matrix E
1980 !    And D_g is the diagonal matrix of values of gamma for this qpoint
1981 
1982 !    Here gkq_sum_bands is indexed with real phonon modes (not atom+idir)
1983 !    turn gkq_sum_bands to atom+cartesian coordinates (instead of normal coordinates for qpoint)
1984 !    This is not a full matrix multiplication, just vector one, by
1985 !    gkq_sum_bands(1,jbranch,jbranch)
1986      tmp_mat2(:,:,:) = zero
1987      do ibranch =1,elph_ds%nbranch
1988        do jbranch =1,elph_ds%nbranch
1989          tmp_mat2(1,ibranch,jbranch) = tmp_mat2(1,ibranch,jbranch) + &
1990 &         eigvec(1,ibranch,jbranch) * &
1991 &         gkq_sum_bands(1,jbranch,jbranch)
1992          tmp_mat2(2,ibranch,jbranch) = tmp_mat2(2,ibranch,jbranch) + &
1993 &         eigvec(2,ibranch,jbranch) * &
1994 &         gkq_sum_bands(1,jbranch,jbranch)
1995        end do
1996      end do
1997 
1998 !    here eigvec is transposed and complex conjugated.
1999      zgemm_tmp_mat=zero
2000      call zgemm('n','c',elph_ds%nbranch,elph_ds%nbranch,elph_ds%nbranch,cone,&
2001 &     tmp_mat2,elph_ds%nbranch,eigvec,elph_ds%nbranch,czero,zgemm_tmp_mat,elph_ds%nbranch)
2002 
2003      gkq_sum_bands = zgemm_tmp_mat
2004 
2005 !    ! gamma matrix contribution in cartesian coordinates (ie interpolatable form)
2006 !    gamma matrix contribution in reduced coordinates (ie interpolatable form)
2007      h1_mat_el_sq(:,1,:,ik_this_proc,isppol) = reshape(gkq_sum_bands(:,:,:),(/2,elph_ds%nbranch*elph_ds%nbranch/))
2008 
2009 !    accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_sum_bands(:,:,:)
2010    end do
2011 !  END loop over kpt_phon
2012  end do
2013 !END loop over sppol
2014 
2015  ABI_DEALLOCATE(zgemm_tmp_mat)
2016 
2017 end subroutine nmsq_gam_sumFS

ABINIT/nmsq_pure_gkk [ Functions ]

[ Top ] [ 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

PARENTS

      normsq_gkq

CHILDREN

      gam_mult_displ

SOURCE

2053 subroutine nmsq_pure_gkk(accum_mat,accum_mat2,displ_red,elph_ds,FSfullpqtofull,&
2054 &   h1_mat_el_sq,iqptirred)
2055 
2056 
2057 !This section has been created automatically by the script Abilint (TD).
2058 !Do not modify the following lines by hand.
2059 #undef ABI_FUNC
2060 #define ABI_FUNC 'nmsq_pure_gkk'
2061 !End of the abilint section
2062 
2063  implicit none
2064 
2065 !Arguments ------------------------------------
2066 !scalars
2067  integer,intent(in) :: iqptirred
2068  type(elph_type),intent(inout) :: elph_ds
2069 !arrays
2070  integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)
2071  real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
2072  real(dp),intent(inout) :: &
2073 & 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)
2074  real(dp),intent(inout) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
2075  real(dp),intent(inout) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
2076 
2077 !Local variables-------------------------------
2078 !scalars
2079  integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ipert1,isppol
2080  integer :: iqpt_fullbz
2081  integer :: ik_this_proc
2082  real(dp) :: sd1,sd2
2083  character(len=500) :: message
2084 !arrays
2085  real(dp) :: gkq_sum_bands(2,elph_ds%nbranch,elph_ds%nbranch)
2086  real(dp) :: tmp_mat2(2,elph_ds%nbranch,elph_ds%nbranch)
2087  real(dp) :: zgemm_tmp_mat(2,elph_ds%nbranch,elph_ds%nbranch)
2088 
2089 ! *************************************************************************
2090 
2091  if (elph_ds%ep_keepbands /= 1) then
2092    message = ' elph_ds%ep_keepbands should be 1 to keep bands!'
2093    MSG_ERROR(message)
2094  end if
2095 
2096  iqpt_fullbz = elph_ds%qirredtofull(iqptirred)
2097 
2098 !h1_mat_el_sq is already fine here - nothing to do
2099 
2100 
2101 !MG20060603 NOTE:
2102 !accum_mat and accum_mat2 are real, the imaginary part is used for debugging purpose
2103 !accum_mat2 is used to store the phonon-linewidhts before interpolation
2104 
2105 !MJV 20070525 NOTE:
2106 !in some of the nmsq routines, in particular this one, the work done to
2107 !calculate accum_mat,accum_mat2 is completely superfluous and will be re-done
2108 !on the interpolated values.
2109 !MG uses them for the QPT output, however, so keep it for consistency for the
2110 !moment.
2111 
2112  do isppol=1,elph_ds%nsppol
2113    do ik_this_proc =1, elph_ds%k_phon%my_nkpt
2114      ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
2115 
2116      ikpt_phonq = FSfullpqtofull(ikpt_phon,iqpt_fullbz)
2117 
2118      gkq_sum_bands(:,:,:) = zero
2119 
2120 !    gkq_sum_bands = \sum_{ib1,ib2} \langle k+q \mid H^{(1)}_{q,\tau_i,\alpha_i} \mid k   \rangle
2121 !    \cdot \langle k   \mid H^{(1)}_{q,\tau_j,\alpha_j} \mid k+q \rangle
2122 !    where ibranch -> \tau_i,\alpha_i  and  jbranch -> \tau_j,\alpha_j
2123 
2124      do ib1=1,elph_ds%nFSband
2125 
2126        sd1 = elph_ds%k_phon%wtk(ib1,ikpt_phon,isppol)      !  weights for distance from the fermi surface
2127 
2128        do ib2=1,elph_ds%nFSband
2129 
2130          sd2 = elph_ds%k_phon%wtk(ib2,ikpt_phonq,isppol)  !  weights for distance from the fermi surface
2131          ibeff = ib2+(ib1-1)*elph_ds%nFSband
2132 
2133          gkq_sum_bands = gkq_sum_bands + &
2134 &         sd1*sd2*reshape(h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol),(/2,elph_ds%nbranch,elph_ds%nbranch/))
2135 
2136        end do !ib2
2137      end do !ib1
2138 !    END loops over bands
2139 
2140 
2141      accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_sum_bands(:,:,:)
2142    end do
2143 !  END loop over kpt_phon
2144 
2145 !  MG20060603
2146 !  do scalar product with the displ_red to calculate the ph lwdth before interpolation (stored in accum_mat2)
2147 
2148    zgemm_tmp_mat = accum_mat(:,:,:,isppol)
2149 
2150    call gam_mult_displ(elph_ds%nbranch, displ_red, zgemm_tmp_mat, tmp_mat2)
2151 
2152    do ipert1=1,elph_ds%nbranch
2153      accum_mat2(1,ipert1,ipert1,isppol) = accum_mat2(1,ipert1,ipert1,isppol) + tmp_mat2(1,ipert1,ipert1)
2154    end do
2155 
2156 !  ENDMG
2157 
2158  end do ! isppol
2159 
2160 end subroutine nmsq_pure_gkk

ABINIT/nmsq_pure_gkk_sumfs [ Functions ]

[ Top ] [ 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

PARENTS

      normsq_gkq

CHILDREN

      gam_mult_displ

SOURCE

2193 subroutine nmsq_pure_gkk_sumfs(accum_mat,accum_mat2,displ_red,elph_ds,FSfullpqtofull,h1_mat_el_sq,iqptirred)
2194 
2195 
2196 !This section has been created automatically by the script Abilint (TD).
2197 !Do not modify the following lines by hand.
2198 #undef ABI_FUNC
2199 #define ABI_FUNC 'nmsq_pure_gkk_sumfs'
2200 !End of the abilint section
2201 
2202  implicit none
2203 
2204 !Arguments ------------------------------------
2205 !scalars
2206  integer,intent(in) :: iqptirred
2207  type(elph_type),intent(in) :: elph_ds
2208 !arrays
2209  integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)
2210  real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
2211  real(dp),intent(inout) :: &
2212 & 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)
2213  real(dp),intent(inout) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
2214  real(dp),intent(inout) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
2215 
2216 !Local variables-------------------------------
2217 !scalars
2218  integer :: ikpt_phon,ikpt_phonq,ib1,ib2,ibeff,ipert1,isppol,iqpt_fullbz
2219  integer :: nbranch,nsppol,nFSband,nkpt_phon
2220  integer :: ik_this_proc
2221  real(dp) :: sd1,sd2
2222  !character(len=500) :: message
2223 !arrays
2224  real(dp) :: gkq_sum_bands(2,elph_ds%nbranch,elph_ds%nbranch)
2225  real(dp) :: tmp_mat2(2,elph_ds%nbranch,elph_ds%nbranch)
2226  real(dp) :: zgemm_tmp_mat(2,elph_ds%nbranch,elph_ds%nbranch)
2227 
2228 ! *************************************************************************
2229 
2230  if (elph_ds%ep_keepbands /= 0) then
2231    MSG_BUG('ep_keepbands should be 0 to average over bands!')
2232  end if
2233 
2234  nbranch   = elph_ds%nbranch
2235  nsppol    = elph_ds%nsppol
2236  nFSband   = elph_ds%nFSband
2237  nkpt_phon = elph_ds%k_phon%nkpt
2238 
2239  iqpt_fullbz = elph_ds%qirredtofull(iqptirred)
2240 
2241 !MG20060603 NOTE:
2242 !accum_mat and accum_mat2 are real, the imaginary part is used for debugging purpose
2243 !accum_mat2 is used to store the phonon-linewidhts before interpolation
2244 
2245  do isppol=1,nsppol
2246    do ik_this_proc =1, elph_ds%k_phon%my_nkpt
2247      ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
2248 
2249 !
2250 !    The index of k+q in the BZ.
2251      ikpt_phonq = FSfullpqtofull(ikpt_phon,iqpt_fullbz)
2252 !
2253 !    gkq_sum_bands =
2254 !    \sum_{ib1,ib2} <k+q| H^{(1)}_{q,\tau_i,\alpha_i} |k> \cdot <k| H^{(1)}_{q,\tau_j,\alpha_j}|k+q>
2255 !
2256 !    where ibranch = (\tau_i,\alpha_i) and  jbranch = (\tau_j,\alpha_j).
2257      gkq_sum_bands(:,:,:) = zero
2258 
2259      do ib1=1,nFSband
2260        sd1 = elph_ds%k_phon%wtk(ib1,ikpt_phon,isppol)      !  weights for distance from the fermi surface
2261 
2262        do ib2=1,nFSband
2263          sd2 = elph_ds%k_phon%wtk(ib2,ikpt_phonq,isppol)  !  weights for distance from the fermi surface
2264          ibeff=ib2+(ib1-1)*nFSband
2265 
2266          gkq_sum_bands = gkq_sum_bands + &
2267 &         sd1*sd2* reshape(h1_mat_el_sq(:,ibeff,:,ik_this_proc,isppol),(/2,nbranch,nbranch/))
2268        end do !ib2
2269      end do !ib1
2270 !
2271 !    gamma matrix contribution in reduced coordinates (ie interpolatable form)
2272 !    The sum over Fermi surface bands is done here, and fed into (ib1,ib2)=(1,1)
2273      h1_mat_el_sq(:,1,:,ik_this_proc,isppol) = reshape(gkq_sum_bands,(/2,nbranch**2/))
2274 
2275      accum_mat(:,:,:,isppol) = accum_mat(:,:,:,isppol) + gkq_sum_bands(:,:,:)
2276    end do ! kpt_phon
2277  end do ! isppol
2278 !
2279 !MG20060603
2280 !do scalar product wit displ_red to calculate the ph lwdth before interpolation (stored in accum_mat2)
2281  do isppol=1,nsppol
2282    zgemm_tmp_mat = accum_mat(:,:,:,isppol)
2283 !
2284    call gam_mult_displ(nbranch, displ_red, zgemm_tmp_mat, tmp_mat2)
2285 
2286    do ipert1=1,nbranch
2287      accum_mat2(1,ipert1,ipert1,isppol) = accum_mat2(1,ipert1,ipert1,isppol) + tmp_mat2(1,ipert1,ipert1)
2288    end do
2289 !
2290  end do
2291 
2292 end subroutine nmsq_pure_gkk_sumfs

ABINIT/normsq_gkq [ Functions ]

[ Top ] [ 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}$.

PARENTS

      read_gkk

CHILDREN

      gam_mult_displ,nmsq_gam,nmsq_gam_sumfs,nmsq_pure_gkk
      nmsq_pure_gkk_sumfs,wrtout,xmpi_sum,zhpev

SOURCE

1507 subroutine normsq_gkq(displ_red,eigvec,elph_ds,FSfullpqtofull,&
1508 &    h1_mat_el_sq,iqptirred,phfrq_tmp,qpt_irred,qdata)
1509 
1510 
1511 !This section has been created automatically by the script Abilint (TD).
1512 !Do not modify the following lines by hand.
1513 #undef ABI_FUNC
1514 #define ABI_FUNC 'normsq_gkq'
1515 !End of the abilint section
1516 
1517  implicit none
1518 
1519 !Arguments ------------------------------------
1520 !scalars
1521  integer,intent(in) :: iqptirred
1522  type(elph_type),intent(inout) :: elph_ds
1523 !arrays
1524  integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)
1525  real(dp),intent(in) :: displ_red(2,elph_ds%nbranch,elph_ds%nbranch)
1526  real(dp),intent(in) :: eigvec(2,elph_ds%nbranch,elph_ds%nbranch)
1527  real(dp),intent(inout) :: &
1528 & 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)
1529  real(dp),intent(in) :: phfrq_tmp(elph_ds%nbranch),qpt_irred(3,elph_ds%nqptirred)
1530  real(dp),intent(out) :: qdata(elph_ds%nbranch,elph_ds%nsppol,3)
1531 
1532 !Local variables-------------------------------
1533 !scalars
1534  integer :: i1,i2,ier,ii,isppol,jbranch,comm
1535  real(dp) :: lambda_tot
1536  character(len=500) :: message
1537 !arrays
1538  real(dp) :: accum_mat(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1539  real(dp) :: accum_mat2(2,elph_ds%nbranch,elph_ds%nbranch,elph_ds%nsppol)
1540  real(dp) :: gam_now2(2,elph_ds%nbranch,elph_ds%nbranch)
1541  real(dp) :: lambda(elph_ds%nsppol)
1542  real(dp),allocatable :: matrx(:,:),val(:),vec(:,:,:)
1543  real(dp),allocatable :: zhpev1(:,:),zhpev2(:)
1544 
1545 ! *************************************************************************
1546 
1547  DBG_ENTER("COLL")
1548 
1549  accum_mat  = zero
1550  accum_mat2 = zero
1551  comm = xmpi_world
1552 
1553  if (elph_ds%ep_scalprod == 1) then
1554 !
1555    if (elph_ds%ep_keepbands == 0) then
1556      call wrtout(std_out,' normsq_gkq : calling nmsq_gam_sumFS',"COLL")
1557      call nmsq_gam_sumFS (accum_mat,accum_mat2,displ_red,eigvec,elph_ds,FSfullpqtofull,&
1558 &     h1_mat_el_sq,iqptirred)
1559 
1560    else if (elph_ds%ep_keepbands == 1) then
1561      call wrtout(std_out,' normsq_gkq : calling nmsq_gam',"COLL")
1562      call nmsq_gam (accum_mat,accum_mat2,displ_red,eigvec,elph_ds,FSfullpqtofull,&
1563 &     h1_mat_el_sq,iqptirred)
1564 
1565    else
1566      write (message,'(a,i0)')' Wrong value for elph_ds%ep_keepbands = ',elph_ds%ep_keepbands
1567      MSG_BUG(message)
1568    end if
1569 !
1570  else if (elph_ds%ep_scalprod == 0) then  ! Interpolate on the pure "matrix of matrix elements" and do the scalar products later.
1571 !
1572    if (elph_ds%ep_keepbands == 0) then
1573      call wrtout(std_out,' normsq_gkq : calling nmsq_pure_gkk_sumFS',"COLL")
1574      call nmsq_pure_gkk_sumFS (accum_mat,accum_mat2,displ_red,elph_ds,FSfullpqtofull,&
1575 &     h1_mat_el_sq,iqptirred)
1576 
1577    else if (elph_ds%ep_keepbands == 1) then
1578      call wrtout(std_out,' normsq_gkq : calling nmsq_pure_gkk',"COLL")
1579 
1580      call nmsq_pure_gkk (accum_mat,accum_mat2,displ_red,elph_ds,FSfullpqtofull,&
1581 &     h1_mat_el_sq,iqptirred)
1582    else
1583      write (message,'(a,i0)')' Wrong value for elph_ds%ep_keepbands = ',elph_ds%ep_keepbands
1584      MSG_BUG(message)
1585    end if
1586 !
1587  else
1588    write (message,'(a,i0)')' Wrong value for elph_ds%ep_scalprod = ',elph_ds%ep_scalprod
1589    MSG_BUG(message)
1590  end if
1591 !end if flag for doing scalar product now.
1592 
1593 
1594 !MG: values without the good prefactor
1595  accum_mat = accum_mat * elph_ds%occ_factor/elph_ds%k_phon%nkpt
1596 
1597 !MG: accum_mat2 contains the line-widhts before the Fourier interpolation
1598  accum_mat2 = accum_mat2 * elph_ds%occ_factor/elph_ds%k_phon%nkpt
1599 
1600 !mpi sum over procs for accum_mat2
1601  call xmpi_sum (accum_mat, comm, ier)
1602  call xmpi_sum (accum_mat2, comm, ier)
1603 
1604 !MG20060531i
1605 !write e-ph quantities before Fourier interpolation
1606 !save e-ph values in the temporary array qdata that will be copied into elph_ds%qgrid_data
1607 
1608  write (message,'(4a,3es16.6,63a)')ch10,                  &
1609 & ' Phonon linewidths before interpolation ',ch10,        &
1610 & ' Q point = ',qpt_irred(:,iqptirred),ch10,('=',ii=1,60),ch10,&
1611 & ' Mode          Frequency (Ha)  Linewidth (Ha)  Lambda '
1612  call wrtout(std_out,message,'COLL')
1613 
1614  lambda_tot = zero
1615  do isppol=1,elph_ds%nsppol
1616    do ii=1,elph_ds%nbranch
1617      lambda(isppol)=zero
1618 !    MG: the tolerance factor is somehow arbitrary
1619      if (abs(phfrq_tmp(ii)) > tol10) lambda(isppol)=accum_mat2(1,ii,ii,isppol)/&
1620 &     (pi*elph_ds%n0(isppol)*phfrq_tmp(ii)**2)
1621      lambda_tot=lambda_tot+lambda(isppol)
1622      write(message,'(i8,es20.6,2es16.6)' )ii,phfrq_tmp(ii),accum_mat2(1,ii,ii,isppol),lambda(isppol)
1623      call wrtout(std_out,message,'COLL')
1624 !    save values
1625      qdata(ii,isppol,1)=phfrq_tmp(ii)
1626      qdata(ii,isppol,2)=accum_mat2(1,ii,ii,isppol)
1627      qdata(ii,isppol,3)=lambda(isppol)
1628    end do !loop over branch
1629  end do !loop over sppol
1630 
1631 !normalize for number of spins
1632  lambda_tot = lambda_tot / elph_ds%nsppol
1633 
1634  write(message,'(61a,44x,es16.6,62a)' )('=',ii=1,60),ch10,lambda_tot,ch10,('=',ii=1,60),ch10
1635  call wrtout(std_out,message,'COLL')
1636 !ENDMG20060531
1637 
1638 !immediately calculate linewidths:
1639  write(std_out,*) 'summed accum_mat = '
1640  write(std_out,'(3(2E18.6,1x))') accum_mat(:,:,:,1)
1641  write(std_out,*) 'summed accum_mat2 = '
1642  write(std_out,'(3(2E18.6,1x))')  (accum_mat2(:,ii,ii,1),ii=1,elph_ds%nbranch)
1643  write(std_out,*) 'displ_red  = '
1644  write(std_out,'(3(2E18.6,1x))') displ_red
1645 
1646  if (elph_ds%ep_scalprod == 1) then
1647    do isppol=1,elph_ds%nsppol
1648 !    Diagonalize gamma matrix at qpoint (complex matrix). Copied from dfpt_phfrq
1649      ier=0
1650      ii=1
1651      ABI_ALLOCATE(matrx,(2,(elph_ds%nbranch*(elph_ds%nbranch+1))/2))
1652      do i2=1,elph_ds%nbranch
1653        do i1=1,i2
1654          matrx(1,ii)=accum_mat2(1,i1,i2,isppol)
1655          matrx(2,ii)=accum_mat2(2,i1,i2,isppol)
1656          ii=ii+1
1657        end do
1658      end do
1659      ABI_ALLOCATE(zhpev1,(2,2*elph_ds%nbranch-1))
1660      ABI_ALLOCATE(zhpev2,(3*elph_ds%nbranch-2))
1661      ABI_ALLOCATE(val,(elph_ds%nbranch))
1662      ABI_ALLOCATE(vec,(2,elph_ds%nbranch,elph_ds%nbranch))
1663      call ZHPEV ('V','U',elph_ds%nbranch,matrx,val,vec,elph_ds%nbranch,zhpev1,zhpev2,ier)
1664 
1665      write (std_out,*) ' normsq_gkq : accumulated eigenvalues isppol ',isppol, ' = '
1666      write (std_out,'(3E18.6)') val
1667      ABI_DEALLOCATE(matrx)
1668      ABI_DEALLOCATE(zhpev1)
1669      ABI_DEALLOCATE(zhpev2)
1670      ABI_DEALLOCATE(vec)
1671      ABI_DEALLOCATE(val)
1672    end do ! isppol
1673 
1674  else if (elph_ds%ep_scalprod == 0) then
1675 
1676 
1677    do isppol=1,elph_ds%nsppol
1678      call gam_mult_displ(elph_ds%nbranch, displ_red, accum_mat(:,:,:,isppol), gam_now2)
1679 
1680      write (std_out,*) ' normsq_gkq : accumulated eigenvalues isppol ', isppol, ' = '
1681      write (std_out,'(3(E14.6,1x))') (gam_now2(1,jbranch,jbranch), jbranch=1,elph_ds%nbranch)
1682      write (std_out,*) ' normsq_gkq : imag part = '
1683      write (std_out,'(3(E14.6,1x))') (gam_now2(2,jbranch,jbranch), jbranch=1,elph_ds%nbranch)
1684    end do ! isppol
1685 
1686  end if
1687 
1688  DBG_EXIT("COLL")
1689 
1690 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

PARENTS

      read_gkk

CHILDREN

      d2sym3

SOURCE

1349 subroutine completeperts(Cryst,nbranch,nFSband,nkpt,nsppol,gkk_flag,h1_mat_el,h1_mat_el_sq,&
1350 &   qpt,symq,qtimrev)
1351 
1352 
1353 !This section has been created automatically by the script Abilint (TD).
1354 !Do not modify the following lines by hand.
1355 #undef ABI_FUNC
1356 #define ABI_FUNC 'completeperts'
1357 !End of the abilint section
1358 
1359  implicit none
1360 
1361 !Arguments ------------------------------------
1362 !scalars
1363  integer,intent(in) :: qtimrev,nbranch,nFSband,nkpt,nsppol
1364  type(crystal_t),intent(in) :: Cryst
1365 !arrays
1366  integer,intent(in) :: symq(4,2,Cryst%nsym)
1367  integer,intent(inout) :: gkk_flag(nbranch,nbranch,nkpt,nsppol)
1368  real(dp),intent(in) :: qpt(3)
1369  real(dp),intent(in) :: h1_mat_el(2,nFSband**2,nbranch,nkpt,nsppol)
1370  real(dp),intent(out) :: h1_mat_el_sq(2,nFSband**2,nbranch**2,nkpt,nsppol)
1371 
1372 !Local variables-------------------------------
1373 !scalars
1374  integer :: ikpt_phon,iatom1,iatom2,ibb,idir1,idir2,ipert1,ipert2,isppol,mpert,natom
1375  real(dp) :: im1,im2,re1,re2,res
1376  character(len=500) :: msg
1377 !arrays
1378  integer,allocatable :: tmpflg(:,:,:,:)
1379  real(dp),allocatable :: tmpval(:,:,:,:,:)
1380 
1381 ! *************************************************************************
1382 
1383 !WARNING! Stupid patch in d2sym3 imposes these matrices to have size natom+2
1384  natom = Cryst%natom
1385  mpert = natom+2
1386 
1387  ABI_ALLOCATE(tmpflg,(3,mpert,3,mpert))
1388  ABI_ALLOCATE(tmpval,(2,3,mpert,3,mpert))
1389 
1390  h1_mat_el_sq = zero
1391 
1392  write (std_out,*) ' completeperts: shape(h1_mat_el_sq) = ', shape(h1_mat_el_sq)
1393 
1394  do isppol=1,nsppol
1395    write(std_out,*)'completeperts: isppol = ', isppol
1396 !
1397    do ikpt_phon=1,nkpt
1398      do ibb=1,nFSband**2
1399 !
1400        tmpval = zero
1401        tmpflg = 0
1402        ! for a fixed k (q) band and sppol construct the gamma matrix for (3 natom)^2 perturbation pairs
1403        do iatom1=1,natom
1404          do idir1=1,3
1405            ipert1 = (iatom1-1)*3+idir1
1406            if (gkk_flag(ipert1,ipert1,ikpt_phon,isppol) < 0) cycle
1407            re1 = h1_mat_el(1,ibb,ipert1,ikpt_phon,isppol)
1408            im1 = h1_mat_el(2,ibb,ipert1,ikpt_phon,isppol)
1409 
1410            do iatom2=1,natom
1411              do idir2=1,3
1412                ipert2 = (iatom2-1)*3+idir2
1413                if (gkk_flag(ipert2,ipert2,ikpt_phon,isppol) < 0) cycle
1414                tmpflg(idir1,iatom1,idir2,iatom2) = 1
1415                re2 = h1_mat_el(1,ibb,ipert2,ikpt_phon,isppol)
1416                im2 = h1_mat_el(2,ibb,ipert2,ikpt_phon,isppol)
1417 !
1418 !              conjg(h1_mat_el_2) * h1_mat_el_1
1419                res =  re1*re2 + im1*im2
1420                tmpval(1,idir1,iatom1,idir2,iatom2) =  res
1421                res =  re1*im2 - im1*re2
1422                tmpval(2,idir1,iatom1,idir2,iatom2) = res
1423 
1424              end do !idir2
1425            end do !iatom2
1426          end do !idir1
1427        end do !iatom1
1428 
1429        ! matrix is symmetrized like a dynamical matrix. No change of band or k
1430        !  in here. This should be checked (if we have to restrict further the symmetry operations)
1431        call d2sym3(tmpflg,tmpval,Cryst%indsym,mpert,natom,Cryst%nsym,qpt,symq,Cryst%symrec,Cryst%symrel,qtimrev,1)
1432        if (sum(tmpflg(:,1:natom,:,1:natom)) /= 3*natom*3*natom) then
1433          write(msg,'(3a,4i0)')&
1434 &         'A perturbation is missing after completion with d2sym3',ch10,&
1435 &         'tmpflg, ikpt_phon, isppol: ',tmpflg,ikpt_phon,isppol
1436          MSG_ERROR(msg)
1437        end if
1438 !
1439 !      Save values for calculation of |gkk|^2
1440        do iatom1=1,natom
1441          do idir1=1,3
1442            ipert1 = (iatom1-1)*3+idir1
1443            do iatom2=1,natom
1444              do idir2=1,3
1445 !
1446 !              mjv 29/10/2007 ipert2 now contains the composite index ip1*nperts+ip2
1447                ipert2 = (iatom2-1)*3 + idir2 + (ipert1-1)*3*natom
1448                h1_mat_el_sq(1,ibb,ipert2,ikpt_phon,isppol) = pi*tmpval(1,idir2,iatom2,idir1,iatom1)
1449                h1_mat_el_sq(2,ibb,ipert2,ikpt_phon,isppol) = pi*tmpval(2,idir2,iatom2,idir1,iatom1)
1450              end do
1451            end do
1452          end do
1453        end do
1454 !
1455      end do !end ibb band dos
1456 !
1457 !    Set flags.
1458      do ipert1=1,3*natom
1459        do ipert2=1,3*natom
1460          if (gkk_flag(ipert2,ipert1,ikpt_phon,isppol) < 0) gkk_flag(ipert2,ipert1,ikpt_phon,isppol) = 1
1461        end do
1462      end do
1463 
1464    end do !end kpt_phon do
1465  end do !end sppol do
1466 
1467  ABI_DEALLOCATE(tmpflg)
1468  ABI_DEALLOCATE(tmpval)
1469 
1470 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

PARENTS

      read_el_veloc

CHILDREN

      hdr_fort_read,hdr_free,wrtout

SOURCE

1222 subroutine inpgkk(eigen1,filegkk,hdr1)
1223 
1224 
1225 !This section has been created automatically by the script Abilint (TD).
1226 !Do not modify the following lines by hand.
1227 #undef ABI_FUNC
1228 #define ABI_FUNC 'inpgkk'
1229 !End of the abilint section
1230 
1231  implicit none
1232 
1233 !Arguments ------------------------------------
1234 !scalars
1235  character(len=fnlen),intent(in) :: filegkk
1236  type(hdr_type), intent(out) :: hdr1
1237 !arrays
1238  real(dp),allocatable,intent(out) :: eigen1(:)
1239 
1240 !Local variables-------------------------------
1241 !scalars
1242  integer :: bantot1
1243  integer :: isppol, ikpt, mband, ikb
1244  integer :: unitgkk, fform, ierr, n1wf, i1wf
1245  type(hdr_type) :: hdr0
1246  real(dp), allocatable :: eigen(:)
1247  character(len=500) :: message
1248 
1249 ! *************************************************************************
1250 
1251  if (open_file(filegkk,message,newunit=unitgkk,form='unformatted',status='old') /= 0) then
1252    MSG_ERROR(message)
1253  end if
1254 
1255 
1256 !read in header of GS file and eigenvalues
1257  call hdr_fort_read(hdr0, unitgkk, fform)
1258  ABI_CHECK(fform /= 0, "hdr_fort_read returned fform == 0")
1259 
1260  mband = maxval(hdr0%nband(:))
1261  ABI_ALLOCATE(eigen,(mband))
1262  call wrtout(std_out,'inpgkk : try to reread GS eigenvalues','COLL')
1263 
1264  do isppol=1,hdr0%nsppol
1265    do ikpt=1,hdr0%nkpt
1266      read (unitgkk,IOSTAT=ierr) eigen(1:hdr0%nband(ikpt))
1267      ABI_CHECK(ierr==0,'reading eigen from gkk file')
1268    end do
1269  end do
1270 
1271  read(unitgkk,IOSTAT=ierr) n1wf
1272  ABI_CHECK(ierr==0,"reading n1wf from gkk file")
1273 
1274  ABI_DEALLOCATE(eigen)
1275  call hdr_free(hdr0)
1276 
1277  if (n1wf > 1) then
1278    write(message,'(3a)')&
1279 &   'several 1wf records were found in the file,',ch10, &
1280 &   'which is not allowed for reading with this routine'
1281    MSG_ERROR(message)
1282  end if
1283 
1284 !read in header of 1WF file
1285  call hdr_fort_read(hdr1, unitgkk, fform)
1286  if (fform == 0) then
1287    write(message,'(a,i0,a)')' 1WF header number ',i1wf,' was mis-read. fform == 0'
1288    MSG_ERROR(message)
1289  end if
1290 
1291  bantot1 = 2*hdr1%nsppol*hdr1%nkpt*mband**2
1292  ABI_ALLOCATE(eigen1, (bantot1))
1293 
1294 
1295 !retrieve 1WF <psi_k+q | H | psi_k> from gkk file and echo to output
1296  ikb = 0
1297  do isppol=1,hdr1%nsppol
1298    do ikpt=1,hdr1%nkpt
1299      read (unitgkk,IOSTAT=ierr) eigen1(ikb+1:ikb+2*hdr1%nband(ikpt)**2)
1300      ikb = ikb + 2*hdr1%nband(ikpt)**2
1301      if (ierr /= 0) then
1302        write(message,'(a,2i0)')'reading eigen1 from gkk file, spin, kpt_idx',isppol,ikpt
1303        MSG_ERROR(message)
1304      end if
1305    end do
1306  end do
1307 
1308  close(unitgkk)
1309 
1310 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

PARENTS

      dfpt_looppert

CHILDREN

      hdr_fort_write,wrtout

SOURCE

738 subroutine outgkk(bantot0,bantot1,outfile,eigen0,eigen1,hdr0,hdr1,mpi_enreg,phasecg)
739 
740 
741 !This section has been created automatically by the script Abilint (TD).
742 !Do not modify the following lines by hand.
743 #undef ABI_FUNC
744 #define ABI_FUNC 'outgkk'
745 !End of the abilint section
746 
747  implicit none
748 
749 !Arguments ------------------------------------
750 !scalars
751  integer,intent(in) :: bantot0,bantot1
752  character(len=fnlen),intent(in) :: outfile
753  type(MPI_type),intent(in) :: mpi_enreg
754  type(hdr_type),intent(inout) :: hdr0,hdr1
755 !arrays
756  real(dp),intent(in) :: eigen0(bantot0),eigen1(2*bantot1)
757  real(dp),intent(in) :: phasecg(2,bantot1)
758 
759 !Local variables-------------------------------
760 !scalars
761  integer :: fform,iband,ikpt,isppol,me,ntot,unitout
762  integer :: iband_off, mband, ierr
763  character(len=500) :: msg
764  real(dp), allocatable :: tmpeig(:)
765 
766 ! *************************************************************************
767 
768 !only master should be writing to disk
769 !Init me
770  me=mpi_enreg%me_kpt
771  if (me /= 0) return
772 
773  call wrtout(std_out,' writing gkk file: '//outfile,"COLL")
774 
775 !initializations
776  fform = 42
777  ntot = 1
778 
779 !open gkk file
780  if (open_file(outfile, msg, newunit=unitout, form='unformatted', status='unknown', action="write") /= 0) then
781    MSG_ERROR(msg)
782  end if
783 
784 !output GS header
785  call hdr_fort_write(hdr0, unitout, fform, ierr)
786  ABI_CHECK(ierr == 0 , "hdr_fort_write returned ierr != 0")
787 
788 !output GS eigenvalues
789  iband=0
790  do isppol=1,hdr0%nsppol
791    do ikpt=1,hdr0%nkpt
792      write (unitout) eigen0(iband+1:iband+hdr0%nband(ikpt))
793      iband=iband+hdr0%nband(ikpt)
794    end do
795  end do
796 
797 !output number of gkk in this file (1)
798  write (unitout) ntot
799 
800 !output RF header
801  call hdr_fort_write(hdr1, unitout, fform, ierr)
802  ABI_CHECK(ierr == 0 , "hdr_fort_write returned ierr != 0")
803 
804 !output RF eigenvalues
805  mband = maxval(hdr1%nband(:))
806  ABI_ALLOCATE(tmpeig,(2*mband**2))
807  iband_off = 0
808  tmpeig(1) = phasecg(1, 1)
809  do isppol = 1, hdr1%nsppol
810    do ikpt = 1, hdr1%nkpt
811      tmpeig = zero
812      do iband = 1, hdr1%nband(ikpt)**2
813        tmpeig (2*(iband-1)+1) = eigen1(2*(iband_off+iband-1)+1)
814        tmpeig (2*(iband-1)+2) = eigen1(2*(iband_off+iband-1)+2)
815      end do
816      write (unitout) tmpeig(1:2*hdr1%nband(ikpt)**2)
817      iband_off = iband_off + hdr1%nband(ikpt)**2
818    end do
819  end do
820  ABI_DEALLOCATE(tmpeig)
821 
822 !close gkk file
823  close (unitout)
824 
825 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

PARENTS

      read_gkk

CHILDREN

SOURCE

 861 subroutine prt_gkk_yambo(displ_cart,displ_red,kpt_phon,h1_mat_el,iqpt,&
 862 &       natom,nFSband,nkpt_phon,phfrq,qptn)
 863 
 864 
 865 !This section has been created automatically by the script Abilint (TD).
 866 !Do not modify the following lines by hand.
 867 #undef ABI_FUNC
 868 #define ABI_FUNC 'prt_gkk_yambo'
 869 !End of the abilint section
 870 
 871  implicit none
 872 
 873 !Arguments ------------------------------------
 874 !scalars
 875  integer,intent(in) :: natom,iqpt
 876  integer,intent(in) :: nFSband,nkpt_phon
 877  !arrays
 878  real(dp),intent(in) :: kpt_phon(3,nkpt_phon)
 879  real(dp),intent(in) :: h1_mat_el(2,nFSband*nFSband,3*natom,nkpt_phon,1)
 880  real(dp),intent(in) :: phfrq(3*natom)
 881  real(dp),intent(in) :: displ_cart(2,3*natom,3*natom)
 882  real(dp),intent(in) :: displ_red(2,3*natom,3*natom)
 883  real(dp),intent(in) :: qptn(3)
 884 
 885 !Local variables-------------------------------
 886  !scalars
 887  integer, save :: firsttime=1
 888  integer :: outunit,ikpt,imode,iband,ibandp,iatom,idir,ibandindex
 889  integer :: jmode, outunit2, outunit3
 890  !arrays
 891  real(dp) :: gkk_mode_dep(2)
 892 ! *************************************************************************
 893 
 894 !if first time round:
 895  if (firsttime==1) then
 896    firsttime=0
 897 !  squash file
 898    outunit=get_unit()
 899    open (unit=outunit,file="yambo_elphon_data",status="REPLACE")
 900    outunit2=get_unit()
 901    open (unit=outunit2,file="yambo_elphon_gkk_bymode",status="replace")
 902    outunit3=get_unit()
 903    open (unit=outunit3,file="yambo_elphon_gkksqtw_bymode",status="replace")
 904 
 905 !  write dimensions
 906    write (outunit,'(a,I6)') 'number of el atoms ', natom
 907    write (outunit2,'(a,I6)') 'number of el atoms ', natom
 908    write (outunit3,'(a,I6)') 'number of el atoms ', natom
 909    write (outunit,'(a,I6)') 'number of ph modes ', 3*natom
 910    write (outunit2,'(a,I6)') 'number of ph modes ', 3*natom
 911    write (outunit3,'(a,I6)') 'number of ph modes ', 3*natom
 912    write (outunit,'(a,I6)') 'number of el bands ', nFSband
 913    write (outunit2,'(a,I6)') 'number of el bands ', nFSband
 914    write (outunit3,'(a,I6)') 'number of el bands ', nFSband
 915 
 916 !  write k-points
 917    write (outunit,'(a,I6)') 'number of k-points ', nkpt_phon
 918    write (outunit2,'(a,I6)') 'number of k-points ', nkpt_phon
 919    write (outunit3,'(a,I6)') 'number of k-points ', nkpt_phon
 920    do ikpt=1,nkpt_phon
 921      write (outunit,'(a,I6,3E20.10)') 'reduced coord kpoint no ', ikpt, kpt_phon(:,ikpt)
 922      write (outunit2,'(a,I6,3E20.10)') 'reduced coord kpoint no ', ikpt, kpt_phon(:,ikpt)
 923      write (outunit3,'(a,I6,3E20.10)') 'reduced coord kpoint no ', ikpt, kpt_phon(:,ikpt)
 924    end do
 925 
 926 !  band energies are not accessible this deep in the code: simpler to get them
 927 !  from elsewhere
 928 
 929    close (outunit)
 930    close (outunit2)
 931    close (outunit3)
 932  end if ! first time round
 933 
 934 !open file
 935  outunit=get_unit()
 936  open (unit=outunit,file="yambo_elphon_data",status="unknown",position="append")
 937 
 938 !qpoint
 939  write (outunit,'(a,I6,3E20.10)') 'reduced coord qpoint no ', iqpt, qptn(:)
 940 
 941 !frequencies
 942  do imode=1,3*natom
 943    write (outunit,'(a,I6,3E20.10)') 'phonon freq no ', imode, phfrq(imode)
 944  end do
 945 
 946 !displacement vector
 947  do imode=1,3*natom
 948    write (outunit,'(a,I6,3E20.10)') 'phonon displ vec no ', imode
 949    do iatom=1,natom
 950      write (outunit,'(3(2E20.10,2x))') displ_cart(:,(iatom-1)*3+1:iatom*3,imode)
 951    end do
 952  end do
 953 
 954 !the beef: matrix elements of the first order hamiltonian for displacement of
 955 !all atoms along all reduced directions
 956  write (outunit,'(a)') ' matrix elements of all perturbations for this q-point'
 957  do ikpt=1,nkpt_phon
 958    write (outunit,'(a,I6)') ' kpoint ', ikpt
 959    imode=0
 960    do iatom=1,natom
 961      do idir=1,3
 962        imode=imode+1
 963        write (outunit,'(a,I6,I6)') ' atom, direction = ', iatom,idir
 964        ibandindex=0
 965        do iband=1,nFSband
 966          do ibandp=1,nFSband
 967            ibandindex=ibandindex+1
 968            write (outunit,'(a,I6,I6,2E20.10)') ' mat el for n,np ', iband,ibandp,&
 969 &           h1_mat_el(:,ibandindex,imode,ikpt,1)
 970          end do !bandp
 971        end do !band
 972      end do !dir
 973    end do !atom
 974  end do
 975 
 976 !blank line
 977  write (outunit,*)
 978  close (outunit)
 979 
 980  outunit2=get_unit()
 981  open (unit=outunit2,file="yambo_elphon_gkk_bymode",status="unknown",position="append")
 982  outunit3=get_unit()
 983  open (unit=outunit3,file="yambo_elphon_gkksqtw_bymode",status="unknown",position="append")
 984 
 985 !qpoint
 986  write (outunit2,'(a,I6,3E20.10)') 'reduced coord qpoint no ', iqpt, qptn(:)
 987  write (outunit3,'(a,I6,3E20.10)') 'reduced coord qpoint no ', iqpt, qptn(:)
 988 
 989 !print out mode-dependent matrix elements
 990  write (outunit2,'(a)') ' matrix elements of all phonon modes for this q-point'
 991  write (outunit3,'(a)') ' 1/w**1/2 times matrix elements of all phonon modes for this q-point'
 992  do ikpt=1,nkpt_phon
 993    write (outunit2,'(a,I6)') ' kpoint ', ikpt
 994    write (outunit3,'(a,I6)') ' kpoint ', ikpt
 995    ibandindex=0
 996    do iband=1,nFSband
 997      do ibandp=1,nFSband
 998        ibandindex=ibandindex+1
 999        write (outunit2,'(a,I6,I6)') ' el bands n,np ', iband,ibandp
1000        write (outunit3,'(a,I6,I6)') ' el bands n,np ', iband,ibandp
1001        do imode=1,3*natom
1002 !        gkk_mode_dep = cg_zdotc(3*natom,displ_red(:,:,imode),h1_mat_el(:,ibandindex,:,ikpt,1))
1003          gkk_mode_dep = zero
1004          do jmode=1,3*natom
1005            gkk_mode_dep(1) = gkk_mode_dep(1) &
1006 &           + displ_red(1,jmode,imode)*h1_mat_el(1,ibandindex,jmode,ikpt,1) &
1007 &           + displ_red(2,jmode,imode)*h1_mat_el(2,ibandindex,jmode,ikpt,1)
1008            gkk_mode_dep(2) = gkk_mode_dep(2) &
1009 &           + displ_red(1,jmode,imode)*h1_mat_el(2,ibandindex,jmode,ikpt,1) &
1010 &           - displ_red(2,jmode,imode)*h1_mat_el(1,ibandindex,jmode,ikpt,1)
1011          end do
1012          write (outunit2,'(a,I6,2E20.10)') ' mat el for phonon mode num = ', imode, gkk_mode_dep
1013          write (outunit3,'(a,I6,2E20.10)') ' 1/w**1/2 * mat el for phonon mode num = ', &
1014 &         imode, gkk_mode_dep/sqrt(two*abs(phfrq(imode))+tol10)
1015        end do !imode
1016      end do !bandp
1017    end do !band
1018  end do
1019 !blank line
1020  write (outunit2,*)
1021  write (outunit3,*)
1022 
1023  close (outunit2)
1024  close (outunit3)
1025 
1026 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-2018 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)

PARENTS

      elphon

CHILDREN

      destroy_kptrank,get_rank_1kpt,hdr_free,inpgkk,mkkptrank

SOURCE

1062 subroutine read_el_veloc(nband_in,nkpt_in,kpt_in,nsppol_in,elph_tr_ds)
1063 
1064 
1065 !This section has been created automatically by the script Abilint (TD).
1066 !Do not modify the following lines by hand.
1067 #undef ABI_FUNC
1068 #define ABI_FUNC 'read_el_veloc'
1069 !End of the abilint section
1070 
1071  implicit none
1072 
1073 !Arguments -----------------------------------
1074 !scalars
1075  integer, intent(in) :: nband_in,nkpt_in,nsppol_in
1076  type(elph_tr_type), intent(inout) :: elph_tr_ds
1077  real(dp), intent(in) :: kpt_in(3,nkpt_in)
1078 
1079 !Local variables-------------------------------
1080 !scalars
1081  integer :: bd2tot_index
1082  integer :: iband,ii,ikpt, ikpt_ddk
1083  integer :: isppol,l1,mband
1084  integer :: bantot1
1085  integer :: unit_ddk
1086  integer :: symrankkpt
1087  character(len=fnlen) :: filnam1,filnam2,filnam3
1088  character(len=500) :: msg
1089  type(hdr_type) :: hdr1
1090  type(kptrank_type) :: kptrank_t
1091 
1092 !arrays
1093  real(dp) :: im_el_veloc(3)
1094  real(dp),allocatable :: eig1_k(:,:)
1095  real(dp),allocatable :: eigen11(:),eigen12(:),eigen13(:)
1096 
1097 ! *********************************************************************************
1098 
1099 !Read data file name
1100 !TODO: this should be standardized and read in anaddb always, not
1101 !conditionally. Otherwise when new files are added to the anaddb files
1102 !file...  Catastrophe!
1103 
1104  write(std_out,*)'enter read_el_veloc '
1105 
1106 !Read data file
1107  if (open_file(elph_tr_ds%ddkfilename,msg,newunit=unit_ddk,form='formatted') /= 0) then
1108    MSG_ERROR(msg)
1109  end if
1110 
1111  rewind(unit_ddk)
1112  read(unit_ddk,'(a)')filnam1       ! first ddk file
1113  read(unit_ddk,'(a)')filnam2       ! second ddk file
1114  read(unit_ddk,'(a)')filnam3       ! third ddk file
1115  close (unit_ddk)
1116 
1117  bantot1 = 2*nband_in**2*nkpt_in*nsppol_in
1118 
1119  call inpgkk(eigen11,filnam1,hdr1)
1120  call hdr_free(hdr1)
1121 
1122  call inpgkk(eigen12,filnam2,hdr1)
1123  call hdr_free(hdr1)
1124 
1125 !we use the hdr1 from the last call - should add some consistency
1126 !testing here, we are trusting users not to mix different ddk files...
1127  call inpgkk(eigen13,filnam3,hdr1)
1128 
1129 !Extract info from the header
1130  if(hdr1%nsppol /= nsppol_in) then
1131    MSG_ERROR('nsspol /= input nsppol')
1132  end if
1133 
1134 !Get mband, as the maximum value of nband(nkpt)
1135  mband=maxval(hdr1%nband(1:hdr1%nkpt))
1136  if (mband /= nband_in) then
1137    MSG_ERROR('nband_in input to read_el_veloc is inconsistent with mband')
1138  end if
1139 
1140  write(std_out,*)
1141  write(std_out,*)                     'readings from read_el_veloc header'
1142  write(std_out,'(a,i8)')              ' natom                =',hdr1%natom
1143  write(std_out,'(a,3i8)')             ' nkpt,nband_in,mband  =',hdr1%nkpt,nband_in,mband
1144  write(std_out,'(a, f10.5,a)' )      ' ecut                 =',hdr1%ecut,' Ha'
1145  write(std_out,'(a,e15.5,a,e15.5,a)' )' fermie               =',hdr1%fermie,' Ha ',hdr1%fermie*Ha_eV,' eV'
1146 
1147  ABI_ALLOCATE(eig1_k,(2*nband_in**2,3))
1148  bd2tot_index = 0
1149  elph_tr_ds%el_veloc=zero
1150 
1151 !need correspondence between the DDK kpoints and the kpt_phon
1152  call mkkptrank (hdr1%kptns,hdr1%nkpt,kptrank_t)
1153 
1154  do isppol=1,nsppol_in
1155    im_el_veloc(:)=zero
1156    do ikpt=1,nkpt_in
1157      call get_rank_1kpt (kpt_in(:,ikpt),symrankkpt, kptrank_t)
1158      ikpt_ddk = kptrank_t%invrank(symrankkpt)
1159      if (ikpt_ddk == -1) then
1160        write(std_out,*)'read_el_veloc ******** error in correspondence between ddk and gkk kpoint sets'
1161        write(std_out,*)' kpt sets in gkk and ddk files must agree.'
1162        MSG_ERROR("Aborting now")
1163      end if
1164      bd2tot_index=2*nband_in**2*(ikpt_ddk-1)
1165 
1166 !    first derivative eigenvalues for k-point
1167      eig1_k(:,1)=eigen11(1+bd2tot_index:2*nband_in**2+bd2tot_index)
1168      eig1_k(:,2)=eigen12(1+bd2tot_index:2*nband_in**2+bd2tot_index)
1169      eig1_k(:,3)=eigen13(1+bd2tot_index:2*nband_in**2+bd2tot_index)
1170 
1171 !    turn el_veloc to cartesian coordinates
1172      do iband=1,nband_in
1173        do l1=1,3
1174          do ii=1,3
1175            elph_tr_ds%el_veloc(ikpt,iband,l1,isppol)=elph_tr_ds%el_veloc(ikpt,iband,l1,isppol)+&
1176 &           hdr1%rprimd(l1,ii)*eig1_k(2*iband-1+(iband-1)*2*nband_in,ii)/two_pi
1177            im_el_veloc(l1)=im_el_veloc(l1)+&
1178 &           hdr1%rprimd(l1,ii)*eig1_k(2*iband+(iband-1)*2*nband_in,ii)/two_pi
1179          end do
1180        end do ! l1
1181      end do
1182    end do
1183  end do ! end isppol
1184 
1185  call destroy_kptrank (kptrank_t)
1186  ABI_DEALLOCATE(eig1_k)
1187  ABI_DEALLOCATE(eigen11)
1188  ABI_DEALLOCATE(eigen12)
1189  ABI_DEALLOCATE(eigen13)
1190 
1191  call hdr_free(hdr1)
1192 
1193  write(std_out,*)'out of read_el_veloc '
1194 
1195 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

PARENTS

      get_all_gkq

CHILDREN

      completeperts,get_rank_1kpt,hdr_bcast,hdr_fort_read,hdr_free,ifc_fourq
      littlegroup_pert,littlegroup_q,mati3inv,normsq_gkq,phdispl_cart2red
      prt_gkk_yambo,wrap2_pmhalf,wrtout,xmpi_bcast

SOURCE

 99 subroutine read_gkk(elph_ds,Cryst,ifc,Bst,FSfullpqtofull,gkk_flag,n1wf,nband,ep_prt_yambo,unitgkk)
100 
101 
102 !This section has been created automatically by the script Abilint (TD).
103 !Do not modify the following lines by hand.
104 #undef ABI_FUNC
105 #define ABI_FUNC 'read_gkk'
106 !End of the abilint section
107 
108  implicit none
109 
110 !Arguments ------------------------------------
111 !scalars
112  integer,intent(in) :: n1wf,nband,unitgkk,ep_prt_yambo
113  type(crystal_t),intent(in) :: Cryst
114  type(ifc_type),intent(in) :: ifc
115  type(ebands_t),intent(in) :: Bst
116  type(elph_type),intent(inout) :: elph_ds
117 !arrays
118  integer,intent(in) :: FSfullpqtofull(elph_ds%k_phon%nkpt,elph_ds%nqpt_full)
119  integer,intent(out) :: gkk_flag(elph_ds%nbranch,elph_ds%nbranch,elph_ds%k_phon%my_nkpt,elph_ds%nsppol,elph_ds%nqpt_full)
120 
121 !Local variables-------------------------------
122 !scalars
123  integer :: nsppol,nbranch,nFSband,minFSband,comm,use_sym
124  integer :: fform,i1wf,ikpt_phon,iatom1,iatom2
125  integer :: ib,ib1,ib2,ibb,ibranch,idir,idir1,idir2,ierr,ii,ikpt1
126  integer :: ipert,ipert1,ipert2,iqptirred,iqptfull,isppol,isym1
127  integer :: itim1,jkpt_phon,new
128  integer :: nsym1,qtimrev,syuse
129  integer :: tdonecompl,test_flag,verify
130  integer :: nqptirred_local
131  integer :: master, me
132  integer :: symrankkpt, ikpt1_phon, ik_this_proc
133  real(dp) :: res,ss,timsign
134  character(len=500) :: msg
135  type(hdr_type) :: hdr1
136 !arrays
137  integer :: FSirrtok(3,elph_ds%k_phon%nkpt)
138  integer :: symaf1(Cryst%nsym),symq(4,2,Cryst%nsym)
139  integer :: symrc1(3,3,Cryst%nsym),symrl1(3,3,Cryst%nsym)
140  integer :: tmpflg(3,Cryst%natom+2,3,Cryst%natom+2)
141  real(dp) :: displ_cart(2,3*Cryst%natom,3*Cryst%natom)
142  real(dp) :: displ_red(2,3*Cryst%natom,3*Cryst%natom)
143  real(dp) :: eigvec(2,3*Cryst%natom,3*Cryst%natom),kpt(3),phfrq_tmp(3*Cryst%natom),redkpt(3)
144  real(dp) :: qptirred_local(3,n1wf)
145  real(dp) :: tnons1(3,Cryst%nsym)
146  real(dp),allocatable :: eigen1(:,:,:),gkk_qpt_tmp(:,:,:,:)
147  real(dp),allocatable :: h1_mat_el(:,:,:,:,:),h1_mat_el_sq(:,:,:,:,:)
148  real(dp),allocatable :: qdata(:,:,:),qdata_tmp(:,:,:,:)
149 
150 ! *************************************************************************
151 
152  ABI_UNUSED(Bst%bantot)
153 
154  use_sym   = 1
155  nsppol    = elph_ds%nsppol
156  nbranch   = elph_ds%nbranch
157  if (ep_prt_yambo==1) then
158    nFSband = nband
159    minFSband = 1
160  else
161    nFSband   = elph_ds%nFSband
162    minFSband = elph_ds%minFSband
163  end if
164 
165 !init values for parallelization
166  comm = xmpi_world
167  me = xmpi_comm_rank(comm)
168  master = 0
169 
170  ABI_STAT_ALLOCATE(h1_mat_el,(2, nFSband**2, nbranch, elph_ds%k_phon%my_nkpt, nsppol), ierr)
171  ABI_CHECK(ierr==0, 'trying to allocate array h1_mat_el')
172  h1_mat_el= zero
173 
174  ABI_STAT_ALLOCATE(h1_mat_el_sq,(2, nFSband**2, nbranch**2,elph_ds%k_phon%my_nkpt, nsppol), ierr)
175  ABI_CHECK(ierr==0, 'trying to allocate array h1_mat_el_sq')
176  h1_mat_el_sq = zero
177 
178  ABI_ALLOCATE(elph_ds%qirredtofull,(elph_ds%nqptirred))
179 
180 !MG array to store the e-ph quantities calculated over the input Q-grid
181  ABI_ALLOCATE(qdata_tmp,(elph_ds%nqptirred,nbranch,nsppol,3))
182  qdata_tmp=zero
183 
184  nqptirred_local=0 !zero number of irred q-points found
185  qptirred_local(:,:)=zero
186 
187  gkk_flag = -1
188 
189  if (elph_ds%gkqwrite ==0) then
190    elph_ds%gkk_qpt = zero
191 
192  else if (elph_ds%gkqwrite == 1) then
193    ABI_STAT_ALLOCATE(gkk_qpt_tmp,(2,elph_ds%ngkkband**2,nbranch**2,nsppol), ierr)
194    ABI_CHECK(ierr==0, 'trying to allocate array gkk_qpt_tmp')
195    gkk_qpt_tmp = zero
196    do iqptirred=1,elph_ds%nqptirred*elph_ds%k_phon%nkpt
197      write (elph_ds%unitgkq,REC=iqptirred) gkk_qpt_tmp
198    end do
199    ABI_DEALLOCATE(gkk_qpt_tmp)
200 
201  else
202    write (msg,'(a,i0)')' Wrong values for gkqwrite = ',elph_ds%gkqwrite
203    MSG_BUG(msg)
204  end if !gkqwrite
205 
206 !===========================================================
207 !Loop over all files we have
208 !read in header for perturbation
209 !should check that all files are complete, have same header
210 !(taking into account the symmetries for the qpoint),
211 !represent the correct qpoints ...
212 !MG: this task should be performed in mrggkk
213 !===========================================================
214 
215  ABI_ALLOCATE(eigen1,(2,nband,nband))
216  do i1wf=1,n1wf
217 
218    if (master == me) then
219      write (msg,'(2a,i4,a,i4)')ch10,' read_gkk : reading 1WF header # ',i1wf,' /',n1wf
220      call wrtout(std_out,msg,'COLL')
221 
222 !    Could check for compatibility of natom, kpt grids, ecut, qpt with DDB grid...
223 !    MG: Also this task should be done in mrggkk
224 
225      call hdr_fort_read(hdr1, unitgkk, fform)
226      if (fform == 0) then
227        write (msg,'(a,i0,a)')' 1WF header number ',i1wf,' was mis-read. fform == 0'
228        MSG_ERROR(msg)
229      end if
230 
231      write(msg,'(a,i4)')' read_gkk : have read 1WF header #',i1wf
232      call wrtout(std_out,msg,'COLL')
233      write (msg,'(2a,i4,a)')ch10,' read_gkk : # of kpt for this perturbation: ',hdr1%nkpt,ch10
234      call wrtout(std_out,msg,'COLL')
235 
236    end if
237 
238 !  broadcast data to all nodes:
239    call hdr_bcast(hdr1, master, me, comm)
240 
241 !  Find qpoint in full grid
242    new=1
243    do iqptfull=1,elph_ds%nqpt_full
244      kpt(:) = hdr1%qptn(:) - elph_ds%qpt_full(:,iqptfull)
245      call wrap2_pmhalf(kpt(1),redkpt(1),res)
246      call wrap2_pmhalf(kpt(2),redkpt(2),res)
247      call wrap2_pmhalf(kpt(3),redkpt(3),res)
248      ss=redkpt(1)**2+redkpt(2)**2+redkpt(3)**2
249      if(ss < tol6) then
250        new = 0
251        exit !exit with iqptfull
252      end if
253    end do !iqptfull
254 
255    if (new == 1) then
256 !    Test should be at the end: dont care if there are additional
257 !    qpts in gkk file which are not on the main grid. Ignore them.
258      write (msg,'(4a,3es16.6,2a)')ch10,&
259 &     ' read_gkk : WARNING-  ',ch10,&
260 &     ' qpoint = ',hdr1%qptn(:),ch10,&
261 &     ' not found in the input q-grid. Ignoring this point '
262      call wrtout(ab_out,msg,'COLL')
263      call wrtout(std_out,msg,'COLL')
264      if (me == master) then
265        do isppol=1,hdr1%nsppol
266          do ikpt1=1,hdr1%nkpt
267            read(unitgkk) ((eigen1(:,ii,ib),ii=1,nband),ib=1,nband)
268          end do
269        end do
270      end if
271 
272      cycle !cycle the loop on i1wf
273    end if !end if (new ==1)
274 
275 
276 !  Check whether other pieces of the DDB have used this qpt already
277    new=1
278    do iqptirred=1,nqptirred_local
279      kpt(:) = qptirred_local(:,iqptirred) - hdr1%qptn(:)
280      call wrap2_pmhalf(kpt(1),redkpt(1),res)
281      call wrap2_pmhalf(kpt(2),redkpt(2),res)
282      call wrap2_pmhalf(kpt(3),redkpt(3),res)
283      ss=redkpt(1)**2+redkpt(2)**2+redkpt(3)**2
284      if(ss < tol6) then
285        new=0
286        exit  !MG We can use this information to avoid recalculating the dynamical matrix
287      end if !but we need to use a fixed format in GKK!
288    end do !iqptirred
289 
290    if (new==1) then  !we have a new valid irreducible qpoint, add it!
291      nqptirred_local = nqptirred_local+1
292      if (nqptirred_local > elph_ds%nqptirred) then
293        write (msg, '(a,a,a,i6,i6)') &
294 &       'found too many qpoints in GKK file wrt anaddb input ', ch10, &
295 &       'nqpt_anaddb nqpt_gkk = ', elph_ds%nqptirred, nqptirred_local
296        MSG_ERROR(msg)
297      end if
298      qptirred_local(:,nqptirred_local) = hdr1%qptn(:)
299      iqptirred = nqptirred_local
300      tdonecompl = 0
301      h1_mat_el = zero
302    end if
303 
304 !  now iqptirred is the index of the present qpoint in the array qptirred_local
305 !  and iqptfull is the index in the full qpt_full array for future reference
306    elph_ds%qirredtofull(iqptirred) = iqptfull
307 
308    write (msg,'(a,i5,a,3es16.8)')&
309 &   ' read_gkk : full zone qpt number ',iqptfull,' is ',elph_ds%qpt_full(:,iqptfull)
310    call wrtout(std_out,msg,'COLL')
311 
312 !  if this perturbation has already been filled (overcomplete gkk)
313 !  check only 1st kpoint and spinpol, then check others
314    verify = 0
315    if (gkk_flag(hdr1%pertcase,hdr1%pertcase,1,1,elph_ds%qirredtofull(iqptirred)) /= -1) then
316 !
317      do isppol=1,nsppol
318        do ik_this_proc=1,elph_ds%k_phon%my_nkpt
319          if (gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) == -1) then
320            write (std_out,*)" hdr1%pertcase,ik_this_proc,iqptirred",hdr1%pertcase,ik_this_proc,iqptirred
321            MSG_ERROR('Partially filled perturbation ')
322          end if
323        end do ! ikpt_phon
324      end do ! isppol
325 !
326      MSG_WARNING(' gkk perturbation is already filled')
327      write(std_out,*)' hdr1%pertcase,iqptirred,iqptfull = ',hdr1%pertcase,iqptirred,iqptfull,&
328 &     gkk_flag(hdr1%pertcase,hdr1%pertcase,1,1,elph_ds%qirredtofull(iqptirred))
329      verify = 1
330      write (125,*) '# matrix elements for symmetric perturbation'
331 !    Instead of reading eigen1 into void, verify == 1 checks them later on wrt values in memory
332    end if !gkk_flag
333 
334 !  Examine the symmetries of the q wavevector
335 !  these will be used to complete the perturbations for other atoms and idir
336    if (ep_prt_yambo==1) then
337      ! If one wants to print GKKs along phonon modes, it mean mixing of
338      ! perturbations with differnt jauge. Symmetries must then be disable.
339      call littlegroup_q(Cryst%nsym,qptirred_local(:,iqptirred),symq,Cryst%symrec,Cryst%symafm,qtimrev,prtvol=0,use_sym=0)
340    else
341      call littlegroup_q(Cryst%nsym,qptirred_local(:,iqptirred),symq,Cryst%symrec,Cryst%symafm,qtimrev,prtvol=0)
342    end if
343 
344    ! Determine dynamical matrix, phonon frequencies and displacement vector for qpoint
345    !call wrtout(std_out,' read_gkk: calling inpphon to calculate the dynamical matrix','COLL')
346 
347    call ifc_fourq(ifc,cryst,qptirred_local(:,iqptirred),phfrq_tmp,displ_cart,out_eigvec=eigvec)
348 
349 !  Get displacement vectors for all branches in reduced coordinates
350 !  used in scalar product with H(1)_atom,idir  matrix elements
351 !  Calculate $displ_red = displ_cart \cdot gprimd$ for each phonon branch
352 
353    call phdispl_cart2red(Cryst%natom,Cryst%gprimd,displ_cart,displ_red)
354 
355 !  prefactors for gk+q,n\prime;k,n matrix element
356 !  COMMENT : in decaft there is a weird term in the mass factor, of M-zval(species)
357 !  dont know why. Not needed to reproduce decaft results, though...
358 !  weight is squared in evaluation of
359 !  gamma_{k,q,j} = 2 \pi omega_{q,j} sum_{nu,nu\prime} |g^{q,j}_{k+q,nu\prime; k,nu}|^2
360 !  normally cancels with the 2 \pi omega_{q,j} factor in front of the sum...
361 
362 !  hdr1%pertcase = idir + (ipert-1)*3 where ipert=iatom in the interesting cases
363    idir = mod (hdr1%pertcase-1,3)+1
364    ipert = int(dble(hdr1%pertcase-idir)/three)+1
365 
366    write (msg,'(4a,i3,a,i3,a,i4,a)')ch10,&
367 &   ' read_gkk : calling littlegroup_pert to examine the symmetries of the full perturbation ',ch10,&
368 &   ' idir = ',idir,' ipert = ',ipert,' and Q point = ',iqptirred,ch10
369    call wrtout(std_out,msg,'COLL')
370 
371 !  Examine the symmetries of the full perturbation these will be used to complete the kpoints
372 !  DOESNT USE TIME REVERSAL IN littlegroup_pert except for gamma
373 
374    syuse=0
375 
376    call littlegroup_pert(Cryst%gprimd,idir,Cryst%indsym,ab_out,ipert,Cryst%natom,Cryst%nsym,nsym1,2,Cryst%symafm,symaf1,&
377 &   symq,Cryst%symrec,Cryst%symrel,symrl1,syuse,Cryst%tnons,tnons1)
378 
379    do isym1=1,nsym1
380      call mati3inv(symrl1(:,:,isym1),symrc1(:,:,isym1))
381    end do
382    FSirrtok = 0
383 
384 !  ========================================================
385 !  Loop over irred kpts in file, and fill the default gkk
386 !  ========================================================
387 
388 !  MG NOTE : in the present implementation, if nsppol /=1 the code stops in rchkGSheader!
389    do isppol=1,hdr1%nsppol !Loop over spins is trivial? Not tested.
390      write (std_out,*) ' read_gkk : isppol = ', isppol
391 
392      do ikpt1=1,hdr1%nkpt   !Loop over irred kpoints, WARNING  nkpt depends on qpoint and symmetry!
393 !
394 !      this is the main read of the gkk matrix elements from the file (eigen1 arrays)
395 !      it has to be done exactly nsppol*nkpt times, and the kpt_phon are completed
396 !      where appropriate in the loop below (normally succeeding only once for each kpt)
397 !
398        if (master == me) then
399          read(unitgkk) ((eigen1(:,ii,ib),ii=1,nband),ib=1,nband)
400        end if
401 
402 !      MPI broadcast data to all nodes:
403        call xmpi_bcast(eigen1, master, comm, ierr)
404 
405 !      find place of irred k in k_phon
406 !      the kpoints in the file (kptns) could be ordered arbitrarily
407        call get_rank_1kpt (hdr1%kptns(:,ikpt1)-qptirred_local(:,iqptirred), &
408 &       symrankkpt, elph_ds%k_phon%kptrank_t)
409        ikpt1_phon = elph_ds%k_phon%kptrank_t%invrank(symrankkpt)
410        if (ikpt1_phon < 0) then
411          write (msg,'(a,3es16.6,a)')&
412 &         ' irred k ',hdr1%kptns(:,ikpt1),' was not found in full grid'
413          MSG_ERROR(msg)
414        end if
415 !      find correspondence between this kpt_phon and the others
416 !      symrc1 conserves perturbation as well as qpoint
417 !      add to FSirrtok list
418        do isym1=1,nsym1
419          do itim1=0,qtimrev
420            timsign=one-two*itim1
421            kpt(:) = timsign*matmul(symrc1(:,:,isym1), elph_ds%k_phon%kpt(:,ikpt1_phon))
422 
423            call get_rank_1kpt (kpt,symrankkpt,elph_ds%k_phon%kptrank_t)
424            jkpt_phon = elph_ds%k_phon%kptrank_t%invrank(symrankkpt)
425 
426            if (jkpt_phon > 0) then
427              FSirrtok(1,jkpt_phon) = ikpt1_phon
428              FSirrtok(2,jkpt_phon) = isym1
429              FSirrtok(3,jkpt_phon) = itim1
430            else
431              write (msg,'(a,3es16.6,a,i5,a,i4,a)')&
432 &             ' sym equivalent of kpt ',hdr1%kptns(:,ikpt1),' by sym ',&
433 &             isym1,' and itime ',itim1,' was not found'
434              MSG_ERROR(msg)
435            end if
436          end do !itim1
437        end do !isim1
438 
439 
440        !
441        !  Here check if the symmetry-copied gkk at new k point is equal to the one found in the file for non-irreducible point
442        !  NB This is DEBUG code
443        !
444        if (verify == 1 .and. elph_ds%k_phon%my_kpt(ikpt1_phon) == me) then
445          do ik_this_proc = 1, elph_ds%k_phon%my_nkpt
446            if (elph_ds%k_phon%my_ikpt(ik_this_proc) == ikpt1_phon) exit
447          end do
448          do ib1=1,nFSband
449            do ib2=1,nFSband
450              ibb = (ib1-1)*nFSband+ib2
451              write (125,'(2(2E16.6,2x))') h1_mat_el(:,ibb,hdr1%pertcase,ik_this_proc,isppol),&
452 &             eigen1(:,minFSband-1+ib2,minFSband-1+ib1)
453            end do
454          end do
455        end if !verify end DEBUG code
456 
457 
458        do ik_this_proc = 1, elph_ds%k_phon%my_nkpt
459 !        should I be dealing with this k-point?
460          jkpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
461 
462 !        does present ikpt1 contribute to this k-point?
463          if (FSirrtok(1,jkpt_phon) /= ikpt1_phon) cycle
464 
465 !        if this kpoint has already been filled (overcomplete gkk)
466          if (gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) /= -1) then
467            MSG_WARNING("gkk element is already filled")
468            write(std_out,*)' hdr1%pertcase,ik_this_proc,isppol,iqptirred = ',&
469 &           hdr1%pertcase,ik_this_proc,isppol,iqptirred,&
470 &           gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred))
471 !           exit
472          end if !gkk_flag
473 
474 !        ===============================================================
475 !        TODO: if there is a phase factor in swapping k-points, insert it here in copy to h1_mat_el
476 !        as a function of symops in FSirrtok
477 !        complete gkk for symmetric ikpt_phon with sym1 which conserve
478 !        the full perturbation+qpoint
479 !        Not tested explicitly, but the results for Pb using reduced kpts look good
480 !        should do same RF calculation with nsym=1 and check
481 !        ===============================================================
482 
483 !        save this kpoint
484          do ib1=1,nFSband
485            do ib2=1,nFSband
486              ibb = (ib1-1)*nFSband+ib2
487 
488 !            real
489              res=eigen1(1,minFSband-1+ib2,minFSband-1+ib1)
490              h1_mat_el(1,ibb,hdr1%pertcase,ik_this_proc,isppol) = res
491 
492 !            imag
493              res=eigen1(2,minFSband-1+ib2,minFSband-1+ib1)
494              h1_mat_el(2,ibb,hdr1%pertcase,ik_this_proc,isppol) = res
495            end do !ib2
496          end do !ib1
497 !        if jkpt is equal to ikpt1_phon (if clause above) flag == 3
498          if (FSirrtok(2,jkpt_phon) == 1) then
499            gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) = 3
500 !          if jkpt_phon comes from ikpt1_phon flag == 2 with some symop
501          else
502            gkk_flag(hdr1%pertcase,hdr1%pertcase,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) = 2
503          end if
504 
505        end do !jkpt_phon
506 
507 !      ===============================================================
508 !      we now have contribution to g(k+q,k; \kappa,\alpha) from one
509 !      kpoint,and one perturbation,
510 !      NB: each perturbation will contribute to all the modes later!
511 !
512 !      SHOULD ONLY DO THIS FOR THE SYMS NEEDED
513 !      TO COMPLETE THE PERTURBATIONS!!!
514 !      ================================================================
515 
516      end do !ikpt1
517    end do !isppol
518 
519 ! 14 Jan 2014 removed test on verify - in new scheme full BZ is read in and should be used to avoid phase errors
520 !   if (verify == 1) cycle
521 
522 !  Checks on irred grid provided and on gkk_flag accumulated up to now
523    if (elph_ds%tuniformgrid == 1) then  ! check if irred kpoints found reconstitute the FS kpts
524      do ikpt_phon=1,elph_ds%k_phon%nkpt
525        if (FSirrtok(1,ikpt_phon) == 0) then
526          write(msg,'(a,3es16.6,2a)')&
527 &         ' kpt = ',elph_ds%k_phon%kpt(:,ikpt_phon),ch10,&
528 &         ' is not the symmetric of one of those found in the GKK file'
529          MSG_ERROR(msg)
530        end if
531      end do !ikpt_phon
532 
533 !    normally at this point we have used all the gkk for all kpoints on the FS
534 !    for the given irred perturbation: check
535      do ik_this_proc = 1, elph_ds%k_phon%my_nkpt
536        ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
537 
538        if (gkk_flag(hdr1%pertcase, hdr1%pertcase, ik_this_proc, 1, elph_ds%qirredtofull(iqptirred)) == -1) then
539          write (msg,'(a,i3,a,3es18.6,2a,i3,a,i3,a,3es18.6,a,a,i4,a,a)')&
540 &         ' For irreducible qpt ', iqptirred,' = ',qptirred_local(:,iqptirred),ch10, &
541 &         ' the gkk element : pertcase = ',hdr1%pertcase,' ik_this_proc = ',ik_this_proc, &
542 &         ' kpt = ',elph_ds%k_phon%kpt(:,ikpt_phon),ch10,&
543 &         ' and isppol ',1,ch10,&
544 &         ' was not found by symmetry operations on the irreducible kpoints given'
545          MSG_ERROR(msg)
546        end if
547      end do !ikpt_phon
548    end if ! end elph_ds%tuniformgrid == 1 checks
549 
550    write(msg,'(a,i0)')' read_gkk : Done completing the kpoints for pertcase ',hdr1%pertcase
551    call wrtout(std_out,msg,'COLL')
552 
553    tmpflg(:,:,:,:) = 0
554 
555    do idir1=1,3
556      do iatom1=1,Cryst%natom
557        ipert1 = (iatom1-1)*3+idir1
558        do idir2=1,3
559          do iatom2=1,Cryst%natom
560            ipert2 = (iatom2-1)*3+idir2
561            if (gkk_flag(ipert1,ipert1,1,1,elph_ds%qirredtofull(iqptirred)) >= 0 .and. &
562 &           gkk_flag(ipert2,ipert2,1,1,elph_ds%qirredtofull(iqptirred)) >= 0) then
563              tmpflg(idir1,iatom1,idir2,iatom2) = 1
564            end if
565          end do
566        end do
567      end do
568    end do
569 
570 
571 !  ===============================================
572 !  Full test: need all perturbations explicitly
573 !  ===============================================
574 
575    test_flag = 0
576    if (sum(tmpflg(:,1:Cryst%natom,:,1:Cryst%natom)) == (3*Cryst%natom)**2 .and. tdonecompl == 0) test_flag = 1
577 
578    write(std_out,*)'read_gkk: tdonecompl = ', tdonecompl
579 
580 !  de-activate completion of perts by symmetry for now.
581 !  Must be called when all irreducible perturbations are in memory!!!!
582    if (test_flag == 1 .and. tdonecompl == 0) then
583 
584 !    write(std_out,*) ' read_gkk : enter fxgkkphase before completeperts'
585 !    call fxgkkphase(elph_ds,gkk_flag,h1_mat_el,iqptirred)
586 
587      if (ep_prt_yambo==1) then
588        if (elph_ds%k_phon%my_nkpt /= elph_ds%k_phon%nkpt) then
589          write (msg, '(a)') 'prt_gkk_yambo can not handle parallel anaddb yet'
590          MSG_ERROR(msg)
591        end if
592        call prt_gkk_yambo(displ_cart,displ_red,elph_ds%k_phon%kpt,h1_mat_el,iqptirred,&
593 &       Cryst%natom,nFSband,elph_ds%k_phon%my_nkpt,phfrq_tmp,hdr1%qptn)
594      end if
595 
596 !    ========================================================================
597 !    Now use more general symops to complete the other equivalent
598 !    perturbations: the kpoints are also shuffled by these symops
599 !    afterwards h1_mat_el_sq contains gamma_\tau\alpha,\tau'\alpha' in reduced coordinates
600 !
601 !    \gamma_{\tau'\alpha',\tau\alpha} =
602 !    <psi_{k+q,ib2}| H(1)_{\tau'\alpha'}| psi_{k,ib1}>* \cdot
603 !    <psi_{k+q,ib2}| H(1)_{\tau \alpha }| psi_{k,ib1}>
604 !
605 !    ========================================================================
606 
607      call completeperts(Cryst,nbranch,nFSband,elph_ds%k_phon%my_nkpt,nsppol,&
608 &     gkk_flag(:,:,:,:,elph_ds%qirredtofull(iqptirred)),h1_mat_el,h1_mat_el_sq,qptirred_local(:,iqptirred),symq,qtimrev)
609 
610      tdonecompl = 1
611    end if
612 
613 !  ==============================================================
614 !  if we have all the perturbations for this qpoint, proceed
615 !  with scalar product, norm squared, and add weight factors
616 !
617 !  SHOULD HAVE A TEST SO h1_mat_el IS NOT OVERWRITTEN
618 !  BEFORE PREVIOUS QPOINT IS FINISHED!!!!!
619 !  ==============================================================
620 
621    test_flag = 1
622    do isppol=1,nsppol
623      do ik_this_proc = 1, elph_ds%k_phon%my_nkpt
624        do ibranch=1,nbranch
625          if (gkk_flag (ibranch,ibranch,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) == -1) then
626            test_flag = 0
627            exit
628          end if
629        end do
630      end do
631    end do
632 
633    if (test_flag /= 0) then
634      call wrtout(std_out,' read_gkk : enter normsq_gkq',"COLL")
635 
636 !    MG temporary array to save ph-linewidths before Fourier interpolation
637      ABI_ALLOCATE(qdata,(nbranch,nsppol,3))
638      qdata(:,:,:)=zero
639 
640      call normsq_gkq(displ_red,eigvec,elph_ds,FSfullpqtofull,&
641 &     h1_mat_el_sq,iqptirred,phfrq_tmp,qptirred_local,qdata)
642 
643 !    save gkk_qpt, eventually to disk, for bands up to ngkkband,
644 !    NB: if the sum over bands has been performed ngkkband is 1 instead of nFSband
645      if (elph_ds%gkqwrite == 0) then
646        elph_ds%gkk_qpt(:,:,:,:,:,iqptirred) = h1_mat_el_sq(:,1:elph_ds%ngkkband*elph_ds%ngkkband,:,:,:)
647      else
648 !      write all kpoints to disk
649        write (std_out,*) 'size of record to be written: ', 8  * 2*elph_ds%ngkkband*elph_ds%ngkkband*&
650 &       elph_ds%nbranch*elph_ds%nbranch*elph_ds%k_phon%my_nkpt*elph_ds%nsppol
651        inquire(unit=elph_ds%unitgkq, recl=isppol)
652        write (std_out,*) 'recl =', isppol
653        write (std_out,*) 'iqptirred ', iqptirred
654        do ik_this_proc = 1, elph_ds%k_phon%my_nkpt
655          write (elph_ds%unitgkq,REC=((iqptirred-1)*elph_ds%k_phon%my_nkpt+ik_this_proc)) &
656 &         h1_mat_el_sq(:,1:elph_ds%ngkkband*elph_ds%ngkkband,:,ik_this_proc,:)
657        end do
658      end if
659 
660      qdata_tmp(iqptirred,:,:,:)=qdata(:,:,:)
661      ABI_DEALLOCATE(qdata)
662    end if
663 
664    call hdr_free(hdr1)
665 
666  end do !of i1wf
667 
668 !got all the gkk perturbations
669 
670  ABI_DEALLOCATE(eigen1)
671  ABI_DEALLOCATE(h1_mat_el)
672  ABI_DEALLOCATE(h1_mat_el_sq)
673 
674  if (nqptirred_local /= elph_ds%nqptirred) then
675    write (msg, '(3a,i0,i0)') &
676 &   ' Found wrong number of qpoints in GKK file wrt anaddb input ', ch10, &
677 &   ' nqpt_anaddb nqpt_gkk = ', elph_ds%nqptirred, nqptirred_local
678    MSG_ERROR(msg)
679  end if
680 
681 !normally at this point we have the gkk for all kpoints on the FS
682 !for all the perturbations. Otherwise a 1WF file is missing.
683 !NOTE: still havent checked the qpoint grid completeness
684  do iqptirred=1,elph_ds%nqptirred
685    do isppol=1,nsppol
686      do ik_this_proc = 1, elph_ds%k_phon%my_nkpt
687        ikpt_phon = elph_ds%k_phon%my_ikpt(ik_this_proc)
688        do ipert=1,nbranch
689          if (gkk_flag(ipert,ipert,ik_this_proc,isppol,elph_ds%qirredtofull(iqptirred)) == -1) then
690            write (msg,'(a,i5,1x,i5,1x,i5,1x,i5,a,a)')&
691 &           ' gkk element',ipert,ikpt_phon,isppol,iqptirred,' was not found by symmetry operations ',&
692 &           ' on the irreducible perturbations and qpoints given'
693            MSG_ERROR(msg)
694          end if
695        end do !ipert
696      end do !ik_this_proc
697    end do !isppol
698  end do !iqptirred
699 
700  call wrtout(std_out,'read_gkk : done completing the perturbations (and checked!)','COLL')
701 
702 !MG save phonon frequencies, ph-linewidths and lambda(q,n) values before Fourier interpolation
703  ABI_ALLOCATE(elph_ds%qgrid_data,(elph_ds%nqptirred,nbranch,nsppol,3))
704 
705  do iqptirred=1,elph_ds%nqptirred
706    elph_ds%qgrid_data(iqptirred,:,:,:)=qdata_tmp(iqptirred,:,:,:)
707  end do
708 
709  ABI_DEALLOCATE(qdata_tmp)
710 
711 end subroutine read_gkk