Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief This module returns additional info on PW grids
10 : !> \par History
11 : !> JGH (09-06-2007) : Created from pw_grids
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE pw_grid_info
15 :
16 : USE fft_tools, ONLY: FFT_RADIX_NEXT,&
17 : FFT_RADIX_NEXT_ODD,&
18 : fft_radix_operations
19 : USE kinds, ONLY: dp
20 : USE mathconstants, ONLY: twopi
21 : USE pw_grid_types, ONLY: pw_grid_type
22 : #include "../base/base_uses.f90"
23 :
24 : IMPLICIT NONE
25 :
26 : PRIVATE
27 : PUBLIC :: pw_find_cutoff, pw_grid_init_setup, pw_grid_bounds_from_n, pw_grid_n_for_fft
28 :
29 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_grid_info'
30 :
31 : CONTAINS
32 :
33 : ! **************************************************************************************************
34 : !> \brief ...
35 : !> \param hmat ...
36 : !> \param cutoff ...
37 : !> \param spherical ...
38 : !> \param odd ...
39 : !> \param fft_usage ...
40 : !> \param ncommensurate ...
41 : !> \param icommensurate ...
42 : !> \param ref_grid ...
43 : !> \param n_orig ...
44 : !> \return ...
45 : ! **************************************************************************************************
46 35774 : FUNCTION pw_grid_init_setup(hmat, cutoff, spherical, odd, fft_usage, ncommensurate, &
47 35774 : icommensurate, ref_grid, n_orig) RESULT(n)
48 :
49 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat
50 : REAL(KIND=dp), INTENT(IN) :: cutoff
51 : LOGICAL, INTENT(IN) :: spherical, odd, fft_usage
52 : INTEGER, INTENT(IN) :: ncommensurate, icommensurate
53 : TYPE(pw_grid_type), INTENT(IN), OPTIONAL :: ref_grid
54 : INTEGER, INTENT(IN), OPTIONAL :: n_orig(3)
55 : INTEGER, DIMENSION(3) :: n
56 :
57 : INTEGER :: my_icommensurate
58 :
59 35774 : IF (ncommensurate > 0) THEN
60 1114 : my_icommensurate = icommensurate
61 1114 : CPASSERT(icommensurate > 0)
62 1114 : CPASSERT(icommensurate <= ncommensurate)
63 : ELSE
64 : my_icommensurate = 0
65 : END IF
66 :
67 1114 : IF (my_icommensurate > 1) THEN
68 840 : CPASSERT(PRESENT(ref_grid))
69 4200 : n = ref_grid%npts/2**(my_icommensurate - 1)
70 4200 : CPASSERT(ALL(ref_grid%npts == n*2**(my_icommensurate - 1)))
71 3360 : CPASSERT(ALL(pw_grid_n_for_fft(n) == n))
72 : ELSE
73 : n = pw_grid_find_n(hmat, cutoff=cutoff, fft_usage=fft_usage, ncommensurate=ncommensurate, &
74 34934 : spherical=spherical, odd=odd, n_orig=n_orig)
75 : END IF
76 :
77 35774 : END FUNCTION pw_grid_init_setup
78 :
79 : ! **************************************************************************************************
80 : !> \brief returns the n needed for the grid with all the given constraints
81 : !> \param hmat ...
82 : !> \param cutoff ...
83 : !> \param fft_usage ...
84 : !> \param spherical ...
85 : !> \param odd ...
86 : !> \param ncommensurate ...
87 : !> \param n_orig ...
88 : !> \return ...
89 : !> \author fawzi
90 : ! **************************************************************************************************
91 34934 : FUNCTION pw_grid_find_n(hmat, cutoff, fft_usage, spherical, odd, ncommensurate, &
92 34934 : n_orig) RESULT(n)
93 :
94 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat
95 : REAL(KIND=dp), INTENT(IN) :: cutoff
96 : LOGICAL, INTENT(IN) :: fft_usage, spherical, odd
97 : INTEGER, INTENT(IN) :: ncommensurate
98 : INTEGER, INTENT(IN), OPTIONAL :: n_orig(3)
99 : INTEGER, DIMENSION(3) :: n
100 :
101 : INTEGER :: idir, my_icommensurate, &
102 : my_ncommensurate, nsubgrid, &
103 : nsubgrid_new, ntest(3), t_icommensurate
104 : LOGICAL :: ftest, subgrid_is_OK
105 :
106 : ! ncommensurate is the number of commensurate grids
107 : ! in order to have non-commensurate grids ncommensurate must be 0
108 : ! icommensurte is the level number of communensurate grids
109 : ! this implies that the number of grid points in each direction
110 : ! is k*2**(ncommensurate-icommensurate)
111 :
112 34934 : my_ncommensurate = ncommensurate
113 34934 : IF (my_ncommensurate > 0) THEN
114 : my_icommensurate = 1
115 : ELSE
116 : my_icommensurate = 0
117 : END IF
118 34660 : CPASSERT(my_icommensurate <= my_ncommensurate)
119 34934 : CPASSERT(my_icommensurate > 0 .OR. my_ncommensurate <= 0)
120 34934 : CPASSERT(my_ncommensurate >= 0)
121 :
122 34934 : IF (PRESENT(n_orig)) THEN
123 9592 : n = n_orig
124 : ELSE
125 32536 : CPASSERT(cutoff > 0.0_dp)
126 32536 : n = pw_grid_n_from_cutoff(hmat, cutoff)
127 : END IF
128 :
129 34934 : IF (fft_usage) THEN
130 139736 : n = pw_grid_n_for_fft(n, odd=odd)
131 :
132 34934 : IF (.NOT. spherical) THEN
133 130416 : ntest = n
134 :
135 32604 : IF (my_ncommensurate > 0) THEN
136 1096 : DO idir = 1, 3
137 274 : DO
138 : ! find valid radix >= ntest
139 2260 : CALL fft_radix_operations(ntest(idir), n(idir), FFT_RADIX_NEXT)
140 : ! check every subgrid of n
141 2260 : subgrid_is_OK = .TRUE.
142 4780 : DO t_icommensurate = 1, my_ncommensurate - 1
143 3958 : nsubgrid = n(idir)/2**(my_ncommensurate - t_icommensurate)
144 3958 : CALL fft_radix_operations(nsubgrid, nsubgrid_new, FFT_RADIX_NEXT)
145 : subgrid_is_OK = (nsubgrid == nsubgrid_new) .AND. &
146 3958 : (MODULO(n(idir), 2**(my_ncommensurate - t_icommensurate)) == 0)
147 822 : IF (.NOT. subgrid_is_OK) EXIT
148 : END DO
149 2260 : IF (subgrid_is_OK) THEN
150 : EXIT
151 : ELSE
152 : ! subgrid wasn't OK, increment ntest and try again
153 1438 : ntest(idir) = n(idir) + 1
154 : END IF
155 : END DO
156 : END DO
157 : END IF
158 : END IF
159 : ELSE
160 : ! without a cutoff and HALFSPACE we have to be sure that there is
161 : ! a negative counterpart to every g vector (-> odd number of grid points)
162 0 : IF (odd) n = n + MOD(n + 1, 2)
163 :
164 : END IF
165 :
166 : ! final check if all went fine ...
167 2604 : IF (my_ncommensurate > 0) THEN
168 1388 : DO my_icommensurate = 1, my_ncommensurate
169 5570 : ftest = ANY(MODULO(n, 2**(my_ncommensurate - my_icommensurate)) .NE. 0)
170 1388 : CPASSERT(.NOT. ftest)
171 : END DO
172 : END IF
173 :
174 34934 : END FUNCTION pw_grid_find_n
175 :
176 : ! **************************************************************************************************
177 : !> \brief returns the closest number of points >= n, on which you can perform
178 : !> ffts
179 : !> \param n the minimum number of points you want
180 : !> \param odd if the number has to be odd
181 : !> \return ...
182 : !> \author fawzi
183 : !> \note
184 : !> result<=n
185 : ! **************************************************************************************************
186 36210 : FUNCTION pw_grid_n_for_fft(n, odd) RESULT(nout)
187 : INTEGER, DIMENSION(3), INTENT(in) :: n
188 : LOGICAL, INTENT(in), OPTIONAL :: odd
189 : INTEGER, DIMENSION(3) :: nout
190 :
191 : LOGICAL :: my_odd
192 :
193 36210 : my_odd = .FALSE.
194 36210 : IF (PRESENT(odd)) my_odd = odd
195 144840 : CPASSERT(ALL(n >= 0))
196 36210 : IF (my_odd) THEN
197 468 : CALL fft_radix_operations(n(1), nout(1), FFT_RADIX_NEXT_ODD)
198 468 : CALL fft_radix_operations(n(2), nout(2), FFT_RADIX_NEXT_ODD)
199 468 : CALL fft_radix_operations(n(3), nout(3), FFT_RADIX_NEXT_ODD)
200 : ELSE
201 35742 : CALL fft_radix_operations(n(1), nout(1), FFT_RADIX_NEXT)
202 35742 : CALL fft_radix_operations(n(2), nout(2), FFT_RADIX_NEXT)
203 35742 : CALL fft_radix_operations(n(3), nout(3), FFT_RADIX_NEXT)
204 : END IF
205 :
206 36210 : END FUNCTION pw_grid_n_for_fft
207 :
208 : ! **************************************************************************************************
209 : !> \brief Find the number of points that give at least the requested cutoff
210 : !> \param hmat ...
211 : !> \param cutoff ...
212 : !> \return ...
213 : !> \par History
214 : !> JGH (21-12-2000) : Simplify parameter list, bounds will be global
215 : !> JGH ( 8-01-2001) : Add check to FFT allowd grids (this now depends
216 : !> on the FFT library.
217 : !> Should the pw_grid_type have a reference to the FFT
218 : !> library ?
219 : !> JGH (28-02-2001) : Only do conditional check for FFT
220 : !> JGH (21-05-2002) : Optimise code, remove orthorhombic special case
221 : !> \author apsi
222 : !> Christopher Mundy
223 : ! **************************************************************************************************
224 32536 : FUNCTION pw_grid_n_from_cutoff(hmat, cutoff) RESULT(n)
225 :
226 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat
227 : REAL(KIND=dp), INTENT(IN) :: cutoff
228 : INTEGER, DIMENSION(3) :: n
229 :
230 : INTEGER :: i
231 : REAL(KIND=dp) :: alat(3)
232 :
233 130144 : DO i = 1, 3
234 422968 : alat(i) = SUM(hmat(:, i)**2)
235 : END DO
236 130144 : CPASSERT(ALL(alat /= 0._dp))
237 130144 : n = 2*FLOOR(SQRT(2.0_dp*cutoff*alat)/twopi) + 1
238 :
239 32536 : END FUNCTION pw_grid_n_from_cutoff
240 :
241 : ! **************************************************************************************************
242 : !> \brief returns the bounds that distribute n points evenly around 0
243 : !> \param npts the number of points in each direction
244 : !> \return ...
245 : !> \author fawzi
246 : ! **************************************************************************************************
247 33426 : FUNCTION pw_grid_bounds_from_n(npts) RESULT(bounds)
248 : INTEGER, DIMENSION(3), INTENT(in) :: npts
249 : INTEGER, DIMENSION(2, 3) :: bounds
250 :
251 133704 : bounds(1, :) = -npts/2
252 133704 : bounds(2, :) = bounds(1, :) + npts - 1
253 :
254 33426 : END FUNCTION pw_grid_bounds_from_n
255 :
256 : ! **************************************************************************************************
257 : !> \brief Given a grid and a box, calculate the corresponding cutoff
258 : !> *** This routine calculates the cutoff in MOMENTUM UNITS! ***
259 : !> \param npts ...
260 : !> \param h_inv ...
261 : !> \return ...
262 : !> \par History
263 : !> JGH (20-12-2000) : Deleted some strange comments
264 : !> \author apsi
265 : !> Christopher Mundy
266 : !> \note
267 : !> This routine is local. It works independent from the distribution
268 : !> of PW on processors.
269 : !> npts is the grid size for the full box.
270 : ! **************************************************************************************************
271 3409 : FUNCTION pw_find_cutoff(npts, h_inv) RESULT(cutoff)
272 :
273 : INTEGER, DIMENSION(:), INTENT(IN) :: npts
274 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_inv
275 : REAL(KIND=dp) :: cutoff
276 :
277 : REAL(KIND=dp) :: gcut, gdum(3), length
278 :
279 : ! compute 2*pi*h_inv^t*g where g = (nmax[1],0,0)
280 :
281 13636 : gdum(:) = twopi*h_inv(1, :)*REAL((npts(1) - 1)/2, KIND=dp)
282 3409 : length = SQRT(gdum(1)**2 + gdum(2)**2 + gdum(3)**2)
283 3409 : gcut = length
284 :
285 : ! compute 2*pi*h_inv^t*g where g = (0,nmax[2],0)
286 13636 : gdum(:) = twopi*h_inv(2, :)*REAL((npts(2) - 1)/2, KIND=dp)
287 3409 : length = SQRT(gdum(1)**2 + gdum(2)**2 + gdum(3)**2)
288 3409 : gcut = MIN(gcut, length)
289 :
290 : ! compute 2*pi*h_inv^t*g where g = (0,0,nmax[3])
291 13636 : gdum(:) = twopi*h_inv(3, :)*REAL((npts(3) - 1)/2, KIND=dp)
292 3409 : length = SQRT(gdum(1)**2 + gdum(2)**2 + gdum(3)**2)
293 3409 : gcut = MIN(gcut, length)
294 :
295 3409 : cutoff = gcut - 1.e-8_dp
296 :
297 3409 : END FUNCTION pw_find_cutoff
298 :
299 : END MODULE pw_grid_info
300 :
|
