Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2020 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_types_fftGrid
8 : use m_constants
9 : USE m_juDFT
10 : TYPE t_fftGrid
11 :
12 : INTEGER :: extent(3) = [-1, -1, -1]
13 : INTEGER :: dimensions(3) = [-1, -1, -1]
14 : INTEGER :: gridLength = -1
15 : COMPLEX, ALLOCATABLE :: grid(:)
16 :
17 : CONTAINS
18 :
19 : PROCEDURE :: t_fftGrid_init,t_fftGrid_init_dims
20 : GENERIC :: init => t_fftGrid_init,t_fftGrid_init_dims
21 : FINAL :: free
22 : PROCEDURE :: putFieldOnGrid
23 : PROCEDURE :: takeFieldFromGrid
24 : PROCEDURE :: getRealPartOfGrid
25 : PROCEDURE :: putStateOnGrid
26 : procedure :: put_state_on_external_grid
27 : PROCEDURE :: putRealStateOnGrid
28 : PROCEDURE :: putComplexStateOnGrid
29 : PROCEDURE :: fillStateIndexArray
30 : PROCEDURE :: fillFieldSphereIndexArray
31 : PROCEDURE :: getElement
32 : procedure :: g2fft => map_g_to_fft_grid
33 : PROCEDURE :: perform_fft
34 : END TYPE t_fftGrid
35 :
36 : PUBLIC :: t_fftGrid, put_real_on_external_grid, put_cmplx_on_external_grid
37 :
38 : CONTAINS
39 75027 : subroutine perform_fft(grid,forward)
40 : use m_types_fft
41 : implicit none
42 : CLASS(t_fftGrid), INTENT(INOUT) :: grid
43 : LOGICAL,INTENT(IN) :: forward
44 :
45 300108 : type(t_fft) :: fft
46 :
47 300108 : if (size(grid%grid) .ne. product(grid%dimensions)) call juDFT_error('array bounds are inconsistent', calledby='perform_fft')
48 :
49 75027 : call fft%init(grid%dimensions, forward)
50 75027 : call fft%exec(grid%grid)
51 75027 : call fft%free()
52 75027 : end subroutine perform_fft
53 :
54 10080 : function map_g_to_fft_grid(grid, g_in) result(g_idx)
55 : implicit none
56 : CLASS(t_fftGrid), INTENT(IN) :: grid
57 : integer, intent(in) :: g_in(3)
58 : integer :: g_idx
59 :
60 : integer :: shifted_g(3)
61 :
62 : ! the fft-grid starts at g=0, not -g_max
63 : ! therefore all negative g's need to be shifted
64 :
65 40320 : shifted_g = merge(g_in + grid%dimensions, g_in, g_in < 0)
66 :
67 : ! map it to 1d
68 : g_idx = shifted_g(1) &
69 : + shifted_g(2) * grid%dimensions(1) &
70 10080 : + shifted_g(3) * grid%dimensions(1) * grid%dimensions(2)
71 10080 : end function map_g_to_fft_grid
72 :
73 18191 : SUBROUTINE t_fftGrid_init(this, cell, sym, gCutoff, gzCutoff)
74 : USE m_constants
75 : USE m_boxdim
76 : USE m_spgrot
77 : USE m_ifft
78 : USE m_types_cell
79 : USE m_types_sym
80 : IMPLICIT NONE
81 : CLASS(t_fftGrid), INTENT(INOUT) :: this
82 : TYPE(t_cell), INTENT(IN) :: cell
83 : TYPE(t_sym), INTENT(IN) :: sym
84 : REAL, INTENT(IN) :: gCutoff
85 : REAL, OPTIONAL, INTENT(IN) :: gzCutoff
86 :
87 :
88 : INTEGER, ALLOCATABLE :: ig(:, :, :)
89 :
90 : INTEGER :: k1, k2, k3, i
91 : INTEGER :: mxx(3), kVec(3), kRot(3, sym%nop), inv_du(sym%nop), tempDim(3)
92 : REAL :: gCutoffSquared, gSquared
93 : REAL :: arltv(3), g(3)
94 :
95 18191 : if (.not.present(gzCutoff)) then
96 72764 : this%extent = calc_extent(cell, sym, gCutoff)
97 : else
98 0 : this%extent = calc_extent(cell, sym, gCutoff, gzCutoff)
99 : end if
100 72764 : this%dimensions = 2*this%extent + 1
101 :
102 72764 : do i = 1,3
103 72764 : this%dimensions(i) = next_optimal_fft_size(this%dimensions(i))
104 : enddo
105 72764 : this%gridLength = product(this%dimensions)
106 :
107 18191 : IF (ALLOCATED(this%grid)) THEN
108 412 : IF (size(this%grid)==this%gridlength) RETURN
109 1 : DEALLOCATE (this%grid)
110 : ENDIF
111 67714447 : ALLOCATE (this%grid(0:this%gridLength - 1), source=CMPLX_NOT_INITALIZED)
112 : END SUBROUTINE t_fftGrid_init
113 :
114 63787 : SUBROUTINE t_fftGrid_init_dims(this,dims)
115 : IMPLICIT NONE
116 : CLASS(t_fftGrid), INTENT(INOUT) :: this
117 : INTEGER,INTENT(IN) :: dims(3)
118 255148 : this%dimensions=dims
119 255148 : this%gridLength = product(this%dimensions)
120 :
121 63787 : IF (ALLOCATED(this%grid)) DEALLOCATE (this%grid)
122 480749340 : ALLOCATE (this%grid(0:this%gridLength - 1), source=CMPLX_NOT_INITALIZED)
123 63787 : END SUBROUTINE t_fftGrid_init_dims
124 :
125 :
126 :
127 52275 : SUBROUTINE putFieldOnGrid(this, stars,field, cell, gCutoff, gzCutoff, firstderiv,secondderiv,l_2D)
128 : USE m_types_stars
129 : USE m_types_cell
130 : IMPLICIT NONE
131 : CLASS(t_fftGrid), INTENT(INOUT) :: this
132 : TYPE(t_stars), INTENT(IN) :: stars
133 : COMPLEX, INTENT(IN) :: field(:) ! length is stars%ng3
134 : TYPE(t_cell),INTENT(IN),OPTIONAL:: cell
135 : REAL, OPTIONAL, INTENT(IN) :: gCutoff, gzCutoff
136 : real,optional,intent(in) :: firstderiv(3),secondderiv(3)
137 : LOGICAL,OPTIONAL,intent(in) :: l_2d
138 :
139 : INTEGER :: x, y, z, iStar
140 : INTEGER :: xGrid, yGrid, zGrid, layerDim
141 : REAL :: gCutoffInternal, gvec(3)
142 : COMPLEX :: fct
143 : LOGICAL :: twoD, l_insph
144 52275 : gCutoffInternal = 1.0e99
145 52275 : IF (PRESENT(gCutoff)) gCutoffInternal = gCutoff
146 :
147 52275 : layerDim = this%dimensions(1)*this%dimensions(2)
148 :
149 425791420 : this%grid(:) = CMPLX(0.0, 0.0)
150 :
151 52275 : twod=.false.
152 52275 : if (present(l_2d)) twoD=l_2d
153 42000 : if (twoD.and.this%dimensions(3)>1) call juDFT_error("Bug in putFieldOnGrid: no two-D grid")
154 :
155 348538 : DO z = merge(0,-stars%mx3,twoD), merge(0,stars%mx3,twoD)
156 296263 : zGrid = MODULO(z, merge(1,this%dimensions(3),twoD)) !always 0 in 2d-case
157 7603107 : DO y = -stars%mx2, stars%mx2
158 7254569 : yGrid = MODULO(y, this%dimensions(2))
159 210485971 : DO x = -stars%mx1, stars%mx1
160 202935139 : IF (twod) THEN
161 48778000 : iStar=stars%i2g(x,y)
162 48778000 : fct=stars%r2gphs(x,y)
163 : ELSE
164 154157139 : iStar = stars%ig(x, y, z)
165 154157139 : fct=stars%rgphs(x, y, z)
166 : endif
167 202935139 : if (present(firstderiv)) THEN
168 1560318000 : fct=fct*cmplx(0.0,-1*dot_product(firstderiv,matmul(real([x,y,z]),cell%bmat)))
169 926952000 : if (present(secondderiv)) fct=fct*cmplx(0.0,-1*dot_product(secondderiv,matmul(real([x,y,z]),cell%bmat)))
170 : endif
171 202935139 : IF (iStar .EQ. 0) CYCLE
172 90902601 : l_insph = stars%sk3(iStar) .LE. gCutoffInternal
173 90902601 : IF (PRESENT(gzCutoff)) gvec = matmul(real([x,y,z]),cell%bmat)
174 0 : IF (PRESENT(gzCutoff)) l_insph = (SQRT(DOT_PRODUCT(gvec(:2),gvec(:2))).LE.gCutoffInternal).AND.(ABS(gvec(3)).LE.gzCutoff)
175 90902601 : IF (.NOT.l_insph) CYCLE
176 83919043 : IF (twod) THEN
177 34472000 : if(stars%sk2(istar)>gCutoffInternal) CYCLE
178 : ENDIF
179 83919043 : xGrid = MODULO(x, this%dimensions(1))
180 210189708 : this%grid(xGrid + this%dimensions(1)*yGrid + layerDim*zGrid) = field(iStar)*fct
181 : END DO
182 : END DO
183 : END DO
184 :
185 52275 : END SUBROUTINE putFieldOnGrid
186 :
187 45562 : SUBROUTINE takeFieldFromGrid(this, stars, field, gCutoff, gzCutoff, l_2d)
188 : USE m_types_stars
189 : IMPLICIT NONE
190 : CLASS(t_fftGrid), INTENT(IN) :: this
191 : TYPE(t_stars), INTENT(IN) :: stars
192 : COMPLEX, INTENT(INOUT) :: field(:)
193 : REAL, OPTIONAL, INTENT(IN) :: gCutoff, gzCutoff
194 : LOGICAL,OPTIONAL,intent(in) :: l_2d
195 :
196 :
197 : INTEGER :: x, y, z, iStar
198 : INTEGER :: xGrid, yGrid, zGrid, layerDim
199 : REAL :: elementWeight, gCutoffInternal
200 : COMPLEX :: fct
201 : LOGICAL :: twod, l_insph
202 :
203 45562 : gCutoffInternal = 1.0e99
204 45562 : IF (PRESENT(gCutoff)) gCutoffInternal = gCutoff
205 :
206 45562 : twod=.false.
207 45562 : if (present(l_2d)) twoD=l_2d
208 26400 : if (twoD.and.this%dimensions(3)>1) call juDFT_error("Bug in takeFieldFromGrid: no two-D grid")
209 :
210 50640062 : field(:) = CMPLX(0.0, 0.0)
211 45562 : layerDim = this%dimensions(1)*this%dimensions(2)
212 :
213 530784 : DO z = merge(0,-stars%mx3,twoD), merge(0,stars%mx3,twoD)
214 485222 : zGrid = MODULO(z, merge(1,this%dimensions(3),twoD)) !always 0 in 2d-case
215 11558834 : DO y = -stars%mx2, stars%mx2
216 11028050 : yGrid = MODULO(y, this%dimensions(2))
217 285623738 : DO x = -stars%mx1, stars%mx1
218 274110466 : IF(twod) THEN
219 23948000 : iStar = stars%i2g(x, y)
220 23948000 : fct = conjg(stars%r2gphs(x,y))
221 : ELSE
222 250162466 : iStar = stars%ig(x, y, z)
223 250162466 : fct=CONJG(stars%rgphs(x, y, z))
224 : ENDIF
225 274110466 : IF (iStar .EQ. 0) CYCLE
226 108861956 : l_insph = stars%sk3(iStar) .LE. gCutoffInternal
227 108861956 : IF (PRESENT(gzCutoff)) l_insph = (stars%sk2(stars%i2g(x, y)).LE.gCutoffInternal).AND.(stars%ab3(iStar).LE.gzCutoff)
228 108861956 : IF (.NOT.l_insph) CYCLE
229 80429794 : xGrid = MODULO(x, this%dimensions(1))
230 285138516 : field(iStar) = field(iStar) + this%grid(xGrid + this%dimensions(1)*yGrid + layerDim*zGrid) * fct
231 : END DO
232 : END DO
233 : END DO
234 :
235 45562 : if (twoD) THEN
236 26400 : elementWeight = 1.0/(this%dimensions(1)*this%dimensions(2))
237 8464000 : field(:) = elementWeight*field(:)/stars%nstr2(:)
238 : ELSE
239 19162 : elementWeight = 1.0/(this%dimensions(1)*this%dimensions(2)*this%dimensions(3))
240 42176062 : field(:) = elementWeight*field(:)/stars%nstr(:)
241 : endif
242 :
243 45562 : END SUBROUTINE takeFieldFromGrid
244 :
245 87992 : SUBROUTINE putStateOnGrid(this, lapw, iSpin, zMat, iState)
246 : USE m_types_lapw
247 : USE m_types_mat
248 : IMPLICIT NONE
249 : CLASS(t_fftGrid), INTENT(INOUT) :: this
250 : TYPE(t_lapw), INTENT(IN) :: lapw
251 : TYPE(t_mat), INTENT(IN) :: zMat
252 : INTEGER, INTENT(IN) :: iSpin
253 : INTEGER, INTENT(IN) :: iState
254 :
255 87992 : IF (zMat%l_real) THEN
256 40172 : CALL putRealStateOnGrid(this, lapw, iSpin, zMat%data_r(:, iState))
257 : ELSE
258 47820 : CALL putComplexStateOnGrid(this, lapw, iSpin, zMat%data_c(:, iState))
259 : END IF
260 87992 : END SUBROUTINE putStateOnGrid
261 :
262 5350 : SUBROUTINE put_state_on_external_grid(this, lapw, iSpin, zMat, iState, ext_grid, l_gpu)
263 : USE m_types_lapw
264 : USE m_types_mat
265 : IMPLICIT NONE
266 : CLASS(t_fftGrid), INTENT(INOUT) :: this
267 : TYPE(t_lapw), INTENT(IN) :: lapw
268 : TYPE(t_mat), INTENT(IN) :: zMat
269 : INTEGER, INTENT(IN) :: iSpin
270 : INTEGER, INTENT(IN) :: iState
271 : complex, intent(inout) :: ext_grid(0:)
272 : logical, intent(in), optional :: l_gpu
273 :
274 5350 : if (zMat%l_real) then
275 3964 : call put_real_on_external_grid(this, lapw, ispin, zMat%data_r(:, iState), ext_grid, l_gpu)
276 : else
277 1386 : call put_cmplx_on_external_grid(this, lapw, ispin, zMat%data_c(:, iState), ext_grid, l_gpu)
278 : endif
279 5350 : end subroutine put_state_on_external_grid
280 :
281 40172 : SUBROUTINE putRealStateOnGrid(this, lapw, iSpin, state)
282 : USE m_types_lapw
283 : IMPLICIT NONE
284 : CLASS(t_fftGrid), INTENT(INOUT) :: this
285 : TYPE(t_lapw), INTENT(IN) :: lapw
286 : REAL, INTENT(IN) :: state(:)
287 : INTEGER, INTENT(IN) :: iSpin
288 :
289 40172 : call put_real_on_external_grid(this, lapw, ispin, state, this%grid)
290 40172 : END SUBROUTINE putRealStateOnGrid
291 :
292 44136 : subroutine put_real_on_external_grid(this, lapw, ispin, state, ext_grid, l_gpu)
293 : USE m_types_lapw
294 : IMPLICIT NONE
295 : CLASS(t_fftGrid), INTENT(INOUT) :: this
296 : TYPE(t_lapw), INTENT(IN) :: lapw
297 : REAL, INTENT(IN) :: state(:)
298 : complex, intent(inout) :: ext_grid(0:)
299 : INTEGER, INTENT(IN) :: iSpin
300 : logical, intent(in), optional :: l_gpu
301 :
302 : logical :: use_gpu
303 : INTEGER :: xGrid, yGrid, zGrid, layerDim, iLAPW
304 :
305 44136 : layerDim = this%dimensions(1)*this%dimensions(2)
306 :
307 44136 : if(present(l_gpu)) then
308 3964 : use_gpu = l_gpu
309 : else
310 : use_gpu = .False.
311 : endif
312 :
313 44136 : if(use_gpu) then
314 : !$acc kernels
315 33247974 : ext_grid = cmplx_0
316 : !$acc end kernels
317 :
318 : !$acc parallel loop default(none) private(xGrid, yGrid, zGrid) &
319 : !$acc present(lapw, lapw%nv, lapw%gvec, this, this%dimensions, ext_grid, state) &
320 : !$acc copyin(layerDim, iSpin)
321 491986 : DO iLAPW = 1, lapw%nv(iSpin)
322 488022 : xGrid = MODULO(lapw%gvec(1, iLAPW, iSpin), this%dimensions(1))
323 488022 : yGrid = MODULO(lapw%gvec(2, iLAPW, iSpin), this%dimensions(2))
324 488022 : zGrid = MODULO(lapw%gvec(3, iLAPW, iSpin), this%dimensions(3))
325 491986 : ext_grid(xGrid + this%dimensions(1)*yGrid + layerDim*zGrid) = state(iLAPW)
326 : END DO
327 : !$acc end parallel loop
328 : else
329 136740670 : ext_grid = cmplx_0
330 :
331 4878286 : DO iLAPW = 1, lapw%nv(iSpin)
332 4838114 : xGrid = MODULO(lapw%gvec(1, iLAPW, iSpin), this%dimensions(1))
333 4838114 : yGrid = MODULO(lapw%gvec(2, iLAPW, iSpin), this%dimensions(2))
334 4838114 : zGrid = MODULO(lapw%gvec(3, iLAPW, iSpin), this%dimensions(3))
335 4878286 : ext_grid(xGrid + this%dimensions(1)*yGrid + layerDim*zGrid) = state(iLAPW)
336 : END DO
337 : endif
338 44136 : end subroutine put_real_on_external_grid
339 :
340 66614 : SUBROUTINE putComplexStateOnGrid(this, lapw, iSpin, state)
341 : USE m_types_lapw
342 : IMPLICIT NONE
343 : CLASS(t_fftGrid), INTENT(INOUT) :: this
344 : TYPE(t_lapw), INTENT(IN) :: lapw
345 : COMPLEX, INTENT(IN) :: state(:)
346 : INTEGER, INTENT(IN) :: iSpin
347 :
348 66614 : call put_cmplx_on_external_grid(this, lapw, ispin, state, this%grid)
349 66614 : END SUBROUTINE putComplexStateOnGrid
350 :
351 68000 : SUBROUTINE put_cmplx_on_external_grid(this, lapw, iSpin, state, ext_grid, l_gpu)
352 : USE m_types_lapw
353 : use m_judft
354 : IMPLICIT NONE
355 : CLASS(t_fftGrid), INTENT(INOUT) :: this
356 : TYPE(t_lapw), INTENT(IN) :: lapw
357 : COMPLEX, INTENT(IN) :: state(:)
358 : complex, intent(inout) :: ext_grid(0:)
359 : INTEGER, INTENT(IN) :: iSpin
360 : logical, intent(in), optional :: l_gpu
361 :
362 : logical :: use_gpu
363 : INTEGER :: xGrid, yGrid, zGrid, iLAPW
364 :
365 68000 : if(present(l_gpu)) then
366 1386 : use_gpu = l_gpu
367 : else
368 : use_gpu = .False.
369 : endif
370 :
371 68000 : if(use_gpu) then
372 : !$acc kernels
373 19161450 : ext_grid = cmplx_0
374 : !$acc end kernels
375 :
376 : !$acc parallel loop default(none) private(xGrid, yGrid, zGrid) &
377 : !$acc present(lapw, lapw%nv, lapw%gvec, this, this%dimensions, ext_grid, state, iSpin)
378 237636 : DO iLAPW = 1, lapw%nv(iSpin)
379 236250 : xGrid = MODULO(lapw%gvec(1, iLAPW, iSpin), this%dimensions(1))
380 236250 : yGrid = MODULO(lapw%gvec(2, iLAPW, iSpin), this%dimensions(2))
381 236250 : zGrid = MODULO(lapw%gvec(3, iLAPW, iSpin), this%dimensions(3))
382 237636 : ext_grid(xGrid + this%dimensions(1)*yGrid + (this%dimensions(1)*this%dimensions(2))*zGrid) = state(iLAPW)
383 : END DO
384 : !$acc end parallel loop
385 : else
386 262132500 : ext_grid = cmplx_0
387 :
388 10005454 : DO iLAPW = 1, lapw%nv(iSpin)
389 9938840 : xGrid = MODULO(lapw%gvec(1, iLAPW, iSpin), this%dimensions(1))
390 9938840 : yGrid = MODULO(lapw%gvec(2, iLAPW, iSpin), this%dimensions(2))
391 9938840 : zGrid = MODULO(lapw%gvec(3, iLAPW, iSpin), this%dimensions(3))
392 10005454 : ext_grid(xGrid + this%dimensions(1)*yGrid + (this%dimensions(1)*this%dimensions(2))*zGrid) = state(iLAPW)
393 : END DO
394 : endif
395 68000 : end SUBROUTINE put_cmplx_on_external_grid
396 :
397 8018 : SUBROUTINE fillStateIndexArray(this, lapw, ispin, indexArray)
398 : USE m_types_lapw
399 : IMPLICIT NONE
400 : CLASS(t_fftGrid), INTENT(INOUT) :: this
401 : TYPE(t_lapw), INTENT(IN) :: lapw
402 : INTEGER, INTENT(IN) :: iSpin
403 : INTEGER, INTENT(INOUT) :: indexArray(lapw%nv(ispin))
404 :
405 : INTEGER :: xGrid, yGrid, zGrid, layerDim, iLAPW
406 :
407 8018 : layerDim = this%dimensions(1)*this%dimensions(2)
408 :
409 1049151 : DO iLAPW = 1, lapw%nv(iSpin)
410 1041133 : xGrid = MODULO(lapw%gvec(1, iLAPW, iSpin), this%dimensions(1))
411 1041133 : yGrid = MODULO(lapw%gvec(2, iLAPW, iSpin), this%dimensions(2))
412 1041133 : zGrid = MODULO(lapw%gvec(3, iLAPW, iSpin), this%dimensions(3))
413 1049151 : indexArray(iLAPW) = xGrid + this%dimensions(1)*yGrid + layerDim*zGrid
414 : END DO
415 8018 : END SUBROUTINE fillStateIndexArray
416 :
417 7324 : SUBROUTINE fillFieldSphereIndexArray(this, stars, gCutoff, indexArray)
418 : USE m_types_stars
419 : IMPLICIT NONE
420 : CLASS(t_fftGrid), INTENT(IN) :: this
421 : TYPE(t_stars), INTENT(IN) :: stars
422 : REAL, INTENT(IN) :: gCutoff
423 : INTEGER, ALLOCATABLE, INTENT(INOUT) :: indexArray(:)
424 :
425 : INTEGER :: x, y, z, iStar
426 : INTEGER :: xGrid, yGrid, zGrid, layerDim, tempArrayIndex
427 7324 : INTEGER :: tempArray((2*stars%mx1+1)*(2*stars%mx2+1)*(2*stars%mx3+1))
428 :
429 7324 : layerDim = this%dimensions(1)*this%dimensions(2)
430 :
431 7324 : tempArrayIndex = 0
432 169900 : DO z = -stars%mx3, stars%mx3
433 162576 : zGrid = MODULO(z, this%dimensions(3))
434 3513976 : DO y = -stars%mx2, stars%mx2
435 3344076 : yGrid = MODULO(y, this%dimensions(2))
436 76111660 : DO x = -stars%mx1, stars%mx1
437 72605008 : iStar = stars%ig(x, y, z)
438 72605008 : IF (iStar .EQ. 0) CYCLE
439 25668520 : IF (stars%sk3(iStar) .GT. gCutoff) CYCLE
440 10123452 : xGrid = MODULO(x, this%dimensions(1))
441 10123452 : tempArrayIndex = tempArrayIndex + 1
442 75949084 : tempArray(tempArrayIndex) = xGrid + this%dimensions(1)*yGrid + layerDim*zGrid
443 : END DO
444 : END DO
445 : END DO
446 :
447 7324 : IF(ALLOCATED(indexArray)) DEALLOCATE (indexArray)
448 21972 : ALLOCATE(indexArray(tempArrayIndex))
449 :
450 10130776 : indexArray(1:tempArrayIndex) = tempArray(1:tempArrayIndex)
451 :
452 7324 : END SUBROUTINE fillFieldSphereIndexArray
453 :
454 0 : COMPLEX FUNCTION getElement(this, x, y, z)
455 : IMPLICIT NONE
456 : CLASS(t_fftGrid), INTENT(IN) :: this
457 : INTEGER, INTENT(IN) :: x, y, z
458 :
459 : INTEGER :: xGrid, yGrid, zGrid, layerDim
460 :
461 0 : layerDim = this%dimensions(1)*this%dimensions(2)
462 :
463 0 : xGrid = MODULO(x, this%dimensions(1))
464 0 : yGrid = MODULO(y, this%dimensions(2))
465 0 : zGrid = MODULO(z, this%dimensions(3))
466 :
467 0 : getElement = this%grid(xGrid + this%dimensions(1)*yGrid + layerDim*zGrid)
468 :
469 0 : END FUNCTION getElement
470 :
471 0 : SUBROUTINE getRealPartOfGrid(this, realGrid)
472 : IMPLICIT NONE
473 : CLASS(t_fftGrid), INTENT(IN) :: this
474 : REAL, INTENT(INOUT) :: realGrid(:)
475 :
476 0 : realGrid(:) = REAL(this%grid(:))
477 :
478 0 : END SUBROUTINE getRealPartOfGrid
479 :
480 137282 : subroutine free(fftGrid)
481 : implicit none
482 : TYPE(t_fftGrid), INTENT(INOUT) :: fftGrid
483 :
484 549128 : fftGrid%extent = -1
485 549128 : fftGrid%dimensions = -1
486 137282 : fftGrid%gridLength = -1
487 :
488 137282 : if(allocated(fftGrid%grid)) deallocate(fftGrid%grid)
489 137282 : end subroutine free
490 :
491 36382 : function calc_extent(cell, sym, gCutoff, gzCutoff) result(mxx)
492 : USE m_constants
493 : USE m_boxdim
494 : USE m_spgrot
495 : USE m_ifft
496 : USE m_types_cell
497 : USE m_types_sym
498 : IMPLICIT NONE
499 :
500 : TYPE(t_cell), INTENT(IN) :: cell
501 : TYPE(t_sym), INTENT(IN) :: sym
502 : REAL, INTENT(IN) :: gCutoff
503 : REAL, OPTIONAL, INTENT(IN):: gzCutoff
504 :
505 18191 : INTEGER, ALLOCATABLE :: ig(:, :, :)
506 :
507 : INTEGER :: k1, k2, k3, i, j
508 18191 : INTEGER :: mxx(3), kVec(3), kRot(3, sym%nop), inv_du(sym%nop), tempDim(3)
509 : REAL :: gCutoffSquared, gSquared
510 : REAL :: arltv(3), g(3)
511 : LOGICAL :: l_insph
512 :
513 18191 : CALL boxdim(cell%bmat, arltv(1), arltv(2), arltv(3))
514 :
515 72764 : tempDim(:) = INT(gCutoff/arltv(:)) + 1
516 18191 : IF (PRESENT(gzCutoff)) tempDim(3) = INT(gzCutoff/arltv(3)) + 1
517 :
518 242567 : DO i = 1, sym%nop
519 242567 : inv_du(i) = i ! dummy array for spgrot
520 : END DO
521 :
522 : ALLOCATE (ig(-tempDim(1):tempDim(1), &
523 : -tempDim(2):tempDim(2), &
524 93151057 : -tempDim(3):tempDim(3)), source=0)
525 :
526 72764 : mxx(:) = 0
527 18191 : gCutoffSquared = gCutoff*gCutoff
528 305324 : DO k1 = tempDim(1), -tempDim(1), -1
529 287133 : kVec(1) = k1
530 5021315 : DO k2 = tempDim(2), -tempDim(2), -1
531 4715991 : kVec(2) = k2
532 92567311 : DO k3 = tempDim(3), -tempDim(3), -1
533 87564187 : IF (ig(k1, k2, k3) .NE. 0) CYCLE
534 :
535 71547115 : kVec(3) = k3
536 :
537 286188460 : DO i = 1, 3
538 930112495 : g(i) = DOT_PRODUCT(cell%bmat(:, i), kVec(:)) ! loop to be replaced by MATMUL call later.
539 : END DO
540 :
541 71547115 : gSquared = g(1)**2 + g(2)**2 + g(3)**2
542 71547115 : l_insph = gSquared .LE. gCutoffSquared
543 71547115 : IF (PRESENT(gzCutoff)) l_insph = (g(1)**2 + g(2)**2.LE.gCutoffSquared).AND.(ABS(g(3)).LE.gzCutoff)
544 :
545 76263106 : IF (l_insph) THEN
546 13433205 : CALL spgrot(sym%nop, .true., sym%mrot, sym%tau, inv_du, kVec, kRot)
547 51590401 : DO i = 1, sym%nop
548 152628784 : do j = 1,3
549 152628784 : mxx(j) = max(mxx(j), kRot(j,i))
550 : enddo
551 51590401 : ig(kRot(1, i), kRot(2, i), kRot(3, i)) = 1
552 : END DO
553 : END IF
554 : END DO
555 : END DO
556 : END DO
557 18191 : END function calc_extent
558 :
559 0 : function calc_fft_dim(cell, sym, gCutoff, gzCutoff) result(dims)
560 : USE m_ifft
561 : USE m_types_cell
562 : USE m_types_sym
563 : implicit none
564 : TYPE(t_cell), INTENT(IN) :: cell
565 : TYPE(t_sym), INTENT(IN) :: sym
566 : REAL, INTENT(IN) :: gCutoff
567 : REAL, OPTIONAL, INTENT(IN):: gzCutoff
568 : integer :: dims(3)
569 : integer :: i
570 :
571 0 : if (.not.present(gzCutoff)) then
572 0 : dims = 2* calc_extent(cell, sym, gCutoff) + 1
573 : else
574 0 : dims = 2* calc_extent(cell, sym, gCutoff, gzCutoff) + 1
575 : end if
576 :
577 0 : do i = 1,3
578 0 : dims(i) = next_optimal_fft_size(dims(i))
579 : enddo
580 0 : end function calc_fft_dim
581 :
582 492388 : END MODULE m_types_fftGrid
|
