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 1066 : 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 716498 : do ig1=1,npw
225 715432 : ish1=g2sh(ig1) ; ss=shlim(ish1) ; ee=shlim(ish1+1)-1
226 2861728 : gbase(:)=gvec(:,ig1)
227 :
228 2141747 : do itim=1,timrev
229 35359841 : do isym=1,nsym
230 531506560 : grot=(3-2*itim)*MATMUL(symrec(:,:,isym),gbase)
231 : found=.FALSE.
232 : ! * Loop on the shell of ig1 to speed up the search.
233 2033871256 : do ig2=ss,ee
234 2296357632 : if (ALL(ABS(grot(:)-gvec(:,ig2))==0)) then
235 33219160 : found=.TRUE.
236 33219160 : grottb (ig1,itim,isym)=ig2
237 33219160 : grottbm1(ig2,itim,isym)=ig1
238 : end if
239 : end do
240 34644409 : 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 1066 : 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 2026 : 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 1066 : integer,allocatable :: shlim(:)
299 1066 : integer,pointer :: symrec(:,:,:),gvec_ptr(:,:)
300 : real(dp) :: kptns1(3,nkpt1)
301 1066 : real(dp),allocatable :: shlen(:)
302 1066 : real(dp),pointer :: tnons(:,:)
303 : !************************************************************************
304 :
305 : DBG_ENTER("COLL")
306 :
307 : ! Copy info on symmetries
308 1066 : nsym = Cryst%nsym
309 1066 : timrev = Cryst%timrev
310 1066 : symrec => Cryst%symrec
311 1066 : tnons => Cryst%tnons
312 :
313 : ! Initialize the object
314 1066 : Gsph%istwfk = 1 ! Time reversal is not used here.
315 1066 : Gsph%nsym = nsym
316 1066 : Gsph%timrev = timrev
317 :
318 13858 : Gsph%gmet = Cryst%gmet
319 13858 : Gsph%gprimd = Cryst%gprimd
320 :
321 1066 : if (PRESENT(gvec)) then
322 960 : if (PRESENT(ecut)) then
323 0 : ABI_BUG("ecut cannot be present when gvec is used")
324 : end if
325 960 : Gsph%ng= ng
326 2880 : ABI_MALLOC(Gsph%gvec,(3,ng))
327 1931848 : Gsph%gvec=gvec
328 : !
329 : ! Calculate cutoff energy of the sphere.
330 : max_ecut=-one
331 483442 : do ig=1,ng
332 482482 : g1=gvec(1,ig)
333 482482 : g2=gvec(2,ig)
334 482482 : 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 482482 : two*(Cryst%gmet(1,2)*g1*g2+Cryst%gmet(1,3)*g1*g3+Cryst%gmet(2,3)*g2*g3)
337 483442 : max_ecut=MAX(max_ecut,gsq)
338 : end do
339 960 : max_ecut=two*max_ecut*pi**2
340 960 : 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 4264 : ABI_MALLOC(Gsph%phmGt, (Gsph%ng, nsym))
357 43674 : do isym=1,nsym
358 16763342 : do ig=1,Gsph%ng
359 66921280 : 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 3198 : ABI_MALLOC(Gsph%phmSGt,(Gsph%ng,nsym))
366 716498 : do ig=1,Gsph%ng
367 17436166 : do isym=1,nsym
368 518309708 : sg=MATMUL(symrec(:,:,isym),Gsph%gvec(:,ig))
369 67594104 : 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 4264 : ABI_CHECK(ALL(Gsph%gvec(1:3,1)==0), 'First G must be 0')
378 :
379 3198 : ABI_MALLOC(Gsph%g2sh,(Gsph%ng))
380 1066 : 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 3198 : ABI_MALLOC(shlim,(Gsph%ng+1))
384 1066 : shlim(1)=1
385 :
386 : ! For each shell, gives the radius of the shell.
387 3198 : ABI_MALLOC(shlen,(Gsph%ng))
388 1066 : shlen(1)=zero
389 :
390 1066 : nsh=1; norm_old=zero
391 715432 : do ig=2,Gsph%ng
392 13572954 : norm=two_pi*SQRT(DOT_PRODUCT(Gsph%gvec(:,ig),MATMUL(Cryst%gmet,Gsph%gvec(:,ig))))
393 714366 : eps=norm*tol8
394 714366 : if (ABS(norm-norm_old)>eps) then
395 19085 : norm_old = norm; nsh = nsh + 1
396 19085 : shlim(nsh)=ig
397 19085 : shlen(nsh)=norm
398 : end if
399 715432 : Gsph%g2sh(ig)=nsh
400 : end do
401 1066 : shlim(nsh+1)=Gsph%ng+1
402 :
403 : ! Save info on the shells
404 1066 : Gsph%nsh = nsh
405 3198 : ABI_MALLOC(Gsph%shlim, (nsh+1))
406 23349 : Gsph%shlim = shlim(1:nsh+1)
407 3198 : ABI_MALLOC(Gsph%shlen, (nsh ))
408 22283 : Gsph%shlen = shlen(1:nsh)
409 1066 : ABI_FREE(shlim)
410 1066 : ABI_FREE(shlen)
411 :
412 : ! Calculate tables for rotated G"s
413 5330 : ABI_MALLOC(Gsph%rottb , (Gsph%ng,timrev,nsym))
414 4264 : ABI_MALLOC(Gsph%rottbm1, (Gsph%ng,timrev,nsym))
415 :
416 : call setup_G_rotation(nsym, symrec, timrev, Gsph%ng, Gsph%gvec,&
417 1066 : 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 3198 : ABI_MALLOC(gsph%g2mg, (gsph%ng))
422 :
423 716498 : do ig=1,gsph%ng
424 715432 : ish=gsph%g2sh(ig)
425 715432 : ss=gsph%shlim(ish); ee=gsph%shlim(ish+1)-1
426 2861728 : gsearch = -gsph%gvec(:,ig)
427 25258959 : img = 0
428 : ! Loop on shells to speed up the search.
429 25258959 : do isearch=ss,ee
430 29266560 : if (all(gsph%gvec(:,isearch) == gsearch)) then
431 : img = isearch; exit
432 : end if
433 : end do
434 715432 : if (img==0) ABI_ERROR("Cannot find -G in G-sphere!")
435 716498 : gsph%g2mg(ig) = img
436 : end do
437 :
438 : !call Gsph%print(unit=std_out,prtvol=1)
439 : DBG_EXIT("COLL")
440 :
441 2026 : 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 104432 : 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 104432 : type(MPI_type) :: MPI_enreg_seq
484 : !arrays
485 104432 : integer,allocatable :: gmg0(:,:)
486 104432 : logical,allocatable :: kg_mask(:)
487 : ! *************************************************************************
488 :
489 417728 : if (mgfft/=MAXVAL(ngfft(1:3))) then
490 0 : ABI_ERROR("mgfft/-MAXVAL(ngfft(1:3)")
491 : end if
492 :
493 104432 : ng = Gsph%ng
494 :
495 104432 : ierr=0; use_padfft=0
496 313296 : ABI_MALLOC(gmg0,(3,ng))
497 6084642 : do ig=1,ng
498 23920840 : gmg0(:,ig) = Gsph%gvec(:,ig)-g0
499 : ! Consider possible wrap around errors.
500 47841680 : 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 8066604 : if (ALL(gmg0(:,ig) == 0)) use_padfft=1
506 : end do
507 :
508 104432 : 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 104432 : if (use_padfft == 1) call sphereboundary(gmg0_gbound,1,gmg0,mgfft,ng)
517 :
518 104432 : call initmpi_seq(MPI_enreg_seq) ! No FFT parallelism.
519 104432 : call MPI_enreg_seq%distribfft%init_seq('c',ngfft(2),ngfft(3),'all')
520 :
521 313296 : ABI_MALLOC(kg_mask, (ng))
522 104432 : call kgindex(gmg0_ifft, gmg0, kg_mask, MPI_enreg_seq, ngfft, ng)
523 :
524 6084642 : ABI_CHECK(ALL(kg_mask),"FFT para not yet implemented")
525 104432 : ABI_FREE(kg_mask)
526 :
527 104432 : ABI_FREE(gmg0)
528 104432 : call destroy_mpi_enreg(MPI_enreg_seq)
529 :
530 104432 : 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 1949 : 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 1949 : ABI_SFREE(Gsph%g2sh)
713 1949 : ABI_SFREE(Gsph%gvec)
714 1949 : ABI_SFREE(Gsph%g2mg)
715 1949 : ABI_SFREE(Gsph%rottb)
716 1949 : ABI_SFREE(Gsph%rottbm1)
717 1949 : ABI_SFREE(Gsph%shlim)
718 :
719 1949 : ABI_SFREE(Gsph%shlen)
720 :
721 : ! complex arrays
722 1949 : ABI_SFREE(Gsph%phmGt)
723 1949 : ABI_SFREE(Gsph%phmSGt)
724 :
725 : DBG_EXIT("COLL")
726 :
727 1949 : 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 469 : 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 469 : type(MPI_type) :: MPI_enreg_seq
945 : !arrays
946 938 : integer :: gcur(3),geq(3),dum_kg(3,0),symrec2t(3,3,nsym2)
947 469 : integer,allocatable :: gbase(:,:),gbasek(:,:,:)
948 469 : integer,allocatable :: gcurr(:,:),gshell(:,:),insort(:),gtmp(:,:)
949 469 : integer,allocatable :: nbasek(:),nshell(:),shlim(:), npwarr(:)
950 : real(dp) :: kpoint(3),gmet(3,3)
951 469 : 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 469 : call initmpi_seq(MPI_enreg_seq)
958 :
959 : !Compute reciprocal space metrics
960 1876 : do ii=1,3
961 : gmet(ii,:)=gprimd(1,ii)*gprimd(1,:)+&
962 : & gprimd(2,ii)*gprimd(2,:)+&
963 6097 : & 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 19207 : do isym=1,nsym2
969 244063 : symrec2t(:,:,isym)=TRANSPOSE(symrel2(:,:,isym))
970 : end do
971 : !
972 : ! ==============================================
973 : ! ==== Find irreducible G-vectors at each k ====
974 : ! ==============================================
975 :
976 1407 : ABI_MALLOC(npwarr,(nkpt))
977 469 : exchn2n3d=0; ikg=0
978 2722 : do ikpt=1,nkpt
979 9012 : kpoint=kptns(:,ikpt); istwf_k=1
980 2722 : call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,dum_kg,kpoint,0,MPI_enreg_seq,0,npwarr(ikpt))
981 : end do
982 2722 : mpw = MAXVAL(npwarr)
983 :
984 938 : ABI_MALLOC(nbasek,(nkpt))
985 1876 : ABI_MALLOC(gbasek,(3,mpw,nkpt))
986 1876 : ABI_MALLOC(cnormk,(mpw,nkpt))
987 2722 : nbasek=0 ! # of irreducible G at each k.
988 904022 : cnormk=zero ! Norm of each irreducible G.
989 3607922 : gbasek=0 ! The set of irreducible G"s at each k.
990 :
991 2722 : do ikpt=1,nkpt
992 :
993 9012 : kpoint = kptns(:,ikpt)
994 2253 : npw_k = npwarr(ikpt)
995 :
996 2253 : exchn2n3d=0; ikg=0; istwf_k=1
997 6759 : ABI_MALLOC(gcurr,(3,npw_k))
998 2253 : call kpgsph(ecut,exchn2n3d,gmet,ikg,0,istwf_k,gcurr,kpoint,mkmem_,MPI_enreg_seq,npw_k,onpw_k)
999 :
1000 9012 : 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 2253 : call get_irredg(npw_k,nsym2,pinv,gprimd,symrec2t,gcurr,nbasek(ikpt),gbasek(:,:,ikpt),cnormk(:,ikpt))
1007 :
1008 2722 : 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 469 : sizepw=2*mpw
1014 1407 : ABI_MALLOC(gbase,(3,sizepw))
1015 1407 : ABI_MALLOC(cnorm,(sizepw))
1016 : nbase=0 ! # of irred G found.
1017 :
1018 469 : call merge_kgirr(nsym2,pinv,nkpt,mpw,sizepw,symrec2t,nbasek,cnormk,gbasek,nbase,gbase,cnorm,ierr)
1019 469 : if (ierr/=0) then
1020 0 : ABI_ERROR('merge_kgirr returned a non-zero status error')
1021 : end if
1022 :
1023 469 : ABI_FREE(nbasek)
1024 469 : ABI_FREE(cnormk)
1025 469 : 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 1407 : ABI_MALLOC(gshell,(3,2*nsym2))
1031 1407 : ABI_MALLOC(shlim,(nbase))
1032 :
1033 938 : 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 938 : ABI_MALLOC(insort,(nbase))
1043 938 : ABI_MALLOC(nshell,(nbase))
1044 37400 : do in=1,nbase
1045 37400 : insort(in)=in
1046 : end do
1047 469 : call sort_dp(nbase,cnorm,insort,tol14)
1048 : !
1049 : !Loop over all different modules of g''s (=shells):
1050 469 : maxpw=0
1051 37400 : do in=1,nbase
1052 36931 : nshell(in)=0
1053 147724 : gcur(:)=gbase(:,insort(in))
1054 :
1055 853453 : do is=1,nsym2 ! Loop over all symmetries:
1056 1669975 : do iinv=pinv,1,2
1057 3266088 : geq(:)=iinv*(symrel2(1,:,is)*gcur(1)+symrel2(2,:,is)*gcur(2)+symrel2(3,:,is)*gcur(3))
1058 :
1059 : found=.FALSE.; ish=1
1060 7050483 : do while ((.not.found) .and. (ish<=nshell(in))) ! Search for symmetric of g and eventually add it:
1061 8227704 : found=ALL(geq(:)==gshell(:,ish))
1062 6730442 : ish=ish+1
1063 : end do
1064 1633044 : if (.not.found) then
1065 496481 : nshell(in)=nshell(in)+1
1066 1985924 : gshell(:,nshell(in))=geq(:)
1067 : end if
1068 : end do
1069 : end do
1070 :
1071 36931 : 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 533412 : do ig=1,nshell(in)
1105 2022855 : gbig(:,ig+maxpw)=gshell(:,ig)
1106 : end do
1107 37400 : maxpw=maxpw+nshell(in)
1108 :
1109 : end do ! End loop over shells
1110 : !
1111 : ! * Compute shell limits
1112 469 : ilim=0
1113 37400 : do in=1,nbase
1114 36931 : ilim=ilim+nshell(in)
1115 37400 : shlim(in)=ilim
1116 : end do
1117 :
1118 469 : if (PRESENT(shlim_p)) then ! Return shlim_p
1119 726 : ABI_MALLOC(shlim_p,(nbase))
1120 20880 : shlim_p = shlim
1121 : end if
1122 :
1123 : ! Re-allocate gbig with correct sizes so that caller can inquire the size
1124 1407 : ABI_MALLOC(gtmp,(3,ilim))
1125 1986862 : gtmp = gbig(:,1:ilim)
1126 469 : ABI_FREE(gbig)
1127 938 : ABI_MALLOC(gbig,(3,ilim))
1128 1986393 : gbig=gtmp
1129 469 : ABI_FREE(gtmp)
1130 :
1131 469 : 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 469 : ABI_FREE(gshell)
1144 469 : ABI_FREE(insort)
1145 469 : ABI_FREE(nshell)
1146 : #endif
1147 :
1148 469 : call destroy_mpi_enreg(MPI_enreg_seq)
1149 469 : ABI_FREE(gbase)
1150 469 : ABI_FREE(shlim)
1151 469 : ABI_FREE(cnorm)
1152 469 : ABI_FREE(npwarr)
1153 :
1154 938 : 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 3526 : 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 3526 : 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 4569461 : nbasek = 0; cnormk(:) = zero; gbasek(:,:) = 0
1356 :
1357 916713 : do ig=1,npw_k
1358 3652748 : gcur(:) = gcurr(:,ig); norm = zero
1359 3652748 : do jj=1,3
1360 2739561 : gcar(jj)=gcur(1)*gprimd(jj,1)+gcur(2)*gprimd(jj,2)+gcur(3)*gprimd(jj,3)
1361 3652748 : norm=norm+gcar(jj)**2
1362 : end do
1363 913187 : eps = tol8 * norm; found = .False.; irr = 1
1364 113280552 : do while (.not.found .and. irr <= nbasek) ! This loop can be optimized by looping inside the shell.
1365 112367365 : if (abs(norm - cnormk(irr)) <= eps) then
1366 7264036 : gbas(:) = gbasek(:,irr); isym = 1
1367 31324101 : do while (.not.found .and. isym <= nsym)
1368 383605196 : geq(:) = matmul(symrec(:,:,isym),gcur)
1369 34650661 : found = all(geq(:) == gbas(:))
1370 29508092 : if (pinv == -1) found = (found .or. all(geq == -gbas)) ! For time-reversal
1371 31324101 : isym = isym + 1
1372 : end do
1373 : end if
1374 112445004 : irr = irr + 1
1375 : end do
1376 916713 : if (.not. found) then
1377 310556 : nbasek = nbasek + 1; cnormk(nbasek) = norm; gbasek(:,nbasek) = gcur(:)
1378 : end if
1379 : end do
1380 :
1381 : DBG_EXIT("COLL")
1382 :
1383 3526 : 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 469 : 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 469 : 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 469 : nbase=0 ; ierr=0
1447 2722 : do ikpt=1,nkpt
1448 78198 : do irgk=1,nbasek(ikpt)
1449 301904 : gcur(:)=gbasek(:,irgk,ikpt)
1450 75476 : norm=cnormk(irgk,ikpt) ; eps=tol8*norm
1451 75476 : found=.FALSE. ; inb=1
1452 9601384 : do while ((.not.found).and.(inb<=nbase))
1453 9525908 : if (ABS(norm-cnorm(inb))<=eps) then
1454 510904 : gbas(:)=gbase(:,inb)
1455 : isym=1
1456 1615525 : do while ((.not.found).and.(isym<=nsym))
1457 19341387 : geq(:)=MATMUL(symrec(:,:,isym),gcur)
1458 1720250 : found=ALL(geq(:)==gbas(:))
1459 1487799 : if (pinv==-1) found= (found.or.ALL(geq(:)==-gbas(:)) ) ! For time-reversal
1460 1615525 : isym=isym+1
1461 : end do
1462 : end if
1463 9562839 : inb=inb+1
1464 : end do
1465 77729 : if (.not.found) then
1466 : ! === Add to the list ===
1467 36931 : nbase=nbase+1
1468 36931 : 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 36931 : cnorm(nbase)=cnormk(irgk,ikpt)
1475 147724 : 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 1396 : 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 1396 : type(MPI_type) :: MPI_enreg_seq
1534 : !arrays
1535 : integer :: geq(3)
1536 1396 : integer,allocatable :: gvec(:,:),gvec_sh(:,:),insort(:),npw_sh(:)
1537 : real(dp) :: gctr(3)
1538 1396 : 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 1396 : if (npw==0.and.nsh==0.and.ecut<=tol6) then
1545 : write(msg,'(8a)')&
1546 176 : 'One of the three variables ecut',TRIM(tag),', npw',TRIM(tag),', or nsh',TRIM(tag),ch10,&
1547 352 : 'must be non-null. Returning.'
1548 176 : ABI_COMMENT(msg)
1549 176 : RETURN
1550 : end if
1551 : ! 2-> one and only one should be non-null
1552 1220 : 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 1220 : 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 1220 : 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 1220 : gctr(:)= [zero,zero,zero]
1574 1220 : 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 1092 : npwave=NINT(ucvol*(two*ecut)**1.5_dp/(six*pi**2))
1579 1092 : npwwrk=NINT(DBLE(npwave)*scale)+pad
1580 1092 : 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 1220 : call initmpi_seq(MPI_enreg_seq)
1592 :
1593 3660 : ABI_MALLOC(gvec,(3,npwwrk))
1594 2440 : ifound=0
1595 2440 : do while (ifound==0)
1596 : !write(msg,'(a,f8.2)')' setshells : ecut_trial = ',ecut_trial
1597 : !call wrtout(std_out,msg,'COLL')
1598 1220 : exchn2n3d=0 ! For the time being, no exchange of n2 and n3
1599 :
1600 1220 : call kpgsph(ecut_trial,exchn2n3d,gmet,0,1,1,gvec,gctr,1,MPI_enreg_seq,npwwrk,npw_found)
1601 :
1602 3660 : ABI_MALLOC(gnorm,(npw_found))
1603 3660 : ABI_MALLOC(insort,(npw_found))
1604 :
1605 359454 : do ig=1,npw_found
1606 358234 : insort(ig)=ig
1607 358234 : gnorm(ig)=zero
1608 1434156 : do ii=1,3
1609 : gnorm(ig)=gnorm(ig)+(gvec(1,ig)*gprimd(ii,1)+&
1610 : gvec(2,ig)*gprimd(ii,2)+&
1611 1432936 : gvec(3,ig)*gprimd(ii,3))**2
1612 : end do
1613 : end do
1614 1220 : call sort_dp(npw_found,gnorm,insort,tol14)
1615 :
1616 2440 : ABI_MALLOC(npw_sh,(npw_found))
1617 2440 : ABI_MALLOC(gnorm_sh,(npw_found))
1618 3660 : ABI_MALLOC(gvec_sh,(3,npw_found))
1619 359454 : npw_sh(:)=0
1620 359454 : gnorm_sh(:)=zero
1621 1434156 : gvec_sh(:,:)=0
1622 : ! Count the number of shells:
1623 : ! (search for the G-vectors generating the others by symmetry)
1624 1220 : nsh_found=0
1625 :
1626 359454 : do ig=1,npw_found
1627 358234 : eps=1.d-8*gnorm(ig)
1628 358234 : found=.FALSE.
1629 358234 : ish=1
1630 43087727 : do while ((.not.found).and.(ish<=nsh_found))
1631 42729493 : if (ABS(gnorm(ig)-gnorm_sh(ish))<=eps) then
1632 : isym=1
1633 11871425 : 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 44360616 : 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 11090154 : (geq(3)==gvec_sh(3,ish)))
1641 11417526 : isym=isym+1
1642 : end do
1643 : end if
1644 42729493 : ish=ish+1
1645 : end do
1646 359454 : if (.not.found) then
1647 30862 : nsh_found=nsh_found+1
1648 30862 : gnorm_sh(nsh_found)=gnorm(ig)
1649 123448 : gvec_sh(:,nsh_found)=gvec(:,insort(ig))
1650 30862 : npw_sh(nsh_found)=1
1651 : else
1652 327372 : ish=ish-1
1653 327372 : npw_sh(ish)=npw_sh(ish)+1
1654 : end if
1655 : end do
1656 :
1657 1220 : ecut_found=two*pi**2*gnorm(npw_found)
1658 :
1659 1220 : 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 1220 : 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 1220 : if (ecut<tol6) then
1736 128 : ecut=ecut_found
1737 : end if
1738 1220 : npw=npw_found
1739 1220 : nsh=nsh_found
1740 : end if
1741 :
1742 : end do ! while(ifound==0)
1743 :
1744 1220 : call destroy_mpi_enreg(MPI_enreg_seq)
1745 :
1746 1220 : ABI_FREE(gnorm)
1747 1220 : ABI_FREE(gnorm_sh)
1748 1220 : ABI_FREE(gvec)
1749 1220 : ABI_FREE(gvec_sh)
1750 1220 : ABI_FREE(insort)
1751 1220 : ABI_FREE(npw_sh)
1752 :
1753 1396 : 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 4252 : 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 4252 : integer,allocatable :: iwork(:,:,:)
1792 : !************************************************************************
1793 :
1794 7443324 : g1_max = maxval(abs(kg1))
1795 7011176 : g2_max = maxval(abs(kg2))
1796 17008 : do ii=1,3
1797 17008 : gmax(ii) = max(g1_max(ii), g2_max(ii))
1798 : end do
1799 17008 : gmax = 2*gmax + 1
1800 4252 : 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 21260 : ABI_MALLOC(iwork, (n1, n2, n3))
1804 :
1805 : ! Insert kg1 into work with extra 0 s around outside:
1806 6627552 : iwork = 0
1807 1859768 : do ipw=1,npw1
1808 1855516 : i1 = kg1(1,ipw); if (i1<0) i1=i1+n1; i1=i1+1
1809 1855516 : i2 = kg1(2,ipw); if (i2<0) i2=i2+n2; i2=i2+1
1810 1855516 : i3 = kg1(3,ipw); if (i3<0) i3=i3+n3; i3=i3+1
1811 1859768 : iwork(i1,i2,i3) = ipw
1812 : end do
1813 :
1814 1751731 : g2g1 = 0; nmiss = 0
1815 1751731 : do ipw=1,npw2
1816 1747479 : i1 = kg2(1,ipw); if (i1<0) i1=i1+n1; i1=i1+1
1817 1747479 : i2 = kg2(2,ipw); if (i2<0) i2=i2+n2; i2=i2+1
1818 1747479 : 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 1747479 : g2g1(ipw) = iwork(i1,i2,i3)
1824 1751731 : if (g2g1(ipw) == 0) nmiss = nmiss + 1
1825 : end do
1826 :
1827 4252 : ABI_FREE(iwork)
1828 :
1829 4252 : 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 717789 : end module m_gsphere
2264 : !!***
|
