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