Branch data 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 : : !!***
|
