1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Several screening methods used in HFX calcualtions 8!> \par History 9!> 04.2008 created [Manuel Guidon] 10!> \author Manuel Guidon 11! ************************************************************************************************** 12MODULE hfx_screening_methods 13 USE ao_util, ONLY: exp_radius_very_extended 14 USE basis_set_types, ONLY: gto_basis_set_type 15 USE hfx_libint_interface, ONLY: evaluate_eri_screen 16 USE hfx_types, ONLY: hfx_basis_type,& 17 hfx_p_kind,& 18 hfx_potential_type,& 19 hfx_screen_coeff_type,& 20 log_zero,& 21 powell_min_log 22 USE kinds, ONLY: dp,& 23 int_8 24 USE libint_wrapper, ONLY: cp_libint_t 25 USE machine, ONLY: default_output_unit 26 USE orbital_pointers, ONLY: ncoset 27 USE powell, ONLY: opt_state_type,& 28 powell_optimize 29 USE qs_environment_types, ONLY: get_qs_env,& 30 qs_environment_type 31 USE qs_kind_types, ONLY: get_qs_kind,& 32 qs_kind_type 33#include "./base/base_uses.f90" 34 35 IMPLICIT NONE 36 PRIVATE 37 38 PUBLIC :: update_pmax_mat, & 39 calc_screening_functions, & 40 calc_pair_dist_radii 41 42 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_screening_methods' 43 44!*** 45 46CONTAINS 47 48! ************************************************************************************************** 49!> \brief calculates max values of two-electron integrals in a quartet/shell 50!> w.r.t. different zetas using the library lib_int 51!> \param lib ... 52!> \param ra position 53!> \param rb position 54!> \param zeta zeta 55!> \param zetb zeta 56!> \param la_min angular momentum 57!> \param la_max angular momentum 58!> \param lb_min angular momentum 59!> \param lb_max angular momentum 60!> \param npgfa number of primitive cartesian gaussian in actual shell 61!> \param npgfb number of primitive cartesian gaussian in actual shell 62!> \param max_val schwarz screening value 63!> \param potential_parameter contains info for libint 64!> \param tmp_R_1 pgf_based screening coefficients 65!> \param rab2 Squared Distance of centers ab 66!> \param p_work ... 67!> \par History 68!> 03.2007 created [Manuel Guidon] 69!> 02.2009 refactored [Manuel Guidon] 70!> \author Manuel Guidon 71! ************************************************************************************************** 72 SUBROUTINE screen4(lib, ra, rb, zeta, zetb, & 73 la_min, la_max, lb_min, lb_max, & 74 npgfa, npgfb, & 75 max_val, potential_parameter, tmp_R_1, & 76 rab2, p_work) 77 78 TYPE(cp_libint_t) :: lib 79 REAL(dp), INTENT(IN) :: ra(3), rb(3) 80 REAL(dp), DIMENSION(:), INTENT(IN) :: zeta, zetb 81 INTEGER, INTENT(IN) :: la_min, la_max, lb_min, lb_max, npgfa, & 82 npgfb 83 REAL(dp), INTENT(INOUT) :: max_val 84 TYPE(hfx_potential_type) :: potential_parameter 85 TYPE(hfx_screen_coeff_type), DIMENSION(:, :), & 86 POINTER :: tmp_R_1 87 REAL(dp) :: rab2 88 REAL(dp), DIMENSION(:), POINTER :: p_work 89 90 INTEGER :: ipgf, jpgf, la, lb 91 REAL(dp) :: max_val_temp, R1 92 93 max_val_temp = max_val 94 DO ipgf = 1, npgfa 95 DO jpgf = 1, npgfb 96 R1 = MAX(0.0_dp, tmp_R_1(jpgf, ipgf)%x(1)*rab2 + tmp_R_1(jpgf, ipgf)%x(2)) 97 DO la = la_min, la_max 98 DO lb = lb_min, lb_max 99 !Build primitives 100 CALL evaluate_eri_screen(lib, ra, rb, ra, rb, & 101 zeta(ipgf), zetb(jpgf), zeta(ipgf), zetb(jpgf), & 102 la, lb, la, lb, & 103 max_val_temp, potential_parameter, R1, R1, p_work) 104 105 max_val = MAX(max_val, max_val_temp) 106 END DO !lb 107 END DO !la 108 END DO !jpgf 109 END DO !ipgf 110 END SUBROUTINE screen4 111 112! ************************************************************************************************** 113!> \brief updates the maximum of the density matrix in compressed form for screening purposes 114!> \param pmax_set buffer to store matrix 115!> \param map_atom_to_kind_atom ... 116!> \param set_offset ... 117!> \param atomic_block_offset ... 118!> \param pmax_atom ... 119!> \param full_density_alpha ... 120!> \param full_density_beta ... 121!> \param natom ... 122!> \param kind_of ... 123!> \param basis_parameter ... 124!> \param nkind ... 125!> \param is_assoc_atomic_block ... 126!> \par History 127!> 09.2007 created [Manuel Guidon] 128!> \author Manuel Guidon 129!> \note 130!> - updates for each pair of shells the maximum absolute value of p 131! ************************************************************************************************** 132 SUBROUTINE update_pmax_mat(pmax_set, map_atom_to_kind_atom, set_offset, atomic_block_offset, & 133 pmax_atom, full_density_alpha, full_density_beta, natom, & 134 kind_of, basis_parameter, & 135 nkind, is_assoc_atomic_block) 136 137 TYPE(hfx_p_kind), DIMENSION(:), POINTER :: pmax_set 138 INTEGER, DIMENSION(:), POINTER :: map_atom_to_kind_atom 139 INTEGER, DIMENSION(:, :, :, :), POINTER :: set_offset 140 INTEGER, DIMENSION(:, :), POINTER :: atomic_block_offset 141 REAL(dp), DIMENSION(:, :), POINTER :: pmax_atom, full_density_alpha, & 142 full_density_beta 143 INTEGER, INTENT(IN) :: natom 144 INTEGER :: kind_of(*) 145 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter 146 INTEGER :: nkind 147 INTEGER, DIMENSION(:, :) :: is_assoc_atomic_block 148 149 CHARACTER(LEN=*), PARAMETER :: routineN = 'update_pmax_mat', & 150 routineP = moduleN//':'//routineN 151 152 INTEGER :: act_atomic_block_offset, act_set_offset, handle, i, img, jatom, jkind, jset, & 153 katom, kind_kind_idx, kkind, kset, mb, mc, nimg, nsetb, nsetc 154 INTEGER, DIMENSION(:), POINTER :: nsgfb, nsgfc 155 REAL(dp) :: pmax_tmp 156 157 CALL timeset(routineN, handle) 158 159 DO i = 1, SIZE(pmax_set) 160 pmax_set(i)%p_kind = 0.0_dp 161 END DO 162 163 pmax_atom = log_zero 164 165 nimg = SIZE(full_density_alpha, 2) 166 167 DO jatom = 1, natom 168 jkind = kind_of(jatom) 169 nsetb = basis_parameter(jkind)%nset 170 nsgfb => basis_parameter(jkind)%nsgf 171 172 DO katom = 1, natom 173 IF (.NOT. is_assoc_atomic_block(jatom, katom) >= 1) CYCLE 174 kkind = kind_of(katom) 175 IF (kkind < jkind) CYCLE 176 nsetc = basis_parameter(kkind)%nset 177 nsgfc => basis_parameter(kkind)%nsgf 178 act_atomic_block_offset = atomic_block_offset(katom, jatom) 179 DO jset = 1, nsetb 180 DO kset = 1, nsetc 181 IF (katom >= jatom) THEN 182 pmax_tmp = 0.0_dp 183 act_set_offset = set_offset(kset, jset, kkind, jkind) 184 DO img = 1, nimg 185 i = act_set_offset + act_atomic_block_offset - 1 186 DO mc = 1, nsgfc(kset) 187 DO mb = 1, nsgfb(jset) 188 pmax_tmp = MAX(pmax_tmp, ABS(full_density_alpha(i, img))) 189 IF (ASSOCIATED(full_density_beta)) THEN 190 pmax_tmp = MAX(pmax_tmp, ABS(full_density_beta(i, img))) 191 END IF 192 i = i + 1 193 ENDDO 194 ENDDO 195 ENDDO 196 IF (pmax_tmp == 0.0_dp) THEN 197 pmax_tmp = log_zero 198 ELSE 199 pmax_tmp = LOG10(pmax_tmp) 200 END IF 201 kind_kind_idx = INT(get_1D_idx(jkind, kkind, INT(nkind, int_8))) 202 pmax_set(kind_kind_idx)%p_kind(jset, & 203 kset, & 204 map_atom_to_kind_atom(jatom), & 205 map_atom_to_kind_atom(katom)) = pmax_tmp 206 ELSE 207 pmax_tmp = 0.0_dp 208 act_set_offset = set_offset(jset, kset, jkind, kkind) 209 DO img = 1, nimg 210 DO mc = 1, nsgfc(kset) 211 i = act_set_offset + act_atomic_block_offset - 1 + mc - 1 212 DO mb = 1, nsgfb(jset) 213 pmax_tmp = MAX(pmax_tmp, ABS(full_density_alpha(i, img))) 214 IF (ASSOCIATED(full_density_beta)) THEN 215 pmax_tmp = MAX(pmax_tmp, ABS(full_density_beta(i, img))) 216 END IF 217 i = i + nsgfc(kset) 218 ENDDO 219 ENDDO 220 ENDDO 221 IF (pmax_tmp == 0.0_dp) THEN 222 pmax_tmp = log_zero 223 ELSE 224 pmax_tmp = LOG10(pmax_tmp) 225 END IF 226 kind_kind_idx = INT(get_1D_idx(jkind, kkind, INT(nkind, int_8))) 227 pmax_set(kind_kind_idx)%p_kind(jset, & 228 kset, & 229 map_atom_to_kind_atom(jatom), & 230 map_atom_to_kind_atom(katom)) = pmax_tmp 231 END IF 232 END DO 233 END DO 234 END DO 235 END DO 236 237 DO jatom = 1, natom 238 jkind = kind_of(jatom) 239 nsetb = basis_parameter(jkind)%nset 240 DO katom = 1, natom 241 IF (.NOT. is_assoc_atomic_block(jatom, katom) >= 1) CYCLE 242 kkind = kind_of(katom) 243 IF (kkind < jkind) CYCLE 244 nsetc = basis_parameter(kkind)%nset 245 pmax_tmp = log_zero 246 DO jset = 1, nsetb 247 DO kset = 1, nsetc 248 kind_kind_idx = INT(get_1D_idx(jkind, kkind, INT(nkind, int_8))) 249 pmax_tmp = MAX(pmax_set(kind_kind_idx)%p_kind(jset, & 250 kset, & 251 map_atom_to_kind_atom(jatom), & 252 map_atom_to_kind_atom(katom)), pmax_tmp) 253 END DO 254 END DO 255 pmax_atom(jatom, katom) = pmax_tmp 256 pmax_atom(katom, jatom) = pmax_tmp 257 END DO 258 END DO 259 260 CALL timestop(handle) 261 262 END SUBROUTINE update_pmax_mat 263 264! ************************************************************************************************** 265!> \brief calculates screening functions for schwarz screening 266!> \param qs_env qs_env 267!> \param basis_parameter ... 268!> \param lib structure to libint 269!> \param potential_parameter contains infos on potential 270!> \param coeffs_set set based coefficients 271!> \param coeffs_kind kind based coefficients 272!> \param coeffs_pgf pgf based coefficients 273!> \param radii_pgf coefficients for long-range screening 274!> \param max_set Maximum Number of basis set sets in the system 275!> \param max_pgf Maximum Number of basis set pgfs in the system 276!> \param n_threads ... 277!> \param i_thread Thread ID of current task 278!> \param p_work ... 279!> \par History 280!> 02.2009 created [Manuel Guidon] 281!> \author Manuel Guidon 282!> \note 283!> This routine calculates (ab|ab) for different distances Rab = |a-b| 284!> and uses the powell optimiztion routines in order to fit the results 285!> in the following form: 286!> 287!> (ab|ab) = (ab|ab)(Rab) = c2*Rab^2 + c0 288!> 289!> The missing linear term assures that the functions is monotonically 290!> decaying such that c2 can be used as upper bound when applying the 291!> Minimum Image Convention in the periodic case. Furthermore 292!> it seems to be a good choice to fit the logarithm of the (ab|ab) 293!> The fitting takes place at several levels: kind, set and pgf within 294!> the corresponding ranges of the prodiuct charge distributions 295!> Doing so, we only need arrays of size nkinds^2*2 instead of big 296!> screening matrices 297! ************************************************************************************************** 298 299 SUBROUTINE calc_screening_functions(qs_env, basis_parameter, lib, potential_parameter, & 300 coeffs_set, coeffs_kind, coeffs_pgf, radii_pgf, & 301 max_set, max_pgf, n_threads, i_thread, & 302 p_work) 303 TYPE(qs_environment_type), POINTER :: qs_env 304 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter 305 TYPE(cp_libint_t) :: lib 306 TYPE(hfx_potential_type) :: potential_parameter 307 TYPE(hfx_screen_coeff_type), & 308 DIMENSION(:, :, :, :), POINTER :: coeffs_set 309 TYPE(hfx_screen_coeff_type), DIMENSION(:, :), & 310 POINTER :: coeffs_kind 311 TYPE(hfx_screen_coeff_type), & 312 DIMENSION(:, :, :, :, :, :), POINTER :: coeffs_pgf, radii_pgf 313 INTEGER, INTENT(IN) :: max_set, max_pgf, n_threads, i_thread 314 REAL(dp), DIMENSION(:), POINTER :: p_work 315 316 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_screening_functions', & 317 routineP = moduleN//':'//routineN 318 319 INTEGER :: handle, i, ikind, ipgf, iset, jkind, & 320 jpgf, jset, la, lb, ncoa, ncob, nkind, & 321 nseta, nsetb, sgfa, sgfb 322 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 323 npgfb, nsgfa, nsgfb 324 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 325 REAL(dp) :: kind_radius_a, kind_radius_b, max_contraction_a, max_contraction_b, max_val, & 326 max_val_temp, R1, ra(3), radius, rb(3), x(2) 327 REAL(dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b 328 REAL(dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb 329 REAL(dp), SAVE :: DATA(2, 0:100) 330 TYPE(gto_basis_set_type), POINTER :: orb_basis 331 TYPE(hfx_screen_coeff_type), DIMENSION(:, :), & 332 POINTER :: tmp_R_1 333 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 334 TYPE(qs_kind_type), POINTER :: qs_kind 335 336!$OMP MASTER 337 CALL timeset(routineN, handle) 338!$OMP END MASTER 339 340 CALL get_qs_env(qs_env=qs_env, & 341 qs_kind_set=qs_kind_set) 342 343 nkind = SIZE(qs_kind_set, 1) 344 345!$OMP MASTER 346 ALLOCATE (coeffs_pgf(max_pgf, max_pgf, max_set, max_set, nkind, nkind)) 347 348 DO ikind = 1, nkind 349 DO jkind = 1, nkind 350 DO iset = 1, max_set 351 DO jset = 1, max_set 352 DO ipgf = 1, max_pgf 353 DO jpgf = 1, max_pgf 354 coeffs_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(:) = 0.0_dp 355 END DO 356 END DO 357 END DO 358 END DO 359 END DO 360 END DO 361 362!$OMP END MASTER 363!$OMP BARRIER 364 ra = 0.0_dp 365 rb = 0.0_dp 366 DO ikind = 1, nkind 367 NULLIFY (qs_kind, orb_basis) 368 qs_kind => qs_kind_set(ikind) 369 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis) 370 NULLIFY (la_max, la_min, npgfa, zeta) 371 372 la_max => basis_parameter(ikind)%lmax 373 la_min => basis_parameter(ikind)%lmin 374 npgfa => basis_parameter(ikind)%npgf 375 nseta = basis_parameter(ikind)%nset 376 zeta => basis_parameter(ikind)%zet 377 set_radius_a => basis_parameter(ikind)%set_radius 378 first_sgfa => basis_parameter(ikind)%first_sgf 379 sphi_a => basis_parameter(ikind)%sphi 380 nsgfa => basis_parameter(ikind)%nsgf 381 rpgfa => basis_parameter(ikind)%pgf_radius 382 383 DO jkind = 1, nkind 384 NULLIFY (qs_kind, orb_basis) 385 qs_kind => qs_kind_set(jkind) 386 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis) 387 NULLIFY (lb_max, lb_min, npgfb, zetb) 388 389 lb_max => basis_parameter(jkind)%lmax 390 lb_min => basis_parameter(jkind)%lmin 391 npgfb => basis_parameter(jkind)%npgf 392 nsetb = basis_parameter(jkind)%nset 393 zetb => basis_parameter(jkind)%zet 394 set_radius_b => basis_parameter(jkind)%set_radius 395 first_sgfb => basis_parameter(jkind)%first_sgf 396 sphi_b => basis_parameter(jkind)%sphi 397 nsgfb => basis_parameter(jkind)%nsgf 398 rpgfb => basis_parameter(jkind)%pgf_radius 399 400 DO iset = 1, nseta 401 ncoa = npgfa(iset)*ncoset(la_max(iset)) 402 sgfa = first_sgfa(1, iset) 403 max_contraction_a = MAXVAL((/(SUM(ABS(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)/)) 404 DO jset = 1, nsetb 405 ncob = npgfb(jset)*ncoset(lb_max(jset)) 406 sgfb = first_sgfb(1, jset) 407 max_contraction_b = MAXVAL((/(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/)) 408 radius = set_radius_a(iset) + set_radius_b(jset) 409 DO ipgf = 1, npgfa(iset) 410 DO jpgf = 1, npgfb(jset) 411 radius = rpgfa(ipgf, iset) + rpgfb(jpgf, jset) 412 DO i = i_thread, 100, n_threads 413 rb(1) = 0.0_dp + REAL(i, dp)*0.01_dp*radius 414 max_val = 0.0_dp 415 R1 = MAX(0.0_dp, radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(1)*rb(1)**2 + & 416 radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(2)) 417 DO la = la_min(iset), la_max(iset) 418 DO lb = lb_min(jset), lb_max(jset) 419 !Build primitives 420 max_val_temp = 0.0_dp 421 CALL evaluate_eri_screen(lib, ra, rb, ra, rb, & 422 zeta(ipgf, iset), zetb(jpgf, jset), zeta(ipgf, iset), zetb(jpgf, jset), & 423 la, lb, la, lb, & 424 max_val_temp, potential_parameter, R1, R1, p_work) 425 max_val = MAX(max_val, max_val_temp) 426 END DO !lb 427 END DO !la 428 max_val = SQRT(max_val) 429 max_val = max_val*max_contraction_a*max_contraction_b 430 DATA(1, i) = rb(1) 431 IF (max_val == 0.0_dp) THEN 432 DATA(2, i) = powell_min_log 433 ELSE 434 DATA(2, i) = LOG10((max_val)) 435 END IF 436 END DO 437!$OMP BARRIER 438!$OMP MASTER 439 CALL optimize_it(DATA, x, powell_min_log) 440 coeffs_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x = x 441!$OMP END MASTER 442!$OMP BARRIER 443 444 END DO 445 END DO 446 END DO 447 END DO 448 END DO 449 END DO 450 451!$OMP MASTER 452 ALLOCATE (coeffs_set(max_set, max_set, nkind, nkind)) 453 454 DO ikind = 1, nkind 455 DO jkind = 1, nkind 456 DO iset = 1, max_set 457 DO jset = 1, max_set 458 coeffs_set(jset, iset, jkind, ikind)%x(:) = 0.0_dp 459 END DO 460 END DO 461 END DO 462 END DO 463!$OMP END MASTER 464!$OMP BARRIER 465 ra = 0.0_dp 466 rb = 0.0_dp 467 DO ikind = 1, nkind 468 NULLIFY (qs_kind, orb_basis) 469 qs_kind => qs_kind_set(ikind) 470 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis) 471 NULLIFY (la_max, la_min, npgfa, zeta) 472! CALL get_gto_basis_set(gto_basis_set=orb_basis,& 473! lmax=la_max,& 474! lmin=la_min,& 475! npgf=npgfa,& 476! nset=nseta,& 477! zet=zeta,& 478! set_radius=set_radius_a,& 479! first_sgf=first_sgfa,& 480! sphi=sphi_a,& 481! nsgf_set=nsgfa) 482 la_max => basis_parameter(ikind)%lmax 483 la_min => basis_parameter(ikind)%lmin 484 npgfa => basis_parameter(ikind)%npgf 485 nseta = basis_parameter(ikind)%nset 486 zeta => basis_parameter(ikind)%zet 487 set_radius_a => basis_parameter(ikind)%set_radius 488 first_sgfa => basis_parameter(ikind)%first_sgf 489 sphi_a => basis_parameter(ikind)%sphi 490 nsgfa => basis_parameter(ikind)%nsgf 491 492 DO jkind = 1, nkind 493 NULLIFY (qs_kind, orb_basis) 494 qs_kind => qs_kind_set(jkind) 495 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis) 496 NULLIFY (lb_max, lb_min, npgfb, zetb) 497 498 lb_max => basis_parameter(jkind)%lmax 499 lb_min => basis_parameter(jkind)%lmin 500 npgfb => basis_parameter(jkind)%npgf 501 nsetb = basis_parameter(jkind)%nset 502 zetb => basis_parameter(jkind)%zet 503 set_radius_b => basis_parameter(jkind)%set_radius 504 first_sgfb => basis_parameter(jkind)%first_sgf 505 sphi_b => basis_parameter(jkind)%sphi 506 nsgfb => basis_parameter(jkind)%nsgf 507 508 DO iset = 1, nseta 509 ncoa = npgfa(iset)*ncoset(la_max(iset)) 510 sgfa = first_sgfa(1, iset) 511 max_contraction_a = MAXVAL((/(SUM(ABS(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)/)) 512 DO jset = 1, nsetb 513 ncob = npgfb(jset)*ncoset(lb_max(jset)) 514 sgfb = first_sgfb(1, jset) 515 max_contraction_b = MAXVAL((/(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/)) 516 radius = set_radius_a(iset) + set_radius_b(jset) 517 tmp_R_1 => radii_pgf(:, :, jset, iset, jkind, ikind) 518 DO i = i_thread, 100, n_threads 519 rb(1) = 0.0_dp + REAL(i, dp)*0.01_dp*radius 520 max_val = 0.0_dp 521 CALL screen4(lib, ra, rb, & 522 zeta(:, iset), zetb(:, jset), & 523 la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), & 524 npgfa(iset), npgfb(jset), & 525 max_val, potential_parameter, tmp_R_1, rb(1)**2, p_work) 526 max_val = SQRT(max_val) 527 max_val = max_val*max_contraction_a*max_contraction_b 528 DATA(1, i) = rb(1) 529 IF (max_val == 0.0_dp) THEN 530 DATA(2, i) = powell_min_log 531 ELSE 532 DATA(2, i) = LOG10((max_val)) 533 END IF 534 END DO 535!$OMP BARRIER 536!$OMP MASTER 537 CALL optimize_it(DATA, x, powell_min_log) 538 coeffs_set(jset, iset, jkind, ikind)%x = x 539!$OMP END MASTER 540!$OMP BARRIER 541 END DO 542 END DO 543 544 END DO 545 END DO 546 547 ! ** now kinds 548!$OMP MASTER 549 ALLOCATE (coeffs_kind(nkind, nkind)) 550 551 DO ikind = 1, nkind 552 DO jkind = 1, nkind 553 coeffs_kind(jkind, ikind)%x(:) = 0.0_dp 554 END DO 555 END DO 556!$OMP END MASTER 557 ra = 0.0_dp 558 rb = 0.0_dp 559 DO ikind = 1, nkind 560 NULLIFY (qs_kind, orb_basis) 561 qs_kind => qs_kind_set(ikind) 562 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis) 563 NULLIFY (la_max, la_min, npgfa, zeta) 564 565 la_max => basis_parameter(ikind)%lmax 566 la_min => basis_parameter(ikind)%lmin 567 npgfa => basis_parameter(ikind)%npgf 568 nseta = basis_parameter(ikind)%nset 569 zeta => basis_parameter(ikind)%zet 570 kind_radius_a = basis_parameter(ikind)%kind_radius 571 first_sgfa => basis_parameter(ikind)%first_sgf 572 sphi_a => basis_parameter(ikind)%sphi 573 nsgfa => basis_parameter(ikind)%nsgf 574 575 DO jkind = 1, nkind 576 NULLIFY (qs_kind, orb_basis) 577 qs_kind => qs_kind_set(jkind) 578 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis) 579 NULLIFY (lb_max, lb_min, npgfb, zetb) 580 581 lb_max => basis_parameter(jkind)%lmax 582 lb_min => basis_parameter(jkind)%lmin 583 npgfb => basis_parameter(jkind)%npgf 584 nsetb = basis_parameter(jkind)%nset 585 zetb => basis_parameter(jkind)%zet 586 kind_radius_b = basis_parameter(jkind)%kind_radius 587 first_sgfb => basis_parameter(jkind)%first_sgf 588 sphi_b => basis_parameter(jkind)%sphi 589 nsgfb => basis_parameter(jkind)%nsgf 590 591 radius = kind_radius_a + kind_radius_b 592 DO iset = 1, nseta 593 ncoa = npgfa(iset)*ncoset(la_max(iset)) 594 sgfa = first_sgfa(1, iset) 595 max_contraction_a = MAXVAL((/(SUM(ABS(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)/)) 596 DO jset = 1, nsetb 597 ncob = npgfb(jset)*ncoset(lb_max(jset)) 598 sgfb = first_sgfb(1, jset) 599 max_contraction_b = MAXVAL((/(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/)) 600 DO i = i_thread, 100, n_threads 601 rb(1) = 0.0_dp + REAL(i, dp)*0.01_dp*radius 602 max_val = 0.0_dp 603 tmp_R_1 => radii_pgf(:, :, jset, iset, jkind, ikind) 604 CALL screen4(lib, ra, rb, & 605 zeta(:, iset), zetb(:, jset), & 606 la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), & 607 npgfa(iset), npgfb(jset), & 608 max_val, potential_parameter, tmp_R_1, rb(1)**2, p_work) 609 DATA(1, i) = rb(1) 610 max_val = SQRT(max_val) 611 max_val = max_val*max_contraction_a*max_contraction_b 612 IF (max_val == 0.0_dp) THEN 613 DATA(2, i) = MAX(powell_min_log, DATA(2, i)) 614 ELSE 615 DATA(2, i) = MAX(LOG10(max_val), DATA(2, i)) 616 END IF 617 END DO 618 END DO 619 END DO 620!$OMP BARRIER 621!$OMP MASTER 622 CALL optimize_it(DATA, x, powell_min_log) 623 coeffs_kind(jkind, ikind)%x = x 624!$OMP END MASTER 625!$OMP BARRIER 626 END DO 627 END DO 628 629!$OMP MASTER 630 CALL timestop(handle) 631!$OMP END MASTER 632 633 END SUBROUTINE calc_screening_functions 634 635! ************************************************************************************************** 636!> \brief calculates radius functions for longrange screening 637!> \param qs_env qs_env 638!> \param basis_parameter ... 639!> \param radii_pgf pgf based coefficients 640!> \param max_set Maximum Number of basis set sets in the system 641!> \param max_pgf Maximum Number of basis set pgfs in the system 642!> \param eps_schwarz ... 643!> \param n_threads ... 644!> \param i_thread Thread ID of current task 645!> \par History 646!> 02.2009 created [Manuel Guidon] 647!> \author Manuel Guidon 648!> \note 649!> This routine calculates the pair-distribution radius of a product 650!> Gaussian and uses the powell optimiztion routines in order to fit 651!> the results in the following form: 652!> 653!> (ab| = (ab(Rab) = c2*Rab^2 + c0 654!> 655! ************************************************************************************************** 656 657 SUBROUTINE calc_pair_dist_radii(qs_env, basis_parameter, & 658 radii_pgf, max_set, max_pgf, eps_schwarz, & 659 n_threads, i_thread) 660 661 TYPE(qs_environment_type), POINTER :: qs_env 662 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter 663 TYPE(hfx_screen_coeff_type), & 664 DIMENSION(:, :, :, :, :, :), POINTER :: radii_pgf 665 INTEGER, INTENT(IN) :: max_set, max_pgf 666 REAL(dp) :: eps_schwarz 667 INTEGER, INTENT(IN) :: n_threads, i_thread 668 669 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_pair_dist_radii', & 670 routineP = moduleN//':'//routineN 671 672 INTEGER :: handle, i, ikind, ipgf, iset, jkind, & 673 jpgf, jset, la, lb, ncoa, ncob, nkind, & 674 nseta, nsetb, sgfa, sgfb 675 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, & 676 npgfb, nsgfa, nsgfb 677 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb 678 REAL(dp) :: cutoff, ff, max_contraction_a, max_contraction_b, prefactor, R1, R_max, ra(3), & 679 rab(3), rab2, radius, rap(3), rb(3), rp(3), x(2), zetp 680 REAL(dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb 681 REAL(dp), SAVE :: DATA(2, 0:100) 682 TYPE(gto_basis_set_type), POINTER :: orb_basis 683 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 684 TYPE(qs_kind_type), POINTER :: qs_kind 685 686!$OMP MASTER 687 CALL timeset(routineN, handle) 688!$OMP END MASTER 689 CALL get_qs_env(qs_env=qs_env, & 690 qs_kind_set=qs_kind_set) 691 692 nkind = SIZE(qs_kind_set, 1) 693 ra = 0.0_dp 694 rb = 0.0_dp 695!$OMP MASTER 696 ALLOCATE (radii_pgf(max_pgf, max_pgf, max_set, max_set, nkind, nkind)) 697 DO ikind = 1, nkind 698 DO jkind = 1, nkind 699 DO iset = 1, max_set 700 DO jset = 1, max_set 701 DO ipgf = 1, max_pgf 702 DO jpgf = 1, max_pgf 703 radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x(:) = 0.0_dp 704 END DO 705 END DO 706 END DO 707 END DO 708 END DO 709 END DO 710 711 DATA = 0.0_dp 712!$OMP END MASTER 713!$OMP BARRIER 714 715 DO ikind = 1, nkind 716 NULLIFY (qs_kind, orb_basis) 717 qs_kind => qs_kind_set(ikind) 718 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis) 719 NULLIFY (la_max, la_min, npgfa, zeta) 720 721 la_max => basis_parameter(ikind)%lmax 722 la_min => basis_parameter(ikind)%lmin 723 npgfa => basis_parameter(ikind)%npgf 724 nseta = basis_parameter(ikind)%nset 725 zeta => basis_parameter(ikind)%zet 726 first_sgfa => basis_parameter(ikind)%first_sgf 727 sphi_a => basis_parameter(ikind)%sphi 728 nsgfa => basis_parameter(ikind)%nsgf 729 rpgfa => basis_parameter(ikind)%pgf_radius 730 731 DO jkind = 1, nkind 732 NULLIFY (qs_kind, orb_basis) 733 qs_kind => qs_kind_set(jkind) 734 CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis) 735 NULLIFY (lb_max, lb_min, npgfb, zetb) 736 737 lb_max => basis_parameter(jkind)%lmax 738 lb_min => basis_parameter(jkind)%lmin 739 npgfb => basis_parameter(jkind)%npgf 740 nsetb = basis_parameter(jkind)%nset 741 zetb => basis_parameter(jkind)%zet 742 first_sgfb => basis_parameter(jkind)%first_sgf 743 sphi_b => basis_parameter(jkind)%sphi 744 nsgfb => basis_parameter(jkind)%nsgf 745 rpgfb => basis_parameter(jkind)%pgf_radius 746 747 DO iset = 1, nseta 748 ncoa = npgfa(iset)*ncoset(la_max(iset)) 749 sgfa = first_sgfa(1, iset) 750 max_contraction_a = MAXVAL((/(SUM(ABS(sphi_a(1:ncoa, i))), i=sgfa, sgfa + nsgfa(iset) - 1)/)) 751 DO jset = 1, nsetb 752 ncob = npgfb(jset)*ncoset(lb_max(jset)) 753 sgfb = first_sgfb(1, jset) 754 max_contraction_b = MAXVAL((/(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/)) 755 DO ipgf = 1, npgfa(iset) 756 DO jpgf = 1, npgfb(jset) 757 radius = rpgfa(ipgf, iset) + rpgfb(jpgf, jset) 758 DO i = i_thread, 100, n_threads 759 rb(1) = 0.0_dp + 0.01_dp*radius*i 760 R_max = 0.0_dp 761 DO la = la_min(iset), la_max(iset) 762 DO lb = lb_min(jset), lb_max(jset) 763 zetp = zeta(ipgf, iset) + zetb(jpgf, jset) 764 ff = zetb(jpgf, jset)/zetp 765 rab = 0.0_dp 766 rab(1) = rb(1) 767 rab2 = rb(1)**2 768 prefactor = EXP(-zeta(ipgf, iset)*ff*rab2) 769 rap(:) = ff*rab(:) 770 rp(:) = ra(:) + rap(:) 771 rb(:) = ra(:) + rab(:) 772 cutoff = 1.0_dp 773 R1 = exp_radius_very_extended( & 774 la, la, lb, lb, ra=ra, rb=rb, rp=rp, & 775 zetp=zetp, eps=eps_schwarz, prefactor=prefactor, cutoff=cutoff, epsin=1.0E-12_dp) 776 R_max = MAX(R_max, R1) 777 END DO 778 END DO 779 DATA(1, i) = rb(1) 780 DATA(2, i) = R_max 781 END DO 782 ! the radius can not be negative, we take that into account in the code as well by using a MAX 783 ! the functional form we use for fitting does not seem particularly accurate 784!$OMP BARRIER 785!$OMP MASTER 786 CALL optimize_it(DATA, x, 0.0_dp) 787 radii_pgf(jpgf, ipgf, jset, iset, jkind, ikind)%x = x 788!$OMP END MASTER 789!$OMP BARRIER 790 END DO !jpgf 791 END DO !ipgf 792 END DO 793 END DO 794 END DO 795 END DO 796!$OMP MASTER 797 CALL timestop(handle) 798!$OMP END MASTER 799 END SUBROUTINE calc_pair_dist_radii 800 801! 802! 803! little driver routine for the powell minimizer 804! data is the data to fit, x is of the form (x(1)*DATA(1)**2+x(2)) 805! only values of DATA(2) larger than fmin are taken into account 806! it constructs an approximate upper bound of the fitted function 807! 808! 809! ************************************************************************************************** 810!> \brief ... 811!> \param DATA ... 812!> \param x ... 813!> \param fmin ... 814! ************************************************************************************************** 815 SUBROUTINE optimize_it(DATA, x, fmin) 816 817 REAL(KIND=dp), INTENT(IN) :: DATA(2, 0:100) 818 REAL(KIND=dp), INTENT(OUT) :: x(2) 819 REAL(KIND=dp), INTENT(IN) :: fmin 820 821 INTEGER :: i, k 822 REAL(KIND=dp) :: f, large_weight, small_weight, weight 823 TYPE(opt_state_type) :: opt_state 824 825! initial values 826 827 x(1) = 0.0_dp 828 x(2) = 0.0_dp 829 830 ! we go in two steps, first we do the symmetric weight to get a good, unique initial guess 831 ! we restart afterwards for the assym. 832 ! the assym function appears to have several local minima, depending on the data to fit 833 ! the loop over k can make the switch gradual, but there is not much need, seemingly 834 DO k = 0, 4, 4 835 836 small_weight = 1.0_dp 837 large_weight = small_weight*(10.0_dp**k) 838 839 ! init opt run 840 opt_state%state = 0 841 opt_state%nvar = 2 842 opt_state%iprint = 3 843 opt_state%unit = default_output_unit 844 opt_state%maxfun = 100000 845 opt_state%rhobeg = 0.1_dp 846 opt_state%rhoend = 0.000001_dp 847 848 DO 849 850 ! compute function 851 IF (opt_state%state == 2) THEN 852 opt_state%f = 0.0_dp 853 DO i = 0, 100 854 f = x(1)*DATA(1, i)**2 + x(2) 855 IF (f > DATA(2, i)) THEN 856 weight = small_weight 857 ELSE 858 weight = large_weight 859 END IF 860 IF (DATA(2, i) > fmin) opt_state%f = opt_state%f + weight*(f - DATA(2, i))**2 861 END DO 862 END IF 863 864 IF (opt_state%state == -1) EXIT 865 CALL powell_optimize(opt_state%nvar, x, opt_state) 866 END DO 867 868 ! dealloc mem 869 opt_state%state = 8 870 CALL powell_optimize(opt_state%nvar, x, opt_state) 871 872 ENDDO 873 874 END SUBROUTINE optimize_it 875 876! ************************************************************************************************** 877!> \brief Given a 2d index pair, this function returns a 1d index pair for 878!> a symmetric upper triangle NxN matrix 879!> The compiler should inline this function, therefore it appears in 880!> several modules 881!> \param i 2d index 882!> \param j 2d index 883!> \param N matrix size 884!> \return ... 885!> \par History 886!> 03.2009 created [Manuel Guidon] 887!> \author Manuel Guidon 888! ************************************************************************************************** 889 PURE FUNCTION get_1D_idx(i, j, N) 890 INTEGER, INTENT(IN) :: i, j 891 INTEGER(int_8), INTENT(IN) :: N 892 INTEGER(int_8) :: get_1D_idx 893 894 INTEGER(int_8) :: min_ij 895 896 min_ij = MIN(i, j) 897 get_1D_idx = min_ij*N + MAX(i, j) - (min_ij - 1)*min_ij/2 - N 898 899 END FUNCTION get_1D_idx 900 901END MODULE hfx_screening_methods 902