Line data Source code
1 : !!****m* ABINIT/m_getshell
2 : !! NAME
3 : !! m_getshell
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1999-2026 ABINIT group (MVeithen)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
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
35 : !!***
36 :
37 : public :: getshell
38 : !!***
39 :
40 : contains
41 : !!***
42 :
43 : !!****f* ABINIT/getshell
44 : !! NAME
45 : !! getshell
46 : !!
47 : !! FUNCTION
48 : !! For each k-point, set up the shells of first neighbours and find
49 : !! the weigths required for the finite difference expression
50 : !! of Marzari and Vanderbilt (see PRB 56, 12847 (1997) [[cite:Marzari1997]]).
51 : !!
52 : !! INPUTS
53 : !! gmet(3,3) = metric tensor of reciprocal space
54 : !! kptopt = option for the generation of k points
55 : !! kptrlatt = k-point lattice specification
56 : !! kpt2(3,nkpt2) = reduced coordinates of the k-points in the
57 : !! reduced part of the BZ (see below)
58 : !! mkmem = number of k points which can fit in memory
59 : !! nkpt2 = number of k-points in the reduced BZ
60 : !! nkpt3 = number of k-points in the full BZ
61 : !! nshiftk = number of kpoint grid shifts
62 : !! rmet(3,3) = metric tensor of real space
63 : !! rprimd(3,3) = dimensional primitive translations (bohr)
64 : !! shiftk = shift vectors for k point generation
65 : !! wtk2 = weight assigned to each k point
66 : !! comm=MPI communicator
67 : !!
68 : !! OUTPUT
69 : !! kneigh(30,nkpt2) = for each k-point in the reduced part of the BZ
70 : !! kneigh stores the index (ikpt) of the neighbouring
71 : !! k-points
72 : !! kg_neigh(30,nkpt2,3) = kg-neigh takes values of -1, 0 or 1,
73 : !! and can be non-zero only for a single k-point,
74 : !! a line of k-points or a plane of k-points.
75 : !! The vector joining the ikpt2-th k-point to its
76 : !! ineigh-th nearest neighbour is :
77 : !! dk(:)-nint(dk(:))+real(kg_neigh(ineigh,ikpt2,:))
78 : !! with dk(:)=kpt2(:,kneigh(ineigh,ikpt2))-kpt2(:,ikpt2)
79 : !! kptindex(2,nkpt3)
80 : !! kptindex(1,ikpt) = ikpt_rbz
81 : !! ikpt_rbz = index of the k-point in the reduced BZ
82 : !! ikpt = index of the k-point in the full BZ
83 : !! kptindex(2,ikpt) = 1: use time-reversal symmetry to transform the
84 : !! wavefunction at ikpt_rbz to the wavefunction at ikpt
85 : !! 0: ikpt belongs already to the reduced BZ
86 : !! (no transformation required)
87 : !! kpt3(3,nkpt3) = reduced coordinates of the k-points in the full BZ
88 : !! mvwtk(30,nkpt2) = weights required to evaluate the finite difference
89 : !! formula of Marzari and Vanderbilt, computed for each
90 : !! k-point in the reduced part of the BZ
91 : !! mkmem_max = maximal number of k-points on each processor (MPI //)
92 : !! nneigh = total number of neighbours required to evaluate the finite
93 : !! difference formula
94 : !!
95 : !! COMMENTS
96 : !! The array kpt2 holds the reduced coordinates of the k-points in the
97 : !! reduced part of the BZ. For example, in case time-reversal symmetry is
98 : !! used (kptopt = 2) kpt2 samples half the BZ. Since some of the neighbours
99 : !! of these k-points may lie outside the reduced BZ, getshell also needs the
100 : !! coordinates of the k-points in the full BZ.
101 : !! The coordinates of the k-points in the full BZ are stored in kpt3.
102 : !! The weights mvwtk are computed for the k-points kpt2.
103 : !!
104 : !! In case no symmetry is used to reduce the number of k-points,
105 : !! the arrays kpt2 and kpt3 are equal.
106 : !!
107 : !! SOURCE
108 :
109 21 : subroutine getshell(gmet,kneigh,kg_neigh,kptindex,kptopt,kptrlatt,kpt2,&
110 21 : & kpt3,mkmem,mkmem_max,mvwtk,&
111 21 : & 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 42 : integer :: neigh(0:6,nkpt2),symafm_dummy(1),vacuum(3)
135 21 : 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 21 : 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 21 : mkmem_cp=mkmem
145 21 : call xmpi_max(mkmem_cp,mkmem_max,comm,ier)
146 : else
147 : mkmem_max = mkmem
148 : end if
149 :
150 63 : 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 21 : if (kptopt == 3) then
157 :
158 0 : ABI_MALLOC(wtk3,(nkpt3))
159 0 : kpt3(:,:) = kpt2(:,:)
160 0 : wtk3(:) = wtk2(:)
161 0 : do ikpt = 1,nkpt3
162 0 : kptindex(1,ikpt) = ikpt
163 0 : kptindex(2,ikpt) = 0
164 : end do
165 :
166 21 : else if (kptopt == 2) then
167 :
168 63 : ABI_MALLOC(wtk3,(nkpt3))
169 21 : ii = 5 ; kptopt_used = 3
170 21 : symafm_dummy(1) = 1
171 21 : shiftk_(:,:) = 0._dp
172 105 : shiftk_(:,1:nshiftk) = shiftk(:,1:nshiftk)
173 :
174 21 : nsym1 = 1
175 21 : ABI_MALLOC(symrel1,(3,3,nsym1))
176 21 : ABI_MALLOC(tnons1,(3,nsym1))
177 273 : symrel1(:,:,1) = 0
178 21 : symrel1(1,1,1) = 1 ; symrel1(2,2,1) = 1 ; symrel1(3,3,1) = 1
179 105 : tnons1(:,:) = 0._dp
180 21 : 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 21 : vacuum,wtk3)
186 :
187 21 : if (nkpt_computed /= nkpt3) then
188 : write(msg,'(a,a,a,a,i0,a,a,i0)')&
189 0 : ' The number of k-points in the whole BZ, nkpt_computed= ',nkpt_computed,ch10,&
190 0 : ' is not twice the number of k-points in half the BZ, nkpt3=',nkpt3
191 0 : ABI_BUG(msg)
192 : end if
193 :
194 3966 : kptindex(:,:) = 0
195 1336 : do ikpt3 = 1, nkpt3
196 :
197 190511 : flag = 1
198 190511 : 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 762044 : dk_(:) = kpt3(:,ikpt3) - kpt2(:,ikpt2)
210 762044 : dk(:) = dk_(:) - nint(dk_(:))
211 190511 : if (dk(1)*dk(1) + dk(2)*dk(2) + dk(3)*dk(3) < tol10) then
212 2632 : do ii = 1, 3
213 2632 : 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 0 : 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 762044 : dk_(:) = kpt3(:,ikpt3) + kpt2(:,ikpt2)
225 762044 : dk(:) = dk_(:) - nint(dk_(:))
226 190511 : if (dk(1)*dk(1) + dk(2)*dk(2) + dk(3)*dk(3) < tol10) then
227 2632 : do ii = 1, 3
228 2632 : 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 302 : 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 762044 : dk(:) = kpt3(:,ikpt3) - kpt2(:,ikpt2)
241 190511 : if (dk(1)*dk(1) + dk(2)*dk(2) + dk(3)*dk(3) < tol10) then
242 658 : kptindex(1,ikpt3) = ikpt2
243 658 : 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 759412 : dk(:) = kpt3(:,ikpt3) + kpt2(:,ikpt2)
252 189853 : if (dk(1)*dk(1) + dk(2)*dk(2) + dk(3)*dk(3) < tol10) then
253 657 : kptindex(1,ikpt3) = ikpt2
254 657 : kptindex(2,ikpt3) = 1 ! use time-reversal symmetry
255 : flag = 0
256 : exit
257 : end if
258 :
259 : end do ! ikpt2
260 :
261 21 : if (flag == 1) then
262 0 : write(msg,'(a,i0)')' Could not find a symmetric k-point for ikpt3= ',ikpt3
263 0 : ABI_BUG(msg)
264 : end if
265 : end do ! ikpt3
266 :
267 : else
268 0 : 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 21 : dist_ = 0._dp
282 84 : do ii = 1,3
283 84 : dist_ = dist_ + gmet(ii,ii)
284 : end do
285 21 : max_dist = 3._dp * dist_ * 2._dp
286 21 : 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 21 : & + rmet(1,2)*rmet(1,2) + rmet(2,3)*rmet(2,3) + rmet(3,1)*rmet(3,1)
291 :
292 : !Initialize shell loop
293 21 : ishell = 0
294 21 : last_dist = 0._dp
295 21 : wtkflg = 0
296 20419 : kneigh(:,:) = 0
297 61278 : kg_neigh(:,:,:) = 0
298 5285 : neigh(:,:) = 0
299 :
300 : !Loop over shells until the residuum is zero
301 115 : do while ((wtkflg == 0).and.(resdm > tol8))
302 : ! Advance shell counter
303 94 : ishell = ishell + 1
304 :
305 : ! Initialize shell radius with upper limit
306 94 : 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 2126 : do ikpt = 1,nkpt3
317 : ! !write(std_out,*)ikpt
318 : ! !write(std_out,*)kpt3(:,ikpt)
319 8128 : 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 6190 : 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 14224 : do coord2 = 0,1
337 28448 : do coord3 = 0,1
338 : ! !! if ((coord1/=0).or.(coord2/=0).or.(coord3/=0)) then
339 65024 : dist_ = 0._dp
340 65024 : dk_(:) = dk(:) - nint(dk(:))
341 16256 : dk_(1) = dk_(1) + real(coord1,dp)
342 16256 : dk_(2) = dk_(2) + real(coord2,dp)
343 16256 : dk_(3) = dk_(3) + real(coord3,dp)
344 65024 : do ii = 1,3
345 211328 : do jj = 1,3
346 195072 : 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 24384 : if ((dist_ < dist(ishell)).and.(dist_ - last_dist > tol8)) then
357 360 : 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 94 : if (max_dist-dist(ishell)<tol8) then
407 0 : write(msg,'(a,i0)')' Cannot find shell number',ishell
408 0 : ABI_BUG(msg)
409 : end if
410 :
411 1111 : 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 1111 : do ikpt2 = 1, nkpt2 ! k-points in half the BZ
417 2393 : orig = sum(neigh(0:ishell-1,ikpt2))
418 : ! !write(std_out,*)'ikpt2, orig', ikpt2,orig
419 : ! !write(std_out,*) kpt2(:,ikpt2)
420 1017 : nneigh = 0
421 384289 : do ikpt3 = 1, nkpt3 ! whole k-point grid
422 : ! !! if(border_flag==0)then
423 1533088 : dk(:) = kpt3(:,ikpt3) - kpt2(:,ikpt2)
424 : ! !! dk_(:) = dk(:) - nint(dk(:))
425 : ! !! dist_ = 0._dp
426 1534105 : do coord1 = -1,1
427 4982536 : do coord2 = -1,1
428 14947608 : do coord3 = -1,1
429 : ! !! if ((coord1/=0).or.(coord2/=0).or.(coord3/=0)) then
430 41393376 : dist_ = 0._dp
431 41393376 : dk_(:) = dk(:) - nint(dk(:))
432 10348344 : dk_(1) = dk_(1) + real(coord1,dp)
433 10348344 : dk_(2) = dk_(2) + real(coord2,dp)
434 10348344 : dk_(3) = dk_(3) + real(coord3,dp)
435 41393376 : do ii = 1,3
436 134528472 : do jj = 1,3
437 124180128 : dist_ = dist_ + dk_(ii)*gmet(ii,jj)*dk_(jj)
438 : end do
439 : end do
440 13797792 : if (abs(dist_ - dist(ishell)) < tol8) then
441 4598 : nneigh = nneigh + 1
442 4598 : kneigh(orig+nneigh,ikpt2) = ikpt3
443 4598 : kg_neigh(orig+nneigh,ikpt2,1) = coord1
444 4598 : kg_neigh(orig+nneigh,ikpt2,2) = coord2
445 4598 : 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 1111 : 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 1111 : do ikpt = 1,nkpt2
504 1111 : if (neigh(ishell,ikpt) /= nneigh) flag = 0
505 : end do
506 :
507 94 : if (flag == 0) then
508 : write(msg,'(a,i0,a,a)')&
509 0 : ' The number of points in shell number',ishell,' is not the same',&
510 0 : ' for each k-point.'
511 0 : ABI_BUG(msg)
512 : end if
513 :
514 94 : if (nneigh == 0) then
515 0 : write(msg,'(a,a,a,a)') ch10,&
516 0 : ' getshell: BUG - ',ch10,&
517 0 : ' Cannot find enough neighbor shells'
518 0 : call wrtout(unts, msg)
519 0 : wtkflg = 1
520 : end if
521 :
522 : ! Calculate the total number of neighbors
523 261 : 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 94 : ikpt = 1
536 :
537 : ! Calculate the coefficients of the linear system to be solved
538 94 : mat(:,:) = 0._dp
539 261 : do is1 = 1, ishell
540 407 : orig = sum(neigh(0:is1-1,ikpt))
541 167 : bis = orig + neigh(is1,ikpt)
542 663 : do ineigh = orig+1, bis
543 1608 : dk_(:) = kpt3(:,kneigh(ineigh,ikpt)) - kpt2(:,ikpt)
544 1608 : dk(:) = dk_(:) - nint(dk_(:))
545 1608 : dk(:) = dk(:) + real(kg_neigh(ineigh,ikpt,:),dp)
546 402 : mat(1,is1) = mat(1,is1) + dk(1)*dk(1)
547 402 : mat(2,is1) = mat(2,is1) + dk(2)*dk(2)
548 402 : mat(3,is1) = mat(3,is1) + dk(3)*dk(3)
549 402 : mat(4,is1) = mat(4,is1) + dk(1)*dk(2)
550 402 : mat(5,is1) = mat(5,is1) + dk(2)*dk(3)
551 569 : mat(6,is1) = mat(6,is1) + dk(3)*dk(1)
552 : end do
553 : end do
554 :
555 94 : rvec(1) = rmet(1,1)
556 94 : rvec(2) = rmet(2,2)
557 94 : rvec(3) = rmet(3,3)
558 94 : rvec(4) = rmet(1,2)
559 94 : rvec(5) = rmet(2,3)
560 94 : 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 94 : call dgelss(6,ishell,1,mat,6,rvec,6,sgval,tol8,irank,work,30,info)
571 :
572 94 : if( info /= 0 ) then
573 : write(msg,'(3a,i0,a)')&
574 0 : ' Singular-value decomposition of the linear system determining the',ch10,&
575 0 : ' weights failed (info).',info,ch10
576 0 : ABI_COMMENT(msg)
577 0 : wtkflg = 1
578 : end if
579 :
580 : ! Check that the system has maximum rank
581 115 : if( irank == ishell ) then
582 : ! System has full rank. Calculate the residuum
583 29 : s1 = resdm
584 29 : resdm = 0._dp
585 166 : do is1 = ishell + 1, 6
586 166 : resdm = resdm + rvec(is1) * rvec(is1)
587 : end do
588 :
589 29 : if( ishell == 6 .and. resdm > tol8 ) then
590 : write(msg,'(4a)')&
591 0 : ' Linear system determining the weights could not be solved',ch10,&
592 0 : ' This should not happen.',ch10
593 0 : ABI_COMMENT(msg)
594 : wtkflg = 1
595 : end if
596 : else
597 : ! The system is rank deficient
598 65 : 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 21 : ikpt=1
613 50 : do is1 = 1, ishell
614 66 : orig = sum(neigh(0:is1-1,ikpt))
615 29 : bis = orig + neigh(is1,ikpt)
616 176 : mvwtk(orig+1:bis,1) = rvec(is1)
617 : end do
618 658 : do ikpt = 2,nkpt2
619 4480 : mvwtk(1:nneigh,ikpt) = mvwtk(1:nneigh,1)
620 : end do ! ikpt
621 :
622 : !Report weights
623 21 : write(std_out,*) 'Neighbors(1:ishell,1) ', neigh(1:ishell,1)
624 21 : write(std_out,*) 'Weights (1:ishell) ', rvec(1:ishell)
625 21 : write(std_out,*) mvwtk(1:nneigh,1)
626 :
627 : !Check the computed weights
628 21 : if (wtkflg == 0) then
629 21 : max_err = zero
630 21 : my_tol = five * tol6
631 679 : do ikpt = 1, nkpt2
632 2653 : do ii = 1,3
633 8554 : do jj = 1,3
634 5922 : s1 = 0._dp
635 41454 : do ineigh = 1, nneigh
636 142128 : dk_(:) = kpt3(:,kneigh(ineigh,ikpt)) - kpt2(:,ikpt)
637 142128 : dk(:) = dk_(:) - nint(dk_(:))
638 142128 : dk(:) = dk(:) + real(kg_neigh(ineigh,ikpt,:),dp)
639 41454 : s1 = s1 + dk(ii)*dk(jj)*mvwtk(ineigh,ikpt)
640 : end do
641 7896 : if (abs(s1 - rmet(ii,jj)) > my_tol) then
642 0 : max_err = max(max_err, abs(s1 - rmet(ii,jj)))
643 0 : wtkflg = 1
644 : end if
645 : end do
646 : end do
647 : end do
648 :
649 21 : if (wtkflg /= 0) then
650 0 : write(msg,'(5a, 2(a, es16.8))') ch10,&
651 0 : ' getshell: BUG -',ch10,&
652 0 : ' The calculated weights do not solve the linear system for all k-points.', ch10, &
653 0 : " max_err: ", max_err, " > tolerance: ", my_tol
654 0 : 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 0 : //' sets of k-points which are not quite part of the same shell. Try rectifying angles and acell.'
663 0 : ABI_BUG(msg)
664 :
665 : else
666 :
667 21 : nshell = ishell
668 :
669 21 : write(msg,'(a,a,a,a,a,a,a,i3,a,a,f16.7)') ch10,&
670 21 : ' getshell : finite difference formula of Marzari and Vanderbilt',ch10,&
671 21 : ' (see Marzari and Vanderbilt, PRB 56, 12847 (1997), Appendix B)',& ! [[cite:Marzari1997]]
672 21 : ch10,ch10,&
673 21 : ' number of first neighbours : ', neigh(1,1),ch10,&
674 42 : ' weight : ',mvwtk(1,1)
675 21 : call wrtout(unts, msg)
676 :
677 21 : if (nshell > 1) then
678 8 : is1 = neigh(1,1) + 1
679 8 : write(msg,'(a,a,i3,a,a,f16.7)')ch10,&
680 8 : ' number of second neighbours : ', neigh(2,1),ch10,&
681 16 : ' weight : ',mvwtk(is1,1)
682 8 : call wrtout(unts, msg)
683 : end if
684 :
685 8 : if (nshell > 2) then
686 0 : is1 = sum(neigh(1:2,1)) + 1
687 0 : write(msg,'(a,a,i3,a,a,f16.7)')ch10,&
688 0 : ' number of third neighbours : ', neigh(3,1),ch10,&
689 0 : ' weight : ',mvwtk(is1,1)
690 0 : call wrtout(unts, msg)
691 : end if
692 :
693 0 : if (nshell > 3) then
694 0 : is1 = sum(neigh(1:3,1)) + 1
695 0 : write(msg,'(a,a,i3,a,a,f16.7)')ch10,&
696 0 : ' number of fourth neighbours : ', neigh(4,1),ch10,&
697 0 : ' weight : ',mvwtk(is1,1)
698 0 : call wrtout(unts, msg)
699 : end if
700 :
701 0 : if (nshell > 4) then
702 0 : is1 = sum(neigh(1:4,1)) + 1
703 0 : write(msg,'(a,a,i3,a,a,f16.7)')ch10,&
704 0 : ' number of fifth neighbours : ', neigh(5,1),ch10,&
705 0 : ' weight : ',mvwtk(is1,1)
706 0 : call wrtout(unts, msg)
707 : end if
708 :
709 0 : if (nshell > 5) then
710 0 : is1 = sum(neigh(1:5,1)) + 1
711 0 : write(msg,'(a,a,i3,a,a,f16.7)')ch10,&
712 0 : ' number of sixth neighbours : ', neigh(6,1),ch10,&
713 0 : ' weight : ',mvwtk(is1,1)
714 0 : call wrtout(unts, msg)
715 : end if
716 :
717 : end if
718 :
719 21 : ABI_SFREE(tnons1)
720 21 : ABI_SFREE(symrel1)
721 :
722 21 : ABI_FREE(wtk3)
723 :
724 21 : end subroutine getshell
725 : !!***
726 :
727 : end module m_getshell
728 : !!***
|
