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