TABLE OF CONTENTS


ABINIT/chern_number [ Functions ]

[ Top ] [ Functions ]

NAME

 chern_number

FUNCTION

 This routine computes the Chern number based on input wavefunctions.
 It is assumed that only completely filled bands are present.

COPYRIGHT

 Copyright (C) 2003-2017 ABINIT  group (JWZ)
 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

 atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
 cg(2,mcg)=planewave coefficients of wavefunctions
 cprj(natom,mcprj*usecrpj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
 dtset <type(dataset_type)>=all input variables in this dataset
 gmet(3,3)=metric in reciprocal space
 gprimd(3,3)=reciprocal space dimensional primitive translations
 kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere
 mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
 mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
 mpi_enreg=information about MPI parallelization
 npwarr(nkpt)=number of planewaves in basis at this k point
 pawang <type(pawang_type)>=paw angular mesh and related data
 pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
 pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
 pwind(pwind_alloc,2,3) = array used to compute
           the overlap matrix smat between k-points (see initberry.f)
 pwind_alloc = first dimension of pwind
 symrec(3,3,nsym) = symmetries in reciprocal space in terms of
   reciprocal space primitive translations
 usecprj=1 if cprj datastructure has been allocated
 usepaw=1 if PAW calculation
 xred(3,natom) = location of atoms in unit cell

OUTPUT

SIDE EFFECTS

 dtorbmag <type(orbmag_type)> = variables related to orbital magnetization

TODO

NOTES

 See Ceresoli et al, PRB 74, 024408 (2006) [[cite:Ceresoli2006]],
 and Gonze and Zwanziger, PRB 84 064445 (2011) [[cite:Gonze2011a]].
 This routine computes the Chern number as
 $C_\alpha = \frac{i}{2\pi}\int_{\mathrm{BZ}} dk \epsilon_{\alpha\beta\gamma}
 \mathrm{Tr}[\rho_k \partial_\beta \rho_k (1 - \rho_k) \partial_gamma\rho_k] $
 The derivative of the density operator is obtained from a discretized formula
 $\partial_\beta \rho_k = \frac{1}{2\Delta}(\rho_{k+b} - \rho_{k-b})$ with
 $\Delta = |b|$. When reduced to wavefunction overlaps the computation amounts to
 multiple calls to smatrix.F90, exactly as in other Berry phase computations, with
 the one additional complication of overlaps like $\langle u_{n,k+b}|u_{n',k+g}\rangle$.
 At this stage mkpwind_k is invoked, which generalizes the code in initberry
 and initorbmag necessary to index plane waves around different k points.
 Direct questions and comments to J Zwanziger

PARENTS

CHILDREN

SOURCE

 938 subroutine chern_number(atindx1,cg,cprj,dtset,dtorbmag,gmet,gprimd,kg,&
 939      &            mcg,mcprj,mpi_enreg,npwarr,pawang,pawrad,pawtab,pwind,pwind_alloc,&
 940      &            symrec,usecprj,usepaw,xred)
 941 
 942 
 943 !This section has been created automatically by the script Abilint (TD).
 944 !Do not modify the following lines by hand.
 945 #undef ABI_FUNC
 946 #define ABI_FUNC 'chern_number'
 947 !End of the abilint section
 948 
 949   implicit none
 950 
 951   !Arguments ------------------------------------
 952   !scalars
 953   integer,intent(in) :: mcg,mcprj,pwind_alloc,usecprj,usepaw
 954   type(dataset_type),intent(in) :: dtset
 955   type(MPI_type), intent(inout) :: mpi_enreg
 956   type(orbmag_type), intent(inout) :: dtorbmag
 957   type(pawang_type),intent(in) :: pawang
 958 
 959   !arrays
 960   integer,intent(in) :: atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem)
 961   integer,intent(in) :: npwarr(dtset%nkpt),pwind(pwind_alloc,2,3),symrec(3,3,dtset%nsym)
 962   real(dp), intent(in) :: cg(2,mcg),gmet(3,3),gprimd(3,3),xred(3,dtset%natom)
 963   type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*usepaw)
 964   type(pawcprj_type),intent(in) ::  cprj(dtset%natom,mcprj*usecprj)
 965   type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*usepaw)
 966 
 967   !Local variables -------------------------
 968   !scalars
 969   integer :: adir,bdir,bdx,bdxc,bfor,bsigma,ddkflag,epsabg,gdir,gdx,gdxc,gfor,gsigma
 970   integer :: icg,icgb,icgg,icprj,icprjb,icprjg
 971   integer :: ikg,ikpt,ikptb,ikptg,isppol,itrs,job
 972   integer :: mcg1_k,my_nspinor,nband_k,ncpgr,nn,n1,n2,n3,npw_k,npw_kb,npw_kg,shiftbd
 973   real(dp) :: deltab,deltag
 974   complex(dpc) :: IA,IB,t1A,t2A,t3A,t1B,t2B,t3B,t4B,tprodA,tprodB
 975   character(len=500) :: message
 976   !arrays
 977   integer,allocatable :: dimlmn(:),nattyp_dum(:),pwind_kb(:),pwind_kg(:),pwind_bg(:),sflag_k(:)
 978   real(dp) :: cnum(2,3),dkb(3),dkg(3),dkbg(3),dtm_k(2)
 979   real(dp),allocatable :: cg1_k(:,:),kk_paw(:,:,:),pwnsfac_k(:,:)
 980   real(dp),allocatable :: smat_all_indx(:,:,:,:,:,:),smat_inv(:,:,:),smat_kk(:,:,:)
 981   logical,allocatable :: has_smat_indx(:,:,:)
 982   type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kb(:,:),cprj_kg(:,:)
 983   type(pawcprj_type),allocatable :: cprj_fkn(:,:),cprj_ikn(:,:)
 984 
 985   ! ***********************************************************************
 986   ! my_nspinor=max(1,dtorbmag%nspinor/mpi_enreg%nproc_spinor)
 987 
 988   ! TODO: generalize to nsppol > 1
 989   isppol = 1
 990   my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
 991 
 992   nband_k = dtorbmag%mband_occ
 993 
 994   if (usepaw == 1) then ! cprj allocation
 995      ncpgr = cprj(1,1)%ncpgr
 996      ABI_ALLOCATE(dimlmn,(dtset%natom))
 997      call pawcprj_getdim(dimlmn,dtset%natom,nattyp_dum,dtset%ntypat,dtset%typat,pawtab,'R')
 998      ABI_DATATYPE_ALLOCATE(cprj_k,(dtset%natom,dtorbmag%nspinor*dtset%mband))
 999      call pawcprj_alloc(cprj_k,ncpgr,dimlmn)
