Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief methods of pw_env that have dependence on qs_env
10 : !> \par History
11 : !> 10.2002 created [fawzi]
12 : !> JGH (22-Feb-03) PW grid options added
13 : !> 04.2003 added rs grid pools [fawzi]
14 : !> 02.2004 added commensurate grids
15 : !> \author Fawzi Mohamed
16 : ! **************************************************************************************************
17 : MODULE pw_env_methods
18 : USE ao_util, ONLY: exp_radius
19 : USE basis_set_types, ONLY: get_gto_basis_set,&
20 : gto_basis_set_type
21 : USE cell_types, ONLY: cell_type
22 : USE cp_control_types, ONLY: dft_control_type
23 : USE cp_log_handling, ONLY: cp_get_default_logger,&
24 : cp_logger_type
25 : USE cp_output_handling, ONLY: cp_p_file,&
26 : cp_print_key_finished_output,&
27 : cp_print_key_should_output,&
28 : cp_print_key_unit_nr
29 : USE cp_realspace_grid_init, ONLY: init_input_type
30 : USE cube_utils, ONLY: destroy_cube_info,&
31 : init_cube_info,&
32 : return_cube_max_iradius
33 : USE d3_poly, ONLY: init_d3_poly_module
34 : USE dct, ONLY: neumannX,&
35 : neumannXY,&
36 : neumannXYZ,&
37 : neumannXZ,&
38 : neumannY,&
39 : neumannYZ,&
40 : neumannZ,&
41 : setup_dct_pw_grids
42 : USE dielectric_types, ONLY: derivative_cd3,&
43 : derivative_cd5,&
44 : derivative_cd7,&
45 : rho_dependent
46 : USE fft_tools, ONLY: init_fft_scratch_pool
47 : USE gaussian_gridlevels, ONLY: destroy_gaussian_gridlevel,&
48 : gaussian_gridlevel,&
49 : init_gaussian_gridlevel
50 : USE input_constants, ONLY: do_method_gapw,&
51 : do_method_gapw_xc,&
52 : do_method_gpw,&
53 : do_method_lrigpw,&
54 : do_method_ofgpw,&
55 : do_method_rigpw,&
56 : xc_vdw_fun_nonloc
57 : USE input_section_types, ONLY: section_get_ival,&
58 : section_vals_get,&
59 : section_vals_get_subs_vals,&
60 : section_vals_type,&
61 : section_vals_val_get
62 : USE kinds, ONLY: dp
63 : USE message_passing, ONLY: mp_para_env_type
64 : USE ps_implicit_types, ONLY: MIXED_BC,&
65 : MIXED_PERIODIC_BC,&
66 : NEUMANN_BC,&
67 : PERIODIC_BC
68 : USE ps_wavelet_types, ONLY: WAVELET0D,&
69 : WAVELET2D,&
70 : WAVELET3D
71 : USE pw_env_types, ONLY: pw_env_type
72 : USE pw_grid_info, ONLY: pw_grid_init_setup
73 : USE pw_grid_types, ONLY: FULLSPACE,&
74 : HALFSPACE,&
75 : pw_grid_type
76 : USE pw_grids, ONLY: do_pw_grid_blocked_false,&
77 : pw_grid_change,&
78 : pw_grid_create,&
79 : pw_grid_release
80 : USE pw_poisson_methods, ONLY: pw_poisson_set
81 : USE pw_poisson_read_input, ONLY: pw_poisson_read_parameters
82 : USE pw_poisson_types, ONLY: pw_poisson_analytic,&
83 : pw_poisson_implicit,&
84 : pw_poisson_mt,&
85 : pw_poisson_multipole,&
86 : pw_poisson_none,&
87 : pw_poisson_parameter_type,&
88 : pw_poisson_periodic,&
89 : pw_poisson_wavelet
90 : USE pw_pool_types, ONLY: pw_pool_create,&
91 : pw_pool_p_type,&
92 : pw_pool_release,&
93 : pw_pools_dealloc
94 : USE qs_cneo_types, ONLY: cneo_potential_type
95 : USE qs_dispersion_types, ONLY: qs_dispersion_type
96 : USE qs_environment_types, ONLY: get_qs_env,&
97 : qs_environment_type
98 : USE qs_kind_types, ONLY: get_qs_kind,&
99 : qs_kind_type
100 : USE qs_rho0_types, ONLY: get_rho0_mpole,&
101 : rho0_mpole_type
102 : USE realspace_grid_types, ONLY: &
103 : realspace_grid_desc_p_type, realspace_grid_desc_type, realspace_grid_input_type, &
104 : realspace_grid_type, rs_grid_create, rs_grid_create_descriptor, rs_grid_print, &
105 : rs_grid_release, rs_grid_release_descriptor
106 : USE xc_input_constants, ONLY: &
107 : xc_deriv_collocate, xc_deriv_nn10_smooth, xc_deriv_nn50_smooth, xc_deriv_pw, &
108 : xc_deriv_spline2, xc_deriv_spline2_smooth, xc_deriv_spline3, xc_deriv_spline3_smooth, &
109 : xc_rho_nn10, xc_rho_nn50, xc_rho_no_smooth, xc_rho_spline2_smooth, xc_rho_spline3_smooth
110 : #include "./base/base_uses.f90"
111 :
112 : IMPLICIT NONE
113 : PRIVATE
114 :
115 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
116 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pw_env_methods'
117 :
118 : PUBLIC :: pw_env_create, pw_env_rebuild
119 :
120 : ! **************************************************************************************************
121 :
122 : CONTAINS
123 :
124 : ! **************************************************************************************************
125 : !> \brief creates a pw_env, if qs_env is given calls pw_env_rebuild
126 : !> \param pw_env the pw_env that gets created
127 : !> \par History
128 : !> 10.2002 created [fawzi]
129 : !> \author Fawzi Mohamed
130 : ! **************************************************************************************************
131 8628 : SUBROUTINE pw_env_create(pw_env)
132 : TYPE(pw_env_type), POINTER :: pw_env
133 :
134 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_env_create'
135 :
136 : INTEGER :: handle
137 :
138 8628 : CALL timeset(routineN, handle)
139 :
140 8628 : CPASSERT(.NOT. ASSOCIATED(pw_env))
141 120792 : ALLOCATE (pw_env)
142 : NULLIFY (pw_env%pw_pools, pw_env%gridlevel_info, pw_env%poisson_env, &
143 : pw_env%cube_info, pw_env%rs_descs, pw_env%rs_grids, &
144 : pw_env%xc_pw_pool, pw_env%vdw_pw_pool, &
145 : pw_env%interp_section)
146 8628 : pw_env%auxbas_grid = -1
147 8628 : pw_env%ref_count = 1
148 :
149 8628 : CALL timestop(handle)
150 :
151 8628 : END SUBROUTINE pw_env_create
152 :
153 : ! **************************************************************************************************
154 : !> \brief rebuilds the pw_env data (necessary if cell or cutoffs change)
155 : !> \param pw_env the environment to rebuild
156 : !> \param qs_env the qs_env where to get the cell, cutoffs,...
157 : !> \param external_para_env ...
158 : !> \par History
159 : !> 10.2002 created [fawzi]
160 : !> \author Fawzi Mohamed
161 : ! **************************************************************************************************
162 10058 : SUBROUTINE pw_env_rebuild(pw_env, qs_env, external_para_env)
163 : TYPE(pw_env_type), POINTER :: pw_env
164 : TYPE(qs_environment_type), POINTER :: qs_env
165 : TYPE(mp_para_env_type), INTENT(IN), OPTIONAL, &
166 : TARGET :: external_para_env
167 :
168 : CHARACTER(len=*), PARAMETER :: routineN = 'pw_env_rebuild'
169 :
170 : CHARACTER(LEN=3) :: string
171 : INTEGER :: blocked_id, blocked_id_input, boundary_condition, grid_span, handle, i, &
172 : igrid_level, iounit, ncommensurate, ngrid_level, xc_deriv_method_id, xc_smooth_method_id
173 : INTEGER, DIMENSION(2) :: distribution_layout
174 : INTEGER, DIMENSION(3) :: higher_grid_layout
175 : LOGICAL :: do_io, efg_present, linres_present, odd, set_vdw_pool, should_output, &
176 : smooth_required, spherical, uf_grid, use_ref_cell
177 : REAL(KIND=dp) :: cutilev, rel_cutoff
178 10058 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radius
179 10058 : REAL(KIND=dp), DIMENSION(:), POINTER :: cutoff
180 : TYPE(cell_type), POINTER :: cell, cell_ref, my_cell
181 : TYPE(cp_logger_type), POINTER :: logger
182 : TYPE(dft_control_type), POINTER :: dft_control
183 : TYPE(mp_para_env_type), POINTER :: para_env
184 : TYPE(pw_grid_type), POINTER :: dct_pw_grid, mt_super_ref_grid, old_pw_grid, pw_grid, &
185 : super_ref_grid, vdw_grid, vdw_ref_grid, xc_super_ref_grid
186 : TYPE(pw_poisson_parameter_type) :: poisson_params
187 10058 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
188 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
189 10058 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
190 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
191 10058 : POINTER :: rs_descs
192 : TYPE(realspace_grid_input_type) :: input_settings
193 10058 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_grids
194 : TYPE(section_vals_type), POINTER :: efg_section, input, linres_section, &
195 : poisson_section, print_section, &
196 : rs_grid_section, xc_section
197 :
198 : ! a very small safety factor might be needed for roundoff issues
199 : ! e.g. radius being computed here as 12.998 (13) and 13.002 (14) during the collocation
200 : ! the latter can happen due to the lower precision in the computation of the radius in collocate
201 : ! parallel cost of rs_pw_transfer goes as safety_factor**3 so it is worthwhile keeping it tight
202 : ! Edit: Safety Factor was unused
203 :
204 10058 : CALL timeset(routineN, handle)
205 :
206 : !
207 : !
208 : ! Part one, deallocate old data if needed
209 : !
210 : !
211 10058 : NULLIFY (cutoff, cell, pw_grid, old_pw_grid, dft_control, qs_kind_set, &
212 10058 : pw_pools, rs_descs, para_env, cell_ref, vdw_ref_grid, &
213 10058 : mt_super_ref_grid, input, poisson_section, xc_super_ref_grid, &
214 10058 : dct_pw_grid, vdw_grid, super_ref_grid, my_cell, rs_grids, dispersion_env)
215 :
216 : CALL get_qs_env(qs_env=qs_env, &
217 : dft_control=dft_control, &
218 : qs_kind_set=qs_kind_set, &
219 : cell_ref=cell_ref, &
220 : cell=cell, &
221 : para_env=para_env, &
222 : input=input, &
223 10058 : dispersion_env=dispersion_env)
224 :
225 10058 : CPASSERT(ASSOCIATED(pw_env))
226 10058 : CPASSERT(pw_env%ref_count > 0)
227 10058 : CALL pw_pool_release(pw_env%vdw_pw_pool)
228 10058 : CALL pw_pool_release(pw_env%xc_pw_pool)
229 10058 : CALL pw_pools_dealloc(pw_env%pw_pools)
230 10058 : IF (ASSOCIATED(pw_env%rs_descs)) THEN
231 3972 : DO i = 1, SIZE(pw_env%rs_descs)
232 3972 : CALL rs_grid_release_descriptor(pw_env%rs_descs(i)%rs_desc)
233 : END DO
234 1452 : DEALLOCATE (pw_env%rs_descs)
235 : END IF
236 10058 : IF (ASSOCIATED(pw_env%rs_grids)) THEN
237 3972 : DO i = 1, SIZE(pw_env%rs_grids)
238 3972 : CALL rs_grid_release(pw_env%rs_grids(i))
239 : END DO
240 1452 : DEALLOCATE (pw_env%rs_grids)
241 : END IF
242 10058 : IF (ASSOCIATED(pw_env%gridlevel_info)) THEN
243 1452 : CALL destroy_gaussian_gridlevel(pw_env%gridlevel_info)
244 : ELSE
245 8606 : ALLOCATE (pw_env%gridlevel_info)
246 : END IF
247 :
248 10058 : IF (ASSOCIATED(pw_env%cube_info)) THEN
249 3972 : DO igrid_level = 1, SIZE(pw_env%cube_info)
250 3972 : CALL destroy_cube_info(pw_env%cube_info(igrid_level))
251 : END DO
252 1452 : DEALLOCATE (pw_env%cube_info)
253 : END IF
254 10058 : NULLIFY (pw_env%pw_pools, pw_env%cube_info)
255 :
256 : ! remove fft scratch pool, as it depends on pw_env mpi group handles
257 10058 : CALL init_fft_scratch_pool()
258 :
259 : !
260 : !
261 : ! Part two, setup the pw_grids
262 : !
263 : !
264 :
265 10058 : do_io = .TRUE.
266 10058 : IF (PRESENT(external_para_env)) THEN
267 1100 : para_env => external_para_env
268 : CPASSERT(ASSOCIATED(para_env))
269 1100 : do_io = .FALSE. !multiple MPI subgroups mess-up the output file
270 : END IF
271 : ! interpolation section
272 10058 : pw_env%interp_section => section_vals_get_subs_vals(input, "DFT%MGRID%INTERPOLATOR")
273 :
274 10058 : CALL get_qs_env(qs_env, use_ref_cell=use_ref_cell)
275 10058 : IF (use_ref_cell) THEN
276 60 : my_cell => cell_ref
277 : ELSE
278 9998 : my_cell => cell
279 : END IF
280 10058 : rel_cutoff = dft_control%qs_control%relative_cutoff
281 10058 : cutoff => dft_control%qs_control%e_cutoff
282 10058 : CALL section_vals_val_get(input, "DFT%XC%XC_GRID%USE_FINER_GRID", l_val=uf_grid)
283 10058 : ngrid_level = SIZE(cutoff)
284 :
285 : ! init gridlevel_info XXXXXXXXX setup mapping to the effective cutoff ?
286 : ! XXXXXXXXX the cutoff array here is more a 'wish-list'
287 : ! XXXXXXXXX same holds for radius
288 : print_section => section_vals_get_subs_vals(input, &
289 10058 : "PRINT%GRID_INFORMATION")
290 : CALL init_gaussian_gridlevel(pw_env%gridlevel_info, &
291 : ngrid_levels=ngrid_level, cutoff=cutoff, rel_cutoff=rel_cutoff, &
292 10058 : print_section=print_section)
293 : ! init pw_grids and pools
294 60338 : ALLOCATE (pw_pools(ngrid_level))
295 :
296 10058 : IF (dft_control%qs_control%commensurate_mgrids) THEN
297 274 : ncommensurate = ngrid_level
298 : ELSE
299 9784 : ncommensurate = 0
300 : END IF
301 : !
302 : ! If Tuckerman is present let's perform the set-up of the super-reference-grid
303 : !
304 10058 : cutilev = cutoff(1)
305 10058 : IF (dft_control%qs_control%pw_grid_opt%spherical) THEN
306 0 : grid_span = HALFSPACE
307 0 : spherical = .TRUE.
308 10058 : ELSE IF (dft_control%qs_control%pw_grid_opt%fullspace) THEN
309 10058 : grid_span = FULLSPACE
310 10058 : spherical = .FALSE.
311 : ELSE
312 0 : grid_span = HALFSPACE
313 0 : spherical = .FALSE.
314 : END IF
315 :
316 10058 : poisson_section => section_vals_get_subs_vals(qs_env%input, "DFT%POISSON")
317 : CALL setup_super_ref_grid(super_ref_grid, mt_super_ref_grid, &
318 : xc_super_ref_grid, cutilev, grid_span, spherical, my_cell, para_env, &
319 : poisson_section, ncommensurate, uf_grid=uf_grid, &
320 10058 : print_section=print_section)
321 10058 : old_pw_grid => super_ref_grid
322 10058 : IF (.NOT. ASSOCIATED(mt_super_ref_grid)) vdw_ref_grid => super_ref_grid
323 : !
324 : ! Setup of the multi-grid pw_grid and pw_pools
325 : !
326 :
327 10058 : IF (do_io) THEN
328 8958 : logger => cp_get_default_logger()
329 8958 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.Log')
330 : ELSE
331 1100 : iounit = 0
332 : END IF
333 :
334 10058 : IF (dft_control%qs_control%pw_grid_opt%spherical) THEN
335 0 : grid_span = HALFSPACE
336 0 : spherical = .TRUE.
337 0 : odd = .TRUE.
338 10058 : ELSE IF (dft_control%qs_control%pw_grid_opt%fullspace) THEN
339 10058 : grid_span = FULLSPACE
340 10058 : spherical = .FALSE.
341 10058 : odd = .FALSE.
342 : ELSE
343 0 : grid_span = HALFSPACE
344 0 : spherical = .FALSE.
345 0 : odd = .TRUE.
346 : END IF
347 :
348 : ! use input suggestion for blocked
349 10058 : blocked_id_input = dft_control%qs_control%pw_grid_opt%blocked
350 :
351 : ! methods that require smoothing or nearest neighbor have to use a plane distributed setup
352 : ! find the xc properties (FIXME this could miss other xc sections that operate on the grid ...)
353 10058 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
354 10058 : xc_deriv_method_id = section_get_ival(xc_section, "XC_GRID%XC_DERIV")
355 10058 : xc_smooth_method_id = section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO")
356 10058 : smooth_required = .FALSE.
357 : SELECT CASE (xc_deriv_method_id)
358 : CASE (xc_deriv_pw, xc_deriv_collocate, xc_deriv_spline3, xc_deriv_spline2)
359 80 : smooth_required = smooth_required .OR. .FALSE.
360 : CASE (xc_deriv_spline2_smooth, &
361 : xc_deriv_spline3_smooth, xc_deriv_nn10_smooth, xc_deriv_nn50_smooth)
362 80 : smooth_required = smooth_required .OR. .TRUE.
363 : CASE DEFAULT
364 10058 : CPABORT("")
365 : END SELECT
366 : SELECT CASE (xc_smooth_method_id)
367 : CASE (xc_rho_no_smooth)
368 42 : smooth_required = smooth_required .OR. .FALSE.
369 : CASE (xc_rho_spline2_smooth, xc_rho_spline3_smooth, xc_rho_nn10, xc_rho_nn50)
370 42 : smooth_required = smooth_required .OR. .TRUE.
371 : CASE DEFAULT
372 10058 : CPABORT("")
373 : END SELECT
374 : ! EPR, NMR, EFG can require splines. If the linres/EFG section is present we assume
375 : ! it could be on and splines might be used (not quite sure if this is due to their use of splines or something else)
376 : linres_section => section_vals_get_subs_vals(section_vals=input, &
377 10058 : subsection_name="PROPERTIES%LINRES")
378 10058 : CALL section_vals_get(linres_section, explicit=linres_present)
379 10058 : IF (linres_present) THEN
380 282 : smooth_required = smooth_required .OR. .TRUE.
381 : END IF
382 :
383 : efg_section => section_vals_get_subs_vals(section_vals=input, &
384 10058 : subsection_name="DFT%PRINT%ELECTRIC_FIELD_GRADIENT")
385 10058 : CALL section_vals_get(efg_section, explicit=efg_present)
386 10058 : IF (efg_present) THEN
387 2 : smooth_required = smooth_required .OR. .TRUE.
388 : END IF
389 :
390 40222 : DO igrid_level = 1, ngrid_level
391 30164 : cutilev = cutoff(igrid_level)
392 :
393 : ! the whole of QS seems to work fine with either blocked/non-blocked distribution in g-space
394 : ! the default choice should be made free
395 30164 : blocked_id = blocked_id_input
396 :
397 90492 : distribution_layout = dft_control%qs_control%pw_grid_opt%distribution_layout
398 :
399 : ! qmmm does not support a ray distribution
400 : ! FIXME ... check if a plane distributed lower grid is sufficient
401 30164 : IF (qs_env%qmmm) THEN
402 3606 : distribution_layout = [para_env%num_pe, 1]
403 : END IF
404 :
405 : ! If splines are required
406 : ! FIXME.... should only be true for the highest grid
407 30164 : IF (smooth_required) THEN
408 4158 : distribution_layout = [para_env%num_pe, 1]
409 : END IF
410 :
411 30164 : IF (igrid_level == 1) THEN
412 10058 : IF (ASSOCIATED(old_pw_grid)) THEN
413 : CALL pw_grid_create(pw_grid, para_env, my_cell%hmat, grid_span=grid_span, &
414 : cutoff=cutilev, &
415 : spherical=spherical, odd=odd, fft_usage=.TRUE., &
416 : ncommensurate=ncommensurate, icommensurate=igrid_level, &
417 : blocked=do_pw_grid_blocked_false, &
418 : ref_grid=old_pw_grid, &
419 : rs_dims=distribution_layout, &
420 1222 : iounit=iounit)
421 1222 : old_pw_grid => pw_grid
422 : ELSE
423 : CALL pw_grid_create(pw_grid, para_env, my_cell%hmat, grid_span=grid_span, &
424 : cutoff=cutilev, &
425 : spherical=spherical, odd=odd, fft_usage=.TRUE., &
426 : ncommensurate=ncommensurate, icommensurate=igrid_level, &
427 : blocked=blocked_id, &
428 : rs_dims=distribution_layout, &
429 8836 : iounit=iounit)
430 8836 : old_pw_grid => pw_grid
431 : END IF
432 : ELSE
433 : CALL pw_grid_create(pw_grid, para_env, my_cell%hmat, grid_span=grid_span, &
434 : cutoff=cutilev, &
435 : spherical=spherical, odd=odd, fft_usage=.TRUE., &
436 : ncommensurate=ncommensurate, icommensurate=igrid_level, &
437 : blocked=do_pw_grid_blocked_false, &
438 : ref_grid=old_pw_grid, &
439 : rs_dims=distribution_layout, &
440 20106 : iounit=iounit)
441 : END IF
442 :
443 : ! init pw_pools
444 30164 : NULLIFY (pw_pools(igrid_level)%pool)
445 30164 : CALL pw_pool_create(pw_pools(igrid_level)%pool, pw_grid=pw_grid)
446 :
447 40222 : CALL pw_grid_release(pw_grid)
448 :
449 : END DO
450 :
451 10058 : pw_env%pw_pools => pw_pools
452 :
453 : ! init auxbas_grid
454 40222 : DO i = 1, ngrid_level
455 40222 : IF (cutoff(i) == dft_control%qs_control%cutoff) pw_env%auxbas_grid = i
456 : END DO
457 :
458 : ! init xc_pool
459 10058 : IF (ASSOCIATED(xc_super_ref_grid)) THEN
460 : CALL pw_pool_create(pw_env%xc_pw_pool, &
461 4 : pw_grid=xc_super_ref_grid)
462 4 : CALL pw_grid_release(xc_super_ref_grid)
463 : ELSE
464 10054 : pw_env%xc_pw_pool => pw_pools(pw_env%auxbas_grid)%pool
465 10054 : CALL pw_env%xc_pw_pool%retain()
466 : END IF
467 :
468 : ! init vdw_pool
469 10058 : set_vdw_pool = .FALSE.
470 10058 : IF (ASSOCIATED(dispersion_env)) THEN
471 9724 : IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN
472 74 : IF (dispersion_env%pw_cutoff > 0._dp) set_vdw_pool = .TRUE.
473 : END IF
474 : END IF
475 : IF (set_vdw_pool) THEN
476 68 : CPASSERT(ASSOCIATED(old_pw_grid))
477 68 : IF (.NOT. ASSOCIATED(vdw_ref_grid)) vdw_ref_grid => old_pw_grid
478 68 : IF (iounit > 0) WRITE (iounit, "(/,T2,A)") "PW_GRID| Grid for non-local vdW functional"
479 : CALL pw_grid_create(vdw_grid, para_env, my_cell%hmat, grid_span=grid_span, &
480 : cutoff=dispersion_env%pw_cutoff, &
481 : spherical=spherical, odd=odd, fft_usage=.TRUE., &
482 : ncommensurate=0, icommensurate=0, &
483 : blocked=do_pw_grid_blocked_false, &
484 : ref_grid=vdw_ref_grid, &
485 : rs_dims=distribution_layout, &
486 68 : iounit=iounit)
487 68 : CALL pw_pool_create(pw_env%vdw_pw_pool, pw_grid=vdw_grid)
488 68 : CALL pw_grid_release(vdw_grid)
489 : ELSE
490 9990 : pw_env%vdw_pw_pool => pw_pools(pw_env%auxbas_grid)%pool
491 9990 : CALL pw_env%vdw_pw_pool%retain()
492 : END IF
493 :
494 10058 : IF (do_io) CALL cp_print_key_finished_output(iounit, logger, print_section, '')
495 :
496 : ! complete init of the poisson_env
497 10058 : IF (.NOT. ASSOCIATED(pw_env%poisson_env)) THEN
498 137696 : ALLOCATE (pw_env%poisson_env)
499 8606 : CALL pw_env%poisson_env%create()
500 : END IF
501 : ! poisson_section => section_vals_get_subs_vals(input, "DFT%POISSON")
502 :
503 10058 : CALL pw_poisson_read_parameters(poisson_section, poisson_params)
504 : CALL pw_poisson_set(pw_env%poisson_env, cell_hmat=my_cell%hmat, pw_pools=pw_env%pw_pools, &
505 : parameters=poisson_params, mt_super_ref_pw_grid=mt_super_ref_grid, &
506 10058 : dct_pw_grid=dct_pw_grid, use_level=pw_env%auxbas_grid)
507 10058 : CALL pw_grid_release(mt_super_ref_grid)
508 10058 : CALL pw_grid_release(dct_pw_grid)
509 : !
510 : ! If reference cell is present, then use pw_grid_change to keep bounds constant...
511 : ! do not re-init the Gaussian grid level (fix the gridlevel on which the pgf should go.
512 : !
513 10058 : IF (use_ref_cell) THEN
514 260 : DO igrid_level = 1, SIZE(pw_pools)
515 260 : CALL pw_grid_change(cell%hmat, pw_pools(igrid_level)%pool%pw_grid)
516 : END DO
517 60 : IF (set_vdw_pool) CALL pw_grid_change(cell%hmat, pw_env%vdw_pw_pool%pw_grid)
518 60 : CALL pw_poisson_read_parameters(poisson_section, poisson_params)
519 : CALL pw_poisson_set(pw_env%poisson_env, cell_hmat=cell%hmat, pw_pools=pw_env%pw_pools, &
520 : parameters=poisson_params, mt_super_ref_pw_grid=mt_super_ref_grid, &
521 60 : dct_pw_grid=dct_pw_grid, use_level=pw_env%auxbas_grid)
522 : END IF
523 :
524 10058 : IF ((poisson_params%ps_implicit_params%boundary_condition == MIXED_PERIODIC_BC) .OR. &
525 : (poisson_params%ps_implicit_params%boundary_condition == MIXED_BC)) THEN
526 : pw_env%poisson_env%parameters%dbc_params%do_dbc_cube = &
527 : BTEST(cp_print_key_should_output(logger%iter_info, input, &
528 38 : "DFT%PRINT%IMPLICIT_PSOLVER%DIRICHLET_BC_CUBE"), cp_p_file)
529 : END IF
530 : ! setup dct_pw_grid (an extended pw_grid) for Discrete Cosine Transformation (DCT)
531 10058 : IF ((poisson_params%ps_implicit_params%boundary_condition == NEUMANN_BC) .OR. &
532 : (poisson_params%ps_implicit_params%boundary_condition == MIXED_BC)) THEN
533 : CALL setup_dct_pw_grids(pw_env%poisson_env%pw_pools(pw_env%poisson_env%pw_level)%pool%pw_grid, &
534 : my_cell%hmat, poisson_params%ps_implicit_params%neumann_directions, &
535 22 : pw_env%poisson_env%dct_pw_grid)
536 : END IF
537 : ! setup real space grid for finite difference derivatives of dielectric constant function
538 10058 : IF (poisson_params%has_dielectric .AND. &
539 : ((poisson_params%dielectric_params%derivative_method == derivative_cd3) .OR. &
540 : (poisson_params%dielectric_params%derivative_method == derivative_cd5) .OR. &
541 : (poisson_params%dielectric_params%derivative_method == derivative_cd7))) THEN
542 :
543 70 : SELECT CASE (poisson_params%ps_implicit_params%boundary_condition)
544 : CASE (NEUMANN_BC, MIXED_BC)
545 : CALL setup_diel_rs_grid(pw_env%poisson_env%diel_rs_grid, &
546 : poisson_params%dielectric_params%derivative_method, input, &
547 20 : pw_env%poisson_env%dct_pw_grid)
548 : CASE (PERIODIC_BC, MIXED_PERIODIC_BC)
549 : CALL setup_diel_rs_grid(pw_env%poisson_env%diel_rs_grid, &
550 : poisson_params%dielectric_params%derivative_method, input, &
551 50 : pw_env%poisson_env%pw_pools(pw_env%poisson_env%pw_level)%pool%pw_grid)
552 : END SELECT
553 :
554 : END IF
555 :
556 : !
557 : !
558 : ! determine the maximum radii for mapped gaussians, needed to
559 : ! set up distributed rs grids
560 : !
561 : !
562 :
563 30174 : ALLOCATE (radius(ngrid_level))
564 :
565 10058 : CALL compute_max_radius(radius, pw_env, qs_env)
566 :
567 : !
568 : !
569 : ! set up the rs_grids and the cubes, requires 'radius' to be set up correctly
570 : !
571 : !
572 60338 : ALLOCATE (rs_descs(ngrid_level))
573 :
574 211208 : ALLOCATE (rs_grids(ngrid_level))
575 :
576 331904 : ALLOCATE (pw_env%cube_info(ngrid_level))
577 10058 : higher_grid_layout = [-1, -1, -1]
578 :
579 40222 : DO igrid_level = 1, ngrid_level
580 30164 : pw_grid => pw_pools(igrid_level)%pool%pw_grid
581 :
582 30164 : CALL init_d3_poly_module() ! a fairly arbitrary but sufficient spot to do this
583 : CALL init_cube_info(pw_env%cube_info(igrid_level), &
584 : pw_grid%dr(:), pw_grid%dh(:, :), pw_grid%dh_inv(:, :), pw_grid%orthorhombic, &
585 30164 : radius(igrid_level))
586 :
587 30164 : rs_grid_section => section_vals_get_subs_vals(input, "DFT%MGRID%RS_GRID")
588 :
589 : CALL init_input_type(input_settings, nsmax=2*MAX(1, return_cube_max_iradius(pw_env%cube_info(igrid_level))) + 1, &
590 : rs_grid_section=rs_grid_section, ilevel=igrid_level, &
591 30164 : higher_grid_layout=higher_grid_layout)
592 :
593 30164 : NULLIFY (rs_descs(igrid_level)%rs_desc)
594 30164 : CALL rs_grid_create_descriptor(rs_descs(igrid_level)%rs_desc, pw_grid, input_settings)
595 :
596 30242 : IF (rs_descs(igrid_level)%rs_desc%distributed) higher_grid_layout = rs_descs(igrid_level)%rs_desc%group_dim
597 :
598 40222 : CALL rs_grid_create(rs_grids(igrid_level), rs_descs(igrid_level)%rs_desc)
599 : END DO
600 10058 : pw_env%rs_descs => rs_descs
601 10058 : pw_env%rs_grids => rs_grids
602 :
603 10058 : DEALLOCATE (radius)
604 :
605 : ! Print grid information
606 :
607 10058 : IF (do_io) THEN
608 8958 : logger => cp_get_default_logger()
609 8958 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.Log')
610 : END IF
611 10058 : IF (iounit > 0) THEN
612 3592 : SELECT CASE (poisson_params%solver)
613 : CASE (pw_poisson_periodic)
614 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,A30)") &
615 1566 : "POISSON| Solver", "PERIODIC"
616 : CASE (pw_poisson_analytic)
617 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,A30)") &
618 16 : "POISSON| Solver", "ANALYTIC"
619 : CASE (pw_poisson_mt)
620 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,A30)") &
621 261 : "POISSON| Solver", ADJUSTR("Martyna-Tuckerman (MT)")
622 : WRITE (UNIT=iounit, FMT="(T2,A,T71,F10.3,/,T2,A,T71,F10.1)") &
623 261 : "POISSON| MT| Alpha", poisson_params%mt_alpha, &
624 522 : "POISSON| MT| Relative cutoff", poisson_params%mt_rel_cutoff
625 : CASE (pw_poisson_multipole)
626 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,A30)") &
627 9 : "POISSON| Solver", "MULTIPOLE (Bloechl)"
628 : CASE (pw_poisson_wavelet)
629 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,A30)") &
630 147 : "POISSON| Solver", "WAVELET"
631 : WRITE (UNIT=iounit, FMT="(T2,A,T71,I10)") &
632 147 : "POISSON| Wavelet| Scaling function", poisson_params%wavelet_scf_type
633 296 : SELECT CASE (poisson_params%wavelet_method)
634 : CASE (WAVELET0D)
635 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") &
636 122 : "POISSON| Periodicity", "NONE"
637 : CASE (WAVELET2D)
638 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") &
639 1 : "POISSON| Periodicity", "XZ"
640 : CASE (WAVELET3D)
641 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") &
642 24 : "POISSON| Periodicity", "XYZ"
643 : CASE DEFAULT
644 147 : CPABORT("Invalid periodicity for wavelet solver selected")
645 : END SELECT
646 : CASE (pw_poisson_implicit)
647 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,A30)") &
648 27 : "POISSON| Solver", "IMPLICIT (GENERALIZED)"
649 27 : boundary_condition = poisson_params%ps_implicit_params%boundary_condition
650 5 : SELECT CASE (boundary_condition)
651 : CASE (PERIODIC_BC)
652 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") &
653 5 : "POISSON| Boundary Condition", "PERIODIC"
654 : CASE (NEUMANN_BC, MIXED_BC)
655 11 : IF (boundary_condition == NEUMANN_BC) THEN
656 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") &
657 3 : "POISSON| Boundary Condition", "NEUMANN"
658 : ELSE
659 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") &
660 8 : "POISSON| Boundary Condition", "MIXED"
661 : END IF
662 30 : SELECT CASE (poisson_params%ps_implicit_params%neumann_directions)
663 : CASE (neumannXYZ)
664 8 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| homogeneous Neumann directions", "X, Y, Z"
665 8 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| periodic directions", "NONE"
666 : CASE (neumannXY)
667 0 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| homogeneous Neumann directions", "X, Y"
668 0 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| periodic directions", "Z"
669 : CASE (neumannXZ)
670 0 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| homogeneous Neumann directions", "X, Z"
671 0 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| periodic directions", "Y"
672 : CASE (neumannYZ)
673 1 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| homogeneous Neumann directions", "Y, Z"
674 1 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| periodic directions", "X"
675 : CASE (neumannX)
676 1 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| homogeneous Neumann directions", "X"
677 1 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| periodic directions", "Y, Z"
678 : CASE (neumannY)
679 0 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| homogeneous Neumann directions", "Y"
680 0 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| periodic directions", "X, Z"
681 : CASE (neumannZ)
682 1 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| homogeneous Neumann directions", "Z"
683 1 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") "POISSON| periodic directions", "X, Y"
684 : CASE DEFAULT
685 11 : CPABORT("Invalid combination of Neumann and periodic conditions.")
686 : END SELECT
687 : CASE (MIXED_PERIODIC_BC)
688 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") &
689 11 : "POISSON| Boundary Condition", "PERIODIC & DIRICHLET"
690 : CASE DEFAULT
691 27 : CPABORT("Invalid boundary conditions for the implicit (generalized) poisson solver.")
692 : END SELECT
693 : WRITE (UNIT=iounit, FMT="(T2,A,T71,I10)") &
694 27 : "POISSON| Maximum number of iterations", poisson_params%ps_implicit_params%max_iter
695 : WRITE (UNIT=iounit, FMT="(T2,A,T51,ES30.2)") &
696 27 : "POISSON| Convergence threshold", poisson_params%ps_implicit_params%tol
697 27 : IF (poisson_params%dielectric_params%dielec_functiontype == rho_dependent) THEN
698 : WRITE (UNIT=iounit, FMT="(T2,A,T51,F30.2)") &
699 25 : "POISSON| Dielectric Constant", poisson_params%dielectric_params%eps0
700 : ELSE
701 : WRITE (UNIT=iounit, FMT="(T2,A,T31,F9.2)", ADVANCE='NO') &
702 2 : "POISSON| Dielectric Constants", poisson_params%dielectric_params%eps0
703 3 : DO i = 1, poisson_params%dielectric_params%n_aa_cuboidal
704 : WRITE (UNIT=iounit, FMT="(F9.2)", ADVANCE='NO') &
705 3 : poisson_params%dielectric_params%aa_cuboidal_eps(i)
706 : END DO
707 4 : DO i = 1, poisson_params%dielectric_params%n_xaa_annular
708 : WRITE (UNIT=iounit, FMT="(F9.2)", ADVANCE='NO') &
709 4 : poisson_params%dielectric_params%xaa_annular_eps(i)
710 : END DO
711 2 : WRITE (UNIT=iounit, FMT='(A1,/)')
712 : END IF
713 : WRITE (UNIT=iounit, FMT="(T2,A,T51,ES30.2)") &
714 27 : "POISSON| Relaxation parameter", poisson_params%ps_implicit_params%omega
715 : CASE (pw_poisson_none)
716 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,A30)") &
717 0 : "POISSON| Solver", "NONE"
718 : CASE default
719 2026 : CPABORT("Invalid Poisson solver selected")
720 : END SELECT
721 2026 : IF ((poisson_params%solver /= pw_poisson_wavelet) .AND. &
722 : (poisson_params%solver /= pw_poisson_implicit)) THEN
723 7408 : IF (SUM(poisson_params%periodic(1:3)) == 0) THEN
724 : WRITE (UNIT=iounit, FMT="(T2,A,T77,A4)") &
725 278 : "POISSON| Periodicity", "NONE"
726 : ELSE
727 1574 : string = ""
728 1574 : IF (poisson_params%periodic(1) == 1) string = TRIM(string)//"X"
729 1574 : IF (poisson_params%periodic(2) == 1) string = TRIM(string)//"Y"
730 1574 : IF (poisson_params%periodic(3) == 1) string = TRIM(string)//"Z"
731 : WRITE (UNIT=iounit, FMT="(T2,A,T78,A3)") &
732 1574 : "POISSON| Periodicity", ADJUSTR(string)
733 : END IF
734 : END IF
735 : END IF
736 :
737 : IF ((dft_control%qs_control%method_id == do_method_gpw) .OR. &
738 : (dft_control%qs_control%method_id == do_method_gapw) .OR. &
739 : (dft_control%qs_control%method_id == do_method_gapw_xc) .OR. &
740 : (dft_control%qs_control%method_id == do_method_lrigpw) .OR. &
741 10058 : (dft_control%qs_control%method_id == do_method_rigpw) .OR. &
742 : (dft_control%qs_control%method_id == do_method_ofgpw)) THEN
743 6762 : IF ((poisson_params%solver /= pw_poisson_wavelet) .AND. &
744 : (poisson_params%solver /= pw_poisson_implicit)) THEN
745 21558 : IF (ANY(my_cell%perd(1:3) /= poisson_params%periodic(1:3))) THEN
746 : CALL cp_warn(__LOCATION__, &
747 636 : "The selected periodicities in the sections &CELL and &POISSON do not match")
748 : END IF
749 : END IF
750 : END IF
751 :
752 10058 : IF (do_io) THEN
753 : should_output = BTEST(cp_print_key_should_output(logger%iter_info, &
754 8958 : print_section, ''), cp_p_file)
755 8958 : IF (should_output) THEN
756 15928 : DO igrid_level = 1, ngrid_level
757 15928 : CALL rs_grid_print(rs_grids(igrid_level), iounit)
758 : END DO
759 : END IF
760 8958 : CALL cp_print_key_finished_output(iounit, logger, print_section, "")
761 : END IF
762 :
763 10058 : CALL timestop(handle)
764 :
765 110638 : END SUBROUTINE pw_env_rebuild
766 :
767 : ! **************************************************************************************************
768 : !> \brief computes the maximum radius
769 : !> \param radius ...
770 : !> \param pw_env ...
771 : !> \param qs_env ...
772 : !> \par History
773 : !> 10.2010 refactored [Joost VandeVondele]
774 : !> 01.2020 igrid_zet0_s initialisation code is reused in rho0_s_grid_create() [Sergey Chulkov]
775 : !> \author Joost VandeVondele
776 : ! **************************************************************************************************
777 10058 : SUBROUTINE compute_max_radius(radius, pw_env, qs_env)
778 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: radius
779 : TYPE(pw_env_type), POINTER :: pw_env
780 : TYPE(qs_environment_type), POINTER :: qs_env
781 :
782 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_max_radius'
783 : CHARACTER(LEN=8), DIMENSION(10), PARAMETER :: sbas = ["ORB ", "AUX ", "RI_AUX ", &
784 : "MAO ", "HARRIS ", "RI_HXC ", "RI_K ", "LRI_AUX ", "RHOIN ", "NUC "]
785 : CHARACTER(LEN=8), DIMENSION(5), PARAMETER :: &
786 : pbas = ["ORB ", "AUX_FIT ", "MAO ", "HARRIS ", "NUC "]
787 : REAL(KIND=dp), PARAMETER :: safety_factor = 1.015_dp
788 :
789 : INTEGER :: handle, ibasis_set_type, igrid_level, igrid_zet0_s, ikind, ipgf, iset, ishell, &
790 : jkind, jpgf, jset, jshell, la, lb, lgrid_level, ngrid_level, nkind, nseta, nsetb
791 10058 : INTEGER, DIMENSION(:), POINTER :: npgfa, npgfb, nshella, nshellb
792 10058 : INTEGER, DIMENSION(:, :), POINTER :: lshella, lshellb
793 : REAL(KIND=dp) :: alpha, core_charge, eps_gvg, eps_rho, &
794 : max_rpgf0_s, maxradius, zet0_h, zetp
795 10058 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb
796 : TYPE(cneo_potential_type), POINTER :: cneo_potential
797 : TYPE(dft_control_type), POINTER :: dft_control
798 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
799 10058 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
800 : TYPE(qs_kind_type), POINTER :: qs_kind
801 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
802 :
803 : ! a very small safety factor might be needed for roundoff issues
804 : ! e.g. radius being computed here as 12.998 (13) and 13.002 (14) during the collocation
805 : ! the latter can happen due to the lower precision in the computation of the radius in collocate
806 : ! parallel cost of rs_pw_transfer goes as safety_factor**3 so it is worthwhile keeping it tight
807 :
808 10058 : CALL timeset(routineN, handle)
809 10058 : NULLIFY (dft_control, qs_kind_set, rho0_mpole)
810 :
811 10058 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
812 :
813 10058 : eps_rho = dft_control%qs_control%eps_rho_rspace
814 10058 : eps_gvg = dft_control%qs_control%eps_gvg_rspace
815 :
816 10058 : IF (dft_control%qs_control%gapw) THEN
817 934 : CALL get_qs_env(qs_env=qs_env, rho0_mpole=rho0_mpole)
818 934 : CPASSERT(ASSOCIATED(rho0_mpole))
819 :
820 934 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, max_rpgf0_s=max_rpgf0_s, zet0_h=zet0_h)
821 934 : igrid_zet0_s = gaussian_gridlevel(pw_env%gridlevel_info, 2.0_dp*zet0_h)
822 934 : rho0_mpole%igrid_zet0_s = igrid_zet0_s
823 : END IF
824 :
825 10058 : ngrid_level = SIZE(radius)
826 10058 : nkind = SIZE(qs_kind_set)
827 :
828 : ! try to predict the maximum radius of the gaussians to be mapped on the grid
829 : ! up to now, it is not yet very good
830 40222 : radius = 0.0_dp
831 40222 : DO igrid_level = 1, ngrid_level
832 :
833 30164 : maxradius = 0.0_dp
834 : ! Take into account the radius of the soft compensation charge rho0_soft1
835 30164 : IF (dft_control%qs_control%gapw) THEN
836 3464 : IF (igrid_zet0_s == igrid_level) maxradius = MAX(maxradius, max_rpgf0_s)
837 : END IF
838 :
839 : ! this is to be sure that the core charge is mapped ok
840 : ! right now, the core is mapped on the auxiliary basis,
841 : ! this should, at a give point be changed
842 : ! so that also for the core a multigrid is used
843 85894 : DO ikind = 1, nkind
844 55730 : NULLIFY (cneo_potential)
845 55730 : CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
846 55730 : IF (ASSOCIATED(cneo_potential)) CYCLE
847 : CALL get_qs_kind(qs_kind_set(ikind), &
848 55716 : alpha_core_charge=alpha, ccore_charge=core_charge)
849 85880 : IF (alpha > 0.0_dp .AND. core_charge /= 0.0_dp) THEN
850 55114 : maxradius = MAX(maxradius, exp_radius(0, alpha, eps_rho, core_charge, rlow=maxradius))
851 : ! forces
852 55114 : maxradius = MAX(maxradius, exp_radius(1, alpha, eps_rho, core_charge, rlow=maxradius))
853 : END IF
854 : END DO
855 :
856 : ! loop over basis sets that are used in grid collocation directly (no product)
857 : ! e.g. for calculating a wavefunction or a RI density
858 331804 : DO ibasis_set_type = 1, SIZE(sbas)
859 889104 : DO ikind = 1, nkind
860 557300 : qs_kind => qs_kind_set(ikind)
861 557300 : NULLIFY (orb_basis_set)
862 : CALL get_qs_kind(qs_kind=qs_kind, &
863 557300 : basis_set=orb_basis_set, basis_type=sbas(ibasis_set_type))
864 557300 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
865 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
866 68284 : npgf=npgfa, nset=nseta, zet=zeta, l=lshella, nshell=nshella)
867 624002 : DO iset = 1, nseta
868 1305876 : DO ipgf = 1, npgfa(iset)
869 1616072 : DO ishell = 1, nshella(iset)
870 867496 : zetp = zeta(ipgf, iset)
871 867496 : la = lshella(ishell, iset)
872 867496 : lgrid_level = gaussian_gridlevel(pw_env%gridlevel_info, zetp)
873 1361994 : IF (lgrid_level == igrid_level) THEN
874 289061 : maxradius = MAX(maxradius, exp_radius(la, zetp, eps_rho, 1.0_dp, rlow=maxradius))
875 : END IF
876 : END DO
877 : END DO
878 : END DO
879 : END DO
880 : END DO
881 : ! loop over basis sets that are used in product function grid collocation
882 180984 : DO ibasis_set_type = 1, SIZE(pbas)
883 459634 : DO ikind = 1, nkind
884 278650 : qs_kind => qs_kind_set(ikind)
885 278650 : NULLIFY (orb_basis_set)
886 : CALL get_qs_kind(qs_kind=qs_kind, &
887 278650 : basis_set=orb_basis_set, basis_type=pbas(ibasis_set_type))
888 278650 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
889 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
890 60582 : npgf=npgfa, nset=nseta, zet=zeta, l=lshella, nshell=nshella)
891 :
892 340202 : DO jkind = 1, nkind
893 128800 : qs_kind => qs_kind_set(jkind)
894 : CALL get_qs_kind(qs_kind=qs_kind, &
895 128800 : basis_set=orb_basis_set, basis_type=pbas(ibasis_set_type))
896 128800 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
897 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
898 128786 : npgf=npgfb, nset=nsetb, zet=zetb, l=lshellb, nshell=nshellb)
899 682088 : DO iset = 1, nseta
900 1228560 : DO ipgf = 1, npgfa(iset)
901 2661042 : DO ishell = 1, nshella(iset)
902 1561282 : la = lshella(ishell, iset)
903 5630892 : DO jset = 1, nsetb
904 16094788 : DO jpgf = 1, npgfb(jset)
905 44019100 : DO jshell = 1, nshellb(jset)
906 29485594 : zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
907 29485594 : lb = lshellb(jshell, jset) + la
908 29485594 : lgrid_level = gaussian_gridlevel(pw_env%gridlevel_info, zetp)
909 40774598 : IF (lgrid_level == igrid_level) THEN
910 : ! density (scale is at most 2)
911 12319659 : maxradius = MAX(maxradius, exp_radius(lb, zetp, eps_rho, 2.0_dp, rlow=maxradius))
912 : ! tau, properties?
913 12319659 : maxradius = MAX(maxradius, exp_radius(lb + 1, zetp, eps_rho, 2.0_dp, rlow=maxradius))
914 : ! potential
915 12319659 : maxradius = MAX(maxradius, exp_radius(lb, zetp, eps_gvg, 2.0_dp, rlow=maxradius))
916 : ! forces
917 12319659 : maxradius = MAX(maxradius, exp_radius(lb + 1, zetp, eps_gvg, 2.0_dp, rlow=maxradius))
918 : END IF
919 : END DO
920 : END DO
921 : END DO
922 : END DO
923 : END DO
924 : END DO
925 : END DO
926 : END DO
927 : END DO
928 :
929 : ! this is a bit of hack, but takes into account numerics and rounding
930 30164 : maxradius = maxradius*safety_factor
931 40222 : radius(igrid_level) = maxradius
932 : END DO
933 :
934 10058 : CALL timestop(handle)
935 :
936 10058 : END SUBROUTINE compute_max_radius
937 :
938 : ! **************************************************************************************************
939 : !> \brief Initialize the super-reference grid for Tuckerman or xc
940 : !> \param super_ref_pw_grid ...
941 : !> \param mt_super_ref_pw_grid ...
942 : !> \param xc_super_ref_pw_grid ...
943 : !> \param cutilev ...
944 : !> \param grid_span ...
945 : !> \param spherical ...
946 : !> \param cell_ref ...
947 : !> \param para_env ...
948 : !> \param poisson_section ...
949 : !> \param my_ncommensurate ...
950 : !> \param uf_grid ...
951 : !> \param print_section ...
952 : !> \author 03-2005 Teodoro Laino [teo]
953 : !> \note
954 : !> move somewere else?
955 : ! **************************************************************************************************
956 20116 : SUBROUTINE setup_super_ref_grid(super_ref_pw_grid, mt_super_ref_pw_grid, &
957 : xc_super_ref_pw_grid, cutilev, grid_span, spherical, &
958 : cell_ref, para_env, poisson_section, my_ncommensurate, uf_grid, print_section)
959 : TYPE(pw_grid_type), POINTER :: super_ref_pw_grid, mt_super_ref_pw_grid, &
960 : xc_super_ref_pw_grid
961 : REAL(KIND=dp), INTENT(IN) :: cutilev
962 : INTEGER, INTENT(IN) :: grid_span
963 : LOGICAL, INTENT(in) :: spherical
964 : TYPE(cell_type), POINTER :: cell_ref
965 : TYPE(mp_para_env_type), POINTER :: para_env
966 : TYPE(section_vals_type), POINTER :: poisson_section
967 : INTEGER, INTENT(IN) :: my_ncommensurate
968 : LOGICAL, INTENT(in) :: uf_grid
969 : TYPE(section_vals_type), POINTER :: print_section
970 :
971 : INTEGER :: iounit, my_val, nn(3), no(3)
972 : LOGICAL :: mt_s_grid
973 : REAL(KIND=dp) :: mt_rel_cutoff, my_cutilev
974 : TYPE(cp_logger_type), POINTER :: logger
975 :
976 10058 : CPASSERT(.NOT. ASSOCIATED(mt_super_ref_pw_grid))
977 10058 : CPASSERT(.NOT. ASSOCIATED(xc_super_ref_pw_grid))
978 10058 : CPASSERT(.NOT. ASSOCIATED(super_ref_pw_grid))
979 10058 : CALL section_vals_val_get(poisson_section, "POISSON_SOLVER", i_val=my_val)
980 : !
981 : ! Check if grids will be the same... In this case we don't use a super-reference grid
982 : !
983 10058 : mt_s_grid = .FALSE.
984 10058 : IF (my_val == pw_poisson_mt) THEN
985 : CALL section_vals_val_get(poisson_section, "MT%REL_CUTOFF", &
986 1224 : r_val=mt_rel_cutoff)
987 1224 : IF (mt_rel_cutoff > 1._dp) mt_s_grid = .TRUE.
988 : END IF
989 :
990 10058 : logger => cp_get_default_logger()
991 : iounit = cp_print_key_unit_nr(logger, print_section, "", &
992 10058 : extension=".Log")
993 :
994 10058 : IF (uf_grid) THEN
995 : CALL pw_grid_create(xc_super_ref_pw_grid, para_env, cell_ref%hmat, grid_span=grid_span, &
996 : cutoff=4._dp*cutilev, spherical=spherical, odd=.FALSE., fft_usage=.TRUE., &
997 : ncommensurate=my_ncommensurate, icommensurate=1, &
998 : blocked=do_pw_grid_blocked_false, rs_dims=[para_env%num_pe, 1], &
999 12 : iounit=iounit)
1000 4 : super_ref_pw_grid => xc_super_ref_pw_grid
1001 : END IF
1002 10058 : IF (mt_s_grid) THEN
1003 1218 : IF (ASSOCIATED(xc_super_ref_pw_grid)) THEN
1004 0 : CPABORT("special grid for mt and fine xc grid not compatible")
1005 : ELSE
1006 1218 : my_cutilev = cutilev*mt_rel_cutoff
1007 :
1008 : no = pw_grid_init_setup(cell_ref%hmat, cutoff=cutilev, spherical=spherical, &
1009 1218 : odd=.FALSE., fft_usage=.TRUE., ncommensurate=0, icommensurate=1)
1010 : nn = pw_grid_init_setup(cell_ref%hmat, cutoff=my_cutilev, spherical=spherical, &
1011 1218 : odd=.FALSE., fft_usage=.TRUE., ncommensurate=0, icommensurate=1)
1012 :
1013 : ! bug appears for nn==no, also in old versions
1014 4872 : CPASSERT(ALL(nn > no))
1015 : CALL pw_grid_create(mt_super_ref_pw_grid, para_env, cell_ref%hmat, &
1016 : cutoff=my_cutilev, spherical=spherical, fft_usage=.TRUE., &
1017 : blocked=do_pw_grid_blocked_false, rs_dims=[para_env%num_pe, 1], &
1018 3654 : iounit=iounit)
1019 1218 : super_ref_pw_grid => mt_super_ref_pw_grid
1020 : END IF
1021 : END IF
1022 : CALL cp_print_key_finished_output(iounit, logger, print_section, &
1023 10058 : "")
1024 10058 : END SUBROUTINE setup_super_ref_grid
1025 :
1026 : ! **************************************************************************************************
1027 : !> \brief sets up a real-space grid for finite difference derivative of dielectric
1028 : !> constant function
1029 : !> \param diel_rs_grid real space grid to be created
1030 : !> \param method preferred finite difference derivative method
1031 : !> \param input input file
1032 : !> \param pw_grid plane-wave grid
1033 : !> \par History
1034 : !> 12.2014 created [Hossein Bani-Hashemian]
1035 : !> \author Mohammad Hossein Bani-Hashemian
1036 : ! **************************************************************************************************
1037 50 : SUBROUTINE setup_diel_rs_grid(diel_rs_grid, method, input, pw_grid)
1038 :
1039 : TYPE(realspace_grid_type), POINTER :: diel_rs_grid
1040 : INTEGER, INTENT(IN) :: method
1041 : TYPE(section_vals_type), POINTER :: input
1042 : TYPE(pw_grid_type), POINTER :: pw_grid
1043 :
1044 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_diel_rs_grid'
1045 :
1046 : INTEGER :: border_points, handle
1047 : TYPE(realspace_grid_desc_type), POINTER :: rs_desc
1048 : TYPE(realspace_grid_input_type) :: input_settings
1049 : TYPE(section_vals_type), POINTER :: rs_grid_section
1050 :
1051 50 : CALL timeset(routineN, handle)
1052 :
1053 50 : NULLIFY (rs_desc)
1054 50 : rs_grid_section => section_vals_get_subs_vals(input, "DFT%MGRID%RS_GRID")
1055 74 : SELECT CASE (method)
1056 : CASE (derivative_cd3)
1057 24 : border_points = 1
1058 : CASE (derivative_cd5)
1059 14 : border_points = 2
1060 : CASE (derivative_cd7)
1061 50 : border_points = 3
1062 : END SELECT
1063 : CALL init_input_type(input_settings, 2*border_points + 1, rs_grid_section, &
1064 50 : 1, [-1, -1, -1])
1065 : CALL rs_grid_create_descriptor(rs_desc, pw_grid, input_settings, &
1066 50 : border_points=border_points)
1067 1050 : ALLOCATE (diel_rs_grid)
1068 50 : CALL rs_grid_create(diel_rs_grid, rs_desc)
1069 50 : CALL rs_grid_release_descriptor(rs_desc)
1070 :
1071 50 : CALL timestop(handle)
1072 :
1073 200 : END SUBROUTINE setup_diel_rs_grid
1074 :
1075 : END MODULE pw_env_methods
|
