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