Branch data Line data Source code
1 : : !!****m* ABINIT/m_gsphere
2 : : !! NAME
3 : : !! m_gsphere
4 : : !!
5 : : !! FUNCTION
6 : : !! The Gsphere data type defines the set of G-vectors
7 : : !! centered on Gamma used to describe (chi0|epsilon|W) in the GW code.
8 : : !! Note that, unlike the kg_k arrays used for wavefunctions, here the
9 : : !! G-vectors are ordered in shells (increasing length). Moreover
10 : : !! the sphere can be enlarged to take into account umklapps for which
11 : : !! one need the knowledge of several quantities at G-G0.
12 : : !!
13 : : !! COPYRIGHT
14 : : !! Copyright (C) 1999-2026 ABINIT group (MG, GMR, VO, LR, RWG, MT, XG)
15 : : !! This file is distributed under the terms of the
16 : : !! GNU General Public License, see ~abinit/COPYING
17 : : !! or http://www.gnu.org/copyleft/gpl.txt .
18 : : !!
19 : : !! SOURCE
20 : :
21 : : #if defined HAVE_CONFIG_H
22 : : #include "config.h"
23 : : #endif
24 : :
25 : : #include "abi_common.h"
26 : :
27 : : module m_gsphere
28 : :
29 : : use defs_basis
30 : : use m_abicore
31 : : use m_errors
32 : : use m_sort
33 : :
34 : : use defs_abitypes, only : MPI_type
35 : : use m_fstrings, only : sjoin, itoa
36 : : use m_numeric_tools, only : bisect
37 : : use m_geometry, only : normv
38 : : use m_crystal, only : crystal_t
39 : : use m_fftcore, only : kpgsph, kgindex, sphereboundary
40 : : use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq
41 : :
42 : : implicit none
43 : :
44 : : private
45 : :
46 : : ! Low-level procedures.
47 : : public :: merge_and_sort_kg ! Merges a set of k-centered G-spheres of cutoff ecut. Return a Gamma-centered G-spheres.
48 : : public :: table_gbig2kg ! Associate the kg_k set of G-vectors with Gamma-centered G-sphere.
49 : : public :: get_irredg ! Given a set of G vectors, find the set of G"s generating the others by symmetry.
50 : : public :: merge_kgirr ! Merge a list of irreducible G vectors (see routine for more info)
51 : : public :: setshells ! Set consistently the number of shells, the number of plane-waves, and the energy cut-off
52 : : public :: kg_map ! Compute the mapping between two lists of g-vectors.
53 : : public :: make_istwfk_table
54 : : public :: getkpgnorm ! Compute the norms of the k+G vectors
55 : : public :: symg
56 : : !!***
57 : :
58 : : !----------------------------------------------------------------------
59 : :
60 : : !!****t* m_gsphere/gsphere_t
61 : : !! NAME
62 : : !! gsphere_t
63 : : !!
64 : : !! FUNCTION
65 : : !! The gsphere_t data type contains information related to the set of G vectors
66 : : !! used during a screening or a GW calculation, as well as symmetry tables relating
67 : : !! these vectors. Presently the following quantities are stored
68 : : !!
69 : : !! 1) The reduced coordinates of the G vectors (arrays gvec)
70 : : !! 2) Tables giving the correspondence between a G-vector and its rotated image
71 : : !! through a symmetry operation in reciprocal space.
72 : : !! 3) List of the irreducible G pairs
73 : : !! 4) Tables giving, for each pair in the full reciprocal space, the corresponding
74 : : !! irreducible pair as well as the symmetry operation in reciprocal space
75 : : !!
76 : : !! Note that, unlike the GS part, the basis set does not depend on the k-point.
77 : : !!
78 : : !! NOTES
79 : : !!
80 : : !! To indicate the indices in the arrays grottb, grottbm1 we use the following notation:
81 : : !!
82 : : !! g defines the index of the reciprocal lattice vector in the array gvec
83 : : !! s indicates the index of the symmetry operation in reciprocal space
84 : : !! i can be one or two. 1 is used to indicate the identity operator
85 : : !!
86 : : !! SOURCE
87 : :
88 : : type,public :: gsphere_t
89 : :
90 : : integer :: ng
91 : : ! Total number of G vectors in the sphere taking into account umklapps
92 : : ! it defines the size of the array gvec and it accounts for possible umklapps for which
93 : : ! one has to shift the sphere.
94 : :
95 : : ! TODO: The sphere should be enlarged in gpairs_init using mg0 and (ecut|ng) as input.
96 : : ! For the time being we keep the old implementation.
97 : : !%integer :: ng_eff
98 : : ! Effective number of G vectors, i.e. the number of G in the smaller sphere without umklapps.
99 : : ! ng_eff<=ng and should be used to loop over the elements of (chi0|epsilon|W).
100 : : !
101 : : ! TODO: Add info on FFT including zero-padding algorithm.
102 : : ! table must be recalculated for each G0 and rho_tw_g should accept gpshere in input.
103 : :
104 : : integer :: nsh ! Number of shells
105 : : integer :: nsym ! The number of symmetry operations
106 : : integer :: timrev ! 2 if time-reversal is used, 1 otherwise
107 : : integer :: istwfk=1 ! Storage mode. At present time-reversal is not used.
108 : :
109 : : !integer :: mg0(3)=0
110 : : ! For each reduced direction gives the max G0 component to account for umklapp processes.
111 : :
112 : : real(dp) :: ecut
113 : : ! Cutoff energy of the sphere.
114 : :
115 : : real(dp) :: gmet(3,3)
116 : : ! Reciprocal space metric ($\textrm{bohr}^{-2}$).
117 : :
118 : : real(dp) :: gprimd(3,3)
119 : : ! Dimensional reciprocal space primitive translations (bohr^{-1})
120 : :
121 : : integer,allocatable :: g2sh(:)
122 : : ! g2sh(ng)
123 : : ! For each G, it gives the index of the shell to which it belongs.
124 : :
125 : : integer,allocatable :: gvec(:,:)
126 : : ! gvec(3,ng)
127 : : ! Reduced coordinates of G vectors.
128 : :
129 : : integer,allocatable :: g2mg(:)
130 : : ! g2mg(ng)
131 : : ! Correspondence G --> -G
132 : :
133 : : integer,allocatable :: rottb(:,:,:)
134 : : ! rottb(ng,timrev,nsym)
135 : : ! rottb(G,I,S) is the index of (SI) G in the array gvec
136 : : ! where I is either the identity or the inversion.
137 : :
138 : : integer,allocatable :: rottbm1(:,:,:)
139 : : ! rottb(ng,timrev,nsym)
140 : : ! rottbm1(G,I,S) is the index of IS{^-1} G in the array gvec
141 : :
142 : : integer,allocatable :: shlim(:)
143 : : ! shlim(nsh+1)
144 : : ! Index of the first G vector in each shell, =ng+1 for nsh+1
145 : :
146 : : real(dp),allocatable :: shlen(:)
147 : : ! shlen(nsh)
148 : : ! Radius of each shell.
149 : :
150 : : !TODO switch to dpc
151 : : complex(gwp),allocatable :: phmGt(:,:)
152 : : ! phmGt(ng,nsym)
153 : : ! Phase factor e^{-i2\pi(G.\tau)} where $\tau$ is the fractional translation associated to isym.
154 : :
155 : : complex(gwp),allocatable :: phmSGt(:,:)
156 : : ! phmSGt(ng,nsym)
157 : : ! Phase factor e^{-i2\pi(SG.\tau)} where S is one of the symmetry properties in reciprocal space.
158 : :
159 : : contains
160 : :
161 : : procedure :: init => gsph_init ! Initialize the G-sphere.
162 : : procedure :: fft_tabs => gsph_fft_tabs ! Returns useful tables for FFT (with or without padding).
163 : : procedure :: in_fftbox => gsph_in_fftbox ! Initialize the largest Gsphere contained in the FFT box.
164 : : procedure :: print => gsph_print ! Printout of basic dimensions.
165 : : procedure :: free => gsph_free ! Free memory allocated in the object.
166 : : procedure :: g_idx => gsph_g_idx ! Returns the index of G from its reduced coordinates.
167 : : procedure :: gmg_idx => gsph_gmg_idx ! Returns the index of G1-G2 from their indices
168 : : procedure :: gmg_fftidx => gsph_gmg_fftidx ! Returns the index of G1-G2 in the FFT mesh defined by ngfft.
169 : : procedure :: extend => gsph_extend ! Construct a new gsphere_t with a larger cutoff energy
170 : :
171 : : end type gsphere_t
172 : : !!***
173 : :
174 : : CONTAINS !=================================================================================
175 : : !!***
176 : :
177 : : !!****f* m_gsphere/setup_G_rotation
178 : : !! NAME
179 : : !! setup_G_rotation
180 : : !!
181 : : !! FUNCTION
182 : : !! Set up tables indicating rotation of G-vectors.
183 : : !!
184 : : !! INPUTS
185 : : !! nsym=Number of symmetry operations.
186 : : !! symrec(3,3,nsym)=Symmetry operations in reciprocal space.
187 : : !! timrev=2 if time reversal can be used, 1 otherwise.
188 : : !! npw=Number of planewaves in the sphere.
189 : : !! gvec(3,npw)=Coordinates of plane waves, supposed to be ordered in increasing modulus
190 : : !! g2sh(npw)=For each G, it gives the index of the shell to which it belongs.
191 : : !! nsh=Number of shells
192 : : !! shlim(nsh+1)=Index of the first G vector in each shell, =npw+1 for nsh+1
193 : : !!
194 : : !! OUTPUT
195 : : !! grottb (npw,2,nsym)= grottb(G,I,S) is the index of (SI) G in the array gvec.
196 : : !! grottbm1(npw,2,nsym)= index of IS^{-1} G.
197 : : !!
198 : : !! NOTES:
199 : : !! I is either the identity or the inversion (time reversal in reciprocal space).
200 : : !! S is one of the symmetry operation in reciprocal space belonging to the Space group.
201 : : !!
202 : : !! SOURCE
203 : :
204 : 1085 : subroutine setup_G_rotation(nsym,symrec,timrev,npw,gvec,g2sh,nsh,shlim,grottb,grottbm1)
205 : :
206 : : !Arguments ------------------------------------
207 : : !scalars
208 : : integer,intent(in) :: npw,nsh,nsym,timrev
209 : : !arrays
210 : : integer,intent(in) :: g2sh(npw),gvec(3,npw),shlim(nsh+1),symrec(3,3,nsym)
211 : : integer,intent(inout) :: grottb (npw,timrev,nsym)
212 : : integer,intent(inout) :: grottbm1(npw,timrev,nsym)
213 : :
214 : : !Local variables ------------------------------
215 : : !scalars
216 : : integer :: ee,ig1,ig2,ish1,isym,itim,ss
217 : : logical :: found
218 : : character(len=500) :: msg
219 : : !arrays
220 : : integer :: gbase(3),grot(3)
221 : : !************************************************************************
222 : :
223 : : ! === Set up G-rotation table ===
224 [ + + ]: 718560 : do ig1=1,npw
225 : 717475 : ish1=g2sh(ig1) ; ss=shlim(ish1) ; ee=shlim(ish1+1)-1
226 [ + + ]: 2869900 : gbase(:)=gvec(:,ig1)
227 : :
228 [ + + ]: 2147895 : do itim=1,timrev
229 [ + + ]: 35464034 : do isym=1,nsym
230 [ + + + + : 533075584 : grot=(3-2*itim)*MATMUL(symrec(:,:,isym),gbase)
+ + ]
231 : : found=.FALSE.
232 : : ! * Loop on the shell of ig1 to speed up the search.
233 [ + + ]: 2036209336 : do ig2=ss,ee
234 [ + + + + : 2299329072 : if (ALL(ABS(grot(:)-gvec(:,ig2))==0)) then
+ + ]
235 : 33317224 : found=.TRUE.
236 : 33317224 : grottb (ig1,itim,isym)=ig2
237 : 33317224 : grottbm1(ig2,itim,isym)=ig1
238 : : end if
239 : : end do
240 [ - + ]: 34746559 : if (.not.found) then
241 : : write(msg,'(3a,i5,a,i5,1x,2(3i5,a),a,i3,a,i3)')&
242 : 0 : 'G-shell not closed',ch10,&
243 : 0 : ' Initial G vector ',ig1,'/',npw,gbase(:),' Rotated G vector ',grot(:),ch10,&
244 : 0 : ' Through sym ',isym,' and itim ',itim
245 : 0 : ABI_ERROR(msg)
246 : : end if
247 : : end do
248 : : end do
249 : :
250 : : end do !ig1
251 : :
252 : 1085 : end subroutine setup_G_rotation
253 : : !!***
254 : :
255 : : !----------------------------------------------------------------------
256 : :
257 : : !!****f* m_gsphere/gsph_init
258 : : !! NAME
259 : : !! gsph_init
260 : : !!
261 : : !! FUNCTION
262 : : !! Main creation method for the Gvectors data type
263 : : !!
264 : : !! INPUTS
265 : : !! Cryst<crystal_t> = Info on unit cell and its symmetries
266 : : !! ng=number of G vectors, needed only if gvec is passed.
267 : : !! [gvec(3,ng)]=coordinates of G vectors
268 : : !! [ecut]=Cutoff energy for G-sphere. gvec and ecut are mutually exclusive.
269 : : !!
270 : : !! OUTPUT
271 : : !! Gsph<gsphere_t>=Data type containing information related to the set of G vectors
272 : : !! completely initialized in output.
273 : : !!
274 : : !! NOTES
275 : : !! gvec are supposed to be ordered with increasing norm.
276 : : !!
277 : : !! SOURCE
278 : :
279 : 2064 : subroutine gsph_init(Gsph, Cryst, ng, gvec, ecut)
280 : :
281 : : !Arguments ------------------------------------
282 : : !scalars
283 : : class(gsphere_t),intent(out) :: Gsph
284 : : integer,intent(in) :: ng
285 : : real(dp),optional,intent(in) :: ecut
286 : : type(crystal_t),target,intent(in) :: Cryst
287 : :
288 : : !arrays
289 : : integer,optional,intent(in) :: gvec(3,ng)
290 : : !Local variables-------------------------------
291 : : !scalars
292 : : integer,parameter :: nkpt1=1
293 : : integer :: ig,isearch,img,ish,isym,nsh,nsym,timrev,pinv,g1,g2,g3,ss,ee
294 : : real(dp) :: eps,norm,norm_old,max_ecut,gsq
295 : : !arrays
296 : : real(dp),parameter :: k_gamma(3)=(/zero,zero,zero/)
297 : : integer :: sg(3),gsearch(3)
298 : 1085 : integer,allocatable :: shlim(:)
299 : 1085 : integer,pointer :: symrec(:,:,:),gvec_ptr(:,:)
300 : : real(dp) :: kptns1(3,nkpt1)
301 : 1085 : real(dp),allocatable :: shlen(:)
302 : 1085 : real(dp),pointer :: tnons(:,:)
303 : : !************************************************************************
304 : :
305 : : DBG_ENTER("COLL")
306 : :
307 : : ! Copy info on symmetries
308 : 1085 : nsym = Cryst%nsym
309 : 1085 : timrev = Cryst%timrev
310 : 1085 : symrec => Cryst%symrec
311 : 1085 : tnons => Cryst%tnons
312 : :
313 : : ! Initialize the object
314 : 1085 : Gsph%istwfk = 1 ! Time reversal is not used here.
315 : 1085 : Gsph%nsym = nsym
316 : 1085 : Gsph%timrev = timrev
317 : :
318 [ + + + + ]: 14105 : Gsph%gmet = Cryst%gmet
319 [ + + + + ]: 14105 : Gsph%gprimd = Cryst%gprimd
320 : :
321 [ + + ]: 1085 : if (PRESENT(gvec)) then
322 [ - + ]: 979 : if (PRESENT(ecut)) then
323 : 0 : ABI_BUG("ecut cannot be present when gvec is used")
324 : : end if
325 : 979 : Gsph%ng= ng
326 [ - + + - : 2937 : ABI_MALLOC(Gsph%gvec,(3,ng))
+ - + - -
+ - + -
+ ]
327 [ + - - + : 1940058 : Gsph%gvec=gvec
- - + + +
+ ]
328 : : !
329 : : ! Calculate cutoff energy of the sphere.
330 : : max_ecut=-one
331 [ + + ]: 485504 : do ig=1,ng
332 : 484525 : g1=gvec(1,ig)
333 : 484525 : g2=gvec(2,ig)
334 : 484525 : g3=gvec(3,ig)
335 : : gsq= Cryst%gmet(1,1)*g1**2+Cryst%gmet(2,2)*g2**2+Cryst%gmet(3,3)*g3**2+ &
336 : 484525 : two*(Cryst%gmet(1,2)*g1*g2+Cryst%gmet(1,3)*g1*g3+Cryst%gmet(2,3)*g2*g3)
337 : 485504 : max_ecut=MAX(max_ecut,gsq)
338 : : end do
339 : 979 : max_ecut=two*max_ecut*pi**2
340 : 979 : Gsph%ecut= max_ecut
341 : :
342 : : else
343 : : ! To be consistent with the previous implementation.
344 : : !ABI_WARNING("Init from ecut has to be tested")
345 : : !call setshells(ecut,npw,nsh,nsym,Cryst%gmet,Cryst%gprimd,Cryst%symrel,tag,Cryst%ucvol)
346 : 106 : Gsph%ecut = ecut
347 : 106 : pinv=+1; kptns1(:,1)=k_gamma
348 : 106 : call merge_and_sort_kg(nkpt1,kptns1,ecut,Cryst%nsym,pinv,Cryst%symrel,Cryst%gprimd,gvec_ptr,0)
349 : 106 : Gsph%ng = SIZE(gvec_ptr,DIM=2)
350 [ - + + - : 318 : ABI_MALLOC(Gsph%gvec, (3,Gsph%ng))
+ - + - -
+ - + -
+ ]
351 [ + - + - : 932012 : Gsph%gvec = gvec_ptr
- + - - -
- - - + +
+ + ]
352 [ - + ]: 106 : ABI_FREE(gvec_ptr)
353 : : end if
354 : :
355 : : ! Calculate phase exp{-i2\pi G.\tau}
356 [ - + + - : 4340 : ABI_MALLOC(Gsph%phmGt, (Gsph%ng, nsym))
- + + - +
- + - - +
- + - + ]
357 [ + + ]: 44149 : do isym=1,nsym
358 [ + + ]: 16812849 : do ig=1,Gsph%ng
359 [ + + ]: 67117864 : Gsph%phmGt(ig, isym) = EXP(-j_dpc*two_pi*DOT_PRODUCT(Gsph%gvec(:,ig), tnons(:,isym)))
360 : : end do
361 : : end do
362 : : !
363 : : ! === Calculate phase phsgt= exp{-i2\pi SG\cdot t} ===
364 : : ! TODO Here we can store only one of this arrays but I have to rewrite screeening!
365 [ - + + - : 3255 : ABI_MALLOC(Gsph%phmSGt,(Gsph%ng,nsym))
- + + - -
+ - + -
+ ]
366 [ + + ]: 718560 : do ig=1,Gsph%ng
367 [ + + ]: 17487260 : do isym=1,nsym
368 [ + + + + : 519829700 : sg=MATMUL(symrec(:,:,isym),Gsph%gvec(:,ig))
+ - + - +
- ]
369 [ + + ]: 67792275 : Gsph%phmSGt(ig,isym)=EXP(-j_dpc*two_pi*DOT_PRODUCT(sg,tnons(:,isym)))
370 : : end do
371 : : end do
372 : : !
373 : : ! === Calculate number of shells and corresponding starting index ===
374 : : ! * Shells are useful to speed up search algorithms see e.g setup_G_rotation.
375 : : ! * The last shell ends at ng+1, thus gvec is supposed to be closed.
376 : :
377 [ + + + - : 4340 : ABI_CHECK(ALL(Gsph%gvec(1:3,1)==0), 'First G must be 0')
- + ]
378 : :
379 [ - + + - : 3255 : ABI_MALLOC(Gsph%g2sh,(Gsph%ng))
+ - + - -
+ - + -
+ ]
380 : 1085 : Gsph%g2sh(1)=1 ! This table is useful if we dont loop over shell
381 : :
382 : : ! For each shell, gives the index of the initial G-vector.
383 [ + - - + : 3255 : ABI_MALLOC(shlim,(Gsph%ng+1))
+ - + - +
- - + -
+ ]
384 : 1085 : shlim(1)=1
385 : :
386 : : ! For each shell, gives the radius of the shell.
387 [ - + + - : 3255 : ABI_MALLOC(shlen,(Gsph%ng))
+ - + - -
+ - + ]
388 : 1085 : shlen(1)=zero
389 : :
390 : 1085 : nsh=1; norm_old=zero
391 [ + + ]: 717475 : do ig=2,Gsph%ng
392 [ + + + + : 13611410 : norm=two_pi*SQRT(DOT_PRODUCT(Gsph%gvec(:,ig),MATMUL(Cryst%gmet,Gsph%gvec(:,ig))))
+ + + - +
+ ]
393 : 716390 : eps=norm*tol8
394 [ + + ]: 716390 : if (ABS(norm-norm_old)>eps) then
395 : 19213 : norm_old = norm; nsh = nsh + 1
396 : 19213 : shlim(nsh)=ig
397 : 19213 : shlen(nsh)=norm
398 : : end if
399 : 717475 : Gsph%g2sh(ig)=nsh
400 : : end do
401 : 1085 : shlim(nsh+1)=Gsph%ng+1
402 : :
403 : : ! Save info on the shells
404 : 1085 : Gsph%nsh = nsh
405 [ - + + - : 3255 : ABI_MALLOC(Gsph%shlim, (nsh+1))
+ - - + -
+ - + ]
406 [ - + + + ]: 23553 : Gsph%shlim = shlim(1:nsh+1)
407 [ - + + - : 3255 : ABI_MALLOC(Gsph%shlen, (nsh ))
+ - - + -
+ - + ]
408 [ - + + + ]: 22468 : Gsph%shlen = shlen(1:nsh)
409 : 1085 : ABI_FREE(shlim)
410 : 1085 : ABI_FREE(shlen)
411 : :
412 : : ! Calculate tables for rotated G"s
413 [ - + + - : 5425 : ABI_MALLOC(Gsph%rottb , (Gsph%ng,timrev,nsym))
- + + - -
+ + - + -
+ - - + -
+ - + ]
414 [ - + + - : 4340 : ABI_MALLOC(Gsph%rottbm1, (Gsph%ng,timrev,nsym))
- + + - -
+ + - - +
- + - + ]
415 : :
416 : : call setup_G_rotation(nsym, symrec, timrev, Gsph%ng, Gsph%gvec,&
417 [ - + ]: 1085 : Gsph%g2sh, Gsph%nsh, Gsph%shlim, Gsph%rottb, Gsph%rottbm1)
418 : :
419 : : ! Store Mapping G --> -G
420 : : ! (we use a specialized table instead of rootb since rottb assumes time-reversal symmetry.
421 [ - + + - : 3255 : ABI_MALLOC(gsph%g2mg, (gsph%ng))
+ - + - -
+ - + -
+ ]
422 : :
423 [ + + ]: 718560 : do ig=1,gsph%ng
424 : 717475 : ish=gsph%g2sh(ig)
425 : 717475 : ss=gsph%shlim(ish); ee=gsph%shlim(ish+1)-1
426 [ + + ]: 2869900 : gsearch = -gsph%gvec(:,ig)
427 : 25283314 : img = 0
428 : : ! Loop on shells to speed up the search.
429 [ + - ]: 25283314 : do isearch=ss,ee
430 [ + + + + : 29300577 : if (all(gsph%gvec(:,isearch) == gsearch)) then
+ + ]
431 : : img = isearch; exit
432 : : end if
433 : : end do
434 [ - + ]: 717475 : if (img==0) ABI_ERROR("Cannot find -G in G-sphere!")
435 : 718560 : gsph%g2mg(ig) = img
436 : : end do
437 : :
438 : : !call Gsph%print(unit=std_out,prtvol=1)
439 : : DBG_EXIT("COLL")
440 : :
441 [ + + + - ]: 2064 : end subroutine gsph_init
442 : : !!***
443 : :
444 : : !----------------------------------------------------------------------
445 : :
446 : : !!****f* m_gsphere/gsph_fft_tabs
447 : : !! NAME
448 : : !! gsph_fft_tabs
449 : : !!
450 : : !! FUNCTION
451 : : !!
452 : : !! INPUTS
453 : : !! Gsph<gsphere_t>=Info on the G-sphere
454 : : !! g0(3)
455 : : !! mgfft=MAXVAL(ngfft(1:3))
456 : : !! ngfftf(18)=Info on the FFT mesh.
457 : : !!
458 : : !! OUTPUT
459 : : !! use_padfft=1 if padded FFT can be used, 0 otherwise.
460 : : !! gmg0_gbound(2*mgfft+8,2)=Tables for improved zero-padded FFTS. Calculated only if use_padfft==1
461 : : !! gmg0_ifft(Gsph%ng)=Index of G-G0 in the FFT mesh defined by ngfft.
462 : : !!
463 : : !! NOTES
464 : : !! The routine will stop if any G-G0 happens to be outside the FFT box.
465 : : !!
466 : : !! SOURCE
467 : :
468 : 104550 : subroutine gsph_fft_tabs(Gsph, g0, mgfft, ngfft, use_padfft, gmg0_gbound, gmg0_ifft)
469 : :
470 : : !Arguments ------------------------------------
471 : : !scalars
472 : : class(gsphere_t),intent(in) :: Gsph
473 : : integer,intent(in) :: mgfft
474 : : integer,intent(out) :: use_padfft
475 : : !arrays
476 : : integer,intent(in) :: g0(3),ngfft(18)
477 : : integer,intent(out) :: gmg0_gbound(2*mgfft+8,2),gmg0_ifft(Gsph%ng)
478 : :
479 : : !Local variables-------------------------------
480 : : !scalars
481 : : integer :: ig,ng,ierr
482 : : character(len=500) :: msg
483 : 104550 : type(MPI_type) :: MPI_enreg_seq
484 : : !arrays
485 : 104550 : integer,allocatable :: gmg0(:,:)
486 : 104550 : logical,allocatable :: kg_mask(:)
487 : : ! *************************************************************************
488 : :
489 [ + + - + ]: 418200 : if (mgfft/=MAXVAL(ngfft(1:3))) then
490 : 0 : ABI_ERROR("mgfft/-MAXVAL(ngfft(1:3)")
491 : : end if
492 : :
493 : 104550 : ng = Gsph%ng
494 : :
495 : 104550 : ierr=0; use_padfft=0
496 [ - + + - : 313650 : ABI_MALLOC(gmg0,(3,ng))
+ - + - -
+ - + ]
497 [ + + ]: 6089462 : do ig=1,ng
498 [ + + ]: 23939648 : gmg0(:,ig) = Gsph%gvec(:,ig)-g0
499 : : ! Consider possible wrap around errors.
500 [ + + + - : 47879296 : if ( ANY(gmg0(:,ig)>ngfft(1:3)/2) .or. ANY(gmg0(:,ig)<-(ngfft(1:3)-1)/2) ) then
+ + + - -
+ ]
501 : : !gmg0_ifft(ig,ig01+mg0(1)+1,ig02+mg0(2)+1,ig03+mg0(3)+1) = 0
502 : 0 : write(std_out,*)" outside FFT box ",gmg0(:,ig)
503 : 0 : ierr=ierr+1
504 : : end if
505 [ + + + + : 8073152 : if (ALL(gmg0(:,ig) == 0)) use_padfft=1
+ + ]
506 : : end do
507 : :
508 [ - + ]: 104550 : if (ierr/=0) then
509 : 0 : write(msg,'(a,i0,a)')'Found ',ierr,' G-G0 vectors falling outside the FFT box. This is not allowed '
510 : 0 : ABI_ERROR(msg)
511 : : end if
512 : : !
513 : : ! Evaluate the tables needed for the padded FFT performed in rhotwg. Note that we have
514 : : ! to pass G-G0 to sphereboundary instead of G as we need FFT results on the shifted G-sphere,
515 : : ! If Gamma is not inside G-G0 one has to disable FFT padding as sphereboundary will give wrong tables.
516 [ + + ]: 104550 : if (use_padfft == 1) call sphereboundary(gmg0_gbound,1,gmg0,mgfft,ng)
517 : :
518 : 104550 : call initmpi_seq(MPI_enreg_seq) ! No FFT parallelism.
519 : 104550 : call MPI_enreg_seq%distribfft%init_seq('c',ngfft(2),ngfft(3),'all')
520 : :
521 [ - + + - : 313650 : ABI_MALLOC(kg_mask, (ng))
+ - + - -
+ - + ]
522 : 104550 : call kgindex(gmg0_ifft, gmg0, kg_mask, MPI_enreg_seq, ngfft, ng)
523 : :
524 [ + + + - : 6089462 : ABI_CHECK(ALL(kg_mask),"FFT para not yet implemented")
- + ]
525 [ - + ]: 104550 : ABI_FREE(kg_mask)
526 : :
527 [ - + ]: 104550 : ABI_FREE(gmg0)
528 : 104550 : call destroy_mpi_enreg(MPI_enreg_seq)
529 : :
530 [ - + - + : 104550 : end subroutine gsph_fft_tabs
- + - + -
+ - + - +
- + - + -
+ - + ]
531 : : !!***
532 : :
533 : : !----------------------------------------------------------------------
534 : :
535 : : !!****f* m_gsphere/gsph_in_fftbox
536 : : !! NAME
537 : : !! gsph_in_fftbox
538 : : !!
539 : : !! FUNCTION
540 : : !! Initialize the largest Gsphere contained in the FFT box.
541 : : !!
542 : : !! INPUTS
543 : : !! Cryst<crystal_t> = Info on unit cell and its symmetries.
544 : : !! ngfft(18)=Info on the FFT box.
545 : : !!
546 : : !! OUTPUT
547 : : !! Gsph<gsphere_t>=Data type containing information related to the set of G vectors
548 : : !! completely initialized in output.
549 : : !!
550 : : !! SOURCE
551 : :
552 : 39 : subroutine gsph_in_fftbox(Gsph, Cryst, ngfft)
553 : :
554 : : !Arguments ------------------------------------
555 : : !scalars
556 : : class(gsphere_t),intent(out) :: Gsph
557 : : type(crystal_t),intent(in) :: Cryst
558 : : !arrays
559 : : integer,intent(in) :: ngfft(18)
560 : :
561 : : !Local variables-------------------------------
562 : : !scalars
563 : : integer :: dir1,dir2,dir3,npw,ig,i_st
564 : : real(dp) :: ecut,trial_ene
565 : : !arrays
566 : : integer :: n1_max(3),n2_max(3),n3_max(3),vec(3)
567 [ + - ]: 39 : integer,allocatable :: gvec(:,:)
568 : : !************************************************************************
569 : :
570 : : ! Find ecut for the largest G-sphere contained in the FFT box.
571 : 39 : n1_max(1) = -(ngfft(1)-1)/2
572 : 39 : n2_max(1) = -(ngfft(2)-1)/2
573 : 39 : n3_max(1) = -(ngfft(3)-1)/2
574 : :
575 : 39 : n1_max(2) = 0
576 : 39 : n2_max(2) = 0
577 : 39 : n3_max(2) = 0
578 : :
579 : 39 : n1_max(3) = ngfft(1)/2
580 : 39 : n2_max(3) = ngfft(2)/2
581 : 39 : n3_max(3) = ngfft(3)/2
582 : :
583 : 39 : ecut = HUGE(one)
584 [ + + ]: 156 : do dir3=1,3
585 : 117 : vec(3) = n1_max(dir3)
586 [ + + ]: 507 : do dir2=1,3
587 : 351 : vec(2) = n2_max(dir2)
588 [ + + ]: 1521 : do dir1=1,3
589 : 1053 : vec(1) = n1_max(dir1)
590 [ + + + + : 1911 : if (ANY(vec/=0)) then
+ + ]
591 : 1014 : trial_ene = half * normv(vec,Cryst%gmet,"G")**2
592 : 1014 : ecut = MIN(ecut,trial_ene)
593 : : !write(std_out,*)vec(:),trial_ene
594 : : end if
595 : : end do
596 : : end do
597 : : end do
598 : : !
599 : : ! Init sphere from ecut.
600 : 39 : call Gsph%init(Cryst, 0, ecut=ecut)
601 : : !
602 : : ! Make sure that Gsph does not contain G vectors outside the FFT box.
603 : : ! kpgsph might return G whose energy is larger than the input ecut.
604 : 39 : npw = Gsph%ng
605 [ + + ]: 3013 : star_loop: do i_st=1,Gsph%nsh-1
606 [ + + ]: 159018 : do ig=Gsph%shlim(i_st),Gsph%shlim(i_st+1)
607 [ + + + - : 1251266 : if ( ANY(Gsph%gvec(:,ig)>ngfft(1:3)/2) .or. ANY(Gsph%gvec(:,ig)<-(ngfft(1:3)-1)/2) ) then
+ + + + +
+ ]
608 : 36 : npw = Gsph%shlim(i_st)-1 ! Gsph exceeds the FFT box. Only the shells up to npw will be used.
609 : 36 : EXIT star_loop
610 : : end if
611 : : end do
612 : : end do star_loop
613 : :
614 [ + + ]: 39 : if (npw<Gsph%ng) then
615 : 36 : ABI_COMMENT("Have to reinit Gpshere")
616 [ - + + - : 108 : ABI_MALLOC(gvec,(3,npw))
+ - + - -
+ - + ]
617 [ + - - + : 610328 : gvec = Gsph%gvec(:,1:npw)
- - + + +
+ ]
618 : 36 : call Gsph%free()
619 : 36 : call Gsph%init(Cryst, npw, gvec=gvec)
620 [ - + ]: 36 : ABI_FREE(gvec)
621 : : end if
622 : :
623 : 39 : end subroutine gsph_in_fftbox
624 : : !!***
625 : :
626 : : !----------------------------------------------------------------------
627 : :
628 : : !!****f* m_gsphere/gsph_print
629 : : !! NAME
630 : : !! gsph_print
631 : : !!
632 : : !! FUNCTION
633 : : !! Print info on object.
634 : : !!
635 : : !! INPUTS
636 : : !! unit=the unit number for output
637 : : !! prtvol = verbosity level
638 : : !!
639 : : !! SOURCE
640 : :
641 : 144 : subroutine gsph_print(Gsph, units, prtvol, header)
642 : :
643 : : !Arguments ------------------------------------
644 : : !scalars
645 : : class(gsphere_t),intent(in) :: Gsph
646 : : integer,intent(in) :: units(:), prtvol
647 : : character(len=*),optional,intent(in) :: header
648 : :
649 : : !Local variables-------------------------------
650 : : !scalars
651 : : integer :: ish, nsc
652 : : real(dp) :: fact, kin
653 : : character(len=500) :: msg
654 : : ! *************************************************************************
655 : :
656 : 72 : msg = ch10 // ' ==== Info on the G-sphere ==== ' // ch10
657 [ + + + - : 72 : if (present(header)) msg = ' ==== '//trim(adjustl(header))//' ==== '
+ - ]
658 : 72 : call wrtout(units, msg)
659 : :
660 : : write(msg,'(2(a,i8,a))')&
661 : 72 : ' Number of G vectors ... ',Gsph%ng,ch10,&
662 : 144 : ' Number of shells ...... ',Gsph%nsh,ch10
663 : 72 : call wrtout(units, msg)
664 : :
665 : 80 : select case (Gsph%timrev)
666 : : case (1)
667 : 8 : call wrtout(units, ' Time reversal symmetry cannot be used')
668 : : case (2)
669 : 64 : call wrtout(units, ' Time reversal symmetry is used')
670 : : case default
671 [ + + - ]: 72 : ABI_BUG("Wrong timrev")
672 : : end select
673 : :
674 [ + + ]: 72 : if (prtvol /= 0) then
675 : 2 : fact = half*two_pi**2
676 : 2 : write(msg,'(a)')
677 : 2 : call wrtout(units, ' Shell Tot no. of Gs Cutoff [Ha]')
678 [ + + ]: 88 : do ish=1,Gsph%nsh
679 : 86 : nsc=Gsph%shlim(ish+1)-1
680 : 86 : kin=half*Gsph%shlen(ish)**2
681 : 86 : write(msg, '(2x,i4,10x,i6,5x,f8.3)')ish,nsc,kin
682 : 88 : call wrtout(units, msg)
683 : : end do
684 : 2 : call wrtout(units, ch10)
685 : : end if
686 : :
687 [ - + ]: 72 : end subroutine gsph_print
688 : : !!***
689 : :
690 : : !----------------------------------------------------------------------
691 : :
692 : : !!****f* m_gsphere/gsph_free
693 : : !! NAME
694 : : !! gsph_free
695 : : !!
696 : : !! FUNCTION
697 : : !! Deallocate the memory in a gsphere_t data type.
698 : : !!
699 : : !! SOURCE
700 : :
701 : 1977 : subroutine gsph_free(Gsph)
702 : :
703 : : !Arguments ------------------------------------
704 : : class(gsphere_t),intent(inout) :: Gsph
705 : : ! *************************************************************************
706 : :
707 : : DBG_ENTER("COLL")
708 : :
709 : : !@gsphere_t
710 : :
711 : : ! integer arrays.
712 [ + + ]: 1977 : ABI_SFREE(Gsph%g2sh)
713 [ + + ]: 1977 : ABI_SFREE(Gsph%gvec)
714 [ + + ]: 1977 : ABI_SFREE(Gsph%g2mg)
715 [ + + ]: 1977 : ABI_SFREE(Gsph%rottb)
716 [ + + ]: 1977 : ABI_SFREE(Gsph%rottbm1)
717 [ + + ]: 1977 : ABI_SFREE(Gsph%shlim)
718 : :
719 [ + + ]: 1977 : ABI_SFREE(Gsph%shlen)
720 : :
721 : : ! complex arrays
722 [ + + ]: 1977 : ABI_SFREE(Gsph%phmGt)
723 [ + + ]: 1977 : ABI_SFREE(Gsph%phmSGt)
724 : :
725 : : DBG_EXIT("COLL")
726 : :
727 : 1977 : end subroutine gsph_free
728 : : !!***
729 : :
730 : : !----------------------------------------------------------------------
731 : :
732 : : !!****f* m_gsphere/gsph_g_idx
733 : : !! NAME
734 : : !! gsph_g_idx
735 : : !!
736 : : !! FUNCTION
737 : : !! Return the index of G in the sphere. zero if not in the sphere
738 : : !!
739 : : !! INPUTS
740 : : !! Gsph<gsphere_t>=Info on the G-sphere
741 : : !! gg(3)=Reduced coordinates of the G-vector.
742 : : !!
743 : : !! NOTES
744 : : !! The function assumes that the G-vectors are ordered with increasing length.
745 : : !!
746 : : !! SOURCE
747 : :
748 : 0 : pure function gsph_g_idx(Gsph, gg) result(g_idx)
749 : :
750 : : !Arguments ------------------------------------
751 : : !scalars
752 : : class(gsphere_t),intent(in) :: Gsph
753 : : integer :: g_idx
754 : : !arrays
755 : : integer,intent(in) :: gg(3)
756 : :
757 : : !Local variables-------------------------------
758 : : !scalars
759 : : integer :: ishbsc,igs,ige
760 : : real(dp) :: glen
761 : : logical :: found
762 : : ! *************************************************************************
763 : :
764 : : ! Use shells and bisection to find the star and stop index thus avoiding the storage of a table (ig1,ig2)
765 [ # # # # : 0 : glen = two_pi*SQRT(DOT_PRODUCT(gg,MATMUL(Gsph%gmet,gg)))
# # # # ]
766 : :
767 : 0 : ishbsc = bisect(Gsph%shlen,glen)
768 [ # # # # : 0 : if ( ANY(ishbsc==(/0,Gsph%nsh/)) ) then ! glen out of range.
# # ]
769 : 0 : g_idx=0; RETURN
770 : : end if
771 : :
772 : 0 : igs = Gsph%shlim(ishbsc)
773 : 0 : ige = Gsph%shlim(MIN(ishbsc+2,Gsph%nsh+1))-1
774 : :
775 : 0 : g_idx=igs-1; found=.FALSE.
776 [ # # ]: 0 : do while (.not.found .and. g_idx<ige)
777 : 0 : g_idx=g_idx+1
778 [ # # # # ]: 0 : found=(ALL(Gsph%gvec(:,g_idx)==gg(:)))
779 : : end do
780 [ # # ]: 0 : if (.not.found) g_idx=0
781 : :
782 : 0 : end function gsph_g_idx
783 : : !!***
784 : :
785 : : !----------------------------------------------------------------------
786 : :
787 : : !!****f* m_gsphere/gsph_gmg_idx
788 : : !! NAME
789 : : !! gsph_gmg_idx
790 : : !!
791 : : !! FUNCTION
792 : : !! Return the index of G1-G2 in the sphere. zero if not in the sphere
793 : : !!
794 : : !! INPUTS
795 : : !! Gsph<gsphere_t>=Info on the G-sphere
796 : : !! ig1,ig2 index of g1 and g2 in the G-sphere.
797 : : !!
798 : : !! NOTES
799 : : !! The function assumes that the G-vectors are ordered with increasing length.
800 : : !!
801 : : !! SOURCE
802 : :
803 : 6422840 : pure function gsph_gmg_idx(Gsph, ig1, ig2) result(ig1mg2)
804 : :
805 : : !Arguments ------------------------------------
806 : : !scalars
807 : : class(gsphere_t),intent(in) :: Gsph
808 : : integer,intent(in) :: ig1,ig2
809 : : integer :: ig1mg2
810 : :
811 : : !Local variables-------------------------------
812 : : !scalars
813 : : integer :: ishbsc,igs,ige
814 : : real(dp) :: difflen
815 : : logical :: found
816 : : !arrays
817 : : integer :: g1mg2(3)
818 : : ! *************************************************************************
819 : :
820 [ + + ]: 25691360 : g1mg2 = Gsph%gvec(:,ig1)-Gsph%gvec(:,ig2)
821 : :
822 : : ! Use shells and bisect to find the star and stop index thus avoiding the storage of a table (ig1,ig2)
823 [ + + + + : 122033960 : difflen = two_pi*SQRT(DOT_PRODUCT(g1mg2,MATMUL(Gsph%gmet,g1mg2)))
+ + + + ]
824 : :
825 : : ! FIXME It seems bisect is not portable, on my laptop test v5/t72 the number of skipped G-vectors is > 0
826 : 6422840 : ishbsc = bisect(Gsph%shlen,difflen)
827 [ + - - + ]: 6422840 : if ( ANY(ishbsc==(/0,Gsph%nsh/)) ) then ! difflen out of range.
828 : 0 : ig1mg2=0; RETURN
829 : : end if
830 : :
831 : 6422840 : igs = Gsph%shlim(ishbsc)
832 : 6422840 : ige = Gsph%shlim(MIN(ishbsc+2,Gsph%nsh+1))-1
833 : :
834 : 6422840 : ig1mg2=igs-1; found=.FALSE.
835 [ + + ]: 78583046 : do while (.not.found .and. ig1mg2<ige)
836 : 72160206 : ig1mg2=ig1mg2+1
837 [ + + + + ]: 108457970 : found=(ALL(Gsph%gvec(:,ig1mg2)==g1mg2(:)))
838 : : end do
839 [ - + ]: 6422840 : if (.not.found) ig1mg2=0
840 : :
841 : 6422840 : end function gsph_gmg_idx
842 : : !!***
843 : :
844 : : !----------------------------------------------------------------------
845 : :
846 : : !!****f* m_gsphere/gsph_gmg_fftidx
847 : : !! NAME
848 : : !! gsph_gmg_fftidx
849 : : !!
850 : : !! FUNCTION
851 : : !! Return the index of G1-G2 in the FFT mesh defined by ngfft. zero if not found.
852 : : !!
853 : : !! INPUTS
854 : : !! Gsph<gsphere_t>=Info on the G-sphere
855 : : !! ig1,ig2 index of g1 and g2 in the G-sphere.
856 : : !! ngfft(18)=Info on the FFT mesh.
857 : : !!
858 : : !! SOURCE
859 : :
860 : 3586176 : pure function gsph_gmg_fftidx(Gsph, ig1, ig2, ngfft) result(fft_idx)
861 : :
862 : : !Arguments ------------------------------------
863 : : !scalars
864 : : class(gsphere_t),intent(in) :: Gsph
865 : : integer,intent(in) :: ig1,ig2
866 : : integer :: fft_idx
867 : : !arrays
868 : : integer,intent(in) :: ngfft(18)
869 : :
870 : : !Local variables-------------------------------
871 : : !scalars
872 : : integer :: id1,id2,id3
873 : : !arrays
874 : : integer :: g1mg2(3)
875 : : ! *************************************************************************
876 : :
877 [ + + ]: 14344704 : g1mg2(:)=Gsph%gvec(:,ig1)-Gsph%gvec(:,ig2)
878 : :
879 : : ! Make sure G1-G2 is still in the FFT mesh.
880 : : ! MODULO wraps G1-G2 in the FFT box but the Fourier components are not periodic!
881 [ + + + - : 28689408 : if (ANY(g1mg2(:)>ngfft(1:3)/2) .or. ANY(g1mg2(:)<-(ngfft(1:3)-1)/2)) then
+ + + - +
- ]
882 : 3586176 : fft_idx=0; RETURN
883 : : end if
884 : :
885 : 3586176 : id1=MODULO(g1mg2(1),ngfft(1))
886 : 3586176 : id2=MODULO(g1mg2(2),ngfft(2))
887 : 3586176 : id3=MODULO(g1mg2(3),ngfft(3))
888 : 3586176 : fft_idx= 1 + id1 + id2*ngfft(1) + id3*ngfft(1)*ngfft(2)
889 : :
890 : 3586176 : end function gsph_gmg_fftidx
891 : : !!***
892 : :
893 : : !----------------------------------------------------------------------
894 : :
895 : : !!****f* m_gsphere/merge_and_sort_kg
896 : : !! NAME
897 : : !! merge_and_sort_kg
898 : : !!
899 : : !! FUNCTION
900 : : !! This routine merges a set of k-centered G-spheres of cutoff energy ecut and
901 : : !! returns a Gamma-centered G-spheres. The elements in the final G-spheres are packed with increasing module.
902 : : !!
903 : : !! INPUTS
904 : : !! nkpt=Number of k-points
905 : : !! kptns(3,nkpt)=The k-points in reduced coordinates defining the k-centered G-spheres.
906 : : !! ecut=Cutoff energy for the k-centered G-spheres.
907 : : !! nsym2=Number of symmetry operations.
908 : : !! pinv=-1 if time-reversal can be used, 1 otherwise
909 : : !! symrel2(3,3,nsym2)=symmetry operations in real space.
910 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
911 : : !! prtvol=Flag defining the verbosity level.
912 : : !!
913 : : !! SIDE EFFECTS
914 : : !! gbig(:,:)
915 : : !! in input : pointer to NULL
916 : : !! in output: gbig(3,1:npw) contains the set of G-vectors ordered by shell obtained by
917 : : !! merging the k-centered sphere.
918 : : !! shlim_p(:)
919 : : !! in input : pointer to NULL
920 : : !! in output: shlim_p(nbase)=Cumulative number of G-vectors for each shell.
921 : : !! where nbase is the number of irreducible G"s found.
922 : : !!
923 : : !! SOURCE
924 : :
925 : 476 : subroutine merge_and_sort_kg(nkpt,kptns,ecut,nsym2,pinv,symrel2,gprimd,gbig,prtvol,shlim_p)
926 : :
927 : : !Arguments ------------------------------------
928 : : !scalars
929 : : integer,intent(in) :: nkpt,nsym2,pinv,prtvol
930 : : real(dp),intent(in) :: ecut
931 : : !arrays
932 : : integer,intent(in) :: symrel2(3,3,nsym2)
933 : : real(dp),intent(in) :: kptns(3,nkpt),gprimd(3,3)
934 : : integer,pointer :: gbig(:,:)
935 : : integer,optional,pointer :: shlim_p(:)
936 : :
937 : : !Local variables-------------------------------
938 : : !scalars
939 : : integer,parameter :: mkmem_=1
940 : : integer :: ikg,ig,ikpt,nbase,sizepw,in,maxpw,is,iinv,ish,ilim,mpw
941 : : integer :: exchn2n3d,istwf_k,onpw_k,ierr,npw_k,ii,isym,sizeold
942 : : logical :: found
943 : : character(len=500) :: msg
944 : 476 : type(MPI_type) :: MPI_enreg_seq
945 : : !arrays
946 : 952 : integer :: gcur(3),geq(3),dum_kg(3,0),symrec2t(3,3,nsym2)
947 : 476 : integer,allocatable :: gbase(:,:),gbasek(:,:,:)
948 : 476 : integer,allocatable :: gcurr(:,:),gshell(:,:),insort(:),gtmp(:,:)
949 : 476 : integer,allocatable :: nbasek(:),nshell(:),shlim(:), npwarr(:)
950 : : real(dp) :: kpoint(3),gmet(3,3)
951 : 476 : real(dp),allocatable :: cnorm(:),cnormk(:,:),ctmp(:)
952 : : ! *********************************************************************
953 : :
954 : : ! Fake MPI_type for the sequential part.
955 : : ! This routine should not be parallelized as communicating gbig and other
956 : : ! tables takes more time than recalculating them in sequential.
957 : 476 : call initmpi_seq(MPI_enreg_seq)
958 : :
959 : : !Compute reciprocal space metrics
960 [ + + ]: 1904 : do ii=1,3
961 : : gmet(ii,:)=gprimd(1,ii)*gprimd(1,:)+&
962 : : & gprimd(2,ii)*gprimd(2,:)+&
963 [ + + ]: 6188 : & gprimd(3,ii)*gprimd(3,:)
964 : : end do
965 : :
966 : : ! * Here we use TRANSPOSE(symrel2) instead of the more intuitive symrel2^{-1t} for historical reasons
967 : : ! It does not affect the results since in the code below we only check the module of G
968 [ + + ]: 19382 : do isym=1,nsym2
969 [ + + + + ]: 246254 : symrec2t(:,:,isym)=TRANSPOSE(symrel2(:,:,isym))
970 : : end do
971 : : !
972 : : ! ==============================================
973 : : ! ==== Find irreducible G-vectors at each k ====
974 : : ! ==============================================
975 : :
976 [ - + + - : 1428 : ABI_MALLOC(npwarr,(nkpt))
+ - + - -
+ - + ]
977 : 476 : exchn2n3d=0; ikg=0
978 [ + + ]: 2749 : do ikpt=1,nkpt
979 [ + + ]: 9092 : kpoint=kptns(:,ikpt); istwf_k=1
980 : 2749 : call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,dum_kg,kpoint,0,MPI_enreg_seq,0,npwarr(ikpt))
981 : : end do
982 [ + + ]: 2749 : mpw = MAXVAL(npwarr)
983 : :
984 [ - + + - : 952 : ABI_MALLOC(nbasek,(nkpt))
- + - + ]
985 [ - + + - : 1904 : ABI_MALLOC(gbasek,(3,mpw,nkpt))
- + + - +
- + - - +
- + ]
986 [ - + + - : 1904 : ABI_MALLOC(cnormk,(mpw,nkpt))
- + + - +
- + - - +
+ - ]
987 [ + + ]: 2749 : nbasek=0 ! # of irreducible G at each k.
988 [ + + + + ]: 908043 : cnormk=zero ! Norm of each irreducible G.
989 [ + + + + : 3623925 : gbasek=0 ! The set of irreducible G"s at each k.
+ + ]
990 : :
991 [ + + ]: 2749 : do ikpt=1,nkpt
992 : :
993 [ + + ]: 9092 : kpoint = kptns(:,ikpt)
994 : 2273 : npw_k = npwarr(ikpt)
995 : :
996 : 2273 : exchn2n3d=0; ikg=0; istwf_k=1
997 [ - + + - : 6819 : ABI_MALLOC(gcurr,(3,npw_k))
+ - + - -
+ - + ]
998 : 2273 : call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,gcurr,kpoint,mkmem_,MPI_enreg_seq,npw_k,onpw_k)
999 : :
1000 [ + + + - : 9092 : if (ANY(gcurr(:,1)/=0)) then
- + ]
1001 : 0 : ABI_BUG("gcurr(:,1)/=0")
1002 : : end if
1003 : : !
1004 : : ! * Search for the G"s generating the others by symmetry.
1005 : : ! NB: Here we use symrec2t=TRANSPOSE(symrel2) for historical reasons, see note above
1006 : 2273 : call get_irredg(npw_k,nsym2,pinv,gprimd,symrec2t,gcurr,nbasek(ikpt),gbasek(:,:,ikpt),cnormk(:,ikpt))
1007 : :
1008 [ - + ]: 2749 : ABI_FREE(gcurr)
1009 : : end do
1010 : : !
1011 : : ! === Reduce info over k-points ===
1012 : : ! * Here symrec2t=TRANSPOSE(symrel2) for historical reasons, see note above
1013 : 476 : sizepw=2*mpw
1014 [ - + + - : 1428 : ABI_MALLOC(gbase,(3,sizepw))
+ - + - -
+ - + ]
1015 [ - + + - : 1428 : ABI_MALLOC(cnorm,(sizepw))
+ - + - -
+ - + ]
1016 : : nbase=0 ! # of irred G found.
1017 : :
1018 : 476 : call merge_kgirr(nsym2,pinv,nkpt,mpw,sizepw,symrec2t,nbasek,cnormk,gbasek,nbase,gbase,cnorm,ierr)
1019 [ - + ]: 476 : if (ierr/=0) then
1020 : 0 : ABI_ERROR('merge_kgirr returned a non-zero status error')
1021 : : end if
1022 : :
1023 [ - + ]: 476 : ABI_FREE(nbasek)
1024 [ - + ]: 476 : ABI_FREE(cnormk)
1025 [ - + ]: 476 : ABI_FREE(gbasek)
1026 : : !
1027 : : !=== Reorder base G-vectors in order of increasing module ===
1028 : : !
1029 : : !Generate all shells of G-vectors: star of a g==set of all symmetrics of this g
1030 [ + - - + : 1428 : ABI_MALLOC(gshell,(3,2*nsym2))
+ - + - +
- - + -
+ ]
1031 [ - + + - : 1428 : ABI_MALLOC(shlim,(nbase))
+ - + - -
+ - + ]
1032 : :
1033 [ - + + - : 952 : ABI_MALLOC(gbig,(3,sizepw))
- + - + ]
1034 : : !
1035 : : !TODO
1036 : : #if 0
1037 : : !* Here symrec2t=TRANSPOSE(symrel2) for historical reasons, see note above
1038 : : call getfullg(nbase,nsym2,pinv,sizepw,gbase,symrec2t,cnorm,maxpw,gbig,shlim,ierr)
1039 : : if (ierr/0) RETURN
1040 : :
1041 : : #else
1042 [ - + + - : 952 : ABI_MALLOC(insort,(nbase))
- + - + ]
1043 [ - + + - : 952 : ABI_MALLOC(nshell,(nbase))
- + - + ]
1044 [ + + ]: 37595 : do in=1,nbase
1045 : 37595 : insort(in)=in
1046 : : end do
1047 : 476 : call sort_dp(nbase,cnorm,insort,tol14)
1048 : : !
1049 : : !Loop over all different modules of g''s (=shells):
1050 : 476 : maxpw=0
1051 [ + + ]: 37595 : do in=1,nbase
1052 : 37119 : nshell(in)=0
1053 [ + + ]: 148476 : gcur(:)=gbase(:,insort(in))
1054 : :
1055 [ + + ]: 858153 : do is=1,nsym2 ! Loop over all symmetries:
1056 [ + - + - ]: 1679187 : do iinv=pinv,1,2
1057 [ + + ]: 3284136 : geq(:)=iinv*(symrel2(1,:,is)*gcur(1)+symrel2(2,:,is)*gcur(2)+symrel2(3,:,is)*gcur(3))
1058 : :
1059 : : found=.FALSE.; ish=1
1060 [ + + + + ]: 7079526 : do while ((.not.found) .and. (ish<=nshell(in))) ! Search for symmetric of g and eventually add it:
1061 [ + + + + ]: 8262240 : found=ALL(geq(:)==gshell(:,ish))
1062 : 6756998 : ish=ish+1
1063 : : end do
1064 [ + + ]: 1642068 : if (.not.found) then
1065 : 498506 : nshell(in)=nshell(in)+1
1066 [ + + ]: 1994024 : gshell(:,nshell(in))=geq(:)
1067 : : end if
1068 : : end do
1069 : : end do
1070 : :
1071 [ - + ]: 37119 : if ((maxpw+nshell(in)) > sizepw) then
1072 : : ! We need to increase the size of the gbase, gbig and cnorm arrays while still keeping their content.
1073 : : ! This is done using two temporary arrays gtmp and ctmp
1074 : 0 : ABI_WARNING("Had to reallocate gbase, gbig, cnorm. Perhaps geometry too inaccurate. Possible fix: correct your input file.")
1075 [ # # # # : 0 : ABI_MALLOC(ctmp,(sizepw))
# # # # #
# # # ]
1076 [ # # # # : 0 : ABI_MALLOC(gtmp,(3,sizepw))
# # # # #
# # # ]
1077 : 0 : sizeold=sizepw
1078 : 0 : sizepw=maxpw+nshell(in)
1079 : :
1080 [ # # ]: 0 : ctmp(:)=cnorm(:)
1081 [ # # # # ]: 0 : gtmp(:,:)=gbase(:,:)
1082 : :
1083 [ # # ]: 0 : ABI_FREE(cnorm)
1084 [ # # # # : 0 : ABI_MALLOC(cnorm,(sizepw))
# # # # #
# # # ]
1085 [ # # ]: 0 : cnorm(1:sizeold)=ctmp(1:sizeold)
1086 [ # # ]: 0 : cnorm(sizeold+1:sizepw)=zero
1087 [ # # ]: 0 : ABI_FREE(ctmp)
1088 : :
1089 : : ! MG why this? gbase should not be changed!
1090 [ # # ]: 0 : ABI_FREE(gbase)
1091 [ # # # # : 0 : ABI_MALLOC(gbase,(3,sizepw))
# # # # #
# # # ]
1092 [ # # # # ]: 0 : gbase(:,:sizeold)=gtmp(:,:sizeold)
1093 [ # # # # ]: 0 : gbase(:,sizeold+1:sizepw)=0
1094 [ # # # # ]: 0 : gtmp(:,:)=gbig(:,:)
1095 : :
1096 [ # # ]: 0 : ABI_FREE(gbig)
1097 [ # # # # : 0 : ABI_MALLOC(gbig,(3,sizepw))
# # # # ]
1098 [ # # # # ]: 0 : gbig(:,:sizeold)=gtmp(:,:sizeold)
1099 [ # # # # ]: 0 : gbig(:,sizeold+1:sizepw)=0
1100 [ # # ]: 0 : ABI_FREE(gtmp)
1101 : : end if
1102 : : !
1103 : : ! Store this shell of g''s in a big array of g (gbig):
1104 [ + + ]: 535625 : do ig=1,nshell(in)
1105 [ + + ]: 2031143 : gbig(:,ig+maxpw)=gshell(:,ig)
1106 : : end do
1107 : 37595 : maxpw=maxpw+nshell(in)
1108 : :
1109 : : end do ! End loop over shells
1110 : : !
1111 : : ! * Compute shell limits
1112 : 476 : ilim=0
1113 [ + + ]: 37595 : do in=1,nbase
1114 : 37119 : ilim=ilim+nshell(in)
1115 : 37595 : shlim(in)=ilim
1116 : : end do
1117 : :
1118 [ + + ]: 476 : if (PRESENT(shlim_p)) then ! Return shlim_p
1119 [ - + + - : 740 : ABI_MALLOC(shlim_p,(nbase))
- + - + ]
1120 [ + + ]: 21075 : shlim_p = shlim
1121 : : end if
1122 : :
1123 : : ! Re-allocate gbig with correct sizes so that caller can inquire the size
1124 [ - + + - : 1428 : ABI_MALLOC(gtmp,(3,ilim))
+ - + - -
+ - + ]
1125 [ + - - + : 1994976 : gtmp = gbig(:,1:ilim)
- - + + +
+ ]
1126 [ - + ]: 476 : ABI_FREE(gbig)
1127 [ - + + - : 952 : ABI_MALLOC(gbig,(3,ilim))
- + - + ]
1128 [ + + + + ]: 1994500 : gbig=gtmp
1129 [ - + ]: 476 : ABI_FREE(gtmp)
1130 : :
1131 [ - + ]: 476 : if (prtvol>10) then ! Print out shell limits
1132 : : write(msg,'(3a)')&
1133 : 0 : & ' Shells found:',ch10,&
1134 : 0 : & ' number of shell number of G vectors cut-off energy [Ha} '
1135 : 0 : call wrtout(std_out,msg)
1136 [ # # ]: 0 : do in=1,nbase
1137 : 0 : write(msg,'(12x,i4,17x,i6,12x,f8.3)')in,shlim(in),2*pi**2*cnorm(in)
1138 : 0 : call wrtout(std_out,msg)
1139 : : end do
1140 : 0 : call wrtout(std_out,ch10)
1141 : : end if
1142 : :
1143 [ - + ]: 476 : ABI_FREE(gshell)
1144 [ - + ]: 476 : ABI_FREE(insort)
1145 [ - + ]: 476 : ABI_FREE(nshell)
1146 : : #endif
1147 : :
1148 : 476 : call destroy_mpi_enreg(MPI_enreg_seq)
1149 [ - + ]: 476 : ABI_FREE(gbase)
1150 [ - + ]: 476 : ABI_FREE(shlim)
1151 [ - + ]: 476 : ABI_FREE(cnorm)
1152 [ - + ]: 476 : ABI_FREE(npwarr)
1153 : :
1154 [ - + - + : 952 : end subroutine merge_and_sort_kg
- + - + -
+ - + - +
- + - + -
+ - + ]
1155 : : !!***
1156 : :
1157 : : !----------------------------------------------------------------------
1158 : :
1159 : : !!****f* m_gsphere/getfullg
1160 : : !! NAME
1161 : : !! getfullg
1162 : : !!
1163 : : !! FUNCTION
1164 : : !! Reconstruct a G-sphere starting from a set of irreducible lattice vectors
1165 : : !!
1166 : : !! INPUTS
1167 : : !! pinv=-1 if time-reversal can be used, 1 otherwise
1168 : : !! nsym=number of symmetry operations
1169 : : !! sizepw=Max expected number of G vectors in the shere
1170 : : !! symrec(3,3,nsym)=symmetry operation in reciprocal space
1171 : : !! nbase=number of irreducible G vectors
1172 : : !! gbase(3,nbase)=irreducible G-vectors
1173 : : !! cnorm(nbase)=norm of the irreducible G vectors (supposed not yet sorted)
1174 : : !!
1175 : : !! OUTPUT
1176 : : !! maxpw=Number of G vectors found
1177 : : !! gbig(3,sizepw)=G vectors in the sphere packed in the first maxpw columns
1178 : : !! shlim(nbase)=number of G vectors within each shell
1179 : : !! ierr= Exit status, if /=0 the number of G vectors found exceeds sizepw
1180 : : !!
1181 : : !! SIDE EFFECTS
1182 : : !!
1183 : : !! NOTES
1184 : : !! cnorm is a bit redundant since it can be calculated from gbase. However this procedure
1185 : : !! is called by outkss in which cnorm is already calculated and we dont want to do it twice
1186 : : !!
1187 : : !! SOURCE
1188 : :
1189 : : subroutine getfullg(nbase,nsym,pinv,sizepw,gbase,symrec,cnorm,maxpw,gbig,shlim,ierr)
1190 : :
1191 : : !Arguments ------------------------------------
1192 : : !scalars
1193 : : integer,intent(in) :: nbase,nsym,pinv,sizepw
1194 : : integer,intent(out) :: ierr,maxpw
1195 : : !arrays
1196 : : integer,intent(in) :: gbase(3,nbase),symrec(3,3,nsym)
1197 : : integer,intent(out) :: gbig(3,sizepw),shlim(nbase)
1198 : : real(dp),intent(inout) :: cnorm(nbase) !sort_dp can change cnorm
1199 : :
1200 : : !Local variables-------------------------------
1201 : : !scalars
1202 : : integer :: ibase,ig,ilim,ish,isym,itim
1203 : : logical :: found
1204 : : character(len=500) :: msg
1205 : : !arrays
1206 : : integer :: gcur(3),geq(3)
1207 : : integer,allocatable :: gshell(:,:),insort(:),nshell(:)
1208 : : ! *************************************************************************
1209 : :
1210 : : if (pinv/=1.and.pinv/=-1) then
1211 : : write(msg,'(a,i6)')&
1212 : : & ' The argument pinv should be -1 or 1, however, pinv =',pinv
1213 : : ABI_BUG(msg)
1214 : : end if
1215 : : !
1216 : : ! === Reorder base g-vectors in order of increasing module ===
1217 : : ABI_MALLOC(insort,(nbase))
1218 : : do ibase=1,nbase
1219 : : insort(ibase)=ibase
1220 : : end do
1221 : : call sort_dp(nbase,cnorm,insort,tol14)
1222 : : !
1223 : : ! === Generate all stars of G-vectors ===
1224 : : ! Star of G is the set of all symmetrical images of the vector
1225 : : ! gshell contains the symmetrical G at fixed gbase. No need to add an additional dimension
1226 : : ! or initialize to zero the array inside the loop over nbase as we loop over (ish<=nshell(ibase))
1227 : : ABI_MALLOC(nshell,(nbase))
1228 : : ABI_MALLOC(gshell,(3,2*nsym))
1229 : : !
1230 : : ! === Start with zero number of G vectors found ===
1231 : : maxpw=0 ; ierr=0
1232 : : do ibase=1,nbase
1233 : : !
1234 : : ! === Loop over all different modules of G ===
1235 : : ! * Start with zero G vectors found in this star
1236 : : nshell(ibase)=0
1237 : : gcur(:)=gbase(:,insort(ibase))
1238 : : !
1239 : : ! === Loop over symmetries ===
1240 : : do isym=1,nsym
1241 : : do itim=pinv,1,2
1242 : : geq(:)=itim*MATMUL(symrec(:,:,isym),gcur)
1243 : : !
1244 : : ! * Search for symmetric of g and eventually add it:
1245 : : found=.FALSE. ; ish=1
1246 : : do while ((.not.found).and. (ish<=nshell(ibase)))
1247 : : found=ALL(geq(:)==gshell(:,ish))
1248 : : ish=ish+1
1249 : : end do
1250 : : if (.not.found) then
1251 : : nshell(ibase)=nshell(ibase)+1
1252 : : gshell(:,nshell(ibase))=geq(:)
1253 : : end if
1254 : : end do
1255 : : end do
1256 : : !
1257 : : ! * Was sizepw large enough?
1258 : : if ((maxpw+nshell(ibase))>sizepw) then
1259 : : write(msg,'(a,i6,2a)')&
1260 : : & ' Number of G in sphere exceeds maximum allowed value =',sizepw,ch10,&
1261 : : & ' check the value of sizepw in calling routine '
1262 : : ABI_WARNING(msg)
1263 : : ierr=1; RETURN
1264 : : end if
1265 : : !
1266 : : ! === Store this shell of Gs in a big array (gbig) ===
1267 : : do ig=1,nshell(ibase)
1268 : : gbig(:,ig+maxpw)=gshell(:,ig)
1269 : : end do
1270 : : maxpw=maxpw+nshell(ibase)
1271 : : end do ! ibase
1272 : : !
1273 : : ! === Compute number of G"s within each shell ===
1274 : : ilim=0
1275 : : do ibase=1,nbase
1276 : : ilim=ilim+nshell(ibase)
1277 : : shlim(ibase)=ilim
1278 : : end do
1279 : : !
1280 : : ! === Print out shell limits ===
1281 : : write(msg,'(3a)')&
1282 : : & ' Shells found:',ch10,&
1283 : : & ' number of shell number of G vectors cut-off energy [Ha] '
1284 : : call wrtout(std_out,msg)
1285 : :
1286 : : do ibase=1,nbase
1287 : : write(msg,'(12x,i4,17x,i6,12x,f8.3)')ibase,shlim(ibase),two*pi**2*cnorm(ibase)
1288 : : call wrtout(std_out,msg)
1289 : : end do
1290 : : write(msg,'(a)')ch10
1291 : : call wrtout(std_out,msg)
1292 : : ABI_FREE(gshell)
1293 : : ABI_FREE(insort)
1294 : : ABI_FREE(nshell)
1295 : :
1296 : : end subroutine getfullg
1297 : : !!***
1298 : :
1299 : : !----------------------------------------------------------------------
1300 : :
1301 : : !!****f* m_gsphere/get_irredg
1302 : : !! NAME
1303 : : !! get_irredg
1304 : : !!
1305 : : !! FUNCTION
1306 : : !! Given a set of reciprocal lattice vectors, find the set of G"s generating the others by symmetry.
1307 : : !!
1308 : : !! INPUTS
1309 : : !! nsym=number of symmetry operations
1310 : : !! pinv=-1 if time-reversal can be used, 1 otherwise
1311 : : !! npw_k=number of G vectors (for this k-point, as the set of G is k-centered)
1312 : : !! gcurr(3,npw_k)=the list of G vectors
1313 : : !! gprimd(3,3)=dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
1314 : : !! symrec(3,3,nsym)=symmetry operations in terms of reciprocal space primitive translations.
1315 : : !!
1316 : : !! OUTPUT
1317 : : !! nbasek=number of irreducible G vectors found
1318 : : !! cnormk(npw_k)=first nbasek elements are the norm of each irreducible G-vector
1319 : : !! gbasek(3,npw_k)=first nbasek elements are the irreducible G vectors
1320 : : !!
1321 : : !! NOTES
1322 : : !! The search can be optimized by looping over shells. See m_skw for a faster algo
1323 : : !!
1324 : : !! SOURCE
1325 : :
1326 : 3727 : subroutine get_irredg(npw_k,nsym,pinv,gprimd,symrec,gcurr,nbasek,gbasek,cnormk)
1327 : :
1328 : : !Arguments ------------------------------------
1329 : : !scalars
1330 : : integer,intent(in) :: npw_k,nsym,pinv
1331 : : integer,intent(out) :: nbasek
1332 : : !arrays
1333 : : integer,intent(in) :: gcurr(3,npw_k),symrec(3,3,nsym)
1334 : : integer,intent(out) :: gbasek(3,npw_k)
1335 : : real(dp),intent(in) :: gprimd(3,3)
1336 : : real(dp),intent(out) :: cnormk(npw_k)
1337 : :
1338 : : !Local variables-------------------------------
1339 : : !scalars
1340 : : integer :: ig,irr,isym,jj
1341 : : real(dp) :: eps,norm
1342 : : logical :: found
1343 : : !arrays
1344 : : integer :: gbas(3),gcur(3),geq(3)
1345 : : real(dp) :: gcar(3)
1346 : : ! *************************************************************************
1347 : :
1348 : : DBG_ENTER("COLL")
1349 : :
1350 [ - + ]: 3727 : if (pinv/=1.and.pinv/=-1) then
1351 : 0 : ABI_BUG(sjoin('pinv should be -1 or 1, however, pinv =', itoa(pinv)))
1352 : : end if
1353 : :
1354 : : ! Zero irred G vectors found, zeroing output arrays.
1355 [ + + + + : 4613627 : nbasek = 0; cnormk(:) = zero; gbasek(:,:) = 0
+ + ]
1356 : :
1357 [ + + ]: 925707 : do ig=1,npw_k
1358 [ + + ]: 3687920 : gcur(:) = gcurr(:,ig); norm = zero
1359 [ + + ]: 3687920 : do jj=1,3
1360 : 2765940 : gcar(jj)=gcur(1)*gprimd(jj,1)+gcur(2)*gprimd(jj,2)+gcur(3)*gprimd(jj,3)
1361 : 3687920 : norm=norm+gcar(jj)**2
1362 : : end do
1363 : 921980 : eps = tol8 * norm; found = .False.; irr = 1
1364 [ + + + + ]: 113366921 : do while (.not.found .and. irr <= nbasek) ! This loop can be optimized by looping inside the shell.
1365 [ + + ]: 112444941 : if (abs(norm - cnormk(irr)) <= eps) then
1366 [ + + ]: 7313164 : gbas(:) = gbasek(:,irr); isym = 1
1367 [ + + + + ]: 31630910 : do while (.not.found .and. isym <= nsym)
1368 [ + + + + ]: 387434047 : geq(:) = matmul(symrec(:,:,isym),gcur)
1369 [ + + + + ]: 35004994 : found = all(geq(:) == gbas(:))
1370 [ - + - - : 29802619 : if (pinv == -1) found = (found .or. all(geq == -gbas)) ! For time-reversal
- - ]
1371 : 31630910 : isym = isym + 1
1372 : : end do
1373 : : end if
1374 : 112523342 : irr = irr + 1
1375 : : end do
1376 [ + + ]: 925707 : if (.not. found) then
1377 [ + + ]: 313604 : nbasek = nbasek + 1; cnormk(nbasek) = norm; gbasek(:,nbasek) = gcur(:)
1378 : : end if
1379 : : end do
1380 : :
1381 : : DBG_EXIT("COLL")
1382 : :
1383 : 3727 : end subroutine get_irredg
1384 : : !!***
1385 : :
1386 : : !----------------------------------------------------------------------
1387 : :
1388 : : !!****f* m_gsphere/merge_kgirr
1389 : : !! NAME
1390 : : !! merge_kgirr
1391 : : !!
1392 : : !! FUNCTION
1393 : : !! Given a list of irreducible reciprocal vectors associated to different k-centered spheres,
1394 : : !! this subroutine finds the minimal set of G vectors needed to reconstruct the union of the spheres
1395 : : !! through symmetry operations.
1396 : : !!
1397 : : !! INPUTS
1398 : : !! nsym=number of symmetry operations
1399 : : !! pinv=-1 if time-reversal can be used, 0 otherwise
1400 : : !! nkpt=number of k-points for k-centered spheres
1401 : : !! mpw=Max number of G vectors for each k-point
1402 : : !! sizepw=Max expected number of G vectors found
1403 : : !! symrec(3,3,nsym)=symmetries in reciprocal space given in reduced coordinates
1404 : : !! nbasek(nkpt)=number of irred G for each k-point
1405 : : !! cnormk(mpw,nkpt)=the norm of each k-centered G (only 1:nbase(ik)) is used
1406 : : !! gbasek(3,mpw,nkpt)
1407 : : !!
1408 : : !! OUTPUT
1409 : : !! nbase=number of irreducible G needed to reconstruct the initial set of spheres
1410 : : !! gbase(3,sizepw)=irreducible G found in reciprocal coordinates
1411 : : !! cnorm(sizepw)=Norm of each irred G vector
1412 : : !! ierr= Exit status, if /=0 the number of G vectors found exceeds sizepw
1413 : : !!
1414 : : !! SOURCE
1415 : :
1416 : 476 : subroutine merge_kgirr(nsym,pinv,nkpt,mpw,sizepw,symrec,nbasek,cnormk,gbasek,nbase,gbase,cnorm,ierr)
1417 : :
1418 : : !Arguments ------------------------------------
1419 : : !scalars
1420 : : integer,intent(in) :: mpw,nkpt,nsym,pinv,sizepw
1421 : : integer,intent(out) :: ierr,nbase
1422 : : !arrays
1423 : : integer,intent(in) :: gbasek(3,mpw,nkpt),nbasek(nkpt),symrec(3,3,nsym)
1424 : : integer,intent(inout) :: gbase(3,sizepw) !vz_i
1425 : : real(dp),intent(in) :: cnormk(mpw,nkpt)
1426 : : real(dp),intent(inout) :: cnorm(sizepw) !vz_i
1427 : :
1428 : : !Local variables-------------------------------
1429 : : !scalars
1430 : : integer :: ikpt,inb,irgk,isym
1431 : : real(dp) :: eps,norm
1432 : : logical :: found
1433 : : character(len=500) :: msg
1434 : : !arrays
1435 : : integer :: gbas(3),gcur(3),geq(3)
1436 : : ! *************************************************************************
1437 : :
1438 : : DBG_ENTER("COLL")
1439 : :
1440 [ - + ]: 476 : if (pinv/=1.and.pinv/=-1) then
1441 : 0 : write(msg,'(a,i6)')' The argument pinv should be -1 or 1, however, pinv =',pinv
1442 : 0 : ABI_BUG(msg)
1443 : : end if
1444 : : !
1445 : : ! === Start with zero number of G found ===
1446 : 476 : nbase=0 ; ierr=0
1447 [ + + ]: 2749 : do ikpt=1,nkpt
1448 [ + + ]: 78702 : do irgk=1,nbasek(ikpt)
1449 [ + + ]: 303812 : gcur(:)=gbasek(:,irgk,ikpt)
1450 : 75953 : norm=cnormk(irgk,ikpt) ; eps=tol8*norm
1451 : 75953 : found=.FALSE. ; inb=1
1452 [ + + + + ]: 9608891 : do while ((.not.found).and.(inb<=nbase))
1453 [ + + ]: 9532938 : if (ABS(norm-cnorm(inb))<=eps) then
1454 [ + + ]: 513072 : gbas(:)=gbase(:,inb)
1455 : : isym=1
1456 [ + + + + ]: 1622889 : do while ((.not.found).and.(isym<=nsym))
1457 [ + + + + ]: 19430073 : geq(:)=MATMUL(symrec(:,:,isym),gcur)
1458 [ + + + + ]: 1729181 : found=ALL(geq(:)==gbas(:))
1459 [ - + - - : 1494621 : if (pinv==-1) found= (found.or.ALL(geq(:)==-gbas(:)) ) ! For time-reversal
- - ]
1460 : 1622889 : isym=isym+1
1461 : : end do
1462 : : end if
1463 : 9570057 : inb=inb+1
1464 : : end do
1465 [ + + ]: 78226 : if (.not.found) then
1466 : : ! === Add to the list ===
1467 : 37119 : nbase=nbase+1
1468 [ - + ]: 37119 : if (nbase>sizepw) then
1469 : : write(msg,'(2(a,i5),a)')&
1470 : 0 : & ' nbase (',nbase,') became greater than sizepw = ',sizepw,' returning ierr=1 '
1471 : 0 : ABI_WARNING(msg)
1472 : 0 : ierr=1; RETURN
1473 : : end if
1474 : 37119 : cnorm(nbase)=cnormk(irgk,ikpt)
1475 [ + + ]: 148476 : gbase(:,nbase)=gcur(:)
1476 : : end if
1477 : : end do
1478 : : end do
1479 : :
1480 : : DBG_EXIT("COLL")
1481 : :
1482 : : end subroutine merge_kgirr
1483 : : !!***
1484 : :
1485 : : !----------------------------------------------------------------------
1486 : :
1487 : : !!****f* m_gpshere/setshells
1488 : : !! NAME
1489 : : !! setshells
1490 : : !!
1491 : : !! FUNCTION
1492 : : !! Set consistently the number of shells, the number of plane-waves, and the energy cut-off
1493 : : !!
1494 : : !! INPUTS
1495 : : !! nsym=number of symmetry operations
1496 : : !! gmet(3,3)=metric tensor in reciprocal space
1497 : : !! gprimd(3,3)=dimensional primitive vectors in reciprocal space
1498 : : !! symrel(3,3,nsym)=symmetry operations in real space
1499 : : !! tag=suffix to account for the different possibilities for these variables (npw, ecut or nsh ..)
1500 : : !! ucvol=unit cell volume
1501 : : !!
1502 : : !! OUTPUT
1503 : : !! (see side effects)
1504 : : !!
1505 : : !! SIDE EFFECTS
1506 : : !! ecut,npw,nsh=one of them is an input, the two others are output
1507 : : !! ecut=cut-off energy for plane wave basis sphere (Ha)
1508 : : !! npw=number of plane waves
1509 : : !! nsh=number of shells
1510 : : !!
1511 : : !! SOURCE
1512 : :
1513 : 1419 : subroutine setshells(ecut,npw,nsh,nsym,gmet,gprimd,symrel,tag,ucvol)
1514 : :
1515 : : !Arguments ------------------------------------
1516 : : !scalars
1517 : : integer,intent(in) :: nsym
1518 : : integer,intent(inout) :: npw,nsh
1519 : : real(dp),intent(in) :: ucvol
1520 : : real(dp),intent(inout) :: ecut
1521 : : character(len=*),intent(in) :: tag
1522 : : !arrays
1523 : : integer,intent(in) :: symrel(3,3,nsym)
1524 : : real(dp),intent(in) :: gmet(3,3),gprimd(3,3)
1525 : :
1526 : : !Local variables-------------------------------
1527 : : !scalars
1528 : : integer :: exchn2n3d,ifound,ig,ii,ish,isym,npw_found,npwave
1529 : : integer :: npwwrk,nsh_found,pad=50
1530 : : real(dp) :: ecut_found,ecut_trial,eps,scale=1.3_dp
1531 : : logical :: found
1532 : : character(len=500) :: msg
1533 : 1419 : type(MPI_type) :: MPI_enreg_seq
1534 : : !arrays
1535 : : integer :: geq(3)
1536 : 1419 : integer,allocatable :: gvec(:,:),gvec_sh(:,:),insort(:),npw_sh(:)
1537 : : real(dp) :: gctr(3)
1538 : 1419 : real(dp),allocatable :: gnorm(:),gnorm_sh(:)
1539 : :
1540 : : !******************************************************************
1541 : :
1542 : : ! Check coherence of input variables ecut, npw, and nsh.
1543 : : ! 1-> one at least should be non-null
1544 [ + + + - : 1419 : if (npw==0.and.nsh==0.and.ecut<=tol6) then
+ + ]
1545 : : write(msg,'(8a)')&
1546 : 181 : 'One of the three variables ecut',TRIM(tag),', npw',TRIM(tag),', or nsh',TRIM(tag),ch10,&
1547 : 362 : 'must be non-null. Returning.'
1548 : 181 : ABI_COMMENT(msg)
1549 : 181 : RETURN
1550 : : end if
1551 : : ! 2-> one and only one should be non-null
1552 [ + + - + ]: 1238 : if (npw/=0.and.nsh/=0) then
1553 : : write(msg,'(6a)')&
1554 : 0 : 'Only one of the two variables npw',TRIM(tag),' and nsh',TRIM(tag),ch10,&
1555 : 0 : 'can be non-null. Modify the value of one of these in input file.'
1556 : 0 : ABI_ERROR(msg)
1557 : : end if
1558 [ + + - + ]: 1238 : if (ecut>tol6.and.npw/=0) then
1559 : : write(msg,'(6a)')&
1560 : 0 : 'Only one of the two variables ecut',TRIM(tag),' and npw',TRIM(tag),ch10,&
1561 : 0 : 'can be non-null. Modify the value of one of these in input file.'
1562 : 0 : ABI_ERROR(msg)
1563 : : end if
1564 [ + + - + ]: 1238 : if (ecut>tol6.and.nsh/=0) then
1565 : : write(msg,'(6a)')&
1566 : 0 : 'Only one of the two variables ecut',TRIM(tag),' and nsh',TRIM(tag),ch10,&
1567 : 0 : 'can be non-null Action : modify the value of one of these in input file.'
1568 : 0 : ABI_ERROR(msg)
1569 : : end if
1570 : :
1571 : : ! Calculate an upper bound for npw.
1572 : : ! gctr is center of the g-vector sphere
1573 : 1238 : gctr(:)= [zero,zero,zero]
1574 [ + + ]: 1238 : if (ecut>tol6) then
1575 : : ! The average number of plane-waves in the cutoff sphere is given by:
1576 : : ! npwave = (2*ecut)**(3/2)*ucvol/(6*pi**2)
1577 : : ! The upper bound is calculated as npwwrk=int(scale * npwave) + pad
1578 : 1110 : npwave=NINT(ucvol*(two*ecut)**1.5_dp/(six*pi**2))
1579 : 1110 : npwwrk=NINT(DBLE(npwave)*scale)+pad
1580 : 1110 : ecut_trial=ecut
1581 [ + - ]: 128 : else if (npw/=0) then
1582 : : ! npw is given in the input
1583 : 128 : npwwrk=NINT(DBLE(npw)*scale)+pad
1584 : 128 : ecut_trial=(six*pi**2*npw/ucvol)**two_thirds/two
1585 : : else
1586 : : ! If nsh is given in the input
1587 : 0 : npwwrk=nsh*18+2*pad
1588 : 0 : ecut_trial=(six*pi**2*nsh*18/ucvol)**two_thirds/two
1589 : : end if
1590 : :
1591 : 1238 : call initmpi_seq(MPI_enreg_seq)
1592 : :
1593 [ - + + - : 3714 : ABI_MALLOC(gvec,(3,npwwrk))
+ - + - -
+ - + ]
1594 : 2476 : ifound=0
1595 [ + + ]: 2476 : do while (ifound==0)
1596 : : !write(msg,'(a,f8.2)')' setshells : ecut_trial = ',ecut_trial
1597 : : !call wrtout(std_out,msg,'COLL')
1598 : 1238 : exchn2n3d=0 ! For the time being, no exchange of n2 and n3
1599 : :
1600 : 1238 : call kpgsph(ecut_trial,exchn2n3d,gmet,0,1,1,gvec,gctr,1,MPI_enreg_seq,npwwrk,npw_found)
1601 : :
1602 [ - + + - : 3714 : ABI_MALLOC(gnorm,(npw_found))
+ - + - -
+ - + -
+ ]
1603 [ - + + - : 3714 : ABI_MALLOC(insort,(npw_found))
+ - + - -
+ - + -
+ ]
1604 : :
1605 [ + + ]: 361732 : do ig=1,npw_found
1606 : 360494 : insort(ig)=ig
1607 : 360494 : gnorm(ig)=zero
1608 [ + + ]: 1443214 : do ii=1,3
1609 : : gnorm(ig)=gnorm(ig)+(gvec(1,ig)*gprimd(ii,1)+&
1610 : : gvec(2,ig)*gprimd(ii,2)+&
1611 : 1441976 : gvec(3,ig)*gprimd(ii,3))**2
1612 : : end do
1613 : : end do
1614 : 1238 : call sort_dp(npw_found,gnorm,insort,tol14)
1615 : :
1616 [ - + + - : 2476 : ABI_MALLOC(npw_sh,(npw_found))
- + - + -
+ ]
1617 [ - + + - : 2476 : ABI_MALLOC(gnorm_sh,(npw_found))
- + - + -
+ ]
1618 [ - + + - : 3714 : ABI_MALLOC(gvec_sh,(3,npw_found))
+ - + - -
+ - + -
+ ]
1619 [ + + ]: 361732 : npw_sh(:)=0
1620 [ + + ]: 361732 : gnorm_sh(:)=zero
1621 [ + + + + ]: 1443214 : gvec_sh(:,:)=0
1622 : : ! Count the number of shells:
1623 : : ! (search for the G-vectors generating the others by symmetry)
1624 : 1238 : nsh_found=0
1625 : :
1626 [ + + ]: 361732 : do ig=1,npw_found
1627 : 360494 : eps=1.d-8*gnorm(ig)
1628 : 360494 : found=.FALSE.
1629 : 360494 : ish=1
1630 [ + + ]: 43121633 : do while ((.not.found).and.(ish<=nsh_found))
1631 [ + + ]: 42761139 : if (ABS(gnorm(ig)-gnorm_sh(ish))<=eps) then
1632 : : isym=1
1633 [ + + + + ]: 11918774 : do while ((.not.found).and.(isym<=nsym))
1634 : : geq(:)=(symrel(1,:,isym)*gvec(1,insort(ig))+&
1635 : : symrel(2,:,isym)*gvec(2,insort(ig))+&
1636 [ + + ]: 44537276 : symrel(3,:,isym)*gvec(3,insort(ig)))
1637 : :
1638 : : found=((geq(1)==gvec_sh(1,ish)).and.&
1639 : : (geq(2)==gvec_sh(2,ish)).and.&
1640 [ + + + + : 11134319 : (geq(3)==gvec_sh(3,ish)))
+ + ]
1641 : 11463703 : isym=isym+1
1642 : : end do
1643 : : end if
1644 : 42761139 : ish=ish+1
1645 : : end do
1646 [ + + ]: 361732 : if (.not.found) then
1647 : 31110 : nsh_found=nsh_found+1
1648 : 31110 : gnorm_sh(nsh_found)=gnorm(ig)
1649 [ + + ]: 124440 : gvec_sh(:,nsh_found)=gvec(:,insort(ig))
1650 : 31110 : npw_sh(nsh_found)=1
1651 : : else
1652 : 329384 : ish=ish-1
1653 : 329384 : npw_sh(ish)=npw_sh(ish)+1
1654 : : end if
1655 : : end do
1656 : :
1657 : 1238 : ecut_found=two*pi**2*gnorm(npw_found)
1658 : :
1659 [ + + ]: 1238 : if(ecut>tol6) then
1660 : : ! ecut is given in the input
1661 : : !if (ecut_found<ecut-0.1) then
1662 : : ! write(msg,'(3a,e14.6,9a,e14.6,3a)')&
1663 : : ! 'The value ecut',TRIM(tag),'=',ecut,' given in the input file leads to',ch10,&
1664 : : ! 'the same values for nsh',TRIM(tag),' and npw',TRIM(tag),' as ecut',TRIM(tag),'=',ecut_found,ch10
1665 : : ! ABI_COMMENT(msg)
1666 : : !end if
1667 : : ifound=1
1668 [ + - ]: 128 : else if (npw/=0) then
1669 : : ! If npw is given in the input
1670 [ + + ]: 128 : if (npw_found==npw) then
1671 : : ecut_found=two*pi**2*gnorm(npw_found)
1672 : : ifound=1
1673 [ + - ]: 11 : else if (npw_found>npw) then
1674 : 11 : npw_found=0
1675 : 11 : nsh_found=0
1676 [ + + ]: 60 : do while (npw_found<npw)
1677 : 49 : nsh_found=nsh_found+1
1678 : 49 : npw_found=npw_found+npw_sh(nsh_found)
1679 : : end do
1680 : : ! check that the shell is closed
1681 [ - + ]: 11 : if(npw_found>npw) then
1682 : : ! shell not closed
1683 : 0 : npw_found=npw_found-npw_sh(nsh_found)
1684 : 0 : nsh_found=nsh_found-1
1685 [ # # ]: 0 : do while (ABS(gnorm_sh(nsh_found)-gnorm_sh(nsh_found+1))<0.000001)
1686 : 0 : npw_found=npw_found-npw_sh(nsh_found)
1687 : 0 : nsh_found=nsh_found-1
1688 : : end do
1689 : : write(msg,'(3a,i6,5a,i6,3a)')&
1690 : 0 : 'The value npw',TRIM(tag),'=',npw,' given in the input file does not close the shell',ch10,&
1691 : 0 : 'The lower closed-shell is obtained for a value npw',TRIM(tag),'=',npw_found,ch10,&
1692 : 0 : 'This value will be adopted for the calculation.',ch10
1693 : 0 : ABI_WARNING(msg)
1694 : : end if
1695 : 11 : ecut_found=two*pi**2*gnorm(npw_found)
1696 : : ifound=1
1697 : : end if
1698 [ # # ]: 0 : else if (nsh/=0) then
1699 : : ! If nsh is given in the input
1700 [ # # ]: 0 : if (nsh_found==nsh) then
1701 : : ecut_found=two*pi**2*gnorm(npw_found)
1702 : : ifound=1
1703 [ # # ]: 0 : else if (nsh_found>nsh) then
1704 : 0 : npw_found=0
1705 : 0 : nsh_found=0
1706 [ # # ]: 0 : do ish=1,nsh
1707 : 0 : npw_found=npw_found+npw_sh(ish)
1708 : 0 : nsh_found=nsh_found+1
1709 : : end do
1710 [ # # ]: 0 : if (ABS(gnorm_sh(nsh_found)-gnorm_sh(nsh_found+1))<0.000001) then
1711 [ # # ]: 0 : do while (ABS(gnorm_sh(nsh_found)-gnorm_sh(nsh_found+1))<0.000001)
1712 : 0 : nsh_found=nsh_found+1
1713 : 0 : npw_found=npw_found+npw_sh(nsh_found)
1714 : : end do
1715 : : write(msg,'(3a,i6,5a,i6,3a)')&
1716 : 0 : 'The value nsh',TRIM(tag),'=',nsh,' given in the input file corresponds to the same',ch10,&
1717 : 0 : 'cut-off energy as for closed-shell upto nsh',TRIM(tag),'=',nsh_found,ch10,&
1718 : 0 : 'This value will be adopted for the calculation.',ch10
1719 : 0 : ABI_WARNING(msg)
1720 : : end if
1721 : 0 : ecut_found=two*pi**2*gnorm(npw_found)
1722 : : ifound=1
1723 : : end if
1724 : : end if
1725 : :
1726 : 1238 : if (ifound==0) then
1727 : 0 : ecut_trial=1.1*ecut_trial
1728 : 0 : ABI_FREE(gnorm)
1729 [ # # ]: 0 : ABI_FREE(gnorm_sh)
1730 [ # # ]: 0 : ABI_FREE(gvec_sh)
1731 [ # # ]: 0 : ABI_FREE(insort)
1732 [ # # ]: 0 : ABI_FREE(npw_sh)
1733 : : else
1734 : : ! ecut was not provided as an input, then set it now!
1735 [ + + ]: 1238 : if (ecut<tol6) then
1736 : 128 : ecut=ecut_found
1737 : : end if
1738 : 1238 : npw=npw_found
1739 : 1238 : nsh=nsh_found
1740 : : end if
1741 : :
1742 : : end do ! while(ifound==0)
1743 : :
1744 : 1238 : call destroy_mpi_enreg(MPI_enreg_seq)
1745 : :
1746 [ - + ]: 1238 : ABI_FREE(gnorm)
1747 [ - + ]: 1238 : ABI_FREE(gnorm_sh)
1748 [ - + ]: 1238 : ABI_FREE(gvec)
1749 [ - + ]: 1238 : ABI_FREE(gvec_sh)
1750 [ - + ]: 1238 : ABI_FREE(insort)
1751 [ - + ]: 1238 : ABI_FREE(npw_sh)
1752 : :
1753 [ - + - + : 1419 : end subroutine setshells
- + - + -
+ - + - +
- + - + -
+ - + ]
1754 : : !!***
1755 : :
1756 : : !----------------------------------------------------------------------
1757 : :
1758 : : !!****f* m_gsphere/kg_map
1759 : : !! NAME
1760 : : !! kg_map
1761 : : !!
1762 : : !! FUNCTION
1763 : : !! Compute the mapping between two lists of g-vectors.
1764 : : !!
1765 : : !! INPUTS
1766 : : !! npw1, kg1(3,npw1)=First list of G-vectors
1767 : : !! npw2, kg2(3,npw2)=Second list of G-vectors
1768 : : !!
1769 : : !! OUTPUT
1770 : : !! g2g1(npw2) = Mapping kg2 index --> kg1 index.
1771 : : !! Set to 0 if kg2(:,ig) not in kg1
1772 : : !! nmiss = Number of G-vectors in kg2 not found in kg1
1773 : : !!
1774 : : !! SOURCE
1775 : :
1776 : 4344 : subroutine kg_map(npw1, kg1, npw2, kg2, g2g1, nmiss)
1777 : :
1778 : : !Arguments ------------------------------------
1779 : : !scalars
1780 : : integer,intent(in) :: npw1,npw2
1781 : : integer,intent(out) :: nmiss
1782 : : !arrays
1783 : : integer,intent(in) :: kg1(3,npw1),kg2(3,npw2)
1784 : : integer,intent(out) :: g2g1(npw2)
1785 : :
1786 : : !Local variables ------------------------------
1787 : : !scalars
1788 : : integer :: ii,ipw,i1,i2,i3,n1,n2,n3
1789 : : !arrays
1790 : : integer :: gmax(3),g1_max(3),g2_max(3)
1791 : 4344 : integer,allocatable :: iwork(:,:,:)
1792 : : !************************************************************************
1793 : :
1794 [ + + + + : 7535520 : g1_max = maxval(abs(kg1))
+ + ]
1795 [ + + + + : 7103372 : g2_max = maxval(abs(kg2))
+ + ]
1796 [ + + ]: 17376 : do ii=1,3
1797 : 17376 : gmax(ii) = max(g1_max(ii), g2_max(ii))
1798 : : end do
1799 [ + + ]: 17376 : gmax = 2*gmax + 1
1800 : 4344 : n1 = gmax(1); n2 = gmax(2); n3 = gmax(3)
1801 : :
1802 : : !print *, "n1, n2, n3", n1, n2, n3; print *, "kg1:", kg1; print *, "kg2:", kg2
1803 [ - + + - : 21720 : ABI_MALLOC(iwork, (n1, n2, n3))
- + + - -
+ + - + -
+ - - + -
+ ]
1804 : :
1805 : : ! Insert kg1 into work with extra 0 s around outside:
1806 [ + + + + : 6722756 : iwork = 0
+ + ]
1807 [ + + ]: 1882794 : do ipw=1,npw1
1808 [ + + ]: 1878450 : i1 = kg1(1,ipw); if (i1<0) i1=i1+n1; i1=i1+1
1809 [ + + ]: 1878450 : i2 = kg1(2,ipw); if (i2<0) i2=i2+n2; i2=i2+1
1810 [ + + ]: 1878450 : i3 = kg1(3,ipw); if (i3<0) i3=i3+n3; i3=i3+1
1811 : 1882794 : iwork(i1,i2,i3) = ipw
1812 : : end do
1813 : :
1814 [ + + ]: 1774757 : g2g1 = 0; nmiss = 0
1815 [ + + ]: 1774757 : do ipw=1,npw2
1816 [ + + ]: 1770413 : i1 = kg2(1,ipw); if (i1<0) i1=i1+n1; i1=i1+1
1817 [ + + ]: 1770413 : i2 = kg2(2,ipw); if (i2<0) i2=i2+n2; i2=i2+1
1818 [ + + ]: 1770413 : i3 = kg2(3,ipw); if (i3<0) i3=i3+n3; i3=i3+1
1819 : : !#if defined FC_NVHPC
1820 : : ! !if (n1 == -1) write(std_out, *)"NVHPC raises SIGSEGV"
1821 : : ! print *, "ipw, i1, i2, i3, n1, n2, n3", ipw, i1, i2, i3, n1, n2, n3
1822 : : !#endif
1823 : 1770413 : g2g1(ipw) = iwork(i1,i2,i3)
1824 [ + + ]: 1774757 : if (g2g1(ipw) == 0) nmiss = nmiss + 1
1825 : : end do
1826 : :
1827 [ - + ]: 4344 : ABI_FREE(iwork)
1828 : :
1829 : 4344 : end subroutine kg_map
1830 : : !!***
1831 : :
1832 : : !----------------------------------------------------------------------
1833 : :
1834 : : !!****f* m_gsphere/make_istwk_table
1835 : : !! NAME
1836 : : !! make_istwfk_table
1837 : : !!
1838 : : !! FUNCTION
1839 : : !!
1840 : : !! INPUTS
1841 : : !! ng1,ng2,ng3
1842 : : !!
1843 : : !! OUTPUT
1844 : : !!
1845 : : !! NOTES
1846 : : !! Useful relations:
1847 : : !! u_k(G) = u_{k+G0}(G-G0); u_{-k}(G) = u_k(G)^*
1848 : : !! and therefore:
1849 : : !! u_{G0/2}(G) = u_{G0/2}(-G-G0)^*.
1850 : : !!
1851 : : !! SOURCE
1852 : :
1853 : 0 : subroutine make_istwfk_table(istwf_k,ng1,ng2,ng3,ig1_inver,ig2_inver,ig3_inver)
1854 : :
1855 : : !Arguments ------------------------------------
1856 : : !scalars
1857 : : integer,intent(in) :: ng1,ng2,ng3,istwf_k
1858 : : !arrays
1859 : : integer,intent(out) :: ig1_inver(ng1),ig2_inver(ng2),ig3_inver(ng3)
1860 : :
1861 : : !Local variables ------------------------------
1862 : : !scalars
1863 : : integer :: i1,i2,i3
1864 : : character(len=500) :: msg
1865 : : !************************************************************************
1866 : :
1867 : : ! Initialize the inverse coordinates
1868 : 0 : select case (istwf_k)
1869 : :
1870 : : case (1)
1871 : 0 : ig1_inver(1)=1
1872 [ # # ]: 0 : do i1=2,ng1
1873 : 0 : ig1_inver(i1)=ng1+2-i1
1874 : : end do
1875 : 0 : ig2_inver(1)=1
1876 [ # # ]: 0 : do i2=2,ng2
1877 : 0 : ig2_inver(i2)=ng2+2-i2
1878 : : end do
1879 : 0 : ig3_inver(1)=1
1880 [ # # ]: 0 : do i3=2,ng3
1881 : 0 : ig3_inver(i3)=ng3+2-i3
1882 : : end do
1883 : :
1884 : : case (2:8)
1885 [ # # ]: 0 : if (istwf_k==2 .or. istwf_k==4 .or. istwf_k==6 .or. istwf_k==8) then
1886 : 0 : ig1_inver(1)=1
1887 [ # # ]: 0 : do i1=2,ng1
1888 : 0 : ig1_inver(i1)=ng1+2-i1
1889 : : end do
1890 : : else
1891 [ # # ]: 0 : do i1=1,ng1
1892 : 0 : ig1_inver(i1)=ng1+1-i1
1893 : : end do
1894 : : end if
1895 [ # # ]: 0 : if (istwf_k>=2 .and. istwf_k<=5) then
1896 : 0 : ig2_inver(1)=1
1897 [ # # ]: 0 : do i2=2,ng2
1898 : 0 : ig2_inver(i2)=ng2+2-i2
1899 : : end do
1900 : : else
1901 [ # # ]: 0 : do i2=1,ng2
1902 : 0 : ig2_inver(i2)=ng2+1-i2
1903 : : end do
1904 : : end if
1905 [ # # # # ]: 0 : if (istwf_k==2 .or. istwf_k==3 .or. istwf_k==6 .or. istwf_k==7) then
1906 : 0 : ig3_inver(1)=1
1907 [ # # ]: 0 : do i3=2,ng3
1908 : 0 : ig3_inver(i3)=ng3+2-i3
1909 : : end do
1910 : : else
1911 [ # # ]: 0 : do i3=1,ng3
1912 : 0 : ig3_inver(i3)=ng3+1-i3
1913 : : end do
1914 : : end if
1915 : :
1916 : : case default
1917 : 0 : write(msg,'(a,i0)')" Wrong value for istwf_k: ",istwf_k
1918 [ # # # ]: 0 : ABI_ERROR(msg)
1919 : : end select
1920 : :
1921 : 0 : end subroutine make_istwfk_table
1922 : : !!***
1923 : :
1924 : : !----------------------------------------------------------------------
1925 : :
1926 : : !!****f* m_gsphere/table_gbig2kg
1927 : : !! NAME
1928 : : !! table_gbig2kg
1929 : : !!
1930 : : !! FUNCTION
1931 : : !! Associate the kg_k set of g-vectors with the big array of gbig
1932 : : !! The array gbig(3,maxpw) contains all g-vectors used for all k-points, in order of
1933 : : !! increasing shells. For a each k-point, the wave-functions are defined only on a particular set
1934 : : !! of g-vectors kg_k (included in gbig). This set is defined by array gamma2k:
1935 : : !! The array gamma2k(ig=1,maxpw) translates the index of the gbig (from 1 to maxpw) into the corresponding
1936 : : !! index in array kg_k. If gbig(ig) does not exist in kg_k, gamma2k(ig) contains npw_k+1.
1937 : : !!
1938 : : !! INPUTS
1939 : : !! npw_k=Number of planewaves in the k-centered basis set
1940 : : !! kg_k(3,npw_k)=The k-centered basis set
1941 : : !! maxpw=Number of G in gbig
1942 : : !! gbig(3,maxpw)=The union of the G-spheres at different k-points.
1943 : : !!
1944 : : !! OUTPUT
1945 : : !! ierr=Status error. It gives the number of G of kg_k not contained in gbig.
1946 : : !! gamma2k(maxpw)=Mapping gbig -> kg_k
1947 : : !!
1948 : : !! SOURCE
1949 : :
1950 : 31 : pure subroutine table_gbig2kg(npw_k,kg_k,maxpw,gbig,gamma2k,ierr)
1951 : :
1952 : : !Arguments ------------------------------------
1953 : : !scalars
1954 : : integer,intent(in) :: npw_k,maxpw
1955 : : integer,intent(out) :: ierr
1956 : : !arrays
1957 : : integer,intent(in) :: kg_k(3,npw_k)
1958 : : integer,intent(in) :: gbig(3,maxpw)
1959 : : integer,intent(out) :: gamma2k(maxpw)
1960 : :
1961 : : !Local variables-------------------------------
1962 : : !scalars
1963 : : integer :: ig,igp
1964 : : logical :: found
1965 : : !arrays
1966 : : integer :: gcur(3)
1967 : : ! *********************************************************************
1968 : :
1969 : 31 : ierr=0
1970 [ + + ]: 30608 : gamma2k(:)=npw_k+1 ! Initialize array gamma2k
1971 : :
1972 [ + + ]: 21877 : do ig=1,npw_k ! Loop over g-vectors, for this k point.
1973 [ + + ]: 87384 : gcur(:)=kg_k(:,ig)
1974 : : igp=0; found=.FALSE.
1975 [ + + - + ]: 11381569 : do while ((.not.found) .and. igp<maxpw) ! Search selected vector in array gbig: TODO this part can be optimized
1976 : 11359723 : igp=igp+1
1977 [ + + + + ]: 12455395 : found=ALL(gcur(:)==gbig(:,igp))
1978 : : end do
1979 [ + - ]: 21877 : if (found) then ! Store it if found:
1980 : 21846 : gamma2k(igp)=ig
1981 : : else
1982 : 0 : ierr=ierr+1
1983 : : end if
1984 : : end do
1985 : :
1986 : 31 : end subroutine table_gbig2kg
1987 : : !!***
1988 : :
1989 : : !----------------------------------------------------------------------
1990 : :
1991 : : !!****f* m_gsphere/gsph_extend
1992 : : !! NAME
1993 : : !! gsph_extend
1994 : : !!
1995 : : !! FUNCTION
1996 : : !! Construct a new gsphere_t with a larger cutoff energy
1997 : : !! while preserving the ordering of the first G-vectors stored in in_Gsph
1998 : : !!
1999 : : !! INPUTS
2000 : : !!
2001 : : !! OUTPUT
2002 : : !!
2003 : : !! SOURCE
2004 : :
2005 : 36 : subroutine gsph_extend(in_Gsph, Cryst, new_ecut, new_Gsph)
2006 : :
2007 : : !Arguments ------------------------------------
2008 : : !scalars
2009 : : class(gsphere_t),intent(in) :: in_Gsph
2010 : : type(crystal_t),intent(in) :: Cryst
2011 : : real(dp),intent(in) :: new_ecut
2012 : : class(gsphere_t),intent(out) :: new_Gsph
2013 : :
2014 : : !Local variables-------------------------------
2015 : : !scalars
2016 : : integer :: new_ng,in_ng,ig,ierr,sh
2017 : : !arrays
2018 [ + - ]: 36 : integer,allocatable :: new_gvec(:,:)
2019 : : ! *********************************************************************
2020 : :
2021 : 36 : call new_Gsph%init(Cryst, 0, ecut=new_ecut)
2022 : :
2023 [ + + ]: 36 : if (new_Gsph%ng > in_Gsph%ng) then
2024 : : ! new_gpsh larger than in_gsph
2025 : 34 : new_ng = new_Gsph%ng
2026 : 34 : in_ng = in_Gsph%ng
2027 : :
2028 : 34 : ierr = 0
2029 [ + + ]: 1800 : do ig=1,in_ng
2030 [ + + + - : 7098 : if (ANY(new_Gsph%gvec(:,ig) /= in_Gsph%gvec(:,ig)) ) then
- + ]
2031 : 0 : ierr = ierr + 1
2032 : 0 : write(std_out,*)" new_gvec, in_gvec",ig,new_Gsph%gvec(:,ig),in_Gsph%gvec(:,ig)
2033 : : end if
2034 : : end do
2035 : :
2036 [ + - ]: 34 : if (ierr == 0) RETURN
2037 : :
2038 : 0 : ierr = 0
2039 [ # # ]: 0 : do sh=1,in_Gsph%nsh
2040 [ # # # # ]: 0 : if (new_Gsph%shlim(sh) /= in_Gsph%shlim(sh) .or. &
2041 : 0 : ABS(new_Gsph%shlen(sh)-in_Gsph%shlen(sh)) > tol12 ) then
2042 : 0 : ierr = ierr + 1
2043 : 0 : write(std_out,*)"new_shlim, in_shlim",sh,new_Gsph%shlim(sh),in_Gsph%shlim(sh)
2044 : 0 : write(std_out,*)"new_shlen, in_shlen",sh,new_Gsph%shlen(sh),in_Gsph%shlen(sh)
2045 : : end if
2046 : : end do
2047 [ # # ]: 0 : ABI_CHECK(ierr == 0,"Wrong shells")
2048 : :
2049 [ # # # # : 0 : ABI_MALLOC(new_gvec,(3, new_ng))
# # # # #
# # # ]
2050 [ # # # # : 0 : new_gvec = new_Gsph%gvec
# # # # #
# # # # #
# # ]
2051 [ # # # # ]: 0 : new_gvec(:,1:in_ng) = in_Gsph%gvec
2052 : :
2053 : 0 : call new_Gsph%free()
2054 : 0 : call new_Gsph%init(Cryst, new_ng, gvec=new_gvec)
2055 [ # # ]: 0 : ABI_FREE(new_gvec)
2056 : :
2057 : : else
2058 : : ! new_gpsh smaller/equal than in_gsph
2059 : 2 : ierr = 0
2060 [ + + ]: 132 : do ig=1,MIN(new_Gsph%ng,in_Gsph%ng)
2061 [ + + + - : 522 : if (ANY(new_Gsph%gvec(:,ig) /= in_Gsph%gvec(:,ig)) ) then
- + ]
2062 : 0 : ierr = ierr + 1
2063 : 0 : write(std_out,*)" new_gvec, in_gvec",ig,new_Gsph%gvec(:,ig),in_Gsph%gvec(:,ig)
2064 : : end if
2065 : : end do
2066 [ - + ]: 2 : ABI_CHECK(ierr==0,"Fatal error")
2067 : : end if
2068 : :
2069 : : end subroutine gsph_extend
2070 : : !!***
2071 : :
2072 : : !----------------------------------------------------------------------
2073 : :
2074 : : !!****f* m_gsphere/getkpgnorm
2075 : : !! NAME
2076 : : !! getkpgnorm
2077 : : !!
2078 : : !! FUNCTION
2079 : : !! compute the norms of the k+G vectors
2080 : : !!
2081 : : !! INPUTS
2082 : : !! gprimd(3,3)=metric tensor
2083 : : !! kg_k(3,npw_k)= G vectors, in reduced coordinates
2084 : : !! kpt(3)=k vector, in reduced coordinates
2085 : : !! npw_k=size of the G-vector set
2086 : : !!
2087 : : !! OUTPUT
2088 : : !! kpgnorm(npw_k)=norms of the k+G vectors
2089 : : !!
2090 : : !! SOURCE
2091 : :
2092 : 306 : subroutine getkpgnorm(gprimd,kpt,kg_k,kpgnorm,npw_k)
2093 : :
2094 : : !Arguments ------------------------------------
2095 : : !scalars
2096 : : integer,intent(in) :: npw_k
2097 : : !arrays
2098 : : integer,intent(in) :: kg_k(3,npw_k)
2099 : : real(dp),intent(in) :: gprimd(3,3),kpt(3)
2100 : : real(dp),intent(out) :: kpgnorm(npw_k)
2101 : :
2102 : : !Local variables-------------------------------
2103 : : !scalars
2104 : : integer :: ipw
2105 : : real(dp) :: g11,g12,g13,g21,g22,g23,g31,g32,g33,k1,k2,k3,kpg1,kpg2,kpg3,rr,xx
2106 : : real(dp) :: yy,zz
2107 : : ! *************************************************************************
2108 : :
2109 : 306 : k1=kpt(1) ; k2=kpt(2) ; k3=kpt(3)
2110 : 306 : g11=gprimd(1,1)
2111 : 306 : g12=gprimd(1,2)
2112 : 306 : g13=gprimd(1,3)
2113 : 306 : g21=gprimd(2,1)
2114 : 306 : g22=gprimd(2,2)
2115 : 306 : g23=gprimd(2,3)
2116 : 306 : g31=gprimd(3,1)
2117 : 306 : g32=gprimd(3,2)
2118 : 306 : g33=gprimd(3,3)
2119 : :
2120 : : !Loop over all k+G
2121 [ + + ]: 68931 : do ipw=1,npw_k
2122 : :
2123 : : ! Load k+G
2124 : 68625 : kpg1=k1+dble(kg_k(1,ipw))
2125 : 68625 : kpg2=k2+dble(kg_k(2,ipw))
2126 : 68625 : kpg3=k3+dble(kg_k(3,ipw))
2127 : :
2128 : : ! Calculate module of k+G
2129 : 68625 : xx=g11*kpg1+g12*kpg2+g13*kpg3
2130 : 68625 : yy=g21*kpg1+g22*kpg2+g23*kpg3
2131 : 68625 : zz=g31*kpg1+g32*kpg2+g33*kpg3
2132 : 68625 : rr=sqrt(xx**2+yy**2+zz**2)
2133 : 68931 : kpgnorm(ipw) = rr
2134 : :
2135 : : end do ! ipw
2136 : :
2137 : 306 : end subroutine getkpgnorm
2138 : : !!***
2139 : :
2140 : : !!****f* m_gsphere/symg
2141 : : !! NAME
2142 : : !! symg
2143 : : !!
2144 : : !! FUNCTION
2145 : : !! Treat symmetries applied to the G vectors, in view of the application
2146 : : !! to symmetrization of the dielectric matrix.
2147 : : !! Generate a list of time-reversed G vectors, as well as a list
2148 : : !! of spatially-symmetric G vectors.
2149 : : !!
2150 : : !! INPUTS
2151 : : !! kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
2152 : : !! npwdiel=number of planewaves for the dielectric matrix
2153 : : !! nsym=number of symmetry
2154 : : !! symrel(3,3,nsym)=symmetry matrices in real space (integers)
2155 : : !! tnons(3,nsym)=reduced nonsymmorphic translations
2156 : : !! (symrel and tnons are in terms of real space primitive translations)
2157 : : !!
2158 : : !! OUTPUT
2159 : : !! phdiel(2,npwdiel,nsym)=phase associated with point symmetries applied to G
2160 : : !! sym_g(npwdiel,nsym)=index list of symmetric G vectors
2161 : : !! (could save a bit of space by suppressing isym=1, since the
2162 : : !! corresponding symmetry is the identity)
2163 : : !! tmrev_g(npwdiel)=index list of inverted G vectors (time-reversed)
2164 : : !!
2165 : : !! SOURCE
2166 : :
2167 : 82 : subroutine symg(kg_diel,npwdiel,nsym,phdiel,sym_g,symrel,tmrev_g,tnons)
2168 : :
2169 : : !Arguments ------------------------------------
2170 : : !scalars
2171 : : integer,intent(in) :: npwdiel,nsym
2172 : : !arrays
2173 : : integer,intent(in) :: kg_diel(3,npwdiel),symrel(3,3,nsym)
2174 : : integer,intent(out) :: sym_g(npwdiel,nsym),tmrev_g(npwdiel)
2175 : : real(dp),intent(in) :: tnons(3,nsym)
2176 : : real(dp),intent(out) :: phdiel(2,npwdiel,nsym)
2177 : :
2178 : : !Local variables-------------------------------
2179 : : !scalars
2180 : : integer :: g1,g2,g3,ipw,isym,j1,j2,j3,m1m,m1p,m2m,m2p,m3m,m3p,symmg,trevg
2181 : : real(dp) :: arg,tau1,tau2,tau3
2182 : : !character(len=500) :: msg
2183 : : !arrays
2184 : 82 : integer,allocatable :: grid(:,:,:)
2185 : : ! *************************************************************************
2186 : :
2187 : : !Determines maximal bounds of the zone spanned by the planewaves
2188 : 82 : m1m=0 ; m2m=0 ; m3m=0 ; m1p=0 ; m2p=0 ; m3p=0
2189 [ + + ]: 12500 : do ipw=1,npwdiel
2190 : 12418 : g1=kg_diel(1,ipw)
2191 : 12418 : g2=kg_diel(2,ipw)
2192 : 12418 : g3=kg_diel(3,ipw)
2193 : 12418 : if(g1<m1m)m1m=g1 ; if(g1>m1p)m1p=g1
2194 : 12418 : if(g2<m2m)m2m=g2 ; if(g2>m2p)m2p=g2
2195 : 12500 : if(g3<m3m)m3m=g3 ; if(g3>m3p)m3p=g3
2196 : : end do
2197 : :
2198 : : !Set up grid, that associate to each point the index of the
2199 : : !corresponding planewave, if there is one
2200 [ - + + - : 410 : ABI_MALLOC(grid, (m1m:m1p,m2m:m2p,m3m:m3p))
- + + - -
+ + - + -
+ - - + -
+ ]
2201 [ + + + + : 28718 : grid(:,:,:)=0
+ + ]
2202 [ + + ]: 12500 : do ipw=1,npwdiel
2203 : 12418 : g1=kg_diel(1,ipw)
2204 : 12418 : g2=kg_diel(2,ipw)
2205 : 12418 : g3=kg_diel(3,ipw)
2206 : 12500 : grid(g1,g2,g3)=ipw
2207 : : end do
2208 : :
2209 : : !Set up tmrev_g and sym_g arrays
2210 [ + + ]: 12500 : do ipw=1,npwdiel
2211 : 12418 : g1=kg_diel(1,ipw)
2212 : 12418 : g2=kg_diel(2,ipw)
2213 : 12418 : g3=kg_diel(3,ipw)
2214 : :
2215 : : ! Treat first time-reversal symmetry
2216 : 12418 : trevg=grid(-g1,-g2,-g3)
2217 [ - + ]: 12418 : if(trevg==0)then
2218 : 0 : ABI_BUG('Do not find the time-reversed symmetric of a G-vector.')
2219 : : end if
2220 : 12418 : tmrev_g(ipw)=trevg
2221 : :
2222 : : ! Treat now spatial symmetries
2223 [ + + ]: 98240 : do isym=1,nsym
2224 : :
2225 : : ! Get rotated G vector Gj for each symmetry element
2226 : : ! -- here we use the TRANSPOSE of symrel; assuming symrel expresses
2227 : : ! the rotation in real space, the transpose is then appropriate
2228 : : ! for G space symmetrization (according to Doug : see routine irrzg.f)
2229 : : j1=symrel(1,1,isym)*g1+&
2230 : 85740 : & symrel(2,1,isym)*g2+symrel(3,1,isym)*g3
2231 : : j2=symrel(1,2,isym)*g1+&
2232 : 85740 : & symrel(2,2,isym)*g2+symrel(3,2,isym)*g3
2233 : : j3=symrel(1,3,isym)*g1+&
2234 : 85740 : & symrel(2,3,isym)*g2+symrel(3,3,isym)*g3
2235 : 85740 : symmg=grid(j1,j2,j3)
2236 [ - + ]: 85740 : if(symmg==0)then
2237 : 0 : ABI_BUG('Do not find the spatially symmetric of a G-vector.')
2238 : : end if
2239 : 85740 : sym_g(ipw,isym)=symmg
2240 : :
2241 : : ! Get associated phase
2242 : 85740 : tau1=tnons(1,isym)
2243 : 85740 : tau2=tnons(2,isym)
2244 : 85740 : tau3=tnons(3,isym)
2245 [ + + + - : 98158 : if (abs(tau1)>tol12.or.abs(tau2)>tol12.or.abs(tau3)>tol12) then
+ + ]
2246 : : ! compute exp(-2*Pi*I*G dot tau) using original G
2247 : 31800 : arg=two_pi*(dble(g1)*tau1+dble(g2)*tau2+dble(g3)*tau3)
2248 : 31800 : phdiel(1,ipw,isym)=cos(arg)
2249 : 31800 : phdiel(2,ipw,isym)=-sin(arg)
2250 : : else
2251 : 53940 : phdiel(1,ipw,isym)=1._dp
2252 : 53940 : phdiel(2,ipw,isym)=0._dp
2253 : : end if
2254 : :
2255 : : end do
2256 : : end do
2257 : :
2258 [ - + ]: 82 : ABI_FREE(grid)
2259 : :
2260 : 82 : end subroutine symg
2261 : : !!***
2262 : :
2263 [ - + - - : 719870 : end module m_gsphere
+ + - + +
- - + - +
+ - - + -
+ - + - +
- + - + -
+ ]
2264 : : !!***
|
