TABLE OF CONTENTS
- ABINIT/m_paw_overlap
- m_paw_overlap/expibi
- m_paw_overlap/overlap_k1k2_paw
- m_paw_overlap/qijb_kk
- m_paw_overlap/smatrix_k_paw
- m_paw_overlap/smatrix_pawinit
ABINIT/m_paw_overlap [ 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-2024 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
18 #if defined HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 22 #include "abi_common.h" 23 24 MODULE m_paw_overlap 25 26 use defs_basis 27 use m_errors 28 use m_abicore 29 use m_xmpi 30 31 use defs_abitypes, only : MPI_type 32 use m_special_funcs, only : sbf8 33 use m_efield, only : efield_type 34 use m_pawang, only : pawang_type 35 use m_pawcprj, only : pawcprj_type 36 use m_pawrad, only : pawrad_type,simp_gen 37 use m_pawtab, only : pawtab_type 38 use m_paw_sphharm, only : initylmr 39 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free, pawcprj_getdim 40 41 implicit none 42 43 private 44 45 !public procedures. 46 public :: overlap_k1k2_paw 47 public :: smatrix_pawinit 48 public :: smatrix_k_paw 49 public :: qijb_kk 50 public :: expibi 51 52 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
SOURCE
842 subroutine expibi(calc_expibi,dkvecs,natom,xred) 843 844 !Arguments--------------------------- 845 !scalars 846 integer,intent(in) :: natom 847 real(dp),intent(out) :: calc_expibi(2,natom) 848 !arrays 849 real(dp),intent(in) :: dkvecs(3),xred(3,natom) 850 851 !Local variables--------------------------- 852 !scalars 853 integer :: iatom 854 real(dp) :: bdotr 855 856 ! ************************************************************************* 857 858 calc_expibi(:,:) = zero 859 860 !calc_expibi(2,natom) 861 !used for PAW field calculations (distributed over atomic sites) 862 !stores the on-site phase factors arising from 863 !$\langle\phi_{i,k}|\phi_{j,k+\sigma_k k_k}\rangle$ 864 !where $\sigma = \pm 1$. These overlaps arise in various Berry 865 !phase calculations of electric and magnetic polarization. The on-site 866 !phase factor is $\exp[-i\sigma_k k_k)\cdot I]$ where 867 !$I$ is the nuclear position. 868 869 do iatom = 1, natom 870 871 ! note the definition used for the k-dependence of the PAW basis functions: 872 !$|\phi_{i,k}\rangle = exp(-i k\cdot r)|\phi_i\rangle 873 ! see Umari, Gonze, and Pasquarello, PRB 69,235102 [[cite:Umari2004]] Eq. 23. 874 bdotr = DOT_PRODUCT(xred(1:3,iatom),-dkvecs(1:3)) 875 ! here is exp(i b.R) for the given site 876 calc_expibi(1,iatom) = cos(two_pi*bdotr) 877 calc_expibi(2,iatom) = sin(two_pi*bdotr) 878 879 end do ! end loop on natom 880 881 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.
SOURCE
95 subroutine overlap_k1k2_paw(cprj_k1,cprj_k2,dk,gprimd,k1k2_paw,lmn2max,lmnsize,& 96 & natom,nband,nband_occ,nspinor,ntypat,pawang,pawrad,pawtab,typat,xred) 97 98 !Arguments--------------------------- 99 !scalars 100 integer,intent(in) :: lmn2max,natom,nband,nband_occ,nspinor,ntypat 101 type(pawang_type),intent(in) :: pawang 102 type(pawcprj_type),intent(in) :: cprj_k1(natom,nband),cprj_k2(natom,nband) 103 104 !arrays 105 integer,intent(in) :: lmnsize(ntypat),typat(natom) 106 real(dp),intent(in) :: dk(3),gprimd(3,3),xred(natom,3) 107 real(dp),intent(out) :: k1k2_paw(2,nband_occ,nband_occ) 108 type(pawrad_type),intent(in) :: pawrad(ntypat) 109 type(pawtab_type),intent(in) :: pawtab(ntypat) 110 111 !Local variables--------------------------- 112 !scalars 113 integer :: iatom,iband,ibs,ilmn,ispinor,itypat 114 integer :: jband,jbs,jlmn,klmn 115 complex(dpc) :: cpk1,cpk2,cterm,paw_onsite 116 117 ! arrays 118 real(dp),allocatable :: calc_expibi(:,:),calc_qijb(:,:,:) 119 120 ! ************************************************************************* 121 122 !initialize k1k2_paw output variable 123 k1k2_paw(:,:,:) = zero 124 125 ! obtain the atomic phase factors for the input k vector shift 126 ABI_MALLOC(calc_expibi,(2,natom)) 127 call expibi(calc_expibi,dk,natom,xred) 128 129 ! obtain the onsite PAW terms for the input k vector shift 130 ABI_MALLOC(calc_qijb,(2,lmn2max,natom)) 131 call qijb_kk(calc_qijb,dk,calc_expibi,gprimd,lmn2max,natom,ntypat,pawang,pawrad,pawtab,typat) 132 ABI_FREE(calc_expibi) 133 134 do iatom = 1, natom 135 itypat = typat(iatom) 136 137 do ilmn=1,lmnsize(itypat) 138 do jlmn=1,lmnsize(itypat) 139 klmn=max(ilmn,jlmn)*(max(ilmn,jlmn)-1)/2 + min(ilmn,jlmn) 140 paw_onsite = cmplx(calc_qijb(1,klmn,iatom),calc_qijb(2,klmn,iatom)) 141 do iband = 1, nband_occ 142 do jband = 1, nband_occ 143 do ispinor = 1, nspinor 144 ibs = nspinor*(iband-1) + ispinor 145 jbs = nspinor*(jband-1) + ispinor 146 cpk1=cmplx(cprj_k1(iatom,ibs)%cp(1,ilmn),cprj_k1(iatom,ibs)%cp(2,ilmn)) 147 cpk2=cmplx(cprj_k2(iatom,jbs)%cp(1,jlmn),cprj_k2(iatom,jbs)%cp(2,jlmn)) 148 cterm = conjg(cpk1)*paw_onsite*cpk2 149 k1k2_paw(1,iband,jband) = k1k2_paw(1,iband,jband)+real(cterm) 150 k1k2_paw(2,iband,jband) = k1k2_paw(2,iband,jband)+aimag(cterm) 151 end do ! end loop over ispinor 152 end do ! end loop over jband 153 end do ! end loop over iband 154 end do ! end loop over ilmn 155 end do ! end loop over jlmn 156 157 end do ! end loop over atoms 158 159 ABI_FREE(calc_qijb) 160 161 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.
SOURCE
706 subroutine qijb_kk(calc_qijb,dkvecs,expibi,gprimd,lmn2max,natom,ntypat,& 707 & pawang,pawrad,pawtab,typat) 708 709 !Arguments--------------------------- 710 !scalars 711 integer,intent(in) :: lmn2max,natom,ntypat 712 type(pawang_type),intent(in) :: pawang 713 real(dp),intent(out) :: calc_qijb(2,lmn2max,natom) 714 !arrays 715 integer,intent(in) :: typat(natom) 716 real(dp),intent(in) :: dkvecs(3),expibi(2,natom),gprimd(3,3) 717 type(pawrad_type),intent(in) :: pawrad(ntypat) 718 type(pawtab_type),intent(in) :: pawtab(ntypat) 719 720 !Local variables--------------------------- 721 !scalars 722 integer :: iatom,ir,isel,itypat 723 integer :: klm,kln,klmn,lbess,lbesslm,lmin,lmax,mbess,mesh_size 724 integer :: ylmr_normchoice,ylmr_npts,ylmr_option 725 real(dp) :: arg,bessg,bnorm,intg,rterm 726 complex(dpc) :: cterm,etb,ifac 727 !arrays 728 real(dp) :: bb(3),bbn(3),bcart(3),ylmgr(1,1,0),ylmr_nrm(1) 729 real(dp),allocatable :: ff(:),j_bessel(:,:),ylmb(:),sb_out(:) 730 ! the following is (i)^L mod 4. 731 complex(dpc),dimension(0:3) :: il(0:3)=(/cone,j_dpc,-cone,-j_dpc/) 732 733 ! ************************************************************************* 734 735 calc_qijb(:,:,:) = zero 736 737 ylmr_normchoice = 0 ! input to initylmr are normalized 738 ylmr_npts = 1 ! only 1 point to compute in initylmr 739 ylmr_nrm(1) = one ! weight of normed point for initylmr 740 ylmr_option = 1 ! compute only ylm's in initylmr 741 742 ABI_MALLOC(sb_out, (pawang%l_size_max)) 743 744 do iatom = 1, natom 745 746 itypat = typat(iatom) 747 mesh_size = pawtab(itypat)%mesh_size 748 749 ABI_MALLOC(j_bessel,(mesh_size,pawang%l_size_max)) 750 ABI_MALLOC(ff,(mesh_size)) 751 ABI_MALLOC(ylmb,(pawang%l_size_max*pawang%l_size_max)) 752 753 ! here is exp(-i b.R) for current atom: recall storage in expibi 754 etb = cmplx(expibi(1,iatom),expibi(2,iatom)) 755 756 ! note the definition used for the k-dependence of the PAW basis functions: 757 !$|\phi_{i,k}\rangle = exp(-i k\cdot r)|\phi_i\rangle 758 ! see Umari, Gonze, and Pasquarello, PRB 69,235102 [[cite:Umari2004]], Eq. 23. Thus the k-vector on the 759 ! bra side enters as k, while on the ket side it enters as -k. 760 bb(:) = -dkvecs(:) 761 762 ! reference bb to cartesian axes 763 bcart(1:3)=MATMUL(gprimd(1:3,1:3),bb(1:3)) 764 765 ! bbn is b-hat (the unit vector in the b direction) 766 bnorm=dsqrt(dot_product(bcart,bcart)) 767 bbn(:) = bcart(:)/bnorm 768 769 ! as an argument to the bessel function, need 2pi*b*r = 1 so b is re-normed to two_pi 770 bnorm = two_pi*bnorm 771 do ir=1,mesh_size 772 arg=bnorm*pawrad(itypat)%rad(ir) 773 call sbf8(pawang%l_size_max,arg,sb_out) ! spherical bessel functions at each mesh point 774 j_bessel(ir,:) = sb_out 775 end do ! end loop over mesh 776 777 ! compute Y_LM(b) here 778 call initylmr(pawang%l_size_max,ylmr_normchoice,ylmr_npts,ylmr_nrm,ylmr_option,bbn,ylmb(:),ylmgr) 779 780 do klmn = 1, pawtab(itypat)%lmn2_size 781 klm =pawtab(itypat)%indklmn(1,klmn) 782 kln =pawtab(itypat)%indklmn(2,klmn) 783 lmin=pawtab(itypat)%indklmn(3,klmn) 784 lmax=pawtab(itypat)%indklmn(4,klmn) 785 do lbess = lmin, lmax, 2 ! only possible choices for L s.t. Gaunt integrals 786 ! will be non-zero 787 ifac = il(mod(lbess,4)) 788 do mbess = -lbess, lbess 789 lbesslm = lbess*lbess+lbess+mbess+1 790 isel=pawang%gntselect(lbesslm,klm) 791 if (isel > 0) then 792 bessg = pawang%realgnt(isel) 793 ff(1:mesh_size)=(pawtab(itypat)%phiphj(1:mesh_size,kln)& 794 & -pawtab(itypat)%tphitphj(1:mesh_size,kln))& 795 & *j_bessel(1:mesh_size,lbess+1) 796 call simp_gen(intg,ff,pawrad(itypat)) 797 rterm = four_pi*bessg*intg*ylmb(lbesslm) 798 cterm = etb*ifac*rterm 799 calc_qijb(1,klmn,iatom) = & 800 & calc_qijb(1,klmn,iatom) + dreal(cterm) 801 calc_qijb(2,klmn,iatom) = & 802 & calc_qijb(2,klmn,iatom) + dimag(cterm) 803 804 end if ! end selection on non-zero Gaunt factors 805 end do ! end loop on mbess = -lbess, lbess 806 end do ! end loop on lmin-lmax bessel l values 807 end do ! end loop on lmn2_size klmn basis pairs 808 809 ABI_FREE(j_bessel) 810 ABI_FREE(ff) 811 ABI_FREE(ylmb) 812 end do ! end loop over atoms 813 814 ABI_FREE(sb_out) 815 816 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.
SOURCE
615 subroutine smatrix_k_paw(cprj_k,cprj_kb,dtefield,kdir,kfor,mband,natom,smat_k_paw,typat) 616 617 !Arguments--------------------------- 618 !scalars 619 integer,intent(in) :: kdir,kfor,mband,natom 620 type(efield_type),intent(in) :: dtefield 621 type(pawcprj_type),intent(in) :: cprj_k(natom,dtefield%nspinor*mband) 622 type(pawcprj_type),intent(in) :: cprj_kb(natom,dtefield%nspinor*mband) 623 624 !arrays 625 integer,intent(in) :: typat(natom) 626 real(dp),intent(out) :: smat_k_paw(2,dtefield%mband_occ,dtefield%mband_occ) 627 628 !Local variables--------------------------- 629 !scalars 630 integer :: iatom,iband,ibs,ilmn,ispinor,itypat 631 integer :: jband,jbs,jlmn,klmn,nspinor 632 complex(dpc) :: cpk,cpkb,cterm,paw_onsite 633 634 ! ************************************************************************* 635 636 !initialize smat_k_paw 637 smat_k_paw(:,:,:) = zero 638 639 nspinor = dtefield%nspinor 640 641 do iatom = 1, natom 642 itypat = typat(iatom) 643 644 do ilmn=1,dtefield%lmn_size(itypat) 645 do jlmn=1,dtefield%lmn_size(itypat) 646 klmn=max(ilmn,jlmn)*(max(ilmn,jlmn)-1)/2 + min(ilmn,jlmn) 647 paw_onsite = cmplx(dtefield%qijb_kk(1,klmn,iatom,kdir),& 648 & dtefield%qijb_kk(2,klmn,iatom,kdir)) 649 if (kfor > 1) paw_onsite = conjg(paw_onsite) 650 do iband = 1, dtefield%mband_occ 651 do jband = 1, dtefield%mband_occ 652 do ispinor = 1, nspinor 653 ibs = nspinor*(iband-1) + ispinor 654 jbs = nspinor*(jband-1) + ispinor 655 cpk=cmplx(cprj_k(iatom,ibs)%cp(1,ilmn),cprj_k(iatom,ibs)%cp(2,ilmn)) 656 cpkb=cmplx(cprj_kb(iatom,jbs)%cp(1,jlmn),cprj_kb(iatom,jbs)%cp(2,jlmn)) 657 cterm = conjg(cpk)*paw_onsite*cpkb 658 smat_k_paw(1,iband,jband) = smat_k_paw(1,iband,jband)+dreal(cterm) 659 smat_k_paw(2,iband,jband) = smat_k_paw(2,iband,jband)+dimag(cterm) 660 end do ! end loop over ispinor 661 end do ! end loop over jband 662 end do ! end loop over iband 663 end do ! end loop over ilmn 664 end do ! end loop over jlmn 665 666 end do ! end loop over atoms 667 668 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
SOURCE
208 subroutine smatrix_pawinit(atindx1,cm2,cprj,ikpt1,ikpt2,isppol,& 209 & g1,gprimd,kpt,mband,mbandw,mkmem,mpi_enreg,& 210 & natom,nband,nkpt,nspinor,nsppol,ntypat,pawang,pawrad,pawtab,rprimd,& 211 & seed_name,typat,xred) 212 213 !Arguments--------------------------- 214 !scalars 215 integer,intent(in) :: ikpt1,ikpt2,isppol,mband,mbandw,mkmem,natom,nkpt,nspinor,nsppol 216 integer,intent(in) :: ntypat 217 character(len=fnlen) :: seed_name !seed names of files containing cg info used in case of MPI 218 type(MPI_type),intent(in) :: mpi_enreg 219 type(pawang_type),intent(in) :: pawang 220 221 !arrays 222 integer,intent(in) :: atindx1(natom),g1(3),nband(nsppol*nkpt),typat(natom) 223 real(dp),intent(in) :: gprimd(3,3),kpt(3,nkpt),rprimd(3,3),xred(3,natom) 224 real(dp),intent(inout) :: cm2(2,mbandw,mbandw) 225 type(pawcprj_type) :: cprj(natom,nspinor*mband*mkmem*nsppol) 226 type(pawrad_type),intent(in) :: pawrad(ntypat) 227 type(pawtab_type),intent(in) :: pawtab(ntypat) 228 229 !Local variables--------------------------- 230 !scalars 231 integer :: dummy 232 integer :: iatom,iband1,iband2,icg1,icg2,idx1,idx2,ii 233 integer :: ilmn,ios,iunit,ir 234 integer :: iorder_cprj,isel,ispinor,itypat,j0lmn,jj,jlmn,klm,klmn,kln,ll,lm0,lmax 235 integer :: lmin,lmn_size,max_lmn,mesh_size,mm,nband_k 236 integer :: nprocs,spaceComm,rank !for mpi 237 real(dp) :: arg,bnorm,delta,intg,ppi,ppr,qijbtemp,qijtot,x1 238 real(dp) :: x2,xsum,xtemp,xx,yy,zz 239 character(len=500) :: message 240 character(len=fnlen) :: cprj_file !file containing cg info used in case of MPI 241 logical::lfile 242 243 !arrays 244 integer,allocatable :: dimcprj(:),nattyp_dum(:) 245 real(dp),parameter :: ili(7)=(/zero,-one,zero,one,zero,-one,zero/) 246 real(dp),parameter :: ilr(7)=(/one,zero,-one,zero,one,zero,-one/) 247 real(dp) :: bb(3),bb1(3),bbn(3),qijb(2),xcart(3,natom) 248 real(dp),allocatable :: ff(:),j_bessel(:,:),ylmb(:),ylmrgr_dum(:,:,:) 249 real(dp),allocatable :: sb_out(:) 250 type(pawcprj_type),allocatable :: cprj_k1(:,:) 251 type(pawcprj_type),allocatable :: cprj_k2(:,:) 252 253 ! ************************************************************************* 254 255 DBG_ENTER("COLL") 256 257 ! 258 !Allocate cprj_k1 and cprj_k2 259 ! 260 ABI_MALLOC(dimcprj,(natom)) 261 call pawcprj_getdim(dimcprj,natom,nattyp_dum,ntypat,typat,pawtab,'R') 262 263 nband_k=nband(ikpt1) 264 ABI_MALLOC(cprj_k1,(natom,nband_k*nspinor)) 265 call pawcprj_alloc(cprj_k1,0,dimcprj) 266 267 nband_k=nband(ikpt2) 268 ABI_MALLOC(cprj_k2,(natom,nband_k*nspinor)) 269 call pawcprj_alloc(cprj_k2,0,dimcprj) 270 ABI_FREE(dimcprj) 271 272 !mpi initialization 273 spaceComm=MPI_enreg%comm_cell 274 nprocs=xmpi_comm_size(spaceComm) 275 rank=MPI_enreg%me_kpt 276 277 lfile=.false. 278 ! 279 !write(std_out,*) "compute PAW overlap for k-points",ikpt1,ikpt2 280 do iatom=1,natom 281 xcart(:,iatom)=rprimd(:,1)*xred(1,iatom)+& 282 & rprimd(:,2)*xred(2,iatom)+& 283 & rprimd(:,3)*xred(3,iatom) 284 end do 285 286 ! 287 !Calculate indices icg1 and icg2 288 ! 289 icg1=0 290 do ii=1,isppol 291 ll=nkpt 292 if(ii==isppol) ll=ikpt1-1 293 do jj=1,ll 294 ! MPI: cycle over kpts not treated by this node 295 if ( ABS(MPI_enreg%proc_distrb(jj,1,ii)-rank)/=0) CYCLE 296 ! write(std_out,'("kpt loop2: ikpt",i3," rank ",i3)') jj,rank 297 icg1=icg1+nspinor*nband(jj+(ii-1)*nkpt) 298 end do 299 end do 300 icg2=0 301 do ii=1,isppol 302 ll=nkpt 303 if(isppol==ii) ll=ikpt2-1 304 do jj=1,ll 305 ! MPI: cycle over kpts not treated by this node 306 if (ABS(MPI_enreg%proc_distrb(jj,1,ii)-rank)/=0) CYCLE 307 ! write(std_out,'("kpt loop2: ikpt",i3," rank ",i3)') jj,rank 308 icg2=icg2+nspinor*nband(jj+(ii-1)*nkpt) 309 end do 310 end do 311 ! 312 !MPI: if ikpt2 not found in this processor then 313 !read info from an unformatted file 314 ! 315 if (nprocs>1) then 316 if (ABS(MPI_enreg%proc_distrb(ikpt2,1,isppol)-rank)/=0) then 317 lfile=.true. 318 ! 319 ! get maximum of lmn_size 320 max_lmn=0 321 do itypat=1,ntypat 322 lmn_size=pawtab(itypat)%lmn_size 323 if(lmn_size>max_lmn) max_lmn=lmn_size 324 end do 325 ! 326 ! get file name and open it 327 ! 328 write(cprj_file,'(a,I5.5,".",I1)') trim(seed_name),ikpt2,isppol 329 iunit=1000 330 ! write(std_out,*)'reading file',trim(cprj_file) 331 open (unit=iunit, file=cprj_file,form='unformatted',status='old',iostat=ios) 332 if(ios /= 0) then 333 write(message,*) " smatrix_pawinit: file",trim(cprj_file), "not found" 334 ABI_ERROR(message) 335 end if 336 ! 337 ! start reading 338 do ii=1,mband*nspinor 339 do iatom=1,natom 340 itypat=typat(iatom) 341 lmn_size=pawtab(itypat)%lmn_size 342 do ilmn=1,lmn_size 343 read(iunit)(cprj_k2(iatom,ii)%cp(jj,ilmn),jj=1,2) 344 end do !ilmn 345 end do 346 end do 347 ! 348 ! close file 349 ! 350 close (unit=iunit,iostat=ios) 351 if(ios /= 0) then 352 write(message,*) " smatrix_pawinit: error closing file ",trim(cprj_file) 353 ABI_ERROR(message) 354 end if 355 ! 356 end if 357 end if !mpi 358 359 !Extract cprj_k1 and cprj_k2 360 !these contain the projectors cprj for just one k-point (ikpt1 or ikpt2) 361 362 !Extract cprj for k-point 1 363 iorder_cprj=0 !do not change the ordering of cprj 364 nband_k=nband(ikpt1) 365 dummy=1000 !index of file not implemented here, mkmem==0 not implemented 366 call pawcprj_get(atindx1,cprj_k1,cprj,natom,1,icg1,ikpt1,iorder_cprj,isppol,& 367 & mband,mkmem,natom,nband_k,nband_k,nspinor,nsppol,dummy,& 368 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 369 370 !Extract cprj for k-point 2 371 if( lfile .eqv. .false. ) then !if it was not already read above 372 iorder_cprj=0 !do not change the ordering of cprj 373 nband_k=nband(ikpt2) 374 dummy=1000 !index of file not implemented here, mkmem==0 not implemented 375 call pawcprj_get(atindx1,cprj_k2,cprj,natom,1,icg2,ikpt2,iorder_cprj,isppol,& 376 & mband,mkmem,natom,nband_k,nband_k,nspinor,nsppol,dummy,& 377 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 378 end if 379 380 !DEBUG 381 !if(ikpt2==2) then 382 !if(ikpt1==1) then 383 !do iband1=1,mbandw 384 !do iatom=1,natom 385 !itypat=typat(atindx1(iatom)) 386 !lmn_size=pawtab(itypat)%lmn_size 387 !do ilmn=1,1!lmn_size 388 !!write(500,'(a,i3,a,i3,a,i3,a,2f13.7)')'iband ',iband1,' iatom ',iatom,' ilmn ',ilmn,' cprj',cprj(iatom,iband1+icg1)%cp(:,ilmn) 389 !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) 390 !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) 391 !end do 392 !end do 393 !end do 394 !end if 395 !end if 396 !NDDEBUG 397 398 !!!!!!!!!!!!!!!!! 399 !--- Compute intermediate quantities: "b" vector=k2-k1 and its 400 !normalized value: bbn (and its norm: bnorm) 401 !compute also Ylm(b). 402 ABI_MALLOC(ylmb,(pawang%l_size_max*pawang%l_size_max)) 403 ABI_MALLOC(ylmrgr_dum,(1,1,0)) 404 bb(:)=kpt(:,ikpt2)-kpt(:,ikpt1)+g1(:) 405 bb1=bb 406 xx=gprimd(1,1)*bb(1)+gprimd(1,2)*bb(2)+gprimd(1,3)*bb(3) 407 yy=gprimd(2,1)*bb(1)+gprimd(2,2)*bb(2)+gprimd(2,3)*bb(3) 408 zz=gprimd(3,1)*bb(1)+gprimd(3,2)*bb(2)+gprimd(3,3)*bb(3) 409 bnorm=two_pi*dsqrt(xx**2+yy**2+zz**2) 410 if(bnorm<tol8) then 411 ! write(std_out,*) "WARNING: bnorm=",bnorm 412 bbn(:)=zero 413 else 414 xx=xx*two_pi 415 yy=yy*two_pi 416 zz=zz*two_pi 417 bb(1)=xx 418 bb(2)=yy 419 bb(3)=zz 420 bbn(1)=xx/bnorm 421 bbn(2)=yy/bnorm 422 bbn(3)=zz/bnorm 423 end if 424 425 !debug bbn=0 426 !debug bnorm=0 427 !bbn has to ne normalized 428 call initylmr(pawang%l_size_max,0,1,(/one/),1,bbn(:),ylmb(:),ylmrgr_dum) 429 !write(std_out,*) "ylmb(:)",ylmb(:) 430 !write(std_out,*) pawang%l_size_max 431 !write(std_out,*) "bbn",bbn(:) 432 !write(std_out,*) "xx,yy,zz",xx,yy,zz 433 !write(std_out,*) "bnorm",bnorm 434 ABI_FREE(ylmrgr_dum) 435 436 !------- First Compute Qij(b)- 437 ABI_MALLOC(sb_out, (pawang%l_size_max)) 438 cm2=zero 439 440 do iatom=1,natom 441 itypat=typat(iatom) 442 lmn_size=pawtab(itypat)%lmn_size 443 ! --- en coordonnnes reelles cartesiennes (espace reel) 444 ! --- first radial part(see pawinit) 445 mesh_size=pawtab(itypat)%mesh_size 446 ABI_MALLOC(j_bessel,(mesh_size,pawang%l_size_max)) 447 448 449 ! --- compute bessel function for (br) for all angular momenta necessary 450 ! --- and for all value of r. 451 ! --- they are needed for radial part 452 ! --- of the integration => j_bessel(ir,:) 453 do ir=1,mesh_size 454 arg=bnorm*pawrad(itypat)%rad(ir) 455 call sbf8(pawang%l_size_max,arg,sb_out) 456 j_bessel(ir,:) = sb_out 457 end do 458 459 ! do jlmn=1,pawang%l_size_max 460 ! write(665,*) "j_bessel",j_bessel(1:mesh_size,jlmn) 461 ! enddo 462 ! write(std_out,*) "bessel function computed" 463 ! --- Compute \Sum b.R=xsum for future use 464 xtemp=zero 465 do mm=1,3 466 xtemp=xtemp+xred(mm,iatom)*bb1(mm) 467 end do 468 xtemp=xtemp*two_pi 469 xsum=zero 470 do mm=1,3 471 xsum=xsum+xcart(mm,iatom)*bbn(mm)*bnorm 472 end do 473 ! write(std_out,*)'xsum',xsum,xtemp,lmn_size 474 475 ! --- Loop on jlmn and ilmn 476 qijtot=zero 477 do jlmn=1,lmn_size 478 j0lmn=jlmn*(jlmn-1)/2 479 do ilmn=1,jlmn 480 481 klmn=j0lmn+ilmn 482 klm=pawtab(itypat)%indklmn(1,klmn);kln=pawtab(itypat)%indklmn(2,klmn) 483 lmin=pawtab(itypat)%indklmn(3,klmn);lmax=pawtab(itypat)%indklmn(4,klmn) 484 ! --- Sum over angular momenta 485 ! --- compute radial part integration for each angular momentum => intg 486 ! --- (3j) symbols follows the rule: l belongs to abs(li-lj), li+lj. 487 qijb=zero 488 do ll=lmin,lmax,2 489 lm0=ll*ll+ll+1 490 ABI_MALLOC(ff,(mesh_size)) 491 ff(1:mesh_size)=(pawtab(itypat)%phiphj(1:mesh_size,kln)& 492 & -pawtab(itypat)%tphitphj(1:mesh_size,kln))& 493 & *j_bessel(1:mesh_size,ll+1) 494 call simp_gen(intg,ff,pawrad(itypat)) 495 ABI_FREE(ff) 496 qijbtemp=zero 497 do mm=-ll,ll 498 isel=pawang%gntselect(lm0+mm,klm) 499 if (isel>0) qijbtemp=qijbtemp& 500 & +pawang%realgnt(isel)*ylmb(lm0+mm) 501 end do ! mm 502 ! --- compute angular part with a summation 503 ! --- qijb =\sum_{lm} intg(lm)*qijbtemp 504 qijb(1)=qijb(1) +intg*qijbtemp*ilr(ll+1) 505 qijb(2)=qijb(2) +intg*qijbtemp*ili(ll+1) 506 ! if(ilmn==jlmn) write(std_out,*) "intg, qij",intg,qijbtemp 507 end do ! ll 508 509 ! --- Add exp(-i.b*R) for each atom. 510 if(ilmn==jlmn) qijtot=qijtot+qijb(1) 511 ! if(ilmn==jlmn) write(std_out,*) "qijtot",qijtot 512 x1=qijb(1)*dcos(-xsum)-qijb(2)*dsin(-xsum) 513 x2=qijb(1)*dsin(-xsum)+qijb(2)*dcos(-xsum) 514 ! x1 x2 necessary to avoid changing qijb(1) before 515 ! computing qijb(2) 516 qijb(1)=x1 517 qijb(2)=x2 ! 518 ! if(ilmn==jlmn) write(std_out,*) "qij",jlmn,ilmn,qijb(1),qijb(2) 519 520 do iband1=1,mbandw ! limite inferieure a preciser 521 do iband2=1,mbandw 522 ppr=0.d0 523 ppi=0.d0 524 do ispinor=1,nspinor 525 idx1=iband1*nspinor-(nspinor-ispinor) 526 idx2=iband2*nspinor-(nspinor-ispinor) !to take into account spinors 527 ! write(std_out,*) "iband2",iband2 528 ! product of (a1+ia2)*(b1-ib2) (minus sign because conjugated) 529 ppr=ppr+& 530 ! real part a_1*b_1+a_2*b_2 531 & cprj_k1(iatom,idx1)%cp(1,ilmn)*cprj_k2(iatom,idx2)%cp(1,jlmn)+& 532 & cprj_k1(iatom,idx1)%cp(2,ilmn)*cprj_k2(iatom,idx2)%cp(2,jlmn)+& 533 ! & cprj(iatom,idx1+icg1)%cp(1,ilmn)*cprj(iatom,idx2+icg2)%cp(1,jlmn)+& 534 ! & cprj(iatom,idx1+icg1)%cp(2,ilmn)*cprj(iatom,idx2+icg2)%cp(2,jlmn)+& 535 ! add term on the other triangle of the matrix 536 ! qij is the same for this part because phi are real. 537 & cprj_k1(iatom,idx1)%cp(1,jlmn)*cprj_k2(iatom,idx2)%cp(1,ilmn)+& 538 & cprj_k1(iatom,idx1)%cp(2,jlmn)*cprj_k2(iatom,idx2)%cp(2,ilmn) 539 ! & cprj(iatom,idx1+icg1)%cp(1,jlmn)*cprj(iatom,idx2+icg2)%cp(1,ilmn)+& 540 ! & cprj(iatom,idx1+icg1)%cp(2,jlmn)*cprj(iatom,idx2+icg2)%cp(2,ilmn) 541 ppi=ppi+& 542 ! imaginary part a_1*b_2-a_2*b_1 543 & cprj_k1(iatom,idx1)%cp(1,ilmn)*cprj_k2(iatom,idx2)%cp(2,jlmn)-& 544 & cprj_k1(iatom,idx1)%cp(2,ilmn)*cprj_k2(iatom,idx2)%cp(1,jlmn)+& 545 ! & cprj(iatom,idx1+icg1)%cp(1,ilmn)*cprj(iatom,idx2+icg2)%cp(2,jlmn)-& 546 ! & cprj(iatom,idx1+icg1)%cp(2,ilmn)*cprj(iatom,idx2+icg2)%cp(1,jlmn)+& 547 ! add term on the other triangle of the matrix 548 & cprj_k1(iatom,idx1)%cp(1,jlmn)*cprj_k2(iatom,idx2)%cp(2,ilmn)-& 549 & cprj_k1(iatom,idx1)%cp(2,jlmn)*cprj_k2(iatom,idx2)%cp(1,ilmn) 550 ! & cprj(iatom,idx1+icg1)%cp(1,jlmn)*cprj(iatom,idx2+icg2)%cp(2,ilmn)-& 551 ! & cprj(iatom,idx1+icg1)%cp(2,jlmn)*cprj(iatom,idx2+icg2)%cp(1,ilmn) 552 end do !ispinor 553 ! 554 ! delta: diagonal terms are counted twice ! so 555 ! we need a 0.5 factor for diagonal elements. 556 delta=one 557 ! write(std_out,*) "ppr and ppi computed",ikpt1,ikpt2,iband1,iband2 558 if(ilmn==jlmn) delta=half 559 cm2(1,iband1,iband2)= cm2(1,iband1,iband2)+ & 560 & (qijb(1)*ppr-qijb(2)*ppi)*delta 561 cm2(2,iband1,iband2)= cm2(2,iband1,iband2)+ & 562 & (qijb(2)*ppr+qijb(1)*ppi)*delta 563 end do ! iband2 564 end do ! iband1 565 566 end do ! ilmn 567 end do ! jlmn 568 ! write(std_out,*) "final qijtot",qijtot 569 ABI_FREE(j_bessel) 570 end do ! iatom 571 572 ABI_FREE(sb_out) 573 ABI_FREE(ylmb) 574 call pawcprj_free(cprj_k1) 575 call pawcprj_free(cprj_k2) 576 ABI_FREE(cprj_k1) 577 ABI_FREE(cprj_k2) 578 579 DBG_EXIT("COLL") 580 581 end subroutine smatrix_pawinit