TABLE OF CONTENTS
ABINIT/initorbmag [ 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,setsymrhoij,smpbz,symatm,timab,wrtout,xmpi_max,xmpi_sum
SOURCE
44 #if defined HAVE_CONFIG_H 45 #include "config.h" 46 #endif 47 48 #include "abi_common.h" 49 50 subroutine initorbmag(dtorbmag,dtset,gmet,gprimd,kg,mpi_enreg,npwarr,occ,& 51 & pawtab,psps,pwind,pwind_alloc,pwnsfac,& 52 & rprimd,symrec,xred) 53 54 use defs_basis 55 use defs_datatypes 56 use defs_abitypes 57 use m_orbmag 58 use m_profiling_abi 59 use m_errors 60 use m_xmpi 61 62 use m_fftcore, only : kpgsph 63 use m_pawtab, only : pawtab_type 64 use m_pawcprj, only : pawcprj_alloc, pawcprj_getdim 65 66 !This section has been created automatically by the script Abilint (TD). 67 !Do not modify the following lines by hand. 68 #undef ABI_FUNC 69 #define ABI_FUNC 'initorbmag' 70 use interfaces_14_hidewrite 71 use interfaces_18_timing 72 use interfaces_32_util 73 use interfaces_41_geometry 74 use interfaces_56_recipspace 75 use interfaces_65_paw 76 !End of the abilint section 77 78 implicit none 79 80 !Arguments ------------------------------------ 81 !scalars 82 integer,intent(out) :: pwind_alloc 83 type(MPI_type),intent(inout) :: mpi_enreg 84 type(dataset_type),intent(inout) :: dtset 85 type(orbmag_type),intent(out) :: dtorbmag 86 type(pseudopotential_type),intent(in) :: psps 87 !arrays 88 integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt) 89 integer,intent(in) :: symrec(3,3,dtset%nsym) 90 integer,pointer :: pwind(:,:,:) 91 real(dp),intent(in) :: gmet(3,3),gprimd(3,3),occ(dtset%mband*dtset%nkpt*dtset%nsppol) 92 real(dp),intent(in) :: rprimd(3,3),xred(3,dtset%natom) 93 real(dp),pointer :: pwnsfac(:,:) 94 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 95 96 !Local variables------------------------------- 97 !scalars 98 integer :: brav,exchn2n3d,fnkpt_computed 99 integer :: iband,icg,icprj,idir,idum,idum1,ierr,ifor,ikg,ikg1 100 integer :: ikpt,ikpt_loc,ikpti,ikpt1,ikpt1f,ikpt1i 101 integer :: index,ipw,ipwnsfac,isign,isppol,istwf_k,isym,isym1,itrs,itypat 102 integer :: jpw,lmax,lmn2_size_max 103 integer :: mband_occ_k,me,me_g0,mkmem_,mkpt,my_nspinor,nband_k,nkptlatt,nproc,npw_k,npw_k1 104 integer :: option,spaceComm 105 real(dp) :: diffk1,diffk2,diffk3,ecut_eff,kpt_shifted1,kpt_shifted2,kpt_shifted3,rdum 106 character(len=500) :: message 107 !arrays 108 integer :: iadum(3),iadum1(3),dg(3) 109 integer,allocatable :: kg1_k(:,:) 110 real(dp) :: diffk(3),dk(3),dum33(3,3),kpt1(3),tsec(2) 111 real(dp),allocatable :: spkpt(:,:) 112 113 ! ************************************************************************* 114 115 DBG_ENTER("COLL") 116 117 call timab(1001,1,tsec) 118 call timab(1002,1,tsec) 119 120 !save the current value of nspinor 121 dtorbmag%nspinor = dtset%nspinor 122 123 !---------------------------------------------------------------------------- 124 !-------------------- Obtain k-point grid in the full BZ -------------------- 125 !---------------------------------------------------------------------------- 126 127 if(dtset%kptopt==1 .or. dtset%kptopt==2 .or. dtset%kptopt==4)then 128 ! Compute the number of k points in the G-space unit cell 129 nkptlatt=dtset%kptrlatt(1,1)*dtset%kptrlatt(2,2)*dtset%kptrlatt(3,3) & 130 & +dtset%kptrlatt(1,2)*dtset%kptrlatt(2,3)*dtset%kptrlatt(3,1) & 131 & +dtset%kptrlatt(1,3)*dtset%kptrlatt(2,1)*dtset%kptrlatt(3,2) & 132 & -dtset%kptrlatt(1,2)*dtset%kptrlatt(2,1)*dtset%kptrlatt(3,3) & 133 & -dtset%kptrlatt(1,3)*dtset%kptrlatt(2,2)*dtset%kptrlatt(3,1) & 134 & -dtset%kptrlatt(1,1)*dtset%kptrlatt(2,3)*dtset%kptrlatt(3,2) 135 136 ! Call smpbz to obtain the list of k-point in the full BZ - without symmetry reduction 137 option = 0 138 brav = 1 139 mkpt=nkptlatt*dtset%nshiftk 140 ABI_ALLOCATE(spkpt,(3,mkpt)) 141 call smpbz(1,ab_out,dtset%kptrlatt,mkpt,fnkpt_computed,dtset%nshiftk,option,dtset%shiftk,spkpt) 142 dtorbmag%fnkpt = fnkpt_computed 143 ABI_ALLOCATE(dtorbmag%fkptns,(3,dtorbmag%fnkpt)) 144 dtorbmag%fkptns(:,:)=spkpt(:,1:dtorbmag%fnkpt) 145 ABI_DEALLOCATE(spkpt) 146 else if(dtset%kptopt==3.or.dtset%kptopt==0)then 147 dtorbmag%fnkpt=dtset%nkpt 148 ABI_ALLOCATE(dtorbmag%fkptns,(3,dtorbmag%fnkpt)) 149 dtorbmag%fkptns(1:3,1:dtorbmag%fnkpt)=dtset%kpt(1:3,1:dtorbmag%fnkpt) 150 if(dtset%kptopt==0)then 151 write(message,'(10a)') ch10,& 152 & ' initorbmag : WARNING -',ch10,& 153 & ' you have defined manually the k-point grid with kptopt = 0',ch10,& 154 & ' the orbital magnetization calculation works only with a regular k-points grid,',ch10,& 155 & ' abinit doesn''t check if your grid is regular...' 156 call wrtout(std_out,message,'PERS') 157 end if 158 end if 159 160 !call listkk to get mapping from FBZ to IBZ 161 rdum=1.0d-5 ! cutoff distance to decide when two k points match 162 ABI_ALLOCATE(dtorbmag%indkk_f2ibz,(dtorbmag%fnkpt,6)) 163 164 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 165 166 !JWZ: The following may need modification in the future 167 !**** no spin-polarization doubling ; do not allow use of time reversal symmetry **** 168 169 call timab(1002,2,tsec) 170 call timab(1003,1,tsec) 171 172 call listkk(rdum,gmet,dtorbmag%indkk_f2ibz,dtset%kptns,dtorbmag%fkptns,dtset%nkpt,& 173 & dtorbmag%fnkpt,dtset%nsym,1,dtset%symafm,symrec,0,use_symrec=.True.) 174 175 call timab(1003,2,tsec) 176 call timab(1004,1,tsec) 177 178 !Construct i2fbz and f2ibz 179 ABI_ALLOCATE(dtorbmag%i2fbz,(dtset%nkpt)) 180 idum=0 181 do ikpt=1,dtorbmag%fnkpt 182 if (dtorbmag%indkk_f2ibz(ikpt,2)==1 .and. & 183 & dtorbmag%indkk_f2ibz(ikpt,6) == 0 .and. & 184 & maxval(abs(dtorbmag%indkk_f2ibz(ikpt,3:5))) == 0 ) then 185 dtorbmag%i2fbz(dtorbmag%indkk_f2ibz(ikpt,1))=ikpt 186 idum=idum+1 187 end if 188 end do 189 if (idum/=dtset%nkpt)then 190 message = ' Found wrong number of k-points in IBZ' 191 MSG_ERROR(message) 192 end if 193 194 !---------------------------------------------------------------------------- 195 !------------- Allocate PAW space as necessary ------------------------------ 196 !---------------------------------------------------------------------------- 197 198 dtorbmag%usepaw = psps%usepaw 199 dtorbmag%natom = dtset%natom 200 dtorbmag%my_natom = mpi_enreg%my_natom 201 202 ABI_ALLOCATE(dtorbmag%lmn_size,(dtset%ntypat)) 203 ABI_ALLOCATE(dtorbmag%lmn2_size,(dtset%ntypat)) 204 do itypat = 1, dtset%ntypat 205 dtorbmag%lmn_size(itypat) = pawtab(itypat)%lmn_size 206 dtorbmag%lmn2_size(itypat) = pawtab(itypat)%lmn2_size 207 end do 208 209 lmn2_size_max = psps%lmnmax*(psps%lmnmax+1)/2 210 dtorbmag%lmn2max = lmn2_size_max 211 212 ABI_ALLOCATE(dtorbmag%cprjindex,(dtset%nkpt,dtset%nsppol)) 213 dtorbmag%cprjindex(:,:) = 0 214 215 if (dtset%kptopt /= 3) then 216 ABI_ALLOCATE(dtorbmag%atom_indsym,(4,dtset%nsym,dtorbmag%natom)) 217 call symatm(dtorbmag%atom_indsym,dtorbmag%natom,dtset%nsym,symrec,dtset%tnons,tol8,dtset%typat,xred) 218 lmax = psps%mpsang - 1 219 ABI_ALLOCATE(dtorbmag%zarot,(2*lmax+1,2*lmax+1,lmax+1,dtset%nsym)) 220 call setsymrhoij(gprimd,lmax,dtset%nsym,1,rprimd,symrec,dtorbmag%zarot) 221 dtorbmag%nsym = dtset%nsym 222 dtorbmag%lmax = lmax 223 dtorbmag%lmnmax = psps%lmnmax 224 end if 225 226 ! !------------------------------------------------------------------------------ 227 ! !------------------- Compute variables related to MPI // ---------------------- 228 ! !------------------------------------------------------------------------------ 229 spaceComm=mpi_enreg%comm_cell 230 nproc=xmpi_comm_size(spaceComm) 231 me=xmpi_comm_rank(spaceComm) 232 233 if (nproc==1) then 234 dtorbmag%fmkmem = dtorbmag%fnkpt 235 dtorbmag%fmkmem_max = dtorbmag%fnkpt 236 dtorbmag%mkmem_max = dtset%nkpt 237 else 238 dtorbmag%fmkmem = 0 239 do ikpt = 1, dtorbmag%fnkpt 240 ikpti = dtorbmag%indkk_f2ibz(ikpt,1) 241 nband_k = dtset%nband(ikpti) 242 if (.not.(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,-1,me))) & 243 & dtorbmag%fmkmem = dtorbmag%fmkmem + 1 244 end do 245 ! Maximum value of mkmem and fmkmem 246 call xmpi_max(dtorbmag%fmkmem,dtorbmag%fmkmem_max,spaceComm,ierr) 247 ! I have to use the dummy variable mkmem_ because 248 ! mkmem is declared as intent(in) while the first 249 ! argument of xmpi_max must be intent(inout) 250 mkmem_ = dtset%mkmem 251 call xmpi_max(mkmem_,dtorbmag%mkmem_max,spaceComm,ierr) 252 end if 253 254 ABI_ALLOCATE(mpi_enreg%kpt_loc2fbz_sp,(0:nproc-1,1:dtorbmag%fmkmem_max*dtset%nsppol, 1:2)) 255 ABI_ALLOCATE(mpi_enreg%kpt_loc2ibz_sp,(0:nproc-1,1:dtorbmag%mkmem_max*dtset%nsppol, 1:2)) 256 ABI_ALLOCATE(mpi_enreg%kptdstrb,(nproc,6,dtorbmag%fmkmem_max*dtset%nsppol*2)) 257 ABI_ALLOCATE(mpi_enreg%mkmem,(0:nproc-1)) 258 mpi_enreg%kpt_loc2fbz_sp(:,:,:) = 0 259 mpi_enreg%kpt_loc2ibz_sp(:,:,:) = 0 260 mpi_enreg%kptdstrb(:,:,:) = 0 261 mpi_enreg%mkmem(:) = 0 262 263 pwind_alloc = dtset%mpw*dtorbmag%fmkmem_max 264 265 ABI_ALLOCATE(pwind,(pwind_alloc,2,3)) 266 ABI_ALLOCATE(pwnsfac,(2,pwind_alloc)) 267 268 ! !------------------------------------------------------------------------------ 269 ! !---------------------- Compute orbmag_type variables ------------------------- 270 ! !------------------------------------------------------------------------------ 271 272 !Initialization of orbmag_type variables 273 dtorbmag%dkvecs(:,:) = zero 274 ABI_ALLOCATE(dtorbmag%ikpt_dk,(dtorbmag%fnkpt,2,3)) 275 ABI_ALLOCATE(dtorbmag%cgindex,(dtset%nkpt,dtset%nsppol)) 276 ABI_ALLOCATE(dtorbmag%kgindex,(dtset%nkpt)) 277 ABI_ALLOCATE(dtorbmag%fkgindex,(dtorbmag%fnkpt)) 278 dtorbmag%ikpt_dk(:,:,:) = 0 279 dtorbmag%cgindex(:,:) = 0 280 dtorbmag%mband_occ = 0 281 ABI_ALLOCATE(dtorbmag%nband_occ,(dtset%nsppol)) 282 dtorbmag%kgindex(:) = 0 283 dtorbmag%fkgindex(:) = 0 284 285 !Compute spin degeneracy 286 if (dtset%nsppol == 1 .and. dtset%nspinor == 1) then 287 dtorbmag%sdeg = two 288 else if (dtset%nsppol == 2 .or. my_nspinor == 2) then 289 dtorbmag%sdeg = one 290 end if 291 292 !Compute the number of occupied bands and check that 293 !it is the same for each k-point 294 295 index = 0 296 do isppol = 1, dtset%nsppol 297 dtorbmag%nband_occ(isppol) = 0 298 do ikpt = 1, dtset%nkpt 299 300 mband_occ_k = 0 301 nband_k = dtset%nband(ikpt + (isppol - 1)*dtset%nkpt) 302 303 do iband = 1, nband_k 304 index = index + 1 305 if (abs(occ(index) - dtorbmag%sdeg) < tol8) mband_occ_k = mband_occ_k + 1 306 end do 307 308 if (ikpt > 1) then 309 if (dtorbmag%nband_occ(isppol) /= mband_occ_k) then 310 message = "The number of valence bands is not the same for every k-point of present spin channel" 311 MSG_ERROR(message) 312 end if 313 else 314 dtorbmag%mband_occ = max(dtorbmag%mband_occ, mband_occ_k) 315 dtorbmag%nband_occ(isppol) = mband_occ_k 316 end if 317 318 end do ! close loop over ikpt 319 end do ! close loop over isppol 320 321 !Compute the location of each wavefunction 322 323 icg = 0 324 icprj = 0 325 !ikg = 0 326 do isppol = 1, dtset%nsppol 327 do ikpt = 1, dtset%nkpt 328 329 nband_k = dtset%nband(ikpt + (isppol-1)*dtset%nkpt) 330 331 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) cycle 332 333 dtorbmag%cgindex(ikpt,isppol) = icg 334 npw_k = npwarr(ikpt) 335 icg = icg + npw_k*dtorbmag%nspinor*nband_k 336 337 if (psps%usepaw == 1) then 338 dtorbmag%cprjindex(ikpt,isppol) = icprj 339 icprj = icprj + dtorbmag%nspinor*nband_k 340 end if 341 342 end do 343 end do 344 345 ikg = 0 346 do ikpt = 1, dtset%nkpt 347 if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,me)).and.& 348 & (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,dtset%nsppol,me))) cycle 349 350 npw_k = npwarr(ikpt) 351 dtorbmag%kgindex(ikpt) = ikg 352 ikg = ikg + npw_k 353 end do 354 355 call timab(1004,2,tsec) 356 357 !------------------------------------------------------------------------------ 358 !---------------------- Compute dk -------------------------------------------- 359 !------------------------------------------------------------------------------ 360 361 call timab(1005,1,tsec) 362 363 do idir = 1, 3 364 365 ! Compute dk(:), the vector between a k-point and its nearest 366 ! neighbour along the direction idir 367 368 dk(:) = zero 369 dk(idir) = 1._dp ! 1 mean there is no other k-point un the direction idir 370 do ikpt = 2, dtorbmag%fnkpt 371 diffk(:) = abs(dtorbmag%fkptns(:,ikpt) - dtorbmag%fkptns(:,1)) 372 if ((diffk(1) < dk(1)+tol8).and.(diffk(2) < dk(2)+tol8).and.& 373 & (diffk(3) < dk(3)+tol8)) dk(:) = diffk(:) 374 end do 375 dtorbmag%dkvecs(:,idir) = dk(:) 376 ! DEBUG 377 ! write(std_out,*)' initorbmag : idir, dk', idir, dk 378 ! ENDDEBUG 379 380 ! For each k point, find k_prim such that k_prim= k + dk mod(G) 381 ! where G is a vector of the reciprocal lattice 382 383 do ikpt = 1, dtorbmag%fnkpt 384 385 ! First k+dk, then k-dk 386 do isign=-1,1,2 387 kpt_shifted1=dtorbmag%fkptns(1,ikpt)- isign*dk(1) 388 kpt_shifted2=dtorbmag%fkptns(2,ikpt)- isign*dk(2) 389 kpt_shifted3=dtorbmag%fkptns(3,ikpt)- isign*dk(3) 390 ! Note that this is still a order fnkpt**2 algorithm. 391 ! It is possible to implement a order fnkpt algorithm, see listkk.F90. 392 do ikpt1 = 1, dtorbmag%fnkpt 393 diffk1=dtorbmag%fkptns(1,ikpt1) - kpt_shifted1 394 if(abs(diffk1-nint(diffk1))>tol8)cycle 395 diffk2=dtorbmag%fkptns(2,ikpt1) - kpt_shifted2 396 if(abs(diffk2-nint(diffk2))>tol8)cycle 397 diffk3=dtorbmag%fkptns(3,ikpt1) - kpt_shifted3 398 if(abs(diffk3-nint(diffk3))>tol8)cycle 399 dtorbmag%ikpt_dk(ikpt,(isign+3)/2,idir) = ikpt1 400 exit 401 end do ! ikpt1 402 end do ! isign 403 404 end do ! ikpt 405 406 end do ! close loop over idir 407 408 call timab(1005,2,tsec) 409 call timab(1006,1,tsec) 410 411 !------------------------------------------------------------------------------ 412 !------------ Build the array pwind that is needed to compute the ------------- 413 !------------ overlap matrices at k +- dk ------------- 414 !------------------------------------------------------------------------------ 415 416 ecut_eff = dtset%ecut*(dtset%dilatmx)**2 417 exchn2n3d = 0 ; istwf_k = 1 ; ikg1 = 0 418 pwind(:,:,:) = 0 419 pwnsfac(1,:) = 1.0_dp 420 pwnsfac(2,:) = 0.0_dp 421 ABI_ALLOCATE(kg1_k,(3,dtset%mpw)) 422 423 ipwnsfac = 0 424 425 do idir = 1, 3 426 427 dk(:) = dtorbmag%dkvecs(:,idir) 428 429 do ifor = 1, 2 430 431 if (ifor == 2) dk(:) = -1._dp*dk(:) 432 433 ! Build pwind and kgindex 434 ! NOTE: The array kgindex is important for parallel execution. 435 ! In case nsppol = 2, it may happen that a particular processor 436 ! treats k-points at different spin polarizations. 437 ! In this case, it is not possible to address the elements of 438 ! pwind correctly without making use of the kgindex array. 439 440 ikg = 0 ; ikpt_loc = 0 ; isppol = 1 441 do ikpt = 1, dtorbmag%fnkpt 442 443 ikpti = dtorbmag%indkk_f2ibz(ikpt,1) 444 nband_k = dtset%nband(ikpti) 445 ikpt1f = dtorbmag%ikpt_dk(ikpt,ifor,idir) 446 ikpt1i = dtorbmag%indkk_f2ibz(ikpt1f,1) 447 448 if ((proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,1,me)).and.& 449 & (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,dtset%nsppol,me))) cycle 450 451 ikpt_loc = ikpt_loc + 1 452 453 ! Build basis sphere of plane waves for the nearest neighbour of 454 ! the k-point (important for MPI //) 455 456 kg1_k(:,:) = 0 457 kpt1(:) = dtset%kptns(:,ikpt1i) 458 call kpgsph(ecut_eff,exchn2n3d,gmet,ikg1,ikpt,istwf_k,kg1_k,kpt1,& 459 & 1,mpi_enreg,dtset%mpw,npw_k1) 460 me_g0=mpi_enreg%me_g0 461 462 463 ! ji: fkgindex is defined here ! 464 dtorbmag%fkgindex(ikpt) = ikg 465 466 ! 467 ! Deal with symmetry transformations 468 ! 469 470 ! bra k-point k(b) and IBZ k-point kIBZ(b) related by 471 ! k(b) = alpha(b) S(b)^t kIBZ(b) + G(b) 472 ! where alpha(b), S(b) and G(b) are given by indkk_f2ibz 473 ! 474 ! For the ket k-point: 475 ! k(k) = alpha(k) S(k)^t kIBZ(k) + G(k) - GBZ(k) 476 ! where GBZ(k) takes k(k) to the BZ 477 ! 478 479 isym = dtorbmag%indkk_f2ibz(ikpt,2) 480 isym1 = dtorbmag%indkk_f2ibz(ikpt1f,2) 481 482 ! Construct transformed G vector that enters the matching condition: 483 ! alpha(k) S(k)^{t,-1} ( -G(b) - GBZ(k) + G(k) ) 484 485 dg(:) = -dtorbmag%indkk_f2ibz(ikpt,3:5) & 486 & -nint(-dtorbmag%fkptns(:,ikpt) - dk(:) - tol10 & 487 & +dtorbmag%fkptns(:,ikpt1f)) & 488 & +dtorbmag%indkk_f2ibz(ikpt1f,3:5) 489 490 iadum(:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),dg(:)) 491 492 dg(:) = iadum(:) 493 494 if ( dtorbmag%indkk_f2ibz(ikpt1f,6) == 1 ) dg(:) = -dg(:) 495 496 ! Construct S(k)^{t,-1} S(b)^{t} 497 498 dum33(:,:) = MATMUL(TRANSPOSE(dtset%symrel(:,:,isym1)),symrec(:,:,isym)) 499 500 ! Construct alpha(k) alpha(b) 501 502 if (dtorbmag%indkk_f2ibz(ikpt,6) == dtorbmag%indkk_f2ibz(ikpt1f,6)) then 503 itrs=0 504 else 505 itrs=1 506 end if 507 508 509 npw_k = npwarr(ikpti) 510 ! npw_k1 = npwarr(ikpt1i) 511 512 ! loop over bra G vectors 513 do ipw = 1, npw_k 514 515 ! NOTE: the bra G vector is taken for the sym-related IBZ k point, 516 ! not for the FBZ k point 517 iadum(:) = kg(:,dtorbmag%kgindex(ikpti) + ipw) 518 519 ! Store non-symmorphic operation phase factor exp[i2\pi \alpha G \cdot t] 520 521 if ( ipwnsfac == 0 ) then 522 rdum=0.0_dp 523 do idum=1,3 524 rdum=rdum+dble(iadum(idum))*dtset%tnons(idum,isym) 525 end do 526 rdum=two_pi*rdum 527 if ( dtorbmag%indkk_f2ibz(ikpt,6) == 1 ) rdum=-rdum 528 pwnsfac(1,ikg+ipw) = cos(rdum) 529 pwnsfac(2,ikg+ipw) = sin(rdum) 530 end if 531 532 ! to determine r.l.v. matchings, we transformed the bra vector 533 ! Rotation 534 iadum1(:)=0 535 do idum1=1,3 536 iadum1(:)=iadum1(:)+dum33(:,idum1)*iadum(idum1) 537 end do 538 iadum(:)=iadum1(:) 539 ! Time reversal 540 if (itrs==1) iadum(:)=-iadum(:) 541 ! Translation 542 iadum(:) = iadum(:) + dg(:) 543 544 do jpw = 1, npw_k1 545 iadum1(1:3) = kg1_k(1:3,jpw) 546 if ( (iadum(1) == iadum1(1)).and. & 547 & (iadum(2) == iadum1(2)).and. & 548 & (iadum(3) == iadum1(3)) ) then 549 pwind(ikg + ipw,ifor,idir) = jpw 550 ! write(std_out,'(a,2x,3i4,2x,i4)') 'Found !:',iadum1(:),jpw 551 exit 552 end if 553 end do 554 end do 555 556 ikg = ikg + npw_k 557 558 end do ! close loop over ikpt 559 560 ipwnsfac = 1 561 562 end do ! close loop over ifor 563 564 end do ! close loop over idir 565 566 567 call timab(1008,2,tsec) 568 call timab(1009,1,tsec) 569 570 !Build mpi_enreg%kptdstrb 571 !array required to communicate the WFs between cpus 572 !(MPI // over k-points) 573 if (nproc>1) then 574 do idir = 1, 3 575 do ifor = 1, 2 576 577 ikpt_loc = 0 578 do isppol = 1, dtset%nsppol 579 580 do ikpt = 1, dtorbmag%fnkpt 581 582 ikpti = dtorbmag%indkk_f2ibz(ikpt,1) 583 nband_k = dtset%nband(ikpti) 584 ikpt1f = dtorbmag%ikpt_dk(ikpt,ifor,idir) 585 ikpt1i = dtorbmag%indkk_f2ibz(ikpt1f,1) 586 587 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,isppol,me)) cycle 588 589 ikpt_loc = ikpt_loc + 1 590 mpi_enreg%kptdstrb(me + 1,ifor+2*(idir-1),ikpt_loc) = & 591 & ikpt1i + (isppol - 1)*dtset%nkpt 592 593 mpi_enreg%kptdstrb(me+1,ifor+2*(idir-1),ikpt_loc+dtorbmag%fmkmem_max*dtset%nsppol) = & 594 & ikpt1f + (isppol - 1)*dtorbmag%fnkpt 595 596 end do ! ikpt 597 end do ! isppol 598 end do ! ifor 599 end do ! idir 600 end if ! nproc>1 601 602 !build mpi_enreg%kpt_loc2fbz_sp 603 ikpt_loc = 0 604 do isppol = 1, dtset%nsppol 605 do ikpt = 1, dtorbmag%fnkpt 606 607 ikpti = dtorbmag%indkk_f2ibz(ikpt,1) 608 nband_k = dtset%nband(ikpti) 609 610 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpti,1,nband_k,isppol,me)) cycle 611 612 ikpt_loc = ikpt_loc + 1 613 614 mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc, 1) = ikpt 615 mpi_enreg%kpt_loc2fbz_sp(me, ikpt_loc, 2) = isppol 616 617 end do 618 end do 619 620 !should be temporary 621 !unassigned mpi_enreg%kpt_loc2fbz_sp are empty ; inform other cpu (there are better ways...) 622 mpi_enreg%mkmem(me) = dtset%mkmem 623 !do ii=ikpt_loc+1,dtefield%fmkmem_max 624 !mpi_enreg%kpt_loc2fbz_sp(me, ii, 1) = -1 625 !end do 626 627 call xmpi_sum(mpi_enreg%kptdstrb,spaceComm,ierr) 628 call xmpi_sum(mpi_enreg%kpt_loc2fbz_sp,spaceComm,ierr) 629 630 ABI_DEALLOCATE(kg1_k) 631 632 call timab(1009,2,tsec) 633 call timab(1001,2,tsec) 634 635 DBG_EXIT("COLL") 636 637 end subroutine initorbmag