TABLE OF CONTENTS
ABINIT/dfptff_initberry [ Functions ]
NAME
dfptff_initberry
FUNCTION
Initialization of response calculations in finite electric field.
COPYRIGHT
Copyright (C) 2004-2018 ABINIT group (XW). 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 kg(3,mpw*mkmem) = reduced (integer) coordinates of G vecs in basis sphere kg1(3,mpw1*mkmem) = reduced (integer) coordinates of G vecs for response wfs mband = maximum number of bands mkmem = maximum number of k-points in core memory mpw = maximum number of plane waves mpw1 = maximum number of plane waves for response wavefunctions nkpt = number of k points npwarr(nkpt) = number of planewaves in basis and boundary at this k point npwar1(nkpt) = number of planewaves in basis and boundary for response wfs nsppol = 1 for unpolarized, 2 for spin-polarized occ(mband*nkpt*nsppol) = occup number for each band at each k point rprimd(3,3) = dimensional primitive vectors
OUTPUT
dtefield=variables related to response Berry-phase calculation pwindall(max(mpw,mpw1)*mkmem,8,3) = array used to compute the overlap matrices pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1> pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1> pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1> pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1> pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1> pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1> pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1> pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
NOTES
this duplicates in part initberry for the init of the dtefield - should be made into a common constructor in m_dtefield or somethin
PARENTS
dfpt_scfcv
CHILDREN
kpgio,wrtout
SOURCE
56 #if defined HAVE_CONFIG_H 57 #include "config.h" 58 #endif 59 60 #include "abi_common.h" 61 62 63 subroutine dfptff_initberry(dtefield,dtset,gmet,kg,kg1,mband,mkmem,mpi_enreg,& 64 & mpw,mpw1,nkpt,npwarr,npwar1,nsppol,occ,pwindall,rprimd) 65 66 use defs_basis 67 use defs_abitypes 68 use m_profiling_abi 69 use m_errors 70 use m_efield 71 72 use m_kg, only : kpgio 73 74 !This section has been created automatically by the script Abilint (TD). 75 !Do not modify the following lines by hand. 76 #undef ABI_FUNC 77 #define ABI_FUNC 'dfptff_initberry' 78 use interfaces_14_hidewrite 79 use interfaces_32_util 80 !End of the abilint section 81 82 implicit none 83 84 !Arguments ---------------------------------------- 85 ! scalars 86 ! arrays 87 !scalars 88 integer,intent(in) :: mband,mkmem,mpw,mpw1,nkpt,nsppol 89 type(MPI_type),intent(inout) :: mpi_enreg 90 type(dataset_type),intent(in) :: dtset 91 type(efield_type),intent(inout) :: dtefield !vz_i needs efield2 92 !arrays 93 integer,intent(in) :: kg(3,mpw*mkmem),kg1(3,mpw1*mkmem),npwar1(nkpt) 94 integer,intent(in) :: npwarr(nkpt) 95 integer,intent(out) :: pwindall(max(mpw,mpw1)*mkmem,8,3) 96 real(dp),intent(in) :: gmet(3,3),occ(mband*nkpt*nsppol),rprimd(3,3) 97 98 !Local variables ---------------------------------- 99 ! scalars 100 ! arrays 101 !scalars 102 integer :: flag,iband,icg,idir,ifor,ikg,ikg1,ikpt,ikpt1,ikpt2,ikstr 103 integer :: index,ipw,isppol,istr,iunmark,jpw,nband_k,mband_occ_k,nkstr,npw_k 104 integer :: npw_k1,orig 105 real(dp) :: ecut_eff,occ_val 106 character(len=500) :: message 107 !arrays 108 integer :: dg(3) 109 integer,allocatable :: kg_tmp(:,:),kpt_mark(:),npwarr_tmp(:),npwtot(:) 110 real(dp) :: diffk(3),dk(3) 111 real(dp),allocatable :: kpt1(:,:) 112 113 ! ************************************************************************* 114 115 !----------------------------------------------------------- 116 !---------------- initilize dtefield //--------------------- 117 !----------------------------------------------------------- 118 119 !Initialization of efield_type variables 120 dtefield%efield_dot(:) = zero 121 dtefield%dkvecs(:,:) = zero 122 dtefield%maxnstr = 0 ; dtefield%maxnkstr = 0 123 dtefield%nstr(:) = 0 ; dtefield%nkstr(:) = 0 124 ABI_ALLOCATE(dtefield%ikpt_dk,(nkpt,9,3)) 125 ABI_ALLOCATE(dtefield%cgindex,(nkpt,nsppol*2)) 126 ABI_ALLOCATE(dtefield%kgindex,(nkpt)) 127 dtefield%ikpt_dk(:,:,:) = 0 128 dtefield%cgindex(:,:) = 0 129 dtefield%mband_occ = 0 130 ABI_ALLOCATE(dtefield%nband_occ,(nsppol)) 131 dtefield%nband_occ = 0 132 pwindall(:,:,:) = 0 133 134 !Compute the number of occupied bands and check that-------- 135 !it is the same for each k-point---------------------------- 136 137 occ_val = two/(dtset%nsppol*one) 138 139 index = 0 140 do isppol = 1, nsppol 141 do ikpt = 1, nkpt 142 143 mband_occ_k = 0 144 nband_k = dtset%nband(ikpt) 145 146 do iband = 1, nband_k 147 index = index + 1 148 if (abs(occ(index) - occ_val) < tol8) mband_occ_k = mband_occ_k + 1 149 end do 150 151 if (ikpt > 1) then 152 if (dtefield%nband_occ(isppol) /= mband_occ_k) then 153 message = ' The number of valence bands is not the same for every k-point for present spin' 154 MSG_ERROR(message) 155 end if 156 else 157 dtefield%mband_occ = max(dtefield%mband_occ,mband_occ_k) 158 dtefield%nband_occ(isppol) = mband_occ_k 159 end if 160 161 end do 162 end do 163 164 !DEBUG 165 !write(std_out,*)'dfptff_initberry:nkpt',nkpt 166 !ENDDEBUG 167 168 !Compute the location of zero-order wavefunction -------------- 169 icg = 0 170 do isppol = 1, nsppol 171 do ikpt = 1, nkpt 172 173 dtefield%cgindex(ikpt,isppol) = icg 174 nband_k = dtset%nband(ikpt) 175 npw_k = npwarr(ikpt) 176 icg = icg + dtset%nspinor*npw_k*nband_k 177 178 end do 179 end do 180 181 !Compute the location of kg -------------- 182 ikg = 0 183 do ikpt = 1, nkpt 184 185 dtefield%kgindex(ikpt) = ikg 186 npw_k = npwarr(ikpt) 187 ikg = ikg + npw_k 188 189 end do 190 191 !Compute the location of first-order wavefunction ------------------- 192 icg = 0 193 do isppol = 1, nsppol 194 do ikpt = 1, nkpt 195 196 dtefield%cgindex(ikpt,isppol+nsppol) = icg 197 nband_k = dtset%nband(ikpt) 198 npw_k = npwar1(ikpt) 199 icg = icg + dtset%nspinor*npw_k*nband_k 200 201 end do 202 end do 203 204 !Compute the reciprocal lattice coordinates of the electric field----- 205 206 dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1)) 207 dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2)) 208 dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3)) 209 210 write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,& 211 & ' initberry: Reciprocal lattice coordinates of the electric field',ch10,& 212 & ' efield_dot(1:3) = ',dtefield%efield_dot(1:3),ch10 213 call wrtout(std_out,message,'COLL') 214 215 !find the related k points to every k point in full BZ----------------- 216 217 !loop over three reciprocal directions 218 do idir = 1, 3 219 220 if (dtset%rfdir(idir) == 1) then 221 222 ! Compute dk(:), the vector between a k-point and its nearest 223 ! neighbour along the direction idir 224 225 dk(:) = zero 226 dk(idir) = 1000_dp ! initialize with a big number 227 do ikpt = 2, nkpt 228 diffk(:) = abs(dtset%kptns(:,ikpt) - dtset%kptns(:,1)) 229 if ((diffk(1) < dk(1)+tol8).and.(diffk(2) < dk(2)+tol8).and.& 230 & (diffk(3) < dk(3)+tol8)) dk(:) = diffk(:) 231 end do 232 dtefield%dkvecs(:,idir) = dk(:) 233 234 ! For each k point, find k_prim such that k_prim= k + dk mod(G) 235 ! where G is a vector of the reciprocal lattice 236 237 do ikpt = 1, nkpt 238 239 ! First: k + dk 240 do ikpt1 = 1, nkpt 241 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dk(:)) 242 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 243 dtefield%ikpt_dk(ikpt,1,idir) = ikpt1 244 exit 245 end if 246 end do 247 248 ! Second: k - dk 249 do ikpt1 = 1, nkpt 250 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dk(:)) 251 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 252 dtefield%ikpt_dk(ikpt,2,idir) = ikpt1 253 exit 254 end if 255 end do 256 257 ! new 258 ! 3rd: k + (n+1)dk 259 260 do ikpt1 = 1, nkpt 261 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dk(:) - dtset%qptn(:)) 262 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 263 dtefield%ikpt_dk(ikpt,3,idir) = ikpt1 264 exit 265 end if 266 end do 267 268 ! 6th: k - (n+1)dk 269 270 do ikpt1 = 1, nkpt 271 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dk(:) + dtset%qptn(:)) 272 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 273 dtefield%ikpt_dk(ikpt,6,idir) = ikpt1 274 exit 275 end if 276 end do 277 278 ! 4th: k + (n-1)dk 279 280 do ikpt1 = 1, nkpt 281 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dk(:) - dtset%qptn(:)) 282 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 283 dtefield%ikpt_dk(ikpt,4,idir) = ikpt1 284 exit 285 end if 286 end do 287 288 ! 5th: k - (n-1)dk 289 290 do ikpt1 = 1, nkpt 291 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dk(:) + dtset%qptn(:)) 292 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 293 dtefield%ikpt_dk(ikpt,5,idir) = ikpt1 294 exit 295 end if 296 end do 297 298 ! 7th: k+n dk 299 300 do ikpt1 = 1, nkpt 301 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dtset%qptn(:)) 302 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 303 dtefield%ikpt_dk(ikpt,7,idir) = ikpt1 304 exit 305 end if 306 end do 307 308 ! 8th: k-n dk 309 310 do ikpt1 = 1, nkpt 311 diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dtset%qptn(:)) 312 if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then 313 dtefield%ikpt_dk(ikpt,8,idir) = ikpt1 314 exit 315 end if 316 end do 317 318 ! 9th: -k 319 320 do ikpt1 = 1, nkpt 321 diffk(:) = abs(dtset%kptns(:,ikpt1) + dtset%kptns(:,ikpt)) 322 if(sum(diffk(:)) < 3*tol8) then 323 dtefield%ikpt_dk(ikpt,9,idir) = ikpt1 324 exit 325 end if 326 end do 327 328 end do ! ikpt 329 330 331 332 ! Find the string length, starting from k point 1 333 ! (all strings must have the same number of points) 334 335 nkstr = 1 336 ikpt1 = 1 337 do ikpt = 1, nkpt 338 ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir) 339 if (ikpt1 == 1) exit 340 nkstr = nkstr + 1 341 end do 342 343 ! Check that the string length is a divisor of nkpt 344 if(mod(nkpt,nkstr) /= 0) then 345 write(message,'(a,i0,a,i0)')' The string length = ',nkstr,', is not a divisor of nkpt =',nkpt 346 MSG_BUG(message) 347 end if 348 dtefield%nkstr(idir) = nkstr 349 dtefield%nstr(idir) = nkpt/nkstr 350 351 end if ! dtset%rfdir(idir) == 1 352 353 write(message,'(a,i1,a,i3,a,i3)')& 354 & ' dfptff_initberry: for direction ',idir,', nkstr = ',dtefield%nkstr(idir),& 355 & ', nstr = ',dtefield%nstr(idir) 356 call wrtout(std_out,message,'COLL') 357 358 end do ! close loop over idir 359 360 dtefield%maxnstr = maxval(dtefield%nstr(:)) 361 dtefield%maxnkstr = maxval(dtefield%nkstr(:)) 362 ABI_ALLOCATE(dtefield%idxkstr,(dtefield%maxnkstr,dtefield%maxnstr,3)) 363 dtefield%idxkstr(:,:,:) = 0 364 365 366 !Build the different strings------------------------------------------ 367 368 ABI_ALLOCATE(kpt_mark,(nkpt)) 369 do idir = 1, 3 370 371 if (dtset%rfdir(idir) == 1) then 372 373 iunmark = 1 374 kpt_mark(:) = 0 375 do istr = 1, dtefield%nstr(idir) 376 377 do while(kpt_mark(iunmark) /= 0) 378 iunmark = iunmark + 1 379 end do 380 dtefield%idxkstr(1,istr,idir) = iunmark 381 kpt_mark(iunmark) = 1 382 do ikstr = 2, dtefield%nkstr(idir) 383 ikpt1 = dtefield%idxkstr(ikstr-1,istr,idir) 384 ikpt2 = dtefield%ikpt_dk(ikpt1,1,idir) 385 dtefield%idxkstr(ikstr,istr,idir) = ikpt2 386 kpt_mark(ikpt2) = 1 387 end do 388 389 end do ! istr 390 391 end if ! rfdir(idir) == 1 392 393 end do ! close loop over idir 394 395 ABI_DEALLOCATE(kpt_mark) 396 397 398 !Build the array pwindall that is needed to compute the different overlap matrices 399 !at k +- dk 400 401 ABI_ALLOCATE(kg_tmp,(3,max(mpw,mpw1)*mkmem)) 402 ABI_ALLOCATE(kpt1,(3,nkpt)) 403 ABI_ALLOCATE(npwarr_tmp,(nkpt)) 404 ABI_ALLOCATE(npwtot,(nkpt)) 405 ecut_eff = dtset%ecut*(dtset%dilatmx)**2 406 407 do idir = 1, 3 408 409 if (dtset%rfdir(idir) == 1) then 410 411 dk(:) = dtefield%dkvecs(:,idir) 412 413 do ifor = 1, 2 414 415 if (ifor == 2) dk(:) = -1_dp*dk(:) 416 417 ! Build the array kpt1 = kptns + dk 418 ! all k-poins of kpt1 must be in the same BZ as those of kptns 419 kpt1(:,:) = zero 420 do ikpt = 1, nkpt 421 ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir) 422 kpt1(:,ikpt) = dtset%kptns(:,ikpt1) 423 end do ! close loop over ikpt 424 425 ! Set up the basis sphere of plane waves at kpt1 426 kg_tmp(:,:) = 0 427 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,& 428 & kpt1,dtset%mkmem,dtset%nband,nkpt,'PERS',mpi_enreg,mpw,& 429 & npwarr_tmp,npwtot,dtset%nsppol) 430 431 ikg = 0 ; ikg1 = 0 432 do ikpt = 1, nkpt 433 434 nband_k = dtset%nband(ikpt) 435 ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir) 436 437 438 dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - & 439 & dtset%kptns(:,ikpt1)) 440 441 flag = 0; orig = 1 442 if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1 443 444 ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir) 445 npw_k = npwarr(ikpt) 446 npw_k1 = npwarr_tmp(ikpt) 447 448 do ipw = 1, npw_k 449 do jpw = orig, npw_k1 450 if ((kg(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. & 451 (kg(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. & 452 (kg(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then 453 pwindall((ikpt-1)*max(mpw,mpw1) + ipw,ifor,idir) = jpw 454 if (flag == 0) orig = jpw 455 exit 456 end if 457 end do 458 end do 459 460 ikg = ikg + npw_k 461 ikg1 = ikg1 + npw_k1 462 463 end do ! close loop over ikpt 464 465 end do ! close loop over ifor 466 467 end if ! rfdir(idir) == 1 468 469 end do ! close loop over idir 470 471 !------------------------------------------------------------------------------ 472 !<u_q|u_q> 473 !at k +- dk 474 475 do idir = 1, 3 476 477 if (dtset%rfdir(idir) == 1) then 478 479 dk(:) = dtefield%dkvecs(:,idir) 480 481 do ifor = 1, 2 482 483 if (ifor == 2) dk(:) = -1_dp*dk(:) 484 485 ! Build the array kpt1 = kptns + qptn + dk 486 ! all k-poins of kpt1 must be in the same BZ as those of kptns 487 kpt1(:,:) = zero 488 do ikpt = 1, nkpt 489 ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir) 490 kpt1(:,ikpt) = dtset%kptns(:,ikpt1)+dtset%qptn(:) 491 end do ! close loop over ikpt 492 493 ! Set UP THE BASIS SPHERE OF PLANE waves at kpt1 494 kg_tmp(:,:) = 0 495 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,& 496 & kpt1,dtset%mkmem,dtset%nband,nkpt,'PERS',mpi_enreg,mpw1,& 497 & npwarr_tmp,npwtot,dtset%nsppol) 498 499 500 ikg = 0 ; ikg1 = 0 501 do ikpt = 1, nkpt 502 503 nband_k = dtset%nband(ikpt) 504 ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir) 505 506 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,mpi_enreg%me)) cycle 507 508 dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - & 509 & dtset%kptns(:,ikpt1)) 510 511 flag = 0; orig = 1 512 if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1 513 514 ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir) 515 npw_k = npwar1(ikpt) 516 npw_k1 = npwarr_tmp(ikpt) 517 518 do ipw = 1, npw_k 519 do jpw = orig, npw_k1 520 if ((kg1(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. & 521 (kg1(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. & 522 (kg1(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then 523 pwindall((ikpt-1)*max(mpw,mpw1) +ipw,ifor+2,idir) = jpw 524 if (flag == 0) orig = jpw 525 exit 526 end if 527 end do 528 end do 529 530 ikg = ikg + npw_k 531 ikg1 = ikg1 + npw_k1 532 533 end do ! close loop over ikpt 534 535 end do ! close loop over ifor 536 537 end if ! rfdir(idir) == 1 538 539 end do ! close loop over idir 540 541 542 543 544 545 546 !--------------------------------------------------------------------------- 547 548 do idir = 1, 3 549 550 if (dtset%rfdir(idir) == 1) then 551 552 dk(:) = dtset%qptn(:) + dtefield%dkvecs(:,idir) 553 554 do ifor = 1, 2 555 556 if (ifor == 2) dk(:) = dtset%qptn(:) - dtefield%dkvecs(:,idir) 557 558 ! Build the array kpt1 = kptns + qptn + dk 559 ! all k-poins of kpt1 must be in the same BZ as those of kptns 560 kpt1(:,:) = zero 561 do ikpt = 1, nkpt 562 ikpt1 = dtefield%ikpt_dk(ikpt,ifor+2,idir) 563 kpt1(:,ikpt) = dtset%kptns(:,ikpt1) 564 end do ! close loop over ikpt 565 566 ! Set UP THE BASIS SPHERE OF PLANE waves at kpt1 567 kg_tmp(:,:) = 0 568 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,& 569 & kpt1,dtset%mkmem,dtset%nband,nkpt,'PERS',mpi_enreg,mpw,& 570 & npwarr_tmp,npwtot,dtset%nsppol) 571 572 ikg = 0 ; ikg1 = 0 573 do ikpt = 1, nkpt 574 575 nband_k = dtset%nband(ikpt) 576 ikpt1 = dtefield%ikpt_dk(ikpt,ifor+2,idir) 577 578 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,mpi_enreg%me)) cycle 579 580 dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - & 581 & dtset%kptns(:,ikpt1)) 582 583 flag = 0; orig = 1 584 if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1 585 586 ikpt1 = dtefield%ikpt_dk(ikpt,ifor+2,idir) 587 npw_k = npwar1(ikpt) 588 npw_k1 = npwarr_tmp(ikpt) 589 590 do ipw = 1, npw_k 591 do jpw = orig, npw_k1 592 if ((kg1(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. & 593 (kg1(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. & 594 (kg1(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then 595 pwindall((ikpt-1)*max(mpw,mpw1)+ipw,ifor+4,idir) = jpw 596 if (flag == 0) orig = jpw 597 exit 598 end if 599 end do 600 end do 601 ikg = ikg + npw_k 602 ikg1 = ikg1 + npw_k1 603 604 end do ! close loop over ikpt 605 606 end do ! close loop over ifor 607 608 end if ! rfdir(idir) == 1 609 610 end do ! close loop over idir=============================================================== 611 612 613 614 615 616 !Build the array pwind3 that is needed to compute the overlap matrices=============================== 617 618 do idir = 1, 3 619 620 if (dtset%rfdir(idir) == 1) then 621 622 dk(:) = - dtset%qptn(:) + dtefield%dkvecs(:,idir) 623 624 do ifor = 1, 2 625 626 if (ifor == 2) dk(:) = - dtset%qptn(:) - dtefield%dkvecs(:,idir) 627 628 ! Build the array kpt1 = kptns + qptn + dk 629 ! all k-poins of kpt1 must be in the same BZ as those of kptns 630 kpt1(:,:) = zero 631 do ikpt = 1, nkpt 632 ikpt1 = dtefield%ikpt_dk(ikpt,ifor+4,idir) 633 kpt1(:,ikpt) = dtset%kptns(:,ikpt1)+ dtset%qptn(:) 634 end do ! close loop over ikpt 635 636 ! Set UP THE BASIS SPHERE OF PLANE waves at kpt1 637 kg_tmp(:,:) = 0 638 call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,& 639 & kpt1,dtset%mkmem,dtset%nband,nkpt,'PERS',mpi_enreg,mpw1,& 640 & npwarr_tmp,npwtot,dtset%nsppol) 641 642 ikg = 0 ; ikg1 = 0 643 do ikpt = 1, nkpt 644 645 nband_k = dtset%nband(ikpt) 646 ikpt1 = dtefield%ikpt_dk(ikpt,ifor+4,idir) 647 648 if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,mpi_enreg%me)) cycle 649 650 dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - & 651 & dtset%kptns(:,ikpt1)) 652 653 flag = 0; orig = 1 654 if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1 655 656 ikpt1 = dtefield%ikpt_dk(ikpt,ifor+4,idir) 657 npw_k = npwarr(ikpt) 658 npw_k1 = npwarr_tmp(ikpt) 659 660 do ipw = 1, npw_k 661 do jpw = orig, npw_k1 662 if ((kg(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. & 663 (kg(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. & 664 (kg(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then 665 pwindall((ikpt-1)*max(mpw,mpw1) + ipw,ifor+6,idir) = jpw 666 if (flag == 0) orig = jpw 667 exit 668 end if 669 end do 670 end do 671 ikg = ikg + npw_k 672 ikg1 = ikg1 + npw_k1 673 674 end do ! close loop over ikpt 675 676 end do ! close loop over ifor 677 678 end if ! rfdir(idir) == 1 679 680 end do ! close loop over idir==================================================================== 681 682 ABI_DEALLOCATE(kg_tmp) 683 ABI_DEALLOCATE(kpt1) 684 ABI_DEALLOCATE(npwarr_tmp) 685 ABI_DEALLOCATE(npwtot) 686 687 688 end subroutine dfptff_initberry