Line data Source code
1 : !!****m* ABINIT/m_fft_mesh
2 : !! NAME
3 : !! m_fft_mesh
4 : !!
5 : !! FUNCTION
6 : !! This module contains routines and helper functions to perform the setup of the FFT mesh
7 : !! It also provides a set of tools to test the grid, rotate the mesh according to the symmetry
8 : !! operations of the space group etc.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2008-2026 ABINIT group (MG, XG, GMR, VO, LR, RWG, YMN, RS, TR, DC)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : MODULE m_fft_mesh
25 :
26 : use defs_basis
27 : use m_errors
28 : use m_abicore
29 : use m_hide_blas
30 : use, intrinsic :: iso_c_binding
31 :
32 : use defs_fftdata, only : size_goed_fft
33 : use m_fstrings, only : sjoin, itoa, ltoa
34 : use m_numeric_tools, only : denominator, mincm, iseven, pfactorize
35 : use m_matrix, only : mati3inv
36 : use m_geometry, only : xred2xcart
37 : use m_crystal, only : crystal_t
38 :
39 : implicit none
40 :
41 : private
42 :
43 : public :: setmesh ! Perform the setup of the FFT mesh for the GW oscillator strengths.
44 : public :: check_rot_fft ! Test whether the mesh is compatible with the rotational part of the space group.
45 : public :: fft_check_rotrans ! Test whether the mesh is compatible with the symmetries of the space group.
46 : public :: rotate_fft_mesh ! Calculate the FFT index of the rotated mesh.
47 : public :: denpot_project ! Compute n(r) + n( $R^{-1}(r-\tau)$ in) / 2
48 : ! Mainly used with R = inversion to select the even/odd part under inversion
49 : public :: cigfft ! Calculate the FFT index of G-G0.
50 : public :: ig2gfft ! Returns the component of a G in the FFT Box from its sequential index.
51 : public :: g2ifft ! Returns the index of the G in the FFT box from its reduced coordinates.
52 : public :: get_gfft ! Calculate the G-vectors in the FFT box from ngfft.
53 : public :: calc_ceigr ! e^{iG.r} on the FFT mesh (complex valued).
54 : public :: calc_eigr ! e^{iG.r} on the FFT mesh (version for real array with RE,IM).
55 : public :: calc_ceikr ! e^{ik.r} on the FFT mesh (complex valued).
56 : public :: times_eigr ! Multiply an array on the real-space mesh by e^{iG0.r}
57 : public :: times_eikr ! Multiply an array on the real-space mesh by e^{ik.r}
58 : public :: ctimes_eikr ! Version for complex array
59 : public :: phase ! Compute ph(ig)=$\exp(\pi\ i \ n/ngfft)$ for n=0,...,ngfft/2,-ngfft/2+1,...,-1
60 : public :: mkgrid_fft ! Sets the grid of fft (or real space) points to be treated.
61 : public :: supercell_fft
62 :
63 : interface calc_ceigr
64 : module procedure calc_ceigr_spc
65 : module procedure calc_ceigr_dpc
66 : end interface calc_ceigr
67 :
68 : interface calc_ceikr
69 : module procedure calc_ceikr_spc
70 : module procedure calc_ceikr_dpc
71 : end interface calc_ceikr
72 :
73 : !interface times_eikr
74 : ! module procedure times_eikr_dp
75 : ! module procedure ctimes_eikr_dpc
76 : !end interface times_eikr
77 : !!***
78 :
79 : !----------------------------------------------------------------------
80 :
81 : !!****t* m_fft_mesh/zpad_t
82 : !! NAME
83 : !! zpad_t
84 : !!
85 : !! FUNCTION
86 : !! Tables used for zero-padded FFTs.
87 : !!
88 : !! SOURCE
89 :
90 : type,public :: zpad_t
91 :
92 : integer :: nlinex
93 : ! Total number of 1D transforms.
94 :
95 : integer :: n_zplanes
96 : ! Number of z-planes intersecting the sphere.
97 :
98 : integer,allocatable :: zplane(:,:)
99 : ! zplane(3,n_zplanes)
100 : ! zplane(1,zpl) : mapping z-plane index -> FFT index_z
101 : ! zplane(2,zpl) : mapping z-plane index -> igb index in array gbound
102 :
103 : integer,allocatable :: linex2ifft_yz(:,:)
104 : ! linex2ifft_yz(2,nlinex)
105 : ! mapping 1D-FFT -> (FFT_index_y, FFT index_z)
106 : contains
107 : procedure :: init => zpad_init
108 : procedure :: free => zpad_free
109 : end type zpad_t
110 : !!***
111 :
112 : CONTAINS !========================================================================================
113 : !!***
114 :
115 : !----------------------------------------------------------------------
116 :
117 : !!****f* m_fft_mesh/zpad_init
118 : !! NAME
119 : !! zpad_init
120 : !!
121 : !! FUNCTION
122 : !! Creation method for zpad_t instance
123 : !!
124 : !! INPUTS
125 : !! mgfft=MAX(nx,ny,nz), only used to dimension gbound
126 : !! gbound(2*mgfft+8,2)= The boundaries of the basis sphere of G vectors at a given k-point.
127 : !! See sphereboundary for more info.
128 : !!
129 : !! SOURCE
130 :
131 5736024 : subroutine zpad_init(zpad, nx, ny, nz, ldx, ldy, ldz, mgfft, gbound)
132 :
133 : !Arguments ------------------------------------
134 : !scalars
135 : class(zpad_t),intent(out) :: zpad
136 : integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,mgfft
137 : !arrays
138 : integer,intent(in) :: gbound(2*mgfft+8,2)
139 :
140 : !Local variables-------------------------------
141 : integer :: jj,g3_max,g3_min,gg3,ifft_g3,igb,g2min,g2max,nlinex
142 : ! *************************************************************************
143 :
144 5736024 : g3_min = gbound(3, 2)
145 5736024 : g3_max = gbound(4, 2)
146 :
147 5736024 : zpad%n_zplanes = g3_max - g3_min + 1
148 :
149 17208072 : ABI_MALLOC(zpad%zplane, (2, nz))
150 17208072 : ABI_MALLOC(zpad%linex2ifft_yz, (2, nx*ny*nz))
151 :
152 : ! Loop over the z-planes intersecting the G-sphere.
153 5736024 : nlinex = 0
154 35502948 : do gg3=1,zpad%n_zplanes
155 : !
156 29766924 : if (gg3<=g3_max+1) then
157 : ifft_g3 = gg3
158 : else
159 12319478 : ifft_g3 = gg3 + nz - zpad%n_zplanes ! Wrap around for negative gg3.
160 : end if
161 : !
162 : ! Select the set of y for this z-plane.
163 29766924 : igb=2*gg3+3
164 29766924 : g2min = gbound(igb ,2)
165 29766924 : g2max = gbound(igb+1,2)
166 :
167 29766924 : zpad%zplane(1,gg3) = ifft_g3
168 29766924 : zpad%zplane(2,gg3) = igb
169 :
170 : !(1:g2max+1,ifft_g3) ! Positive g_y.
171 : !(g2min+ny+1:ny,ifft_g3) ! Negative g_y.
172 :
173 99287341 : do jj=1,g2max+1
174 69520417 : nlinex = nlinex + 1
175 69520417 : zpad%linex2ifft_yz(1,nlinex) = jj
176 99287341 : zpad%linex2ifft_yz(2,nlinex) = ifft_g3
177 : end do
178 :
179 78049234 : do jj=g2min+ny+1,ny
180 42546286 : nlinex = nlinex + 1
181 42546286 : zpad%linex2ifft_yz(1,nlinex) = jj
182 72313210 : zpad%linex2ifft_yz(2,nlinex) = ifft_g3
183 : end do
184 : end do
185 :
186 5736024 : zpad%nlinex = nlinex
187 :
188 : RETURN
189 : ABI_UNUSED((/ldx,ldy,ldz/))
190 :
191 : end subroutine zpad_init
192 : !!***
193 :
194 : !----------------------------------------------------------------------
195 :
196 : !!****f* m_fft_mesh/zpad_free
197 : !! NAME
198 : !! zpad_free
199 : !!
200 : !! FUNCTION
201 : !! Free dynamic memory
202 : !!
203 : !! SOURCE
204 :
205 5736024 : subroutine zpad_free(zpad)
206 :
207 : !Arguments ------------------------------------
208 : class(zpad_t),intent(inout) :: zpad
209 : ! *************************************************************************
210 :
211 5736024 : ABI_SFREE(zpad%zplane)
212 5736024 : ABI_SFREE(zpad%linex2ifft_yz)
213 :
214 5736024 : end subroutine zpad_free
215 : !!***
216 :
217 : !----------------------------------------------------------------------
218 :
219 : !!****f* m_fft_mesh/setmesh
220 : !!
221 : !! NAME
222 : !! setmesh
223 : !!
224 : !! FUNCTION
225 : !! Calculate the size of the FFT grid for the GW calculation.
226 : !!
227 : !! INPUTS
228 : !! gmet(3,3)=Reciprocal space metric.
229 : !! gvec(3,npwvec)=G-vectors in reduced coordinates.
230 : !! npwvec=Number of G vectors in the array gvec max(npwwfn,npwsigx)
231 : !! npwsigx=Size of the dielectric or self-energy matrix.
232 : !! npwwfn=Number of G-vectors in the wavefunctions.
233 : !! method=Integer flag for FFT grid (see below)
234 : !! mG0=Number of shells that must be added to take into account umklapp processes.
235 : !! Cryst<crystal_t>=Data type gathering information on unit cell and symmetries
236 : !! %nsym=Number of symmetry operations in the SG.
237 : !! %symrel(3,3,nsym)=Symmetry operations in real space.
238 : !! %tnons(3,nsym)=Fractional translations.
239 : !! enforce_sym=Flag to enforce a FFT which fulfils all symmetry operations, both the
240 : !! rotational part and fractional translations.
241 : !! [unit]=Output unit, defaults to std_out
242 : !!
243 : !! OUTPUT
244 : !! ngfft(18)=contain all needed information about 3D FFT,
245 : !! see also ~abinit/doc/variables/vargs.htm#ngfft
246 : !! nfftot= ngfft(1)*ngfft(2)*ngfft(3)=Total number of points in the FFT grid.
247 : !!
248 : !! NOTES
249 : !! Four methods are implemented for the calculation of the mesh:
250 : !! method=0 --> FFT mesh defined by the user, useful for debugging.
251 : !! method=1 Roughly takes the FFT box which encloses the larger of the two spheres of radius
252 : !! aliasing_factor * rwfn and rsigx, where rwfn and rsigx are the radius of the spheres
253 : !! with npwwfn and npwsigx planewaves respectively. The default aliasing_factor is 1.
254 : !! method=2 --> Calculates the optimal FFT grid which allows aliasing only outside the sphere of the
255 : !! npwsigx planewaves (finer than method=1 with aliasing_factor=1).
256 : !! method=3 --> Calculates the FFT grid needed to expand the density.
257 : !! (even finer than method=2, roughly corresponds to method=1 with aliasing_factor=2).
258 : !!
259 : !! See defs_fftdata for a list of allowed sizes of FFT.
260 : !!
261 : !! SOURCE
262 :
263 579 : subroutine setmesh(gmet, gvec, ngfft, npwvec, npwsigx, npwwfn, nfftot, method, mG0, Cryst, enforce_sym, unit)
264 :
265 : !Arguments ------------------------------------
266 : !scalars
267 : integer,intent(in) :: enforce_sym,method,npwsigx,npwvec,npwwfn
268 : integer,intent(out) :: nfftot
269 : integer,optional,intent(in) :: unit
270 : type(crystal_t),target,intent(in) :: Cryst
271 : !arrays
272 : integer,intent(in) :: gvec(3,npwvec),mG0(3)
273 : integer,intent(inout) :: ngfft(18)
274 : real(dp),intent(in) :: gmet(3,3)
275 :
276 : !Local variables ------------------------------
277 : !scalars
278 : integer :: aliasing_factor,fftalg,fftalga,fftalgc,ig,ig1,ig1max,ig2,ig2max,ig3,ig3max,ii,idx,ierr
279 : integer :: is,m1,m2,m3,mm1,mm2,mm3,n1,n2,n3,nsym,nt,ount
280 : real(dp) :: ecuteff,ecutsigx,ecutwfn,g1,g2,g3,gsq,gsqmax,reff,rsigx,rwfn
281 : logical :: fft_ok
282 : character(len=500) :: msg, tnons_warn
283 : !arrays
284 : integer :: fftnons(3),fftsym(3),mdum(3)
285 : !integer,allocatable :: pfactors(:),powers(:)
286 579 : integer,pointer :: symrel(:,:,:)
287 579 : real(dp),pointer :: tnons(:,:)
288 : !************************************************************************
289 :
290 : DBG_ENTER("COLL")
291 :
292 2316 : if (any(mg0 < 0)) then
293 0 : ABI_BUG(sjoin('Wrong mG0:', trim(ltoa(mG0))))
294 : end if
295 :
296 : tnons_warn = "Check your fractional translations tnons. "//ch10//&
297 : "Components should be a rational fraction in 1/8th or in 1/12th."//ch10//&
298 : "You may need to polish the structure by running AbiPy `abistruct.py abisanitize` on the input file."//ch10//&
299 579 : "to get rid of spurious tnons"
300 :
301 579 : ount = std_out; if (present(unit)) ount = unit
302 :
303 579 : nsym = Cryst%nsym
304 579 : symrel => Cryst%symrel
305 579 : tnons => Cryst%tnons
306 :
307 : ! Calculate the limits of the sphere of npwwfn G-vectors in each direction.
308 331814 : m1 = MAXVAL(ABS(gvec(1,1:npwwfn)))
309 331814 : m2 = MAXVAL(ABS(gvec(2,1:npwwfn)))
310 331814 : m3 = MAXVAL(ABS(gvec(3,1:npwwfn)))
311 :
312 : ! Calculate the limits of the sphere of npsigx G-vectors in each direction.
313 : ! Ensure that G+G0 will fit into the FFT grid, where G is any of the npwsigx/npweps vectors
314 : ! and G0 is (i,j,k) [-nG0shell<i,j,k<nG0shell]. This is required when npwsigx>npwwfn since
315 : ! we have to take into account umklapp G0 vectors to evaluate the oscillator matrix elements
316 : ! (see rho_tw_g) or to symmetrize these quantities (see also cigfft).
317 117182 : mm1 = MAXVAL(ABS(gvec(1,1:npwsigx)))
318 117182 : mm2 = MAXVAL(ABS(gvec(2,1:npwsigx)))
319 117182 : mm3 = MAXVAL(ABS(gvec(3,1:npwsigx)))
320 :
321 579 : mm1=mm1+mG0(1)
322 579 : mm2=mm2+mG0(2)
323 579 : mm3=mm3+mG0(3)
324 :
325 : ! To avoid possible wrap-around errors in cigfft, it is safe to start
326 : ! with odd divisions so that the FFT box is centered on Gamma
327 : ! This holds only if npwsigx > npwwfn.
328 : !if (iseven(mm1)) mm1=mm1+1
329 : !if (iseven(mm2)) mm2=mm2+1
330 : !if (iseven(mm3)) mm3=mm3+1
331 :
332 579 : write(msg,'(2(2a,i8,a,3i6),2a,3i3)')ch10,&
333 579 : ' setmesh: npwwfn = ',npwwfn, '; Max (m1,m2,m3) = ',m1,m2,m3,ch10,&
334 579 : ' npweps/npwsigx= ',npwsigx,'; Max (mm1,mm2,mm3)= ',mm1,mm2,mm3,ch10,&
335 1158 : ' mG0 added = ',mG0(:)
336 579 : call wrtout(ount, msg)
337 : !
338 : ! === Different FFT grids according to method ==
339 579 : select case (method)
340 :
341 : case (0)
342 : ! * FFT mesh defined by user, useful for testing.
343 0 : n1=ngfft(1)
344 0 : n2=ngfft(2)
345 0 : n3=ngfft(3)
346 0 : write(msg,'(3(a,i3))')' Mesh size enforced by user = ',n1,'x',n2,'x',n3
347 0 : ABI_COMMENT(msg)
348 :
349 0 : ngfft(1)=n1
350 0 : ngfft(2)=n2
351 0 : ngfft(3)=n3
352 0 : ngfft(4)=2*(ngfft(1)/2)+1
353 0 : ngfft(5)=2*(ngfft(2)/2)+1
354 0 : ngfft(6)= ngfft(3)
355 : !ngfft(4:6)=ngfft(1:3)
356 0 : nfftot=n1*n2*n3
357 0 : RETURN
358 :
359 : case (1)
360 38 : aliasing_factor=1
361 38 : write(msg,'(2a,i3)')ch10,' using method 1 with aliasing_factor = ',aliasing_factor
362 38 : call wrtout(ount, msg)
363 : m1=m1*aliasing_factor
364 : m2=m2*aliasing_factor
365 38 : m3=m3*aliasing_factor
366 :
367 : case (2,3)
368 :
369 541 : ecutwfn=-one ! Calculate the radius of the sphere of npwwfn G-vectors.
370 315490 : do ig=1,npwwfn
371 314949 : g1=REAL(gvec(1,ig))
372 314949 : g2=REAL(gvec(2,ig))
373 314949 : g3=REAL(gvec(3,ig))
374 : gsq= gmet(1,1)*g1**2+gmet(2,2)*g2**2+gmet(3,3)*g3**2+ &
375 314949 : two*(gmet(1,2)*g1*g2+gmet(1,3)*g1*g3+gmet(2,3)*g2*g3)
376 315490 : ecutwfn=MAX(ecutwfn,gsq)
377 : end do
378 541 : rwfn=SQRT(ecutwfn); ecutwfn=two*ecutwfn*pi**2
379 :
380 : ! * Calculate the radius of the sphere of (npwsigx|npweps) G-vectors.
381 541 : ecutsigx=-one
382 111318 : do ig=1,npwsigx
383 110777 : g1=REAL(gvec(1,ig))
384 110777 : g2=REAL(gvec(2,ig))
385 110777 : g3=REAL(gvec(3,ig))
386 : gsq= gmet(1,1)*g1**2+gmet(2,2)*g2**2+gmet(3,3)*g3**2+ &
387 110777 : two*(gmet(1,2)*g1*g2+gmet(1,3)*g1*g3+gmet(2,3)*g2*g3)
388 111318 : ecutsigx=MAX(ecutsigx,gsq)
389 : end do
390 541 : rsigx=SQRT(ecutsigx); ecutsigx=two*ecutsigx*pi**2
391 :
392 : write(msg,'(a,f7.3,3a,f7.3,a)')&
393 541 : ' calculated ecutwfn = ',ecutwfn, ' [Ha] ',ch10,&
394 1082 : ' calculated ecutsigx/ecuteps = ',ecutsigx,' [Ha]'
395 541 : call wrtout(ount, msg)
396 : !
397 : ! In the calculation of the GW self-energy or of the RPA dielectric matrix,
398 : ! we have products $ \rho_{12}(r)=u_1*(r) u_2(r) $ of wavefunctions whose Fourier
399 : ! coefficients lie in the sphere of radius rwfn. Such products will have non
400 : ! vanishing Fourier coefficients in the whole sphere of radius 2*rwfn since:
401 : ! $ rho_{12}(G) = \sum_T u_1*(T) u_2(T+G) $.
402 : ! However, we only need the Fourier coefficients of $rho_{12}$ that lie in the sphere
403 : ! of radius rsigx. We can thus allow aliasing outside that sphere, so that the FFT box
404 : ! will only enclose a sphere of radius reff given by:
405 :
406 541 : reff=rsigx+rwfn
407 541 : if (method==3) reff=two*rwfn ! Yields back the GS FFT grid if full wavefunctions are considered.
408 541 : ecuteff=two*(pi*reff)**2
409 541 : gsqmax=reff**2
410 :
411 541 : write(msg,'(a,i2,a,f7.3,a)')' using method = ',method,' with ecuteff = ',ecuteff,' [Ha]'
412 541 : call wrtout(ount, msg)
413 : !
414 : ! === Search the limits of the reff sphere in each direction ===
415 : !ig1max=2*m1+1
416 : !ig2max=2*m2+1
417 : !ig3max=2*m3+1
418 541 : if (method==2) then
419 537 : ig1max=mm1+m1+1
420 537 : ig2max=mm2+m2+1
421 537 : ig3max=mm3+m3+1
422 : else if (method==3) then
423 4 : ig1max=MAX(2*m1+1,2*mm1+1,mm1+m1+1)
424 4 : ig2max=MAX(2*m2+1,2*mm2+1,mm2+m2+1)
425 4 : ig3max=MAX(2*m3+1,2*mm3+1,mm3+m3+1)
426 : else
427 : ABI_BUG(sjoin("Wrong method:", itoa(method)))
428 : end if
429 :
430 541 : m1=-1; m2=-1; m3=-1
431 5607 : do ig1=0,ig1max
432 60809 : do ig2=0,ig2max
433 731128 : do ig3=0,ig3max
434 670860 : g1=REAL(ig1)
435 670860 : g2=REAL(ig2)
436 670860 : g3=REAL(ig3)
437 : gsq= gmet(1,1)*g1**2+gmet(2,2)*g2**2+gmet(3,3)*g3**2+ &
438 670860 : two*(gmet(1,2)*g1*g2+gmet(1,3)*g1*g3+gmet(2,3)*g2*g3)
439 670860 : if (gsq>gsqmax+tol6) CYCLE ! tol6 to improve portability
440 342171 : m1=MAX(m1,ig1)
441 342171 : m2=MAX(m2,ig2)
442 726062 : m3=MAX(m3,ig3)
443 : end do
444 : end do
445 : end do
446 :
447 : case default
448 579 : ABI_BUG(sjoin('Method > 3 or < 0 not allowed in setmesh while method:', itoa(method)))
449 : end select
450 : !
451 : ! * Warning if low npwwfn.
452 579 : if (m1<mm1 .or. m2<mm2 .or. m3<mm3) then
453 : write(msg,'(5a)')&
454 5 : 'Note that npwwfn is small with respect to npweps or with respect to npwsigx. ',ch10,&
455 5 : 'Such a small npwwfn is a waste: ',ch10,&
456 10 : 'You could raise npwwfn without loss in cpu time. '
457 5 : ABI_COMMENT(msg)
458 : end if
459 : !
460 : ! Keep the largest of the m/mm and and find the FFT grid which is compatible
461 : ! with the library and, if required, with the symmetry operations.
462 579 : m1=MAX(m1,mm1)
463 579 : m2=MAX(m2,mm2)
464 579 : m3=MAX(m3,mm3)
465 :
466 579 : if (enforce_sym==0) then
467 : ! === Determine the best size for the FFT grid *without* considering the symm ops ===
468 : ! * Ideally n=2*m+1 but this could not be allowed by the FFT library.
469 0 : call size_goed_fft(m1, n1, ierr)
470 0 : ABI_CHECK(ierr == 0, sjoin("size_goed_fft failed", ch10, tnons_warn))
471 0 : call size_goed_fft(m2, n2, ierr)
472 0 : ABI_CHECK(ierr == 0, sjoin("size_goed_fft failed", ch10, tnons_warn))
473 0 : call size_goed_fft(m3, n3, ierr)
474 0 : ABI_CHECK(ierr == 0, sjoin("size_goed_fft failed", ch10, tnons_warn))
475 :
476 0 : nfftot=n1*n2*n3
477 :
478 : ! * Check if the FFT is compatible, write ONLY a warning if it breaks the symmetry
479 0 : fftnons(1)=n1
480 0 : fftnons(2)=n2
481 0 : fftnons(3)=n3
482 0 : fft_ok=.TRUE.
483 0 : rd: do ii=1,3
484 0 : do is=1,nsym
485 0 : nt=denominator(tnons(ii,is), ierr)
486 0 : if (((fftnons(ii)/nt)*nt) /= fftnons(ii)) then
487 : fft_ok=.FALSE.; EXIT rd
488 : end if
489 : end do
490 : end do rd
491 : !
492 : ! Warn if not compatible with tnons or rotational part.
493 0 : if (.not.fft_ok) then
494 0 : ABI_WARNING('FFT mesh is not compatible with non-symmorphic translations')
495 : end if
496 0 : if (.not.(check_rot_fft(nsym,symrel,n1,n2,n3))) then
497 0 : ABI_WARNING('FFT mesh is not compatible with rotations')
498 : end if
499 :
500 : else
501 : ! === Determine the best size for the FFT grid considering symm ops ===
502 : ! * Ideally n=2*m+1 but this could not be allowed by the FFT library (at present only Goedecker)
503 579 : call wrtout(ount,' Finding a FFT mesh compatible with all the symmetries')
504 :
505 : ! 1) Find a FFT mesh compatible with the non-symmorphic operations
506 2316 : fftnons(:)=1
507 2316 : do ii=1,3
508 1737 : fftnons(ii)=1
509 73968 : do is=1,nsym
510 71652 : nt=denominator(tnons(ii,is), ierr)
511 73389 : if (((fftnons(ii)/nt)*nt)/=fftnons(ii)) fftnons(ii)=mincm(fftnons(ii),nt)
512 : end do
513 : end do
514 579 : write(msg,'(a,3(i0,1x))')' setmesh: divisor mesh ',fftnons(:)
515 579 : call wrtout(ount, msg)
516 : !
517 : ! 2) Check if also rotations preserve the grid.
518 : ! * Use previous m values as Initial guess.
519 579 : call size_goed_fft(m1,fftsym(1),ierr)
520 579 : ABI_CHECK(ierr == 0, sjoin("size_goed_fft failed", ch10, tnons_warn))
521 579 : call size_goed_fft(m2,fftsym(2),ierr)
522 579 : ABI_CHECK(ierr == 0, sjoin("size_goed_fft failed", ch10, tnons_warn))
523 579 : call size_goed_fft(m3,fftsym(3),ierr)
524 579 : ABI_CHECK(ierr == 0, sjoin("size_goed_fft failed", ch10, tnons_warn))
525 :
526 579 : mdum(1)=m1
527 579 : mdum(2)=m2
528 579 : mdum(3)=m3
529 :
530 579 : idx=0
531 1759 : do ! If a FFT division gets too large the code stops in size_goed_fft.
532 : if ( check_rot_fft(nsym,symrel,fftsym(1),fftsym(2),fftsym(3)) .and. &
533 : (MOD(fftsym(1),fftnons(1))==0) .and. &
534 4097 : (MOD(fftsym(2),fftnons(2))==0) .and. &
535 : (MOD(fftsym(3),fftnons(3))==0) &
536 : ) EXIT
537 1759 : ii=MOD(idx,3)+1
538 1759 : mdum(ii)=mdum(ii)+1
539 1759 : call size_goed_fft(mdum(ii),fftsym(ii),ierr)
540 1759 : ABI_CHECK(ierr == 0, sjoin("size_goed_fft failed", ch10, tnons_warn))
541 1759 : idx=idx+1
542 : end do
543 : !
544 : ! Got a good FFT grid, Calculate the number of FFT grid points
545 579 : n1=fftsym(1)
546 579 : n2=fftsym(2)
547 579 : n3=fftsym(3); nfftot=n1*n2*n3
548 :
549 : if (.not.( check_rot_fft(nsym,symrel,n1,n2,n3)) &
550 579 : .or.( MOD(fftsym(1),fftnons(1))/=0) .and. &
551 : ( MOD(fftsym(2),fftnons(2))/=0) .and. &
552 : ( MOD(fftsym(3),fftnons(3))/=0) &
553 : ) then
554 0 : ABI_BUG('Not able to generate a symmetric FFT')
555 : end if
556 : end if ! enforce_sym
557 :
558 : write(msg,'(3(a,i5),2a,i12,a)')&
559 579 : ' setmesh: FFT mesh size selected = ',n1,'x',n2,'x',n3,ch10,&
560 1158 : ' total number of points = ',nfftot,ch10
561 579 : call wrtout(ount, msg)
562 579 : if (ount /= dev_null) call wrtout(ab_out, msg)
563 :
564 579 : ngfft(1)=n1
565 579 : ngfft(2)=n2
566 579 : ngfft(3)=n3
567 579 : ngfft(4)=2*(ngfft(1)/2)+1
568 579 : ngfft(5)=2*(ngfft(2)/2)+1
569 579 : ngfft(6)= ngfft(3)
570 : !ngfft(4:6) = ngfft(1:3)
571 : !
572 : ! === Check the value of fftalg i.e ngfft(7) ===
573 : ! * Presently only Goedecker"s library or FFTW3 are allowed, see size_goed_fft.F90
574 579 : fftalg=ngfft(7); fftalga=fftalg/100; fftalgc=MOD(fftalg,10)
575 :
576 579 : if (all(fftalga /= [FFT_SG, FFT_FFTW3, FFT_DFTI]) ) then
577 0 : write(msg,'(6a)')ch10,&
578 0 : "Only Goedecker's routines with fftalg=1xx or FFTW3/DFTI routines are allowed in GW calculations. ",ch10,&
579 0 : "Action : check the value of fftalg in your input file, ",ch10,&
580 0 : "or modify setmesh.F90 to make sure the FFT mesh is compatible with the FFT library. "
581 0 : ABI_ERROR(msg)
582 : end if
583 :
584 : ! TODO Had to change setmesh to avoid bad values for FFTW3
585 : ! if (fftalga==3) then ! check whether mesh is optimal for FFTW3
586 : ! ABI_MALLOC(pfactors,(5))
587 : ! ABI_MALLOC(powers,(6))
588 : ! pfactors = (/2, 3, 5, 7, 11/)
589 : ! do ii=1,3
590 : ! call pfactorize(ngfft(ii),5,pfactors,powers)
591 : ! if (powers(6)/=1 .or. powers(4)/=0 .or. powers(5)/=0) then
592 : ! write(msg,'(a,i0,a)')&
593 : !& "ngfft(ii) ",ngfft(ii)," contains powers of 7-11 or greater; FFTW3 is not optimal "
594 : ! ABI_WARNING(msg)
595 : ! end if
596 : ! end do
597 : ! ABI_FREE(pfactors)
598 : ! ABI_FREE(powers)
599 : ! end if
600 :
601 : DBG_EXIT("COLL")
602 :
603 579 : end subroutine setmesh
604 : !!***
605 :
606 : !----------------------------------------------------------------------
607 :
608 : !!****f* m_fft_mesh/check_rot_fft
609 : !! NAME
610 : !! check_rot_fft
611 : !!
612 : !! FUNCTION
613 : !! Return .TRUE. if the given grid in real space is compatible
614 : !! with the rotational part of the space group symmetries.
615 : !!
616 : !! INPUTS
617 : !! nsym=Number of symmetry operations
618 : !! symrel(3,3,nsym)=Symmetry operations in real space.
619 : !! nr1,nr2,nr3=FFT divisions.
620 : !!
621 : !! SOURCE
622 :
623 4043 : pure function check_rot_fft(nsym,symrel,nr1,nr2,nr3)
624 :
625 : !Arguments
626 : !Scalar
627 : integer,intent(in) :: nr1,nr2,nr3,nsym
628 : logical :: check_rot_fft
629 : !Arrays
630 : integer,intent(in) :: symrel(3,3,nsym)
631 :
632 : !local variables
633 : integer :: is
634 : !************************************************************************
635 :
636 : ! The grid is compatible with the symmetries (only rotational part) if
637 : ! for each symmetry, each n_i and n_j ==> $n_i*R_{ij}/n_j$ is an integer
638 4043 : check_rot_fft=.TRUE.
639 138631 : do is=1,nsym
640 : if (MOD(symrel(2,1,is)*nr2, nr1) /=0 .or. &
641 : MOD(symrel(3,1,is)*nr3, nr1) /=0 .or. &
642 : MOD(symrel(1,2,is)*nr1, nr2) /=0 .or. &
643 : MOD(symrel(3,2,is)*nr3, nr2) /=0 .or. &
644 135536 : MOD(symrel(1,3,is)*nr1, nr3) /=0 .or. &
645 : MOD(symrel(2,3,is)*nr2, nr3) /=0 &
646 3095 : ) then
647 : check_rot_fft=.FALSE.; EXIT
648 : end if
649 : end do
650 :
651 4043 : end function check_rot_fft
652 : !!***
653 :
654 : !----------------------------------------------------------------------
655 :
656 : !!****f* m_fft_mesh/fft_check_rotrans
657 : !! NAME
658 : !! fft_check_rotrans
659 : !!
660 : !! FUNCTION
661 : !! Checks if the real space FFT mesh is compatible both with the rotational
662 : !! and the translational part of space group of the crystal.
663 : !!
664 : !! INPUTS
665 : !! nsym=Number of symmetries.
666 : !! symrel(3,3,nsym)=Symmetries in real space in reduced coordinates.
667 : !! tnons(3,nsym)=Fractional translations.
668 : !! ngfft(18)=Information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
669 : !!
670 : !! OUTPUT
671 : !! err(3,nsym)=The max error for each symmetry. (given in terms of the FFT vectors)
672 : !! isok=.FALSE. if the FFT mesh does not fulfil all symmetry properties of the crystal.
673 : !!
674 : !! SOURCE
675 :
676 0 : function fft_check_rotrans(nsym,symrel,tnons,ngfft,err) result(isok)
677 :
678 : !Arguments ------------------------------------
679 : !scalars
680 : integer,intent(in) :: nsym
681 : logical :: isok
682 : !arrays
683 : integer,intent(in) :: symrel(3,3,nsym)
684 : integer,intent(in) :: ngfft(18)
685 : real(dp),intent(in) :: tnons(3,nsym)
686 : real(dp),intent(out) :: err(3,nsym)
687 :
688 : !Local variables-------------------------------
689 : !scalars
690 : integer :: isym,ix,iy,iz,jx,jy,jz,ngfft1,ngfft2,ngfft3
691 : !character(len=500) :: msg
692 : !arrays
693 0 : integer :: Rm1(3,3,nsym),r1_FFT(3),red2fft(3,3)
694 0 : real(dp) :: Rm1_FFT(3,3,nsym),fft2red(3,3),r2_FFT(3),tnons_FFT(3,nsym)
695 : ! *************************************************************************
696 :
697 : ! Precalculate R^-1 and fractional translations in FFT coordinates
698 0 : ngfft1=ngfft(1)
699 0 : ngfft2=ngfft(2)
700 0 : ngfft3=ngfft(3)
701 :
702 0 : red2fft=RESHAPE([ngfft1,0,0,0,ngfft2,0,0,0,ngfft3], [3,3])
703 0 : fft2red=RESHAPE((/(one/ngfft1),zero,zero,zero,(one/ngfft2),zero,zero,zero, (one/ngfft3)/),(/3,3/))
704 : !
705 : ! === For a fully compatible mesh, each Rm1_FFT should be integer ===
706 0 : do isym=1,nsym
707 0 : call mati3inv(symrel(:,:,isym),Rm1(:,:,isym))
708 0 : Rm1(:,:,isym)=TRANSPOSE(Rm1(:,:,isym))
709 0 : Rm1_FFT(:,:,isym)=MATMUL(MATMUL(red2fft,Rm1(:,:,isym)),fft2red)
710 0 : tnons_FFT(:,isym)=MATMUL(red2fft,tnons(:,isym))
711 : end do
712 :
713 0 : err(:,:)=smallest_real
714 0 : do iz=0,ngfft3-1
715 0 : R1_FFT(3)=DBLE(iz)
716 0 : do iy=0,ngfft2-1
717 0 : R1_FFT(2)=DBLE(iy)
718 0 : do ix=0,ngfft1-1
719 0 : R1_FFT(1)=DBLE(ix)
720 0 : do isym=1,nsym ! Form R^-1 (r-\tau) in the FFT basis.
721 0 : R2_FFT(:)=MATMUL(Rm1_FFT(:,:,isym),R1_FFT(:)-tnons_FFT(:,isym))
722 0 : jx=NINT(R2_FFT(1)); err(1,isym)=MAX(err(1,isym),ABS(R2_FFT(1)-jx)/ngfft1)
723 0 : jy=NINT(R2_FFT(2)); err(2,isym)=MAX(err(2,isym),ABS(R2_FFT(2)-jy)/ngfft2)
724 0 : jz=NINT(R2_FFT(3)); err(3,isym)=MAX(err(3,isym),ABS(R2_FFT(3)-jz)/ngfft3)
725 : end do
726 : end do
727 : end do
728 : end do
729 :
730 : isok=.TRUE.
731 0 : do isym=1,nsym
732 0 : if (ANY(err(:,isym)>tol6)) then
733 0 : isok=.FALSE.
734 : !write(msg,'(a,i3,a,3es14.6)')' symmetry ',isym,') not compatible with FFT grid, error ',err(:,isym)
735 : !ABI_WARNING(msg)
736 : end if
737 : end do
738 :
739 0 : end function fft_check_rotrans
740 : !!***
741 :
742 : !----------------------------------------------------------------------
743 :
744 : !!****f* m_fft_mesh/rotate_FFT_mesh
745 : !! NAME
746 : !! rotate_FFT_mesh
747 : !!
748 : !! FUNCTION
749 : !! Find the FFT index of $ R{-1}(r-\tau) $ for each point in the FFT box.
750 : !! $R$ is a symmetry operation in real space, $\tau$ is the associated
751 : !! fractional translation.
752 : !!
753 : !! INPUTS
754 : !! nsym=Number of symmetries.
755 : !! symrel(3,3,nsym)=Symmetries in real space in reduced coordinates.
756 : !! tnons(3,nsym)=Fractional translations.
757 : !! ngfft(18)=Information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
758 : !!
759 : !! OUTPUT
760 : !! irottb(ngfftot,nsym)=Indices of $R^{-1}(r-\tau)$ in the FFT box.
761 : !! preserve=.FALSE. if the FFT mesh does not fulfil all symmetry properties of the crystal.
762 : !!
763 : !! NOTES
764 : !! The evaluation of the rotated point $R^{-1}(r-\tau)$ is done using real arithmetic.
765 : !! As a consequence, if the FFT mesh does not fulfil the symmetry properties
766 : !! of the crystal, the array irottb will contain the index of the FFT point which
767 : !! is the closest one to $R^{-1}(r-\tau)$. This might lead to inaccuracies in the
768 : !! final results, in particular in the description of degenerate states.
769 : !!
770 : !! SOURCE
771 :
772 2423 : subroutine rotate_fft_mesh(nsym, symrel, tnons, ngfft, irottb, preserve)
773 :
774 : !Arguments ------------------------------------
775 : !scalars
776 : integer,intent(in) :: nsym
777 : logical,intent(out) :: preserve
778 : !arrays
779 : integer,intent(in) :: symrel(3,3,nsym)
780 : integer,intent(in) :: ngfft(18)
781 : integer,intent(out) :: irottb(ngfft(1)*ngfft(2)*ngfft(3),nsym)
782 : real(dp),intent(in) :: tnons(3,nsym)
783 :
784 : !Local variables-------------------------------
785 : !scalars
786 : integer :: ir1,isym,ix,iy,iz,jx,jy,jz,ngfft1,ngfft2,ngfft3
787 : !character(len=500) :: msg
788 : !arrays
789 4846 : integer :: Rm1(3,3,nsym),r1_FFT(3),red2fft(3,3)
790 4846 : real(dp) :: Rm1_FFT(3,3,nsym),err(3,nsym),fft2red(3,3),r2_FFT(3), tnons_FFT(3,nsym)
791 : ! *************************************************************************
792 :
793 : ! Precalculate R^-1 and fractional translations in FFT coordinates.
794 2423 : ngfft1 = ngfft(1); ngfft2 = ngfft(2); ngfft3 = ngfft(3)
795 :
796 24230 : red2fft = reshape([ngfft1, 0, 0, 0, ngfft2, 0, 0, 0, ngfft3], [3, 3])
797 24230 : fft2red = reshape([(one/ngfft1), zero, zero, zero, one/ngfft2, zero, zero, zero, one/ngfft3], [3, 3])
798 :
799 : ! For a fully compatible mesh, each Rm1_FFT should be integer ===
800 104839 : do isym=1,nsym
801 102416 : call mati3inv(symrel(:,:,isym), Rm1(:,:,isym))
802 2560400 : Rm1(:,:,isym) = transpose(Rm1(:,:,isym))
803 10548848 : Rm1_FFT(:,:,isym) = matmul(matmul(red2fft, Rm1(:,:,isym)), fft2red)
804 2870071 : tnons_FFT(:,isym) = matmul(red2fft, tnons(:,isym))
805 : end do
806 :
807 412087 : err(:,:) = zero
808 :
809 : !$OMP PARALLEL DO PRIVATE(R1_FFT,ir1,R2_FFT,jx,jy,jz) reduction(MAX:err)
810 43992 : do iz=0,ngfft3-1
811 41569 : R1_FFT(3) = dble(iz)
812 865654 : do iy=0,ngfft2-1
813 821662 : R1_FFT(2) = dble(iy)
814 20037571 : do ix=0,ngfft1-1
815 19174340 : R1_FFT(1) = dble(ix)
816 19174340 : ir1 = 1+ix+iy*ngfft1+iz*ngfft1*ngfft2
817 687292558 : do isym=1,nsym
818 : ! Form R^-1 (r-\tau) in the FFT basis.
819 10676744896 : R2_FFT(:) = MATMUL(Rm1_FFT(:,:,isym),R1_FFT(:)-tnons_FFT(:,isym))
820 667296556 : jx = NINT(R2_FFT(1)); err(1,isym) = MAX(err(1,isym), ABS(R2_FFT(1)-jx)/ngfft1)
821 667296556 : jy = NINT(R2_FFT(2)); err(2,isym) = MAX(err(2,isym), ABS(R2_FFT(2)-jy)/ngfft2)
822 667296556 : jz = NINT(R2_FFT(3)); err(3,isym) = MAX(err(3,isym), ABS(R2_FFT(3)-jz)/ngfft3)
823 667296556 : jx = MODULO(jx, ngfft1)
824 667296556 : jy = MODULO(jy, ngfft2)
825 667296556 : jz = MODULO(jz, ngfft3)
826 686470896 : irottb(ir1,isym) = 1+jx+jy*ngfft1+jz*ngfft1*ngfft2
827 : end do
828 : end do
829 : end do
830 : end do
831 :
832 2423 : preserve = .TRUE.
833 104839 : do isym=1,nsym
834 403431 : if (any(err(:,isym) > tol6)) then
835 2912 : preserve = .FALSE.
836 : !write(msg,'(a,i0,a,3es14.6)')' symmetry ',isym,') not compatible with FFT grid, error ',err(:,isym)
837 : !ABI_WARNING(msg)
838 : end if
839 : end do
840 :
841 2423 : end subroutine rotate_fft_mesh
842 : !!***
843 :
844 : !----------------------------------------------------------------------
845 :
846 : !!****f* m_numeric_tools/denpot_project
847 : !! NAME
848 : !!
849 : !! FUNCTION
850 : !! Compute n(r) + n( $R^{-1}(r-\tau)$ in) / 2
851 : !! Mainly used with R = inversion to select the even/odd part under inversion
852 : !!
853 : !! INPUTS
854 : !! cplex=1 for real, 2 for complex data.
855 : !! ngfft(3)=Mesh divisions of input array
856 : !! nspden=Number of density components.
857 : !! in_rhor(cplex * nfftot * nspden)=Input array
858 : !! one_symrel(3,3)= R operation
859 : !! tau(3)=Fractional translation.
860 : !!
861 : !! OUTPUT
862 : !! out_rhor(cplex * nfftot * nspden)=Output array
863 : !!
864 : !! SOURCE
865 :
866 0 : subroutine denpot_project(cplex, ngfft, nspden, in_rhor, one_symrel, one_tnons, out_rhor)
867 :
868 : !Arguments-------------------------------------------------------------
869 : !scalars
870 : integer,intent(in) :: cplex, nspden
871 : !arrays
872 : integer,intent(in) :: ngfft(18), one_symrel(3,3)
873 : real(dp),intent(in) :: in_rhor(cplex, product(ngfft(1:3)), nspden)
874 : real(dp),intent(in) :: one_tnons(3)
875 : real(dp),intent(out) :: out_rhor(cplex, product(ngfft(1:3)), nspden)
876 :
877 : !Local variables--------------------------------------------------------
878 : !scalars
879 : integer,parameter :: nsym1 = 1, isgn = 1
880 : integer :: ispden, ii, ifft, ifft_rot, nfft
881 : logical :: preserve
882 : !arrays
883 0 : integer,allocatable :: irottb(:)
884 : ! *************************************************************************
885 :
886 0 : nfft = product(ngfft(1:3))
887 0 : ABI_MALLOC(irottb, (nfft))
888 :
889 0 : call rotate_fft_mesh(nsym1, one_symrel, one_tnons, ngfft, irottb, preserve)
890 0 : ABI_CHECK(preserve, "FFT mesh is not compatible with {R, tau}")
891 :
892 0 : do ispden=1,nspden
893 0 : do ifft=1,nfft
894 0 : ifft_rot = irottb(ifft)
895 0 : do ii=1,cplex
896 0 : out_rhor(cplex, ifft, ispden) = (in_rhor(cplex, ifft, ispden) + isgn * in_rhor(cplex, ifft_rot, ispden)) * half
897 : end do
898 : end do
899 : end do
900 :
901 0 : ABI_FREE(irottb)
902 :
903 0 : end subroutine denpot_project
904 : !!***
905 :
906 : !----------------------------------------------------------------------
907 :
908 : !!****f* m_fft_mesh/cigfft
909 : !! NAME
910 : !! cigfft
911 : !!
912 : !! FUNCTION
913 : !! For each of the (2*nG0sh+1)**3 vectors G0 around the origin,
914 : !! calculate G-G0 and its FFT index number for all the NPWVEC vectors G.
915 : !!
916 : !! INPUTS
917 : !! mG0(3)= For each reduced direction gives the max G0 component to account for umklapp processes.
918 : !! npwvec=Number of plane waves
919 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
920 : !! gvec(3,npwvec)=Reduced coordinates of G vectors.
921 : !!
922 : !! OUTPUT
923 : !! igfft(npwvec,2*mG0(1)+1,2*mG0(2)+1,2*mG0(3)+1)=For each G, and each G0 vector,
924 : !! it gives the FFT grid index of the G-G0 vector.
925 : !! ierr=Number of G-G0 vectors falling outside the inout FFT box.
926 : !!
927 : !! SOURCE
928 :
929 0 : subroutine cigfft(mG0,npwvec,ngfft,gvec,igfft,ierr)
930 :
931 : !Arguments ------------------------------------
932 : !scalars
933 : integer,intent(in) :: npwvec
934 : integer,intent(out) :: ierr
935 : !arrays
936 : integer,intent(in) :: gvec(3,npwvec)
937 : integer,intent(in) :: mg0(3),ngfft(18)
938 : integer,intent(out) :: igfft(npwvec,2*mg0(1)+1,2*mg0(2)+1,2*mg0(3)+1)
939 :
940 : !Local variables ------------------------------
941 : !scalars
942 : integer :: gmg01,gmg02,gmg03,ig,ig01,ig02,ig03,n1,n2,n3
943 : character(len=500) :: msg
944 : !arrays
945 : integer :: gmg0(3)
946 : !************************************************************************
947 :
948 : DBG_ENTER("COLL")
949 :
950 0 : if (ANY(mg0<0)) then
951 0 : ABI_BUG(sjoin('Found negative value of mg0:', trim(ltoa(mg0))))
952 : end if
953 :
954 0 : n1=ngfft(1)
955 0 : n2=ngfft(2)
956 0 : n3=ngfft(3)
957 0 : ierr=0
958 :
959 0 : do ig=1,npwvec
960 0 : do ig01=-mg0(1),mG0(1)
961 0 : gmg0(1) = gvec(1,ig)-ig01
962 0 : do ig02=-mg0(2),mg0(2)
963 0 : gmg0(2) = gvec(2,ig)-ig02
964 0 : do ig03=-mg0(3),mg0(3)
965 0 : gmg0(3) = gvec(3,ig)-ig03
966 : ! Calculate FFT index of G-G0. Consider possible wrap around errors.
967 0 : gmg01=MODULO(gmg0(1),n1)
968 0 : gmg02=MODULO(gmg0(2),n2)
969 0 : gmg03=MODULO(gmg0(3),n3)
970 0 : igfft(ig,ig01+mg0(1)+1,ig02+mg0(2)+1,ig03+mg0(3)+1) = 1+gmg01+gmg02*n1+gmg03*n1*n2
971 0 : if ( ANY(gmg0>ngfft(1:3)/2) .or. ANY(gmg0<-(ngfft(1:3)-1)/2) ) then
972 0 : igfft(ig,ig01+mg0(1)+1,ig02+mg0(2)+1,ig03+mg0(3)+1) = 0
973 0 : ierr=ierr+1
974 : end if
975 : end do
976 : end do
977 : end do
978 : end do !ig
979 :
980 0 : if (ierr/=0) then
981 : write(msg,'(a,i0,3a)')&
982 0 : 'Found ',ierr,' G-G0 vectors falling outside the FFT box. ',ch10,&
983 0 : 'igfft will be set to zero for these particular G-G0 '
984 0 : ABI_WARNING(msg)
985 : end if
986 :
987 : DBG_EXIT("COLL")
988 :
989 0 : end subroutine cigfft
990 : !!***
991 :
992 : !----------------------------------------------------------------------
993 :
994 : !!****f* m_fft_mesh/ig2gfft
995 : !! NAME
996 : !! ig2gfft
997 : !!
998 : !! FUNCTION
999 : !! Return the reduced component of a G-vector in the FFT mesh starting from is index.
1000 : !!
1001 : !! INPUTS
1002 : !! ig = The index >=1, <=ng
1003 : !! ng = The number of FFT points along this direction.
1004 : !!
1005 : !! OUTPUT
1006 : !! gc = The reduced component
1007 : !!
1008 : !! SOURCE
1009 :
1010 117932501 : elemental integer function ig2gfft(ig, ng) result (gc)
1011 :
1012 : !Arguments ------------------------------------
1013 : !scalars
1014 : integer,intent(in) :: ig,ng
1015 : !************************************************************************
1016 :
1017 : ! Use the following indexing (N means ngfft of the adequate direction)
1018 : ! 0 1 2 3 ... N/2 -(N-1)/2 ... -1 <= gc
1019 : ! 1 2 3 4 ....N/2+1 N/2+2 ... N <= index ig
1020 : !
1021 117932501 : if (ig <= 0 .or. ig > ng) then
1022 : ! Wrong ig, returns huge. Parent code will likely crash with SIGSEV.
1023 117932501 : gc = huge(1)
1024 : return
1025 : end if
1026 :
1027 117932501 : if ( ig > ng/2 + 1) then
1028 53327880 : gc = ig - ng -1
1029 : else
1030 64604621 : gc = ig -1
1031 : end if
1032 :
1033 : end function ig2gfft
1034 : !!***
1035 :
1036 : !----------------------------------------------------------------------
1037 :
1038 : !!****f* m_fft_mesh/g2ifft
1039 : !! NAME
1040 : !! g2ifft
1041 : !!
1042 : !! FUNCTION
1043 : !! Returns the index of G in the FFT box from its reduced coordinates. 0 if not in the BOX.
1044 : !!
1045 : !! INPUTS
1046 : !! gg(3)=Reduced coordinated of the G vector.
1047 : !! ngfft(18) = Info on the FFT box.
1048 : !!
1049 : !! OUTPUT
1050 : !! gidx=Index in the FFT box. 0 if G is outside the box.
1051 : !!
1052 : !! SOURCE
1053 :
1054 52524 : pure integer function g2ifft(gg,ngfft) result (gidx)
1055 :
1056 : !Arguments ------------------------------------
1057 : integer,intent(in) :: gg(3),ngfft(3)
1058 :
1059 : !Local variables-------------------------------
1060 : integer :: n1,n2,n3,ig1,ig2,ig3
1061 : !************************************************************************
1062 :
1063 : ! Use the following indexing (N means ngfft of the adequate direction)
1064 : ! 0 1 2 3 ... N/2 -(N-1)/2 ... -1 <= gg
1065 : ! 1 2 3 4 ....N/2+1 N/2+2 ... N <= index
1066 : !
1067 420192 : if (any(gg > ngfft(1:3)/2) .or. any(gg < -(ngfft(1:3)-1)/2)) then ! out of the box.
1068 : gidx = 0
1069 : else
1070 52524 : n1 = ngfft(1)
1071 52524 : n2 = ngfft(2)
1072 52524 : n3 = ngfft(3)
1073 52524 : ig1 = MODULO(gg(1),n1)
1074 52524 : ig2 = MODULO(gg(2),n2)
1075 52524 : ig3 = MODULO(gg(3),n3)
1076 52524 : gidx = 1 + ig1 + n1*(ig2+ig3*n2)
1077 : end if
1078 :
1079 52524 : end function g2ifft
1080 : !!***
1081 :
1082 : !----------------------------------------------------------------------
1083 :
1084 : !!****f* m_fft_mesh/get_gfft
1085 : !! NAME
1086 : !! get_gfft
1087 : !!
1088 : !! FUNCTION
1089 : !! Returns the set of G-vectors in the FFT mesh and the maximal kinetic energy of k+G.
1090 : !!
1091 : !! INPUTS
1092 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1093 : !! kpt(3)=input k vector (reduced coordinates --in terms of reciprocal lattice primitive translations)
1094 : !! gmet(3,3)=reciprocal space metric (bohr^-2)
1095 : !!
1096 : !! OUTPUT
1097 : !! gsq_max=Max value of (k+G)^2 for G in the FFT box
1098 : !! gfft(3,nfft_tot) = The reduced components of the G in the FFT mesh (nfft_tot=PRODUCT(ngfft(1:3))
1099 : !!
1100 : !! SOURCE
1101 :
1102 22248 : pure subroutine get_gfft(ngfft, kpt, gmet, gsq_max, gfft)
1103 :
1104 : !Arguments ------------------------------------
1105 : !scalars
1106 : real(dp),intent(out) :: gsq_max
1107 : !arrays
1108 : integer,intent(in) :: ngfft(18)
1109 : integer,intent(out) :: gfft(3,ngfft(1)*ngfft(2)*ngfft(3))
1110 : real(dp),intent(in) :: kpt(3), gmet(3,3)
1111 :
1112 : !Local variables-------------------------------
1113 : integer :: ifft,g1,g2,g3,i1,i2,i3
1114 : real(dp) :: dsq
1115 : !************************************************************************
1116 :
1117 22248 : ifft=0; gsq_max=smallest_real
1118 333022 : do i3=1,ngfft(3)
1119 310774 : g3 = ig2gfft(i3,ngfft(3))
1120 5702980 : do i2=1,ngfft(2)
1121 5369958 : g2 = ig2gfft(i2,ngfft(2))
1122 117930158 : do i1=1,ngfft(1)
1123 112249426 : g1 = ig2gfft(i1,ngfft(1))
1124 112249426 : ifft = ifft+1
1125 448997704 : gfft(:,ifft) = [g1, g2, g3]
1126 : dsq=gmet(1,1)*(kpt(1)+dble(i1))**2 &
1127 : +gmet(2,2)*(kpt(2)+dble(i2))**2 &
1128 : +gmet(3,3)*(kpt(3)+dble(i3))**2 &
1129 : +2._dp*(gmet(1,2)*(kpt(1)+dble(i1))*(kpt(2)+dble(i2)) &
1130 : +gmet(2,3)*(kpt(2)+dble(i2))*(kpt(3)+dble(i3)) &
1131 112249426 : +gmet(3,1)*(kpt(3)+dble(i3))*(kpt(1)+dble(i1)))
1132 117619384 : gsq_max = MAX(dsq,gsq_max)
1133 : end do
1134 : end do
1135 : end do
1136 :
1137 22248 : end subroutine get_gfft
1138 : !!***
1139 :
1140 : !----------------------------------------------------------------------
1141 :
1142 : !!****f* m_fft_mesh/calc_ceigr_spc
1143 : !! NAME
1144 : !! calc_ceigr_spc
1145 : !!
1146 : !! FUNCTION
1147 : !! Helper function to calculate e^{iG.r} on the FFT mesh.
1148 : !!
1149 : !! INPUTS
1150 : !! gg(3)=G vector in reduced coordinates.
1151 : !! nfft=Total number of points in the FFT mesh.
1152 : !! nspinor=Number of spinors
1153 : !! ngfft(18)=information about 3D FFT,
1154 : !!
1155 : !! OUTPUT
1156 : !! ceigr(nfft*nspinor)=e^{ik.r} on the FFT mesh.
1157 : !!
1158 : !! SOURCE
1159 :
1160 0 : subroutine calc_ceigr_spc(gg, nfft, nspinor, ngfft, ceigr)
1161 :
1162 : !Arguments ------------------------------------
1163 : !scalars
1164 : integer,intent(in) :: nfft,nspinor
1165 : !arrays
1166 : integer,intent(in) :: gg(3)
1167 : integer,intent(in) :: ngfft(18)
1168 : complex(sp),intent(out) :: ceigr(nfft*nspinor)
1169 :
1170 : !Local variables-------------------------------
1171 : integer :: ix,iy,iz,ifft,base,isp
1172 : real(dp) :: gdotr
1173 : ! *************************************************************************
1174 :
1175 0 : if (ALL(gg==0)) then
1176 0 : ceigr=(1._sp,0._sp)
1177 : RETURN
1178 : end if
1179 :
1180 0 : do iz=0,ngfft(3)-1
1181 0 : do iy=0,ngfft(2)-1
1182 0 : do ix=0,ngfft(1)-1
1183 : gdotr= two_pi*( gg(1)*(ix/DBLE(ngfft(1))) &
1184 : +gg(2)*(iy/DBLE(ngfft(2))) &
1185 0 : +gg(3)*(iz/DBLE(ngfft(3))) )
1186 0 : ifft = ix + ngfft(1) * (iy + ngfft(2) * iz) + 1
1187 0 : ceigr(ifft)=CMPLX(DCOS(gdotr),DSIN(gdotr), KIND=sp)
1188 : end do
1189 : end do
1190 : end do
1191 :
1192 0 : if (nspinor > 1) then
1193 0 : do isp=2,nspinor
1194 0 : base = 1 + (isp-1)*nfft
1195 0 : call xcopy(nfft,ceigr,1,ceigr(base:),1)
1196 : end do
1197 : end if
1198 :
1199 : end subroutine calc_ceigr_spc
1200 : !!***
1201 :
1202 : !----------------------------------------------------------------------
1203 :
1204 : !!****f* m_fft_mesh/calc_ceigr_dpc
1205 : !! NAME
1206 : !! calc_ceigr_dpc
1207 : !!
1208 : !! FUNCTION
1209 : !! Helper function to calculate e^{iG.r} on the FFT mesh.
1210 : !!
1211 : !! INPUTS
1212 : !! gg(3)=G vector in reduced coordinates.
1213 : !! nfft=Total number of points in the FFT mesh.
1214 : !! nspinor=Number of spinors
1215 : !! ngfft(18)=information about 3D FFT,
1216 : !!
1217 : !! OUTPUT
1218 : !! ceigr(nfft*nspinor)=e^{ik.r} on the FFT mesh.
1219 : !!
1220 : !! SOURCE
1221 :
1222 901 : subroutine calc_ceigr_dpc(gg, nfft, nspinor, ngfft, ceigr)
1223 :
1224 : !Arguments ------------------------------------
1225 : !scalars
1226 : integer,intent(in) :: nfft,nspinor
1227 : !arrays
1228 : integer,intent(in) :: gg(3)
1229 : integer,intent(in) :: ngfft(18)
1230 : complex(dp),intent(out) :: ceigr(nfft*nspinor)
1231 :
1232 : !Local variables-------------------------------
1233 : integer :: ix,iy,iz,ifft,base,isp
1234 : real(dp) :: gdotr
1235 : ! *************************************************************************
1236 :
1237 1621 : if (ALL(gg==0)) then
1238 0 : ceigr=cone; RETURN
1239 : end if
1240 :
1241 16279 : do iz=0,ngfft(3)-1
1242 348907 : do iy=0,ngfft(2)-1
1243 13462254 : do ix=0,ngfft(1)-1
1244 : gdotr= two_pi*( gg(1)*(ix/DBLE(ngfft(1))) &
1245 : +gg(2)*(iy/DBLE(ngfft(2))) &
1246 13114248 : +gg(3)*(iz/DBLE(ngfft(3))) )
1247 13114248 : ifft = ix + ngfft(1) * (iy + ngfft(2) * iz) + 1
1248 13446876 : ceigr(ifft)=DCMPLX(DCOS(gdotr),DSIN(gdotr))
1249 : end do
1250 : end do
1251 : end do
1252 :
1253 901 : if (nspinor > 1) then
1254 0 : do isp=2,nspinor
1255 0 : base = 1 + (isp-1)*nfft
1256 0 : call xcopy(nfft,ceigr,1,ceigr(base:),1)
1257 : end do
1258 : end if
1259 :
1260 : end subroutine calc_ceigr_dpc
1261 : !!***
1262 :
1263 : !----------------------------------------------------------------------
1264 :
1265 : !!****f* m_fft_mesh/calc_eigr
1266 : !! NAME
1267 : !! calc_eigr
1268 : !!
1269 : !! FUNCTION
1270 : !! Helper function to calculate e^{iG.r} on the FFT mesh.
1271 : !!
1272 : !! INPUTS
1273 : !! gg(3)=G vector in reduced coordinates.
1274 : !! nfft=Total number of points in the FFT mesh.
1275 : !! ngfft(18)=information about 3D FFT,
1276 : !!
1277 : !! OUTPUT
1278 : !! eigr(2*nfft)=e^{ig.r} on the FFT mesh.
1279 : !!
1280 : !! SOURCE
1281 :
1282 3 : pure subroutine calc_eigr(gg, nfft, ngfft, eigr)
1283 :
1284 : !Arguments ------------------------------------
1285 : !scalars
1286 : integer,intent(in) :: nfft
1287 : !arrays
1288 : integer,intent(in) :: gg(3)
1289 : integer,intent(in) :: ngfft(18)
1290 : real(dp),intent(out) :: eigr(2*nfft)
1291 :
1292 : !Local variables-------------------------------
1293 : integer :: ix,iy,iz,ifft
1294 : real(dp) :: gdotr
1295 : ! *************************************************************************
1296 :
1297 3 : if (ALL(gg==0)) then
1298 0 : eigr(1:2*nfft:2)=one
1299 0 : eigr(2:2*nfft:2)=zero
1300 : RETURN
1301 : end if
1302 :
1303 3 : ifft=1
1304 273 : do iz=0,ngfft(3)-1
1305 24573 : do iy=0,ngfft(2)-1
1306 2211570 : do ix=0,ngfft(1)-1
1307 : gdotr= two_pi*( gg(1)*(ix/DBLE(ngfft(1))) &
1308 : +gg(2)*(iy/DBLE(ngfft(2))) &
1309 2187000 : +gg(3)*(iz/DBLE(ngfft(3))) )
1310 2187000 : eigr(ifft )=DCOS(gdotr)
1311 2187000 : eigr(ifft+1)=DSIN(gdotr)
1312 2211300 : ifft = ifft + 2
1313 : end do
1314 : end do
1315 : end do
1316 :
1317 : end subroutine calc_eigr
1318 : !!***
1319 :
1320 : !----------------------------------------------------------------------
1321 :
1322 : !!****f* m_fft_mesh/calc_ceikr_dpc
1323 : !! NAME
1324 : !! calc_ceikr_dpc
1325 : !!
1326 : !! FUNCTION
1327 : !! calculate e^{ik.r} on the fft mesh.
1328 : !!
1329 : !! INPUTS
1330 : !! kk(3)=k-point in reduced coordinates.
1331 : !! nfft=total number of points in the fft mesh.
1332 : !! ngfft(18)=information about 3d fft,
1333 : !! nspinor=number of spinor components.
1334 : !!
1335 : !! OUTPUT
1336 : !! ceikr(nfft*nspinor) = e^{ik.r} on the fft mesh.
1337 : !!
1338 : !! SOURCE
1339 :
1340 0 : subroutine calc_ceikr_dpc(kk, ngfft, nfft, nspinor, ceikr)
1341 :
1342 : !arguments ------------------------------------
1343 : !scalars
1344 : integer,intent(in) :: nfft, nspinor
1345 : !arrays
1346 : real(dp),intent(in) :: kk(3)
1347 : integer,intent(in) :: ngfft(18)
1348 : complex(dp),intent(out) :: ceikr(nfft*nspinor)
1349 :
1350 : !local variables-------------------------------
1351 : integer :: ix, iy, iz, ifft
1352 : real(dp) :: kdotr
1353 : ! *************************************************************************
1354 :
1355 0 : if (all(abs(kk) < tol12)) then
1356 0 : ceikr = cone; return
1357 : end if
1358 :
1359 : !$OMP PARALLEL DO PRIVATE(kdotr, ifft)
1360 0 : do iz=0,ngfft(3)-1
1361 0 : do iy=0,ngfft(2)-1
1362 0 : do ix=0,ngfft(1)-1
1363 : kdotr = two_pi*( kk(1) * (ix / dble(ngfft(1))) &
1364 : +kk(2) * (iy / dble(ngfft(2))) &
1365 0 : +kk(3) * (iz / dble(ngfft(3))) )
1366 0 : ifft = ix + ngfft(1) * (iy + ngfft(2) * iz) + 1
1367 0 : ceikr(ifft) = dcmplx(cos(kdotr), sin(kdotr))
1368 : end do
1369 : end do
1370 : end do
1371 :
1372 0 : if (nspinor > 1) ceikr(nfft+1:) = ceikr(1:nfft)
1373 :
1374 : end subroutine calc_ceikr_dpc
1375 : !!***
1376 :
1377 : !----------------------------------------------------------------------
1378 :
1379 : !!****f* m_fft_mesh/calc_ceikr_spc
1380 : !! name
1381 : !! calc_ceikr_spc
1382 : !!
1383 : !! function
1384 : !! calculate e^{ik.r} on the fft mesh.
1385 : !!
1386 : !! inputs
1387 : !! kk(3)=k-point in reduced coordinates.
1388 : !! nfft=total number of points in the fft mesh.
1389 : !! ngfft(18)=information about 3d fft,
1390 : !! nspinor=number of spinor components.
1391 : !!
1392 : !! output
1393 : !! ceikr(nfft*nspinor) = e^{ik.r} on the fft mesh.
1394 : !!
1395 : !! source
1396 :
1397 0 : subroutine calc_ceikr_spc(kk, ngfft, nfft, nspinor, ceikr)
1398 :
1399 : !arguments ------------------------------------
1400 : !scalars
1401 : integer,intent(in) :: nfft, nspinor
1402 : !arrays
1403 : real(dp),intent(in) :: kk(3)
1404 : integer,intent(in) :: ngfft(18)
1405 : complex(sp),intent(out) :: ceikr(nfft*nspinor)
1406 :
1407 : !local variables-------------------------------
1408 : integer :: ix, iy, iz, ifft
1409 : real(dp) :: kdotr
1410 : ! *************************************************************************
1411 :
1412 0 : if (all(abs(kk) < tol12)) then
1413 0 : ceikr = cone; return
1414 : end if
1415 :
1416 : !$OMP PARALLEL DO PRIVATE(kdotr, ifft)
1417 0 : do iz=0,ngfft(3)-1
1418 0 : do iy=0,ngfft(2)-1
1419 0 : do ix=0,ngfft(1)-1
1420 : kdotr = two_pi*( kk(1) * (ix / dble(ngfft(1))) &
1421 : +kk(2) * (iy / dble(ngfft(2))) &
1422 0 : +kk(3) * (iz / dble(ngfft(3))) )
1423 0 : ifft = ix + ngfft(1) * (iy + ngfft(2) * iz) + 1
1424 0 : ceikr(ifft) = cmplx(cos(kdotr), sin(kdotr), kind=sp)
1425 : end do
1426 : end do
1427 : end do
1428 :
1429 0 : if (nspinor > 1) ceikr(nfft+1:) = ceikr(1:nfft)
1430 :
1431 : end subroutine calc_ceikr_spc
1432 : !!***
1433 :
1434 : !----------------------------------------------------------------------
1435 :
1436 : !!****f* m_fft_mesh/times_eigr
1437 : !! NAME
1438 : !! times_eigr
1439 : !!
1440 : !! FUNCTION
1441 : !! Multiply an array on the real-space mesh by e^{iG0.r} where G0 is a reciprocal lattice vector.
1442 : !!
1443 : !! INPUTS
1444 : !! gg(3)=G vector in reduced coordinates.
1445 : !! ngfft(18)=information about 3D FFT,
1446 : !! nfft=Number of points in the FFT mesh.
1447 : !! ndat=Number of arrays
1448 : !!
1449 : !! SIDE EFFECTS
1450 : !! ur(2,nfft,ndat)= contains u(r) in input. output: u(r) e^{ig.r} on the real-space FFT mesh.
1451 : !!
1452 : !! SOURCE
1453 :
1454 13583 : subroutine times_eigr(gg, ngfft, nfft, ndat, ur)
1455 :
1456 : !Arguments ------------------------------------
1457 : !scalars
1458 : integer,intent(in) :: nfft,ndat
1459 : !arrays
1460 : integer,intent(in) :: gg(3)
1461 : integer,intent(in) :: ngfft(18)
1462 : real(dp),intent(inout) :: ur(2,nfft,ndat)
1463 :
1464 : !Local variables-------------------------------
1465 : !scalars
1466 : integer :: ix,iy,iz,ifft,idat
1467 : real(dp) :: gr
1468 : !arrays
1469 : real(dp) :: ph(2),val(2)
1470 : ! *************************************************************************
1471 :
1472 40628 : if (all(gg == 0)) return
1473 :
1474 : !$OMP PARALLEL DO PRIVATE(ifft, gr, ph, val)
1475 10552 : do idat=1,ndat
1476 110118 : do iz=0,ngfft(3)-1
1477 2287818 : do iy=0,ngfft(2)-1
1478 56664066 : do ix=0,ngfft(1)-1
1479 54381524 : ifft = ix + ngfft(1) * (iy + ngfft(2) * iz) + 1
1480 : gr = two_pi*(gg(1)*(ix/dble(ngfft(1))) &
1481 : +gg(2)*(iy/dble(ngfft(2))) &
1482 54381524 : +gg(3)*(iz/dble(ngfft(3))) )
1483 54381524 : ph(1) = cos(gr); ph(2) = sin(gr)
1484 54381524 : val(1) = ur(1,ifft,idat); val(2) = ur(2,ifft,idat)
1485 :
1486 54381524 : ur(1,ifft,idat) = ph(1) * val(1) - ph(2) * val(2)
1487 56564500 : ur(2,ifft,idat) = ph(1) * val(2) + ph(2) * val(1)
1488 : end do
1489 : end do
1490 : end do
1491 : end do ! idat
1492 :
1493 : end subroutine times_eigr
1494 : !!***
1495 :
1496 : !----------------------------------------------------------------------
1497 :
1498 : !!****f* m_fft_mesh/times_eikr
1499 : !! NAME
1500 : !! times_eikr
1501 : !!
1502 : !! FUNCTION
1503 : !! Multiply an array on the real-space mesh by e^{ik.r} where k
1504 : !! is a real(dp) vector in reduced coordinates
1505 : !!
1506 : !! INPUTS
1507 : !! kk(3)=k-vector in reduced coordinates.
1508 : !! ngfft(18)=information about 3D FFT,
1509 : !! nfft=Number of points in the FFT mesh.
1510 : !! ndat=Number of arrays to transform
1511 : !!
1512 : !! SIDE EFFECTS
1513 : !! ur(2,nfft,ndat)= contains u(r) in input. output: u(r) e^{ig.r} on the real-space FFT mesh.
1514 : !!
1515 : !! SOURCE
1516 :
1517 44036 : subroutine times_eikr(kk, ngfft, nfft, ndat, ur)
1518 :
1519 : !Arguments ------------------------------------
1520 : !scalars
1521 : integer,intent(in) :: nfft,ndat
1522 : !arrays
1523 : real(dp),intent(in) :: kk(3)
1524 : integer,intent(in) :: ngfft(18)
1525 : real(dp),intent(inout) :: ur(2,nfft,ndat)
1526 :
1527 : !Local variables-------------------------------
1528 : integer :: ix,iy,iz,ifft,idat
1529 : real(dp) :: kr, ph(2),val(2)
1530 : ! *************************************************************************
1531 :
1532 53949 : if (all(abs(kk) < tol12)) return
1533 :
1534 : !$OMP PARALLEL DO PRIVATE(ifft, kr, ph, val) IF (ndat > 1)
1535 93338 : do idat=1,ndat
1536 854164 : do iz=0,ngfft(3)-1
1537 15123966 : do iy=0,ngfft(2)-1
1538 341073634 : do ix=0,ngfft(1)-1
1539 326000018 : ifft = ix + ngfft(1) * (iy + ngfft(2) * iz) + 1
1540 : kr = two_pi*(kk(1)*(ix/dble(ngfft(1))) &
1541 : +kk(2)*(iy/dble(ngfft(2))) &
1542 326000018 : +kk(3)*(iz/dble(ngfft(3))) )
1543 326000018 : ph(1) = cos(kr); ph(2) = sin(kr)
1544 326000018 : val(1) = ur(1,ifft,idat); val(2) = ur(2,ifft,idat)
1545 :
1546 326000018 : ur(1,ifft,idat) = ph(1) * val(1) - ph(2) * val(2)
1547 340312808 : ur(2,ifft,idat) = ph(1) * val(2) + ph(2) * val(1)
1548 : end do
1549 : end do
1550 : end do
1551 : end do
1552 :
1553 : end subroutine times_eikr
1554 : !!***
1555 :
1556 : ! Version for double precision complex arrays.
1557 0 : subroutine ctimes_eikr(kk, ngfft, nfft, ndat, ur)
1558 :
1559 : !Arguments ------------------------------------
1560 : !scalars
1561 : integer,intent(in) :: nfft,ndat
1562 : !arrays
1563 : real(dp),intent(in) :: kk(3)
1564 : integer,intent(in) :: ngfft(18)
1565 : complex(dp),target,intent(inout) :: ur(nfft,ndat)
1566 :
1567 : !Local variables-------------------------------
1568 : real(dp),contiguous,pointer :: ur_ptr(:,:,:)
1569 : ! *************************************************************************
1570 :
1571 0 : call c_f_pointer(c_loc(ur), ur_ptr, shape=[2, nfft, ndat])
1572 0 : call times_eikr(kk, ngfft, nfft, ndat, ur_ptr)
1573 :
1574 0 : end subroutine ctimes_eikr
1575 : !!***
1576 :
1577 : !!****f* m_fft_mesh/phase
1578 : !! NAME
1579 : !! phase
1580 : !!
1581 : !! FUNCTION
1582 : !! Compute ph(ig)=$\exp(\pi\ i \ n/ngfft)$ for n=0,...,ngfft/2,-ngfft/2+1,...,-1
1583 : !! while ig runs from 1 to ngfft.
1584 : !!
1585 : !! INPUTS
1586 : !! ngfft=number of points
1587 : !!
1588 : !! OUTPUT
1589 : !! ph(2*ngfft)=phase array (complex)
1590 : !!
1591 : !! NOTES
1592 : !! XG 990504: changed the formulation, in order to preserve
1593 : !! the invariance between n and -n, that was broken for n=ngfft/2 if ngfft even.
1594 : !! Simply suppresses the corresponding sine.
1595 : !!
1596 : !! SOURCE
1597 :
1598 25218 : subroutine phase(ngfft, ph)
1599 :
1600 : !Arguments ------------------------------------
1601 : !scalars
1602 : integer,intent(in) :: ngfft
1603 : !arrays
1604 : real(dp),intent(out) :: ph(2*ngfft)
1605 :
1606 : !Local variables-------------------------------
1607 : integer :: id,ig,nn
1608 : real(dp) :: arg,fac
1609 : ! *************************************************************************
1610 :
1611 25218 : id=ngfft/2+2
1612 25218 : fac=pi/dble(ngfft)
1613 490086 : do ig=1,ngfft
1614 464868 : nn=ig-1-(ig/id)*ngfft
1615 464868 : arg=fac*dble(nn)
1616 464868 : ph(2*ig-1)=cos(arg)
1617 490086 : ph(2*ig) =sin(arg)
1618 : end do
1619 :
1620 : ! XG 990504 Here zero the corresponding sine
1621 25218 : if((ngfft/2)*2==ngfft) ph(2*(id-1))=zero
1622 :
1623 25218 : end subroutine phase
1624 : !!***
1625 :
1626 : !!****f* ABINIT/mkgrid_fft
1627 : !! NAME
1628 : !! mkgrid_fft
1629 : !!
1630 : !! FUNCTION
1631 : !! Sets the grid of fft (or real space) points to be treated.
1632 : !!
1633 : !! INPUTS
1634 : !!
1635 : !! OUTPUT
1636 : !!
1637 : !! SOURCE
1638 :
1639 0 : subroutine mkgrid_fft(ffti3_local,fftn3_distrib,gridcart,nfft,ngfft,rprimd)
1640 :
1641 : !Arguments ------------------------------------
1642 : integer, intent(in) :: nfft
1643 : integer,intent(in) :: ngfft(18)
1644 : integer, dimension(*), intent(in) :: ffti3_local,fftn3_distrib
1645 : real(dp), dimension(3,nfft), intent(out) :: gridcart
1646 : real(dp),intent(in) :: rprimd(3,3)
1647 :
1648 : !Local variables-------------------------------
1649 : integer :: ind,i1,i2,i3,i3loc,me,nproc
1650 : integer :: n1,n2,n3
1651 : real(dp), dimension(3) :: coord
1652 0 : real(dp), dimension(3,nfft) :: gridred
1653 : ! *************************************************************************
1654 :
1655 0 : n1 = ngfft(1)
1656 0 : n2 = ngfft(2)
1657 0 : n3 = ngfft(3)
1658 0 : nproc = ngfft(10)
1659 0 : me = ngfft(11)
1660 :
1661 0 : do i3 = 1, n3, 1
1662 0 : if(fftn3_distrib(i3) == me) then !MPI
1663 0 : i3loc=ffti3_local(i3)
1664 0 : coord(3) = real(i3 - 1, dp) / real(n3, dp)
1665 0 : do i2 = 1, n2, 1
1666 0 : coord(2) = real(i2 - 1, dp) / real(n2, dp)
1667 0 : do i1 = 1, n1, 1
1668 0 : ind=i1+(i2-1)*n1+(i3loc-1)*n1*n2
1669 0 : coord(1) = real(i1 - 1, dp) / real(n1, dp)
1670 0 : gridred(:, ind) = coord(:)
1671 : end do
1672 : end do
1673 : end if
1674 : end do
1675 0 : call xred2xcart(nfft, rprimd, gridcart, gridred)
1676 :
1677 0 : end subroutine mkgrid_fft
1678 : !!***
1679 :
1680 : !!****f* ABINIT/supercell_fft
1681 : !! NAME
1682 : !! supercell_fft
1683 : !!
1684 : !! FUNCTION
1685 : !! Build table for supercell calculations
1686 : !! Note that this version is only used in the BSE code as memory scales badly with the supercell size.
1687 : !!
1688 : !! INPUTS
1689 : !! ncells(3)= Number of cells along the three reduced directions
1690 : !!
1691 : !! OUTPUT
1692 : !! sc_nfft=The total number of points in the supercell.
1693 : !! sc2uc(sc_fft): The image of the point in the small box.
1694 : !! scred(3,sc_nfft): The reduced coordinates of the point in the supercell in terms of rprimd.
1695 : !!
1696 : !! SOURCE
1697 :
1698 0 : subroutine supercell_fft(ncells, ngfft, sc_nfft, sc_ngfft, sc2uc, scred)
1699 :
1700 : !Arguments ------------------------------------
1701 : integer,intent(in) :: ncells(3), ngfft(18)
1702 : integer,intent(out) :: sc_nfft, sc_ngfft(18)
1703 : integer,allocatable,intent(out) :: sc2uc(:)
1704 : real(dp),allocatable,intent(out) :: scred(:,:)
1705 :
1706 : !Local variables-------------------------------
1707 : integer :: irc, ir1, ir2, ir3, wp1, wp2, wp3, wp_idx
1708 : ! *************************************************************************
1709 :
1710 0 : sc_ngfft = ngfft
1711 0 : sc_ngfft(1:3) = ncells(1:3) * ngfft(1:3)
1712 0 : sc_ngfft(4:6) = sc_ngfft(1:3)
1713 : !sc_ngfft(4) = 2*(sc_ngfft(1)/2)+1
1714 : !sc_ngfft(5) = 2*(sc_ngfft(2)/2)+1
1715 : !sc_ngfft(6) = sc_ngfft(3)
1716 0 : sc_nfft = product(sc_ngfft(1:3)) ! Total number of points in the supercell
1717 :
1718 0 : ABI_MALLOC(sc2uc, (sc_nfft))
1719 0 : ABI_MALLOC(scred, (3, sc_nfft))
1720 :
1721 0 : irc = 0
1722 0 : do ir3=0,sc_ngfft(3)-1 ! Loop over the points in the supercell.
1723 0 : do ir2=0,sc_ngfft(2)-1
1724 0 : do ir1=0,sc_ngfft(1)-1
1725 0 : irc = 1+irc
1726 0 : wp1=MODULO(ir1, ngfft(1)) ! The FFT index of the point wrapped into the unit cell.
1727 0 : wp2=MODULO(ir2, ngfft(2))
1728 0 : wp3=MODULO(ir3, ngfft(3))
1729 0 : wp_idx = 1 + wp1 + wp2*ngfft(1) + wp3*ngfft(1)*ngfft(2)
1730 0 : sc2uc(irc) = wp_idx
1731 0 : scred(1,irc) = DBLE(ir1)/ngfft(1) ! Reduced coordinates in terms of the unit cell lattice vectors
1732 0 : scred(2,irc) = DBLE(ir2)/ngfft(2)
1733 0 : scred(3,irc) = DBLE(ir3)/ngfft(3)
1734 : end do
1735 : end do
1736 : end do
1737 :
1738 0 : end subroutine supercell_fft
1739 :
1740 11472048 : END MODULE m_fft_mesh
1741 : !!***
|
