TABLE OF CONTENTS


ABINIT/m_paw_overlap [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_overlap

FUNCTION

  This module contains several routines used to compute the overlap between
  2 wave-functions (PAW only), and associated tools.
  Mainly used in Berry phase formalism.

COPYRIGHT

 Copyright (C) 2018-2018 ABINIT group (JWZ,TRangel,BA,FJ,PHermet)
 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

19 #if defined HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 
23 #include "abi_common.h"
24 
25 MODULE m_paw_overlap
26 
27  use defs_basis
28  use defs_abitypes
29  use m_errors
30  use m_abicore
31  use m_xmpi
32 
33  use m_special_funcs, only : sbf8
34  use m_efield, only : efield_type
35 
36  use m_pawang, only : pawang_type
37  use m_pawcprj, only : pawcprj_type
38  use m_pawrad, only : pawrad_type,simp_gen
39  use m_pawtab, only : pawtab_type
40  use m_paw_sphharm, only : initylmr
41  use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free, pawcprj_getdim
42 
43  implicit none
44 
45  private
46 
47 !public procedures.
48  public :: overlap_k1k2_paw
49  public :: smatrix_pawinit
50  public :: smatrix_k_paw
51  public :: qijb_kk
52  public :: expibi
53 
54 CONTAINS  !========================================================================================

m_paw_overlap/expibi [ Functions ]

[ Top ] [ m_paw_overlap ] [ Functions ]

NAME

 expibi

FUNCTION

 Routine that computes exp(i (-b_ket).R) at each site.

INPUTS

  dkvecs(3) :: $\Delta k$ increment
  natom :: number of atoms in unit cell
  xred(natom,3) :: reduced coordinates of atoms in unit cell

OUTPUT

  calc_expibi(2,natom) :: phase factors at each atom for vector shift

SIDE EFFECTS

NOTES

PARENTS

      initberry,overlap_k1k2_paw

CHILDREN

SOURCE

909  subroutine expibi(calc_expibi,dkvecs,natom,xred)
910 
911 
912 !This section has been created automatically by the script Abilint (TD).
913 !Do not modify the following lines by hand.
914 #undef ABI_FUNC
915 #define ABI_FUNC 'expibi'
916 !End of the abilint section
917 
918  implicit none
919 
920 !Arguments---------------------------
921 !scalars
922  integer,intent(in) :: natom
923  real(dp),intent(out) :: calc_expibi(2,natom)
924 !arrays
925  real(dp),intent(in) :: dkvecs(3),xred(3,natom)
926 
927 !Local variables---------------------------
928 !scalars
929  integer :: iatom
930  real(dp) :: bdotr
931 
932 ! *************************************************************************
933 
934  calc_expibi(:,:) = zero
935 
936  !calc_expibi(2,natom)
937  !used for PAW field calculations (distributed over atomic sites)
938  !stores the on-site phase factors arising from
939  !$\langle\phi_{i,k}|\phi_{j,k+\sigma_k k_k}\rangle$
940  !where $\sigma = \pm 1$. These overlaps arise in various Berry
941  !phase calculations of electric and magnetic polarization. The on-site
942  !phase factor is $\exp[-i\sigma_k k_k)\cdot I]$ where
943  !$I$ is the nuclear position.
944 
945  do iatom = 1, natom
946 
947     !    note the definition used for the k-dependence of the PAW basis functions:
948     !$|\phi_{i,k}\rangle = exp(-i k\cdot r)|\phi_i\rangle
949     !    see Umari, Gonze, and Pasquarello, PRB 69,235102 [[cite:Umari2004]] Eq. 23.
950    bdotr = DOT_PRODUCT(xred(1:3,iatom),-dkvecs(1:3))
951     !    here is exp(i b.R) for the given site
952    calc_expibi(1,iatom) = cos(two_pi*bdotr)
953    calc_expibi(2,iatom) = sin(two_pi*bdotr)
954 
955  end do ! end loop on natom
956 
957  end subroutine expibi

m_paw_overlap/overlap_k1k2_paw [ Functions ]

[ Top ] [ m_paw_overlap ] [ Functions ]

NAME

 overlap_k1k2_paw

FUNCTION

 compute PAW overlap between two k points,
 similar to smatrix_k_paw.F90 but more generic

INPUTS

  cprj_k1 (pawcprj_type) :: cprj for occupied bands at point k1
  cprj_k2 :: cprj for occupied bands at point k2
  dk(3) :: vector k2 - k1
  gprimd(3,3)=dimensioned primitive translations of reciprocal lattice
  lmn2max :: lmnmax*(lmnmax+1)/2
  lmnsize(ntypat) :: lmnsize for each atom type
  mband :: number of bands
  natom=number of atoms in unit cell
  nspinor :: number of spinors (1 or 2)
  ntypat=number of types of atoms in unit cell
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  typat=typat(natom) list of atom types
  xred(natom,3) :: locations of atoms in cell

OUTPUT

 k1k2_paw(2,mband,mband) :: array of the on-site PAW parts of the overlaps between Bloch states at points
   k1 and k2, for the various pairs of bands, that is, the on-site part of
   <u_nk1|u_mk2>

SIDE EFFECTS

NOTES

 This routine assumes that the cprj are not explicitly ordered by
 atom type.

PARENTS

      chern_number

CHILDREN

      expibi,qijb_kk

SOURCE

