1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief pw_types 8!> \author CJM 9! ************************************************************************************************** 10MODULE ewald_pw_types 11 USE ao_util, ONLY: exp_radius 12 USE cell_types, ONLY: cell_type 13 USE cp_log_handling, ONLY: cp_get_default_logger,& 14 cp_logger_type 15 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 16 cp_print_key_unit_nr 17 USE cp_para_types, ONLY: cp_para_env_type 18 USE cp_realspace_grid_init, ONLY: init_input_type 19 USE dg_types, ONLY: dg_create,& 20 dg_release,& 21 dg_retain,& 22 dg_type 23 USE dgs, ONLY: dg_pme_grid_setup 24 USE ewald_environment_types, ONLY: ewald_env_get,& 25 ewald_environment_type 26 USE input_section_types, ONLY: section_vals_get_subs_vals,& 27 section_vals_type 28 USE kinds, ONLY: dp 29 USE mathconstants, ONLY: pi 30 USE message_passing, ONLY: mp_comm_self 31 USE pw_grid_types, ONLY: HALFSPACE,& 32 pw_grid_type 33 USE pw_grids, ONLY: pw_grid_create,& 34 pw_grid_release,& 35 pw_grid_setup 36 USE pw_poisson_methods, ONLY: pw_poisson_set 37 USE pw_poisson_read_input, ONLY: pw_poisson_read_parameters 38 USE pw_poisson_types, ONLY: & 39 do_ewald_ewald, do_ewald_none, do_ewald_pme, do_ewald_spme, pw_poisson_create, & 40 pw_poisson_parameter_type, pw_poisson_release, pw_poisson_retain, pw_poisson_type 41 USE pw_pool_types, ONLY: pw_pool_create,& 42 pw_pool_p_type,& 43 pw_pool_release,& 44 pw_pool_retain,& 45 pw_pool_type 46 USE realspace_grid_types, ONLY: & 47 realspace_grid_desc_type, realspace_grid_input_type, realspace_grid_type, rs_grid_create, & 48 rs_grid_create_descriptor, rs_grid_print, rs_grid_release, rs_grid_release_descriptor, & 49 rs_grid_retain_descriptor 50#include "./base/base_uses.f90" 51 52 IMPLICIT NONE 53 54 PRIVATE 55 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewald_pw_types' 56 INTEGER, PRIVATE, SAVE :: last_ewald_pw_id_nr = 0 57 PUBLIC :: ewald_pw_type, ewald_pw_release, & 58 ewald_pw_retain, ewald_pw_create, & 59 ewald_pw_get, ewald_pw_set 60 61! ************************************************************************************************** 62 TYPE ewald_pw_type 63 PRIVATE 64 INTEGER :: ref_count, id_nr 65 TYPE(pw_pool_type), POINTER :: pw_small_pool 66 TYPE(pw_pool_type), POINTER :: pw_big_pool 67 TYPE(realspace_grid_desc_type), POINTER :: rs_desc 68 TYPE(pw_poisson_type), POINTER :: poisson_env 69 TYPE(dg_type), POINTER :: dg 70 END TYPE ewald_pw_type 71 72CONTAINS 73 74! ************************************************************************************************** 75!> \brief retains the structure ewald_pw_type 76!> \param ewald_pw ... 77! ************************************************************************************************** 78 SUBROUTINE ewald_pw_retain(ewald_pw) 79 TYPE(ewald_pw_type), POINTER :: ewald_pw 80 81 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_pw_retain', & 82 routineP = moduleN//':'//routineN 83 84 CPASSERT(ASSOCIATED(ewald_pw)) 85 CPASSERT(ewald_pw%ref_count > 0) 86 ewald_pw%ref_count = ewald_pw%ref_count + 1 87 END SUBROUTINE ewald_pw_retain 88 89! ************************************************************************************************** 90!> \brief creates the structure ewald_pw_type 91!> \param ewald_pw ... 92!> \param ewald_env ... 93!> \param cell ... 94!> \param cell_ref ... 95!> \param print_section ... 96! ************************************************************************************************** 97 SUBROUTINE ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section) 98 TYPE(ewald_pw_type), POINTER :: ewald_pw 99 TYPE(ewald_environment_type), POINTER :: ewald_env 100 TYPE(cell_type), POINTER :: cell, cell_ref 101 TYPE(section_vals_type), POINTER :: print_section 102 103 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_pw_create', & 104 routineP = moduleN//':'//routineN 105 106 TYPE(dg_type), POINTER :: dg 107 108 NULLIFY (dg) 109 ALLOCATE (ewald_pw) 110 NULLIFY (ewald_pw%pw_big_pool) 111 NULLIFY (ewald_pw%pw_small_pool) 112 NULLIFY (ewald_pw%rs_desc) 113 NULLIFY (ewald_pw%poisson_env) 114 CALL dg_create(dg) 115 ewald_pw%dg => dg 116 ewald_pw%ref_count = 1 117 last_ewald_pw_id_nr = last_ewald_pw_id_nr + 1 118 ewald_pw%id_nr = last_ewald_pw_id_nr 119 CALL ewald_pw_init(ewald_pw, ewald_env, cell, cell_ref, print_section) 120 END SUBROUTINE ewald_pw_create 121 122!****f* ewald_pw_types/ewald_pw_release [1.0] * 123 124! ************************************************************************************************** 125!> \brief releases the memory used by the ewald_pw 126!> \param ewald_pw ... 127! ************************************************************************************************** 128 SUBROUTINE ewald_pw_release(ewald_pw) 129 TYPE(ewald_pw_type), POINTER :: ewald_pw 130 131 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_pw_release', & 132 routineP = moduleN//':'//routineN 133 134 INTEGER :: handle 135 136 CALL timeset(routineN, handle) 137 IF (ASSOCIATED(ewald_pw)) THEN 138 CPASSERT(ewald_pw%ref_count > 0) 139 ewald_pw%ref_count = ewald_pw%ref_count - 1 140 IF (ewald_pw%ref_count < 1) THEN 141 CALL pw_pool_release(ewald_pw%pw_small_pool) 142 CALL pw_pool_release(ewald_pw%pw_big_pool) 143 CALL rs_grid_release_descriptor(ewald_pw%rs_desc) 144 CALL pw_poisson_release(ewald_pw%poisson_env) 145 CALL dg_release(ewald_pw%dg) 146 DEALLOCATE (ewald_pw) 147 END IF 148 END IF 149 NULLIFY (ewald_pw) 150 CALL timestop(handle) 151 END SUBROUTINE ewald_pw_release 152 153!****** ewald_pw_types/ewald_pw_init [1.0] * 154 155! ************************************************************************************************** 156!> \brief ... 157!> \param ewald_pw ... 158!> \param ewald_env ... 159!> \param cell ... 160!> \param cell_ref ... 161!> \param print_section ... 162!> \par History 163!> JGH (12-Jan-2001): Added SPME part 164!> JGH (15-Mar-2001): Work newly distributed between initialize, setup, 165!> and force routine 166!> \author CJM 167! ************************************************************************************************** 168 SUBROUTINE ewald_pw_init(ewald_pw, ewald_env, cell, cell_ref, print_section) 169 TYPE(ewald_pw_type), POINTER :: ewald_pw 170 TYPE(ewald_environment_type), POINTER :: ewald_env 171 TYPE(cell_type), POINTER :: cell, cell_ref 172 TYPE(section_vals_type), POINTER :: print_section 173 174 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_pw_init', routineP = moduleN//':'//routineN 175 176 INTEGER :: bo(2, 3), ewald_type, gmax(3), handle, & 177 npts_s(3), ns_max, o_spline, & 178 output_unit 179 REAL(KIND=dp) :: alpha, alphasq, cutoff_radius, epsilon, & 180 norm 181 TYPE(cp_logger_type), POINTER :: logger 182 TYPE(cp_para_env_type), POINTER :: para_env 183 TYPE(pw_grid_type), POINTER :: pw_big_grid, pw_small_grid 184 TYPE(pw_poisson_parameter_type) :: poisson_params 185 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 186 TYPE(pw_pool_type), POINTER :: pw_pool 187 TYPE(realspace_grid_desc_type), POINTER :: rs_desc 188 TYPE(realspace_grid_input_type) :: input_settings 189 TYPE(realspace_grid_type), POINTER :: rs 190 TYPE(section_vals_type), POINTER :: poisson_section, rs_grid_section 191 192 CALL timeset(routineN, handle) 193 194 NULLIFY (pw_big_grid) 195 NULLIFY (pw_small_grid, poisson_section) 196 197 CPASSERT(ASSOCIATED(ewald_pw)) 198 CPASSERT(ASSOCIATED(ewald_env)) 199 CPASSERT(ASSOCIATED(cell)) 200 CPASSERT(ewald_pw%ref_count > 0) 201 CALL ewald_env_get(ewald_env=ewald_env, & 202 para_env=para_env, & 203 gmax=gmax, alpha=alpha, & 204 ns_max=ns_max, & 205 ewald_type=ewald_type, & 206 o_spline=o_spline, & 207 poisson_section=poisson_section, & 208 epsilon=epsilon) 209 210 rs_grid_section => section_vals_get_subs_vals(poisson_section, "EWALD%RS_GRID") 211 212 SELECT CASE (ewald_type) 213 CASE (do_ewald_ewald) 214 ! set up Classic EWALD sum 215 logger => cp_get_default_logger() 216 output_unit = cp_print_key_unit_nr(logger, print_section, "", extension=".Log") 217 CALL pw_grid_create(pw_big_grid, mp_comm_self) 218 219 IF (ANY(gmax == 2*(gmax/2))) THEN 220 CPABORT("gmax has to be odd.") 221 END IF 222 bo(1, :) = -gmax/2 223 bo(2, :) = +gmax/2 224 CALL pw_grid_setup(cell_ref%hmat, pw_big_grid, grid_span=HALFSPACE, bounds=bo, spherical=.TRUE., & 225 fft_usage=.FALSE., iounit=output_unit) 226 NULLIFY (pw_pool) 227 CALL pw_pool_create(pw_pool, pw_grid=pw_big_grid) 228 ewald_pw%pw_big_pool => pw_pool 229 CALL pw_pool_retain(ewald_pw%pw_big_pool) 230 CALL pw_pool_release(pw_pool) 231 CALL pw_grid_release(pw_big_grid) 232 CALL cp_print_key_finished_output(output_unit, logger, print_section, "") 233 234 CASE (do_ewald_pme) 235 ! set up Particle-Mesh EWALD sum 236 logger => cp_get_default_logger() 237 output_unit = cp_print_key_unit_nr(logger, print_section, "", extension=".Log") 238 IF (.NOT. ASSOCIATED(ewald_pw%poisson_env)) THEN 239 CALL pw_poisson_create(ewald_pw%poisson_env) 240 END IF 241 CALL pw_grid_create(pw_small_grid, mp_comm_self) 242 CALL pw_grid_create(pw_big_grid, para_env%group) 243 IF (ns_max == 2*(ns_max/2)) THEN 244 CPABORT("ns_max has to be odd.") 245 END IF 246 npts_s(:) = ns_max 247 ! compute cut-off radius 248 alphasq = alpha**2 249 norm = (2.0_dp*alphasq/pi)**(1.5_dp) 250 cutoff_radius = exp_radius(0, 2.0_dp*alphasq, epsilon, norm) 251 252 CALL dg_pme_grid_setup(cell_ref%hmat, npts_s, cutoff_radius, & 253 pw_small_grid, pw_big_grid, rs_dims=(/para_env%num_pe, 1/), & 254 iounit=output_unit, fft_usage=.TRUE.) 255 ! Write some useful info 256 IF (output_unit > 0) THEN 257 WRITE (output_unit, '( A,T71,E10.4 )') & 258 ' EWALD| Gaussian tolerance (effective) ', epsilon 259 WRITE (output_unit, '( A,T63,3I6 )') & 260 ' EWALD| Small box grid ', pw_small_grid%npts 261 WRITE (output_unit, '( A,T63,3I6 )') & 262 ' EWALD| Full box grid ', pw_big_grid%npts 263 END IF 264 265 ! pw pools initialized 266 NULLIFY (pw_pool) 267 CALL pw_pool_create(pw_pool, pw_grid=pw_big_grid) 268 ewald_pw%pw_big_pool => pw_pool 269 CALL pw_pool_retain(ewald_pw%pw_big_pool) 270 CALL pw_pool_release(pw_pool) 271 272 NULLIFY (pw_pool) 273 CALL pw_pool_create(pw_pool, pw_grid=pw_small_grid) 274 ewald_pw%pw_small_pool => pw_pool 275 CALL pw_pool_retain(ewald_pw%pw_small_pool) 276 CALL pw_pool_release(pw_pool) 277 278 NULLIFY (rs_desc) 279 CALL init_input_type(input_settings, nsmax=MAXVAL(pw_small_grid%npts(1:3)), & 280 rs_grid_section=rs_grid_section, ilevel=1, & 281 higher_grid_layout=(/-1, -1, -1/)) 282 CALL rs_grid_create_descriptor(rs_desc, pw_big_grid, input_settings) 283 284 CALL rs_grid_create(rs, rs_desc) 285 CALL rs_grid_print(rs, output_unit) 286 CALL rs_grid_release(rs) 287 288 CALL cp_print_key_finished_output(output_unit, logger, print_section, "") 289 290 ewald_pw%rs_desc => rs_desc 291 292 CALL rs_grid_retain_descriptor(ewald_pw%rs_desc) 293 CALL rs_grid_release_descriptor(rs_desc) 294 295 CALL pw_grid_release(pw_small_grid) 296 CALL pw_grid_release(pw_big_grid) 297 298 CASE (do_ewald_spme) 299 ! set up the Smooth-Particle-Mesh EWALD sum 300 logger => cp_get_default_logger() 301 output_unit = cp_print_key_unit_nr(logger, print_section, "", extension=".Log") 302 IF (.NOT. ASSOCIATED(ewald_pw%poisson_env)) THEN 303 CALL pw_poisson_create(ewald_pw%poisson_env) 304 END IF 305 CALL pw_grid_create(pw_big_grid, para_env%group) 306 npts_s = gmax 307 CALL pw_grid_setup(cell_ref%hmat, pw_big_grid, grid_span=HALFSPACE, npts=npts_s, spherical=.TRUE., & 308 rs_dims=(/para_env%num_pe, 1/), iounit=output_unit, fft_usage=.TRUE.) 309 310 ! pw pools initialized 311 NULLIFY (pw_pool) 312 CALL pw_pool_create(pw_pool, pw_grid=pw_big_grid) 313 ewald_pw%pw_big_pool => pw_pool 314 CALL pw_pool_retain(ewald_pw%pw_big_pool) 315 CALL pw_pool_release(pw_pool) 316 317 NULLIFY (rs_desc) 318 CALL init_input_type(input_settings, nsmax=o_spline, & 319 rs_grid_section=rs_grid_section, ilevel=1, & 320 higher_grid_layout=(/-1, -1, -1/)) 321 CALL rs_grid_create_descriptor(rs_desc, pw_big_grid, input_settings) 322 323 CALL rs_grid_create(rs, rs_desc) 324 CALL rs_grid_print(rs, output_unit) 325 CALL rs_grid_release(rs) 326 CALL cp_print_key_finished_output(output_unit, logger, print_section, "") 327 328 ewald_pw%rs_desc => rs_desc 329 330 CALL rs_grid_retain_descriptor(ewald_pw%rs_desc) 331 CALL rs_grid_release_descriptor(rs_desc) 332 333 CALL pw_grid_release(pw_big_grid) 334 CASE (do_ewald_none) 335 ! No EWALD sums.. 336 CASE default 337 CPABORT("") 338 END SELECT 339 ! Poisson Environment 340 IF (ASSOCIATED(ewald_pw%poisson_env)) THEN 341 ALLOCATE (pw_pools(1)) 342 pw_pools(1)%pool => ewald_pw%pw_big_pool 343 CALL pw_poisson_read_parameters(poisson_section, poisson_params) 344 poisson_params%ewald_type = ewald_type 345 poisson_params%ewald_o_spline = o_spline 346 poisson_params%ewald_alpha = alpha 347 CALL pw_poisson_set(ewald_pw%poisson_env, cell_hmat=cell%hmat, parameters=poisson_params, & 348 use_level=1, pw_pools=pw_pools) 349 DEALLOCATE (pw_pools) 350 END IF 351 CALL timestop(handle) 352 END SUBROUTINE ewald_pw_init 353 354!****** ewald_pw_types/ewald_pw_get [1.0] * 355 356! ************************************************************************************************** 357!> \brief get the ewald_pw environment to the correct program. 358!> \param ewald_pw ... 359!> \param pw_big_pool ... 360!> \param pw_small_pool ... 361!> \param rs_desc ... 362!> \param poisson_env ... 363!> \param dg ... 364!> \author CJM 365! ************************************************************************************************** 366 SUBROUTINE ewald_pw_get(ewald_pw, pw_big_pool, pw_small_pool, rs_desc, poisson_env, dg) 367 368 TYPE(ewald_pw_type), POINTER :: ewald_pw 369 TYPE(pw_pool_type), OPTIONAL, POINTER :: pw_big_pool, pw_small_pool 370 TYPE(realspace_grid_desc_type), OPTIONAL, POINTER :: rs_desc 371 TYPE(pw_poisson_type), OPTIONAL, POINTER :: poisson_env 372 TYPE(dg_type), OPTIONAL, POINTER :: dg 373 374 CHARACTER(LEN=*), PARAMETER :: routineN = 'ewald_pw_get', routineP = moduleN//':'//routineN 375 376 IF (PRESENT(poisson_env)) poisson_env => ewald_pw%poisson_env 377 IF (PRESENT(pw_big_pool)) pw_big_pool => ewald_pw%pw_big_pool 378 IF (PRESENT(pw_small_pool)) pw_small_pool => ewald_pw%pw_small_pool 379 IF (PRESENT(rs_desc)) rs_desc => ewald_pw%rs_desc 380 IF (PRESENT(dg)) dg => ewald_pw%dg 381 382 END SUBROUTINE ewald_pw_get 383 384!****** ewald_pw_types/ewald_pw_set [1.0] * 385 386! ************************************************************************************************** 387!> \brief set the ewald_pw environment to the correct program. 388!> \param ewald_pw ... 389!> \param pw_big_pool ... 390!> \param pw_small_pool ... 391!> \param rs_desc ... 392!> \param dg ... 393!> \param poisson_env ... 394!> \author CJM 395! ************************************************************************************************** 396 SUBROUTINE ewald_pw_set(ewald_pw, pw_big_pool, pw_small_pool, rs_desc, dg, & 397 poisson_env) 398 399 TYPE(ewald_pw_type), POINTER :: ewald_pw 400 TYPE(pw_pool_type), OPTIONAL, POINTER :: pw_big_pool, pw_small_pool 401 TYPE(realspace_grid_desc_type), OPTIONAL, POINTER :: rs_desc 402 TYPE(dg_type), OPTIONAL, POINTER :: dg 403 TYPE(pw_poisson_type), OPTIONAL, POINTER :: poisson_env 404 405 CHARACTER(LEN=*), PARAMETER :: routineN = 'ewald_pw_set', routineP = moduleN//':'//routineN 406 407 IF (PRESENT(pw_big_pool)) THEN 408 CALL pw_pool_retain(pw_big_pool) 409 CALL pw_pool_release(ewald_pw%pw_big_pool) 410 ewald_pw%pw_big_pool => pw_big_pool 411 ENDIF 412 IF (PRESENT(pw_small_pool)) THEN 413 CALL pw_pool_retain(pw_small_pool) 414 CALL pw_pool_release(ewald_pw%pw_small_pool) 415 ewald_pw%pw_small_pool => pw_small_pool 416 ENDIF 417 IF (PRESENT(rs_desc)) THEN 418 CALL rs_grid_retain_descriptor(rs_desc) 419 CALL rs_grid_release_descriptor(ewald_pw%rs_desc) 420 ewald_pw%rs_desc => rs_desc 421 ENDIF 422 IF (PRESENT(dg)) THEN 423 CALL dg_retain(dg) 424 CALL dg_release(ewald_pw%dg) 425 ewald_pw%dg => dg 426 ENDIF 427 IF (PRESENT(poisson_env)) THEN 428 IF (ASSOCIATED(poisson_env)) & 429 CALL pw_poisson_retain(poisson_env) 430 CALL pw_poisson_release(ewald_pw%poisson_env) 431 ewald_pw%poisson_env => poisson_env 432 END IF 433 434 END SUBROUTINE ewald_pw_set 435 436END MODULE ewald_pw_types 437