TABLE OF CONTENTS


ABINIT/getshell [ Functions ]

[ Top ] [ 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.

PARENTS

      nonlinear

CHILDREN

      dgelss,getkgrid,wrtout,xmpi_max

SOURCE

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

ABINIT/m_getshell [ Modules ]

[ Top ] [ Modules ]

NAME

  m_getshell

FUNCTION

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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_getshell
28 
29  use defs_basis
30  use m_abicore
31  use m_xmpi
32  use m_errors
33  use m_linalg_interfaces
34 
35  use m_kpts,            only : getkgrid
36 
37  implicit none
38 
39  private