1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> JGH FEB-13-2007 : Distributed/replicated realspace grids 9!> Teodoro Laino [tlaino] - University of Zurich - 12.2007 10!> \author CJM NOV-30-2003 11! ************************************************************************************************** 12MODULE ewald_environment_types 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_env, ONLY: cp_para_env_release,& 18 cp_para_env_retain 19 USE cp_para_types, ONLY: cp_para_env_type 20 USE cp_units, ONLY: cp_unit_from_cp2k 21 USE input_cp2k_poisson, ONLY: create_ewald_section 22 USE input_enumeration_types, ONLY: enum_i2c,& 23 enumeration_type 24 USE input_keyword_types, ONLY: keyword_get,& 25 keyword_type 26 USE input_section_types, ONLY: section_get_keyword,& 27 section_release,& 28 section_type,& 29 section_vals_get_subs_vals,& 30 section_vals_release,& 31 section_vals_retain,& 32 section_vals_type,& 33 section_vals_val_get 34 USE kinds, ONLY: dp 35 USE mathconstants, ONLY: twopi 36 USE pw_grid_info, ONLY: pw_grid_n_for_fft 37 USE pw_poisson_types, ONLY: do_ewald_ewald,& 38 do_ewald_none,& 39 do_ewald_pme,& 40 do_ewald_spme 41#include "./base/base_uses.f90" 42 43 IMPLICIT NONE 44 PRIVATE 45 46! ************************************************************************************************** 47!> \brief to build arrays of pointers 48!> \param ewald_env the pointer to the ewald_env 49!> \par History 50!> 11/03 51!> \author CJM 52! ************************************************************************************************** 53 TYPE ewald_environment_type 54 PRIVATE 55 INTEGER :: id_nr, ref_count 56 LOGICAL :: do_multipoles ! Flag for using the multipole code 57 INTEGER :: do_ipol ! Solver for induced dipoles 58 INTEGER :: max_multipole ! max expansion in the multipoles 59 INTEGER :: max_ipol_iter ! max number of interation for induced dipoles 60 INTEGER :: ewald_type ! type of ewald 61 INTEGER :: gmax(3) ! max Miller index 62 INTEGER :: ns_max ! # grid points for small grid (PME) 63 INTEGER :: o_spline ! order of spline (SPME) 64 REAL(KIND=dp) :: precs ! precision achieved when evaluating the real-space part 65 REAL(KIND=dp) :: alpha, rcut ! ewald alpha and real-space cutoff 66 REAL(KIND=dp) :: epsilon ! tolerance for small grid (PME) 67 REAL(KIND=dp) :: eps_pol ! tolerance for convergence of induced dipoles 68 TYPE(cp_para_env_type), POINTER :: para_env 69 TYPE(section_vals_type), POINTER :: poisson_section 70 ! interaction cutoff is required to make the electrostatic interaction 71 ! continuous at a pair distance equal to rcut. this is ignored by the 72 ! multipole code and is otherwise only active when SHIFT_CUTOFF is used. 73 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: interaction_cutoffs 74 ! store current cell, used to rebuild lazily. 75 REAL(KIND=dp), DIMENSION(3, 3) :: cell_hmat = -1.0_dp 76 END TYPE ewald_environment_type 77 78! ************************************************************************************************** 79!> \brief to build arrays of pointers 80!> \param ewald_env the pointer to the ewald_env 81!> \par History 82!> 11/03 83!> \author CJM 84! ************************************************************************************************** 85 TYPE ewald_environment_p_type 86 TYPE(ewald_environment_type), POINTER :: ewald_env 87 END TYPE ewald_environment_p_type 88 89! *** Public data types *** 90 PUBLIC :: ewald_environment_type 91 92! *** Public subroutines *** 93 PUBLIC :: ewald_env_get, & 94 ewald_env_set, & 95 ewald_env_create, & 96 ewald_env_retain, & 97 ewald_env_release, & 98 read_ewald_section, & 99 read_ewald_section_tb 100 101 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewald_environment_types' 102 INTEGER, PRIVATE, SAVE :: last_ewald_env_id_nr = 0 103 104CONTAINS 105 106! ************************************************************************************************** 107!> \brief Purpose: Get the EWALD environment. 108!> \param ewald_env the pointer to the ewald_env 109!> \param ewald_type ... 110!> \param alpha ... 111!> \param eps_pol ... 112!> \param epsilon ... 113!> \param gmax ... 114!> \param ns_max ... 115!> \param o_spline ... 116!> \param group ... 117!> \param para_env ... 118!> \param id_nr ... 119!> \param poisson_section ... 120!> \param precs ... 121!> \param rcut ... 122!> \param do_multipoles ... 123!> \param max_multipole ... 124!> \param do_ipol ... 125!> \param max_ipol_iter ... 126!> \param interaction_cutoffs ... 127!> \param cell_hmat ... 128!> \par History 129!> 11/03 130!> \author CJM 131! ************************************************************************************************** 132 SUBROUTINE ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, & 133 gmax, ns_max, o_spline, group, para_env, id_nr, poisson_section, precs, & 134 rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, & 135 interaction_cutoffs, cell_hmat) 136 TYPE(ewald_environment_type), POINTER :: ewald_env 137 INTEGER, OPTIONAL :: ewald_type 138 REAL(KIND=dp), OPTIONAL :: alpha, eps_pol, epsilon 139 INTEGER, OPTIONAL :: gmax(3), ns_max, o_spline, group 140 TYPE(cp_para_env_type), OPTIONAL, POINTER :: para_env 141 INTEGER, INTENT(OUT), OPTIONAL :: id_nr 142 TYPE(section_vals_type), OPTIONAL, POINTER :: poisson_section 143 REAL(KIND=dp), OPTIONAL :: precs, rcut 144 LOGICAL, INTENT(OUT), OPTIONAL :: do_multipoles 145 INTEGER, INTENT(OUT), OPTIONAL :: max_multipole, do_ipol, max_ipol_iter 146 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 147 POINTER :: interaction_cutoffs 148 REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: cell_hmat 149 150 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_env_get', routineP = moduleN//':'//routineN 151 152 CPASSERT(ASSOCIATED(ewald_env)) 153 154 IF (PRESENT(id_nr)) id_nr = ewald_env%id_nr 155 IF (PRESENT(ewald_type)) ewald_type = ewald_env%ewald_type 156 IF (PRESENT(do_multipoles)) do_multipoles = ewald_env%do_multipoles 157 IF (PRESENT(do_ipol)) do_ipol = ewald_env%do_ipol 158 IF (PRESENT(max_multipole)) max_multipole = ewald_env%max_multipole 159 IF (PRESENT(max_ipol_iter)) max_ipol_iter = ewald_env%max_ipol_iter 160 IF (PRESENT(alpha)) alpha = ewald_env%alpha 161 IF (PRESENT(precs)) precs = ewald_env%precs 162 IF (PRESENT(rcut)) rcut = ewald_env%rcut 163 IF (PRESENT(epsilon)) epsilon = ewald_env%epsilon 164 IF (PRESENT(eps_pol)) eps_pol = ewald_env%eps_pol 165 IF (PRESENT(gmax)) gmax = ewald_env%gmax 166 IF (PRESENT(ns_max)) ns_max = ewald_env%ns_max 167 IF (PRESENT(o_spline)) o_spline = ewald_env%o_spline 168 IF (PRESENT(group)) group = ewald_env%para_env%group 169 IF (PRESENT(para_env)) para_env => ewald_env%para_env 170 IF (PRESENT(poisson_section)) poisson_section => ewald_env%poisson_section 171 IF (PRESENT(interaction_cutoffs)) interaction_cutoffs => & 172 ewald_env%interaction_cutoffs 173 IF (PRESENT(cell_hmat)) cell_hmat = ewald_env%cell_hmat 174 END SUBROUTINE ewald_env_get 175 176! ************************************************************************************************** 177!> \brief Purpose: Set the EWALD environment. 178!> \param ewald_env the pointer to the ewald_env 179!> \param ewald_type ... 180!> \param alpha ... 181!> \param epsilon ... 182!> \param eps_pol ... 183!> \param gmax ... 184!> \param ns_max ... 185!> \param precs ... 186!> \param o_spline ... 187!> \param para_env ... 188!> \param id_nr ... 189!> \param poisson_section ... 190!> \param interaction_cutoffs ... 191!> \param cell_hmat ... 192!> \par History 193!> 11/03 194!> \author CJM 195! ************************************************************************************************** 196 SUBROUTINE ewald_env_set(ewald_env, ewald_type, alpha, epsilon, eps_pol, & 197 gmax, ns_max, precs, o_spline, para_env, id_nr, poisson_section, & 198 interaction_cutoffs, cell_hmat) 199 200 TYPE(ewald_environment_type), POINTER :: ewald_env 201 INTEGER, OPTIONAL :: ewald_type 202 REAL(KIND=dp), OPTIONAL :: alpha, epsilon, eps_pol 203 INTEGER, OPTIONAL :: gmax(3), ns_max 204 REAL(KIND=dp), OPTIONAL :: precs 205 INTEGER, OPTIONAL :: o_spline 206 TYPE(cp_para_env_type), OPTIONAL, POINTER :: para_env 207 INTEGER, INTENT(IN), OPTIONAL :: id_nr 208 TYPE(section_vals_type), OPTIONAL, POINTER :: poisson_section 209 REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, & 210 POINTER :: interaction_cutoffs 211 REAL(KIND=dp), DIMENSION(3, 3), OPTIONAL :: cell_hmat 212 213 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_env_set', routineP = moduleN//':'//routineN 214 215 CPASSERT(ASSOCIATED(ewald_env)) 216 217 IF (PRESENT(id_nr)) ewald_env%id_nr = id_nr 218 IF (PRESENT(ewald_type)) ewald_env%ewald_type = ewald_type 219 IF (PRESENT(alpha)) ewald_env%alpha = alpha 220 IF (PRESENT(precs)) ewald_env%precs = precs 221 IF (PRESENT(epsilon)) ewald_env%epsilon = epsilon 222 IF (PRESENT(eps_pol)) ewald_env%eps_pol = eps_pol 223 IF (PRESENT(gmax)) ewald_env%gmax = gmax 224 IF (PRESENT(ns_max)) ewald_env%ns_max = ns_max 225 IF (PRESENT(o_spline)) ewald_env%o_spline = o_spline 226 IF (PRESENT(para_env)) ewald_env%para_env => para_env 227 IF (PRESENT(poisson_section)) THEN 228 CALL section_vals_retain(poisson_section) 229 CALL section_vals_release(ewald_env%poisson_section) 230 ewald_env%poisson_section => poisson_section 231 END IF 232 IF (PRESENT(interaction_cutoffs)) ewald_env%interaction_cutoffs => & 233 interaction_cutoffs 234 IF (PRESENT(cell_hmat)) ewald_env%cell_hmat = cell_hmat 235 END SUBROUTINE ewald_env_set 236 237! ************************************************************************************************** 238!> \brief allocates and intitializes a ewald_env 239!> \param ewald_env the object to create 240!> \param para_env ... 241!> \par History 242!> 12.2002 created [fawzi] 243!> \author Fawzi Mohamed 244! ************************************************************************************************** 245 SUBROUTINE ewald_env_create(ewald_env, para_env) 246 TYPE(ewald_environment_type), POINTER :: ewald_env 247 TYPE(cp_para_env_type), POINTER :: para_env 248 249 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_env_create', & 250 routineP = moduleN//':'//routineN 251 252 ALLOCATE (ewald_env) 253 ewald_env%ref_count = 1 254 last_ewald_env_id_nr = last_ewald_env_id_nr + 1 255 ewald_env%id_nr = last_ewald_env_id_nr 256 NULLIFY (ewald_env%poisson_section) 257 CALL cp_para_env_retain(para_env) 258 ewald_env%para_env => para_env 259 NULLIFY (ewald_env%interaction_cutoffs) ! allocated and initialized later 260 END SUBROUTINE ewald_env_create 261 262! ************************************************************************************************** 263!> \brief retains the given ewald_env (see doc/ReferenceCounting.html) 264!> \param ewald_env the object to retain 265!> \par History 266!> 12.2002 created [fawzi] 267!> \author Fawzi Mohamed 268! ************************************************************************************************** 269 SUBROUTINE ewald_env_retain(ewald_env) 270 TYPE(ewald_environment_type), POINTER :: ewald_env 271 272 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_env_retain', & 273 routineP = moduleN//':'//routineN 274 275 CPASSERT(ASSOCIATED(ewald_env)) 276 CPASSERT(ewald_env%ref_count > 0) 277 ewald_env%ref_count = ewald_env%ref_count + 1 278 END SUBROUTINE ewald_env_retain 279 280! ************************************************************************************************** 281!> \brief releases the given ewald_env (see doc/ReferenceCounting.html) 282!> \param ewald_env the object to release 283!> \par History 284!> 12.2002 created [fawzi] 285!> \author Fawzi Mohamed 286! ************************************************************************************************** 287 SUBROUTINE ewald_env_release(ewald_env) 288 TYPE(ewald_environment_type), POINTER :: ewald_env 289 290 CHARACTER(len=*), PARAMETER :: routineN = 'ewald_env_release', & 291 routineP = moduleN//':'//routineN 292 293 IF (ASSOCIATED(ewald_env)) THEN 294 CPASSERT(ewald_env%ref_count > 0) 295 ewald_env%ref_count = ewald_env%ref_count - 1 296 IF (ewald_env%ref_count < 1) THEN 297 CALL cp_para_env_release(ewald_env%para_env) 298 CALL section_vals_release(ewald_env%poisson_section) 299 IF (ASSOCIATED(ewald_env%interaction_cutoffs)) THEN 300 DEALLOCATE (ewald_env%interaction_cutoffs) 301 ENDIF 302 DEALLOCATE (ewald_env) 303 ENDIF 304 END IF 305 NULLIFY (ewald_env) 306 END SUBROUTINE ewald_env_release 307 308! ************************************************************************************************** 309!> \brief Purpose: read the EWALD section 310!> \param ewald_env the pointer to the ewald_env 311!> \param ewald_section ... 312!> \author Teodoro Laino [tlaino] -University of Zurich - 2005 313! ************************************************************************************************** 314 SUBROUTINE read_ewald_section(ewald_env, ewald_section) 315 TYPE(ewald_environment_type), POINTER :: ewald_env 316 TYPE(section_vals_type), POINTER :: ewald_section 317 318 CHARACTER(len=*), PARAMETER :: routineN = 'read_ewald_section', & 319 routineP = moduleN//':'//routineN 320 321 INTEGER :: iw 322 INTEGER, DIMENSION(:), POINTER :: gmax_read 323 LOGICAL :: explicit 324 REAL(KIND=dp) :: dummy 325 TYPE(cp_logger_type), POINTER :: logger 326 TYPE(enumeration_type), POINTER :: enum 327 TYPE(keyword_type), POINTER :: keyword 328 TYPE(section_type), POINTER :: section 329 TYPE(section_vals_type), POINTER :: multipole_section 330 331 NULLIFY (enum, keyword, section, multipole_section) 332 logger => cp_get_default_logger() 333 CALL section_vals_val_get(ewald_section, "EWALD_TYPE", i_val=ewald_env%ewald_type) 334 CALL section_vals_val_get(ewald_section, "ALPHA", r_val=ewald_env%alpha) 335 CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs) 336 337 IF (ewald_env%ewald_type == do_ewald_none) THEN 338 ewald_env%rcut = 0.0_dp 339 ELSE 340 CALL section_vals_val_get(ewald_section, "RCUT", explicit=explicit) 341 IF (explicit) THEN 342 CALL section_vals_val_get(ewald_section, "RCUT", r_val=ewald_env%rcut) 343 ELSE 344 ewald_env%rcut = find_ewald_optimal_value(ewald_env%precs)/ewald_env%alpha 345 ENDIF 346 END IF 347 ! we have no defaults for gmax, gmax is only needed for ewald and spme 348 SELECT CASE (ewald_env%ewald_type) 349 CASE (do_ewald_ewald, do_ewald_spme) 350 CALL section_vals_val_get(ewald_section, "GMAX", i_vals=gmax_read) 351 SELECT CASE (SIZE(gmax_read, 1)) 352 CASE (1) 353 ewald_env%gmax = gmax_read(1) 354 CASE (3) 355 ewald_env%gmax = gmax_read 356 CASE DEFAULT 357 CPABORT("") 358 END SELECT 359 IF (ewald_env%ewald_type == do_ewald_spme) THEN 360 CALL section_vals_val_get(ewald_section, "O_SPLINE", i_val=ewald_env%o_spline) 361 END IF 362 CASE (do_ewald_pme) 363 CALL section_vals_val_get(ewald_section, "NS_MAX", i_val=ewald_env%ns_max) 364 CALL section_vals_val_get(ewald_section, "EPSILON", r_val=ewald_env%epsilon) 365 CASE DEFAULT 366 ! this should not be used for do_ewald_none 367 ewald_env%gmax = HUGE(0) 368 ewald_env%ns_max = HUGE(0) 369 END SELECT 370 371 ! Multipoles 372 multipole_section => section_vals_get_subs_vals(ewald_section, "MULTIPOLES") 373 CALL section_vals_val_get(multipole_section, "_SECTION_PARAMETERS_", l_val=ewald_env%do_multipoles) 374 CALL section_vals_val_get(multipole_section, "POL_SCF", i_val=ewald_env%do_ipol) 375 CALL section_vals_val_get(multipole_section, "EPS_POL", r_val=ewald_env%eps_pol) 376 IF (ewald_env%do_multipoles) THEN 377 SELECT CASE (ewald_env%ewald_type) 378 CASE (do_ewald_ewald) 379 CALL section_vals_val_get(multipole_section, "MAX_MULTIPOLE_EXPANSION", & 380 i_val=ewald_env%max_multipole) 381 CALL section_vals_val_get(multipole_section, "MAX_IPOL_ITER", i_val=ewald_env%max_ipol_iter) 382 CASE DEFAULT 383 CPABORT("Multipole code works at the moment only with standard EWALD sums.") 384 END SELECT 385 END IF 386 387 iw = cp_print_key_unit_nr(logger, ewald_section, "PRINT%PROGRAM_RUN_INFO", & 388 extension=".log") 389 IF (iw > 0) THEN 390 NULLIFY (keyword, enum) 391 CALL create_ewald_section(section) 392 IF (ewald_env%ewald_type /= do_ewald_none) THEN 393 keyword => section_get_keyword(section, "EWALD_TYPE") 394 CALL keyword_get(keyword, enum=enum) 395 WRITE (iw, '(/,T2,"EWALD| ",A,T67,A14 )') 'Summation is done by:', & 396 ADJUSTR(TRIM(enum_i2c(enum, ewald_env%ewald_type))) 397 IF (ewald_env%do_multipoles) THEN 398 NULLIFY (keyword, enum) 399 keyword => section_get_keyword(section, "MULTIPOLES%MAX_MULTIPOLE_EXPANSION") 400 CALL keyword_get(keyword, enum=enum) 401 WRITE (iw, '( T2,"EWALD| ",A )') 'Enabled Multipole Method' 402 WRITE (iw, '( T2,"EWALD| ",A,T67,A14 )') 'Max Term in Multipole Expansion :', & 403 ADJUSTR(TRIM(enum_i2c(enum, ewald_env%max_multipole))) 404 WRITE (iw, '( T2,"EWALD| ",A,T67,3I10 )') 'Max number Iterations for IPOL :', & 405 ewald_env%max_ipol_iter 406 END IF 407 dummy = cp_unit_from_cp2k(ewald_env%alpha, "angstrom^-1") 408 WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') & 409 'Alpha parameter [', 'ANGSTROM^-1', ']', dummy 410 dummy = cp_unit_from_cp2k(ewald_env%rcut, "angstrom") 411 WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') & 412 'Real Space Cutoff [', 'ANGSTROM', ']', dummy 413 414 SELECT CASE (ewald_env%ewald_type) 415 CASE (do_ewald_ewald) 416 WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') & 417 'G-space max. Miller index', ewald_env%gmax 418 CASE (do_ewald_pme) 419 WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') & 420 'Max small-grid points (input) ', ewald_env%ns_max 421 WRITE (iw, '( T2,"EWALD| ",A,T71,E10.4 )') & 422 'Gaussian tolerance (input) ', ewald_env%epsilon 423 CASE (do_ewald_spme) 424 WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') & 425 'G-space max. Miller index', ewald_env%gmax 426 WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') & 427 'Spline interpolation order ', ewald_env%o_spline 428 CASE DEFAULT 429 CPABORT("") 430 END SELECT 431 ELSE 432 WRITE (iw, '( T2,"EWALD| ",T73, A )') 'not used' 433 END IF 434 CALL section_release(section) 435 END IF 436 CALL cp_print_key_finished_output(iw, logger, ewald_section, & 437 "PRINT%PROGRAM_RUN_INFO") 438 439 END SUBROUTINE read_ewald_section 440 441! ************************************************************************************************** 442!> \brief Purpose: read the EWALD section for TB methods 443!> \param ewald_env the pointer to the ewald_env 444!> \param ewald_section ... 445!> \param hmat ... 446!> \author JGH 447! ************************************************************************************************** 448 SUBROUTINE read_ewald_section_tb(ewald_env, ewald_section, hmat) 449 TYPE(ewald_environment_type), POINTER :: ewald_env 450 TYPE(section_vals_type), POINTER :: ewald_section 451 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat 452 453 CHARACTER(len=*), PARAMETER :: routineN = 'read_ewald_section_tb', & 454 routineP = moduleN//':'//routineN 455 456 INTEGER :: i, iw, n(3) 457 INTEGER, DIMENSION(:), POINTER :: gmax_read 458 LOGICAL :: explicit 459 REAL(KIND=dp) :: alat, cutoff, dummy 460 TYPE(cp_logger_type), POINTER :: logger 461 462 logger => cp_get_default_logger() 463 464 ewald_env%do_multipoles = .FALSE. 465 ewald_env%do_ipol = 0 466 ewald_env%eps_pol = 1.e-12_dp 467 ewald_env%max_multipole = 0 468 ewald_env%max_ipol_iter = 0 469 ewald_env%epsilon = 1.e-12_dp 470 ewald_env%ns_max = HUGE(0) 471 472 CALL section_vals_val_get(ewald_section, "EWALD_TYPE", explicit=explicit) 473 IF (explicit) THEN 474 CALL section_vals_val_get(ewald_section, "EWALD_TYPE", i_val=ewald_env%ewald_type) 475 IF (ewald_env%ewald_type /= do_ewald_spme) THEN 476 CPABORT("TB needs EWALD_TYPE SPME") 477 END IF 478 ELSE 479 ewald_env%ewald_type = do_ewald_spme 480 ENDIF 481 482 CALL section_vals_val_get(ewald_section, "ALPHA", explicit=explicit) 483 IF (explicit) THEN 484 CALL section_vals_val_get(ewald_section, "ALPHA", r_val=ewald_env%alpha) 485 ELSE 486 ewald_env%alpha = 1.0_dp 487 ENDIF 488 489 CALL section_vals_val_get(ewald_section, "EWALD_ACCURACY", r_val=ewald_env%precs) 490 CALL section_vals_val_get(ewald_section, "O_SPLINE", i_val=ewald_env%o_spline) 491 492 CALL section_vals_val_get(ewald_section, "RCUT", explicit=explicit) 493 IF (explicit) THEN 494 CALL section_vals_val_get(ewald_section, "RCUT", r_val=ewald_env%rcut) 495 ELSE 496 ewald_env%rcut = find_ewald_optimal_value(ewald_env%precs)/ewald_env%alpha 497 ENDIF 498 499 CALL section_vals_val_get(ewald_section, "GMAX", explicit=explicit) 500 IF (explicit) THEN 501 CALL section_vals_val_get(ewald_section, "GMAX", i_vals=gmax_read) 502 SELECT CASE (SIZE(gmax_read, 1)) 503 CASE (1) 504 ewald_env%gmax = gmax_read(1) 505 CASE (3) 506 ewald_env%gmax = gmax_read 507 CASE DEFAULT 508 CPABORT("") 509 END SELECT 510 ELSE 511 ! set GMAX using ECUT=alpha*45 Ry 512 cutoff = 45._dp*ewald_env%alpha 513 DO i = 1, 3 514 alat = SUM(hmat(:, i)**2) 515 CPASSERT(alat /= 0._dp) 516 ewald_env%gmax(i) = 2*FLOOR(SQRT(2.0_dp*cutoff*alat)/twopi) + 1 517 ENDDO 518 ENDIF 519 n = ewald_env%gmax 520 ewald_env%gmax = pw_grid_n_for_fft(n, odd=.TRUE.) 521 522 iw = cp_print_key_unit_nr(logger, ewald_section, "PRINT%PROGRAM_RUN_INFO", & 523 extension=".log") 524 IF (iw > 0) THEN 525 WRITE (iw, '(/,T2,"EWALD| ",A,T67,A14 )') 'Summation is done by:', ADJUSTR("SPME") 526 dummy = cp_unit_from_cp2k(ewald_env%alpha, "angstrom^-1") 527 WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') & 528 'Alpha parameter [', 'ANGSTROM^-1', ']', dummy 529 dummy = cp_unit_from_cp2k(ewald_env%rcut, "angstrom") 530 WRITE (iw, '( T2,"EWALD| ",A,A18,A,T71,F10.4 )') & 531 'Real Space Cutoff [', 'ANGSTROM', ']', dummy 532 WRITE (iw, '( T2,"EWALD| ",A,T51,3I10 )') & 533 'G-space max. Miller index', ewald_env%gmax 534 WRITE (iw, '( T2,"EWALD| ",A,T71,I10 )') & 535 'Spline interpolation order ', ewald_env%o_spline 536 END IF 537 CALL cp_print_key_finished_output(iw, logger, ewald_section, & 538 "PRINT%PROGRAM_RUN_INFO") 539 540 END SUBROUTINE read_ewald_section_tb 541 542! ************************************************************************************************** 543!> \brief triggers (by bisection) the optimal value for EWALD parameter x 544!> EXP(-x^2)/x^2 = EWALD_ACCURACY 545!> \param precs ... 546!> \return ... 547!> \author Teodoro Laino [tlaino] - University of Zurich - 12.2007 548! ************************************************************************************************** 549 FUNCTION find_ewald_optimal_value(precs) RESULT(value) 550 REAL(KIND=dp) :: precs, value 551 552 CHARACTER(len=*), PARAMETER :: routineN = 'find_ewald_optimal_value', & 553 routineP = moduleN//':'//routineN 554 555 REAL(KIND=dp) :: func, func1, func2, s, s1, s2 556 557 s = 0.1_dp 558 func = EXP(-s**2)/s**2 - precs 559 CPASSERT(func > 0.0_dp) 560 DO WHILE (func > 0.0_dp) 561 s = s + 0.1_dp 562 func = EXP(-s**2)/s**2 - precs 563 END DO 564 s2 = s 565 s1 = s - 0.1_dp 566 ! Start bisection 567 DO WHILE (.TRUE.) 568 func2 = EXP(-s2**2)/s2**2 - precs 569 func1 = EXP(-s1**2)/s1**2 - precs 570 CPASSERT(func1 >= 0) 571 CPASSERT(func2 <= 0) 572 s = 0.5_dp*(s1 + s2) 573 func = EXP(-s**2)/s**2 - precs 574 IF (func > 0.0_dp) THEN 575 s1 = s 576 ELSE IF (func < 0.0_dp) THEN 577 s2 = s 578 END IF 579 IF (ABS(func) < 100.0_dp*EPSILON(0.0_dp)) EXIT 580 END DO 581 value = s 582 END FUNCTION find_ewald_optimal_value 583 584END MODULE ewald_environment_types 585 586