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)).
COPYRIGHT
Copyright (C) 1999-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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 mpi_enreg = information about MPI parallelization 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
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.
PARENTS
nonlinear
CHILDREN
dgelss,getkgrid,wrtout,xmpi_max
SOURCE
81 #if defined HAVE_CONFIG_H 82 #include "config.h" 83 #endif 84 85 #include "abi_common.h" 86 87 88 subroutine getshell(gmet,kneigh,kg_neigh,kptindex,kptopt,kptrlatt,kpt2,& 89 & kpt3,mkmem,mkmem_max,mpi_enreg,mvwtk,& 90 & nkpt2,nkpt3,nneigh,nshiftk,rmet,rprimd,shiftk,wtk2) 91 92 use defs_basis 93 use defs_abitypes 94 use m_profiling_abi 95 use m_xmpi 96 use m_errors 97 use m_linalg_interfaces 98 99 !This section has been created automatically by the script Abilint (TD). 100 !Do not modify the following lines by hand. 101 #undef ABI_FUNC 102 #define ABI_FUNC 'getshell' 103 use interfaces_14_hidewrite 104 use interfaces_56_recipspace 105 !End of the abilint section 106 107 implicit none 108 109 !Arguments ------------------------------------ 110 !scalars 111 integer,intent(in) :: kptopt,mkmem,nkpt2,nkpt3 112 integer,intent(inout) :: nshiftk 113 integer,intent(out) :: mkmem_max,nneigh 114 type(MPI_type),intent(in) :: mpi_enreg 115 !arrays 116 integer,intent(inout) :: kptrlatt(3,3) 117 integer,intent(out) :: kneigh(30,nkpt2),kptindex(2,nkpt3),kg_neigh(30,nkpt2,3) 118 real(dp),intent(in) :: gmet(3,3),kpt2(3,nkpt2),rmet(3,3),rprimd(3,3) 119 real(dp),intent(in) :: shiftk(3,nshiftk),wtk2(nkpt2) 120 real(dp),intent(out) :: kpt3(3,nkpt3),mvwtk(30,nkpt2) 121 122 !Local variables------------------------------- 123 !scalars 124 integer :: bis,flag,ier,ii,ikpt,ikpt2,ikpt3,ineigh,info,irank,is1,ishell 125 integer :: jj,kptopt_used,mkmem_cp,nkpt_computed,nshell,nsym1,orig 126 integer :: spaceComm,wtkflg, coord1, coord2, coord3 127 real(dp) :: dist_,kptrlen,last_dist,max_dist,resdm,s1 128 character(len=500) :: message 129 !arrays 130 integer :: neigh(0:6,nkpt2),symafm_dummy(1),vacuum(3) 131 integer,allocatable :: symrel1(:,:,:) 132 real(dp) :: dist(6),dk(3),dk_(3),mat(6,6),rvec(6),sgval(6) 133 real(dp) :: shiftk_(3,210),work(30) 134 real(dp),allocatable :: tnons1(:,:),wtk3(:) 135 136 !************************************************************************ 137 138 !In case of MPI //: compute maximum number of k-points per processor 139 if (xmpi_paral == 1) then 140 spaceComm=mpi_enreg%comm_cell 141 mkmem_cp=mkmem 142 call xmpi_max(mkmem_cp,mkmem_max,spaceComm,ier) 143 else 144 mkmem_max = mkmem 145 end if 146 147 !------------- In case kptopt = 2 set up the whole k-point grid ------------- 148 149 !kpt3(3,nkpt3) = reduced coordinates of k-points in the full BZ 150 151 if (kptopt == 3) then 152 153 ABI_ALLOCATE(wtk3,(nkpt3)) 154 kpt3(:,:) = kpt2(:,:) 155 wtk3(:) = wtk2(:) 156 do ikpt = 1,nkpt3 157 kptindex(1,ikpt) = ikpt 158 kptindex(2,ikpt) = 0 159 end do 160 161 else if (kptopt == 2) then 162 163 ABI_ALLOCATE(wtk3,(nkpt3)) 164 ii = 5 ; kptopt_used = 3 165 symafm_dummy(1) = 1 166 shiftk_(:,:) = 0._dp 167 shiftk_(:,1:nshiftk) = shiftk(:,1:nshiftk) 168 169 nsym1 = 1 170 ABI_ALLOCATE(symrel1,(3,3,nsym1)) 171 ABI_ALLOCATE(tnons1,(3,nsym1)) 172 symrel1(:,:,1) = 0 173 symrel1(1,1,1) = 1 ; symrel1(2,2,1) = 1 ; symrel1(3,3,1) = 1 174 tnons1(:,:) = 0._dp 175 vacuum(:) = 0 176 177 call getkgrid(0,0,ii,kpt3,kptopt_used,kptrlatt,& 178 & kptrlen,nsym1,nkpt3,nkpt_computed,nshiftk,nsym1,& 179 & rprimd,shiftk_,symafm_dummy,symrel1,& 180 & vacuum,wtk3) 181 182 if (nkpt_computed /= nkpt3) then 183 write(message,'(a,a,a,a,i4,a,a,i4)')& 184 & ' The number of k-points in the whole BZ, nkpt_computed= ',nkpt_computed,ch10,& 185 & ' is not twice the number of k-points in half the BZ, nkpt3=',nkpt3 186 MSG_BUG(message) 187 end if 188 189 kptindex(:,:) = 0 190 do ikpt3 = 1, nkpt3 191 192 flag = 1 193 do ikpt2 = 1, nkpt2 194 195 ! In case, the k-points differ only by one reciprocal lattice 196 ! vector, apply shift of one g-vector to kpt(:,ikpt3) 197 198 dk_(:) = kpt3(:,ikpt3) - kpt2(:,ikpt2) 199 dk(:) = dk_(:) - nint(dk_(:)) 200 if (dk(1)*dk(1) + dk(2)*dk(2) + dk(3)*dk(3) < tol10) then 201 do ii = 1, 3 202 if ((dk(ii)*dk(ii) < tol10).and.(dk_(ii)*dk_(ii) > tol10)) then 203 kpt3(ii,ikpt3) = -1._dp*kpt3(ii,ikpt3) 204 end if 205 end do 206 end if 207 208 dk_(:) = kpt3(:,ikpt3) + kpt2(:,ikpt2) 209 dk(:) = dk_(:) - nint(dk_(:)) 210 if (dk(1)*dk(1) + dk(2)*dk(2) + dk(3)*dk(3) < tol10) then 211 do ii = 1, 3 212 if ((dk(ii)*dk(ii) < tol10).and.(dk_(ii)*dk_(ii) > tol10)) then 213 kpt3(ii,ikpt3) = -1._dp*kpt3(ii,ikpt3) 214 end if 215 end do 216 end if 217 218 dk(:) = kpt3(:,ikpt3) - kpt2(:,ikpt2) 219 if (dk(1)*dk(1) + dk(2)*dk(2) + dk(3)*dk(3) < tol10) then 220 kptindex(1,ikpt3) = ikpt2 221 kptindex(2,ikpt3) = 0 ! no use of time-reversal symmetry 222 flag = 0 223 exit 224 end if 225 226 dk(:) = kpt3(:,ikpt3) + kpt2(:,ikpt2) 227 if (dk(1)*dk(1) + dk(2)*dk(2) + dk(3)*dk(3) < tol10) then 228 kptindex(1,ikpt3) = ikpt2 229 kptindex(2,ikpt3) = 1 ! use time-reversal symmetry 230 flag = 0 231 exit 232 end if 233 234 end do ! ikpt2 235 236 if (flag == 1) then 237 write(message,'(a,i0)')' Could not find a symmetric k-point for ikpt3= ',ikpt3 238 MSG_BUG(message) 239 end if 240 end do ! ikpt3 241 242 else 243 message = ' the only values for kptopt that are allowed are 2 and 3 ' 244 MSG_ERROR(message) 245 end if ! condition on kptopt 246 247 248 !--------- Compute the weights required for the Marzari-Vanderbilt --------- 249 !--------- finite difference formula --------------------------------------- 250 251 252 !Initialize distance between k-points 253 !The trace of gmet is an upper limit for its largest eigenvalue. Since the 254 !components of the distance vectors do not exceed 1, 3. * Tr[gmet] is 255 !an upper limit for the squared shell radius. 256 !we take something two times larger to make a bug checking. 257 dist_ = 0._dp 258 do ii = 1,3 259 dist_ = dist_ + gmet(ii,ii) 260 end do 261 max_dist = 3._dp * dist_ * 2._dp 262 write(std_out,*)'max_dist',max_dist 263 264 !Calculate an upper limit for the residuum 265 resdm = rmet(1,1)*rmet(1,1) + rmet(2,2)*rmet(2,2) + rmet(3,3)*rmet(3,3)& 266 & + rmet(1,2)*rmet(1,2) + rmet(2,3)*rmet(2,3) + rmet(3,1)*rmet(3,1) 267 268 !Initialize shell loop 269 ishell = 0 270 last_dist = 0._dp 271 wtkflg = 0 272 kneigh(:,:) = 0 273 kg_neigh(:,:,:) = 0 274 neigh(:,:) = 0 275 276 !Loop over shells until the residuum is zero 277 do while ((wtkflg == 0).and.(resdm > tol8)) 278 ! Advance shell counter 279 ishell = ishell + 1 280 281 ! Initialize shell radius with upper limit 282 dist(ishell) = max_dist 283 ! !! border_flag = 1 284 285 ! Find the (squared) radius of the next shell 286 ! !write(std_out,*)'gmet' 287 ! !do ikpt=1,3 288 ! !write(std_out,*)gmet(ikpt,:) 289 ! !enddo 290 ! !write(std_out,*)kpt3(:,1) 291 do ikpt = 1,nkpt3 292 ! !write(std_out,*)ikpt 293 ! !write(std_out,*)kpt3(:,ikpt) 294 dk(:) = kpt3(:,1) - kpt3(:,ikpt) 295 ! !!dk_(:) = dk(:) - nint(dk(:)) 296 ! !!dist_ = 0._dp 297 ! !!do ii = 1,3 298 ! !! do jj = 1,3 299 ! !! dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj) 300 ! !! end do 301 ! !!end do 302 ! !!write(std_out,*)'dist_1', dist_ 303 ! !! dist_ = 0._dp 304 ! !! do ii = 1,3 305 ! !! do jj = 1,3 306 ! !! dist_ = dist_ + dk(ii)*gmet(ii,jj)*dk(jj) 307 ! !! end do 308 ! !! end do 309 ! !!write(std_out,*)'dist_2',dist_ 310 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 311 do coord2 = 0,1 312 do coord3 = 0,1 313 ! !! if ((coord1/=0).or.(coord2/=0).or.(coord3/=0)) then 314 dist_ = 0._dp 315 dk_(:) = dk(:) - nint(dk(:)) 316 dk_(1) = dk_(1) + real(coord1,dp) 317 dk_(2) = dk_(2) + real(coord2,dp) 318 dk_(3) = dk_(3) + real(coord3,dp) 319 do ii = 1,3 320 do jj = 1,3 321 dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj) 322 end do 323 end do 324 ! Note : for ipkt3 = 1, coord1 = coord2 = coord3 = 0, the distance is 0 ; 325 ! but the next "if" statement is false with the tol8 criteria and the k-point 326 ! should be ignored even for ishell = 1 and last_dist= 0. 327 ! !$write(std_out,*)ikpt,coord1,coord2,coord3 328 ! !$write(std_out,*)dk_ 329 ! !$write(std_out,*)'dist_2', dist_ 330 ! !! end if 331 if ((dist_ < dist(ishell)).and.(dist_ - last_dist > tol8)) then 332 dist(ishell) = dist_ 333 end if 334 end do 335 end do 336 end do 337 338 ! !! if ((dist_ < dist(ishell)).and.(dist_ - last_dist > tol8)) then 339 ! !! dist(ishell) = dist_ 340 ! !! border_flag = 0 341 ! !! end if 342 end do 343 344 ! !! if (border_flag==1) then !we haven't found any shell in the interior of the BZ, we need to search on the border 345 ! !!write(std_out,*)ch10 346 ! !!write(std_out,*)'search on the border' 347 ! !! do ikpt = 1,nkpt3 348 ! !! dk(:) = kpt3(:,1) - kpt3(:,ikpt) 349 ! !! do coord1 = 0,1 350 ! !! do coord2 = 0,1 351 ! !! do coord3 = 0,1 352 ! !! if ((coord1/=0).or.(coord2/=0).or.(coord3/=0)) then 353 ! !! dist_ = 0._dp 354 ! !! dk_(:) = dk(:) - nint(dk(:)) 355 ! !! dk_(1) = dk_(1) + real(coord1,dp) 356 ! !! dk_(2) = dk_(2) + real(coord2,dp) 357 ! !! dk_(3) = dk_(3) + real(coord3,dp) 358 ! !! do ii = 1,3 359 ! !! do jj = 1,3 360 ! !! dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj) 361 ! !! end do 362 ! !! end do 363 ! !!write(std_out,*)ikpt,coord1,coord2,coord3 364 ! !!write(std_out,*)dk_ 365 ! !!write(std_out,*)'dist_2', dist_ 366 ! !! end if 367 ! !! if ((dist_ < dist(ishell)).and.(dist_ - last_dist > tol8)) then 368 ! !! dist(ishell) = dist_ 369 ! !! end if 370 ! !! end do 371 ! !! end do 372 ! !! end do 373 ! !! end do 374 ! !! endif 375 376 ! DEBUG 377 ! !write(std_out,*)ch10 378 ! write(std_out,*)'ishell, dist = ',ishell,dist(ishell) 379 ! ENDDEBUG 380 381 if (max_dist-dist(ishell)<tol8) then 382 write(message,'(a,i0)')' Cannot find shell number',ishell 383 MSG_BUG(message) 384 end if 385 386 last_dist = dist(ishell) 387 388 ! For each k-point in halft the BZ get the shells of nearest neighbours. 389 ! These neighbours can be out of the zone sampled by kpt2. 390 ! !$write(std_out,*)'nkpt2', nkpt2, 'nkpt3', nkpt3 391 do ikpt2 = 1, nkpt2 ! k-points in half the BZ 392 orig = sum(neigh(0:ishell-1,ikpt2)) 393 ! !write(std_out,*)'ikpt2, orig', ikpt2,orig 394 ! !write(std_out,*) kpt2(:,ikpt2) 395 nneigh = 0 396 do ikpt3 = 1, nkpt3 ! whole k-point grid 397 ! !! if(border_flag==0)then 398 dk(:) = kpt3(:,ikpt3) - kpt2(:,ikpt2) 399 ! !! dk_(:) = dk(:) - nint(dk(:)) 400 ! !! dist_ = 0._dp 401 do coord1 = -1,1 402 do coord2 = -1,1 403 do coord3 = -1,1 404 ! !! if ((coord1/=0).or.(coord2/=0).or.(coord3/=0)) then 405 dist_ = 0._dp 406 dk_(:) = dk(:) - nint(dk(:)) 407 dk_(1) = dk_(1) + real(coord1,dp) 408 dk_(2) = dk_(2) + real(coord2,dp) 409 dk_(3) = dk_(3) + real(coord3,dp) 410 do ii = 1,3 411 do jj = 1,3 412 dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj) 413 end do 414 end do 415 if (abs(dist_ - dist(ishell)) < tol8) then 416 nneigh = nneigh + 1 417 kneigh(orig+nneigh,ikpt2) = ikpt3 418 kg_neigh(orig+nneigh,ikpt2,1) = coord1 419 kg_neigh(orig+nneigh,ikpt2,2) = coord2 420 kg_neigh(orig+nneigh,ikpt2,3) = coord3 421 end if 422 ! !! end if 423 end do 424 end do 425 end do 426 ! !write(std_out,*)'ikpt3', ikpt3 427 ! !write(std_out,*) kpt3(:,ikpt3) 428 ! write(std_out,*) kpt2(:,ikpt2) 429 ! !write(std_out,*) dk 430 ! write(std_out,*) dk_ 431 ! !! do ii = 1,3 432 ! !! do jj = 1,3 433 ! !! dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj) 434 ! !! end do 435 ! !! end do 436 ! !write(std_out,*)'dist_', dist_ 437 ! !! if (abs(dist_ - dist(ishell)) < tol8) then 438 ! !! nneigh = nneigh + 1 439 ! !! kneigh(orig+nneigh,ikpt2) = ikpt3 440 ! !! end if 441 ! !! else !search on the border 442 ! !! dk(:) = kpt3(:,ikpt3) - kpt2(:,ikpt2) 443 ! !! do coord1 = -1,1 444 ! !! do coord2 = -1,1 445 ! !! do coord3 = -1,1 446 ! !! if ((coord1/=0).or.(coord2/=0).or.(coord3/=0)) then 447 ! !! dist_ = 0._dp 448 ! !! dk_(:) = dk(:) - nint(dk(:)) 449 ! !! dk_(1) = dk_(1) + real(coord1,dp) 450 ! !! dk_(2) = dk_(2) + real(coord2,dp) 451 ! !! dk_(3) = dk_(3) + real(coord3,dp) 452 ! !! do ii = 1,3 453 ! !! do jj = 1,3 454 ! !! dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj) 455 ! !! end do 456 ! !! end do 457 ! !! if (abs(dist_ - dist(ishell)) < tol8) then 458 ! !! nneigh = nneigh + 1 459 ! !! kneigh(orig+nneigh,ikpt2) = ikpt3 460 ! !! kneigh_border(orig+nneigh,ikpt2,1) = real(coord1,dp) 461 ! !! kneigh_border(orig+nneigh,ikpt2,2) = real(coord2,dp) 462 ! !! kneigh_border(orig+nneigh,ikpt2,3) = real(coord3,dp) 463 ! !! end if 464 ! !! end if 465 ! !! end do 466 ! !! end do 467 ! !! end do 468 ! !! end if 469 end do 470 neigh(ishell,ikpt2) = nneigh 471 end do 472 473 474 ! Check if the number of points in shell number ishell 475 ! is the same for each k-point 476 477 flag = 1 478 do ikpt = 1,nkpt2 479 if (neigh(ishell,ikpt) /= nneigh) flag = 0 480 end do 481 482 if (flag == 0) then 483 write(message,'(a,i0,a,a)')& 484 & ' The number of points in shell number',ishell,' is not the same',& 485 & ' for each k-point.' 486 MSG_BUG(message) 487 end if 488 489 if (nneigh == 0) then 490 write(message,'(a,a,a,a)') ch10,& 491 & ' getshell: BUG - ',ch10,& 492 & ' Cannot find enough neighbor shells' 493 call wrtout(ab_out,message,'COLL') 494 call wrtout(std_out, message,'COLL') 495 wtkflg = 1 496 end if 497 498 ! Calculate the total number of neighbors 499 nneigh = sum(neigh(1:ishell,1)) 500 ! DEBUG 501 write(std_out,*)'ishell = ',ishell,'nneigh = ',nneigh 502 ! ENDDEBUG 503 504 ! Find the weights needed to compute the finite difference expression 505 ! of the ddk 506 ! ********************************************************************** 507 508 ! mvwtk(:,:) = 0._dp 509 510 ! The weights are calculated for ikpt=1. The results are copied later 511 ikpt = 1 512 513 ! Calculate the coefficients of the linear system to be solved 514 mat(:,:) = 0._dp 515 do is1 = 1, ishell 516 orig = sum(neigh(0:is1-1,ikpt)) 517 bis = orig + neigh(is1,ikpt) 518 do ineigh = orig+1, bis 519 dk_(:) = kpt3(:,kneigh(ineigh,ikpt)) - kpt2(:,ikpt) 520 dk(:) = dk_(:) - nint(dk_(:)) 521 dk(:) = dk(:) + real(kg_neigh(ineigh,ikpt,:),dp) 522 mat(1,is1) = mat(1,is1) + dk(1)*dk(1) 523 mat(2,is1) = mat(2,is1) + dk(2)*dk(2) 524 mat(3,is1) = mat(3,is1) + dk(3)*dk(3) 525 mat(4,is1) = mat(4,is1) + dk(1)*dk(2) 526 mat(5,is1) = mat(5,is1) + dk(2)*dk(3) 527 mat(6,is1) = mat(6,is1) + dk(3)*dk(1) 528 end do 529 end do 530 531 rvec(1) = rmet(1,1) 532 rvec(2) = rmet(2,2) 533 rvec(3) = rmet(3,3) 534 rvec(4) = rmet(1,2) 535 rvec(5) = rmet(2,3) 536 rvec(6) = rmet(3,1) 537 538 ! DEBUG 539 do ii = 1, 6 540 write(std_out,*)mat(ii,1:ishell), ' : ', rvec(ii) 541 end do 542 ! ENDDEBUG 543 544 ! Solve the linear least square problem 545 call dgelss(6,ishell,1,mat,6,rvec,6,sgval,tol8,irank,work,30,info) 546 547 if( info /= 0 ) then 548 write(message,'(3a,i0,a)')& 549 & ' Singular-value decomposition of the linear system determining the',ch10,& 550 & ' weights failed (info).',info,ch10 551 MSG_COMMENT(message) 552 wtkflg = 1 553 end if 554 555 ! Check that the system has maximum rank 556 if( irank == ishell ) then 557 ! System has full rank. Calculate the residuum 558 s1 = resdm 559 resdm = 0._dp 560 do is1 = ishell + 1, 6 561 resdm = resdm + rvec(is1) * rvec(is1) 562 end do 563 564 if( ishell == 6 .and. resdm > tol8 ) then 565 write(message,'(4a)')& 566 & ' Linear system determining the weights could not be solved',ch10,& 567 & ' This should not happen.',ch10 568 MSG_COMMENT(message) 569 wtkflg = 1 570 end if 571 else 572 ! The system is rank deficient 573 ishell = ishell - 1 574 ! DEBUG 575 write(std_out,*) 'Shell not linear independent from previous shells. Skipped.' 576 ! ENDDEBUG 577 end if 578 579 ! DEBUG 580 write(std_out,*) ishell, nneigh, irank, resdm 581 ! ENDDEBUG 582 583 ! end of loop over shells 584 end do 585 586 !Copy weights 587 ikpt=1 588 do is1 = 1, ishell 589 orig = sum(neigh(0:is1-1,ikpt)) 590 bis = orig + neigh(is1,ikpt) 591 mvwtk(orig+1:bis,1) = rvec(is1) 592 end do 593 do ikpt = 2,nkpt2 594 mvwtk(1:nneigh,ikpt) = mvwtk(1:nneigh,1) 595 end do ! ikpt 596 597 !Report weights 598 write(std_out,*) 'Neighbors', neigh(1:ishell,1) 599 write(std_out,*) 'Weights', rvec(1:ishell) 600 write(std_out,*) mvwtk(1:nneigh,1) 601 602 !Check the computed weights 603 if (wtkflg == 0) then 604 do ikpt = 1, nkpt2 605 do ii = 1,3 606 do jj = 1,3 607 s1 = 0._dp 608 do ineigh = 1, nneigh 609 dk_(:) = kpt3(:,kneigh(ineigh,ikpt)) - kpt2(:,ikpt) 610 dk(:) = dk_(:) - nint(dk_(:)) 611 dk(:) = dk(:) + real(kg_neigh(ineigh,ikpt,:),dp) 612 s1 = s1 + dk(ii)*dk(jj)*mvwtk(ineigh,ikpt) 613 end do 614 if (abs(s1 - rmet(ii,jj)) > tol6) wtkflg = 1 615 end do 616 end do 617 end do 618 619 if (wtkflg /= 0) then 620 write(message,'(a,a,a,a)') ch10,& 621 & ' getshell : BUG -',ch10,& 622 & ' The calculated weights do not solve the linear system for all k-points.' 623 call wrtout(ab_out,message,'COLL') 624 call wrtout(std_out, message,'COLL') 625 end if 626 end if 627 628 if (wtkflg /= 0) then 629 630 message = ' There is a problem with the finite difference expression of the ddk ' 631 MSG_BUG(message) 632 633 else 634 635 nshell = ishell 636 637 write(message,'(a,a,a,a,a,a,a,i3,a,a,f16.7)') ch10,& 638 & ' getshell : finite difference formula of Marzari and Vanderbilt',ch10,& 639 & ' (see Marzari and Vanderbilt, PRB 56, 12847 (1997), Appendix B)',& 640 & ch10,ch10,& 641 & ' number of first neighbours : ', neigh(1,1),ch10,& 642 & ' weight : ',mvwtk(1,1) 643 call wrtout(ab_out,message,'COLL') 644 call wrtout(std_out, message,'COLL') 645 646 if (nshell > 1) then 647 is1 = neigh(1,1) + 1 648 write(message,'(a,a,i3,a,a,f16.7)')ch10,& 649 & ' number of second neighbours : ', neigh(2,1),ch10,& 650 & ' weight : ',mvwtk(is1,1) 651 call wrtout(ab_out,message,'COLL') 652 call wrtout(std_out, message,'COLL') 653 end if 654 655 if (nshell > 2) then 656 is1 = sum(neigh(1:2,1)) + 1 657 write(message,'(a,a,i3,a,a,f16.7)')ch10,& 658 & ' number of third neighbours : ', neigh(3,1),ch10,& 659 & ' weight : ',mvwtk(is1,1) 660 call wrtout(ab_out,message,'COLL') 661 call wrtout(std_out, message,'COLL') 662 end if 663 664 if (nshell > 3) then 665 is1 = sum(neigh(1:3,1)) + 1 666 write(message,'(a,a,i3,a,a,f16.7)')ch10,& 667 & ' number of fourth neighbours : ', neigh(4,1),ch10,& 668 & ' weight : ',mvwtk(is1,1) 669 call wrtout(ab_out,message,'COLL') 670 call wrtout(std_out, message,'COLL') 671 end if 672 673 if (nshell > 4) then 674 is1 = sum(neigh(1:4,1)) + 1 675 write(message,'(a,a,i3,a,a,f16.7)')ch10,& 676 & ' number of fifth neighbours : ', neigh(5,1),ch10,& 677 & ' weight : ',mvwtk(is1,1) 678 call wrtout(ab_out,message,'COLL') 679 call wrtout(std_out, message,'COLL') 680 end if 681 682 if (nshell > 5) then 683 is1 = sum(neigh(1:5,1)) + 1 684 write(message,'(a,a,i3,a,a,f16.7)')ch10,& 685 & ' number of sixth neighbours : ', neigh(6,1),ch10,& 686 & ' weight : ',mvwtk(is1,1) 687 call wrtout(ab_out,message,'COLL') 688 call wrtout(std_out, message,'COLL') 689 end if 690 691 end if 692 693 694 695 !---------------------------------------------------------------------------- 696 697 if (allocated(tnons1)) then 698 ABI_DEALLOCATE(tnons1) 699 end if 700 if (allocated(symrel1)) then 701 ABI_DEALLOCATE(symrel1) 702 end if 703 704 ABI_DEALLOCATE(wtk3) 705 706 707 end subroutine getshell