103  subroutine overlap_k1k2_paw(cprj_k1,cprj_k2,dk,gprimd,k1k2_paw,lmn2max,lmnsize,mband,&
104 &                           natom,nspinor,ntypat,pawang,pawrad,pawtab,typat,xred)
105 
106 
107 !This section has been created automatically by the script Abilint (TD).
108 !Do not modify the following lines by hand.
109 #undef ABI_FUNC
110 #define ABI_FUNC 'overlap_k1k2_paw'
111 !End of the abilint section
112 
113  implicit none
114 
115 !Arguments---------------------------
116 !scalars
117  integer,intent(in) :: lmn2max,mband,natom,nspinor,ntypat
118  type(pawang_type),intent(in) :: pawang
119  type(pawcprj_type),intent(in) :: cprj_k1(natom,mband),cprj_k2(natom,mband)
120 
121 !arrays
122  integer,intent(in) :: lmnsize(ntypat),typat(natom)
123  real(dp),intent(in) :: dk(3),gprimd(3,3),xred(natom,3)
124  real(dp),intent(out) :: k1k2_paw(2,mband,mband)
125  type(pawrad_type),intent(in) :: pawrad(ntypat)
126  type(pawtab_type),intent(in) :: pawtab(ntypat)
127 
128 !Local variables---------------------------
129 !scalars
130  integer :: iatom,iband,ibs,ilmn,ispinor,itypat
131  integer :: jband,jbs,jlmn,klmn
132  complex(dpc) :: cpk1,cpk2,cterm,paw_onsite
133 
134  ! arrays
135  real(dp),allocatable :: calc_expibi(:,:),calc_qijb(:,:,:)
136 
137 ! *************************************************************************
138 
139 !initialize k1k2_paw output variable
140  k1k2_paw(:,:,:) = zero
141 
142  ! obtain the atomic phase factors for the input k vector shift
143  ABI_ALLOCATE(calc_expibi,(2,natom))
144  call expibi(calc_expibi,dk,natom,xred)
145 
146  ! obtain the onsite PAW terms for the input k vector shift
147  ABI_ALLOCATE(calc_qijb,(2,lmn2max,natom))
148  call qijb_kk(calc_qijb,dk,calc_expibi,gprimd,lmn2max,natom,ntypat,pawang,pawrad,pawtab,typat)
149  ABI_DEALLOCATE(calc_expibi)
150 
151  do iatom = 1, natom
152    itypat = typat(iatom)
153 
154    do ilmn=1,lmnsize(itypat)
155      do jlmn=1,lmnsize(itypat)
156        klmn=max(ilmn,jlmn)*(max(ilmn,jlmn)-1)/2 + min(ilmn,jlmn)
157        paw_onsite = cmplx(calc_qijb(1,klmn,iatom),calc_qijb(2,klmn,iatom))
158        do iband = 1, mband
159          do jband = 1, mband
160            do ispinor = 1, nspinor
161              ibs = nspinor*(iband-1) + ispinor
162              jbs = nspinor*(jband-1) + ispinor
163              cpk1=cmplx(cprj_k1(iatom,ibs)%cp(1,ilmn),cprj_k1(iatom,ibs)%cp(2,ilmn))
164              cpk2=cmplx(cprj_k2(iatom,jbs)%cp(1,jlmn),cprj_k2(iatom,jbs)%cp(2,jlmn))
165              cterm = conjg(cpk1)*paw_onsite*cpk2
166              k1k2_paw(1,iband,jband) = k1k2_paw(1,iband,jband)+real(cterm)
167              k1k2_paw(2,iband,jband) = k1k2_paw(2,iband,jband)+aimag(cterm)
168            end do ! end loop over ispinor
169          end do ! end loop over jband
170        end do ! end loop over iband
171      end do ! end loop over ilmn
172    end do ! end loop over jlmn
173 
174  end do ! end loop over atoms
175 
176  ABI_DEALLOCATE(calc_qijb)
177 
178  end subroutine overlap_k1k2_paw

m_paw_overlap/qijb_kk [ Functions ]

[ Top ] [ m_paw_overlap ] [ Functions ]

NAME

 qijb_kk

FUNCTION

 Routine which computes PAW onsite part of wavefunction overlap for Bloch
 functions at two k-points k and k+b. These
 quantities are used in PAW-based computations of polarization and magnetization.

INPUTS

  dkvecs(3) :: $\Delta k$ input vector
  expibi(2,my_natom,3) :: phase factors at each atomic site for given k offset
  gprimd(3,3)=dimensioned primitive translations of reciprocal lattice
  lmn2max :: lmnmax*(lmnmax+1)/2
  natom=number of atoms in unit cell
  ntypat=number of types of atoms in unit cell
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  typat=typat(natom) list of atom types

OUTPUT

  calc_qijb(2,lmn2max,natom) :: PAW on-site overlaps of wavefunctions at neighboring
                                   k point

SIDE EFFECTS

NOTES

 this function computes the on-site data for the PAW version of
 <u_nk|u_mk+b>, that is, two Bloch vectors at two different k points.

PARENTS

      initberry,overlap_k1k2_paw

CHILDREN

      initylmr,sbf8,simp_gen

SOURCE