1000      ABI_DATATYPE_ALLOCATE(cprj_kb,(dtset%natom,dtorbmag%nspinor*dtset%mband))
1001      call pawcprj_alloc(cprj_kb,ncpgr,dimlmn)
1002      ABI_DATATYPE_ALLOCATE(cprj_kg,(dtset%natom,dtorbmag%nspinor*dtset%mband))
1003      call pawcprj_alloc(cprj_kg,ncpgr,dimlmn)
1004      if (dtset%kptopt /= 3) then
1005         ABI_DATATYPE_ALLOCATE(cprj_ikn,(dtset%natom,dtorbmag%nspinor*dtset%mband))
1006         ABI_DATATYPE_ALLOCATE(cprj_fkn,(dtset%natom,dtorbmag%nspinor*dtset%mband))
1007         call pawcprj_alloc(cprj_ikn,ncpgr,dimlmn)
1008         call pawcprj_alloc(cprj_fkn,ncpgr,dimlmn)
1009      end if
1010   else
1011      message = ' usepaw /= 1 but Chern number calculation requires PAW '
1012      MSG_ERROR(message)
1013   end if
1014 
1015   ABI_ALLOCATE(kk_paw,(2,dtset%mband,dtset%mband))
1016   ABI_ALLOCATE(pwind_kb,(dtset%mpw))
1017   ABI_ALLOCATE(pwind_kg,(dtset%mpw))
1018   ABI_ALLOCATE(pwind_bg,(dtset%mpw))
1019   ABI_ALLOCATE(pwnsfac_k,(4,dtset%mpw))
1020   pwnsfac_k(1,:) = 1.0_dp ! bra real
1021   pwnsfac_k(2,:) = 0.0_dp ! bra imag
1022   pwnsfac_k(3,:) = 1.0_dp ! ket real
1023   pwnsfac_k(4,:) = 0.0_dp ! ket imag
1024 
1025   mcg1_k = dtset%mpw*nband_k
1026   ABI_ALLOCATE(cg1_k,(2,mcg1_k))
1027   ABI_ALLOCATE(sflag_k,(nband_k))
1028   ABI_ALLOCATE(smat_inv,(2,nband_k,nband_k))
1029   ABI_ALLOCATE(smat_kk,(2,nband_k,nband_k))
1030 
1031   ddkflag = 1
1032 
1033   !itrs = 0 means do not invoke time reversal symmetry in smatrix.F90
1034   itrs = 0
1035 
1036   job = 1
1037   shiftbd = 1
1038 
1039   ABI_ALLOCATE(has_smat_indx,(dtorbmag%fnkpt,0:6,0:6))
1040   ABI_ALLOCATE(smat_all_indx,(2,nband_k,nband_k,dtorbmag%fnkpt,0:6,0:6))
1041   has_smat_indx(:,:,:)=.FALSE.
1042 
1043   do adir = 1, 3
1044 
1045      IA = czero
1046      IB = czero
1047 
1048      do epsabg = 1, -1, -2
1049 
1050         if (epsabg .EQ. 1) then
1051            bdir = modulo(adir,3)+1
1052            gdir = modulo(adir+1,3)+1
1053         else
1054            bdir = modulo(adir+1,3)+1
1055            gdir = modulo(adir,3)+1
1056         end if
1057 
1058         ! loop over kpts, assuming for now kptopt 3, nsppol = 1, nspinor = 1
1059         ! and no parallelism, no symmorphic symmetry elements
1060 
1061 
1062         do ikpt = 1, dtorbmag%fnkpt
1063 
1064            icprj = dtorbmag%cprjindex(ikpt,isppol)
1065 
1066            npw_k = npwarr(ikpt)
1067            icg = dtorbmag%cgindex(ikpt,dtset%nsppol)
1068 
1069            ikg = dtorbmag%fkgindex(ikpt)
1070 
1071            call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,icprj,ikpt,0,isppol,dtset%mband,&
1072                 &       dtset%mkmem,dtset%natom,nband_k,nband_k,my_nspinor,dtset%nsppol,0)
1073 
1074            do bfor = 1, 2
1075               if (bfor .EQ. 1) then
1076                  bsigma = 1
1077               else
1078                  bsigma = -1
1079               end if
1080               ! index of neighbor 1..6
1081               bdx = 2*bdir-2+bfor
1082               ! index of ikpt viewed from neighbor
1083               bdxc = 2*bdir-2+bfor+bsigma
1084 
1085               dkb(1:3) = bsigma*dtorbmag%dkvecs(1:3,bdir)
1086               deltab = sqrt(DOT_PRODUCT(dkb,MATMUL(gmet,dkb)))
1087 
1088               ikptb = dtorbmag%ikpt_dk(ikpt,bfor,bdir)
1089               icprjb = dtorbmag%cprjindex(ikptb,isppol)
1090 
1091               npw_kb = npwarr(ikptb)
1092               icgb = dtorbmag%cgindex(ikptb,dtset%nsppol)
1093 
1094               pwind_kb(1:npw_k) = pwind(ikg+1:ikg+npw_k,bfor,bdir)
1095 
1096               call pawcprj_get(atindx1,cprj_kb,cprj,dtset%natom,1,icprjb,&
1097                    &         ikptb,0,isppol,dtset%mband,dtset%mkmem,dtset%natom,nband_k,nband_k,&
1098                    &         my_nspinor,dtset%nsppol,0)
1099 
1100               if (.NOT. has_smat_indx(ikpt,bdx,0)) then
1101 
1102                  call overlap_k1k2_paw(cprj_k,cprj_kb,dkb,gprimd,kk_paw,dtorbmag%lmn2max,&
1103                       &           dtorbmag%lmn_size,dtset%mband,&
1104                       &           dtset%natom,my_nspinor,dtset%ntypat,pawang,pawrad,pawtab,dtset%typat,xred)
1105 
1106                  sflag_k=0
1107                  call smatrix(cg,cg,cg1_k,ddkflag,dtm_k,icg,icgb,itrs,job,nband_k,&
1108                       &           mcg,mcg,mcg1_k,1,dtset%mpw,nband_k,nband_k,npw_k,npw_kb,my_nspinor,&
1109                       &           pwind_kb,pwnsfac_k,sflag_k,shiftbd,smat_inv,smat_kk,kk_paw,usepaw)
1110 
1111                  do nn = 1, nband_k
1112                     do n1 = 1, nband_k
1113                        smat_all_indx(1,nn,n1,ikpt,bdx,0) =  smat_kk(1,nn,n1)
1114                        smat_all_indx(2,nn,n1,ikpt,bdx,0) =  smat_kk(2,nn,n1)
1115                        smat_all_indx(1,n1,nn,ikptb,bdxc,0) =  smat_kk(1,nn,n1)
1116                        smat_all_indx(2,n1,nn,ikptb,bdxc,0) = -smat_kk(2,nn,n1)
1117                     end do
1118                  end do
1119 
1120                  has_smat_indx(ikpt,bdx,0) = .TRUE.
1121                  has_smat_indx(ikptb,bdxc,0) = .TRUE.
1122 
1123               end if
1124 
1125               do gfor = 1, 2
1126                  if (gfor .EQ. 1) then
1127                     gsigma = 1
1128                  else
1129                     gsigma = -1
1130                  end if
1131                  ! index of neighbor 1..6
1132                  gdx = 2*gdir-2+gfor
1133                  ! index of ikpt viewed from neighbor
1134                  gdxc = 2*gdir-2+gfor+gsigma
1135 
1136                  dkg(1:3) = gsigma*dtorbmag%dkvecs(1:3,gdir)
1137                  deltag = sqrt(DOT_PRODUCT(dkg,MATMUL(gmet,dkg)))
1138 
1139                  ikptg = dtorbmag%ikpt_dk(ikpt,gfor,gdir)
1140                  icprjg = dtorbmag%cprjindex(ikptg,isppol)
1141 
1142                  npw_kg = npwarr(ikptg)
1143                  icgg = dtorbmag%cgindex(ikptg,dtset%nsppol)
1144 
1145                  pwind_kg(1:npw_k) = pwind(ikg+1:ikg+npw_k,gfor,gdir)
1146 
1147                  call pawcprj_get(atindx1,cprj_kg,cprj,dtset%natom,1,icprjg,&
1148                       &           ikptg,0,isppol,dtset%mband,dtset%mkmem,dtset%natom,nband_k,nband_k,&
1149                       &           my_nspinor,dtset%nsppol,0)
1150 
1151                  if (.NOT. has_smat_indx(ikpt,gdx,0)) then
1152 
1153                     call overlap_k1k2_paw(cprj_k,cprj_kg,dkg,gprimd,kk_paw,dtorbmag%lmn2max,&
1154                          &             dtorbmag%lmn_size,dtset%mband,&
1155                          &             dtset%natom,my_nspinor,dtset%ntypat,pawang,pawrad,pawtab,dtset%typat,xred)
1156 
1157                     sflag_k=0
1158                     call smatrix(cg,cg,cg1_k,ddkflag,dtm_k,icg,icgg,itrs,job,nband_k,&
1159                          &             mcg,mcg,mcg1_k,1,dtset%mpw,nband_k,nband_k,npw_k,npw_kg,my_nspinor,&
1160                          &             pwind_kg,pwnsfac_k,sflag_k,shiftbd,smat_inv,smat_kk,kk_paw,usepaw)
1161 
1162                     do nn = 1, nband_k
1163                        do n1 = 1, nband_k
1164                           smat_all_indx(1,nn,n1,ikpt,gdx,0) =  smat_kk(1,nn,n1)
1165                           smat_all_indx(2,nn,n1,ikpt,gdx,0) =  smat_kk(2,nn,n1)
1166                           smat_all_indx(1,n1,nn,ikptg,gdxc,0) =  smat_kk(1,nn,n1)
1167                           smat_all_indx(2,n1,nn,ikptg,gdxc,0) = -smat_kk(2,nn,n1)
1168                        end do
1169                     end do
1170 
1171                     has_smat_indx(ikpt,gdx,0) = .TRUE.
1172                     has_smat_indx(ikptg,gdxc,0) = .TRUE.
1173 
1174                  end if
1175 
1176                  dkbg = dkg - dkb
1177 
1178                  if (.NOT. has_smat_indx(ikpt,bdx,gdx)) then
1179 
1180                     call overlap_k1k2_paw(cprj_kb,cprj_kg,dkbg,gprimd,kk_paw,dtorbmag%lmn2max,&
1181                          &             dtorbmag%lmn_size,dtset%mband,&
1182                          &             dtset%natom,my_nspinor,dtset%ntypat,pawang,pawrad,pawtab,dtset%typat,xred)
1183 
1184                     call mkpwind_k(dkbg,dtset,dtorbmag%fnkpt,dtorbmag%fkptns,gmet,&
1185                          &             dtorbmag%indkk_f2ibz,ikptb,ikptg,&
1186                          &             kg,dtorbmag%kgindex,mpi_enreg,npw_kb,pwind_bg,symrec)
1187 
1188                     sflag_k=0
1189                     call smatrix(cg,cg,cg1_k,ddkflag,dtm_k,icgb,icgg,itrs,job,nband_k,&
1190                          &             mcg,mcg,mcg1_k,1,dtset%mpw,nband_k,nband_k,npw_kb,npw_kg,my_nspinor,&
1191                          &             pwind_bg,pwnsfac_k,sflag_k,shiftbd,smat_inv,smat_kk,kk_paw,usepaw)
1192 
1193                     do nn = 1, nband_k
1194                        do n1 = 1, nband_k
1195                           smat_all_indx(1,nn,n1,ikpt,bdx,gdx) =  smat_kk(1,nn,n1)
1196                           smat_all_indx(2,nn,n1,ikpt,bdx,gdx) =  smat_kk(2,nn,n1)
1197                           smat_all_indx(1,n1,nn,ikpt,gdx,bdx) =  smat_kk(1,nn,n1)
1198                           smat_all_indx(2,n1,nn,ikpt,gdx,bdx) = -smat_kk(2,nn,n1)
1199                        end do
1200                     end do
1201 
1202                     has_smat_indx(ikpt,bdx,gdx) = .TRUE.
1203                     has_smat_indx(ikpt,gdx,bdx) = .TRUE.
1204 
1205                  end if
1206 
1207                  do nn = 1, nband_k
1208                     do n1 = 1, nband_k
1209 
1210                        t1A = cmplx(smat_all_indx(1,nn,n1,ikpt,bdx,0),smat_all_indx(2,nn,n1,ikpt,bdx,0))
1211                        t1B = t1A
1212 
1213                        do n2 = 1, nband_k
1214 
1215                           t2A = cmplx(smat_all_indx(1,n1,n2,ikpt,bdx,gdx),smat_all_indx(2,n1,n2,ikpt,bdx,gdx))
1216                           t3A = conjg(cmplx(smat_all_indx(1,nn,n2,ikpt,gdx,0),smat_all_indx(2,nn,n2,ikpt,gdx,0)))
1217 
1218                           t2B = conjg(cmplx(smat_all_indx(1,n2,n1,ikpt,bdx,0),smat_all_indx(2,n2,n1,ikpt,bdx,0)))
1219 
1220                           do n3 = 1, nband_k
1221 
1222                              t3B = cmplx(smat_all_indx(1,n2,n3,ikpt,gdx,0),smat_all_indx(2,n2,n3,ikpt,gdx,0))
1223                              t4B=conjg(cmplx(smat_all_indx(1,nn,n3,ikpt,gdx,0),smat_all_indx(2,nn,n3,ikpt,gdx,0)))
1224 
1225                              tprodB = t1B*t2B*t3B*t4B
1226                              IB = IB - epsabg*bsigma*gsigma*tprodB/(2.0*deltab*2.0*deltag)
1227                           end do ! end loop over n3
1228 
1229                           tprodA = t1A*t2A*t3A
1230                           IA = IA + epsabg*bsigma*gsigma*tprodA/(2.0*deltab*2.0*deltag)
1231 
1232 
1233                        end do ! end loop over n2
1234                     end do ! end loop over n1
1235                  end do ! end loop over nn
1236 
1237               end do ! end loop over gfor
1238 
1239            end do ! end loop over bfor
1240 
1241         end do ! end loop over fnkpt
1242 
1243      end do ! end loop over epsabg
1244 
1245      cnum(1,adir) = real(IA+IB)
1246      cnum(2,adir) = aimag(IA+IB)
1247 
1248   end do ! end loop over adir
1249 
1250   cnum(1,1:3) = MATMUL(gprimd,cnum(1,1:3))
1251   cnum(2,1:3) = MATMUL(gprimd,cnum(2,1:3))
1252   dtorbmag%chern(1,1:3) = -cnum(2,1:3)/(two_pi*dtorbmag%fnkpt)
1253   dtorbmag%chern(2,1:3) =  cnum(1,1:3)/(two_pi*dtorbmag%fnkpt)
1254 
1255   write(message,'(a,a,a)')ch10,'====================================================',ch10
1256   call wrtout(ab_out,message,'COLL')
1257 
1258   write(message,'(a)')' Chern number C from orbital magnetization '
1259   call wrtout(ab_out,message,'COLL')
1260   write(message,'(a,a)')'----C is a real vector, given along Cartesian directions----',ch10
1261   call wrtout(ab_out,message,'COLL')
1262 
1263   do adir = 1, 3
1264      write(message,'(a,i4,a,2es16.8)')' C(',adir,') : real, imag ',&
1265           &   dtorbmag%chern(1,adir),dtorbmag%chern(2,adir)
1266      call wrtout(ab_out,message,'COLL')
1267   end do
1268 
1269   write(message,'(a,a,a)')ch10,'====================================================',ch10
1270   call wrtout(ab_out,message,'COLL')
1271 
1272   if (usepaw == 1) then
1273      ABI_DEALLOCATE(dimlmn)
1274      call pawcprj_free(cprj_k)
1275      ABI_DATATYPE_DEALLOCATE(cprj_k)
1276      call pawcprj_free(cprj_kb)
1277      ABI_DATATYPE_DEALLOCATE(cprj_kb)
1278      call pawcprj_free(cprj_kg)
1279      ABI_DATATYPE_DEALLOCATE(cprj_kg)
1280      if (dtset%kptopt /= 3) then
1281         call pawcprj_free(cprj_ikn)
1282         call pawcprj_free(cprj_fkn)
1283         ABI_DATATYPE_DEALLOCATE(cprj_ikn)
1284         ABI_DATATYPE_DEALLOCATE(cprj_fkn)
1285      end if
1286   end if
1287 
1288   ABI_DEALLOCATE(kk_paw)
1289   ABI_DEALLOCATE(cg1_k)
1290   ABI_DEALLOCATE(sflag_k)
1291   ABI_DEALLOCATE(smat_inv)
1292   ABI_DEALLOCATE(smat_kk)
1293   ABI_DEALLOCATE(pwind_kb)
1294   ABI_DEALLOCATE(pwind_kg)
1295   ABI_DEALLOCATE(pwind_bg)
1296   ABI_DEALLOCATE(pwnsfac_k)
1297 
1298   ABI_DEALLOCATE(has_smat_indx)
1299   ABI_DEALLOCATE(smat_all_indx)
1300 
1301 end subroutine chern_number

ABINIT/ctocprjb [ Functions ]

[ Top ] [ Functions ]

NAME

 ctocprjb

FUNCTION

 Compute <p_k+b|u_k> cprj's as needed by orbital magnetization,
 at all k points and all bands

COPYRIGHT

 Copyright (C) 2003-2017 ABINIT  group
 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

OUTPUT

SIDE EFFECTS

TODO

NOTES

PARENTS

CHILDREN

SOURCE

2127 subroutine ctocprjb(atindx1,cg,cprj_kb_k,dtorbmag,dtset,gmet,gprimd,&
2128      & istwf_k,kg,mcg,mpi_enreg,nattyp,ncpgr,npwarr,pawtab,psps,rmet,rprimd,ucvol,xred)
2129 
2130 
2131 !This section has been created automatically by the script Abilint (TD).
2132 !Do not modify the following lines by hand.
2133 #undef ABI_FUNC
2134 #define ABI_FUNC 'ctocprjb'
2135 !End of the abilint section
2136 
2137   implicit none
2138 
2139   !Arguments ------------------------------------
2140   !scalars
2141   integer,intent(in) :: istwf_k,mcg,ncpgr
2142   real(dp),intent(in) :: ucvol
2143   type(dataset_type),intent(in) :: dtset
2144   type(MPI_type), intent(inout) :: mpi_enreg
2145   type(orbmag_type), intent(inout) :: dtorbmag
2146   type(pseudopotential_type),intent(in) :: psps
2147 
2148   !arrays
2149   integer,intent(in) :: atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem)
2150   integer,intent(in) :: nattyp(dtset%ntypat),npwarr(dtset%nkpt)
2151   real(dp),intent(in) :: cg(2,mcg),gmet(3,3),gprimd(3,3),rmet(3,3),rprimd(3,3),xred(3,dtset%natom)
2152   type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
2153   type(pawcprj_type),intent(inout) :: cprj_kb_k(dtorbmag%fnkpt,6,dtset%natom,dtorbmag%nspinor*dtset%mband)
2154 
2155   !Locals------------------------------------
2156   !scalars
2157   integer :: bdir,bdx,bfor,bsigma,choice,cpopt,dimffnl,dimph1d,ia,iband,icg,ider,idir,ikg,ikpt
2158   integer :: n1,n2,n3,nband_k,nkpg,npw_k,optder
2159   real(dp) :: arg
2160 
2161   !arrays
2162   integer :: nband_dum(1),npwarr_dum(1)
2163   integer,allocatable :: dimlmn(:),kg_k(:,:)
2164   real(dp) :: dkb(3),kpoint(3),kpointb(3),kptns(3,1)
2165   real(dp),allocatable :: cwavef(:,:),ffnl(:,:,:,:),kpg_k_dummy(:,:)
2166   real(dp),allocatable :: ph1d(:,:),ph3d(:,:,:),phkxred(:,:)
2167   real(dp),allocatable :: ylm_k(:,:),ylm_k_gr(:,:,:)
2168   type(pawcprj_type),allocatable :: cwaveprj(:,:)
2169 
2170   ! ***********************************************************************
2171 
2172   nband_k = dtorbmag%mband_occ
2173 
2174   ABI_ALLOCATE(dimlmn,(dtset%natom))
2175   call pawcprj_getdim(dimlmn,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'R')
2176 
2177   ABI_DATATYPE_ALLOCATE(cwaveprj,(dtset%natom,1))
2178   call pawcprj_alloc(cwaveprj,ncpgr,dimlmn)
2179 
2180   ABI_ALLOCATE(phkxred,(2,dtset%natom))
2181 
2182   n1=dtset%ngfft(1); n2=dtset%ngfft(2); n3=dtset%ngfft(3)
2183   dimph1d=dtset%natom*(2*(n1+n2+n3)+3)
2184   ABI_ALLOCATE(ph1d,(2,dimph1d))
2185   call getph(atindx1,dtset%natom,n1,n2,n3,ph1d,xred)
2186 
2187   do ikpt=1,dtorbmag%fnkpt
2188 
2189      kpoint(:)=dtorbmag%fkptns(:,ikpt)
2190      npw_k = npwarr(ikpt)
2191      ABI_ALLOCATE(cwavef,(2,npw_k))
2192 
2193      dimffnl=1
2194      ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,dtset%ntypat))
2195 
2196      icg = dtorbmag%cgindex(ikpt,dtset%nsppol)
2197 
2198      ikg = dtorbmag%fkgindex(ikpt)
2199      ABI_ALLOCATE(kg_k,(3,npw_k))
2200      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
2201      nkpg = 0
2202      ABI_ALLOCATE(kpg_k_dummy,(npw_k,nkpg))
2203 
2204      ABI_ALLOCATE(ph3d,(2,npw_k,dtset%natom))
2205 
2206      ! data for initylmg call below
2207      optder=0
2208      nband_dum(1) = nband_k
2209      npwarr_dum(1) = npw_k
2210      ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang))
2211      ABI_ALLOCATE(ylm_k_gr,(npw_k,3+6*(optder/2),psps%mpsang*psps%mpsang))
2212 
2213      do bdir=1, 3
2214         do bfor=1, 2
2215 
2216            if (bfor .EQ. 1) then
2217               bsigma = 1
2218            else
2219               bsigma = -1
2220            end if
2221 
2222            bdx = 2*bdir-2+bfor
2223 
2224            dkb(1:3) = bsigma*dtorbmag%dkvecs(1:3,bdir)
2225 
2226            kpointb(:) = kpoint(:) + dkb(:)
2227 
2228            do ia=1, dtset%natom
2229               arg=two_pi*(kpointb(1)*xred(1,ia)+kpointb(2)*xred(2,ia)+kpointb(3)*xred(3,ia))
2230               phkxred(1,ia)=cos(arg);phkxred(2,ia)=sin(arg)
2231            end do
2232 
2233            call ph1d3d(1,dtset%natom,kg_k,dtset%natom,dtset%natom,npw_k,n1,n2,n3,phkxred,ph1d,ph3d)
2234 
2235            kptns(:,1) = kpointb(:)
2236            call initylmg(gprimd,kg_k,kptns,1,mpi_enreg,psps%mpsang,npw_k,&
2237                 & nband_dum,1,npwarr_dum,dtset%nsppol,optder,rprimd,ylm_k,ylm_k_gr)
2238 
2239            !      Compute nonlocal form factors ffnl at all (k+b+G_k):
2240            ider=0;idir=0
2241            call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
2242                 & gmet,gprimd,ider,idir,psps%indlmn,kg_k,kpg_k_dummy,kpointb,psps%lmnmax,&
2243                 & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
2244                 & npw_k,dtset%ntypat,psps%pspso,psps%qgrid_ff,rmet,&
2245                 & psps%usepaw,psps%useylm,ylm_k,ylm_k_gr)
2246 
2247            choice=1;cpopt=0;idir=0
2248            do iband = 1, nband_k
2249               cwavef(1,1:npw_k) = cg(1,icg+(iband-1)*npw_k+1:icg+iband*npw_k)
2250               cwavef(2,1:npw_k) = cg(2,icg+(iband-1)*npw_k+1:icg+iband*npw_k)
2251 
2252               call getcprj(choice,cpopt,cwavef,cwaveprj,ffnl,&
2253                    & idir,psps%indlmn,istwf_k,kg_k,kpg_k_dummy,kpointb,psps%lmnmax,&
2254                    & dtset%mgfft,mpi_enreg,&
2255                    & dtset%natom,nattyp,dtset%ngfft,dtset%nloalg,npw_k,dtset%nspinor,dtset%ntypat,&
2256                    & phkxred,ph1d,ph3d,ucvol,psps%useylm)
2257 
2258               call pawcprj_put(atindx1,cwaveprj,cprj_kb_k(ikpt,bdx,:,:),dtset%natom,&
2259                    & iband,0,ikpt,0,1,nband_k,1,dtset%natom,1,nband_k,dimlmn,dtset%nspinor,dtset%nsppol,0)
2260 
2261            end do ! end loop over bands
2262 
2263         end do ! end loop over bfor
2264      end do ! end loop over bdir
2265 
2266      ABI_DEALLOCATE(kg_k)
2267      ABI_DEALLOCATE(ph3d)
2268      ABI_DEALLOCATE(kpg_k_dummy)
2269      ABI_DEALLOCATE(cwavef)
2270      ABI_DEALLOCATE(ffnl)
2271      ABI_DEALLOCATE(ylm_k)
2272      ABI_DEALLOCATE(ylm_k_gr)
2273 
2274   end do ! end loop over nkpt
2275 
2276   ABI_DEALLOCATE(dimlmn)
2277   call pawcprj_free(cwaveprj)
2278   ABI_DATATYPE_DEALLOCATE(cwaveprj)
2279 
2280   ABI_DEALLOCATE(phkxred)
2281   ABI_DEALLOCATE(ph1d)
2282 
2283 end subroutine ctocprjb

