1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Types and initialization / release routines for Minimax-Ewald method for electron 8!> repulsion integrals. 9!> \par History 10!> 2015 09 created 11!> \author Patrick Seewald 12! ************************************************************************************************** 13 14MODULE eri_mme_types 15 16 USE cp_para_types, ONLY: cp_para_env_type 17 USE eri_mme_error_control, ONLY: calibrate_cutoff,& 18 cutoff_minimax_error,& 19 minimax_error 20 USE eri_mme_gaussian, ONLY: eri_mme_coulomb,& 21 eri_mme_longrange,& 22 eri_mme_yukawa,& 23 get_minimax_coeff_v_gspace 24 USE eri_mme_util, ONLY: G_abs_min,& 25 R_abs_min 26 USE kinds, ONLY: dp 27 USE mathlib, ONLY: det_3x3,& 28 inv_3x3 29 USE orbital_pointers, ONLY: init_orbital_pointers 30#include "../base/base_uses.f90" 31 32 IMPLICIT NONE 33 34 PRIVATE 35 36 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE. 37 38 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_types' 39 40 INTEGER, PARAMETER, PUBLIC :: n_minimax_max = 53 41 42 PUBLIC :: eri_mme_param, & 43 eri_mme_init, & 44 eri_mme_release, & 45 eri_mme_set_params, & 46 eri_mme_print_grid_info, & 47 get_minimax_from_cutoff, & 48 eri_mme_coulomb, & 49 eri_mme_yukawa, & 50 eri_mme_longrange, & 51 eri_mme_set_potential 52 53 TYPE minimax_grid 54 REAL(KIND=dp) :: cutoff 55 INTEGER :: n_minimax 56 REAL(KIND=dp), POINTER, & 57 DIMENSION(:) :: minimax_aw => NULL() 58 REAL(KIND=dp) :: error 59 END TYPE 60 61 TYPE eri_mme_param 62 INTEGER :: n_minimax 63 REAL(KIND=dp), DIMENSION(3, 3) :: hmat, h_inv 64 REAL(KIND=dp) :: vol 65 LOGICAL :: is_ortho 66 REAL(KIND=dp) :: cutoff 67 LOGICAL :: do_calib_cutoff 68 LOGICAL :: do_error_est 69 LOGICAL :: print_calib 70 REAL(KIND=dp) :: cutoff_min, cutoff_max, cutoff_delta, & 71 cutoff_eps 72 REAL(KIND=dp) :: err_mm, err_c 73 REAL(KIND=dp) :: mm_delta 74 REAL(KIND=dp) :: G_min, R_min 75 LOGICAL :: is_valid 76 LOGICAL :: debug 77 REAL(KIND=dp) :: debug_delta 78 INTEGER :: debug_nsum 79 REAL(KIND=dp) :: C_mm 80 INTEGER :: unit_nr 81 REAL(KIND=dp) :: sum_precision 82 INTEGER :: n_grids 83 TYPE(minimax_grid), DIMENSION(:), & 84 ALLOCATABLE :: minimax_grid 85 REAL(KIND=dp) :: zet_max, zet_min 86 INTEGER :: l_mm, l_max_zet 87 INTEGER :: potential 88 REAL(KIND=dp) :: pot_par 89 END TYPE eri_mme_param 90 91CONTAINS 92 93! ************************************************************************************************** 94!> \brief ... 95!> \param param ... 96!> \param n_minimax ... 97!> \param cutoff ... 98!> \param do_calib_cutoff ... 99!> \param do_error_est ... 100!> \param cutoff_min ... 101!> \param cutoff_max ... 102!> \param cutoff_eps ... 103!> \param cutoff_delta ... 104!> \param sum_precision ... 105!> \param debug ... 106!> \param debug_delta ... 107!> \param debug_nsum ... 108!> \param unit_nr ... 109!> \param print_calib ... 110! ************************************************************************************************** 111 SUBROUTINE eri_mme_init(param, n_minimax, cutoff, do_calib_cutoff, do_error_est, & 112 cutoff_min, cutoff_max, cutoff_eps, cutoff_delta, sum_precision, & 113 debug, debug_delta, debug_nsum, unit_nr, print_calib) 114 TYPE(eri_mme_param), INTENT(OUT) :: param 115 INTEGER, INTENT(IN) :: n_minimax 116 REAL(KIND=dp), INTENT(IN) :: cutoff 117 LOGICAL, INTENT(IN) :: do_calib_cutoff, do_error_est 118 REAL(KIND=dp), INTENT(IN) :: cutoff_min, cutoff_max, cutoff_eps, & 119 cutoff_delta, sum_precision 120 LOGICAL, INTENT(IN) :: debug 121 REAL(KIND=dp), INTENT(IN) :: debug_delta 122 INTEGER, INTENT(IN) :: debug_nsum, unit_nr 123 LOGICAL, INTENT(IN) :: print_calib 124 125 CHARACTER(len=2) :: string 126 127 WRITE (string, '(I2)') n_minimax_max 128 IF (n_minimax .GT. n_minimax_max) & 129 CPABORT("The maximum allowed number of minimax points N_MINIMAX is "//TRIM(string)) 130 131 param%n_minimax = n_minimax 132 param%n_grids = 1 133 param%cutoff = cutoff 134 param%do_calib_cutoff = do_calib_cutoff 135 param%do_error_est = do_error_est 136 param%cutoff_min = cutoff_min 137 param%cutoff_max = cutoff_max 138 param%cutoff_eps = cutoff_eps 139 param%cutoff_delta = cutoff_delta 140 param%sum_precision = sum_precision 141 param%debug = debug 142 param%debug_delta = debug_delta 143 param%debug_nsum = debug_nsum 144 param%print_calib = print_calib 145 param%unit_nr = unit_nr 146 param%err_mm = -1.0_dp 147 param%err_c = -1.0_dp 148 149 param%is_valid = .FALSE. 150 ALLOCATE (param%minimax_grid(param%n_grids)) 151 END SUBROUTINE eri_mme_init 152 153! ************************************************************************************************** 154!> \brief Set parameters for MME method with manual specification of basis parameters. 155!> Takes care of cutoff calibration if requested. 156!> \param param ... 157!> \param hmat ... 158!> \param is_ortho ... 159!> \param zet_min Exponent used to estimate error of minimax approximation. 160!> \param zet_max Exponent used to estimate error of finite cutoff. 161!> \param l_max_zet Total ang. mom. quantum numbers l to be combined with exponents in 162!> zet_max. 163!> \param l_max Maximum total angular momentum quantum number 164!> \param para_env ... 165!> \param potential potential to use. Accepts the following values: 166!> 1: coulomb potential V(r)=1/r 167!> 2: yukawa potential V(r)=e(-a*r)/r 168!> 3: long-range coulomb erf(a*r)/r 169!> \param pot_par potential parameter a for yukawa V(r)=e(-a*r)/r or long-range coulomb V(r)=erf(a*r)/r 170! ************************************************************************************************** 171 SUBROUTINE eri_mme_set_params(param, hmat, is_ortho, zet_min, zet_max, l_max_zet, l_max, para_env, & 172 potential, pot_par) 173 TYPE(eri_mme_param), INTENT(INOUT) :: param 174 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat 175 LOGICAL, INTENT(IN) :: is_ortho 176 REAL(KIND=dp), INTENT(IN) :: zet_min, zet_max 177 INTEGER, INTENT(IN) :: l_max_zet, l_max 178 TYPE(cp_para_env_type), INTENT(IN), OPTIONAL, & 179 POINTER :: para_env 180 INTEGER, INTENT(IN), OPTIONAL :: potential 181 REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par 182 183 CHARACTER(LEN=*), PARAMETER :: routineN = 'eri_mme_set_params', & 184 routineP = moduleN//':'//routineN 185 186 INTEGER :: handle, l_mm, n_grids 187 LOGICAL :: s_only 188 REAL(KIND=dp) :: cutoff 189 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw 190 191 CALL timeset(routineN, handle) 192 193 ! Note: in MP2 default logger hacked and does not use global default print level 194 s_only = l_max .EQ. 0 195 196 CALL init_orbital_pointers(3*l_max) ! allow for orbital pointers of combined index 197 198 ! l values for minimax error estimate (l_mm) and for cutoff error estimate (l_max_zet) 199 l_mm = MERGE(0, 1, s_only) 200 201 ! cell info 202 ! Note: we recompute basic quantities from hmat to avoid dependency on cp2k cell type 203 param%hmat = hmat 204 param%h_inv = inv_3x3(hmat) 205 param%vol = ABS(det_3x3(hmat)) 206 param%is_ortho = is_ortho 207 208 ! Minimum lattice vectors 209 param%G_min = G_abs_min(param%h_inv) 210 param%R_min = R_abs_min(param%hmat) 211 212 ! Minimum and maximum exponents 213 param%zet_max = zet_max 214 param%zet_min = zet_min 215 param%l_max_zet = l_max_zet 216 param%l_mm = l_mm 217 218 ! cutoff calibration not yet implemented for general cell 219 IF (.NOT. param%is_ortho) THEN 220 param%do_calib_cutoff = .FALSE. 221 param%do_error_est = .FALSE. 222 ENDIF 223 224 n_grids = param%n_grids 225 226 ! Cutoff calibration and error estimate for orthorhombic cell 227 ! Here we assume Coulomb potential which should give an upper bound error also for the other 228 ! potentials 229 IF (param%do_calib_cutoff) THEN 230 CALL calibrate_cutoff(param%hmat, param%h_inv, param%G_min, param%vol, & 231 zet_min, l_mm, zet_max, l_max_zet, param%n_minimax, & 232 param%cutoff_min, param%cutoff_max, param%cutoff_eps, & 233 param%cutoff_delta, cutoff, param%err_mm, param%err_c, & 234 param%C_mm, para_env, param%print_calib, param%unit_nr) 235 236 param%cutoff = cutoff 237 ELSE IF (param%do_error_est) THEN 238 ALLOCATE (minimax_aw(2*param%n_minimax)) 239 CALL cutoff_minimax_error(param%cutoff, param%hmat, param%h_inv, param%vol, param%G_min, & 240 zet_min, l_mm, zet_max, l_max_zet, param%n_minimax, & 241 minimax_aw, param%err_mm, param%err_c, param%C_mm, para_env) 242 DEALLOCATE (minimax_aw) 243 ENDIF 244 245 param%is_valid = .TRUE. 246 247 CALL eri_mme_set_potential(param, potential=potential, pot_par=pot_par) 248 249 CALL timestop(handle) 250 END SUBROUTINE eri_mme_set_params 251 252! ************************************************************************************************** 253!> \brief ... 254!> \param param ... 255!> \param potential potential to use. Accepts the following values: 256!> 1: coulomb potential V(r)=1/r 257!> 2: yukawa potential V(r)=e(-a*r)/r 258!> 3: long-range coulomb erf(a*r)/r 259!> \param pot_par potential parameter a for yukawa V(r)=e(-a*r)/r or long-range coulomb V(r)=erf(a*r)/r 260! ************************************************************************************************** 261 SUBROUTINE eri_mme_set_potential(param, potential, pot_par) 262 TYPE(eri_mme_param), INTENT(INOUT) :: param 263 INTEGER, INTENT(IN), OPTIONAL :: potential 264 REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par 265 266 REAL(KIND=dp) :: cutoff_max, cutoff_min, cutoff_rel 267 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw 268 269 CPASSERT(param%is_valid) 270 271 IF (PRESENT(potential)) THEN 272 param%potential = potential 273 ELSE 274 param%potential = eri_mme_coulomb 275 ENDIF 276 277 IF (PRESENT(pot_par)) THEN 278 param%pot_par = pot_par 279 ELSE 280 param%pot_par = 0.0_dp 281 ENDIF 282 283 ALLOCATE (minimax_aw(2*param%n_minimax)) 284 285 CALL minimax_error(param%cutoff, param%hmat, param%vol, param%G_min, param%zet_min, param%l_mm, & 286 param%n_minimax, minimax_aw, param%err_mm, param%mm_delta, potential=potential, pot_par=pot_par) 287 288 DEALLOCATE (minimax_aw) 289 290 CPASSERT(param%zet_max + 1.0E-12 .GT. param%zet_min) 291 CPASSERT(param%n_grids .GE. 1) 292 293 cutoff_max = param%cutoff 294 cutoff_rel = param%cutoff/param%zet_max 295 cutoff_min = param%zet_min*cutoff_rel 296 297 CALL eri_mme_destroy_minimax_grids(param%minimax_grid) 298 ALLOCATE (param%minimax_grid(param%n_grids)) 299 300 CALL eri_mme_create_minimax_grids(param%n_grids, param%minimax_grid, param%n_minimax, & 301 cutoff_max, cutoff_min, param%G_min, & 302 param%mm_delta, potential=potential, pot_par=pot_par) 303 304 END SUBROUTINE 305 306! ************************************************************************************************** 307!> \brief ... 308!> \param n_grids ... 309!> \param minimax_grids ... 310!> \param n_minimax ... 311!> \param cutoff_max ... 312!> \param cutoff_min ... 313!> \param G_min ... 314!> \param target_error ... 315!> \param potential ... 316!> \param pot_par ... 317! ************************************************************************************************** 318 SUBROUTINE eri_mme_create_minimax_grids(n_grids, minimax_grids, n_minimax, & 319 cutoff_max, cutoff_min, G_min, & 320 target_error, potential, pot_par) 321 INTEGER, INTENT(IN) :: n_grids 322 TYPE(minimax_grid), DIMENSION(n_grids), & 323 INTENT(OUT) :: minimax_grids 324 INTEGER, INTENT(IN) :: n_minimax 325 REAL(KIND=dp), INTENT(IN) :: cutoff_max, cutoff_min, G_min, & 326 target_error 327 INTEGER, INTENT(IN), OPTIONAL :: potential 328 REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par 329 330 INTEGER :: i_grid, n_mm 331 REAL(KIND=dp) :: cutoff, cutoff_delta, err_mm, err_mm_prev 332 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: minimax_aw, minimax_aw_prev 333 334 cutoff_delta = (cutoff_max/cutoff_min)**(1.0_dp/(n_grids)) 335 cutoff = cutoff_max 336 337 ALLOCATE (minimax_aw(2*n_minimax)) 338 ! for first grid (for max. cutoff) always use default n_minimax 339 CALL get_minimax_coeff_v_gspace(n_minimax, cutoff, G_min, minimax_aw, err_minimax=err_mm, & 340 potential=potential, pot_par=pot_par) 341 CPASSERT(err_mm .LT. 1.1_dp*target_error + 1.0E-12) 342 CALL create_minimax_grid(minimax_grids(n_grids), cutoff, n_minimax, minimax_aw, err_mm) 343 DEALLOCATE (minimax_aw) 344 345 DO i_grid = n_grids - 1, 1, -1 346 DO n_mm = n_minimax, 1, -1 347 ALLOCATE (minimax_aw(2*n_mm)) 348 CALL get_minimax_coeff_v_gspace(n_mm, cutoff, G_min, minimax_aw, err_minimax=err_mm, & 349 potential=potential, pot_par=pot_par) 350 351 IF (err_mm .GT. 1.1_dp*target_error) THEN 352 CPASSERT(n_mm .NE. n_minimax) 353 CALL create_minimax_grid(minimax_grids(i_grid), cutoff, n_mm + 1, minimax_aw_prev, err_mm_prev) 354 355 DEALLOCATE (minimax_aw) 356 EXIT 357 ENDIF 358 359 IF (ALLOCATED(minimax_aw_prev)) DEALLOCATE (minimax_aw_prev) 360 ALLOCATE (minimax_aw_prev(2*n_mm)) 361 minimax_aw_prev(:) = minimax_aw(:) 362 DEALLOCATE (minimax_aw) 363 err_mm_prev = err_mm 364 ENDDO 365 cutoff = cutoff/cutoff_delta 366 ENDDO 367 END SUBROUTINE 368 369! ************************************************************************************************** 370!> \brief ... 371!> \param minimax_grids ... 372! ************************************************************************************************** 373 SUBROUTINE eri_mme_destroy_minimax_grids(minimax_grids) 374 TYPE(minimax_grid), ALLOCATABLE, DIMENSION(:), & 375 INTENT(INOUT) :: minimax_grids 376 377 INTEGER :: igrid 378 379 IF (ALLOCATED(minimax_grids)) THEN 380 DO igrid = 1, SIZE(minimax_grids) 381 IF (ASSOCIATED(minimax_grids(igrid)%minimax_aw)) THEN 382 DEALLOCATE (minimax_grids(igrid)%minimax_aw) 383 ENDIF 384 ENDDO 385 DEALLOCATE (minimax_grids) 386 ENDIF 387 END SUBROUTINE 388 389! ************************************************************************************************** 390!> \brief ... 391!> \param grid ... 392!> \param cutoff ... 393!> \param n_minimax ... 394!> \param minimax_aw ... 395!> \param error ... 396! ************************************************************************************************** 397 SUBROUTINE create_minimax_grid(grid, cutoff, n_minimax, minimax_aw, error) 398 TYPE(minimax_grid), INTENT(OUT) :: grid 399 REAL(KIND=dp), INTENT(IN) :: cutoff 400 INTEGER, INTENT(IN) :: n_minimax 401 REAL(KIND=dp), DIMENSION(2*n_minimax), INTENT(IN) :: minimax_aw 402 REAL(KIND=dp), INTENT(IN) :: error 403 404 grid%cutoff = cutoff 405 grid%n_minimax = n_minimax 406 ALLOCATE (grid%minimax_aw(2*n_minimax)) 407 grid%minimax_aw(:) = minimax_aw(:) 408 grid%error = error 409 410 END SUBROUTINE 411 412! ************************************************************************************************** 413!> \brief ... 414!> \param grids ... 415!> \param cutoff ... 416!> \param n_minimax ... 417!> \param minimax_aw ... 418!> \param grid_no ... 419! ************************************************************************************************** 420 SUBROUTINE get_minimax_from_cutoff(grids, cutoff, n_minimax, minimax_aw, grid_no) 421 TYPE(minimax_grid), DIMENSION(:), INTENT(IN) :: grids 422 REAL(KIND=dp), INTENT(IN) :: cutoff 423 INTEGER, INTENT(OUT) :: n_minimax 424 REAL(KIND=dp), DIMENSION(:), INTENT(OUT), POINTER :: minimax_aw 425 INTEGER, INTENT(OUT) :: grid_no 426 427 INTEGER :: igrid 428 429 grid_no = 0 430 DO igrid = 1, SIZE(grids) 431 IF (grids(igrid)%cutoff .GE. cutoff/2) THEN 432 n_minimax = grids(igrid)%n_minimax 433 minimax_aw => grids(igrid)%minimax_aw 434 grid_no = igrid 435 EXIT 436 ENDIF 437 ENDDO 438 IF (grid_no == 0) THEN 439 igrid = SIZE(grids) 440 n_minimax = grids(igrid)%n_minimax 441 minimax_aw => grids(igrid)%minimax_aw 442 grid_no = igrid 443 ENDIF 444 END SUBROUTINE 445 446! ************************************************************************************************** 447!> \brief ... 448!> \param grid ... 449!> \param grid_no ... 450!> \param unit_nr ... 451! ************************************************************************************************** 452 SUBROUTINE eri_mme_print_grid_info(grid, grid_no, unit_nr) 453 TYPE(minimax_grid), INTENT(IN) :: grid 454 INTEGER, INTENT(IN) :: grid_no, unit_nr 455 456 IF (unit_nr > 0) THEN 457 WRITE (unit_nr, '(T2, A, 1X, I2)') "ERI_MME | Info for grid no.", grid_no 458 WRITE (unit_nr, '(T2, A, 1X, ES9.2)') "ERI_MME | Cutoff", grid%cutoff 459 WRITE (unit_nr, '(T2, A, 1X, I2)') "ERI_MME | Number of minimax points", grid%n_minimax 460 WRITE (unit_nr, '(T2, A, 1X, 2ES9.2)') "ERI_MME | minimax error", grid%error 461 WRITE (unit_nr, *) 462 ENDIF 463 464 END SUBROUTINE 465 466! ************************************************************************************************** 467!> \brief ... 468!> \param param ... 469! ************************************************************************************************** 470 SUBROUTINE eri_mme_release(param) 471 TYPE(eri_mme_param), INTENT(INOUT) :: param 472 473 IF (ALLOCATED(param%minimax_grid)) THEN 474 CALL eri_mme_destroy_minimax_grids(param%minimax_grid) 475 ENDIF 476 END SUBROUTINE eri_mme_release 477 478END MODULE eri_mme_types 479