759  subroutine qijb_kk(calc_qijb,dkvecs,expibi,gprimd,lmn2max,natom,ntypat,&
760 &                   pawang,pawrad,pawtab,typat)
761 
762 
763 !This section has been created automatically by the script Abilint (TD).
764 !Do not modify the following lines by hand.
765 #undef ABI_FUNC
766 #define ABI_FUNC 'qijb_kk'
767 !End of the abilint section
768 
769  implicit none
770 
771 !Arguments---------------------------
772 !scalars
773  integer,intent(in) :: lmn2max,natom,ntypat
774  type(pawang_type),intent(in) :: pawang
775  real(dp),intent(out) :: calc_qijb(2,lmn2max,natom)
776 !arrays
777  integer,intent(in) :: typat(natom)
778  real(dp),intent(in) :: dkvecs(3),expibi(2,natom),gprimd(3,3)
779  type(pawrad_type),intent(in) :: pawrad(ntypat)
780  type(pawtab_type),intent(in) :: pawtab(ntypat)
781 
782 !Local variables---------------------------
783 !scalars
784  integer :: iatom,ir,isel,itypat
785  integer :: klm,kln,klmn,lbess,lbesslm,lmin,lmax,mbess,mesh_size
786  integer :: ylmr_normchoice,ylmr_npts,ylmr_option
787  real(dp) :: arg,bessg,bnorm,intg,rterm
788  complex(dpc) :: cterm,etb,ifac
789 !arrays
790  real(dp) :: bb(3),bbn(3),bcart(3),ylmgr(1,1,0),ylmr_nrm(1)
791  real(dp),allocatable :: ff(:),j_bessel(:,:),ylmb(:),sb_out(:)
792 ! the following is (i)^L mod 4.
793  complex(dpc),dimension(0:3) :: il(0:3)=(/cone,j_dpc,-cone,-j_dpc/)
794 
795 ! *************************************************************************
796 
797  calc_qijb(:,:,:) = zero
798 
799  ylmr_normchoice = 0 ! input to initylmr are normalized
800  ylmr_npts = 1 ! only 1 point to compute in initylmr
801  ylmr_nrm(1) = one ! weight of normed point for initylmr
802  ylmr_option = 1 ! compute only ylm's in initylmr
803 
804  ABI_ALLOCATE(sb_out, (pawang%l_size_max))
805 
806  do iatom = 1, natom
807 
808    itypat = typat(iatom)
809    mesh_size = pawtab(itypat)%mesh_size
810 
811    ABI_ALLOCATE(j_bessel,(mesh_size,pawang%l_size_max))
812    ABI_ALLOCATE(ff,(mesh_size))
813    ABI_ALLOCATE(ylmb,(pawang%l_size_max*pawang%l_size_max))
814 
815    !    here is exp(-i b.R) for current atom: recall storage in expibi
816    etb = cmplx(expibi(1,iatom),expibi(2,iatom))
817 
818    !    note the definition used for the k-dependence of the PAW basis functions:
819    !$|\phi_{i,k}\rangle = exp(-i k\cdot r)|\phi_i\rangle
820    !    see Umari, Gonze, and Pasquarello, PRB 69,235102 [[cite:Umari2004]], Eq. 23. Thus the k-vector on the
821    !    bra side enters as k, while on the ket side it enters as -k.
822    bb(:) = -dkvecs(:)
823 
824    !    reference bb to cartesian axes
825    bcart(1:3)=MATMUL(gprimd(1:3,1:3),bb(1:3))
826 
827    !    bbn is b-hat (the unit vector in the b direction)
828    bnorm=dsqrt(dot_product(bcart,bcart))
829    bbn(:) = bcart(:)/bnorm
830 
831    !    as an argument to the bessel function, need 2pi*b*r = 1 so b is re-normed to two_pi
832    bnorm = two_pi*bnorm
833    do ir=1,mesh_size
834      arg=bnorm*pawrad(itypat)%rad(ir)
835      call sbf8(pawang%l_size_max,arg,sb_out) ! spherical bessel functions at each mesh point
836      j_bessel(ir,:) = sb_out
837    end do ! end loop over mesh
838 
839    !    compute Y_LM(b) here
840    call initylmr(pawang%l_size_max,ylmr_normchoice,ylmr_npts,ylmr_nrm,ylmr_option,bbn,ylmb(:),ylmgr)
841 
842    do klmn = 1, pawtab(itypat)%lmn2_size
843      klm =pawtab(itypat)%indklmn(1,klmn)
844      kln =pawtab(itypat)%indklmn(2,klmn)
845      lmin=pawtab(itypat)%indklmn(3,klmn)
846      lmax=pawtab(itypat)%indklmn(4,klmn)
847      do lbess = lmin, lmax, 2    ! only possible choices for L s.t. Gaunt integrals
848                                   !        will be non-zero
849        ifac = il(mod(lbess,4))
850        do mbess = -lbess, lbess
851          lbesslm = lbess*lbess+lbess+mbess+1
852          isel=pawang%gntselect(lbesslm,klm)
853          if (isel > 0) then
854            bessg = pawang%realgnt(isel)
855            ff(1:mesh_size)=(pawtab(itypat)%phiphj(1:mesh_size,kln)&
856 &           -pawtab(itypat)%tphitphj(1:mesh_size,kln))&
857 &           *j_bessel(1:mesh_size,lbess+1)
858            call simp_gen(intg,ff,pawrad(itypat))
859            rterm = four_pi*bessg*intg*ylmb(lbesslm)
860            cterm = etb*ifac*rterm
861            calc_qijb(1,klmn,iatom) = &
862 &           calc_qijb(1,klmn,iatom) + dreal(cterm)
863            calc_qijb(2,klmn,iatom) = &
864 &           calc_qijb(2,klmn,iatom) + dimag(cterm)
865 
866          end if ! end selection on non-zero Gaunt factors
867        end do ! end loop on mbess = -lbess, lbess
868      end do ! end loop on lmin-lmax bessel l values
869    end do ! end loop on lmn2_size klmn basis pairs
870 
871    ABI_DEALLOCATE(j_bessel)
872    ABI_DEALLOCATE(ff)
873    ABI_DEALLOCATE(ylmb)
874  end do ! end loop over atoms
875 
876  ABI_DEALLOCATE(sb_out)
877 
878  end subroutine qijb_kk