ABINIT/initorbmag [ Functions ]

[ Top ] [ Functions ]

NAME

 initorbmag

FUNCTION

 Initialization of orbital magnetization calculation; similar to initberry

COPYRIGHT

 Copyright (C) 2004-2017 ABINIT group.
 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

  dtset <type(dataset_type)> = all input variables in this dataset
  gmet(3,3) = reciprocal space metric tensor in bohr**-2
  gprimd(3,3) = primitive translations in recip space
  kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere
  npwarr(nkpt) = number of planewaves in basis and boundary at this k point
  occ(mband*nkpt*nsppol) = occup number for each band at each k point
  pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
  psps <type(pseudopotential_type)>=variables related to pseudopotentials
  rprimd(3,3) = dimensional primitive vectors
  symrec(3,3,nsym) = symmetries in reciprocal space in terms of
    reciprocal space primitive translations
  xred(3,natom) = location of atoms in reduced units

OUTPUT

  dtorbmag <type(orbmag_type)> = variables related to orbital magnetization

SIDE EFFECTS

  mpi_enreg = information about MPI parallelization

PARENTS

      gstate

CHILDREN

      kpgsph,listkk,setsym_ylm,smpbz,symatm,timab,wrtout,xmpi_max,xmpi_sum

SOURCE

