TABLE OF CONTENTS
ABINIT/getshell [ Functions ]
NAME
getshell
FUNCTION
For each k-point, set up the shells of first neighbours and find the weigths required for the finite difference expression of Marzari and Vanderbilt (see PRB 56, 12847 (1997) [[cite:Marzari1997]]).
INPUTS
gmet(3,3) = metric tensor of reciprocal space kptopt = option for the generation of k points kptrlatt = k-point lattice specification kpt2(3,nkpt2) = reduced coordinates of the k-points in the reduced part of the BZ (see below) mkmem = number of k points which can fit in memory nkpt2 = number of k-points in the reduced BZ nkpt3 = number of k-points in the full BZ nshiftk = number of kpoint grid shifts rmet(3,3) = metric tensor of real space rprimd(3,3) = dimensional primitive translations (bohr) shiftk = shift vectors for k point generation wtk2 = weight assigned to each k point comm=MPI communicator
OUTPUT
kneigh(30,nkpt2) = for each k-point in the reduced part of the BZ
kneigh stores the index (ikpt) of the neighbouring
k-points
kg_neigh(30,nkpt2,3) = kg-neigh takes values of -1, 0 or 1,
and can be non-zero only for a single k-point,
a line of k-points or a plane of k-points.
The vector joining the ikpt2-th k-point to its
ineigh-th nearest neighbour is :
dk(:)-nint(dk(:))+real(kg_neigh(ineigh,ikpt2,:))
with dk(:)=kpt2(:,kneigh(ineigh,ikpt2))-kpt2(:,ikpt2)
kptindex(2,nkpt3)
kptindex(1,ikpt) = ikpt_rbz
ikpt_rbz = index of the k-point in the reduced BZ
ikpt = index of the k-point in the full BZ
kptindex(2,ikpt) = 1: use time-reversal symmetry to transform the
wavefunction at ikpt_rbz to the wavefunction at ikpt
0: ikpt belongs already to the reduced BZ
(no transformation required)
kpt3(3,nkpt3) = reduced coordinates of the k-points in the full BZ
mvwtk(30,nkpt2) = weights required to evaluate the finite difference
formula of Marzari and Vanderbilt, computed for each
k-point in the reduced part of the BZ
mkmem_max = maximal number of k-points on each processor (MPI //)
nneigh = total number of neighbours required to evaluate the finite
difference formula
COMMENTS
The array kpt2 holds the reduced coordinates of the k-points in the
reduced part of the BZ. For example, in case time-reversal symmetry is
used (kptopt = 2) kpt2 samples half the BZ. Since some of the neighbours
of these k-points may lie outside the reduced BZ, getshell also needs the
coordinates of the k-points in the full BZ.
The coordinates of the k-points in the full BZ are stored in kpt3.
The weights mvwtk are computed for the k-points kpt2.
In case no symmetry is used to reduce the number of k-points,
the arrays kpt2 and kpt3 are equal.
SOURCE
109 subroutine getshell(gmet,kneigh,kg_neigh,kptindex,kptopt,kptrlatt,kpt2,& 110 & kpt3,mkmem,mkmem_max,mvwtk,& 111 & nkpt2,nkpt3,nneigh,nshiftk,rmet,rprimd,shiftk,wtk2, comm) 112 113 !Arguments ------------------------------------ 114 !scalars 115 integer,intent(in) :: kptopt,mkmem,nkpt2,nkpt3,comm 116 integer,intent(inout) :: nshiftk 117 integer,intent(out) :: mkmem_max,nneigh 118 !arrays 119 integer,intent(inout) :: kptrlatt(3,3) 120 integer,intent(out) :: kneigh(30,nkpt2),kptindex(2,nkpt3),kg_neigh(30,nkpt2,3) 121 real(dp),intent(in) :: gmet(3,3),kpt2(3,nkpt2),rmet(3,3),rprimd(3,3) 122 real(dp),intent(in) :: shiftk(3,nshiftk),wtk2(nkpt2) 123 real(dp),intent(out) :: kpt3(3,nkpt3),mvwtk(30,nkpt2) 124 125 !Local variables------------------------------- 126 !scalars 127 integer :: bis,flag,ier,ii,ikpt,ikpt2,ikpt3,ineigh,info,irank,is1,ishell 128 integer :: jj,kptopt_used,mkmem_cp,nkpt_computed,nshell,nsym1,orig 129 integer :: wtkflg, coord1, coord2, coord3 130 real(dp) :: dist_,kptrlen,last_dist,max_dist,resdm,s1, max_err, my_tol 131 character(len=500) :: msg 132 !arrays 133 integer :: unts(2) 134 integer :: neigh(0:6,nkpt2),symafm_dummy(1),vacuum(3) 135 integer,allocatable :: symrel1(:,:,:) 136 real(dp) :: dist(6),dk(3),dk_(3),mat(6,6),rvec(6),sgval(6) 137 real(dp) :: shiftk_(3,MAX_NSHIFTK),work(30) 138 real(dp),allocatable :: tnons1(:,:),wtk3(:) 139 140 !************************************************************************ 141 142 !In case of MPI //: compute maximum number of k-points per processor 143 if (xmpi_paral == 1) then 144 mkmem_cp=mkmem 145 call xmpi_max(mkmem_cp,mkmem_max,comm,ier) 146 else 147 mkmem_max = mkmem 148 end if 149 150 unts = [std_out, ab_out] 151 152 !------------- In case kptopt = 2 set up the whole k-point grid ------------- 153 154 !kpt3(3,nkpt3) = reduced coordinates of k-points in the full BZ 155 156 if (kptopt == 3) then 157 158 ABI_MALLOC(wtk3,(nkpt3)) 159 kpt3(:,:) = kpt2(:,:) 160 wtk3(:) = wtk2(:) 161 do ikpt = 1,nkpt3 162 kptindex(1,ikpt) = ikpt 163 kptindex(2,ikpt) = 0 164 end do 165 166 else if (kptopt == 2) then 167 168 ABI_MALLOC(wtk3,(nkpt3)) 169 ii = 5 ; kptopt_used = 3 170 symafm_dummy(1) = 1 171 shiftk_(:,:) = 0._dp 172 shiftk_(:,1:nshiftk) = shiftk(:,1:nshiftk) 173 174 nsym1 = 1 175 ABI_MALLOC(symrel1,(3,3,nsym1)) 176 ABI_MALLOC(tnons1,(3,nsym1)) 177 symrel1(:,:,1) = 0 178 symrel1(1,1,1) = 1 ; symrel1(2,2,1) = 1 ; symrel1(3,3,1) = 1 179 tnons1(:,:) = 0._dp 180 vacuum(:) = 0 181 182 call getkgrid(0,0,ii,kpt3,kptopt_used,kptrlatt,& 183 kptrlen,nsym1,nkpt3,nkpt_computed,nshiftk,nsym1,& 184 rprimd,shiftk_,symafm_dummy,symrel1,& 185 vacuum,wtk3) 186 187 if (nkpt_computed /= nkpt3) then 188 write(msg,'(a,a,a,a,i0,a,a,i0)')& 189 ' The number of k-points in the whole BZ, nkpt_computed= ',nkpt_computed,ch10,& 190 ' is not twice the number of k-points in half the BZ, nkpt3=',nkpt3 191 ABI_BUG(msg) 192 end if 193 194 kptindex(:,:) = 0 195 do ikpt3 = 1, nkpt3 196 197 flag = 1 198 do ikpt2 = 1, nkpt2 199 200 ! In case the k-points differ only by one reciprocal lattice 201 ! vector, apply shift of one g-vector to kpt(:,ikpt3) 202 ! MJV 10/2019: this appears to be using time reversal sym, instead of the G vector... 203 ! could be equivalent to keep points inside the 1st BZ, but the code below is not consistent 204 ! with this comment 205 206 ! 207 ! here k3 = k2 + G 208 ! 209 dk_(:) = kpt3(:,ikpt3) - kpt2(:,ikpt2) 210 dk(:) = dk_(:) - nint(dk_(:)) 211 if (dk(1)*dk(1) + dk(2)*dk(2) + dk(3)*dk(3) < tol10) then 212 do ii = 1, 3 213 if ((dk(ii)*dk(ii) < tol10).and.(dk_(ii)*dk_(ii) > tol10)) then 214 ! transform k3 to -k3 215 ! TODO: I suspect this should be k3 -= G!! 216 kpt3(ii,ikpt3) = -1._dp*kpt3(ii,ikpt3) 217 end if 218 end do 219 end if 220 221 ! 222 ! here k3 = -k2 + G 223 ! 224 dk_(:) = kpt3(:,ikpt3) + kpt2(:,ikpt2) 225 dk(:) = dk_(:) - nint(dk_(:)) 226 if (dk(1)*dk(1) + dk(2)*dk(2) + dk(3)*dk(3) < tol10) then 227 do ii = 1, 3 228 if ((dk(ii)*dk(ii) < tol10).and.(dk_(ii)*dk_(ii) > tol10)) then 229 ! transform k3 to -k3 230 ! TODO: I suspect this should be k3 -= G!! 231 kpt3(ii,ikpt3) = -1._dp*kpt3(ii,ikpt3) 232 end if 233 end do 234 end if 235 236 237 ! 238 ! here k3 = k2 239 ! 240 dk(:) = kpt3(:,ikpt3) - kpt2(:,ikpt2) 241 if (dk(1)*dk(1) + dk(2)*dk(2) + dk(3)*dk(3) < tol10) then 242 kptindex(1,ikpt3) = ikpt2 243 kptindex(2,ikpt3) = 0 ! no use of time-reversal symmetry 244 flag = 0 245 exit 246 end if 247 248 ! 249 ! here k3 = -k2 250 ! 251 dk(:) = kpt3(:,ikpt3) + kpt2(:,ikpt2) 252 if (dk(1)*dk(1) + dk(2)*dk(2) + dk(3)*dk(3) < tol10) then 253 kptindex(1,ikpt3) = ikpt2 254 kptindex(2,ikpt3) = 1 ! use time-reversal symmetry 255 flag = 0 256 exit 257 end if 258 259 end do ! ikpt2 260 261 if (flag == 1) then 262 write(msg,'(a,i0)')' Could not find a symmetric k-point for ikpt3= ',ikpt3 263 ABI_BUG(msg) 264 end if 265 end do ! ikpt3 266 267 else 268 ABI_ERROR(' the only values for kptopt that are allowed are 2 and 3 ') 269 end if ! condition on kptopt 270 271 272 !--------- Compute the weights required for the Marzari-Vanderbilt --------- 273 !--------- finite difference formula --------------------------------------- 274 275 276 !Initialize distance between k-points 277 !The trace of gmet is an upper limit for its largest eigenvalue. Since the 278 !components of the distance vectors do not exceed 1, 3. * Tr[gmet] is 279 !an upper limit for the squared shell radius. 280 !we take something two times larger to make a bug checking. 281 dist_ = 0._dp 282 do ii = 1,3 283 dist_ = dist_ + gmet(ii,ii) 284 end do 285 max_dist = 3._dp * dist_ * 2._dp 286 write(std_out,*)'max_dist',max_dist 287 288 !Calculate an upper limit for the residuum 289 resdm = rmet(1,1)*rmet(1,1) + rmet(2,2)*rmet(2,2) + rmet(3,3)*rmet(3,3)& 290 & + rmet(1,2)*rmet(1,2) + rmet(2,3)*rmet(2,3) + rmet(3,1)*rmet(3,1) 291 292 !Initialize shell loop 293 ishell = 0 294 last_dist = 0._dp 295 wtkflg = 0 296 kneigh(:,:) = 0 297 kg_neigh(:,:,:) = 0 298 neigh(:,:) = 0 299 300 !Loop over shells until the residuum is zero 301 do while ((wtkflg == 0).and.(resdm > tol8)) 302 ! Advance shell counter 303 ishell = ishell + 1 304 305 ! Initialize shell radius with upper limit 306 dist(ishell) = max_dist 307 ! !! border_flag = 1 308 309 ! !write(std_out,*)'gmet' 310 ! !do ikpt=1,3 311 ! !write(std_out,*)gmet(ikpt,:) 312 ! !enddo 313 ! !write(std_out,*)kpt3(:,1) 314 315 ! Find the (squared) radius of the next shell 316 do ikpt = 1,nkpt3 317 ! !write(std_out,*)ikpt 318 ! !write(std_out,*)kpt3(:,ikpt) 319 dk(:) = kpt3(:,1) - kpt3(:,ikpt) 320 ! !!dk_(:) = dk(:) - nint(dk(:)) 321 ! !!dist_ = 0._dp 322 ! !!do ii = 1,3 323 ! !! do jj = 1,3 324 ! !! dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj) 325 ! !! end do 326 ! !!end do 327 ! !!write(std_out,*)'dist_1', dist_ 328 ! !! dist_ = 0._dp 329 ! !! do ii = 1,3 330 ! !! do jj = 1,3 331 ! !! dist_ = dist_ + dk(ii)*gmet(ii,jj)*dk(jj) 332 ! !! end do 333 ! !! end do 334 ! !!write(std_out,*)'dist_2',dist_ 335 do coord1 = 0,1 !three loop to search also on the border of the BZ, ie for the k-points (1,k2,k3) and the likes 336 do coord2 = 0,1 337 do coord3 = 0,1 338 ! !! if ((coord1/=0).or.(coord2/=0).or.(coord3/=0)) then 339 dist_ = 0._dp 340 dk_(:) = dk(:) - nint(dk(:)) 341 dk_(1) = dk_(1) + real(coord1,dp) 342 dk_(2) = dk_(2) + real(coord2,dp) 343 dk_(3) = dk_(3) + real(coord3,dp) 344 do ii = 1,3 345 do jj = 1,3 346 dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj) 347 end do 348 end do 349 ! Note : for ipkt3 = 1, coord1 = coord2 = coord3 = 0, the distance is 0 ; 350 ! but the next "if" statement is false with the tol8 criteria and the k-point 351 ! should be ignored even for ishell = 1 and last_dist= 0. 352 ! !$write(std_out,*)ikpt,coord1,coord2,coord3 353 ! !$write(std_out,*)dk_ 354 ! !$write(std_out,*)'dist_2', dist_ 355 ! !! end if 356 if ((dist_ < dist(ishell)).and.(dist_ - last_dist > tol8)) then 357 dist(ishell) = dist_ 358 end if 359 end do 360 end do 361 end do 362 363 ! !! if ((dist_ < dist(ishell)).and.(dist_ - last_dist > tol8)) then 364 ! !! dist(ishell) = dist_ 365 ! !! border_flag = 0 366 ! !! end if 367 end do 368 369 ! !! if (border_flag==1) then !we haven't found any shell in the interior of the BZ, we need to search on the border 370 ! !!write(std_out,*)ch10 371 ! !!write(std_out,*)'search on the border' 372 ! !! do ikpt = 1,nkpt3 373 ! !! dk(:) = kpt3(:,1) - kpt3(:,ikpt) 374 ! !! do coord1 = 0,1 375 ! !! do coord2 = 0,1 376 ! !! do coord3 = 0,1 377 ! !! if ((coord1/=0).or.(coord2/=0).or.(coord3/=0)) then 378 ! !! dist_ = 0._dp 379 ! !! dk_(:) = dk(:) - nint(dk(:)) 380 ! !! dk_(1) = dk_(1) + real(coord1,dp) 381 ! !! dk_(2) = dk_(2) + real(coord2,dp) 382 ! !! dk_(3) = dk_(3) + real(coord3,dp) 383 ! !! do ii = 1,3 384 ! !! do jj = 1,3 385 ! !! dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj) 386 ! !! end do 387 ! !! end do 388 ! !!write(std_out,*)ikpt,coord1,coord2,coord3 389 ! !!write(std_out,*)dk_ 390 ! !!write(std_out,*)'dist_2', dist_ 391 ! !! end if 392 ! !! if ((dist_ < dist(ishell)).and.(dist_ - last_dist > tol8)) then 393 ! !! dist(ishell) = dist_ 394 ! !! end if 395 ! !! end do 396 ! !! end do 397 ! !! end do 398 ! !! end do 399 ! !! endif 400 401 ! DEBUG 402 ! !write(std_out,*)ch10 403 ! write(std_out,*)'ishell, dist = ',ishell,dist(ishell) 404 ! ENDDEBUG 405 406 if (max_dist-dist(ishell)<tol8) then 407 write(msg,'(a,i0)')' Cannot find shell number',ishell 408 ABI_BUG(msg) 409 end if 410 411 last_dist = dist(ishell) 412 413 ! For each k-point in half the BZ get the shells of nearest neighbours. 414 ! These neighbours can be out of the zone sampled by kpt2. 415 ! !$write(std_out,*)'nkpt2', nkpt2, 'nkpt3', nkpt3 416 do ikpt2 = 1, nkpt2 ! k-points in half the BZ 417 orig = sum(neigh(0:ishell-1,ikpt2)) 418 ! !write(std_out,*)'ikpt2, orig', ikpt2,orig 419 ! !write(std_out,*) kpt2(:,ikpt2) 420 nneigh = 0 421 do ikpt3 = 1, nkpt3 ! whole k-point grid 422 ! !! if(border_flag==0)then 423 dk(:) = kpt3(:,ikpt3) - kpt2(:,ikpt2) 424 ! !! dk_(:) = dk(:) - nint(dk(:)) 425 ! !! dist_ = 0._dp 426 do coord1 = -1,1 427 do coord2 = -1,1 428 do coord3 = -1,1 429 ! !! if ((coord1/=0).or.(coord2/=0).or.(coord3/=0)) then 430 dist_ = 0._dp 431 dk_(:) = dk(:) - nint(dk(:)) 432 dk_(1) = dk_(1) + real(coord1,dp) 433 dk_(2) = dk_(2) + real(coord2,dp) 434 dk_(3) = dk_(3) + real(coord3,dp) 435 do ii = 1,3 436 do jj = 1,3 437 dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj) 438 end do 439 end do 440 if (abs(dist_ - dist(ishell)) < tol8) then 441 nneigh = nneigh + 1 442 kneigh(orig+nneigh,ikpt2) = ikpt3 443 kg_neigh(orig+nneigh,ikpt2,1) = coord1 444 kg_neigh(orig+nneigh,ikpt2,2) = coord2 445 kg_neigh(orig+nneigh,ikpt2,3) = coord3 446 end if 447 ! !! end if 448 end do 449 end do 450 end do 451 ! !write(std_out,*)'ikpt3', ikpt3 452 ! !write(std_out,*) kpt3(:,ikpt3) 453 ! write(std_out,*) kpt2(:,ikpt2) 454 ! !write(std_out,*) dk 455 ! write(std_out,*) dk_ 456 ! !! do ii = 1,3 457 ! !! do jj = 1,3 458 ! !! dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj) 459 ! !! end do 460 ! !! end do 461 ! !write(std_out,*)'dist_', dist_ 462 ! !! if (abs(dist_ - dist(ishell)) < tol8) then 463 ! !! nneigh = nneigh + 1 464 ! !! kneigh(orig+nneigh,ikpt2) = ikpt3 465 ! !! end if 466 ! !! else !search on the border 467 ! !! dk(:) = kpt3(:,ikpt3) - kpt2(:,ikpt2) 468 ! !! do coord1 = -1,1 469 ! !! do coord2 = -1,1 470 ! !! do coord3 = -1,1 471 ! !! if ((coord1/=0).or.(coord2/=0).or.(coord3/=0)) then 472 ! !! dist_ = 0._dp 473 ! !! dk_(:) = dk(:) - nint(dk(:)) 474 ! !! dk_(1) = dk_(1) + real(coord1,dp) 475 ! !! dk_(2) = dk_(2) + real(coord2,dp) 476 ! !! dk_(3) = dk_(3) + real(coord3,dp) 477 ! !! do ii = 1,3 478 ! !! do jj = 1,3 479 ! !! dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj) 480 ! !! end do 481 ! !! end do 482 ! !! if (abs(dist_ - dist(ishell)) < tol8) then 483 ! !! nneigh = nneigh + 1 484 ! !! kneigh(orig+nneigh,ikpt2) = ikpt3 485 ! !! kneigh_border(orig+nneigh,ikpt2,1) = real(coord1,dp) 486 ! !! kneigh_border(orig+nneigh,ikpt2,2) = real(coord2,dp) 487 ! !! kneigh_border(orig+nneigh,ikpt2,3) = real(coord3,dp) 488 ! !! end if 489 ! !! end if 490 ! !! end do 491 ! !! end do 492 ! !! end do 493 ! !! end if 494 end do 495 neigh(ishell,ikpt2) = nneigh 496 end do 497 498 499 ! Check if the number of points in shell number ishell 500 ! is the same for each k-point 501 502 flag = 1 503 do ikpt = 1,nkpt2 504 if (neigh(ishell,ikpt) /= nneigh) flag = 0 505 end do 506 507 if (flag == 0) then 508 write(msg,'(a,i0,a,a)')& 509 ' The number of points in shell number',ishell,' is not the same',& 510 ' for each k-point.' 511 ABI_BUG(msg) 512 end if 513 514 if (nneigh == 0) then 515 write(msg,'(a,a,a,a)') ch10,& 516 ' getshell: BUG - ',ch10,& 517 ' Cannot find enough neighbor shells' 518 call wrtout(unts, msg) 519 wtkflg = 1 520 end if 521 522 ! Calculate the total number of neighbors 523 nneigh = sum(neigh(1:ishell,1)) 524 ! DEBUG 525 !write(std_out,*)'ishell = ',ishell,'nneigh = ',nneigh 526 ! ENDDEBUG 527 528 ! Find the weights needed to compute the finite difference expression 529 ! of the ddk 530 ! ********************************************************************** 531 532 ! mvwtk(:,:) = 0._dp 533 534 ! The weights are calculated for ikpt=1. The results are copied later 535 ikpt = 1 536 537 ! Calculate the coefficients of the linear system to be solved 538 mat(:,:) = 0._dp 539 do is1 = 1, ishell 540 orig = sum(neigh(0:is1-1,ikpt)) 541 bis = orig + neigh(is1,ikpt) 542 do ineigh = orig+1, bis 543 dk_(:) = kpt3(:,kneigh(ineigh,ikpt)) - kpt2(:,ikpt) 544 dk(:) = dk_(:) - nint(dk_(:)) 545 dk(:) = dk(:) + real(kg_neigh(ineigh,ikpt,:),dp) 546 mat(1,is1) = mat(1,is1) + dk(1)*dk(1) 547 mat(2,is1) = mat(2,is1) + dk(2)*dk(2) 548 mat(3,is1) = mat(3,is1) + dk(3)*dk(3) 549 mat(4,is1) = mat(4,is1) + dk(1)*dk(2) 550 mat(5,is1) = mat(5,is1) + dk(2)*dk(3) 551 mat(6,is1) = mat(6,is1) + dk(3)*dk(1) 552 end do 553 end do 554 555 rvec(1) = rmet(1,1) 556 rvec(2) = rmet(2,2) 557 rvec(3) = rmet(3,3) 558 rvec(4) = rmet(1,2) 559 rvec(5) = rmet(2,3) 560 rvec(6) = rmet(3,1) 561 562 ! DEBUG 563 !write(std_out,*) " mat(1:6, 1:ishell) : rmet(1:6) for all 6 products dx^2... dxdy..." 564 !do ii = 1, 6 565 ! write(std_out,*) mat(ii,1:ishell), ' : ', rvec(ii) 566 !end do 567 ! ENDDEBUG 568 569 ! Solve the linear least square problem 570 call dgelss(6,ishell,1,mat,6,rvec,6,sgval,tol8,irank,work,30,info) 571 572 if( info /= 0 ) then 573 write(msg,'(3a,i0,a)')& 574 ' Singular-value decomposition of the linear system determining the',ch10,& 575 ' weights failed (info).',info,ch10 576 ABI_COMMENT(msg) 577 wtkflg = 1 578 end if 579 580 ! Check that the system has maximum rank 581 if( irank == ishell ) then 582 ! System has full rank. Calculate the residuum 583 s1 = resdm 584 resdm = 0._dp 585 do is1 = ishell + 1, 6 586 resdm = resdm + rvec(is1) * rvec(is1) 587 end do 588 589 if( ishell == 6 .and. resdm > tol8 ) then 590 write(msg,'(4a)')& 591 ' Linear system determining the weights could not be solved',ch10,& 592 ' This should not happen.',ch10 593 ABI_COMMENT(msg) 594 wtkflg = 1 595 end if 596 else 597 ! The system is rank deficient 598 ishell = ishell - 1 599 ! DEBUG 600 !write(std_out,*) 'Shell not linear independent from previous shells. Skipped.' 601 ! ENDDEBUG 602 end if 603 604 ! DEBUG 605 !write(std_out,*) "ishell, nneigh, irank, resdm ", ishell, nneigh, irank, resdm 606 ! ENDDEBUG 607 608 ! end of loop over shells 609 end do 610 611 !Copy weights 612 ikpt=1 613 do is1 = 1, ishell 614 orig = sum(neigh(0:is1-1,ikpt)) 615 bis = orig + neigh(is1,ikpt) 616 mvwtk(orig+1:bis,1) = rvec(is1) 617 end do 618 do ikpt = 2,nkpt2 619 mvwtk(1:nneigh,ikpt) = mvwtk(1:nneigh,1) 620 end do ! ikpt 621 622 !Report weights 623 write(std_out,*) 'Neighbors(1:ishell,1) ', neigh(1:ishell,1) 624 write(std_out,*) 'Weights (1:ishell) ', rvec(1:ishell) 625 write(std_out,*) mvwtk(1:nneigh,1) 626 627 !Check the computed weights 628 if (wtkflg == 0) then 629 max_err = zero 630 my_tol = five * tol6 631 do ikpt = 1, nkpt2 632 do ii = 1,3 633 do jj = 1,3 634 s1 = 0._dp 635 do ineigh = 1, nneigh 636 dk_(:) = kpt3(:,kneigh(ineigh,ikpt)) - kpt2(:,ikpt) 637 dk(:) = dk_(:) - nint(dk_(:)) 638 dk(:) = dk(:) + real(kg_neigh(ineigh,ikpt,:),dp) 639 s1 = s1 + dk(ii)*dk(jj)*mvwtk(ineigh,ikpt) 640 end do 641 if (abs(s1 - rmet(ii,jj)) > my_tol) then 642 max_err = max(max_err, abs(s1 - rmet(ii,jj))) 643 wtkflg = 1 644 end if 645 end do 646 end do 647 end do 648 649 if (wtkflg /= 0) then 650 write(msg,'(5a, 2(a, es16.8))') ch10,& 651 ' getshell: BUG -',ch10,& 652 ' The calculated weights do not solve the linear system for all k-points.', ch10, & 653 " max_err: ", max_err, " > tolerance: ", my_tol 654 call wrtout(unts, msg) 655 end if 656 end if 657 658 if (wtkflg /= 0) then 659 660 msg = ' There is a problem with the finite difference expression of the ddk '//ch10& 661 //' If you are very close to a symmetric structure, you might be confusing the algorithm with'//ch10& 662 //' sets of k-points which are not quite part of the same shell. Try rectifying angles and acell.' 663 ABI_BUG(msg) 664 665 else 666 667 nshell = ishell 668 669 write(msg,'(a,a,a,a,a,a,a,i3,a,a,f16.7)') ch10,& 670 ' getshell : finite difference formula of Marzari and Vanderbilt',ch10,& 671 ' (see Marzari and Vanderbilt, PRB 56, 12847 (1997), Appendix B)',& ! [[cite:Marzari1997]] 672 ch10,ch10,& 673 ' number of first neighbours : ', neigh(1,1),ch10,& 674 ' weight : ',mvwtk(1,1) 675 call wrtout(unts, msg) 676 677 if (nshell > 1) then 678 is1 = neigh(1,1) + 1 679 write(msg,'(a,a,i3,a,a,f16.7)')ch10,& 680 ' number of second neighbours : ', neigh(2,1),ch10,& 681 ' weight : ',mvwtk(is1,1) 682 call wrtout(unts, msg) 683 end if 684 685 if (nshell > 2) then 686 is1 = sum(neigh(1:2,1)) + 1 687 write(msg,'(a,a,i3,a,a,f16.7)')ch10,& 688 ' number of third neighbours : ', neigh(3,1),ch10,& 689 ' weight : ',mvwtk(is1,1) 690 call wrtout(unts, msg) 691 end if 692 693 if (nshell > 3) then 694 is1 = sum(neigh(1:3,1)) + 1 695 write(msg,'(a,a,i3,a,a,f16.7)')ch10,& 696 ' number of fourth neighbours : ', neigh(4,1),ch10,& 697 ' weight : ',mvwtk(is1,1) 698 call wrtout(unts, msg) 699 end if 700 701 if (nshell > 4) then 702 is1 = sum(neigh(1:4,1)) + 1 703 write(msg,'(a,a,i3,a,a,f16.7)')ch10,& 704 ' number of fifth neighbours : ', neigh(5,1),ch10,& 705 ' weight : ',mvwtk(is1,1) 706 call wrtout(unts, msg) 707 end if 708 709 if (nshell > 5) then 710 is1 = sum(neigh(1:5,1)) + 1 711 write(msg,'(a,a,i3,a,a,f16.7)')ch10,& 712 ' number of sixth neighbours : ', neigh(6,1),ch10,& 713 ' weight : ',mvwtk(is1,1) 714 call wrtout(unts, msg) 715 end if 716 717 end if 718 719 ABI_SFREE(tnons1) 720 ABI_SFREE(symrel1) 721 722 ABI_FREE(wtk3) 723 724 end subroutine getshell
ABINIT/m_getshell [ Modules ]
NAME
m_getshell
FUNCTION
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_getshell 23 24 use defs_basis 25 use m_abicore 26 use m_xmpi 27 use m_errors 28 use m_linalg_interfaces 29 30 use m_kpts, only : getkgrid 31 32 implicit none 33 34 private