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 9570 : 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 9570 : CALL timeset(routineN, handle)
140 :
141 9570 : CPASSERT(.NOT. ASSOCIATED(pw_env))
142 133980 : 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 9570 : pw_env%auxbas_grid = -1
148 9570 : pw_env%ref_count = 1
149 :
150 9570 : CALL timestop(handle)
151 :
152 9570 : 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 12418 : 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 12418 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: radius
180 12418 : 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 12418 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
189 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
190 12418 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
191 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
192 12418 : POINTER :: rs_descs
193 : TYPE(realspace_grid_input_type) :: input_settings
194 12418 : 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 12418 : CALL timeset(routineN, handle)
206 :
207 : !
208 : !
209 : ! Part one, deallocate old data if needed
210 : !
211 : !
212 12418 : NULLIFY (cutoff, cell, pw_grid, old_pw_grid, dft_control, qs_kind_set, &
213 12418 : pw_pools, rs_descs, para_env, cell_ref, vdw_ref_grid, &
214 12418 : mt_super_ref_grid, input, poisson_section, xc_super_ref_grid, &
215 12418 : 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 12418 : dispersion_env=dispersion_env)
225 :
226 12418 : CPASSERT(ASSOCIATED(pw_env))
227 12418 : CPASSERT(pw_env%ref_count > 0)
228 12418 : CALL pw_pool_release(pw_env%vdw_pw_pool)
229 12418 : CALL pw_pool_release(pw_env%xc_pw_pool)
230 12418 : CALL pw_pools_dealloc(pw_env%pw_pools)
231 12418 : IF (ASSOCIATED(pw_env%rs_descs)) THEN
232 7550 : DO i = 1, SIZE(pw_env%rs_descs)
233 7550 : CALL rs_grid_release_descriptor(pw_env%rs_descs(i)%rs_desc)
234 : END DO
235 2872 : DEALLOCATE (pw_env%rs_descs)
236 : END IF
237 12418 : IF (ASSOCIATED(pw_env%rs_grids)) THEN
238 7550 : DO i = 1, SIZE(pw_env%rs_grids)
239 7550 : CALL rs_grid_release(pw_env%rs_grids(i))
240 : END DO
241 2872 : DEALLOCATE (pw_env%rs_grids)
242 : END IF
243 12418 : IF (ASSOCIATED(pw_env%gridlevel_info)) THEN
244 2872 : CALL destroy_gaussian_gridlevel(pw_env%gridlevel_info)
245 : ELSE
246 9546 : ALLOCATE (pw_env%gridlevel_info)
247 : END IF
248 :
249 12418 : IF (ASSOCIATED(pw_env%cube_info)) THEN
250 7550 : DO igrid_level = 1, SIZE(pw_env%cube_info)
251 7550 : CALL destroy_cube_info(pw_env%cube_info(igrid_level))
252 : END DO
253 2872 : DEALLOCATE (pw_env%cube_info)
254 : END IF
255 12418 : 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 12418 : CALL init_fft_scratch_pool()
259 :
260 : !
261 : !
262 : ! Part two, setup the pw_grids
263 : !
264 : !
265 :
266 12418 : do_io = .TRUE.
267 12418 : IF (PRESENT(external_para_env)) THEN
268 1116 : para_env => external_para_env
269 : CPASSERT(ASSOCIATED(para_env))
270 1116 : do_io = .FALSE. !multiple MPI subgroups mess-up the output file
271 : END IF
272 : ! interpolation section
273 12418 : pw_env%interp_section => section_vals_get_subs_vals(input, "DFT%MGRID%INTERPOLATOR")
274 :
275 12418 : CALL get_qs_env(qs_env, use_ref_cell=use_ref_cell)
276 12418 : IF (use_ref_cell) THEN
277 64 : my_cell => cell_ref
278 : ELSE
279 12354 : my_cell => cell
280 : END IF
281 12418 : rel_cutoff = dft_control%qs_control%relative_cutoff
282 12418 : cutoff => dft_control%qs_control%e_cutoff
283 12418 : CALL section_vals_val_get(input, "DFT%XC%XC_GRID%USE_FINER_GRID", l_val=uf_grid)
284 12418 : CALL section_vals_val_get(input, "DFT%XC%XC_GRID%FINE_GRID_FACTOR", r_val=fine_grid_factor)
285 12418 : 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 12418 : 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 12418 : "PRINT%GRID_INFORMATION")
295 : CALL init_gaussian_gridlevel(pw_env%gridlevel_info, &
296 : ngrid_levels=ngrid_level, cutoff=cutoff, rel_cutoff=rel_cutoff, &
297 12418 : print_section=print_section)
298 : ! init pw_grids and pools
299 72686 : ALLOCATE (pw_pools(ngrid_level))
300 :
301 12418 : IF (dft_control%qs_control%commensurate_mgrids) THEN
302 274 : ncommensurate = ngrid_level
303 : ELSE
304 12144 : ncommensurate = 0
305 : END IF
306 : !
307 : ! If Tuckerman is present let's perform the set-up of the super-reference-grid
308 : !
309 12418 : cutilev = cutoff(1)
310 12418 : IF (dft_control%qs_control%pw_grid_opt%spherical) THEN
311 0 : grid_span = HALFSPACE
312 0 : spherical = .TRUE.
313 12418 : ELSE IF (dft_control%qs_control%pw_grid_opt%fullspace) THEN
314 12418 : grid_span = FULLSPACE
315 12418 : spherical = .FALSE.
316 : ELSE
317 0 : grid_span = HALFSPACE
318 0 : spherical = .FALSE.
319 : END IF
320 :
321 12418 : 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 12418 : print_section=print_section)
327 12418 : old_pw_grid => super_ref_grid
328 12418 : 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 12418 : IF (do_io) THEN
334 11302 : logger => cp_get_default_logger()
335 11302 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.Log')
336 : ELSE
337 1116 : iounit = 0
338 : END IF
339 :
340 12418 : IF (dft_control%qs_control%pw_grid_opt%spherical) THEN
341 0 : grid_span = HALFSPACE
342 0 : spherical = .TRUE.
343 0 : odd = .TRUE.
344 12418 : ELSE IF (dft_control%qs_control%pw_grid_opt%fullspace) THEN
345 12418 : grid_span = FULLSPACE
346 12418 : spherical = .FALSE.
347 12418 : 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 12418 : 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 12418 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
360 12418 : xc_deriv_method_id = section_get_ival(xc_section, "XC_GRID%XC_DERIV")
361 12418 : xc_smooth_method_id = section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO")
362 12418 : 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 12418 : 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 12418 : 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 12418 : subsection_name="PROPERTIES%LINRES")
384 12418 : CALL section_vals_get(linres_section, explicit=linres_present)
385 12418 : 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 12418 : subsection_name="DFT%PRINT%ELECTRIC_FIELD_GRADIENT")
391 12418 : CALL section_vals_get(efg_section, explicit=efg_present)
392 12418 : IF (efg_present) THEN
393 2 : smooth_required = smooth_required .OR. .TRUE.
394 : END IF
395 :
396 47850 : DO igrid_level = 1, ngrid_level
397 35432 : 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 35432 : blocked_id = blocked_id_input
402 :
403 106296 : 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 35432 : 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 35432 : IF (smooth_required) THEN
414 4326 : distribution_layout = [para_env%num_pe, 1]
415 : END IF
416 :
417 35432 : IF (igrid_level == 1) THEN
418 12418 : 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 1546 : iounit=iounit)
427 1546 : 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 10872 : iounit=iounit)
436 10872 : 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 23014 : iounit=iounit)
447 : END IF
448 :
449 : ! init pw_pools
450 35432 : NULLIFY (pw_pools(igrid_level)%pool)
451 35432 : CALL pw_pool_create(pw_pools(igrid_level)%pool, pw_grid=pw_grid)
452 :
453 47850 : CALL pw_grid_release(pw_grid)
454 :
455 : END DO
456 :
457 12418 : pw_env%pw_pools => pw_pools
458 :
459 : ! init auxbas_grid
460 47850 : DO i = 1, ngrid_level
461 47850 : IF (cutoff(i) == dft_control%qs_control%cutoff) pw_env%auxbas_grid = i
462 : END DO
463 :
464 : ! init xc_pool
465 12418 : 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 12370 : pw_env%xc_pw_pool => pw_pools(pw_env%auxbas_grid)%pool
471 12370 : CALL pw_env%xc_pw_pool%retain()
472 : END IF
473 :
474 : ! init vdw_pool
475 12418 : set_vdw_pool = .FALSE.
476 12418 : IF (ASSOCIATED(dispersion_env)) THEN
477 11614 : 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 12348 : pw_env%vdw_pw_pool => pw_pools(pw_env%auxbas_grid)%pool
497 12348 : CALL pw_env%vdw_pw_pool%retain()
498 : END IF
499 :
500 12418 : IF (do_io) CALL cp_print_key_finished_output(iounit, logger, print_section, '')
501 :
502 : ! complete init of the poisson_env
503 12418 : IF (.NOT. ASSOCIATED(pw_env%poisson_env)) THEN
504 152736 : ALLOCATE (pw_env%poisson_env)
505 9546 : CALL pw_env%poisson_env%create()
506 : END IF
507 : ! poisson_section => section_vals_get_subs_vals(input, "DFT%POISSON")
508 :
509 12418 : 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 12418 : dct_pw_grid=dct_pw_grid, use_level=pw_env%auxbas_grid)
513 12418 : CALL pw_grid_release(mt_super_ref_grid)
514 12418 : 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 12418 : 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 12418 : 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 12418 : 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 12418 : 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 37254 : ALLOCATE (radius(ngrid_level))
570 :
571 12418 : 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 72686 : ALLOCATE (rs_descs(ngrid_level))
579 :
580 258956 : ALLOCATE (rs_grids(ngrid_level))
581 :
582 407972 : ALLOCATE (pw_env%cube_info(ngrid_level))
583 12418 : higher_grid_layout = [-1, -1, -1]
584 :
585 47850 : DO igrid_level = 1, ngrid_level
586 35432 : pw_grid => pw_pools(igrid_level)%pool%pw_grid
587 :
588 35432 : 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 35432 : radius(igrid_level))
592 :
593 35432 : 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 35432 : higher_grid_layout=higher_grid_layout)
598 :
599 35432 : NULLIFY (rs_descs(igrid_level)%rs_desc)
600 35432 : CALL rs_grid_create_descriptor(rs_descs(igrid_level)%rs_desc, pw_grid, input_settings)
601 :
602 35510 : IF (rs_descs(igrid_level)%rs_desc%distributed) higher_grid_layout = rs_descs(igrid_level)%rs_desc%group_dim
603 :
604 47850 : CALL rs_grid_create(rs_grids(igrid_level), rs_descs(igrid_level)%rs_desc)
605 : END DO
606 12418 : pw_env%rs_descs => rs_descs
607 12418 : pw_env%rs_grids => rs_grids
608 :
609 12418 : DEALLOCATE (radius)
610 :
611 : ! Print grid information
612 :
613 12418 : IF (do_io) THEN
614 11302 : logger => cp_get_default_logger()
615 11302 : iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.Log')
616 : END IF
617 12418 : IF (iounit > 0) THEN
618 4794 : SELECT CASE (poisson_params%solver)
619 : CASE (pw_poisson_periodic)
620 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,A30)") &
621 2132 : "POISSON| Solver", "PERIODIC"
622 : CASE (pw_poisson_analytic)
623 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,A30)") &
624 19 : "POISSON| Solver", "ANALYTIC"
625 : CASE (pw_poisson_mt)
626 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,A30)") &
627 320 : "POISSON| Solver", ADJUSTR("Martyna-Tuckerman (MT)")
628 : WRITE (UNIT=iounit, FMT="(T2,A,T71,F10.3,/,T2,A,T71,F10.1)") &
629 320 : "POISSON| MT| Alpha", poisson_params%mt_alpha, &
630 640 : "POISSON| MT| Relative cutoff", poisson_params%mt_rel_cutoff
631 : CASE (pw_poisson_multipole)
632 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,A30)") &
633 10 : "POISSON| Solver", "MULTIPOLE (Bloechl)"
634 : CASE (pw_poisson_wavelet)
635 : WRITE (UNIT=iounit, FMT="(/,T2,A,T51,A30)") &
636 154 : "POISSON| Solver", "WAVELET"
637 : WRITE (UNIT=iounit, FMT="(T2,A,T71,I10)") &
638 154 : "POISSON| Wavelet| Scaling function", poisson_params%wavelet_scf_type
639 308 : SELECT CASE (poisson_params%wavelet_method)
640 : CASE (WAVELET0D)
641 : WRITE (UNIT=iounit, FMT="(T2,A,T51,A30)") &
642 127 : "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 154 : 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 2662 : CPABORT("Invalid Poisson solver selected")
730 : END SELECT
731 2662 : IF ((poisson_params%solver /= pw_poisson_wavelet) .AND. &
732 : (poisson_params%solver /= pw_poisson_implicit)) THEN
733 9924 : IF (SUM(poisson_params%periodic(1:3)) == 0) THEN
734 : WRITE (UNIT=iounit, FMT="(T2,A,T77,A4)") &
735 341 : "POISSON| Periodicity", "NONE"
736 : ELSE
737 2140 : string = ""
738 2140 : IF (poisson_params%periodic(1) == 1) string = TRIM(string)//"X"
739 2140 : IF (poisson_params%periodic(2) == 1) string = TRIM(string)//"Y"
740 2140 : IF (poisson_params%periodic(3) == 1) string = TRIM(string)//"Z"
741 : WRITE (UNIT=iounit, FMT="(T2,A,T78,A3)") &
742 2140 : "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 12418 : (dft_control%qs_control%method_id == do_method_rigpw) .OR. &
752 : (dft_control%qs_control%method_id == do_method_ofgpw)) THEN
753 7734 : IF (poisson_params%solver /= pw_poisson_implicit) THEN
754 28768 : IF (ANY(my_cell%perd(1:3) /= poisson_params%periodic(1:3))) THEN
755 : CALL cp_warn(__LOCATION__, &
756 666 : "The selected periodicities in the sections &CELL and &POISSON do not match")
757 : END IF
758 : END IF
759 : END IF
760 :
761 12418 : IF (do_io) THEN
762 : should_output = BTEST(cp_print_key_should_output(logger%iter_info, &
763 11302 : print_section, ''), cp_p_file)
764 11302 : IF (should_output) THEN
765 19250 : DO igrid_level = 1, ngrid_level
766 19250 : CALL rs_grid_print(rs_grids(igrid_level), iounit)
767 : END DO
768 : END IF
769 11302 : CALL cp_print_key_finished_output(iounit, logger, print_section, "")
770 : END IF
771 :
772 12418 : CALL timestop(handle)
773 :
774 136598 : 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 12418 : 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 12418 : INTEGER, DIMENSION(:), POINTER :: npgfa, npgfb, nshella, nshellb
801 12418 : 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 12418 : 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 12418 : 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 12418 : CALL timeset(routineN, handle)
818 12418 : NULLIFY (dft_control, qs_kind_set, rho0_mpole)
819 :
820 12418 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, dft_control=dft_control)
821 :
822 12418 : eps_rho = dft_control%qs_control%eps_rho_rspace
823 12418 : eps_gvg = dft_control%qs_control%eps_gvg_rspace
824 :
825 12418 : IF (dft_control%qs_control%gapw) THEN
826 1254 : CALL get_qs_env(qs_env=qs_env, rho0_mpole=rho0_mpole)
827 1254 : CPASSERT(ASSOCIATED(rho0_mpole))
828 :
829 1254 : CALL get_rho0_mpole(rho0_mpole=rho0_mpole, max_rpgf0_s=max_rpgf0_s, zet0_h=zet0_h)
830 1254 : igrid_zet0_s = gaussian_gridlevel(pw_env%gridlevel_info, 2.0_dp*zet0_h)
831 1254 : rho0_mpole%igrid_zet0_s = igrid_zet0_s
832 : END IF
833 :
834 12418 : ngrid_level = SIZE(radius)
835 12418 : 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 47850 : radius = 0.0_dp
840 47850 : DO igrid_level = 1, ngrid_level
841 :
842 35432 : maxradius = 0.0_dp
843 : ! Take into account the radius of the soft compensation charge rho0_soft1
844 35432 : IF (dft_control%qs_control%gapw) THEN
845 4740 : 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 99146 : DO ikind = 1, nkind
853 63714 : NULLIFY (cneo_potential)
854 63714 : CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
855 63714 : IF (ASSOCIATED(cneo_potential)) CYCLE
856 : CALL get_qs_kind(qs_kind_set(ikind), &
857 63700 : alpha_core_charge=alpha, ccore_charge=core_charge)
858 99132 : IF (alpha > 0.0_dp .AND. core_charge /= 0.0_dp) THEN
859 63098 : maxradius = MAX(maxradius, exp_radius(0, alpha, eps_rho, core_charge, rlow=maxradius))
860 : ! forces
861 63098 : 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 389752 : DO ibasis_set_type = 1, SIZE(sbas)
868 1026892 : DO ikind = 1, nkind
869 637140 : qs_kind => qs_kind_set(ikind)
870 637140 : NULLIFY (orb_basis_set)
871 : CALL get_qs_kind(qs_kind=qs_kind, &
872 637140 : basis_set=orb_basis_set, basis_type=sbas(ibasis_set_type))
873 637140 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
874 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
875 75470 : npgf=npgfa, nset=nseta, zet=zeta, l=lshella, nshell=nshella)
876 698524 : DO iset = 1, nseta
877 1446490 : DO ipgf = 1, npgfa(iset)
878 1769956 : DO ishell = 1, nshella(iset)
879 960606 : zetp = zeta(ipgf, iset)
880 960606 : la = lshella(ishell, iset)
881 960606 : lgrid_level = gaussian_gridlevel(pw_env%gridlevel_info, zetp)
882 1501222 : IF (lgrid_level == igrid_level) THEN
883 320881 : 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 212592 : DO ibasis_set_type = 1, SIZE(pbas)
892 531162 : DO ikind = 1, nkind
893 318570 : qs_kind => qs_kind_set(ikind)
894 318570 : NULLIFY (orb_basis_set)
895 : CALL get_qs_kind(qs_kind=qs_kind, &
896 318570 : basis_set=orb_basis_set, basis_type=pbas(ibasis_set_type))
897 318570 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
898 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
899 68222 : npgf=npgfa, nset=nseta, zet=zeta, l=lshella, nshell=nshella)
900 :
901 387046 : DO jkind = 1, nkind
902 141664 : qs_kind => qs_kind_set(jkind)
903 : CALL get_qs_kind(qs_kind=qs_kind, &
904 141664 : basis_set=orb_basis_set, basis_type=pbas(ibasis_set_type))
905 141664 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
906 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
907 141650 : npgf=npgfb, nset=nsetb, zet=zetb, l=lshellb, nshell=nshellb)
908 759378 : DO iset = 1, nseta
909 1338522 : DO ipgf = 1, npgfa(iset)
910 2912964 : DO ishell = 1, nshella(iset)
911 1716106 : la = lshella(ishell, iset)
912 6138148 : DO jset = 1, nsetb
913 17425688 : DO jpgf = 1, npgfb(jset)
914 47487164 : DO jshell = 1, nshellb(jset)
915 31777582 : zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
916 31777582 : lb = lshellb(jshell, jset) + la
917 31777582 : lgrid_level = gaussian_gridlevel(pw_env%gridlevel_info, zetp)
918 43962822 : IF (lgrid_level == igrid_level) THEN
919 : ! density (scale is at most 2)
920 12956973 : maxradius = MAX(maxradius, exp_radius(lb, zetp, eps_rho, 2.0_dp, rlow=maxradius))
921 : ! tau, properties?
922 12956973 : maxradius = MAX(maxradius, exp_radius(lb + 1, zetp, eps_rho, 2.0_dp, rlow=maxradius))
923 : ! potential
924 12956973 : maxradius = MAX(maxradius, exp_radius(lb, zetp, eps_gvg, 2.0_dp, rlow=maxradius))
925 : ! forces
926 12956973 : 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 35432 : maxradius = maxradius*safety_factor
940 47850 : radius(igrid_level) = maxradius
941 : END DO
942 :
943 12418 : CALL timestop(handle)
944 :
945 12418 : 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 24836 : 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 12418 : CPASSERT(.NOT. ASSOCIATED(mt_super_ref_pw_grid))
989 12418 : CPASSERT(.NOT. ASSOCIATED(xc_super_ref_pw_grid))
990 12418 : CPASSERT(.NOT. ASSOCIATED(super_ref_pw_grid))
991 12418 : 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 12418 : mt_s_grid = .FALSE.
996 12418 : IF (my_val == pw_poisson_mt) THEN
997 : CALL section_vals_val_get(poisson_section, "MT%REL_CUTOFF", &
998 1506 : r_val=mt_rel_cutoff)
999 1506 : IF (mt_rel_cutoff > 1._dp) mt_s_grid = .TRUE.
1000 : END IF
1001 :
1002 12418 : logger => cp_get_default_logger()
1003 : iounit = cp_print_key_unit_nr(logger, print_section, "", &
1004 12418 : extension=".Log")
1005 :
1006 12418 : 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 12418 : IF (mt_s_grid) THEN
1017 1500 : my_cutilev = cutilev*mt_rel_cutoff
1018 1500 : 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 1500 : odd=.FALSE., fft_usage=.TRUE., ncommensurate=0, icommensurate=1)
1022 : nn = pw_grid_init_setup(cell_ref%hmat, cutoff=my_cutilev, spherical=spherical, &
1023 1500 : odd=.FALSE., fft_usage=.TRUE., ncommensurate=0, icommensurate=1)
1024 :
1025 : ! bug appears for nn==no, also in old versions
1026 6000 : CPASSERT(ALL(nn > no))
1027 1500 : 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 4494 : iounit=iounit)
1035 1498 : 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 12418 : "")
1040 12418 : 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
|