295 subroutine initorbmag(dtorbmag,dtset,gmet,gprimd,kg,mpi_enreg,npwarr,occ,&
296      &                     pawtab,psps,pwind,pwind_alloc,pwnsfac,&
297      &                     rprimd,symrec,xred)
298 
299 
300 !This section has been created automatically by the script Abilint (TD).
301 !Do not modify the following lines by hand.
302 #undef ABI_FUNC
303 #define ABI_FUNC 'initorbmag'
304 !End of the abilint section
305 
306   implicit none
307 
308   !Arguments ------------------------------------
309   !scalars
310   integer,intent(out) :: pwind_alloc
311   type(MPI_type),intent(inout) :: mpi_enreg
312   type(dataset_type),intent(inout) :: dtset
313   type(orbmag_type),intent(out) :: dtorbmag
314   type(pseudopotential_type),intent(in) :: psps
315   !arrays
316   integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt)
317   integer,intent(in) :: symrec(3,3,dtset%nsym)
318   integer,pointer :: pwind(:,:,:)
319   real(dp),intent(in) :: gmet(3,3),gprimd(3,3),occ(dtset%mband*dtset%nkpt*dtset%nsppol)
320   real(dp),intent(in) :: rprimd(3,3),xred(3,dtset%natom)
321   real(dp),pointer :: pwnsfac(:,:)
322   type(pawtab_type),intent(in) :: pawtab(dtset%ntypat)
323 
324   !Local variables-------------------------------
325   !scalars
326   integer :: brav,exchn2n3d,fnkpt_computed
327   integer :: iband,icg,icprj,idir,idum,idum1,ierr,ifor,ikg,ikg1
328   integer :: ikpt,ikpt_loc,ikpti,ikpt1,ikpt1f,ikpt1i
329   integer :: index,ipw,ipwnsfac,isign,isppol,istwf_k,isym,isym1,itrs,itypat
330   integer :: jpw,lmax,lmn2_size_max
331   integer :: mband_occ_k,me,me_g0,mkmem_,mkpt,my_nspinor,nband_k,nkptlatt,nproc,npw_k,npw_k1
332   integer :: option,spaceComm
333   real(dp) :: diffk1,diffk2,diffk3,ecut_eff
334   real(dp) :: kpt_shifted1,kpt_shifted2,kpt_shifted3,rdum
335   character(len=500) :: message
336   !arrays
337   integer :: iadum(3),iadum1(3),dg(3)
338   integer,allocatable :: kg1_k(:,:)
339   real(dp) :: diffk(3),dk(3),dum33(3,3),kpt1(3),tsec(2)
340   real(dp),allocatable :: spkpt(:,:)
341 
342   ! *************************************************************************
343 
344   DBG_ENTER("COLL")
345 
346   call timab(1001,1,tsec)
347   call timab(1002,1,tsec)
348 
349   !save the current value of nspinor
350   dtorbmag%nspinor = dtset%nspinor
351 
352   !----------------------------------------------------------------------------
353   !-------------------- Obtain k-point grid in the full BZ --------------------
354   !----------------------------------------------------------------------------
355 
356   if(dtset%kptopt==1 .or. dtset%kptopt==2 .or. dtset%kptopt==4)then
357      !  Compute the number of k points in the G-space unit cell
358      nkptlatt=dtset%kptrlatt(1,1)*dtset%kptrlatt(2,2)*dtset%kptrlatt(3,3) &
359           &   +dtset%kptrlatt(1,2)*dtset%kptrlatt(2,3)*dtset%kptrlatt(3,1) &
360           &   +dtset%kptrlatt(1,3)*dtset%kptrlatt(2,1)*dtset%kptrlatt(3,2) &
361           &   -dtset%kptrlatt(1,2)*dtset%kptrlatt(2,1)*dtset%kptrlatt(3,3) &
362           &   -dtset%kptrlatt(1,3)*dtset%kptrlatt(2,2)*dtset%kptrlatt(3,1) &
363           &   -dtset%kptrlatt(1,1)*dtset%kptrlatt(2,3)*dtset%kptrlatt(3,2)
364 
365      !  Call smpbz to obtain the list of k-point in the full BZ - without symmetry reduction
366      option = 0
367      brav = 1
368      mkpt=nkptlatt*dtset%nshiftk
369      ABI_ALLOCATE(spkpt,(3,mkpt))
370      call smpbz(1,ab_out,dtset%kptrlatt,mkpt,fnkpt_computed,dtset%nshiftk,option,dtset%shiftk,spkpt)
371      dtorbmag%fnkpt = fnkpt_computed
372      ABI_ALLOCATE(dtorbmag%fkptns,(3,dtorbmag%fnkpt))
373      dtorbmag%fkptns(:,:)=spkpt(:,1:dtorbmag%fnkpt)
374      ABI_DEALLOCATE(spkpt)
375   else if(dtset%kptopt==3.or.dtset%kptopt==0)then
376      dtorbmag%fnkpt=dtset%nkpt
377      ABI_ALLOCATE(dtorbmag%fkptns,(3,dtorbmag%fnkpt))
378      dtorbmag%fkptns(1:3,1:dtorbmag%fnkpt)=dtset%kpt(1:3,1:dtorbmag%fnkpt)
379      if(dtset%kptopt==0)then
380         write(message,'(10a)') ch10,&
381              &     ' initorbmag : WARNING -',ch10,&
382              &     '  you have defined manually the k-point grid with kptopt = 0',ch10,&
383              &     '  the orbital magnetization calculation works only with a regular k-points grid,',ch10,&
384              &     '  abinit doesn''t check if your grid is regular...'
385         call wrtout(std_out,message,'PERS')
386      end if
387   end if
388 
389   !call listkk to get mapping from FBZ to IBZ
390   rdum=1.0d-5  ! cutoff distance to decide when two k points match
391   ABI_ALLOCATE(dtorbmag%indkk_f2ibz,(dtorbmag%fnkpt,6))
392 
393   my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
394 
395   !JWZ: The following may need modification in the future
396   !**** no spin-polarization doubling ; do not allow use of time reversal symmetry ****
397 
398   call timab(1002,2,tsec)
399   call timab(1003,1,tsec)
400 
401   call listkk(rdum,gmet,dtorbmag%indkk_f2ibz,dtset%kptns,dtorbmag%fkptns,dtset%nkpt,&
402        & dtorbmag%fnkpt,dtset%nsym,1,dtset%symafm,symrec,0,use_symrec=.True.)
403 
404   call timab(1003,2,tsec)
405   call timab(1004,1,tsec)
406 
407   !Construct i2fbz and f2ibz
408   ABI_ALLOCATE(dtorbmag%i2fbz,(dtset%nkpt))
409   idum=0
410   do ikpt=1,dtorbmag%fnkpt
411      if (dtorbmag%indkk_f2ibz(ikpt,2)==1 .and. &
412           &   dtorbmag%indkk_f2ibz(ikpt,6) == 0 .and. &
413           &   maxval(abs(dtorbmag%indkk_f2ibz(ikpt,3:5))) == 0 ) then
414         dtorbmag%i2fbz(dtorbmag%indkk_f2ibz(ikpt,1))=ikpt
415         idum=idum+1
416      end if
417   end do
418   if (idum/=dtset%nkpt)then
419      message = ' Found wrong number of k-points in IBZ'
420      MSG_ERROR(message)
421   end if
422 
423   !----------------------------------------------------------------------------
424   !------------- Allocate PAW space as necessary ------------------------------
425   !----------------------------------------------------------------------------
426 
427   dtorbmag%usepaw   = psps%usepaw
428   dtorbmag%natom    = dtset%natom
429   dtorbmag%my_natom = mpi_enreg%my_natom
430 
431   ABI_ALLOCATE(dtorbmag%lmn_size,(dtset%ntypat))
432   ABI_ALLOCATE(dtorbmag%lmn2_size,(dtset%ntypat))
433   do itypat = 1, dtset%ntypat
434      dtorbmag%lmn_size(itypat) = pawtab(itypat)%lmn_size
435      dtorbmag%lmn2_size(itypat) = pawtab(itypat)%lmn2_size
436   end do
437 
438   lmn2_size_max = psps%lmnmax*(psps%lmnmax+1)/2
439   dtorbmag%lmn2max = lmn2_size_max
440 
441   ABI_ALLOCATE(dtorbmag%cprjindex,(dtset%nkpt,dtset%nsppol))
442   dtorbmag%cprjindex(:,:) = 0
443 
444   if (dtset%kptopt /= 3) then
445      ABI_ALLOCATE(dtorbmag%atom_indsym,(4,dtset%nsym,dtorbmag%natom))
446      call symatm(dtorbmag%atom_indsym,dtorbmag%natom,dtset%nsym,symrec,dtset%tnons,tol8,dtset%typat,xred)
447      lmax = psps%mpsang - 1
448      ABI_ALLOCATE(dtorbmag%zarot,(2*lmax+1,2*lmax+1,lmax+1,dtset%nsym))
449      call setsym_ylm(gprimd,lmax,dtset%nsym,1,rprimd,symrec,dtorbmag%zarot)
450      dtorbmag%nsym = dtset%nsym
451      dtorbmag%lmax = lmax
452      dtorbmag%lmnmax = psps%lmnmax
453   end if
454 
455   ! !------------------------------------------------------------------------------
456   ! !------------------- Compute variables related to MPI // ----------------------
457   ! !------------------------------------------------------------------------------
458   spaceComm=mpi_enreg%comm_cell
459   nproc=xmpi_comm_size(spaceComm)
460   me=xmpi_comm_rank(spaceComm)
461 
462   if (nproc==1) then
463      dtorbmag%fmkmem = dtorbmag%fnkpt
464      dtorbmag%fmkmem_max = dtorbmag%fnkpt
465      dtorbmag%mkmem_max = dtset%nkpt
466   else
467      dtorbmag%fmkmem = 0
468      do ikpt = 1, dtorbmag%fnkpt
469         ikpti = dtorbmag%indkk_f2ibz(ikpt,1)
470         nband_k = dtset%nband(ikpti)
471         if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,-1,me))) &
472              &     dtorbmag%fmkmem = dtorbmag%fmkmem + 1
473      end do
474      !  Maximum value of mkmem and fmkmem
475      call xmpi_max(dtorbmag%fmkmem,dtorbmag%fmkmem_max,spaceComm,ierr)
476      !  I have to use the dummy variable mkmem_ because
477      !  mkmem is declared as intent(in) while the first
478      !  argument of xmpi_max must be intent(inout)
479      mkmem_ = dtset%mkmem
480      call xmpi_max(mkmem_,dtorbmag%mkmem_max,spaceComm,ierr)
481   end if
482 
483   ABI_ALLOCATE(mpi_enreg%kpt_loc2fbz_sp,(0:nproc-1,1:dtorbmag%fmkmem_max*dtset%nsppol, 1:2))
484   ABI_ALLOCATE(mpi_enreg%kpt_loc2ibz_sp,(0:nproc-1,1:dtorbmag%mkmem_max*dtset%nsppol, 1:2))
485   ABI_ALLOCATE(mpi_enreg%kptdstrb,(nproc,6,dtorbmag%fmkmem_max*dtset%nsppol*2))
486   ABI_ALLOCATE(mpi_enreg%mkmem,(0:nproc-1))
487   mpi_enreg%kpt_loc2fbz_sp(:,:,:) = 0
488   mpi_enreg%kpt_loc2ibz_sp(:,:,:) = 0
489   mpi_enreg%kptdstrb(:,:,:)       = 0
490   mpi_enreg%mkmem(:)              = 0
491 
492   pwind_alloc = dtset%mpw*dtorbmag%fmkmem_max
493 
494   ABI_ALLOCATE(pwind,(pwind_alloc,2,3))
495   ABI_ALLOCATE(pwnsfac,(2,pwind_alloc))
496 
497   ! !------------------------------------------------------------------------------
498   ! !---------------------- Compute orbmag_type variables -------------------------
499   ! !------------------------------------------------------------------------------
500 
501   !Initialization of orbmag_type variables
502   dtorbmag%dkvecs(:,:) = zero
503   ABI_ALLOCATE(dtorbmag%ikpt_dk,(dtorbmag%fnkpt,2,3))
504   ABI_ALLOCATE(dtorbmag%cgindex,(dtset%nkpt,dtset%nsppol))
505   ABI_ALLOCATE(dtorbmag%kgindex,(dtset%nkpt))
506   ABI_ALLOCATE(dtorbmag%fkgindex,(dtorbmag%fnkpt))
507   dtorbmag%ikpt_dk(:,:,:) = 0
508   dtorbmag%cgindex(:,:) = 0
509   dtorbmag%mband_occ = 0
510   ABI_ALLOCATE(dtorbmag%nband_occ,(dtset%nsppol))
511   dtorbmag%kgindex(:) = 0
512   dtorbmag%fkgindex(:) = 0
513 
514   !Compute spin degeneracy
515   if (dtset%nsppol == 1 .and. dtset%nspinor == 1) then
516      dtorbmag%sdeg = two
517   else if (dtset%nsppol == 2 .or. my_nspinor == 2) then
518      dtorbmag%sdeg = one
519   end if
520 
521   !Compute the number of occupied bands and check that
522   !it is the same for each k-point
523 
524   index = 0
525   do isppol = 1, dtset%nsppol
526      dtorbmag%nband_occ(isppol) = 0
527      do ikpt = 1, dtset%nkpt
528 
529         mband_occ_k = 0
530         nband_k = dtset%nband(ikpt + (isppol - 1)*dtset%nkpt)
531 
532         do iband = 1, nband_k
533            index = index + 1
534            if (abs(occ(index) - dtorbmag%sdeg) < tol8) mband_occ_k = mband_occ_k + 1
535         end do
536 
537         if (ikpt > 1) then
538            if (dtorbmag%nband_occ(isppol) /= mband_occ_k) then
539               message = "The number of valence bands is not the same for every k-point of present spin channel"
540               MSG_ERROR(message)
541            end if
542         else
543            dtorbmag%mband_occ         = max(dtorbmag%mband_occ, mband_occ_k)
544            dtorbmag%nband_occ(isppol) = mband_occ_k
545         end if
546 
547      end do                ! close loop over ikpt
548   end do                ! close loop over isppol
549 
550   !Compute the location of each wavefunction
551 
552   icg = 0
553   icprj = 0
554   !ikg = 0
555   do isppol = 1, dtset%nsppol
556      do ikpt = 1, dtset%nkpt
557 
558         nband_k = dtset%nband(ikpt + (isppol-1)*dtset%nkpt)
559 
560         if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle
561 
562         dtorbmag%cgindex(ikpt,isppol) = icg
563         npw_k = npwarr(ikpt)
564         icg = icg + npw_k*dtorbmag%nspinor*nband_k
565 
566         if (psps%usepaw == 1) then
567            dtorbmag%cprjindex(ikpt,isppol) = icprj
568            icprj = icprj + dtorbmag%nspinor*nband_k
569         end if
570 
571      end do
572   end do
573 
574   ikg = 0
575   do ikpt = 1, dtset%nkpt
576      if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,me)).and.&
577           &   (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,dtset%nsppol,me))) cycle
578 
579      npw_k = npwarr(ikpt)
580      dtorbmag%kgindex(ikpt) = ikg
581      ikg = ikg + npw_k
582   end do
583 
584   call timab(1004,2,tsec)
585 
586   !------------------------------------------------------------------------------
587   !---------------------- Compute dk --------------------------------------------
588   !------------------------------------------------------------------------------
589 
590   call timab(1005,1,tsec)
591 
592   do idir = 1, 3
593 
594      !    Compute dk(:), the vector between a k-point and its nearest
595      !    neighbour along the direction idir
596 
597      dk(:) = zero
598      dk(idir) = 1._dp   ! 1 mean there is no other k-point un the direction idir
599      do ikpt = 2, dtorbmag%fnkpt
600         diffk(:) = abs(dtorbmag%fkptns(:,ikpt) - dtorbmag%fkptns(:,1))
601         if ((diffk(1) < dk(1)+tol8).and.(diffk(2) < dk(2)+tol8).and.&
602              &     (diffk(3) < dk(3)+tol8)) dk(:) = diffk(:)
603      end do
604      dtorbmag%dkvecs(:,idir) = dk(:)
605      !    DEBUG
606      !    write(std_out,*)' initorbmag : idir, dk', idir, dk
607      !    ENDDEBUG
608 
609      !    For each k point, find k_prim such that k_prim= k + dk mod(G)
610      !    where G is a vector of the reciprocal lattice
611 
612      do ikpt = 1, dtorbmag%fnkpt
613 
614         !      First k+dk, then k-dk
615         do isign=-1,1,2
616            kpt_shifted1=dtorbmag%fkptns(1,ikpt)- isign*dk(1)
617            kpt_shifted2=dtorbmag%fkptns(2,ikpt)- isign*dk(2)
618            kpt_shifted3=dtorbmag%fkptns(3,ikpt)- isign*dk(3)
619            !        Note that this is still a order fnkpt**2 algorithm.
620            !        It is possible to implement a order fnkpt algorithm, see listkk.F90.
621            do ikpt1 = 1, dtorbmag%fnkpt
622               diffk1=dtorbmag%fkptns(1,ikpt1) - kpt_shifted1
623               if(abs(diffk1-nint(diffk1))>tol8)cycle
624               diffk2=dtorbmag%fkptns(2,ikpt1) - kpt_shifted2
625               if(abs(diffk2-nint(diffk2))>tol8)cycle
626               diffk3=dtorbmag%fkptns(3,ikpt1) - kpt_shifted3
627               if(abs(diffk3-nint(diffk3))>tol8)cycle
628               dtorbmag%ikpt_dk(ikpt,(isign+3)/2,idir) = ikpt1
629               exit
630            end do   ! ikpt1
631         end do     ! isign
632 
633      end do     ! ikpt
634 
635   end do     ! close loop over idir
636 
637   call timab(1005,2,tsec)
638   call timab(1006,1,tsec)
639 
640   !------------------------------------------------------------------------------
641   !------------ Build the array pwind that is needed to compute the -------------
642   !------------ overlap matrices at k +- dk                         -------------
643   !------------------------------------------------------------------------------
644 
645   ecut_eff = dtset%ecut*(dtset%dilatmx)**2
646   exchn2n3d = 0 ; istwf_k = 1 ; ikg1 = 0
647   pwind(:,:,:) = 0
648   pwnsfac(1,:) = 1.0_dp
649   pwnsfac(2,:) = 0.0_dp
650   ABI_ALLOCATE(kg1_k,(3,dtset%mpw))
651 
652   ipwnsfac = 0
653 
654   do idir = 1, 3
655 
656      dk(:) = dtorbmag%dkvecs(:,idir)
657 
658      do ifor = 1, 2
659 
660         if (ifor == 2) dk(:) = -1._dp*dk(:)
661 
662         !      Build pwind and kgindex
663         !      NOTE: The array kgindex is important for parallel execution.
664         !      In case nsppol = 2, it may happen that a particular processor
665         !      treats k-points at different spin polarizations.
666         !      In this case, it is not possible to address the elements of
667         !      pwind correctly without making use of the kgindex array.
668 
669         ikg = 0 ; ikpt_loc = 0 ; isppol = 1
670         do ikpt = 1, dtorbmag%fnkpt
671 
672            ikpti = dtorbmag%indkk_f2ibz(ikpt,1)
673            nband_k = dtset%nband(ikpti)
674            ikpt1f = dtorbmag%ikpt_dk(ikpt,ifor,idir)
675            ikpt1i = dtorbmag%indkk_f2ibz(ikpt1f,1)
676 
677            if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,1,me)).and.&
678                 &       (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,dtset%nsppol,me))) cycle
679 
680            ikpt_loc = ikpt_loc + 1
681 
682            !        Build basis sphere of plane waves for the nearest neighbour of
683            !        the k-point (important for MPI //)
684 
685            kg1_k(:,:) = 0
686            kpt1(:) = dtset%kptns(:,ikpt1i)
687            call kpgsph(ecut_eff,exchn2n3d,gmet,ikg1,ikpt,istwf_k,kg1_k,kpt1,&
688                 &       1,mpi_enreg,dtset%mpw,npw_k1)
689            me_g0=mpi_enreg%me_g0
690 
691 
692            !        ji: fkgindex is defined here !
693            dtorbmag%fkgindex(ikpt) = ikg
694 
695            !
696            !        Deal with symmetry transformations
697            !
698 
699            !        bra k-point k(b) and IBZ k-point kIBZ(b) related by
700            !        k(b) = alpha(b) S(b)^t kIBZ(b) + G(b)
701            !        where alpha(b), S(b) and G(b) are given by indkk_f2ibz
702            !
703            !        For the ket k-point:
704            !        k(k) = alpha(k) S(k)^t kIBZ(k) + G(k) - GBZ(k)
705            !        where GBZ(k) takes k(k) to the BZ
706            !
707 
708            isym  = dtorbmag%indkk_f2ibz(ikpt,2)
709            isym1 = dtorbmag%indkk_f2ibz(ikpt1f,2)
710 
711            !        Construct transformed G vector that enters the matching condition:
712            !        alpha(k) S(k)^{t,-1} ( -G(b) - GBZ(k) + G(k) )
713 
714            dg(:) = -dtorbmag%indkk_f2ibz(ikpt,3:5) &
715                 &       -nint(-dtorbmag%fkptns(:,ikpt) - dk(:) - tol10 &
716                 &       +dtorbmag%fkptns(:,ikpt1f)) &
717                 &       +dtorbmag%indkk_f2ibz(ikpt1f,3:5)
718 
719            iadum(:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),dg(:))
720 
721            dg(:) = iadum(:)
722 
723            if ( dtorbmag%indkk_f2ibz(ikpt1f,6) == 1 ) dg(:) = -dg(:)
724 
725            !        Construct S(k)^{t,-1} S(b)^{t}
726 
727            dum33(:,:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),symrec(:,:,isym))
728 
729            !        Construct alpha(k) alpha(b)
730 
731            if (dtorbmag%indkk_f2ibz(ikpt,6) == dtorbmag%indkk_f2ibz(ikpt1f,6)) then
732               itrs=0
733            else
734               itrs=1
735            end if
736 
737 
738            npw_k  = npwarr(ikpti)
739            !        npw_k1 = npwarr(ikpt1i)
740 
741            !        loop over bra G vectors
742            do ipw = 1, npw_k
743 
744               !          NOTE: the bra G vector is taken for the sym-related IBZ k point,
745               !          not for the FBZ k point
746               iadum(:) = kg(:,dtorbmag%kgindex(ikpti) + ipw)
747 
748               !          Store non-symmorphic operation phase factor exp[i2\pi \alpha G \cdot t]
749 
750               if ( ipwnsfac == 0 ) then
751                  rdum=0.0_dp
752                  do idum=1,3
753                     rdum=rdum+dble(iadum(idum))*dtset%tnons(idum,isym)
754                  end do
755                  rdum=two_pi*rdum
756                  if ( dtorbmag%indkk_f2ibz(ikpt,6) == 1 ) rdum=-rdum
757                  pwnsfac(1,ikg+ipw) = cos(rdum)
758                  pwnsfac(2,ikg+ipw) = sin(rdum)
759               end if
760 
761               !          to determine r.l.v. matchings, we transformed the bra vector
762               !          Rotation
763               iadum1(:)=0
764               do idum1=1,3
765                  iadum1(:)=iadum1(:)+dum33(:,idum1)*iadum(idum1)
766               end do
767               iadum(:)=iadum1(:)
768               !          Time reversal
769               if (itrs==1) iadum(:)=-iadum(:)
770               !          Translation
771               iadum(:) = iadum(:) + dg(:)
772 
773               do jpw = 1, npw_k1
774                  iadum1(1:3) = kg1_k(1:3,jpw)
775                  if ( (iadum(1) == iadum1(1)).and. &
776                       &           (iadum(2) == iadum1(2)).and. &
777                       &           (iadum(3) == iadum1(3)) ) then
778                     pwind(ikg + ipw,ifor,idir) = jpw
779                     !              write(std_out,'(a,2x,3i4,2x,i4)') 'Found !:',iadum1(:),jpw
780                     exit
781                  end if
782               end do
783            end do
784 
785            ikg  = ikg + npw_k
786 
787         end do    ! close loop over ikpt
788 
789         ipwnsfac = 1
790 
791      end do    ! close loop over ifor
792 
793   end do        ! close loop over idir
794 
795 
796   call timab(1008,2,tsec)
797   call timab(1009,1,tsec)
798 
799   !Build mpi_enreg%kptdstrb
800   !array required to communicate the WFs between cpus
801   !(MPI // over k-points)
802   if (nproc>1) then
803      do idir = 1, 3
804         do ifor = 1, 2
805 
806            ikpt_loc = 0
807            do isppol = 1, dtset%nsppol
808 
809               do ikpt = 1, dtorbmag%fnkpt
810 
811                  ikpti = dtorbmag%indkk_f2ibz(ikpt,1)
812                  nband_k = dtset%nband(ikpti)
813                  ikpt1f = dtorbmag%ikpt_dk(ikpt,ifor,idir)
814                  ikpt1i = dtorbmag%indkk_f2ibz(ikpt1f,1)
815 
816                  if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,isppol,me)) cycle
817 
818                  ikpt_loc = ikpt_loc + 1
819                  mpi_enreg%kptdstrb(me + 1,ifor+2*(idir-1),ikpt_loc) = &
820                       &           ikpt1i + (isppol - 1)*dtset%nkpt
821 
822                  mpi_enreg%kptdstrb(me+1,ifor+2*(idir-1),ikpt_loc+dtorbmag%fmkmem_max*dtset%nsppol) = &
823                       &           ikpt1f + (isppol - 1)*dtorbmag%fnkpt
824 
825               end do   ! ikpt
826            end do     ! isppol
827         end do       ! ifor
828      end do           ! idir
829   end if             ! nproc>1
830 
831   !build mpi_enreg%kpt_loc2fbz_sp
832   ikpt_loc = 0
833   do isppol = 1, dtset%nsppol
834      do ikpt = 1, dtorbmag%fnkpt
835 
836         ikpti = dtorbmag%indkk_f2ibz(ikpt,1)
837         nband_k = dtset%nband(ikpti)
838 
839         if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,isppol,me)) cycle
840 
841         ikpt_loc = ikpt_loc + 1
842 
843         mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc, 1) = ikpt
844         mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc, 2) = isppol
845 
846      end do
847   end do
848 
849   !should be temporary
850   !unassigned mpi_enreg%kpt_loc2fbz_sp are empty ; inform other cpu (there are better ways...)
851   mpi_enreg%mkmem(me) = dtset%mkmem
852   !do ii=ikpt_loc+1,dtefield%fmkmem_max
853   !mpi_enreg%kpt_loc2fbz_sp(me, ii, 1) = -1
854   !end do
855 
856   call xmpi_sum(mpi_enreg%kptdstrb,spaceComm,ierr)
857   call xmpi_sum(mpi_enreg%kpt_loc2fbz_sp,spaceComm,ierr)
858 
859   ABI_DEALLOCATE(kg1_k)
860 
861 
862   call timab(1009,2,tsec)
863   call timab(1001,2,tsec)
864 
865   DBG_EXIT("COLL")
866 
867 end subroutine initorbmag