m_paw_overlap/smatrix_k_paw [ Functions ]

[ Top ] [ m_paw_overlap ] [ Functions ]

NAME

 smatrix_k_paw

FUNCTION

INPUTS

  cprj_k (pawcprj_type) :: cprj for occupied bands at point k
  cprj_kb :: cprj for occupied bands at point k+b
  dtefield :: structure referring to all efield and berry's phase variables
  kdir :: integer giving direction along which overlap is computed for ket
  kfor :: integer indicating whether to compute forward (1) or backward (2)
    along kpt string
  natom :: number of atoms in cell
  typat :: typat(natom) type of each atom

OUTPUT

 smat_k_paw :: array of the on-site PAW parts of the overlaps between Bloch states at points
   k and k+b, for the various pairs of bands, that is, the on-site part of
   <u_nk|u_mk+b>

SIDE EFFECTS

NOTES

 This routine assumes that the cprj are not explicitly ordered by
 atom type.

PARENTS

      berryphase_new,cgwf,make_grad_berry

CHILDREN

SOURCE

653  subroutine smatrix_k_paw(cprj_k,cprj_kb,dtefield,kdir,kfor,mband,natom,smat_k_paw,typat)
654 
655 
656 !This section has been created automatically by the script Abilint (TD).
657 !Do not modify the following lines by hand.
658 #undef ABI_FUNC
659 #define ABI_FUNC 'smatrix_k_paw'
660 !End of the abilint section
661 
662  implicit none
663 
664 !Arguments---------------------------
665 !scalars
666  integer,intent(in) :: kdir,kfor,mband,natom
667  type(efield_type),intent(in) :: dtefield
668  type(pawcprj_type),intent(in) :: cprj_k(natom,dtefield%nspinor*mband)
669  type(pawcprj_type),intent(in) :: cprj_kb(natom,dtefield%nspinor*mband)
670 
671 !arrays
672  integer,intent(in) :: typat(natom)
673  real(dp),intent(out) :: smat_k_paw(2,dtefield%mband_occ,dtefield%mband_occ)
674 
675 !Local variables---------------------------
676 !scalars
677  integer :: iatom,iband,ibs,ilmn,ispinor,itypat
678  integer :: jband,jbs,jlmn,klmn,nspinor
679  complex(dpc) :: cpk,cpkb,cterm,paw_onsite
680 
681 ! *************************************************************************
682 
683 !initialize smat_k_paw
684  smat_k_paw(:,:,:) = zero
685 
686  nspinor = dtefield%nspinor
687 
688  do iatom = 1, natom
689    itypat = typat(iatom)
690 
691    do ilmn=1,dtefield%lmn_size(itypat)
692      do jlmn=1,dtefield%lmn_size(itypat)
693        klmn=max(ilmn,jlmn)*(max(ilmn,jlmn)-1)/2 + min(ilmn,jlmn)
694        paw_onsite = cmplx(dtefield%qijb_kk(1,klmn,iatom,kdir),&
695 &       dtefield%qijb_kk(2,klmn,iatom,kdir))
696        if (kfor > 1) paw_onsite = conjg(paw_onsite)
697        do iband = 1, dtefield%mband_occ
698          do jband = 1, dtefield%mband_occ
699            do ispinor = 1, nspinor
700              ibs = nspinor*(iband-1) + ispinor
701              jbs = nspinor*(jband-1) + ispinor
702              cpk=cmplx(cprj_k(iatom,ibs)%cp(1,ilmn),cprj_k(iatom,ibs)%cp(2,ilmn))
703              cpkb=cmplx(cprj_kb(iatom,jbs)%cp(1,jlmn),cprj_kb(iatom,jbs)%cp(2,jlmn))
704              cterm = conjg(cpk)*paw_onsite*cpkb
705              smat_k_paw(1,iband,jband) = smat_k_paw(1,iband,jband)+dreal(cterm)
706              smat_k_paw(2,iband,jband) = smat_k_paw(2,iband,jband)+dimag(cterm)
707            end do ! end loop over ispinor
708          end do ! end loop over jband
709        end do ! end loop over iband
710      end do ! end loop over ilmn
711    end do ! end loop over jlmn
712 
713  end do ! end loop over atoms
714 
715  end subroutine smatrix_k_paw

m_paw_overlap/smatrix_pawinit [ Functions ]

[ Top ] [ m_paw_overlap ] [ Functions ]

NAME

 smatrix_pawinit

FUNCTION

 Routine which computes paw part of the overlap used to compute LMWF wannier
  functions and berryphase

INPUTS

  atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
  cprj(natom,nspinor*mband*mkmem*nsppol)= <p_lmn|Cnk> coefficients for each WF |Cnk>
                                          and each |p_lmn> non-local projector
  dimcprj(natom)=array of dimensions of array cprj (not ordered)
  g1(3)= reciprocal vector to put k1+b inside the BZ. bb=k2-k1=b-G1
  ("b" is the true b, so we have to correct bb with G1).
  gprimd(3,3)=dimensional reciprocal space primitive translations
  ikpt1(3)=cartesian coordinates of k1
  ikpt2(3)=cartesian coordinates of k2
  isppol  = spin polarization
  mband=maximum number of bands
  mkmem =number of k points treated by this node.
  mpi_enreg=information about MPI parallelization
  natom=number of atoms in cell.
  nkpt=number of k points.
  nspinor=number of spinorial components of the wavefunctions (on current proc)
  nsppol=1 for unpolarized, 2 for spin-polarized
  ntypat=number of types of atoms in unit cell.
  rprimd(3,3)=dimensional primitive translations for real space (bohr)
  seed_name= seed_name of files containing cg for all k-points to be used with MPI
  xred(3,natom)=reduced dimensionless atomic coordinates

