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