ABINIT/m_orbmag [ Modules ]

[ Top ] [ Modules ]

NAME

  m_orbmag

FUNCTION

  This module contains the declaration of data types and methods
  used to handle orbital magnetization

COPYRIGHT

 Copyright (C) 2011-2017 ABINIT group (JWZ)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

OUTPUT

NOTES

PARENTS

CHILDREN

SOURCE

28 #if defined HAVE_CONFIG_H
29 #include "config.h"
30 #endif
31 
32 #include "abi_common.h"
33 
34 module m_orbmag
35 
36   use defs_abitypes
37   use defs_basis
38   use defs_datatypes
39   use m_errors
40   use m_abicore
41   use m_xmpi
42 
43   use m_berrytk,          only : smatrix
44   use m_cgprj,            only : getcprj
45   use m_fft,              only : fftpac
46   use m_fftcore,          only : kpgsph
47   use m_fourier_interpol, only : transgrid
48   use m_geometry,         only : metric
49   use m_getghc,           only : getghc
50   use m_hamiltonian,      only : init_hamiltonian,destroy_hamiltonian,&
51        &                         load_spin_hamiltonian,load_k_hamiltonian,gs_hamiltonian_type
52   use m_initylmg,         only : initylmg
53   use m_kg,               only : getph,mkkin,mkpwind_k,mknucdipmom_k,ph1d3d
54   use m_kpts,             only : listkk, smpbz
55   use m_mkffnl,           only : mkffnl
56   use m_mpinfo,           only : proc_distrb_cycle
57   use m_pawang,           only : pawang_type
58   use m_pawfgr,           only : pawfgr_type
59   use m_paw_ij,           only : paw_ij_type
60   use m_paw_overlap,      only : overlap_k1k2_paw
61   use m_pawrad,           only : pawrad_type
62   use m_paw_sphharm,      only : initylmr,setsym_ylm
63   use m_pawtab,           only : pawtab_type
64   use m_pawcprj,          only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_free,&
65        &                         pawcprj_get, pawcprj_put, pawcprj_getdim, pawcprj_set_zero
66   use m_symtk,            only : symatm
67   use m_time,             only : timab
68 
69 
70   implicit none
71 
72   private

ABINIT/orbmag [ Functions ]

[ Top ] [ Functions ]

NAME

 orbmag

FUNCTION

 This routine computes the orbital magnetization based on input wavefunctions.
 It is assumed that only completely filled bands are present.

COPYRIGHT

 Copyright (C) 2003-2017 ABINIT  group
 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

 atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
 cg(2,mcg)=planewave coefficients of wavefunctions
 cprj(natom,mcprj*usecrpj)=<p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
 dtset <type(dataset_type)>=all input variables in this dataset
 kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere
 mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
 mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol
 mpi_enreg=information about MPI parallelization
 nfftf= - PAW only - number of FFT grid points for the "fine" grid (see NOTES at beginning of scfcv)
 npwarr(nkpt)=number of planewaves in basis at this k point
 paw_ij(my_natom*usepaw) <type(paw_ij_type)>=paw arrays given on (i,j) channels
 pawang <type(pawang_type)>=paw angular mesh and related data
 pawfgr <type(pawfgr_type)>=fine grid parameters and related data
 pawrad(ntypat*psps%usepaw) <type(pawrad_type)>=paw radial mesh and related data
 pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
 psps <type(pseudopotential_type)>=variables related to pseudopotentials
 pwind(pwind_alloc,2,3) = array used to compute
           the overlap matrix smat between k-points (see initberry.f)
 pwind_alloc = first dimension of pwind
 rprimd(3,3)=dimensional primitive translations in real space (bohr)
 symrec(3,3,nsym) = symmetries in reciprocal space in terms of
   reciprocal space primitive translations
 usecprj=1 if cprj datastructure has been allocated
 vhartr(nfftf)=Hartree potential
 vpsp(nfftf)=array for holding local psp
 vxc(nfftf,nspden)=exchange-correlation potential (hartree) in real space
 xred(3,natom) = location of atoms in unit cell
 ylm(mpw*mkmem,mpsang*mpsang*useylm)= real spherical harmonics for each G and k point
 ylmgr(mpw*mkmem,3,mpsang*mpsang*useylm)= gradients of real spherical harmonics

OUTPUT

SIDE EFFECTS

 dtorbmag <type(orbmag_type)> = variables related to orbital magnetization

TODO