OUTPUT

  cm2: Inside sphere part of the overlap needed for constructing wannier function

SIDE EFFECTS

  (only writing, printing)

NOTES

  The mpi part will work with mlwfovlp but not for berryphase_new

PARENTS

      mlwfovlp

CHILDREN

      initylmr,pawcprj_alloc,pawcprj_free,pawcprj_get,pawcprj_getdim,sbf8
      simp_gen

SOURCE

232  subroutine smatrix_pawinit(atindx1,cm2,cprj,ikpt1,ikpt2,isppol,&
233 & g1,gprimd,kpt,mband,mbandw,mkmem,mpi_enreg,&
234 & natom,nband,nkpt,nspinor,nsppol,ntypat,pawang,pawrad,pawtab,rprimd,&
235 & seed_name,typat,xred)
236 
237 
238 !This section has been created automatically by the script Abilint (TD).
239 !Do not modify the following lines by hand.
240 #undef ABI_FUNC
241 #define ABI_FUNC 'smatrix_pawinit'
242 !End of the abilint section
243 
244  implicit none
245 
246 !Arguments---------------------------
247 !scalars
248  integer,intent(in) :: ikpt1,ikpt2,isppol,mband,mbandw,mkmem,natom,nkpt,nspinor,nsppol
249  integer,intent(in) :: ntypat
250  character(len=fnlen) ::  seed_name  !seed names of files containing cg info used in case of MPI
251  type(MPI_type),intent(in) :: mpi_enreg
252  type(pawang_type),intent(in) :: pawang
253 
254 !arrays
255  integer,intent(in) :: atindx1(natom),g1(3),nband(nsppol*nkpt),typat(natom)
256  real(dp),intent(in) :: gprimd(3,3),kpt(3,nkpt),rprimd(3,3),xred(3,natom)
257  real(dp),intent(inout) :: cm2(2,mbandw,mbandw)
258  type(pawcprj_type) :: cprj(natom,nspinor*mband*mkmem*nsppol)
259  type(pawrad_type),intent(in) :: pawrad(ntypat)
260  type(pawtab_type),intent(in) :: pawtab(ntypat)
261 
262 !Local variables---------------------------
263 !scalars
264  integer :: dummy
265  integer :: iatom,iband1,iband2,icg1,icg2,idx1,idx2,ii
266  integer :: ilmn,ios,iunit,ir
267  integer :: iorder_cprj,isel,ispinor,itypat,j0lmn,jj,jlmn,klm,klmn,kln,ll,lm0,lmax
268  integer :: lmin,lmn_size,max_lmn,mesh_size,mm,nband_k
269  integer :: nprocs,spaceComm,rank !for mpi
270  real(dp) :: arg,bnorm,delta,intg,ppi,ppr,qijbtemp,qijtot,x1
271  real(dp) :: x2,xsum,xtemp,xx,yy,zz
272  character(len=500) :: message
273  character(len=fnlen) ::  cprj_file  !file containing cg info used in case of MPI
274  logical::lfile
275 
276 !arrays
277  integer,allocatable :: dimcprj(:),nattyp_dum(:)
278  real(dp),parameter :: ili(7)=(/zero,-one,zero,one,zero,-one,zero/)
279  real(dp),parameter :: ilr(7)=(/one,zero,-one,zero,one,zero,-one/)
280  real(dp) :: bb(3),bb1(3),bbn(3),qijb(2),xcart(3,natom)
281  real(dp),allocatable :: ff(:),j_bessel(:,:),ylmb(:),ylmrgr_dum(:,:,:)
282  real(dp),allocatable :: sb_out(:)
283  type(pawcprj_type),allocatable :: cprj_k1(:,:)
284  type(pawcprj_type),allocatable :: cprj_k2(:,:)
285 
286 ! *************************************************************************
287 
288  DBG_ENTER("COLL")
289 
290 !
291 !Allocate cprj_k1 and cprj_k2
292 !
293  ABI_ALLOCATE(dimcprj,(natom))
294  call pawcprj_getdim(dimcprj,natom,nattyp_dum,ntypat,typat,pawtab,'R')
295 
296  nband_k=nband(ikpt1)
297  ABI_DATATYPE_ALLOCATE(cprj_k1,(natom,nband_k*nspinor))
298  call pawcprj_alloc(cprj_k1,0,dimcprj)
299 
300  nband_k=nband(ikpt2)
301  ABI_DATATYPE_ALLOCATE(cprj_k2,(natom,nband_k*nspinor))
302  call pawcprj_alloc(cprj_k2,0,dimcprj)
303  ABI_DEALLOCATE(dimcprj)
304 
305 !mpi initialization
306  spaceComm=MPI_enreg%comm_cell
307  nprocs=xmpi_comm_size(spaceComm)
308  rank=MPI_enreg%me_kpt
309 
310  lfile=.false.
311 !
312 !write(std_out,*) "compute PAW overlap for k-points",ikpt1,ikpt2
313  do iatom=1,natom
314    xcart(:,iatom)=rprimd(:,1)*xred(1,iatom)+&
315 &   rprimd(:,2)*xred(2,iatom)+&
316 &   rprimd(:,3)*xred(3,iatom)
317  end do
318 
319 !
320 !Calculate indices icg1 and icg2
321 !
322  icg1=0
323  do ii=1,isppol
324    ll=nkpt
325    if(ii==isppol) ll=ikpt1-1
326    do jj=1,ll
327 !    MPI: cycle over kpts not treated by this node
328      if ( ABS(MPI_enreg%proc_distrb(jj,1,ii)-rank)/=0) CYCLE
329 !    write(std_out,'("kpt loop2: ikpt",i3," rank ",i3)') jj,rank
330      icg1=icg1+nspinor*nband(jj+(ii-1)*nkpt)
331    end do
332  end do
333  icg2=0
334  do ii=1,isppol
335    ll=nkpt
336    if(isppol==ii) ll=ikpt2-1
337    do jj=1,ll
338 !    MPI: cycle over kpts not treated by this node
339      if (ABS(MPI_enreg%proc_distrb(jj,1,ii)-rank)/=0) CYCLE
340 !    write(std_out,'("kpt loop2: ikpt",i3," rank ",i3)') jj,rank
341      icg2=icg2+nspinor*nband(jj+(ii-1)*nkpt)
342    end do
343  end do
344 !
345 !MPI: if ikpt2 not found in this processor then
346 !read info from an unformatted file
347 !
348  if (nprocs>1) then
349    if (ABS(MPI_enreg%proc_distrb(ikpt2,1,isppol)-rank)/=0) then
350      lfile=.true.
351 !
352 !    get maximum of lmn_size
353      max_lmn=0
354      do itypat=1,ntypat
355        lmn_size=pawtab(itypat)%lmn_size
356        if(lmn_size>max_lmn) max_lmn=lmn_size
357      end do
358 !
359 !    get file name and open it
360 !
361      write(cprj_file,'(a,I5.5,".",I1)') trim(seed_name),ikpt2,isppol
362      iunit=1000
363 !    write(std_out,*)'reading file',trim(cprj_file)
364      open (unit=iunit, file=cprj_file,form='unformatted',status='old',iostat=ios)
365      if(ios /= 0) then
366        write(message,*) " smatrix_pawinit: file",trim(cprj_file), "not found"
367        MSG_ERROR(message)
368      end if
369 !
370 !    start reading
371      do ii=1,mband*nspinor
372        do iatom=1,natom
373          itypat=typat(iatom)
374          lmn_size=pawtab(itypat)%lmn_size
375          do ilmn=1,lmn_size
376            read(iunit)(cprj_k2(iatom,ii)%cp(jj,ilmn),jj=1,2)
377          end do !ilmn
378        end do
379      end do
380 !
381 !    close file
382 !
383      close (unit=iunit,iostat=ios)
384      if(ios /= 0) then
385        write(message,*) " smatrix_pawinit: error closing file ",trim(cprj_file)
386        MSG_ERROR(message)
387      end if
388 !
389    end if
390  end if !mpi
391 
392 !Extract cprj_k1 and cprj_k2
393 !these contain the projectors cprj for just one k-point (ikpt1 or ikpt2)
394 
395 !Extract cprj for k-point 1
396  iorder_cprj=0 !do not change the ordering of cprj
397  nband_k=nband(ikpt1)
398  dummy=1000 !index of file not implemented here, mkmem==0 not implemented
399  call pawcprj_get(atindx1,cprj_k1,cprj,natom,1,icg1,ikpt1,iorder_cprj,isppol,&
400 & mband,mkmem,natom,nband_k,nband_k,nspinor,nsppol,dummy,&
401 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
402 
403 !Extract cprj for k-point 2
404  if( lfile .eqv. .false. ) then !if it was not already read above
405    iorder_cprj=0 !do not change the ordering of cprj
406    nband_k=nband(ikpt2)
407    dummy=1000 !index of file not implemented here, mkmem==0 not implemented
408    call pawcprj_get(atindx1,cprj_k2,cprj,natom,1,icg2,ikpt2,iorder_cprj,isppol,&
409 &   mband,mkmem,natom,nband_k,nband_k,nspinor,nsppol,dummy,&
410 &   mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
411  end if
412 
413 !DEBUG
414 !if(ikpt2==2) then
415 !if(ikpt1==1) then
416 !do iband1=1,mbandw
417 !do iatom=1,natom
418 !itypat=typat(atindx1(iatom))
419 !lmn_size=pawtab(itypat)%lmn_size
420 !do ilmn=1,1!lmn_size
421 !!write(500,'(a,i3,a,i3,a,i3,a,2f13.7)')'iband ',iband1,' iatom ',iatom,' ilmn ',ilmn,' cprj',cprj(iatom,iband1+icg1)%cp(:,ilmn)
422 !write(500,'(a,i3,a,i3,a,i3,a,2f13.7)')'iband ',iband1,' iatom ',atindx1(iatom),' ilmn ',ilmn,' cprj',cprj(atindx1(iatom),iband1+icg1)%cp(:,ilmn)
423 !write(500,'(a,i3,a,i3,a,i3,a,2f13.7)')'iband ',iband1,' iatom ',iatom,' ilmn ',ilmn,' cprj_k1 ',cprj_k1(iatom,iband1)%cp(:,ilmn)
424 !end do
425 !end do
426 !end do
427 !end if
428 !end if
429 !NDDEBUG
430 
431 !!!!!!!!!!!!!!!!!
432 !--- Compute intermediate quantities: "b" vector=k2-k1 and its
433 !normalized value: bbn (and its norm: bnorm)
434 !compute also Ylm(b).
435  ABI_ALLOCATE(ylmb,(pawang%l_size_max*pawang%l_size_max))
436  ABI_ALLOCATE(ylmrgr_dum,(1,1,0))
437  bb(:)=kpt(:,ikpt2)-kpt(:,ikpt1)+g1(:)
438  bb1=bb
439  xx=gprimd(1,1)*bb(1)+gprimd(1,2)*bb(2)+gprimd(1,3)*bb(3)
440  yy=gprimd(2,1)*bb(1)+gprimd(2,2)*bb(2)+gprimd(2,3)*bb(3)
441  zz=gprimd(3,1)*bb(1)+gprimd(3,2)*bb(2)+gprimd(3,3)*bb(3)
442  bnorm=two_pi*dsqrt(xx**2+yy**2+zz**2)
443  if(bnorm<tol8) then
444 !  write(std_out,*) "WARNING: bnorm=",bnorm
445    bbn(:)=zero
446  else
447    xx=xx*two_pi
448    yy=yy*two_pi
449    zz=zz*two_pi
450    bb(1)=xx
451    bb(2)=yy
452    bb(3)=zz
453    bbn(1)=xx/bnorm
454    bbn(2)=yy/bnorm
455    bbn(3)=zz/bnorm
456  end if
457 
458 !debug  bbn=0
459 !debug  bnorm=0
460 !bbn has to ne normalized
461  call initylmr(pawang%l_size_max,0,1,(/one/),1,bbn(:),ylmb(:),ylmrgr_dum)
462 !write(std_out,*) "ylmb(:)",ylmb(:)
463 !write(std_out,*) pawang%l_size_max
464 !write(std_out,*) "bbn",bbn(:)
465 !write(std_out,*) "xx,yy,zz",xx,yy,zz
466 !write(std_out,*) "bnorm",bnorm
467  ABI_DEALLOCATE(ylmrgr_dum)
468 
469 !------- First Compute Qij(b)-
470  ABI_ALLOCATE(sb_out, (pawang%l_size_max))
471  cm2=zero
472 
473  do iatom=1,natom
474    itypat=typat(iatom)
475    lmn_size=pawtab(itypat)%lmn_size
476 !  ---  en coordonnnes reelles cartesiennes (espace reel)
477 !  ---  first radial part(see pawinit)
478    mesh_size=pawtab(itypat)%mesh_size
479    ABI_ALLOCATE(j_bessel,(mesh_size,pawang%l_size_max))
480 
481 
482 !  ---  compute bessel function for (br) for all angular momenta necessary
483 !  ---  and for all value of r.
484 !  ---  they are needed for radial part
485 !  ---  of the integration => j_bessel(ir,:)
486    do ir=1,mesh_size
487      arg=bnorm*pawrad(itypat)%rad(ir)
488      call sbf8(pawang%l_size_max,arg,sb_out)
489      j_bessel(ir,:) = sb_out
490    end do
491 
492 !  do jlmn=1,pawang%l_size_max
493 !  write(665,*) "j_bessel",j_bessel(1:mesh_size,jlmn)
494 !  enddo
495 !  write(std_out,*) "bessel function computed"
496 !  ---  Compute \Sum b.R=xsum for future use
497    xtemp=zero
498    do mm=1,3
499      xtemp=xtemp+xred(mm,iatom)*bb1(mm)
500    end do
501    xtemp=xtemp*two_pi
502    xsum=zero
503    do mm=1,3
504      xsum=xsum+xcart(mm,iatom)*bbn(mm)*bnorm
505    end do
506 !  write(std_out,*)'xsum',xsum,xtemp,lmn_size
507 
508 !  ---  Loop on jlmn and ilmn
509    qijtot=zero
510    do jlmn=1,lmn_size
511      j0lmn=jlmn*(jlmn-1)/2
512      do ilmn=1,jlmn
513 
514        klmn=j0lmn+ilmn
515        klm=pawtab(itypat)%indklmn(1,klmn);kln=pawtab(itypat)%indklmn(2,klmn)
516        lmin=pawtab(itypat)%indklmn(3,klmn);lmax=pawtab(itypat)%indklmn(4,klmn)
517 !      ---  Sum over angular momenta
518 !      ---  compute radial part integration for each angular momentum => intg
519 !      ---  (3j) symbols follows the rule: l belongs to abs(li-lj), li+lj.
520        qijb=zero
521        do ll=lmin,lmax,2
522          lm0=ll*ll+ll+1
523          ABI_ALLOCATE(ff,(mesh_size))
524          ff(1:mesh_size)=(pawtab(itypat)%phiphj(1:mesh_size,kln)&
525 &         -pawtab(itypat)%tphitphj(1:mesh_size,kln))&
526 &         *j_bessel(1:mesh_size,ll+1)
527          call simp_gen(intg,ff,pawrad(itypat))
528          ABI_DEALLOCATE(ff)
529          qijbtemp=zero
530          do mm=-ll,ll
531            isel=pawang%gntselect(lm0+mm,klm)
532            if (isel>0) qijbtemp=qijbtemp&
533 &           +pawang%realgnt(isel)*ylmb(lm0+mm)
534          end do ! mm
535 !        ---     compute angular part with a summation
536 !        ---     qijb =\sum_{lm} intg(lm)*qijbtemp
537          qijb(1)=qijb(1) +intg*qijbtemp*ilr(ll+1)
538          qijb(2)=qijb(2) +intg*qijbtemp*ili(ll+1)
539 !        if(ilmn==jlmn) write(std_out,*) "intg, qij",intg,qijbtemp
540        end do ! ll
541 
542 !      ---  Add exp(-i.b*R) for each atom.
543        if(ilmn==jlmn) qijtot=qijtot+qijb(1)
544 !      if(ilmn==jlmn) write(std_out,*) "qijtot",qijtot
545        x1=qijb(1)*dcos(-xsum)-qijb(2)*dsin(-xsum)
546        x2=qijb(1)*dsin(-xsum)+qijb(2)*dcos(-xsum)
547 !      x1 x2 necessary to avoid changing qijb(1) before
548 !      computing qijb(2)
549        qijb(1)=x1
550        qijb(2)=x2 !
551 !      if(ilmn==jlmn) write(std_out,*) "qij",jlmn,ilmn,qijb(1),qijb(2)
552 
553        do iband1=1,mbandw ! limite inferieure a preciser
554          do iband2=1,mbandw
555            ppr=0.d0
556            ppi=0.d0
557            do ispinor=1,nspinor
558              idx1=iband1*nspinor-(nspinor-ispinor)
559              idx2=iband2*nspinor-(nspinor-ispinor) !to take into account spinors
560 !            write(std_out,*) "iband2",iband2
561 !            product of (a1+ia2)*(b1-ib2) (minus sign because conjugated)
562              ppr=ppr+&
563 !            real part a_1*b_1+a_2*b_2
564 &             cprj_k1(iatom,idx1)%cp(1,ilmn)*cprj_k2(iatom,idx2)%cp(1,jlmn)+&
565 &             cprj_k1(iatom,idx1)%cp(2,ilmn)*cprj_k2(iatom,idx2)%cp(2,jlmn)+&
566 !            &     cprj(iatom,idx1+icg1)%cp(1,ilmn)*cprj(iatom,idx2+icg2)%cp(1,jlmn)+&
567 !            &     cprj(iatom,idx1+icg1)%cp(2,ilmn)*cprj(iatom,idx2+icg2)%cp(2,jlmn)+&
568 !            add term on the other triangle  of the matrix
569 !            qij is the same for this part because phi are real.
570 &             cprj_k1(iatom,idx1)%cp(1,jlmn)*cprj_k2(iatom,idx2)%cp(1,ilmn)+&
571 &             cprj_k1(iatom,idx1)%cp(2,jlmn)*cprj_k2(iatom,idx2)%cp(2,ilmn)
572 !            &     cprj(iatom,idx1+icg1)%cp(1,jlmn)*cprj(iatom,idx2+icg2)%cp(1,ilmn)+&
573 !            &     cprj(iatom,idx1+icg1)%cp(2,jlmn)*cprj(iatom,idx2+icg2)%cp(2,ilmn)
574              ppi=ppi+&
575 !            imaginary part a_1*b_2-a_2*b_1
576 &             cprj_k1(iatom,idx1)%cp(1,ilmn)*cprj_k2(iatom,idx2)%cp(2,jlmn)-&
577 &             cprj_k1(iatom,idx1)%cp(2,ilmn)*cprj_k2(iatom,idx2)%cp(1,jlmn)+&
578 !            &     cprj(iatom,idx1+icg1)%cp(1,ilmn)*cprj(iatom,idx2+icg2)%cp(2,jlmn)-&
579 !            &     cprj(iatom,idx1+icg1)%cp(2,ilmn)*cprj(iatom,idx2+icg2)%cp(1,jlmn)+&
580 !            add term on the other triangle  of the matrix
581 &             cprj_k1(iatom,idx1)%cp(1,jlmn)*cprj_k2(iatom,idx2)%cp(2,ilmn)-&
582 &             cprj_k1(iatom,idx1)%cp(2,jlmn)*cprj_k2(iatom,idx2)%cp(1,ilmn)
583 !            &     cprj(iatom,idx1+icg1)%cp(1,jlmn)*cprj(iatom,idx2+icg2)%cp(2,ilmn)-&
584 !            &     cprj(iatom,idx1+icg1)%cp(2,jlmn)*cprj(iatom,idx2+icg2)%cp(1,ilmn)
585            end do !ispinor
586 !
587 !          delta: diagonal terms are counted twice ! so
588 !          we need a 0.5 factor for diagonal elements.
589            delta=one
590 !          write(std_out,*) "ppr and ppi computed",ikpt1,ikpt2,iband1,iband2
591            if(ilmn==jlmn) delta=half
592            cm2(1,iband1,iband2)= cm2(1,iband1,iband2)+ &
593 &           (qijb(1)*ppr-qijb(2)*ppi)*delta
594            cm2(2,iband1,iband2)= cm2(2,iband1,iband2)+ &
595 &           (qijb(2)*ppr+qijb(1)*ppi)*delta
596          end do ! iband2
597        end do ! iband1
598 
599      end do ! ilmn
600    end do ! jlmn
601 !  write(std_out,*) "final qijtot",qijtot
602    ABI_DEALLOCATE(j_bessel)
603  end do ! iatom
604 
605  ABI_DEALLOCATE(sb_out)
606  ABI_DEALLOCATE(ylmb)
607  call pawcprj_free(cprj_k1)
608  call pawcprj_free(cprj_k2)
609  ABI_DATATYPE_DEALLOCATE(cprj_k1)
610  ABI_DATATYPE_DEALLOCATE(cprj_k2)
611 
612  DBG_EXIT("COLL")
613 
614  end subroutine smatrix_pawinit