TABLE OF CONTENTS
ABINIT/chern_number [ 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 (JZwanziger, MVeithen) 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), and Gonze and Zwanziger, PRB 84
064446 (2011). 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
afterscfloop
CHILDREN
mkpwind_k,overlap_k1k2_paw,pawcprj_alloc,pawcprj_free,pawcprj_get pawcprj_getdim,smatrix,wrtout
SOURCE
71 #if defined HAVE_CONFIG_H 72 #include "config.h" 73 #endif 74 75 #include "abi_common.h" 76 77 subroutine chern_number(atindx1,cg,cprj,dtset,dtorbmag,gmet,gprimd,kg,& 78 & mcg,mcprj,mpi_enreg,npwarr,pawang,pawrad,pawtab,pwind,pwind_alloc,& 79 & symrec,usecprj,usepaw,xred) 80 81 use defs_basis 82 use defs_abitypes 83 use defs_datatypes 84 use m_xmpi 85 use m_errors 86 use m_profiling_abi 87 88 use m_orbmag 89 90 use m_fftcore, only : kpgsph 91 use m_pawang, only : pawang_type 92 use m_pawrad, only : pawrad_type 93 use m_pawtab, only : pawtab_type 94 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_copy, pawcprj_free,& 95 & pawcprj_get, pawcprj_getdim, pawcprj_set_zero 96 97 !This section has been created automatically by the script Abilint (TD). 98 !Do not modify the following lines by hand. 99 #undef ABI_FUNC 100 #define ABI_FUNC 'chern_number' 101 use interfaces_14_hidewrite 102 use interfaces_32_util 103 use interfaces_56_recipspace 104 use interfaces_65_paw 105 !End of the abilint section 106 107 implicit none 108 109 !Arguments ------------------------------------ 110 !scalars 111 integer,intent(in) :: mcg,mcprj,pwind_alloc,usecprj,usepaw 112 type(dataset_type),intent(in) :: dtset 113 type(MPI_type), intent(inout) :: mpi_enreg 114 type(orbmag_type), intent(inout) :: dtorbmag 115 type(pawang_type),intent(in) :: pawang 116 117 !arrays 118 integer,intent(in) :: atindx1(dtset%natom),kg(3,dtset%mpw*dtset%mkmem) 119 integer,intent(in) :: npwarr(dtset%nkpt),pwind(pwind_alloc,2,3),symrec(3,3,dtset%nsym) 120 real(dp), intent(in) :: cg(2,mcg),gmet(3,3),gprimd(3,3),xred(3,dtset%natom) 121 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*usepaw) 122 type(pawcprj_type),intent(in) :: cprj(dtset%natom,mcprj*usecprj) 123 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*usepaw) 124 125 !Local variables ------------------------- 126 !scalars 127 integer :: adir,bdir,bfor,bsigma,ddkflag,epsabg,gdir,gfor,gsigma 128 integer :: icg,icgb,icgg,icprj,icprjb,icprjg 129 integer :: ikg,ikpt,ikptb,ikptg,isppol,itrs,job 130 integer :: mcg1_k,my_nspinor,nband_k,ncpgr,nn,n1,n2,n3,npw_k,npw_kb,npw_kg,shiftbd 131 real(dp) :: deltab,deltag 132 complex(dpc) :: IA,IB,t1A,t2A,t3A,t1B,t2B,t3B,t4B,tprodA,tprodB 133 character(len=500) :: message 134 !arrays 135 integer,allocatable :: dimlmn(:),nattyp_dum(:),pwind_kb(:),pwind_kg(:),pwind_bg(:),sflag_k(:) 136 real(dp) :: cnum(2,3),dkb(3),dkg(3),dkbg(3),dtm_k(2) 137 real(dp),allocatable :: cg1_k(:,:),kk_paw(:,:,:),pwnsfac_k(:,:),smat_all(:,:,:,:,:),smat_inv(:,:,:) 138 real(dp),allocatable :: smat_kk(:,:,:) 139 logical,allocatable :: has_smat(:,:) 140 type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_kb(:,:),cprj_kg(:,:) 141 type(pawcprj_type),allocatable :: cprj_fkn(:,:),cprj_ikn(:,:) 142 143 ! *********************************************************************** 144 ! my_nspinor=max(1,dtorbmag%nspinor/mpi_enreg%nproc_spinor) 145 146 ! TODO: generalize to nsppol > 1 147 isppol = 1 148 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 149 150 nband_k = dtorbmag%mband_occ 151 152 if (usepaw == 1) then ! cprj allocation 153 ncpgr = cprj(1,1)%ncpgr 154 ABI_ALLOCATE(dimlmn,(dtset%natom)) 155 call pawcprj_getdim(dimlmn,dtset%natom,nattyp_dum,dtset%ntypat,dtset%typat,pawtab,'R') 156 ABI_DATATYPE_ALLOCATE(cprj_k,(dtset%natom,dtorbmag%nspinor*dtset%mband)) 157 call pawcprj_alloc(cprj_k,ncpgr,dimlmn) 158 ABI_DATATYPE_ALLOCATE(cprj_kb,(dtset%natom,dtorbmag%nspinor*dtset%mband)) 159 call pawcprj_alloc(cprj_kb,ncpgr,dimlmn) 160 ABI_DATATYPE_ALLOCATE(cprj_kg,(dtset%natom,dtorbmag%nspinor*dtset%mband)) 161 call pawcprj_alloc(cprj_kg,ncpgr,dimlmn) 162 if (dtset%kptopt /= 3) then 163 ABI_DATATYPE_ALLOCATE(cprj_ikn,(dtset%natom,dtorbmag%nspinor*dtset%mband)) 164 ABI_DATATYPE_ALLOCATE(cprj_fkn,(dtset%natom,dtorbmag%nspinor*dtset%mband)) 165 call pawcprj_alloc(cprj_ikn,ncpgr,dimlmn) 166 call pawcprj_alloc(cprj_fkn,ncpgr,dimlmn) 167 end if 168 else 169 message = ' usepaw /= 1 but Chern number calculation requires PAW ' 170 MSG_ERROR(message) 171 end if 172 173 ABI_ALLOCATE(kk_paw,(2,dtset%mband,dtset%mband)) 174 ABI_ALLOCATE(pwind_kb,(dtset%mpw)) 175 ABI_ALLOCATE(pwind_kg,(dtset%mpw)) 176 ABI_ALLOCATE(pwind_bg,(dtset%mpw)) 177 ABI_ALLOCATE(pwnsfac_k,(4,dtset%mpw)) 178 pwnsfac_k(1,:) = 1.0_dp ! bra real 179 pwnsfac_k(2,:) = 0.0_dp ! bra imag 180 pwnsfac_k(3,:) = 1.0_dp ! ket real 181 pwnsfac_k(4,:) = 0.0_dp ! ket imag 182 183 mcg1_k = dtset%mpw*nband_k 184 ABI_ALLOCATE(cg1_k,(2,mcg1_k)) 185 ABI_ALLOCATE(sflag_k,(nband_k)) 186 ABI_ALLOCATE(smat_inv,(2,nband_k,nband_k)) 187 ABI_ALLOCATE(smat_kk,(2,nband_k,nband_k)) 188 189 ddkflag = 1 190 191 !itrs = 0 means do not invoke time reversal symmetry in smatrix.F90 192 itrs = 0 193 194 job = 1 195 shiftbd = 1 196 197 ABI_ALLOCATE(has_smat,(dtorbmag%fnkpt,dtorbmag%fnkpt)) 198 ABI_ALLOCATE(smat_all,(2,nband_k,nband_k,dtorbmag%fnkpt,dtorbmag%fnkpt)) 199 has_smat(:,:)=.FALSE. 200 201 do adir = 1, 3 202 203 IA = czero 204 IB = czero 205 206 do epsabg = 1, -1, -2 207 208 if (epsabg .EQ. 1) then 209 bdir = modulo(adir,3)+1 210 gdir = modulo(adir+1,3)+1 211 else 212 bdir = modulo(adir+1,3)+1 213 gdir = modulo(adir,3)+1 214 end if 215 216 ! loop over kpts, assuming for now kptopt 3, nsppol = 1, nspinor = 1 217 ! and no parallelism, no symmorphic symmetry elements 218 do ikpt = 1, dtorbmag%fnkpt 219 220 icprj = dtorbmag%cprjindex(ikpt,isppol) 221 222 npw_k = npwarr(ikpt) 223 icg = dtorbmag%cgindex(ikpt,dtset%nsppol) 224 225 ikg = dtorbmag%fkgindex(ikpt) 226 227 call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,icprj,ikpt,0,isppol,dtset%mband,& 228 & dtset%mkmem,dtset%natom,nband_k,nband_k,my_nspinor,dtset%nsppol,0) 229 230 do bfor = 1, 2 231 if (bfor .EQ. 1) then 232 bsigma = 1 233 else 234 bsigma = -1 235 end if 236 237 dkb(1:3) = bsigma*dtorbmag%dkvecs(1:3,bdir) 238 deltab = sqrt(DOT_PRODUCT(dkb,MATMUL(gmet,dkb))) 239 240 ikptb = dtorbmag%ikpt_dk(ikpt,bfor,bdir) 241 icprjb = dtorbmag%cprjindex(ikptb,isppol) 242 243 npw_kb = npwarr(ikptb) 244 icgb = dtorbmag%cgindex(ikptb,dtset%nsppol) 245 246 pwind_kb(1:npw_k) = pwind(ikg+1:ikg+npw_k,bfor,bdir) 247 248 call pawcprj_get(atindx1,cprj_kb,cprj,dtset%natom,1,icprjb,& 249 & ikptb,0,isppol,dtset%mband,dtset%mkmem,dtset%natom,nband_k,nband_k,& 250 & my_nspinor,dtset%nsppol,0) 251 252 if (.NOT. has_smat(ikpt,ikptb)) then 253 254 call overlap_k1k2_paw(cprj_k,cprj_kb,dkb,gprimd,kk_paw,dtorbmag%lmn2max,dtorbmag%lmn_size,dtset%mband,& 255 & dtset%natom,my_nspinor,dtset%ntypat,pawang,pawrad,pawtab,dtset%typat,xred) 256 257 sflag_k=0 258 call smatrix(cg,cg,cg1_k,ddkflag,dtm_k,icg,icgb,itrs,job,nband_k,& 259 & mcg,mcg,mcg1_k,1,dtset%mpw,nband_k,nband_k,npw_k,npw_kb,my_nspinor,& 260 & pwind_kb,pwnsfac_k,sflag_k,shiftbd,smat_inv,smat_kk,kk_paw,usepaw) 261 262 do nn = 1, nband_k 263 do n1 = 1, nband_k 264 smat_all(1,nn,n1,ikpt,ikptb) = smat_kk(1,nn,n1) 265 smat_all(2,nn,n1,ikpt,ikptb) = smat_kk(2,nn,n1) 266 smat_all(1,n1,nn,ikptb,ikpt) = smat_kk(1,nn,n1) 267 smat_all(2,n1,nn,ikptb,ikpt) = -smat_kk(2,nn,n1) 268 end do 269 end do 270 271 has_smat(ikpt,ikptb) = .TRUE. 272 has_smat(ikptb,ikpt) = .TRUE. 273 274 end if 275 276 do gfor = 1, 2 277 if (gfor .EQ. 1) then 278 gsigma = 1 279 else 280 gsigma = -1 281 end if 282 283 dkg(1:3) = gsigma*dtorbmag%dkvecs(1:3,gdir) 284 deltag = sqrt(DOT_PRODUCT(dkg,MATMUL(gmet,dkg))) 285 286 ikptg = dtorbmag%ikpt_dk(ikpt,gfor,gdir) 287 icprjg = dtorbmag%cprjindex(ikptg,isppol) 288 289 npw_kg = npwarr(ikptg) 290 icgg = dtorbmag%cgindex(ikptg,dtset%nsppol) 291 292 pwind_kg(1:npw_k) = pwind(ikg+1:ikg+npw_k,gfor,gdir) 293 294 call pawcprj_get(atindx1,cprj_kg,cprj,dtset%natom,1,icprjg,& 295 & ikptg,0,isppol,dtset%mband,dtset%mkmem,dtset%natom,nband_k,nband_k,& 296 & my_nspinor,dtset%nsppol,0) 297 298 if (.NOT. has_smat(ikpt,ikptg)) then 299 300 call overlap_k1k2_paw(cprj_k,cprj_kg,dkg,gprimd,kk_paw,dtorbmag%lmn2max,dtorbmag%lmn_size,dtset%mband,& 301 & dtset%natom,my_nspinor,dtset%ntypat,pawang,pawrad,pawtab,dtset%typat,xred) 302 303 sflag_k=0 304 call smatrix(cg,cg,cg1_k,ddkflag,dtm_k,icg,icgg,itrs,job,nband_k,& 305 & mcg,mcg,mcg1_k,1,dtset%mpw,nband_k,nband_k,npw_k,npw_kg,my_nspinor,& 306 & pwind_kg,pwnsfac_k,sflag_k,shiftbd,smat_inv,smat_kk,kk_paw,usepaw) 307 308 do nn = 1, nband_k 309 do n1 = 1, nband_k 310 smat_all(1,nn,n1,ikpt,ikptg) = smat_kk(1,nn,n1) 311 smat_all(2,nn,n1,ikpt,ikptg) = smat_kk(2,nn,n1) 312 smat_all(1,n1,nn,ikptg,ikpt) = smat_kk(1,nn,n1) 313 smat_all(2,n1,nn,ikptg,ikpt) = -smat_kk(2,nn,n1) 314 end do 315 end do 316 317 has_smat(ikpt,ikptg) = .TRUE. 318 has_smat(ikptg,ikpt) = .TRUE. 319 320 end if 321 322 dkbg = dkg - dkb 323 324 if (.NOT. has_smat(ikptb,ikptg)) then 325 326 call overlap_k1k2_paw(cprj_kb,cprj_kg,dkbg,gprimd,kk_paw,dtorbmag%lmn2max,dtorbmag%lmn_size,dtset%mband,& 327 & dtset%natom,my_nspinor,dtset%ntypat,pawang,pawrad,pawtab,dtset%typat,xred) 328 329 call mkpwind_k(dkbg,dtset,dtorbmag%fnkpt,dtorbmag%fkptns,gmet,dtorbmag%indkk_f2ibz,ikptb,ikptg,& 330 & kg,dtorbmag%kgindex,mpi_enreg,npw_kb,pwind_bg,symrec) 331 332 sflag_k=0 333 call smatrix(cg,cg,cg1_k,ddkflag,dtm_k,icgb,icgg,itrs,job,nband_k,& 334 & mcg,mcg,mcg1_k,1,dtset%mpw,nband_k,nband_k,npw_kb,npw_kg,my_nspinor,& 335 & pwind_bg,pwnsfac_k,sflag_k,shiftbd,smat_inv,smat_kk,kk_paw,usepaw) 336 337 do nn = 1, nband_k 338 do n1 = 1, nband_k 339 smat_all(1,nn,n1,ikptb,ikptg) = smat_kk(1,nn,n1) 340 smat_all(2,nn,n1,ikptb,ikptg) = smat_kk(2,nn,n1) 341 smat_all(1,n1,nn,ikptg,ikptb) = smat_kk(1,nn,n1) 342 smat_all(2,n1,nn,ikptg,ikptb) = -smat_kk(2,nn,n1) 343 end do 344 end do 345 346 has_smat(ikptb,ikptg) = .TRUE. 347 has_smat(ikptg,ikptb) = .TRUE. 348 349 end if 350 351 do nn = 1, nband_k 352 do n1 = 1, nband_k 353 354 t1A = cmplx(smat_all(1,nn,n1,ikpt,ikptb),smat_all(2,nn,n1,ikpt,ikptb)) 355 t1B = t1A 356 357 do n2 = 1, nband_k 358 359 t2A = cmplx(smat_all(1,n1,n2,ikptb,ikptg),smat_all(2,n1,n2,ikptb,ikptg)) 360 t3A = conjg(cmplx(smat_all(1,nn,n2,ikpt,ikptg),smat_all(2,nn,n2,ikpt,ikptg))) 361 362 t2B = conjg(cmplx(smat_all(1,n2,n1,ikpt,ikptb),smat_all(2,n2,n1,ikpt,ikptb))) 363 364 do n3 = 1, nband_k 365 366 t3B = cmplx(smat_all(1,n2,n3,ikpt,ikptg),smat_all(2,n2,n3,ikpt,ikptg)) 367 t4B=conjg(cmplx(smat_all(1,nn,n3,ikpt,ikptg),smat_all(2,nn,n3,ikpt,ikptg))) 368 369 tprodB = t1B*t2B*t3B*t4B 370 IB = IB - epsabg*bsigma*gsigma*tprodB/(2.0*deltab*2.0*deltag) 371 end do 372 373 tprodA = t1A*t2A*t3A 374 IA = IA + epsabg*bsigma*gsigma*tprodA/(2.0*deltab*2.0*deltag) 375 376 end do 377 end do 378 end do 379 380 end do 381 382 383 end do 384 385 end do ! end loop over fnkpt 386 end do ! end loop over epsabg 387 388 cnum(1,adir) = real(IA+IB) 389 cnum(2,adir) = aimag(IA+IB) 390 end do ! end loop over adir 391 392 cnum(1,1:3) = MATMUL(gprimd,cnum(1,1:3)) 393 cnum(2,1:3) = MATMUL(gprimd,cnum(2,1:3)) 394 dtorbmag%chern(1,1:3) = -cnum(2,1:3)/(two_pi*dtorbmag%fnkpt) 395 dtorbmag%chern(2,1:3) = cnum(1,1:3)/(two_pi*dtorbmag%fnkpt) 396 397 write(message,'(a,a,a)')ch10,'====================================================',ch10 398 call wrtout(ab_out,message,'COLL') 399 400 write(message,'(a)')' Chern number C from orbital magnetization ' 401 call wrtout(ab_out,message,'COLL') 402 write(message,'(a,a)')'----C is a real vector, given along Cartesian directions----',ch10 403 call wrtout(ab_out,message,'COLL') 404 405 do adir = 1, 3 406 write(message,'(a,i4,a,2es16.8)')' C(',adir,') : real, imag ',& 407 & dtorbmag%chern(1,adir),dtorbmag%chern(2,adir) 408 call wrtout(ab_out,message,'COLL') 409 end do 410 411 write(message,'(a,a,a)')ch10,'====================================================',ch10 412 call wrtout(ab_out,message,'COLL') 413 414 ! ! ======================================= 415 ! ! code to test orthonormality of cg_k 416 ! ! ======================================= 417 418 ! ikpt = 2 419 ! npw_k = npwarr(ikpt) 420 ! isppol = 1 421 ! nband_k = dtorbmag%mband_occ 422 ! ABI_ALLOCATE(bra,(2,npw_k*my_nspinor)) 423 ! ABI_ALLOCATE(ket,(2,npw_k*my_nspinor)) 424 ! max_err_ovlp=0.0 425 426 ! call pawcprj_get(atindx1,cprj_k,cprj,dtset%natom,1,dtorbmag%cprjindex(ikpt,isppol),ikpt,0,isppol,dtset%mband,& 427 ! & dtset%mkmem,dtset%natom,nband_k,nband_k,my_nspinor,dtset%nsppol,0) 428 ! do bband = 1, nband_k 429 ! bra_start = dtorbmag%cgindex(ikpt,dtset%nsppol)+1+(bband-1)*npw_k*my_nspinor 430 ! bra_end = bra_start + npw_k*my_nspinor - 1 431 ! bra(1:2,1:npw_k*my_nspinor) = cg(1:2,bra_start:bra_end) 432 ! do kband = 1, nband_k 433 ! ket_start = dtorbmag%cgindex(ikpt,dtset%nsppol)+1+(kband-1)*npw_k*my_nspinor 434 ! ket_end = ket_start + npw_k*my_nspinor - 1 435 ! ket(1:2,1:npw_k*my_nspinor) = cg(1:2,ket_start:ket_end) 436 437 ! tot_r = 0.0; tot_i = 0.0 438 ! do ispinor = 1, my_nspinor 439 ! ovlp_r = 0.0; ovlp_i = 0.0 440 ! spnshft = (ispinor-1)*npw_k 441 ! do ipw = 1, npw_k 442 ! spnipw = ipw + spnshft 443 ! ovlp_r = ovlp_r + bra(1,spnipw)*ket(1,spnipw)+bra(2,spnipw)*ket(2,spnipw) 444 ! ovlp_i = ovlp_i - bra(2,spnipw)*ket(1,spnipw)+bra(1,spnipw)*ket(2,spnipw) 445 ! end do ! end loop over ipw 446 ! paw_r = 0.0; paw_i = 0.0 447 ! do iatom = 1, dtset%natom 448 ! itypat = dtset%typat(iatom) 449 ! do ilmn = 1, dtorbmag%lmn_size(itypat) 450 ! do jlmn = 1, dtorbmag%lmn_size(itypat) 451 ! klmn=max(ilmn,jlmn)*(max(ilmn,jlmn)-1)/2 + min(ilmn,jlmn) 452 ! bbs = my_nspinor*(bband-1)+ispinor 453 ! kbs = my_nspinor*(kband-1)+ispinor 454 ! cpb=cmplx(cprj_k(iatom,bbs)%cp(1,ilmn),cprj_k(iatom,bbs)%cp(2,ilmn)) 455 ! cpk=cmplx(cprj_k(iatom,kbs)%cp(1,jlmn),cprj_k(iatom,kbs)%cp(2,jlmn)) 456 ! cterm = conjg(cpb)*pawtab(itypat)%sij(klmn)*cpk 457 ! paw_r = paw_r + real(cterm) 458 ! paw_i = paw_i + aimag(cterm) 459 ! end do ! end loop over jlmn 460 ! end do ! end loop over ilmn 461 ! end do ! end loop over iatom 462 ! tot_r = tot_r + ovlp_r + paw_r 463 ! tot_i = tot_i + ovlp_i + paw_i 464 ! end do ! end loop over ispinor 465 466 ! write(std_out,'(a,2i4,2es16.8)')' JWZ Debug: chern_number bband kband ovlp : ',& 467 ! & bband,kband,tot_r,tot_i 468 ! mag_ovlp = tot_r*tot_r + tot_i*tot_i 469 ! if(bband==kband) then 470 ! err_ovlp=abs(mag_ovlp-1.0) 471 ! else 472 ! err_ovlp=abs(mag_ovlp) 473 ! end if 474 ! max_err_ovlp=MAX(max_err_ovlp,err_ovlp) 475 ! end do ! end loop over kband 476 ! end do ! end loop over bband 477 ! write(std_out,'(a,i4,es16.8)')' JWZ Debug: chern_number ikpt ovlp err : ',& 478 ! & ikpt,max_err_ovlp 479 ! ABI_DEALLOCATE(bra) 480 ! ABI_DEALLOCATE(ket) 481 482 ! ! ========================================= 483 ! ! end code to test orthonormality of cg_k 484 ! ! ========================================= 485 486 if (usepaw == 1) then 487 ABI_DEALLOCATE(dimlmn) 488 call pawcprj_free(cprj_k) 489 ABI_DATATYPE_DEALLOCATE(cprj_k) 490 call pawcprj_free(cprj_kb) 491 ABI_DATATYPE_DEALLOCATE(cprj_kb) 492 call pawcprj_free(cprj_kg) 493 ABI_DATATYPE_DEALLOCATE(cprj_kg) 494 if (dtset%kptopt /= 3) then 495 call pawcprj_free(cprj_ikn) 496 call pawcprj_free(cprj_fkn) 497 ABI_DATATYPE_DEALLOCATE(cprj_ikn) 498 ABI_DATATYPE_DEALLOCATE(cprj_fkn) 499 end if 500 end if 501 502 ABI_DEALLOCATE(kk_paw) 503 ABI_DEALLOCATE(cg1_k) 504 ABI_DEALLOCATE(sflag_k) 505 ABI_DEALLOCATE(smat_inv) 506 ABI_DEALLOCATE(smat_kk) 507 ABI_DEALLOCATE(pwind_kb) 508 ABI_DEALLOCATE(pwind_kg) 509 ABI_DEALLOCATE(pwind_bg) 510 ABI_DEALLOCATE(pwnsfac_k) 511 512 ABI_DEALLOCATE(has_smat) 513 ABI_DEALLOCATE(smat_all) 514 515 end subroutine chern_number