NOTES

 See Ceresoli et al, PRB 74, 024408 (2006) [[cite:Ceresoli2006]],
 and Gonze and Zwanziger, PRB 84, 064445 (2011) [[cite:Gonze2011a]].
 The derivative of the density operator is obtained from a discretized formula
 $\partial_\beta \rho_k = \frac{1}{2\Delta}(\rho_{k+b} - \rho_{k-b})$ with
 $\Delta = |b|$. When reduced to wavefunction overlaps the computation amounts to
 multiple calls to smatrix.F90, exactly as in other Berry phase computations, with
 the one additional complication of overlaps like $\langle u_{n,k+b}|u_{n',k+g}\rangle$.
 At this stage mkpwind_k is invoked, which generalizes the code in initberry
 and initorbmag necessary to index plane waves around different k points.
 Direct questions and comments to J Zwanziger

PARENTS

CHILDREN

SOURCE

1376 subroutine orbmag(atindx1,cg,cprj,dtset,dtorbmag,kg,&
1377      &            mcg,mcprj,mpi_enreg,nattyp,nfftf,npwarr,paw_ij,pawang,pawfgr,pawrad,pawtab,psps,&
1378      &            pwind,pwind_alloc,rprimd,symrec,usecprj,vhartr,vpsp,vxc,xred,ylm,ylmgr)
1379 
1380 
1381 !This section has been created automatically by the script Abilint (TD).
1382 !Do not modify the following lines by hand.
1383 #undef ABI_FUNC
1384 #define ABI_FUNC 'orbmag'
1385 !End of the abilint section
1386 
1387  implicit none
1388 
1389  !Arguments ------------------------------------
1390  !scalars
1391  integer,intent(in) :: mcg,mcprj,nfftf,pwind_alloc,usecprj
1392  type(dataset_type),intent(in) :: dtset
1393  type(MPI_type), intent(inout) :: mpi_enreg
1394  type(orbmag_type), intent(inout) :: dtorbmag
1395  type(pawang_type),intent(in) :: pawang
1396  type(pawfgr_type),intent(in) :: pawfgr
1397  type(pseudopotential_type),intent(in) :: psps
1398 
1399  !arrays
1400  integer,intent(in) :: atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem),nattyp(dtset%ntypat)
1401  integer,intent(in) :: npwarr(dtset%nkpt),pwind(pwind_alloc,2,3),symrec(3,3,dtset%nsym)
1402  real(dp),intent(in) :: cg(2,mcg),rprimd(3,3)
1403  real(dp),intent(in) :: vhartr(nfftf),vpsp(nfftf),vxc(nfftf,dtset%nspden),xred(3,dtset%natom)
1404  real(dp),intent(in) :: ylm(dtset%mpw*dtset%mkmem,psps%mpsang*psps%mpsang*psps%useylm)
1405  real(dp),intent(in) :: ylmgr(dtset%mpw*dtset%mkmem,3,psps%mpsang*psps%mpsang*psps%useylm)
1406  type(paw_ij_type),intent(inout) :: paw_ij(dtset%natom*psps%usepaw)
1407  type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*psps%usepaw)
1408  type(pawcprj_type),intent(in) ::  cprj(dtset%natom,mcprj*usecprj)
1409  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
1410 
1411  !Local variables -------------------------
1412  !scalars
1413  integer :: adir,bdir,bdx,bdxc,bfor,bsigma,cpopt,ddkflag,dimffnl,epsabg
1414  integer :: gdir,gdx,gdxc,gfor,gsigma
1415  integer :: iatom,icg,icgb,icgg,icprj,icprjb,icprjg,ider,idir
1416  integer :: ikg,ikgb,ikgg,ikpt,ikptb,ikptg,ilm,ilmn,ipw,isppol,istwf_k,itrs,itypat
1417  integer :: jlmn,job,jpw
1418  integer :: klmn,mcg1_k,my_cpopt,my_nspinor,nband_k,ncpgr,ndat,nkpg,nn,n1,n2
1419  integer :: ngfft1,ngfft2,ngfft3,ngfft4,ngfft5,ngfft6,npw_k,npw_kb,npw_kg
1420  integer :: prtvol,shiftbd,sij_opt,tim_getghc,type_calc,type_calc_123
1421  real(dp) :: deltab,deltag,dkg2,dotr,doti,htpisq,keg,lambda,ucvol
1422  complex(dpc) :: cdij,cgdijcb,cpb,cpg
1423  complex(dpc) :: IIA,IIA1,IIA2,IIA3,IIIA,IIIA1,IIIA2,IIIA3,tprodIIA,tprodIIIA
1424  character(len=500) :: message
1425  type(gs_hamiltonian_type) :: gs_hamk,gs_hamk123
1426  !arrays
1427  integer,allocatable :: dimlmn(:),kg_k(:,:),kg_kb(:,:),kg_kg(:,:)
1428  integer,allocatable :: pwind_kb(:),pwind_kg(:),pwind_bg(:),sflag_k(:)
1429  real(dp) :: dkb(3),dkg(3),dkbg(3),dtm_k(2),gmet(3,3),gprimd(3,3)
1430  real(dp) :: kpoint(3),kpointb(3),kpointg(3)
1431  real(dp) :: orbmagvec(2,3),rhodum(1),rmet(3,3)
1432  real(dp),allocatable :: bra(:,:),cg1_k(:,:),cgrvtrial(:,:),cwavef(:,:),ffnl(:,:,:,:),ghc(:,:),gsc(:,:),gvnlc(:,:)
1433  real(dp),allocatable :: hmat(:,:,:,:,:,:),kinpw(:),kk_paw(:,:,:),kpg_k_dummy(:,:)
1434  real(dp),allocatable :: my_nucdipmom(:,:),ph3d(:,:,:),pwnsfac_k(:,:),smat_all(:,:,:,:,:,:),smat_inv(:,:,:)
1435  real(dp),allocatable :: smat_kk(:,:,:),vlocal(:,:,:,:),vtrial(:,:),ylm_k(:,:)
1436  complex(dpc),allocatable :: nucdipmom_k(:)
1437  logical,allocatable :: has_hmat(:,:,:),has_smat(:,:,:)
1438  type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kb(:,:),cprj_kb_k(:,:,:,:)
1439  type(pawcprj_type),allocatable :: cprj_kg(:,:),cwaveprj(:,:)
1440 
1441  ! ***********************************************************************
1442  ! my_nspinor=max(1,dtorbmag%nspinor/mpi_enreg%nproc_spinor)
1443 
1444  ! TODO: generalize to nsppol > 1
1445  isppol = 1
1446  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
1447 
1448  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1449 
1450  nband_k = dtorbmag%mband_occ
1451 
1452  if (psps%usepaw == 1) then ! cprj allocation
1453    ncpgr = cprj(1,1)%ncpgr
1454    ABI_ALLOCATE(dimlmn,(dtset%natom))
1455    call pawcprj_getdim(dimlmn,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'R')
1456 
1457    ABI_DATATYPE_ALLOCATE(cprj_k,(dtset%natom,dtorbmag%nspinor*dtset%mband))
1458    call pawcprj_alloc(cprj_k,ncpgr,dimlmn)
1459 
1460    ABI_DATATYPE_ALLOCATE(cprj_kb,(dtset%natom,dtorbmag%nspinor*dtset%mband))
1461    call pawcprj_alloc(cprj_kb,ncpgr,dimlmn)
1462 
1463    ABI_DATATYPE_ALLOCATE(cprj_kg,(dtset%natom,dtorbmag%nspinor*dtset%mband))
1464    call pawcprj_alloc(cprj_kg,ncpgr,dimlmn)
1465 
1466    ABI_DATATYPE_ALLOCATE(cwaveprj,(dtset%natom,1))
1467    call pawcprj_alloc(cwaveprj,ncpgr,dimlmn)
1468 
1469    ABI_DATATYPE_ALLOCATE(cprj_kb_k,(dtorbmag%fnkpt,6,dtset%natom,dtorbmag%nspinor*dtset%mband))
1470    do ikpt=1,dtorbmag%fnkpt
1471       do bdx=1, 6
1472          call pawcprj_alloc(cprj_kb_k(ikpt,bdx,:,:),ncpgr,dimlmn)
1473       end do
1474    end do
1475  else
1476    message = ' usepaw /= 1 but orbital magnetization calculation requires PAW '
1477    MSG_ERROR(message)
1478  end if
1479 
1480  ABI_ALLOCATE(kk_paw,(2,dtset%mband,dtset%mband))
1481  ABI_ALLOCATE(pwind_kb,(dtset%mpw))
1482  ABI_ALLOCATE(pwind_kg,(dtset%mpw))
1483  ABI_ALLOCATE(pwind_bg,(dtset%mpw))
1484  ABI_ALLOCATE(pwnsfac_k,(4,dtset%mpw))
1485  pwnsfac_k(1,:) = 1.0_dp ! bra real
1486  pwnsfac_k(2,:) = 0.0_dp ! bra imag
1487  pwnsfac_k(3,:) = 1.0_dp ! ket real
1488  pwnsfac_k(4,:) = 0.0_dp ! ket imag
1489 
1490  mcg1_k = dtset%mpw*nband_k
1491  ABI_ALLOCATE(cg1_k,(2,mcg1_k))
1492  ABI_ALLOCATE(sflag_k,(nband_k))
1493  ABI_ALLOCATE(smat_inv,(2,nband_k,nband_k))
1494  ABI_ALLOCATE(smat_kk,(2,nband_k,nband_k))
1495 
1496  ABI_ALLOCATE(kinpw,(dtset%mpw))
1497 
1498  ABI_ALLOCATE(my_nucdipmom,(3,dtset%natom))
1499  my_nucdipmom(:,:) = dtset%nucdipmom(:,:)
1500 
1501  ! input parameters for calls to smatrix.F90
1502  ddkflag = 1
1503  istwf_k = 1
1504  ! itrs = 0 means do not invoke time reversal symmetry in smatrix.F90
1505  itrs = 0
1506  job = 1
1507  shiftbd = 1
1508 
1509  ngfft1=dtset%ngfft(1) ; ngfft2=dtset%ngfft(2) ; ngfft3=dtset%ngfft(3)
1510  ngfft4=dtset%ngfft(4) ; ngfft5=dtset%ngfft(5) ; ngfft6=dtset%ngfft(6)
1511 
1512  ! input parameters for calls to getghc at ikpt
1513  cpopt = -1
1514  ndat = 1
1515  prtvol = 0
1516  sij_opt = 0
1517  tim_getghc = 0
1518  ! getghc: type_calc 0 means kinetic, local, nonlocal
1519  type_calc = 0
1520  lambda = zero
1521 
1522  htpisq = 0.5_dp*(two_pi)**2
1523 
1524  ABI_ALLOCATE(has_smat,(dtorbmag%fnkpt,0:6,0:6))
1525  ABI_ALLOCATE(smat_all,(2,nband_k,nband_k,dtorbmag%fnkpt,0:6,0:6))
1526  has_smat(:,:,:)=.FALSE.
1527  ABI_ALLOCATE(has_hmat,(dtorbmag%fnkpt,0:6,0:6))
1528  ABI_ALLOCATE(hmat,(2,nband_k,nband_k,dtorbmag%fnkpt,0:6,0:6))
1529  has_hmat(:,:,:) = .FALSE.
1530 
1531  !==== Initialize most of the Hamiltonian ====
1532  !Allocate all arrays and initialize quantities that do not depend on k and spin.
1533  !gs_hamk is the normal hamiltonian at k, needed for computing E_nk
1534 !  call init_hamiltonian(gs_hamk,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%natom,&
1535 ! & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,nucdipmom=dtset%nucdipmom,&
1536 ! & paw_ij=paw_ij)
1537  call init_hamiltonian(gs_hamk,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%natom,&
1538       & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,nucdipmom=my_nucdipmom,&
1539       & paw_ij=paw_ij)
1540 
1541  !gs_hamk123 is used to apply vlocal in <u_nk1|Hk2|u_mk3>
1542  ! my_nucdipmom can be used to override the input nuclear dipoles
1543  ! call init_hamiltonian(gs_hamk123,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%natom,&
1544  !      & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,nucdipmom=dtset%nucdipmom)
1545  call init_hamiltonian(gs_hamk123,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%natom,&
1546       & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,nucdipmom=my_nucdipmom)
1547 
1548  !---------construct local potential------------------
1549  ABI_ALLOCATE(vtrial,(nfftf,dtset%nspden))
1550  ! nspden=1 is essentially hard-coded in the following line
1551  vtrial(1:nfftf,1)=vhartr(1:nfftf)+vxc(1:nfftf,1)+vpsp(1:nfftf)
1552 
1553  ABI_ALLOCATE(cgrvtrial,(dtset%nfft,dtset%nspden))
1554  call transgrid(1,mpi_enreg,dtset%nspden,-1,0,0,dtset%paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial)
1555 
1556  ABI_ALLOCATE(vlocal,(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc))
1557  call fftpac(isppol,mpi_enreg,dtset%nspden,ngfft1,ngfft2,ngfft3,ngfft4,ngfft5,ngfft6,dtset%ngfft,cgrvtrial,vlocal,2)
1558 
1559  ABI_DEALLOCATE(cgrvtrial)
1560  ABI_DEALLOCATE(vtrial)
1561 
1562  ! same vlocal in both Hamiltonians
1563  call load_spin_hamiltonian(gs_hamk,isppol,vlocal=vlocal,with_nonlocal=.true.)
1564  call load_spin_hamiltonian(gs_hamk123,isppol,vlocal=vlocal,with_nonlocal=.false.)
1565 
1566  !------- now local potential is attached to gs_hamk and gs_hamk123 -------------------------
1567 
1568  ! compute the shifted cprj's <p_k+b|u_k>
1569  call ctocprjb(atindx1,cg,cprj_kb_k,dtorbmag,dtset,gmet,gprimd,&
1570       & istwf_k,kg,mcg,mpi_enreg,nattyp,ncpgr,npwarr,pawtab,psps,rmet,rprimd,ucvol,xred)
1571 
1572  do adir = 1, 3
1573 
1574     IIA  = czero
1575     IIIA = czero
1576 
1577     do epsabg = 1, -1, -2
1578 
1579        if (epsabg .EQ. 1) then
1580           bdir = modulo(adir,3)+1
1581           gdir = modulo(adir+1,3)+1
1582        else
1583           bdir = modulo(adir+1,3)+1
1584           gdir = modulo(adir,3)+1
1585        end if
1586 
1587        ! loop over kpts, assuming for now kptopt 3, nsppol = 1, nspinor = 1
1588        ! and no parallelism, no symmorphic symmetry elements
1589        do ikpt = 1, dtorbmag%fnkpt
1590 
1591           kpoint(:)=dtorbmag%fkptns(:,ikpt)
1592           icprj = dtorbmag%cprjindex(ikpt,isppol)
1593 
1594           npw_k = npwarr(ikpt)
1595           icg = dtorbmag%cgindex(ikpt,dtset%nsppol)
1596 
1597           ikg = dtorbmag%fkgindex(ikpt)
1598           ABI_ALLOCATE(kg_k,(3,npw_k))
1599           kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
1600 
1601           call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,icprj,ikpt,0,isppol,dtset%mband,&
1602                &       dtset%mkmem,dtset%natom,nband_k,nband_k,my_nspinor,dtset%nsppol,0)
1603 
1604           ! Set up remainder of normal Hamiltonian at k if necessary
1605 
1606           if (.NOT. has_hmat(ikpt,0,0) ) then
1607              ! kpoint(:)=dtorbmag%fkptns(:,ikpt)
1608 
1609              ABI_ALLOCATE(ylm_k,(npw_k,psps%mpsang*psps%mpsang*psps%useylm))
1610              if (psps%useylm==1) then
1611                 do ilm=1,psps%mpsang*psps%mpsang
1612                    ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm)
1613                 end do
1614              end if
1615 
1616              !      Compute (1/2) (2 Pi)**2 (k+G)**2:
1617              kinpw(:) = zero
1618              call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k,0,0)
1619 
1620              !  Compute (k+G) vectors (only if useylm=1)
1621              ! original code from vtorho.F90
1622              ! nkpg=3*optforces*dtset%nloalg(3)
1623              ! ABI_ALLOCATE(kpg_k,(npw_k,nkpg))
1624              ! if ((mpi_enreg%paral_kgb/=1.or.istep<=1).and.nkpg>0) then
1625              !    call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k)
1626              ! end if
1627              ! pretty sure do not need k+g vectors, use dummy kpg_k
1628              ! may eventually need them for nucdipmom_k hamiltonian if
1629              ! generalize to k,k'
1630              nkpg = 0
1631              ABI_ALLOCATE(kpg_k_dummy,(npw_k,nkpg))
1632 
1633              !      Compute nonlocal form factors ffnl at all (k+G):
1634              ider=0;idir=0;dimffnl=1
1635              ABI_ALLOCATE(ffnl,(npw_k,dimffnl,psps%lmnmax,dtset%ntypat))
1636              call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl,psps%ffspl,&
1637                   &         gmet,gprimd,ider,idir,psps%indlmn,kg_k,kpg_k_dummy,kpoint,psps%lmnmax,&
1638                   &         psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,&
1639                   &         npw_k,dtset%ntypat,psps%pspso,psps%qgrid_ff,rmet,&
1640                   &         psps%usepaw,psps%useylm,ylm_k,ylmgr)
1641 
1642              !     compute and load nuclear dipole Hamiltonian at current k point
1643              if(any(abs(gs_hamk%nucdipmom)>0.0)) then
1644                 if(allocated(nucdipmom_k)) then
1645                    ABI_DEALLOCATE(nucdipmom_k)
1646                 end if
1647                 ABI_ALLOCATE(nucdipmom_k,(npw_k*(npw_k+1)/2))
1648                 call mknucdipmom_k(gmet,kg_k,kpoint,dtset%natom,gs_hamk%nucdipmom,&
1649                      &           nucdipmom_k,npw_k,rprimd,ucvol,xred)
1650                 if(allocated(gs_hamk%nucdipmom_k)) then
1651                    ABI_DEALLOCATE(gs_hamk%nucdipmom_k)
1652                 end if
1653                 ABI_ALLOCATE(gs_hamk%nucdipmom_k,(npw_k*(npw_k+1)/2))
1654                 call load_k_hamiltonian(gs_hamk,nucdipmom_k=nucdipmom_k)
1655                 ABI_DEALLOCATE(nucdipmom_k)
1656              end if
1657 
1658 
1659              !      Load k-dependent part in the Hamiltonian datastructure
1660              !       - Compute 3D phase factors
1661              !       - Prepare various tabs in case of band-FFT parallelism
1662              !       - Load k-dependent quantities in the Hamiltonian
1663              ABI_ALLOCATE(ph3d,(2,npw_k,gs_hamk%matblk))
1664 
1665              call load_k_hamiltonian(gs_hamk,kpt_k=kpoint(:),istwf_k=istwf_k,npw_k=npw_k,&
1666                   &         kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k_dummy,ffnl_k=ffnl,ph3d_k=ph3d,&
1667                   &         compute_ph3d=.TRUE.,compute_gbound=.TRUE.)
1668 
1669              ! apply gs_hamk to wavefunctions at k to compute E_nk eigenvalues
1670              ABI_ALLOCATE(cwavef,(2,npw_k))
1671              ABI_ALLOCATE(ghc,(2,npw_k))
1672              ABI_ALLOCATE(gsc,(2,npw_k))
1673              ABI_ALLOCATE(gvnlc,(2,npw_k))
1674              hmat(:,:,:,ikpt,0,0) = zero
1675              do nn = 1, nband_k
1676                 cwavef(1,1:npw_k) = cg(1,icg+(nn-1)*npw_k+1:icg+nn*npw_k)
1677                 cwavef(2,1:npw_k) = cg(2,icg+(nn-1)*npw_k+1:icg+nn*npw_k)
1678                 call pawcprj_get(atindx1,cwaveprj,cprj_k,dtset%natom,nn,0,ikpt,0,isppol,dtset%mband,&
1679                      &           dtset%mkmem,dtset%natom,1,nband_k,my_nspinor,dtset%nsppol,0)
1680                 my_cpopt=2
1681                 call getghc(my_cpopt,cwavef,cwaveprj,ghc,gsc,gs_hamk,gvnlc,lambda,mpi_enreg,ndat,&
1682                      &           prtvol,sij_opt,tim_getghc,type_calc)
1683                 hmat(1,nn,nn,ikpt,0,0)= DOT_PRODUCT(cwavef(1,1:npw_k),ghc(1,1:npw_k)) &
1684                      &           + DOT_PRODUCT(cwavef(2,1:npw_k),ghc(2,1:npw_k))
1685              end do
1686              has_hmat(ikpt,0,0) = .TRUE.
1687 
1688              ABI_DEALLOCATE(cwavef)
1689              ABI_DEALLOCATE(ghc)
1690              ABI_DEALLOCATE(gsc)
1691              ABI_DEALLOCATE(gvnlc)
1692 
1693              ABI_DEALLOCATE(ylm_k)
1694              ABI_DEALLOCATE(kpg_k_dummy)
1695              ABI_DEALLOCATE(ffnl)
1696              if(any(abs(gs_hamk%nucdipmom)>0.0)) then
1697                 if(allocated(nucdipmom_k)) then
1698                    ABI_DEALLOCATE(nucdipmom_k)
1699                 end if
1700              end if
1701              ABI_DEALLOCATE(ph3d)
1702 
1703           end if ! end check on has_hmat
1704 
1705           do bfor = 1, 2
1706              if (bfor .EQ. 1) then
1707                 bsigma = 1
1708              else
1709                 bsigma = -1
1710              end if
1711              ! index of neighbor 1..6
1712              bdx = 2*bdir-2+bfor
1713              ! index of ikpt viewed from neighbor
1714              bdxc = 2*bdir-2+bfor+bsigma
1715 
1716              dkb(1:3) = bsigma*dtorbmag%dkvecs(1:3,bdir)
1717              deltab = sqrt(DOT_PRODUCT(dkb,MATMUL(gmet,dkb)))
1718 
1719              ikptb = dtorbmag%ikpt_dk(ikpt,bfor,bdir)
1720              kpointb(:)=dtorbmag%fkptns(:,ikptb)
1721 
1722              icprjb = dtorbmag%cprjindex(ikptb,isppol)
1723 
1724              ikgb = dtorbmag%fkgindex(ikptb)
1725              npw_kb = npwarr(ikptb)
1726              ABI_ALLOCATE(kg_kb,(3,npw_kb))
1727              kg_kb(:,1:npw_kb)=kg(:,ikgb+1:ikgb+npw_kb)
1728 
1729              icgb = dtorbmag%cgindex(ikptb,dtset%nsppol)
1730 
1731              pwind_kb(1:npw_k) = pwind(ikg+1:ikg+npw_k,bfor,bdir)
1732 
1733              call pawcprj_get(atindx1,cprj_kb,cprj,dtset%natom,1,icprjb,&
1734                   &         ikptb,0,isppol,dtset%mband,dtset%mkmem,dtset%natom,nband_k,nband_k,&
1735                   &         my_nspinor,dtset%nsppol,0)
1736 
1737              if (.NOT. has_smat(ikpt,bdx,0)) then
1738 
1739                 call overlap_k1k2_paw(cprj_k,cprj_kb,dkb,gprimd,kk_paw,dtorbmag%lmn2max,dtorbmag%lmn_size,dtset%mband,&
1740                      &           dtset%natom,my_nspinor,dtset%ntypat,pawang,pawrad,pawtab,dtset%typat,xred)
1741 
1742                 sflag_k=0
1743                 call smatrix(cg,cg,cg1_k,ddkflag,dtm_k,icg,icgb,itrs,job,nband_k,&
1744                      &           mcg,mcg,mcg1_k,1,dtset%mpw,nband_k,nband_k,npw_k,npw_kb,my_nspinor,&
1745                      &           pwind_kb,pwnsfac_k,sflag_k,shiftbd,smat_inv,smat_kk,kk_paw,psps%usepaw)
1746 
1747                 do nn = 1, nband_k
1748                    do n1 = 1, nband_k
1749                       smat_all(1,nn,n1,ikpt,bdx,0) =  smat_kk(1,nn,n1)
1750                       smat_all(2,nn,n1,ikpt,bdx,0) =  smat_kk(2,nn,n1)
1751                       smat_all(1,n1,nn,ikptb,bdxc,0) =  smat_kk(1,nn,n1)
1752                       smat_all(2,n1,nn,ikptb,bdxc,0) = -smat_kk(2,nn,n1)
1753                    end do
1754                 end do
1755 
1756                 has_smat(ikpt,bdx,0) = .TRUE.
1757                 has_smat(ikptb,bdxc,0) = .TRUE.
1758 
1759              end if
1760 
1761              do gfor = 1, 2
1762                 if (gfor .EQ. 1) then
1763                    gsigma = 1
1764                 else
1765                    gsigma = -1
1766                 end if
1767                 ! index of neighbor 1..6
1768                 gdx = 2*gdir-2+gfor
1769                 ! index of ikpt viewed from neighbor
1770                 gdxc = 2*gdir-2+gfor+gsigma
1771 
1772                 dkg(1:3) = gsigma*dtorbmag%dkvecs(1:3,gdir)
1773                 deltag = sqrt(DOT_PRODUCT(dkg,MATMUL(gmet,dkg)))
1774 
1775                 ikptg = dtorbmag%ikpt_dk(ikpt,gfor,gdir)
1776                 kpointg(:)=dtorbmag%fkptns(:,ikptg)
1777 
1778                 icprjg = dtorbmag%cprjindex(ikptg,isppol)
1779 
1780                 npw_kg = npwarr(ikptg)
1781                 ABI_ALLOCATE(kg_kg,(3,npw_kg))
1782                 ikgg = dtorbmag%fkgindex(ikptg)
1783                 kg_kg(:,1:npw_kg)=kg(:,ikgg+1:ikgg+npw_kg)
1784 
1785                 icgg = dtorbmag%cgindex(ikptg,dtset%nsppol)
1786 
1787                 pwind_kg(1:npw_k) = pwind(ikg+1:ikg+npw_k,gfor,gdir)
1788 
1789                 call pawcprj_get(atindx1,cprj_kg,cprj,dtset%natom,1,icprjg,&
1790                      &           ikptg,0,isppol,dtset%mband,dtset%mkmem,dtset%natom,nband_k,nband_k,&
1791                      &           my_nspinor,dtset%nsppol,0)
1792 
1793                 if (.NOT. has_smat(ikpt,gdx,0)) then
1794 
1795                    call overlap_k1k2_paw(cprj_k,cprj_kg,dkg,gprimd,kk_paw,dtorbmag%lmn2max,dtorbmag%lmn_size,dtset%mband,&
1796                         &             dtset%natom,my_nspinor,dtset%ntypat,pawang,pawrad,pawtab,dtset%typat,xred)
1797 
1798                    sflag_k=0
1799                    call smatrix(cg,cg,cg1_k,ddkflag,dtm_k,icg,icgg,itrs,job,nband_k,&
1800                         &             mcg,mcg,mcg1_k,1,dtset%mpw,nband_k,nband_k,npw_k,npw_kg,my_nspinor,&
1801                         &             pwind_kg,pwnsfac_k,sflag_k,shiftbd,smat_inv,smat_kk,kk_paw,psps%usepaw)
1802 
1803                    do nn = 1, nband_k
1804                       do n1 = 1, nband_k
1805                          smat_all(1,nn,n1,ikpt,gdx,0) =  smat_kk(1,nn,n1)
1806                          smat_all(2,nn,n1,ikpt,gdx,0) =  smat_kk(2,nn,n1)
1807                          smat_all(1,n1,nn,ikptg,gdxc,0) =  smat_kk(1,nn,n1)
1808                          smat_all(2,n1,nn,ikptg,gdxc,0) = -smat_kk(2,nn,n1)
1809                       end do
1810                    end do
1811 
1812                    has_smat(ikpt,gdx,0) = .TRUE.
1813                    has_smat(ikptg,gdxc,0) = .TRUE.
1814 
1815                 end if
1816 
1817                 dkbg = dkg - dkb
1818 
1819                 if (.NOT. has_smat(ikpt,bdx,gdx)) then
1820 
1821                    call overlap_k1k2_paw(cprj_kb,cprj_kg,dkbg,gprimd,kk_paw,dtorbmag%lmn2max,dtorbmag%lmn_size,dtset%mband,&
1822                         &             dtset%natom,my_nspinor,dtset%ntypat,pawang,pawrad,pawtab,dtset%typat,xred)
1823 
1824                    call mkpwind_k(dkbg,dtset,dtorbmag%fnkpt,dtorbmag%fkptns,gmet,&
1825                         &             dtorbmag%indkk_f2ibz,ikptb,ikptg,&
1826                         &             kg,dtorbmag%kgindex,mpi_enreg,npw_kb,pwind_bg,symrec)
1827 
1828                    sflag_k=0
1829                    call smatrix(cg,cg,cg1_k,ddkflag,dtm_k,icgb,icgg,itrs,job,nband_k,&
1830                         &             mcg,mcg,mcg1_k,1,dtset%mpw,nband_k,nband_k,npw_kb,npw_kg,my_nspinor,&
1831                         &             pwind_bg,pwnsfac_k,sflag_k,shiftbd,smat_inv,smat_kk,kk_paw,psps%usepaw)
1832 
1833                    do nn = 1, nband_k
1834                       do n1 = 1, nband_k
1835                          smat_all(1,nn,n1,ikpt,bdx,gdx) =  smat_kk(1,nn,n1)
1836                          smat_all(2,nn,n1,ikpt,bdx,gdx) =  smat_kk(2,nn,n1)
1837                          smat_all(1,n1,nn,ikpt,gdx,bdx) =  smat_kk(1,nn,n1)
1838                          smat_all(2,n1,nn,ikpt,gdx,bdx) = -smat_kk(2,nn,n1)
1839                       end do
1840                    end do
1841 
1842                    has_smat(ikpt,bdx,gdx) = .TRUE.
1843                    has_smat(ikpt,gdx,bdx) = .TRUE.
1844 
1845                 end if
1846 
1847                 if (.NOT. has_hmat(ikpt,gdx,bdx)) then
1848 
1849                    call mkpwind_k(-dkbg,dtset,dtorbmag%fnkpt,dtorbmag%fkptns,gmet,&
1850                         &             dtorbmag%indkk_f2ibz,ikptg,ikptb,&
1851                         &             kg,dtorbmag%kgindex,mpi_enreg,npw_kg,pwind_bg,symrec)
1852 
1853                    nkpg = 0
1854                    ABI_ALLOCATE(kpg_k_dummy,(npw_kb,nkpg))
1855 
1856                    !     compute and load nuclear dipole Hamiltonian at current k point
1857                    ! this may need to be modified to take into account "twist"
1858                    if(any(abs(gs_hamk123%nucdipmom)>0.0)) then
1859                       if(allocated(nucdipmom_k)) then
1860                          ABI_DEALLOCATE(nucdipmom_k)
1861                       end if
1862                       ABI_ALLOCATE(nucdipmom_k,(npw_kb*(npw_kb+1)/2))
1863                       call mknucdipmom_k(gmet,kg_kb,kpointb,dtset%natom,gs_hamk123%nucdipmom,&
1864                            &           nucdipmom_k,npw_kb,rprimd,ucvol,xred)
1865                       if(allocated(gs_hamk123%nucdipmom_k)) then
1866                          ABI_DEALLOCATE(gs_hamk123%nucdipmom_k)
1867                       end if
1868                       ABI_ALLOCATE(gs_hamk123%nucdipmom_k,(npw_kb*(npw_kb+1)/2))
1869                       call load_k_hamiltonian(gs_hamk123,nucdipmom_k=nucdipmom_k)
1870                       ABI_DEALLOCATE(nucdipmom_k)
1871                    end if
1872 
1873                    ! this is minimal Hamiltonian information, to apply vlocal (and only vlocal) to |u_kb>
1874                    call load_k_hamiltonian(gs_hamk123,kpt_k=kpointb(:),istwf_k=istwf_k,npw_k=npw_kb,&
1875                         &             kg_k=kg_kb,kpg_k=kpg_k_dummy,compute_gbound=.TRUE.)
1876 
1877 
1878                    ! apply gs_hamk123 to wavefunctions at kb
1879                    ABI_ALLOCATE(cwavef,(2,npw_kb))
1880                    ABI_ALLOCATE(ghc,(2,npw_kb))
1881                    ABI_ALLOCATE(gsc,(2,npw_kb))
1882                    ABI_ALLOCATE(gvnlc,(2,npw_kb))
1883 
1884                    ABI_ALLOCATE(bra,(2,npw_kg))
1885 
1886                    ! getghc: type_calc_123 1 means local only, 3 means kinetic, local only
1887                    type_calc_123 = 1
1888                    dkg2=DOT_PRODUCT(dkg(:),MATMUL(gmet(:,:),dkg(:)))
1889 
1890                    do nn = 1, nband_k
1891                       cwavef(1,1:npw_kb) = cg(1,icgb+(nn-1)*npw_kb+1:icgb+nn*npw_kb)
1892                       cwavef(2,1:npw_kb) = cg(2,icgb+(nn-1)*npw_kb+1:icgb+nn*npw_kb)
1893                       ! apply only vlocal
1894                       call getghc(cpopt,cwavef,cwaveprj,ghc,gsc,gs_hamk123,gvnlc,lambda,mpi_enreg,ndat,&
1895                            &               prtvol,sij_opt,tim_getghc,type_calc_123)
1896                       do n1 = 1, nband_k
1897                          bra(1,1:npw_kg) = cg(1,icgg+(n1-1)*npw_kg+1:icgg+n1*npw_kg)
1898                          bra(2,1:npw_kg) = cg(2,icgg+(n1-1)*npw_kg+1:icgg+n1*npw_kg)
1899                          dotr=zero;doti=zero
1900 
1901                          ! apply local potential through ghc
1902                          do ipw=1,npw_kg
1903                             jpw=pwind_bg(ipw)
1904                             if(jpw .GT. 0) then
1905                                dotr=dotr+bra(1,ipw)*ghc(1,jpw)+bra(2,ipw)*ghc(2,jpw)
1906                                doti=doti+bra(1,ipw)*ghc(2,jpw)-bra(2,ipw)*ghc(1,jpw)
1907 
1908                                ! twisted kinetic energy: here we are computing
1909                                ! -\frac{1}{2}<u_kg|e^{-i.k.r}\nabla^2 e^{i.k.r}|u_kb>,
1910                                ! that is, kinetic energy at k between wavefunctions at kg and kb. The correct
1911                                ! formula is htpisq*(ikpt + G_right)^2\delta(G_left,G_right) but it's hard to apply
1912                                ! because the G's have wrap-around shifts (output of mkpwind_k)
1913                                ! for the kpts near and at the edge of the IBZ.
1914                                ! The following approach is based on the bra <u_kg| expansion, because
1915                                ! these G vectors are unshifted (indexed by ipw, not jpw). So we are using
1916                                ! k+G_left = (k-kg) + (kg+G_left) = -dkg + (kg+G_left). When squared we obtain
1917                                ! |kg+G_left|^2 - 2*dkg*(kg+G_left) + |dkg|^2. In this way we only use the G_left
1918                                ! expansion vectors, with no shift, for each k point.
1919 
1920                                ! normal kinetic energy for bra
1921                                keg=htpisq*dot_product((kpointg(:)+kg_kg(:,ipw)),MATMUL(gmet,(kpointg(:)+kg_kg(:,ipw))))
1922 
1923                                ! addition of |dkg|^2
1924                                keg=keg+htpisq*dkg2
1925 
1926                                ! addition of -2*dkg*(kg+G_left)
1927                                keg=keg-2.0*htpisq*DOT_PRODUCT(dkg(:),MATMUL(gmet,(kpointg(:)+kg_kg(:,ipw))))
1928 
1929                                ! application of ecut filter and wavefunction
1930                                if (keg < dtset%ecut) then
1931                                   dotr=dotr+bra(1,ipw)*keg*cwavef(1,jpw)+bra(2,ipw)*keg*cwavef(2,jpw)
1932                                   doti=doti+bra(1,ipw)*keg*cwavef(2,jpw)-bra(2,ipw)*keg*cwavef(1,jpw)
1933                                end if ! end keg filter
1934                             end if ! end check on jpw > 0
1935                          end do ! end loop over ipw
1936 
1937                          ! apply onsite terms through cprjk+b
1938                          cgdijcb = czero
1939                          do iatom = 1, dtset%natom
1940                             itypat = dtset%typat(iatom)
1941                             do ilmn = 1, pawtab(itypat)%lmn_size
1942                                cpg=cmplx(cprj_kb_k(ikptg,gdxc,iatom,n1)%cp(1,ilmn),&
1943                                     & cprj_kb_k(ikptg,gdxc,iatom,n1)%cp(2,ilmn))
1944                                do jlmn = 1, pawtab(itypat)%lmn_size
1945                                   cpb=cmplx(cprj_kb_k(ikptb,bdxc,iatom,nn)%cp(1,jlmn),&
1946                                        & cprj_kb_k(ikptb,bdxc,iatom,nn)%cp(2,jlmn))
1947                                   if (jlmn .LE. ilmn) then
1948                                      klmn = (ilmn-1)*ilmn/2 + jlmn
1949                                   else
1950                                      klmn = (jlmn-1)*jlmn/2 + ilmn
1951                                   end if
1952                                   if (paw_ij(iatom)%cplex_dij .EQ. 2) then
1953                                      cdij=cmplx(paw_ij(iatom)%dij(2*klmn-1,1),paw_ij(iatom)%dij(2*klmn,1))
1954                                      if (jlmn .GT. ilmn) cdij=conjg(cdij)
1955                                   else
1956                                      cdij=cmplx(paw_ij(iatom)%dij(klmn,1),zero)
1957                                   end if
1958                                   cgdijcb = cgdijcb + conjg(cpg)*cdij*cpb
1959                                end do
1960                             end do
1961                          end do
1962                          hmat(1,n1,nn,ikpt,gdx,bdx) = dotr + real(cgdijcb)
1963                          hmat(2,n1,nn,ikpt,gdx,bdx) = doti + aimag(cgdijcb)
1964                          ! hmat(1,n1,nn,ikpt,gdx,bdx) = dotr
1965                          ! hmat(2,n1,nn,ikpt,gdx,bdx) = doti
1966                          ! hmat(1,nn,n1,ikpt,bdx,gdx) = dotr
1967                          ! hmat(2,nn,n1,ikpt,bdx,gdx) = -doti
1968                       end do ! end loop over n1
1969                    end do ! end loop over nn
1970                    has_hmat(ikpt,gdx,bdx) = .TRUE.
1971                    ! has_hmat(ikpt,bdx,gdx) = .TRUE.
1972 
1973                    ABI_DEALLOCATE(cwavef)
1974                    ABI_DEALLOCATE(bra)
1975                    ABI_DEALLOCATE(ghc)
1976                    ABI_DEALLOCATE(gsc)
1977                    ABI_DEALLOCATE(gvnlc)
1978 
1979                    ABI_DEALLOCATE(kpg_k_dummy)
1980 
1981                 end if
1982 
1983                 do nn = 1, nband_k
1984                    do n1 = 1, nband_k
1985 
1986                       IIA1  = cmplx(smat_all(1,nn,n1,ikpt,gdx,0),smat_all(2,nn,n1,ikpt,gdx,0))
1987 
1988                       IIIA1 = cmplx(smat_all(1,nn,n1,ikpt,bdx,0),smat_all(2,nn,n1,ikpt,bdx,0))
1989 
1990                       do n2 = 1, nband_k
1991 
1992                          IIA2 = cmplx(hmat(1,n1,n2,ikpt,gdx,bdx),hmat(2,n1,n2,ikpt,gdx,bdx))
1993                          IIA3 = cmplx(smat_all(1,n2,nn,ikptb,bdxc,0),smat_all(2,n2,nn,ikptb,bdxc,0))
1994 
1995                          IIIA2 = cmplx(smat_all(1,n1,n2,ikpt,bdx,gdx),smat_all(2,n1,n2,ikpt,bdx,gdx))
1996                          IIIA3 = conjg(cmplx(smat_all(1,nn,n2,ikpt,gdx,0),smat_all(2,nn,n2,ikpt,gdx,0)))
1997 
1998                          tprodIIA  = IIA1*IIA2*IIA3
1999                          IIA = IIA + epsabg*bsigma*gsigma*tprodIIA/(2.0*deltab*2.0*deltag)
2000 
2001                          tprodIIIA = hmat(1,nn,nn,ikpt,0,0)*IIIA1*IIIA2*IIIA3
2002                          IIIA = IIIA - epsabg*bsigma*gsigma*tprodIIIA/(2.0*deltab*2.0*deltag)
2003 
2004                       end do ! end n2
2005                    end do ! end n1
2006                 end do ! end nn
2007 
2008                 ABI_DEALLOCATE(kg_kg)
2009 
2010              end do ! end gfor
2011 
2012              ABI_DEALLOCATE(kg_kb)
2013 
2014           end do ! end bfor
2015 
2016           ABI_DEALLOCATE(kg_k)
2017 
2018        end do ! end loop over fnkpt
2019 
2020     end do ! end loop over epsabg
2021 
2022     orbmagvec(1,adir) = real(IIA+IIIA)
2023     orbmagvec(2,adir) = aimag(IIA+IIIA)
2024     ! orbmagvec(1,adir) = real(IIIA)
2025     ! orbmagvec(2,adir) = aimag(IIIA)
2026  end do ! end loop over adir
2027 
2028  orbmagvec(1,1:3) = MATMUL(gprimd,orbmagvec(1,1:3))
2029  orbmagvec(2,1:3) = MATMUL(gprimd,orbmagvec(2,1:3))
2030  dtorbmag%orbmagvec(1,1:3) =  orbmagvec(2,1:3)/(two*dtorbmag%fnkpt)
2031  dtorbmag%orbmagvec(2,1:3) = -orbmagvec(1,1:3)/(two*dtorbmag%fnkpt)
2032 
2033  write(message,'(a,a,a)')ch10,'====================================================',ch10
2034  call wrtout(ab_out,message,'COLL')
2035 
2036  write(message,'(a)')' Orbital magnetization '
2037  call wrtout(ab_out,message,'COLL')
2038  write(message,'(a,a)')'----Orbital magnetization is a real vector, given along Cartesian directions----',ch10
2039  call wrtout(ab_out,message,'COLL')
2040 
2041  do adir = 1, 3
2042    write(message,'(a,i4,a,2es16.8)')' Orb Mag(',adir,') : real, imag ',&
2043 &   dtorbmag%orbmagvec(1,adir),dtorbmag%orbmagvec(2,adir)
2044    call wrtout(ab_out,message,'COLL')
2045  end do
2046 
2047  write(message,'(a,a,a)')ch10,'====================================================',ch10
2048  call wrtout(ab_out,message,'COLL')
2049 
2050 
2051  if (psps%usepaw == 1) then
2052    ABI_DEALLOCATE(dimlmn)
2053    call pawcprj_free(cprj_k)
2054    ABI_DATATYPE_DEALLOCATE(cprj_k)
2055    call pawcprj_free(cprj_kb)
2056    ABI_DATATYPE_DEALLOCATE(cprj_kb)
2057    call pawcprj_free(cprj_kg)
2058    ABI_DATATYPE_DEALLOCATE(cprj_kg)
2059    call pawcprj_free(cwaveprj)
2060    ABI_DATATYPE_DEALLOCATE(cwaveprj)
2061    do ikpt=1,dtorbmag%fnkpt
2062       do bdx = 1, 6
2063          call pawcprj_free(cprj_kb_k(ikpt,bdx,:,:))
2064       end do
2065    end do
2066    ABI_DATATYPE_DEALLOCATE(cprj_kb_k)
2067  end if
2068 
2069  ABI_DEALLOCATE(kk_paw)
2070  ABI_DEALLOCATE(cg1_k)
2071  ABI_DEALLOCATE(sflag_k)
2072  ABI_DEALLOCATE(smat_inv)
2073  ABI_DEALLOCATE(smat_kk)
2074  ABI_DEALLOCATE(pwind_kb)
2075  ABI_DEALLOCATE(pwind_kg)
2076  ABI_DEALLOCATE(pwind_bg)
2077  ABI_DEALLOCATE(pwnsfac_k)
2078 
2079  ABI_DEALLOCATE(kinpw)
2080 
2081  ABI_DEALLOCATE(has_smat)
2082  ABI_DEALLOCATE(smat_all)
2083  ABI_DEALLOCATE(has_hmat)
2084  ABI_DEALLOCATE(hmat)
2085 
2086  ABI_DEALLOCATE(my_nucdipmom)
2087 
2088  ABI_DEALLOCATE(vlocal)
2089  call destroy_hamiltonian(gs_hamk)
2090  call destroy_hamiltonian(gs_hamk123)
2091 
2092 end subroutine orbmag

