Branch data 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 [ + - ]: 5739117 : 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 : 5739117 : g3_min = gbound(3, 2)
145 : 5739117 : g3_max = gbound(4, 2)
146 : :
147 : 5739117 : zpad%n_zplanes = g3_max - g3_min + 1
148 : :
149 [ - + + - : 17217351 : ABI_MALLOC(zpad%zplane, (2, nz))
+ - + - -
+ - + -
+ ]
150 [ - + + - : 17217351 : ABI_MALLOC(zpad%linex2ifft_yz, (2, nx*ny*nz))
+ - + - -
+ - + -
+ ]
151 : :
152 : : ! Loop over the z-planes intersecting the G-sphere.
153 : 5739117 : nlinex = 0
154 [ + + ]: 35521404 : do gg3=1,zpad%n_zplanes
155 : : !
156 [ + + ]: 29782287 : if (gg3<=g3_max+1) then
157 : : ifft_g3 = gg3
158 : : else
159 : 12325421 : 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 : 29782287 : igb=2*gg3+3
164 : 29782287 : g2min = gbound(igb ,2)
165 : 29782287 : g2max = gbound(igb+1,2)
166 : :
167 : 29782287 : zpad%zplane(1,gg3) = ifft_g3
168 : 29782287 : 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 [ + + ]: 99339905 : do jj=1,g2max+1
174 : 69557618 : nlinex = nlinex + 1
175 : 69557618 : zpad%linex2ifft_yz(1,nlinex) = jj
176 : 99339905 : zpad%linex2ifft_yz(2,nlinex) = ifft_g3
177 : : end do
178 : :
179 [ + + ]: 78086664 : do jj=g2min+ny+1,ny
180 : 42565260 : nlinex = nlinex + 1
181 : 42565260 : zpad%linex2ifft_yz(1,nlinex) = jj
182 : 72347547 : zpad%linex2ifft_yz(2,nlinex) = ifft_g3
183 : : end do
184 : : end do
185 : :
186 : 5739117 : 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 : 5739117 : subroutine zpad_free(zpad)
206 : :
207 : : !Arguments ------------------------------------
208 : : class(zpad_t),intent(inout) :: zpad
209 : : ! *************************************************************************
210 : :
211 [ + - ]: 5739117 : ABI_SFREE(zpad%zplane)
212 [ + - ]: 5739117 : ABI_SFREE(zpad%linex2ifft_yz)
213 : :
214 : 5739117 : 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 : 591 : 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 : 591 : integer,pointer :: symrel(:,:,:)
287 : 591 : real(dp),pointer :: tnons(:,:)
288 : : !************************************************************************
289 : :
290 : : DBG_ENTER("COLL")
291 : :
292 [ + + + - : 2364 : 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 : 591 : "to get rid of spurious tnons"
300 : :
301 [ + + ]: 591 : ount = std_out; if (present(unit)) ount = unit
302 : :
303 : 591 : nsym = Cryst%nsym
304 : 591 : symrel => Cryst%symrel
305 : 591 : tnons => Cryst%tnons
306 : :
307 : : ! Calculate the limits of the sphere of npwwfn G-vectors in each direction.
308 [ + + ]: 333936 : m1 = MAXVAL(ABS(gvec(1,1:npwwfn)))
309 [ + + ]: 333936 : m2 = MAXVAL(ABS(gvec(2,1:npwwfn)))
310 [ + + ]: 333936 : 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 [ + + ]: 117984 : mm1 = MAXVAL(ABS(gvec(1,1:npwsigx)))
318 [ + + ]: 117984 : mm2 = MAXVAL(ABS(gvec(2,1:npwsigx)))
319 [ + + ]: 117984 : mm3 = MAXVAL(ABS(gvec(3,1:npwsigx)))
320 : :
321 : 591 : mm1=mm1+mG0(1)
322 : 591 : mm2=mm2+mG0(2)
323 : 591 : 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 : 591 : write(msg,'(2(2a,i8,a,3i6),2a,3i3)')ch10,&
333 : 591 : ' setmesh: npwwfn = ',npwwfn, '; Max (m1,m2,m3) = ',m1,m2,m3,ch10,&
334 : 591 : ' npweps/npwsigx= ',npwsigx,'; Max (mm1,mm2,mm3)= ',mm1,mm2,mm3,ch10,&
335 : 1182 : ' mG0 added = ',mG0(:)
336 : 591 : call wrtout(ount, msg)
337 : : !
338 : : ! === Different FFT grids according to method ==
339 : 591 : 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 : 553 : ecutwfn=-one ! Calculate the radius of the sphere of npwwfn G-vectors.
370 [ + + ]: 317612 : do ig=1,npwwfn
371 : 317059 : g1=REAL(gvec(1,ig))
372 : 317059 : g2=REAL(gvec(2,ig))
373 : 317059 : g3=REAL(gvec(3,ig))
374 : : gsq= gmet(1,1)*g1**2+gmet(2,2)*g2**2+gmet(3,3)*g3**2+ &
375 : 317059 : two*(gmet(1,2)*g1*g2+gmet(1,3)*g1*g3+gmet(2,3)*g2*g3)
376 : 317612 : ecutwfn=MAX(ecutwfn,gsq)
377 : : end do
378 : 553 : rwfn=SQRT(ecutwfn); ecutwfn=two*ecutwfn*pi**2
379 : :
380 : : ! * Calculate the radius of the sphere of (npwsigx|npweps) G-vectors.
381 : 553 : ecutsigx=-one
382 [ + + ]: 112120 : do ig=1,npwsigx
383 : 111567 : g1=REAL(gvec(1,ig))
384 : 111567 : g2=REAL(gvec(2,ig))
385 : 111567 : g3=REAL(gvec(3,ig))
386 : : gsq= gmet(1,1)*g1**2+gmet(2,2)*g2**2+gmet(3,3)*g3**2+ &
387 : 111567 : two*(gmet(1,2)*g1*g2+gmet(1,3)*g1*g3+gmet(2,3)*g2*g3)
388 : 112120 : ecutsigx=MAX(ecutsigx,gsq)
389 : : end do
390 : 553 : rsigx=SQRT(ecutsigx); ecutsigx=two*ecutsigx*pi**2
391 : :
392 : : write(msg,'(a,f7.3,3a,f7.3,a)')&
393 : 553 : ' calculated ecutwfn = ',ecutwfn, ' [Ha] ',ch10,&
394 : 1106 : ' calculated ecutsigx/ecuteps = ',ecutsigx,' [Ha]'
395 : 553 : 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 : 553 : reff=rsigx+rwfn
407 [ + + ]: 553 : if (method==3) reff=two*rwfn ! Yields back the GS FFT grid if full wavefunctions are considered.
408 : 553 : ecuteff=two*(pi*reff)**2
409 : 553 : gsqmax=reff**2
410 : :
411 : 553 : write(msg,'(a,i2,a,f7.3,a)')' using method = ',method,' with ecuteff = ',ecuteff,' [Ha]'
412 : 553 : 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 [ + + ]: 553 : if (method==2) then
419 : 549 : ig1max=mm1+m1+1
420 : 549 : ig2max=mm2+m2+1
421 : 549 : 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 : 553 : m1=-1; m2=-1; m3=-1
431 [ + + ]: 5715 : do ig1=0,ig1max
432 [ + + ]: 61701 : do ig2=0,ig2max
433 [ + + ]: 738656 : do ig3=0,ig3max
434 : 677508 : g1=REAL(ig1)
435 : 677508 : g2=REAL(ig2)
436 : 677508 : g3=REAL(ig3)
437 : : gsq= gmet(1,1)*g1**2+gmet(2,2)*g2**2+gmet(3,3)*g3**2+ &
438 : 677508 : two*(gmet(1,2)*g1*g2+gmet(1,3)*g1*g3+gmet(2,3)*g2*g3)
439 [ + + ]: 677508 : if (gsq>gsqmax+tol6) CYCLE ! tol6 to improve portability
440 : 344811 : m1=MAX(m1,ig1)
441 : 344811 : m2=MAX(m2,ig2)
442 : 733494 : m3=MAX(m3,ig3)
443 : : end do
444 : : end do
445 : : end do
446 : :
447 : : case default
448 [ - + + - ]: 591 : 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 [ + + + + : 591 : 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 : 591 : m1=MAX(m1,mm1)
463 : 591 : m2=MAX(m2,mm2)
464 : 591 : m3=MAX(m3,mm3)
465 : :
466 [ - + ]: 591 : 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 : 591 : 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 [ + + ]: 2364 : fftnons(:)=1
507 [ + + ]: 2364 : do ii=1,3
508 : 1773 : fftnons(ii)=1
509 [ + + ]: 74880 : do is=1,nsym
510 : 72516 : nt=denominator(tnons(ii,is), ierr)
511 [ + + ]: 74289 : if (((fftnons(ii)/nt)*nt)/=fftnons(ii)) fftnons(ii)=mincm(fftnons(ii),nt)
512 : : end do
513 : : end do
514 : 591 : write(msg,'(a,3(i0,1x))')' setmesh: divisor mesh ',fftnons(:)
515 : 591 : call wrtout(ount, msg)
516 : : !
517 : : ! 2) Check if also rotations preserve the grid.
518 : : ! * Use previous m values as Initial guess.
519 : 591 : call size_goed_fft(m1,fftsym(1),ierr)
520 [ - + ]: 591 : ABI_CHECK(ierr == 0, sjoin("size_goed_fft failed", ch10, tnons_warn))
521 : 591 : call size_goed_fft(m2,fftsym(2),ierr)
522 [ - + ]: 591 : ABI_CHECK(ierr == 0, sjoin("size_goed_fft failed", ch10, tnons_warn))
523 : 591 : call size_goed_fft(m3,fftsym(3),ierr)
524 [ - + ]: 591 : ABI_CHECK(ierr == 0, sjoin("size_goed_fft failed", ch10, tnons_warn))
525 : :
526 : 591 : mdum(1)=m1
527 : 591 : mdum(2)=m2
528 : 591 : mdum(3)=m3
529 : :
530 : 591 : 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 [ + + + + : 4109 : (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 : 591 : n1=fftsym(1)
546 : 591 : n2=fftsym(2)
547 : 591 : n3=fftsym(3); nfftot=n1*n2*n3
548 : :
549 : : if (.not.( check_rot_fft(nsym,symrel,n1,n2,n3)) &
550 [ - + + - : 591 : .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 : 591 : ' setmesh: FFT mesh size selected = ',n1,'x',n2,'x',n3,ch10,&
560 : 1182 : ' total number of points = ',nfftot,ch10
561 : 591 : call wrtout(ount, msg)
562 [ + + ]: 591 : if (ount /= dev_null) call wrtout(ab_out, msg)
563 : :
564 : 591 : ngfft(1)=n1
565 : 591 : ngfft(2)=n2
566 : 591 : ngfft(3)=n3
567 : 591 : ngfft(4)=2*(ngfft(1)/2)+1
568 : 591 : ngfft(5)=2*(ngfft(2)/2)+1
569 : 591 : 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 : 591 : fftalg=ngfft(7); fftalga=fftalg/100; fftalgc=MOD(fftalg,10)
575 : :
576 [ + - - + ]: 591 : 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 : 591 : 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 : 4093 : 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 : 4093 : check_rot_fft=.TRUE.
639 [ + + ]: 139929 : 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 [ + + + + : 136784 : MOD(symrel(1,3,is)*nr1, nr3) /=0 .or. &
+ + + - +
- + - ]
645 : : MOD(symrel(2,3,is)*nr2, nr3) /=0 &
646 : 3145 : ) then
647 : : check_rot_fft=.FALSE.; EXIT
648 : : end if
649 : : end do
650 : :
651 : 4093 : 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 : 2461 : 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 : 4922 : integer :: Rm1(3,3,nsym),r1_FFT(3),red2fft(3,3)
790 : 4922 : 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 : 2461 : ngfft1 = ngfft(1); ngfft2 = ngfft(2); ngfft3 = ngfft(3)
795 : :
796 [ + + ]: 24610 : red2fft = reshape([ngfft1, 0, 0, 0, ngfft2, 0, 0, 0, ngfft3], [3, 3])
797 [ + + ]: 24610 : 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 [ + + ]: 105837 : do isym=1,nsym
801 : 103376 : call mati3inv(symrel(:,:,isym), Rm1(:,:,isym))
802 [ + + + + : 2584400 : Rm1(:,:,isym) = transpose(Rm1(:,:,isym))
+ + + + ]
803 [ + + + + : 10647728 : Rm1_FFT(:,:,isym) = matmul(matmul(red2fft, Rm1(:,:,isym)), fft2red)
+ + + + +
+ + + + +
+ + + + +
+ ]
804 [ + + + + : 2896989 : tnons_FFT(:,isym) = matmul(red2fft, tnons(:,isym))
+ + + + +
+ ]
805 : : end do
806 : :
807 [ + + + + ]: 415965 : err(:,:) = zero
808 : :
809 : : !$OMP PARALLEL DO PRIVATE(R1_FFT,ir1,R2_FFT,jx,jy,jz) reduction(MAX:err)
810 [ + + ]: 44577 : do iz=0,ngfft3-1
811 : 42116 : R1_FFT(3) = dble(iz)
812 [ + + ]: 875382 : do iy=0,ngfft2-1
813 : 830805 : R1_FFT(2) = dble(iy)
814 [ + + ]: 20225126 : do ix=0,ngfft1-1
815 : 19352205 : R1_FFT(1) = dble(ix)
816 : 19352205 : ir1 = 1+ix+iy*ngfft1+iz*ngfft1*ngfft2
817 [ + + ]: 692132326 : do isym=1,nsym
818 : : ! Form R^-1 (r-\tau) in the FFT basis.
819 [ + + + + : 10751189056 : R2_FFT(:) = MATMUL(Rm1_FFT(:,:,isym),R1_FFT(:)-tnons_FFT(:,isym))
+ + ]
820 : 671949316 : jx = NINT(R2_FFT(1)); err(1,isym) = MAX(err(1,isym), ABS(R2_FFT(1)-jx)/ngfft1)
821 : 671949316 : jy = NINT(R2_FFT(2)); err(2,isym) = MAX(err(2,isym), ABS(R2_FFT(2)-jy)/ngfft2)
822 : 671949316 : jz = NINT(R2_FFT(3)); err(3,isym) = MAX(err(3,isym), ABS(R2_FFT(3)-jz)/ngfft3)
823 : 671949316 : jx = MODULO(jx, ngfft1)
824 : 671949316 : jy = MODULO(jy, ngfft2)
825 : 671949316 : jz = MODULO(jz, ngfft3)
826 : 691301521 : irottb(ir1,isym) = 1+jx+jy*ngfft1+jz*ngfft1*ngfft2
827 : : end do
828 : : end do
829 : : end do
830 : : end do
831 : :
832 : 2461 : preserve = .TRUE.
833 [ + + ]: 105837 : do isym=1,nsym
834 [ + + + + : 407309 : 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 : 2461 : 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 : 117943346 : 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 [ + - + - ]: 117943346 : if (ig <= 0 .or. ig > ng) then
1022 : : ! Wrong ig, returns huge. Parent code will likely crash with SIGSEV.
1023 : 117943346 : gc = huge(1)
1024 : : return
1025 : : end if
1026 : :
1027 [ + + ]: 117943346 : if ( ig > ng/2 + 1) then
1028 : 53332941 : gc = ig - ng -1
1029 : : else
1030 : 64610405 : 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 : 22251 : 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 : 22251 : ifft=0; gsq_max=smallest_real
1118 [ + + ]: 333070 : do i3=1,ngfft(3)
1119 : 310819 : g3 = ig2gfft(i3,ngfft(3))
1120 [ + + ]: 5703703 : do i2=1,ngfft(2)
1121 : 5370633 : g2 = ig2gfft(i2,ngfft(2))
1122 [ + + ]: 117941003 : do i1=1,ngfft(1)
1123 : 112259551 : g1 = ig2gfft(i1,ngfft(1))
1124 : 112259551 : ifft = ifft+1
1125 [ + + ]: 449038204 : 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 : 112259551 : +gmet(3,1)*(kpt(3)+dble(i3))*(kpt(1)+dble(i1)))
1132 : 117630184 : gsq_max = MAX(dsq,gsq_max)
1133 : : end do
1134 : : end do
1135 : : end do
1136 : :
1137 : 22251 : 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 [ - + - - : 11478234 : END MODULE m_fft_mesh
+ + - + +
- - + -
+ ]
1741 : : !!***
|
