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