m_orbmag/destroy_orbmag [ Functions ]

[ Top ] [ m_orbmag ] [ Functions ]

NAME

FUNCTION

   deallocate fields in orbmag structure

INPUTS

OUTPUT

SOURCE

189   subroutine destroy_orbmag(dtorbmag)
190 
191 
192 !This section has been created automatically by the script Abilint (TD).
193 !Do not modify the following lines by hand.
194 #undef ABI_FUNC
195 #define ABI_FUNC 'destroy_orbmag'
196 !End of the abilint section
197 
198     implicit none
199 
200     !Arguments ------------------------------------
201     !array
202     type(orbmag_type),intent(inout) :: dtorbmag
203 
204     ! ************************************************************************
205 
206     ! Integer pointers
207     if(allocated(dtorbmag%atom_indsym))  then
208        ABI_DEALLOCATE(dtorbmag%atom_indsym)
209     end if
210     if(allocated(dtorbmag%cgindex))  then
211        ABI_DEALLOCATE(dtorbmag%cgindex)
212     end if
213     if(allocated(dtorbmag%cprjindex))  then
214        ABI_DEALLOCATE(dtorbmag%cprjindex)
215     end if
216     if(allocated(dtorbmag%fkgindex))  then
217        ABI_DEALLOCATE(dtorbmag%fkgindex)
218     end if
219     if(allocated(dtorbmag%ikpt_dk))  then
220        ABI_DEALLOCATE(dtorbmag%ikpt_dk)
221     end if
222     if(allocated(dtorbmag%indkk_f2ibz))  then
223        ABI_DEALLOCATE(dtorbmag%indkk_f2ibz)
224     end if
225     if(allocated(dtorbmag%i2fbz))  then
226        ABI_DEALLOCATE(dtorbmag%i2fbz)
227     end if
228     if(allocated(dtorbmag%kgindex))  then
229        ABI_DEALLOCATE(dtorbmag%kgindex)
230     end if
231     if(allocated(dtorbmag%lmn_size))  then
232        ABI_DEALLOCATE(dtorbmag%lmn_size)
233     end if
234     if(allocated(dtorbmag%lmn2_size))  then
235        ABI_DEALLOCATE(dtorbmag%lmn2_size)
236     end if
237     if(allocated(dtorbmag%nband_occ))  then
238        ABI_DEALLOCATE(dtorbmag%nband_occ)
239     end if
240     ! Real(dp) pointers
241 
242     if(allocated(dtorbmag%fkptns))  then
243        ABI_DEALLOCATE(dtorbmag%fkptns)
244     end if
245     if(allocated(dtorbmag%zarot))  then
246        ABI_DEALLOCATE(dtorbmag%zarot)
247     end if
248 
249   end subroutine destroy_orbmag

m_orbmag/orbmag_type [ Types ]

[ Top ] [ m_orbmag ] [ Types ]

NAME

 orbmag_type

FUNCTION

 variables used in orbital magnetism calculation

SOURCE

 85   type, public :: orbmag_type
 86 
 87 ! WARNING : if you modify this datatype, please check whether there might be creation/destruction/copy routines,
 88 ! declared in another part of ABINIT, that might need to take into account your modification.
 89 
 90 ! Integer variables
 91      integer :: orbmag              ! value of orbmag input variable in use
 92      integer :: fmkmem              ! number of k-points in the FBZ per cpu
 93      integer :: fmkmem_max          ! max of fmkmem
 94      integer :: fnkpt               ! number of k-points in the FBZ
 95      integer :: lmax
 96      integer :: lmnmax
 97      integer :: lmn2max
 98      integer :: mkmem_max           ! max of mkmem
 99      integer :: natom               ! number of atoms in unit cell
100      integer :: my_natom            ! number of atoms treated by current proc
101      integer :: mband_occ           ! max number of occupied bands (over spin)
102      ! this number must be the same for every k
103      integer :: nspinor             ! nspinor input from data set
104      integer :: nsym
105      integer :: usepaw              ! 1 if a PAW calculation, 0 else
106 
107      ! Real(dp) scalars
108      real(dp) :: sdeg               ! spin degeneracy: sdeg = 2 if nsppol = 1
109 
110      ! Real(dp) arrays
111      real(dp) :: chern(2,3)           ! result of chern number calculation
112 
113      real(dp) :: dkvecs(3,3)        ! dkvec(:,idir) = vector between a k-point and its nearest neighbour along idir
114 
115      real(dp) :: orbmagvec(2,3)     ! result of orbital magnetization calculation
116 
117      ! Integer pointers
118      integer, allocatable :: atom_indsym(:,:,:) ! atom_indsym(4,nsym,natom)
119      ! this is data on how the symmetries map the atoms in the cell
120      ! see symatm.F90 for full description
121      integer, allocatable :: cgindex(:,:)    ! cgindex(nkpt,nsppol)
122      ! for each k-point, stores the location
123      ! of the WF in the cg array
124      integer, allocatable :: cprjindex(:,:)  ! cprjindex(nkpt,nsppol)
125      ! for each k-point, stores the location
126      ! of the cprj in the cprj array (used only
127      ! for PAW calculations)
128      integer, allocatable :: fkgindex(:)     ! same as kgindex, but defined
129      ! for the FBZ and intended to use
130      ! with pwindf
131      integer, allocatable :: ikpt_dk(:,:,:)  ! ikpt_dk(nkpt,2,3)
132      ! ikpt_dp(ikpt,ii,idir) = index of the
133      ! k-point at k+dk (ii=1) and k-dk (ii=2)
134      integer, allocatable :: indkk_f2ibz(:,:)   ! indkk_f2ibz(1:dtorbmag%fnkpt,1:6)
135      ! information needed to fold a
136      ! k-point in the FBZ into the IBZ;
137      ! the second index (1:6)
138      ! is as described in listkk
139      integer, allocatable :: i2fbz(:)           ! i2fbz(1:nkpt) gives index of IBZ
140      ! k-points in the FBZ k-point list
141 
142      integer, allocatable :: kgindex(:)      ! kgind(nkpt)
143      ! kgind(ikpt) = ikg
144 
145      integer, allocatable :: lmn_size(:)        ! lmn_size(ntypat)
146      integer, allocatable :: lmn2_size(:)       ! lmn2_size(ntypat)
147 
148      integer, allocatable :: nband_occ(:)       ! nband_occ(nsppol) = actual number of occupied bands
149      !  can be different for spin up and down!!!
150      ! Real(dp) allocatables
151 
152      real(dp), allocatable :: fkptns(:,:)       ! fkptns(3,1:dtorbmag%fnkpt) k-points in FBZ
153 
154      real(dp), allocatable :: zarot(:,:,:,:)
155      !  zarot(l_size_max,l_size_max,l_max,nsym)
156      !  Coeffs of the transformation of real spherical
157      !  harmonics under the symmetry operations. These are needed when the
158      ! cprj's need to be computed in the full BZ, that is,
159      ! in the PAW case with kptopt /= 3.
160 
161      ! complex(dpc) allocatable
162 
163   end type orbmag_type
164 
165   ! Bound methods:
166   public :: destroy_orbmag
167   public :: initorbmag
168   public :: chern_number
169   public :: orbmag
170   public :: ctocprjb
171 
172 CONTAINS  !========================================================================================