1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Rountines to calculate the 3 and 2 center ERI's needed in the RI 8!> approximation using libint 9!> \par History 10!> 08.2013 created [Mauro Del Ben] 11!> \author Mauro Del Ben 12! ************************************************************************************************** 13MODULE mp2_ri_libint 14 USE atomic_kind_types, ONLY: atomic_kind_type,& 15 get_atomic_kind_set 16 USE basis_set_types, ONLY: get_gto_basis_set,& 17 gto_basis_set_type,& 18 init_aux_basis_set 19 USE cell_types, ONLY: cell_type 20 USE cp_blacs_env, ONLY: cp_blacs_env_create,& 21 cp_blacs_env_release,& 22 cp_blacs_env_type 23 USE cp_control_types, ONLY: dft_control_type 24 USE cp_files, ONLY: close_file,& 25 open_file 26 USE cp_fm_basic_linalg, ONLY: cp_fm_triangular_invert 27 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose 28 USE cp_fm_diag, ONLY: choose_eigv_solver 29 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 30 cp_fm_struct_release,& 31 cp_fm_struct_type 32 USE cp_fm_types, ONLY: cp_fm_create,& 33 cp_fm_get_submatrix,& 34 cp_fm_release,& 35 cp_fm_set_all,& 36 cp_fm_set_submatrix,& 37 cp_fm_to_fm,& 38 cp_fm_type 39 USE cp_para_types, ONLY: cp_para_env_type 40 USE gamma, ONLY: init_md_ftable 41 USE hfx_energy_potential, ONLY: coulomb4 42 USE hfx_pair_list_methods, ONLY: build_pair_list_mp2 43 USE hfx_screening_methods, ONLY: calc_pair_dist_radii,& 44 calc_screening_functions 45 USE hfx_types, ONLY: & 46 hfx_basis_info_type, hfx_basis_type, hfx_create_neighbor_cells, hfx_pgf_list, & 47 hfx_pgf_product_list, hfx_potential_type, hfx_screen_coeff_type, hfx_type, log_zero, & 48 pair_set_list_type 49 USE input_constants, ONLY: do_potential_TShPSC,& 50 do_potential_long 51 USE kinds, ONLY: dp,& 52 int_8 53 USE libint_wrapper, ONLY: cp_libint_t 54 USE message_passing, ONLY: mp_sum 55 USE mp2_types, ONLY: init_TShPSC_lmax,& 56 mp2_biel_type,& 57 mp2_type,& 58 pair_list_type_mp2 59 USE orbital_pointers, ONLY: nco,& 60 ncoset,& 61 nso 62 USE particle_types, ONLY: particle_type 63 USE qs_environment_types, ONLY: get_qs_env,& 64 qs_environment_type 65 USE qs_kind_types, ONLY: get_qs_kind,& 66 get_qs_kind_set,& 67 qs_kind_type 68 USE t_sh_p_s_c, ONLY: free_C0,& 69 init 70#include "./base/base_uses.f90" 71 72 IMPLICIT NONE 73 PRIVATE 74 75 PUBLIC :: libint_ri_mp2, read_RI_basis_set, release_RI_basis_set, prepare_integral_calc 76 77 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_ri_libint' 78 79!*** 80 81CONTAINS 82 83! ************************************************************************************************** 84!> \brief ... 85!> \param dimen ... 86!> \param RI_dimen ... 87!> \param occupied ... 88!> \param natom ... 89!> \param mp2_biel ... 90!> \param mp2_env ... 91!> \param C ... 92!> \param kind_of ... 93!> \param RI_basis_parameter ... 94!> \param RI_basis_info ... 95!> \param basis_S0 ... 96!> \param RI_index_table ... 97!> \param qs_env ... 98!> \param para_env ... 99!> \param Lai ... 100! ************************************************************************************************** 101 SUBROUTINE libint_ri_mp2(dimen, RI_dimen, occupied, natom, mp2_biel, mp2_env, C, & 102 kind_of, & 103 RI_basis_parameter, RI_basis_info, basis_S0, RI_index_table, & 104 qs_env, para_env, & 105 Lai) 106 INTEGER :: dimen, RI_dimen, occupied, natom 107 TYPE(mp2_biel_type) :: mp2_biel 108 TYPE(mp2_type), POINTER :: mp2_env 109 REAL(KIND=dp), DIMENSION(dimen, dimen) :: C 110 INTEGER, DIMENSION(:) :: kind_of 111 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: RI_basis_parameter 112 TYPE(hfx_basis_info_type) :: RI_basis_info 113 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_S0 114 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: RI_index_table 115 TYPE(qs_environment_type), POINTER :: qs_env 116 TYPE(cp_para_env_type), POINTER :: para_env 117 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Lai 118 119 CHARACTER(LEN=*), PARAMETER :: routineN = 'libint_ri_mp2', routineP = moduleN//':'//routineN 120 121 INTEGER :: handle, nkind 122 123 CALL timeset(routineN, handle) 124 125 ! Get the RI basis set and store in to a nice form 126 IF (.NOT. (ASSOCIATED(RI_basis_parameter))) THEN 127 IF (ALLOCATED(RI_index_table)) DEALLOCATE (RI_index_table) 128 IF (ASSOCIATED(basis_S0)) DEALLOCATE (basis_S0) 129 CALL read_RI_basis_set(qs_env, RI_basis_parameter, RI_basis_info, & 130 natom, nkind, kind_of, RI_index_table, RI_dimen, & 131 basis_S0) 132 END IF 133 134 CALL calc_lai_libint(mp2_env, qs_env, para_env, & 135 mp2_biel, dimen, C, occupied, & 136 RI_basis_parameter, RI_basis_info, RI_index_table, RI_dimen, basis_S0, & 137 Lai) 138 139 CALL timestop(handle) 140 141 END SUBROUTINE libint_ri_mp2 142 143! ************************************************************************************************** 144!> \brief Read the auxiliary basis set for RI approxiamtion 145!> \param qs_env ... 146!> \param RI_basis_parameter ... 147!> \param RI_basis_info ... 148!> \param natom ... 149!> \param nkind ... 150!> \param kind_of ... 151!> \param RI_index_table ... 152!> \param RI_dimen ... 153!> \param basis_S0 ... 154!> \par History 155!> 08.2013 created [Mauro Del Ben] 156!> \author Mauro Del Ben 157! ************************************************************************************************** 158 SUBROUTINE read_RI_basis_set(qs_env, RI_basis_parameter, RI_basis_info, & 159 natom, nkind, kind_of, RI_index_table, RI_dimen, & 160 basis_S0) 161 TYPE(qs_environment_type), POINTER :: qs_env 162 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: RI_basis_parameter 163 TYPE(hfx_basis_info_type) :: RI_basis_info 164 INTEGER :: natom, nkind 165 INTEGER, DIMENSION(:) :: kind_of 166 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: RI_index_table 167 INTEGER :: RI_dimen 168 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_S0 169 170 CHARACTER(LEN=*), PARAMETER :: routineN = 'read_RI_basis_set', & 171 routineP = moduleN//':'//routineN 172 173 INTEGER :: co_counter, counter, i, iatom, ikind, ipgf, iset, j, k, la, max_am_kind, & 174 max_coeff, max_nsgfl, max_pgf, max_pgf_kind, max_set, nl_count, nset, nseta, offset_a, & 175 offset_a1, s_offset_nl_a, sgfa, so_counter 176 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nshell 177 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, nl_a 178 REAL(dp), DIMENSION(:, :), POINTER :: sphi_a 179 TYPE(gto_basis_set_type), POINTER :: orb_basis_a 180 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 181 TYPE(qs_kind_type), POINTER :: atom_kind 182 183 NULLIFY (RI_basis_parameter) 184 185 NULLIFY (qs_kind_set) 186 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set) 187 188 nkind = SIZE(qs_kind_set, 1) 189 ALLOCATE (RI_basis_parameter(nkind)) 190 ALLOCATE (basis_S0(nkind)) 191 max_set = 0 192 DO ikind = 1, nkind 193 NULLIFY (atom_kind) 194 atom_kind => qs_kind_set(ikind) 195 ! here we reset the initial RI basis such that we can 196 ! work with non-normalized auxiliary basis functions 197 CALL get_qs_kind(qs_kind=atom_kind, & 198 basis_set=orb_basis_a, basis_type="RI_AUX") 199 IF (.NOT. (ASSOCIATED(orb_basis_a))) THEN 200 CPABORT("Initial RI auxiliary basis not specified.") 201 END IF 202 orb_basis_a%gcc = 1.0_dp 203 orb_basis_a%norm_type = 1 204 CALL init_aux_basis_set(orb_basis_a) 205 206 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, & 207 maxsgf=RI_basis_info%max_sgf, & 208 maxnset=RI_basis_info%max_set, & 209 maxlgto=RI_basis_info%max_am, & 210 basis_type="RI_AUX") 211 CALL get_gto_basis_set(gto_basis_set=orb_basis_a, & 212 lmax=RI_basis_parameter(ikind)%lmax, & 213 lmin=RI_basis_parameter(ikind)%lmin, & 214 npgf=RI_basis_parameter(ikind)%npgf, & 215 nset=RI_basis_parameter(ikind)%nset, & 216 zet=RI_basis_parameter(ikind)%zet, & 217 nsgf_set=RI_basis_parameter(ikind)%nsgf, & 218 first_sgf=RI_basis_parameter(ikind)%first_sgf, & 219 sphi=RI_basis_parameter(ikind)%sphi, & 220 nsgf=RI_basis_parameter(ikind)%nsgf_total, & 221 l=RI_basis_parameter(ikind)%nl, & 222 nshell=RI_basis_parameter(ikind)%nshell, & 223 set_radius=RI_basis_parameter(ikind)%set_radius, & 224 pgf_radius=RI_basis_parameter(ikind)%pgf_radius, & 225 kind_radius=RI_basis_parameter(ikind)%kind_radius) 226 227 max_set = MAX(max_set, RI_basis_parameter(ikind)%nset) 228 229 basis_S0(ikind)%kind_radius = RI_basis_parameter(ikind)%kind_radius 230 basis_S0(ikind)%nset = 1 231 basis_S0(ikind)%nsgf_total = 1 232 ALLOCATE (basis_S0(ikind)%set_radius(1)) 233 basis_S0(ikind)%set_radius(1) = RI_basis_parameter(ikind)%kind_radius 234 ALLOCATE (basis_S0(ikind)%lmax(1)) 235 basis_S0(ikind)%lmax(1) = 0 236 ALLOCATE (basis_S0(ikind)%lmin(1)) 237 basis_S0(ikind)%lmin(1) = 0 238 ALLOCATE (basis_S0(ikind)%npgf(1)) 239 basis_S0(ikind)%npgf(1) = 1 240 ALLOCATE (basis_S0(ikind)%nsgf(1)) 241 basis_S0(ikind)%nsgf(1) = 1 242 ALLOCATE (basis_S0(ikind)%nshell(1)) 243 basis_S0(ikind)%nshell(1) = 1 244 ALLOCATE (basis_S0(ikind)%pgf_radius(1, 1)) 245 basis_S0(ikind)%pgf_radius(1, 1) = RI_basis_parameter(ikind)%kind_radius 246 ALLOCATE (basis_S0(ikind)%sphi(1, 1)) 247 basis_S0(ikind)%sphi(1, 1) = 1.0_dp 248 ALLOCATE (basis_S0(ikind)%zet(1, 1)) 249 basis_S0(ikind)%zet(1, 1) = 0.0_dp 250 ALLOCATE (basis_S0(ikind)%first_sgf(1, 1)) 251 basis_S0(ikind)%first_sgf(1, 1) = 1 252 ALLOCATE (basis_S0(ikind)%nl(1, 1)) 253 basis_S0(ikind)%nl(1, 1) = 0 254 255 ALLOCATE (basis_S0(ikind)%nsgfl(0:0, 1)) 256 basis_S0(ikind)%nsgfl = 1 257 ALLOCATE (basis_S0(ikind)%sphi_ext(1, 0:0, 1, 1)) 258 basis_S0(ikind)%sphi_ext(1, 0, 1, 1) = 1.0_dp 259 260 END DO 261 RI_basis_info%max_set = max_set 262 263 DO ikind = 1, nkind 264 ALLOCATE (RI_basis_parameter(ikind)%nsgfl(0:RI_basis_info%max_am, max_set)) 265 RI_basis_parameter(ikind)%nsgfl = 0 266 nset = RI_basis_parameter(ikind)%nset 267 nshell => RI_basis_parameter(ikind)%nshell 268 DO iset = 1, nset 269 DO i = 0, RI_basis_info%max_am 270 nl_count = 0 271 DO j = 1, nshell(iset) 272 IF (RI_basis_parameter(ikind)%nl(j, iset) == i) nl_count = nl_count + 1 273 END DO 274 RI_basis_parameter(ikind)%nsgfl(i, iset) = nl_count 275 END DO 276 END DO 277 END DO 278 279 max_nsgfl = 0 280 max_pgf = 0 281 DO ikind = 1, nkind 282 max_coeff = 0 283 max_am_kind = 0 284 max_pgf_kind = 0 285 npgfa => RI_basis_parameter(ikind)%npgf 286 nseta = RI_basis_parameter(ikind)%nset 287 nl_a => RI_basis_parameter(ikind)%nsgfl 288 la_max => RI_basis_parameter(ikind)%lmax 289 la_min => RI_basis_parameter(ikind)%lmin 290 DO iset = 1, nseta 291 max_pgf_kind = MAX(max_pgf_kind, npgfa(iset)) 292 max_pgf = MAX(max_pgf, npgfa(iset)) 293 DO la = la_min(iset), la_max(iset) 294 max_nsgfl = MAX(max_nsgfl, nl_a(la, iset)) 295 max_coeff = MAX(max_coeff, nso(la)*nl_a(la, iset)*nco(la)) 296 max_am_kind = MAX(max_am_kind, la) 297 END DO 298 END DO 299 ALLOCATE (RI_basis_parameter(ikind)%sphi_ext(max_coeff, 0:max_am_kind, max_pgf_kind, nseta)) 300 RI_basis_parameter(ikind)%sphi_ext = 0.0_dp 301 END DO 302 303 DO ikind = 1, nkind 304 sphi_a => RI_basis_parameter(ikind)%sphi 305 nseta = RI_basis_parameter(ikind)%nset 306 la_max => RI_basis_parameter(ikind)%lmax 307 la_min => RI_basis_parameter(ikind)%lmin 308 npgfa => RI_basis_parameter(ikind)%npgf 309 first_sgfa => RI_basis_parameter(ikind)%first_sgf 310 nl_a => RI_basis_parameter(ikind)%nsgfl 311 DO iset = 1, nseta 312 sgfa = first_sgfa(1, iset) 313 DO ipgf = 1, npgfa(iset) 314 offset_a1 = (ipgf - 1)*ncoset(la_max(iset)) 315 s_offset_nl_a = 0 316 DO la = la_min(iset), la_max(iset) 317 offset_a = offset_a1 + ncoset(la - 1) 318 co_counter = 0 319 co_counter = co_counter + 1 320 so_counter = 0 321 DO k = sgfa + s_offset_nl_a, sgfa + s_offset_nl_a + nso(la)*nl_a(la, iset) - 1 322 DO i = offset_a + 1, offset_a + nco(la) 323 so_counter = so_counter + 1 324 RI_basis_parameter(ikind)%sphi_ext(so_counter, la, ipgf, iset) = sphi_a(i, k) 325 END DO 326 END DO 327 s_offset_nl_a = s_offset_nl_a + nso(la)*(nl_a(la, iset)) 328 END DO 329 END DO 330 END DO 331 END DO 332 333 ALLOCATE (RI_index_table(natom, max_set)) 334 RI_index_table = -HUGE(0) 335 counter = 0 336 RI_dimen = 0 337 DO iatom = 1, natom 338 ikind = kind_of(iatom) 339 nset = RI_basis_parameter(ikind)%nset 340 DO iset = 1, nset 341 RI_index_table(iatom, iset) = counter + 1 342 counter = counter + RI_basis_parameter(ikind)%nsgf(iset) 343 RI_dimen = RI_dimen + RI_basis_parameter(ikind)%nsgf(iset) 344 END DO 345 END DO 346 347 END SUBROUTINE read_RI_basis_set 348 349! ************************************************************************************************** 350!> \brief Release the auxiliary basis set for RI approxiamtion (to be used 351!> only in the case of basis optimization) 352!> \param RI_basis_parameter ... 353!> \param basis_S0 ... 354!> \par History 355!> 08.2013 created [Mauro Del Ben] 356!> \author Mauro Del Ben 357! ************************************************************************************************** 358 SUBROUTINE release_RI_basis_set(RI_basis_parameter, basis_S0) 359 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: RI_basis_parameter, basis_S0 360 361 CHARACTER(LEN=*), PARAMETER :: routineN = 'release_RI_basis_set', & 362 routineP = moduleN//':'//routineN 363 364 INTEGER :: i 365 366! RI basis 367 368 DO i = 1, SIZE(RI_basis_parameter) 369 DEALLOCATE (RI_basis_parameter(i)%nsgfl) 370 DEALLOCATE (RI_basis_parameter(i)%sphi_ext) 371 END DO 372 DEALLOCATE (RI_basis_parameter) 373 374 ! S0 basis 375 DO i = 1, SIZE(basis_S0) 376 DEALLOCATE (basis_S0(i)%set_radius) 377 DEALLOCATE (basis_S0(i)%lmax) 378 DEALLOCATE (basis_S0(i)%lmin) 379 DEALLOCATE (basis_S0(i)%npgf) 380 DEALLOCATE (basis_S0(i)%nsgf) 381 DEALLOCATE (basis_S0(i)%nshell) 382 DEALLOCATE (basis_S0(i)%pgf_radius) 383 DEALLOCATE (basis_S0(i)%sphi) 384 DEALLOCATE (basis_S0(i)%zet) 385 DEALLOCATE (basis_S0(i)%first_sgf) 386 DEALLOCATE (basis_S0(i)%nl) 387 DEALLOCATE (basis_S0(i)%nsgfl) 388 DEALLOCATE (basis_S0(i)%sphi_ext) 389 END DO 390 DEALLOCATE (basis_S0) 391 392 END SUBROUTINE release_RI_basis_set 393 394! ************************************************************************************************** 395!> \brief ... 396!> \param mp2_env ... 397!> \param qs_env ... 398!> \param para_env ... 399!> \param mp2_biel ... 400!> \param dimen ... 401!> \param C ... 402!> \param occupied ... 403!> \param RI_basis_parameter ... 404!> \param RI_basis_info ... 405!> \param RI_index_table ... 406!> \param RI_dimen ... 407!> \param basis_S0 ... 408!> \param Lai ... 409!> \par History 410!> 08.2013 created [Mauro Del Ben] 411!> \author Mauro Del Ben 412! ************************************************************************************************** 413 SUBROUTINE calc_lai_libint(mp2_env, qs_env, para_env, & 414 mp2_biel, dimen, C, occupied, & 415 RI_basis_parameter, RI_basis_info, RI_index_table, RI_dimen, basis_S0, & 416 Lai) 417 418 TYPE(mp2_type), POINTER :: mp2_env 419 TYPE(qs_environment_type), POINTER :: qs_env 420 TYPE(cp_para_env_type), POINTER :: para_env 421 TYPE(mp2_biel_type) :: mp2_biel 422 INTEGER :: dimen 423 REAL(KIND=dp), DIMENSION(dimen, dimen) :: C 424 INTEGER :: occupied 425 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: RI_basis_parameter 426 TYPE(hfx_basis_info_type) :: RI_basis_info 427 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: RI_index_table 428 INTEGER :: RI_dimen 429 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_S0 430 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: Lai 431 432 CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_lai_libint', & 433 routineP = moduleN//':'//routineN 434 435 INTEGER :: counter_L_blocks, handle, i, i_list_kl, i_set_list_kl, i_set_list_kl_start, & 436 i_set_list_kl_stop, iatom, iatom_end, iatom_start, iiB, ikind, info_chol, iset, jatom, & 437 jatom_end, jatom_start, jjB, jkind, jset, katom, katom_end, katom_start, kkB, kkind, & 438 kset, kset_start, L_B_i_end, L_B_i_start, L_B_k_end, L_B_k_start, latom, latom_end, & 439 latom_start, lkind, llB, lset, max_pgf, max_set, natom, ncob, nints, nseta, nsetb, nsetc, & 440 nsetd, nsgf_max, nspins, orb_k_end, orb_k_start, orb_l_end, orb_l_start, & 441 primitive_counter, sphi_a_u1, sphi_a_u2, sphi_a_u3, sphi_b_u1, sphi_b_u2, sphi_b_u3 442 INTEGER :: sphi_c_u1, sphi_c_u2, sphi_c_u3, sphi_d_u1, sphi_d_u2, sphi_d_u3, virtual 443 INTEGER(int_8) :: estimate_to_store_int, neris_tmp, & 444 neris_total, nprim_ints 445 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, nimages 446 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, lc_max, & 447 lc_min, ld_max, ld_min, npgfa, npgfb, & 448 npgfc, npgfd, nsgfa, nsgfb, nsgfc, & 449 nsgfd 450 INTEGER, DIMENSION(:, :), POINTER :: nsgfl_a, nsgfl_b, nsgfl_c, nsgfl_d 451 LOGICAL :: do_periodic 452 LOGICAL, DIMENSION(:, :), POINTER :: shm_atomic_pair_list 453 REAL(KIND=dp) :: cartesian_estimate, coeffs_kind_max0, eps_schwarz, ln_10, & 454 log10_eps_schwarz, log10_pmax, max_contraction_val, pmax_atom, pmax_entry, ra(3), rb(3), & 455 rc(3), rd(3) 456 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ee_buffer1, ee_buffer2, & 457 ee_primitives_tmp, ee_work, ee_work2, & 458 primitive_integrals 459 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: L_block, L_full_matrix, max_contraction 460 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: BI1, MNRS 461 REAL(KIND=dp), DIMENSION(:), POINTER :: p_work 462 REAL(KIND=dp), DIMENSION(:, :), POINTER :: shm_pmax_block, zeta, zetb, zetc, zetd 463 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: sphi_a_ext_set, sphi_b_ext_set, & 464 sphi_c_ext_set, sphi_d_ext_set 465 REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: sphi_a_ext, sphi_b_ext, sphi_c_ext, & 466 sphi_d_ext 467 TYPE(cell_type), POINTER :: cell 468 TYPE(cp_blacs_env_type), POINTER :: blacs_env 469 TYPE(cp_fm_struct_type), POINTER :: fm_struct_L 470 TYPE(cp_fm_type), POINTER :: fm_matrix_L 471 TYPE(cp_libint_t) :: private_lib 472 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter 473 TYPE(hfx_pgf_list), ALLOCATABLE, DIMENSION(:) :: pgf_list_ij, pgf_list_kl 474 TYPE(hfx_pgf_product_list), ALLOCATABLE, & 475 DIMENSION(:) :: pgf_product_list 476 TYPE(hfx_potential_type) :: mp2_potential_parameter 477 TYPE(hfx_screen_coeff_type), ALLOCATABLE, & 478 DIMENSION(:, :), TARGET :: radii_pgf_large, screen_coeffs_pgf_large 479 TYPE(hfx_screen_coeff_type), DIMENSION(:, :), & 480 POINTER :: screen_coeffs_kind, tmp_R_1, tmp_R_2, & 481 tmp_screen_pgf1, tmp_screen_pgf2 482 TYPE(hfx_screen_coeff_type), & 483 DIMENSION(:, :, :, :), POINTER :: screen_coeffs_set 484 TYPE(hfx_screen_coeff_type), & 485 DIMENSION(:, :, :, :, :, :), POINTER :: radii_pgf, screen_coeffs_pgf 486 TYPE(hfx_type), POINTER :: actual_x_data 487 TYPE(pair_list_type_mp2) :: list_kl 488 TYPE(pair_set_list_type), ALLOCATABLE, & 489 DIMENSION(:) :: set_list_kl 490 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 491 492 CALL timeset(routineN, handle) 493 494 ln_10 = LOG(10.0_dp) 495 496 !! initalize some counters 497 neris_tmp = 0_int_8 498 neris_total = 0_int_8 499 nprim_ints = 0_int_8 500 501 CALL prepare_integral_calc(cell, qs_env, mp2_env, para_env, mp2_potential_parameter, actual_x_data, & 502 do_periodic, basis_parameter, max_set, particle_set, natom, kind_of, & 503 nsgf_max, primitive_integrals, ee_work, ee_work2, ee_buffer1, ee_buffer2, & 504 ee_primitives_tmp, nspins, max_contraction, max_pgf, pgf_list_ij, & 505 pgf_list_kl, pgf_product_list, nimages, eps_schwarz, log10_eps_schwarz, & 506 private_lib, p_work, screen_coeffs_set, screen_coeffs_kind, screen_coeffs_pgf, & 507 radii_pgf, RI_basis_parameter, RI_basis_info) 508 509 ALLOCATE (radii_pgf_large(SIZE(radii_pgf, 1), SIZE(radii_pgf, 2))) 510 ALLOCATE (screen_coeffs_pgf_large(SIZE(screen_coeffs_pgf, 1), SIZE(screen_coeffs_pgf, 2))) 511 DO iiB = 1, SIZE(radii_pgf, 1) 512 DO jjB = 1, SIZE(radii_pgf, 2) 513 radii_pgf_large(iiB, jjB)%x = 100_dp 514 END DO 515 END DO 516 DO iiB = 1, SIZE(screen_coeffs_pgf, 1) 517 DO jjB = 1, SIZE(screen_coeffs_pgf, 2) 518 screen_coeffs_pgf_large(iiB, jjB)%x = 5000_dp 519 END DO 520 END DO 521 tmp_R_1 => radii_pgf_large(:, :) 522 tmp_R_2 => radii_pgf_large(:, :) 523 tmp_screen_pgf1 => screen_coeffs_pgf_large(:, :) 524 tmp_screen_pgf2 => screen_coeffs_pgf_large(:, :) 525 526 ! start computing the L matrix 527 ALLOCATE (L_full_matrix(RI_dimen, RI_dimen)) 528 L_full_matrix = 0.0_dp 529 530 counter_L_blocks = 0 531 DO iatom = 1, natom 532 jatom = iatom 533 534 ikind = kind_of(iatom) 535 jkind = kind_of(jatom) 536 ra = particle_set(iatom)%r(:) 537 rb = particle_set(jatom)%r(:) 538 539 la_max => RI_basis_parameter(ikind)%lmax 540 la_min => RI_basis_parameter(ikind)%lmin 541 npgfa => RI_basis_parameter(ikind)%npgf 542 nseta = RI_basis_parameter(ikind)%nset 543 zeta => RI_basis_parameter(ikind)%zet 544 nsgfa => RI_basis_parameter(ikind)%nsgf 545 sphi_a_ext => RI_basis_parameter(ikind)%sphi_ext(:, :, :, :) 546 nsgfl_a => RI_basis_parameter(ikind)%nsgfl 547 sphi_a_u1 = UBOUND(sphi_a_ext, 1) 548 sphi_a_u2 = UBOUND(sphi_a_ext, 2) 549 sphi_a_u3 = UBOUND(sphi_a_ext, 3) 550 551 lb_max => basis_S0(jkind)%lmax 552 lb_min => basis_S0(jkind)%lmin 553 npgfb => basis_S0(jkind)%npgf 554 nsetb = basis_S0(jkind)%nset 555 zetb => basis_S0(jkind)%zet 556 nsgfb => basis_S0(jkind)%nsgf 557 sphi_b_ext => basis_S0(jkind)%sphi_ext(:, :, :, :) 558 nsgfl_b => basis_S0(jkind)%nsgfl 559 sphi_b_u1 = UBOUND(sphi_b_ext, 1) 560 sphi_b_u2 = UBOUND(sphi_b_ext, 2) 561 sphi_b_u3 = UBOUND(sphi_b_ext, 3) 562 563 DO katom = iatom, natom 564 latom = katom 565 566 kkind = kind_of(katom) 567 lkind = kind_of(latom) 568 rc = particle_set(katom)%r(:) 569 rd = particle_set(latom)%r(:) 570 571 lc_max => RI_basis_parameter(kkind)%lmax 572 lc_min => RI_basis_parameter(kkind)%lmin 573 npgfc => RI_basis_parameter(kkind)%npgf 574 nsetc = RI_basis_parameter(kkind)%nset 575 zetc => RI_basis_parameter(kkind)%zet 576 nsgfc => RI_basis_parameter(kkind)%nsgf 577 sphi_c_ext => RI_basis_parameter(kkind)%sphi_ext(:, :, :, :) 578 nsgfl_c => RI_basis_parameter(kkind)%nsgfl 579 sphi_c_u1 = UBOUND(sphi_c_ext, 1) 580 sphi_c_u2 = UBOUND(sphi_c_ext, 2) 581 sphi_c_u3 = UBOUND(sphi_c_ext, 3) 582 583 ld_max => basis_S0(lkind)%lmax 584 ld_min => basis_S0(lkind)%lmin 585 npgfd => basis_S0(lkind)%npgf 586 nsetd = RI_basis_parameter(lkind)%nset 587 zetd => basis_S0(lkind)%zet 588 nsgfd => basis_S0(lkind)%nsgf 589 sphi_d_ext => basis_S0(lkind)%sphi_ext(:, :, :, :) 590 nsgfl_d => basis_S0(lkind)%nsgfl 591 sphi_d_u1 = UBOUND(sphi_d_ext, 1) 592 sphi_d_u2 = UBOUND(sphi_d_ext, 2) 593 sphi_d_u3 = UBOUND(sphi_d_ext, 3) 594 595 jset = 1 596 lset = 1 597 DO iset = 1, nseta 598 599 ncob = npgfb(jset)*ncoset(lb_max(jset)) 600 sphi_a_ext_set => sphi_a_ext(:, :, :, iset) 601 sphi_b_ext_set => sphi_b_ext(:, :, :, jset) 602 603 L_B_i_start = RI_index_table(iatom, iset) 604 L_B_i_end = RI_index_table(iatom, iset) + nsgfa(iset) - 1 605 606 kset_start = 1 607 IF (iatom == katom) kset_start = iset 608 DO kset = kset_start, nsetc 609 610 counter_L_blocks = counter_L_blocks + 1 611 IF (MOD(counter_L_blocks, para_env%num_pe) /= para_env%mepos) CYCLE 612 613 sphi_c_ext_set => sphi_c_ext(:, :, :, kset) 614 sphi_d_ext_set => sphi_d_ext(:, :, :, lset) 615 616 L_B_k_start = RI_index_table(katom, kset) 617 L_B_k_end = RI_index_table(katom, kset) + nsgfc(kset) - 1 618 619 pmax_entry = 0.0_dp 620 log10_pmax = pmax_entry 621 log10_eps_schwarz = log_zero 622 623 max_contraction_val = 1.0000_dp 624 625 ALLOCATE (L_block(nsgfa(iset), nsgfc(kset))) 626 627 CALL coulomb4(private_lib, ra, rb, rc, rd, npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), & 628 la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), & 629 lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), & 630 nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & 631 sphi_a_u1, sphi_a_u2, sphi_a_u3, & 632 sphi_b_u1, sphi_b_u2, sphi_b_u3, & 633 sphi_c_u1, sphi_c_u2, sphi_c_u3, & 634 sphi_d_u1, sphi_d_u2, sphi_d_u3, & 635 zeta(1:npgfa(iset), iset), zetb(1:npgfb(jset), jset), & 636 zetc(1:npgfc(kset), kset), zetd(1:npgfd(lset), lset), & 637 primitive_integrals, & 638 mp2_potential_parameter, & 639 actual_x_data%neighbor_cells, screen_coeffs_set(1, 1, 1, 1)%x, & 640 screen_coeffs_set(1, 1, 1, 1)%x, eps_schwarz, & 641 max_contraction_val, cartesian_estimate, cell, neris_tmp, & 642 log10_pmax, log10_eps_schwarz, & 643 tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2, & 644 pgf_list_ij, pgf_list_kl, pgf_product_list, & 645 nsgfl_a(:, iset), nsgfl_b(:, jset), & 646 nsgfl_c(:, kset), nsgfl_d(:, lset), & 647 sphi_a_ext_set, & 648 sphi_b_ext_set, & 649 sphi_c_ext_set, & 650 sphi_d_ext_set, & 651 ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, & 652 nimages, do_periodic, p_work) 653 654 primitive_counter = 0 655 DO llB = 1, nsgfd(lset) 656 DO kkB = 1, nsgfc(kset) 657 DO jjB = 1, nsgfb(jset) 658 DO iiB = 1, nsgfa(iset) 659 primitive_counter = primitive_counter + 1 660 L_block(iiB, kkB) = primitive_integrals(primitive_counter) 661 END DO 662 END DO 663 END DO 664 END DO 665 666 L_full_matrix(L_B_i_start:L_B_i_end, L_B_k_start:L_B_k_end) = L_block 667 L_full_matrix(L_B_k_start:L_B_k_end, L_B_i_start:L_B_i_end) = TRANSPOSE(L_block) 668 669 DEALLOCATE (L_block) 670 671 END DO ! kset 672 END DO ! iset 673 674 END DO !katom 675 END DO !iatom 676 677 ! now create a fm_matrix for the L matrix 678 CALL mp_sum(L_full_matrix, para_env%group) 679 680 ! create a sub blacs_env 681 NULLIFY (blacs_env) 682 CALL cp_blacs_env_create(blacs_env=blacs_env, para_env=para_env) 683 684 NULLIFY (fm_matrix_L) 685 NULLIFY (fm_struct_L) 686 CALL cp_fm_struct_create(fm_struct_L, context=blacs_env, nrow_global=RI_dimen, & 687 ncol_global=RI_dimen, para_env=para_env) 688 CALL cp_fm_create(fm_matrix_L, fm_struct_L, name="fm_matrix_L") 689 CALL cp_fm_struct_release(fm_struct_L) 690 CALL cp_blacs_env_release(blacs_env) 691 692 CALL cp_fm_set_submatrix(fm=fm_matrix_L, new_values=L_full_matrix, start_row=1, start_col=1, & 693 n_rows=RI_dimen, n_cols=RI_dimen) 694 695 info_chol = 0 696 CALL cp_fm_cholesky_decompose(matrix=fm_matrix_L, n=RI_dimen, info_out=info_chol) 697 CPASSERT(info_chol == 0) 698 699 ! triangual invert 700 CALL cp_fm_triangular_invert(matrix_a=fm_matrix_L, uplo_tr='U') 701 702 ! replicate L matrix to each proc 703 L_full_matrix = 0.0_dp 704 CALL cp_fm_get_submatrix(fm_matrix_L, L_full_matrix, 1, 1, RI_dimen, RI_dimen, .FALSE.) 705 CALL cp_fm_release(fm_matrix_L) 706 707 ! clean lower part 708 DO iiB = 1, RI_dimen 709 L_full_matrix(iiB + 1:RI_dimen, iiB) = 0.0_dp 710 END DO 711 712 ALLOCATE (list_kl%elements(natom**2)) 713 714 coeffs_kind_max0 = MAXVAL(screen_coeffs_kind(:, :)%x(2)) 715 ALLOCATE (set_list_kl((max_set*natom)**2)) 716 717 !! precalculate maximum density matrix elements in blocks 718 actual_x_data%pmax_block = 0.0_dp 719 shm_pmax_block => actual_x_data%pmax_block 720 721 shm_atomic_pair_list => actual_x_data%atomic_pair_list 722 723 iatom_start = 1 724 iatom_end = natom 725 jatom_start = 1 726 jatom_end = natom 727 katom_start = 1 728 katom_end = natom 729 latom_start = 1 730 latom_end = natom 731 732 CALL build_pair_list_mp2(natom, list_kl, set_list_kl, katom_start, katom_end, & 733 latom_start, latom_end, & 734 kind_of, basis_parameter, particle_set, & 735 do_periodic, screen_coeffs_set, screen_coeffs_kind, & 736 coeffs_kind_max0, log10_eps_schwarz, cell, 0.D+00, & 737 shm_atomic_pair_list) 738 739 virtual = dimen - occupied 740 741 ALLOCATE (Lai(RI_dimen, virtual, occupied)) 742 Lai = 0.0_dp 743 744 DO iatom = 1, natom 745 jatom = iatom 746 747 ikind = kind_of(iatom) 748 jkind = kind_of(jatom) 749 ra = particle_set(iatom)%r(:) 750 rb = particle_set(jatom)%r(:) 751 752 la_max => RI_basis_parameter(ikind)%lmax 753 la_min => RI_basis_parameter(ikind)%lmin 754 npgfa => RI_basis_parameter(ikind)%npgf 755 nseta = RI_basis_parameter(ikind)%nset 756 zeta => RI_basis_parameter(ikind)%zet 757 nsgfa => RI_basis_parameter(ikind)%nsgf 758 sphi_a_ext => RI_basis_parameter(ikind)%sphi_ext(:, :, :, :) 759 nsgfl_a => RI_basis_parameter(ikind)%nsgfl 760 sphi_a_u1 = UBOUND(sphi_a_ext, 1) 761 sphi_a_u2 = UBOUND(sphi_a_ext, 2) 762 sphi_a_u3 = UBOUND(sphi_a_ext, 3) 763 764 lb_max => basis_S0(jkind)%lmax 765 lb_min => basis_S0(jkind)%lmin 766 npgfb => basis_S0(jkind)%npgf 767 nsetb = basis_S0(jkind)%nset 768 zetb => basis_S0(jkind)%zet 769 nsgfb => basis_S0(jkind)%nsgf 770 sphi_b_ext => basis_S0(jkind)%sphi_ext(:, :, :, :) 771 nsgfl_b => basis_S0(jkind)%nsgfl 772 sphi_b_u1 = UBOUND(sphi_b_ext, 1) 773 sphi_b_u2 = UBOUND(sphi_b_ext, 2) 774 sphi_b_u3 = UBOUND(sphi_b_ext, 3) 775 776 jset = 1 777 DO iset = 1, nseta 778 779 counter_L_blocks = counter_L_blocks + 1 780 IF (MOD(counter_L_blocks, para_env%num_pe) /= para_env%mepos) CYCLE 781 782 ncob = npgfb(jset)*ncoset(lb_max(jset)) 783 sphi_a_ext_set => sphi_a_ext(:, :, :, iset) 784 sphi_b_ext_set => sphi_b_ext(:, :, :, jset) 785 786 L_B_i_start = RI_index_table(iatom, iset) 787 L_B_i_end = RI_index_table(iatom, iset) + nsgfa(iset) - 1 788 789 ALLOCATE (BI1(dimen, dimen, nsgfa(iset))) 790 BI1 = 0.0_dp 791 792 DO i_list_kl = 1, list_kl%n_element 793 794 katom = list_kl%elements(i_list_kl)%pair(1) 795 latom = list_kl%elements(i_list_kl)%pair(2) 796 797 i_set_list_kl_start = list_kl%elements(i_list_kl)%set_bounds(1) 798 i_set_list_kl_stop = list_kl%elements(i_list_kl)%set_bounds(2) 799 kkind = list_kl%elements(i_list_kl)%kind_pair(1) 800 lkind = list_kl%elements(i_list_kl)%kind_pair(2) 801 rc = list_kl%elements(i_list_kl)%r1 802 rd = list_kl%elements(i_list_kl)%r2 803 804 pmax_atom = 0.0_dp 805 806 lc_max => basis_parameter(kkind)%lmax 807 lc_min => basis_parameter(kkind)%lmin 808 npgfc => basis_parameter(kkind)%npgf 809 zetc => basis_parameter(kkind)%zet 810 nsgfc => basis_parameter(kkind)%nsgf 811 sphi_c_ext => basis_parameter(kkind)%sphi_ext(:, :, :, :) 812 nsgfl_c => basis_parameter(kkind)%nsgfl 813 sphi_c_u1 = UBOUND(sphi_c_ext, 1) 814 sphi_c_u2 = UBOUND(sphi_c_ext, 2) 815 sphi_c_u3 = UBOUND(sphi_c_ext, 3) 816 817 ld_max => basis_parameter(lkind)%lmax 818 ld_min => basis_parameter(lkind)%lmin 819 npgfd => basis_parameter(lkind)%npgf 820 zetd => basis_parameter(lkind)%zet 821 nsgfd => basis_parameter(lkind)%nsgf 822 sphi_d_ext => basis_parameter(lkind)%sphi_ext(:, :, :, :) 823 nsgfl_d => basis_parameter(lkind)%nsgfl 824 sphi_d_u1 = UBOUND(sphi_d_ext, 1) 825 sphi_d_u2 = UBOUND(sphi_d_ext, 2) 826 sphi_d_u3 = UBOUND(sphi_d_ext, 3) 827 828 DO i_set_list_kl = i_set_list_kl_start, i_set_list_kl_stop 829 kset = set_list_kl(i_set_list_kl)%pair(1) 830 lset = set_list_kl(i_set_list_kl)%pair(2) 831 832 IF (katom == latom .AND. lset < kset) CYCLE 833 834 orb_k_start = mp2_biel%index_table(katom, kset) 835 orb_k_end = orb_k_start + nsgfc(kset) - 1 836 orb_l_start = mp2_biel%index_table(latom, lset) 837 orb_l_end = orb_l_start + nsgfd(lset) - 1 838 839 !! get max_vals if we screen on initial density 840 pmax_entry = 0.0_dp 841 842 sphi_c_ext_set => sphi_c_ext(:, :, :, kset) 843 sphi_d_ext_set => sphi_d_ext(:, :, :, lset) 844 845 log10_pmax = pmax_entry 846 log10_eps_schwarz = log_zero 847 848 IF (ALLOCATED(MNRS)) DEALLOCATE (MNRS) 849 ALLOCATE (MNRS(nsgfd(lset), nsgfc(kset), nsgfa(iset))) 850 851 MNRS = 0.0_dp 852 853 max_contraction_val = max_contraction(kset, katom)*max_contraction(lset, latom) 854 855 CALL coulomb4(private_lib, ra, rb, rc, rd, npgfa(iset), npgfb(jset), npgfc(kset), npgfd(lset), & 856 la_min(iset), la_max(iset), lb_min(jset), lb_max(jset), & 857 lc_min(kset), lc_max(kset), ld_min(lset), ld_max(lset), & 858 nsgfa(iset), nsgfb(jset), nsgfc(kset), nsgfd(lset), & 859 sphi_a_u1, sphi_a_u2, sphi_a_u3, & 860 sphi_b_u1, sphi_b_u2, sphi_b_u3, & 861 sphi_c_u1, sphi_c_u2, sphi_c_u3, & 862 sphi_d_u1, sphi_d_u2, sphi_d_u3, & 863 zeta(1:npgfa(iset), iset), zetb(1:npgfb(jset), jset), & 864 zetc(1:npgfc(kset), kset), zetd(1:npgfd(lset), lset), & 865 primitive_integrals, & 866 mp2_potential_parameter, & 867 actual_x_data%neighbor_cells, screen_coeffs_set(kset, kset, kkind, kkind)%x, & 868 screen_coeffs_set(kset, kset, kkind, kkind)%x, eps_schwarz, & 869 max_contraction_val, cartesian_estimate, cell, neris_tmp, & 870 log10_pmax, log10_eps_schwarz, & 871 tmp_R_1, tmp_R_2, tmp_screen_pgf1, tmp_screen_pgf2, & 872 pgf_list_ij, pgf_list_kl, pgf_product_list, & 873 nsgfl_a(:, iset), nsgfl_b(:, jset), & 874 nsgfl_c(:, kset), nsgfl_d(:, lset), & 875 sphi_a_ext_set, & 876 sphi_b_ext_set, & 877 sphi_c_ext_set, & 878 sphi_d_ext_set, & 879 ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp, & 880 nimages, do_periodic, p_work) 881 882 nints = nsgfa(iset)*nsgfb(jset)*nsgfc(kset)*nsgfd(lset) 883 neris_total = neris_total + nints 884 nprim_ints = nprim_ints + neris_tmp 885 IF (cartesian_estimate == 0.0_dp) cartesian_estimate = TINY(cartesian_estimate) 886 estimate_to_store_int = EXPONENT(cartesian_estimate) 887 estimate_to_store_int = MAX(estimate_to_store_int, -15_int_8) 888 cartesian_estimate = SET_EXPONENT(1.0_dp, estimate_to_store_int + 1) 889 890 primitive_counter = 0 891 DO llB = 1, nsgfd(lset) 892 DO kkB = 1, nsgfc(kset) 893 DO jjB = 1, nsgfb(jset) 894 DO iiB = 1, nsgfa(iset) 895 primitive_counter = primitive_counter + 1 896 MNRS(llB, kkB, iiB) = primitive_integrals(primitive_counter) 897 END DO 898 END DO 899 END DO 900 END DO 901 902 DO iiB = 1, nsgfa(iset) 903 BI1(orb_l_start:orb_l_end, orb_k_start:orb_k_end, iiB) = MNRS(:, :, iiB) 904 BI1(orb_k_start:orb_k_end, orb_l_start:orb_l_end, iiB) = TRANSPOSE(MNRS(:, :, iiB)) 905 END DO 906 907 END DO ! i_set_list_kl 908 END DO ! i_list_kl 909 910 DO iiB = 1, nsgfa(iset) 911 BI1(1:virtual, 1:occupied, iiB) = MATMUL(TRANSPOSE(C(1:dimen, occupied + 1:dimen)), & 912 MATMUL(BI1(1:dimen, 1:dimen, iiB), C(1:dimen, 1:occupied))) 913 Lai(L_B_i_start + iiB - 1, 1:virtual, 1:occupied) = BI1(1:virtual, 1:occupied, iiB) 914 END DO 915 916 DEALLOCATE (BI1) 917 918 END DO 919 END DO 920 921 CALL mp_sum(Lai, para_env%group) 922 923 DO iiB = 1, occupied 924 IF (MOD(iiB, para_env%num_pe) == para_env%mepos) THEN 925 Lai(1:RI_dimen, 1:virtual, iiB) = MATMUL(TRANSPOSE(L_full_matrix), Lai(1:RI_dimen, 1:virtual, iiB)) 926 ELSE 927 Lai(:, :, iiB) = 0.0_dp 928 END IF 929 END DO 930 931 CALL mp_sum(Lai, para_env%group) 932 933 DEALLOCATE (set_list_kl) 934 935 DO i = 1, max_pgf**2 936 DEALLOCATE (pgf_list_ij(i)%image_list) 937 DEALLOCATE (pgf_list_kl(i)%image_list) 938 END DO 939 940 DEALLOCATE (pgf_list_ij) 941 DEALLOCATE (pgf_list_kl) 942 DEALLOCATE (pgf_product_list) 943 944 DEALLOCATE (max_contraction, kind_of) 945 946 DEALLOCATE (ee_work, ee_work2, ee_buffer1, ee_buffer2, ee_primitives_tmp) 947 948 DEALLOCATE (nimages) 949 950 IF (mp2_env%potential_parameter%potential_type == do_potential_TShPSC) THEN 951 init_TShPSC_lmax = -1 952 CALL free_C0() 953 END IF 954 955 CALL timestop(handle) 956 957 END SUBROUTINE calc_lai_libint 958 959! ************************************************************************************************** 960!> \brief ... 961!> \param cell ... 962!> \param qs_env ... 963!> \param mp2_env ... 964!> \param para_env ... 965!> \param mp2_potential_parameter ... 966!> \param actual_x_data ... 967!> \param do_periodic ... 968!> \param basis_parameter ... 969!> \param max_set ... 970!> \param particle_set ... 971!> \param natom ... 972!> \param kind_of ... 973!> \param nsgf_max ... 974!> \param primitive_integrals ... 975!> \param ee_work ... 976!> \param ee_work2 ... 977!> \param ee_buffer1 ... 978!> \param ee_buffer2 ... 979!> \param ee_primitives_tmp ... 980!> \param nspins ... 981!> \param max_contraction ... 982!> \param max_pgf ... 983!> \param pgf_list_ij ... 984!> \param pgf_list_kl ... 985!> \param pgf_product_list ... 986!> \param nimages ... 987!> \param eps_schwarz ... 988!> \param log10_eps_schwarz ... 989!> \param private_lib ... 990!> \param p_work ... 991!> \param screen_coeffs_set ... 992!> \param screen_coeffs_kind ... 993!> \param screen_coeffs_pgf ... 994!> \param radii_pgf ... 995!> \param RI_basis_parameter ... 996!> \param RI_basis_info ... 997! ************************************************************************************************** 998 SUBROUTINE prepare_integral_calc(cell, qs_env, mp2_env, para_env, mp2_potential_parameter, actual_x_data, & 999 do_periodic, basis_parameter, max_set, particle_set, natom, kind_of, & 1000 nsgf_max, primitive_integrals, ee_work, ee_work2, ee_buffer1, ee_buffer2, & 1001 ee_primitives_tmp, nspins, max_contraction, max_pgf, pgf_list_ij, & 1002 pgf_list_kl, pgf_product_list, nimages, eps_schwarz, log10_eps_schwarz, & 1003 private_lib, p_work, screen_coeffs_set, screen_coeffs_kind, screen_coeffs_pgf, & 1004 radii_pgf, RI_basis_parameter, RI_basis_info) 1005 1006 TYPE(cell_type), POINTER :: cell 1007 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env 1008 TYPE(mp2_type), INTENT(IN), POINTER :: mp2_env 1009 TYPE(cp_para_env_type), INTENT(IN), POINTER :: para_env 1010 TYPE(hfx_potential_type), INTENT(OUT) :: mp2_potential_parameter 1011 TYPE(hfx_type), POINTER :: actual_x_data 1012 LOGICAL, INTENT(OUT) :: do_periodic 1013 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter 1014 INTEGER, INTENT(OUT) :: max_set 1015 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1016 INTEGER, INTENT(OUT) :: natom 1017 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: kind_of 1018 INTEGER, INTENT(OUT) :: nsgf_max 1019 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), & 1020 INTENT(OUT) :: primitive_integrals, ee_work, ee_work2, & 1021 ee_buffer1, ee_buffer2, & 1022 ee_primitives_tmp 1023 INTEGER, INTENT(OUT) :: nspins 1024 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 1025 INTENT(OUT) :: max_contraction 1026 INTEGER, INTENT(OUT) :: max_pgf 1027 TYPE(hfx_pgf_list), ALLOCATABLE, DIMENSION(:), & 1028 INTENT(OUT) :: pgf_list_ij, pgf_list_kl 1029 TYPE(hfx_pgf_product_list), ALLOCATABLE, & 1030 DIMENSION(:), INTENT(OUT) :: pgf_product_list 1031 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: nimages 1032 REAL(KIND=dp), INTENT(OUT) :: eps_schwarz, log10_eps_schwarz 1033 TYPE(cp_libint_t), INTENT(OUT) :: private_lib 1034 REAL(KIND=dp), DIMENSION(:), POINTER :: p_work 1035 TYPE(hfx_screen_coeff_type), & 1036 DIMENSION(:, :, :, :), POINTER :: screen_coeffs_set 1037 TYPE(hfx_screen_coeff_type), DIMENSION(:, :), & 1038 POINTER :: screen_coeffs_kind 1039 TYPE(hfx_screen_coeff_type), & 1040 DIMENSION(:, :, :, :, :, :), POINTER :: screen_coeffs_pgf, radii_pgf 1041 TYPE(hfx_basis_type), DIMENSION(:), OPTIONAL, & 1042 POINTER :: RI_basis_parameter 1043 TYPE(hfx_basis_info_type), INTENT(IN), OPTIONAL :: RI_basis_info 1044 1045 CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_integral_calc', & 1046 routineP = moduleN//':'//routineN 1047 1048 INTEGER :: handle, i_thread, ii, ikind, irep, & 1049 l_max, n_threads, ncos_max, nkind, & 1050 nneighbors 1051 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1052 TYPE(dft_control_type), POINTER :: dft_control 1053 TYPE(hfx_basis_info_type), POINTER :: basis_info 1054 TYPE(hfx_type), POINTER :: shm_master_x_data 1055 1056 CALL timeset(routineN, handle) 1057 1058 NULLIFY (dft_control) 1059 1060 irep = 1 1061 1062 CALL get_qs_env(qs_env, & 1063 atomic_kind_set=atomic_kind_set, & 1064 cell=cell, & 1065 dft_control=dft_control) 1066 1067 !! Calculate l_max used in fgamma , because init_md_ftable is definitely not thread safe 1068 nkind = SIZE(atomic_kind_set, 1) 1069 l_max = 0 1070 DO ikind = 1, nkind 1071 l_max = MAX(l_max, MAXVAL(qs_env%x_data(1, 1)%basis_parameter(ikind)%lmax)) 1072 ENDDO 1073 IF (PRESENT(RI_basis_parameter)) THEN 1074 DO ikind = 1, nkind 1075 l_max = MAX(l_max, MAXVAL(RI_basis_parameter(ikind)%lmax)) 1076 ENDDO 1077 END IF 1078 l_max = 4*l_max 1079 CALL init_md_ftable(l_max) 1080 1081 CALL get_potential_parameters(mp2_env, l_max, para_env, mp2_potential_parameter) 1082 1083 n_threads = 1 1084 1085 i_thread = 0 1086 1087 actual_x_data => qs_env%x_data(irep, i_thread + 1) 1088 1089 shm_master_x_data => qs_env%x_data(irep, 1) 1090 1091 do_periodic = actual_x_data%periodic_parameter%do_periodic 1092 1093 IF (do_periodic) THEN 1094 ! ** Rebuild neighbor lists in case the cell has changed (i.e. NPT MD) 1095 actual_x_data%periodic_parameter%number_of_shells = actual_x_data%periodic_parameter%mode 1096 CALL hfx_create_neighbor_cells(actual_x_data, actual_x_data%periodic_parameter%number_of_shells_from_input, & 1097 cell, i_thread) 1098 END IF 1099 1100 basis_parameter => actual_x_data%basis_parameter 1101 basis_info => actual_x_data%basis_info 1102 1103 max_set = basis_info%max_set 1104 IF (PRESENT(RI_basis_info)) max_set = MAX(max_set, RI_basis_info%max_set) 1105 1106 CALL get_qs_env(qs_env=qs_env, & 1107 atomic_kind_set=atomic_kind_set, & 1108 particle_set=particle_set) 1109 1110 natom = SIZE(particle_set, 1) 1111 1112 ALLOCATE (kind_of(natom)) 1113 1114 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, & 1115 kind_of=kind_of) 1116 1117 !! precompute maximum nco and allocate scratch 1118 ncos_max = 0 1119 nsgf_max = 0 1120 CALL get_ncos_and_ncsgf(natom, kind_of, basis_parameter, ncos_max, nsgf_max) 1121 IF (PRESENT(RI_basis_parameter)) THEN 1122 CALL get_ncos_and_ncsgf(natom, kind_of, RI_basis_parameter, ncos_max, nsgf_max) 1123 END IF 1124 1125 !! Allocate the arrays for the integrals. 1126 ALLOCATE (primitive_integrals(nsgf_max**4)) 1127 primitive_integrals = 0.0_dp 1128 1129 ALLOCATE (ee_work(ncos_max**4)) 1130 ALLOCATE (ee_work2(ncos_max**4)) 1131 ALLOCATE (ee_buffer1(ncos_max**4)) 1132 ALLOCATE (ee_buffer2(ncos_max**4)) 1133 ALLOCATE (ee_primitives_tmp(nsgf_max**4)) 1134 1135 nspins = dft_control%nspins 1136 1137 CALL get_max_contraction(max_contraction, max_set, natom, max_pgf, kind_of, basis_parameter) 1138 1139 ! ** Allocate buffers for pgf_lists 1140 nneighbors = SIZE(actual_x_data%neighbor_cells) 1141 ALLOCATE (pgf_list_ij(max_pgf**2)) 1142 ALLOCATE (pgf_list_kl(max_pgf**2)) 1143 ALLOCATE (pgf_product_list(nneighbors**3)) 1144 ALLOCATE (nimages(max_pgf**2)) 1145 1146 DO ii = 1, max_pgf**2 1147 ALLOCATE (pgf_list_ij(ii)%image_list(nneighbors)) 1148 ALLOCATE (pgf_list_kl(ii)%image_list(nneighbors)) 1149 END DO 1150 1151 !!! Skipped part 1152 1153 !! Get screening parameter 1154 eps_schwarz = actual_x_data%screening_parameter%eps_schwarz 1155 IF (eps_schwarz <= 0.0_dp) THEN 1156 log10_eps_schwarz = log_zero 1157 ELSE 1158 log10_eps_schwarz = LOG10(eps_schwarz) 1159 END IF 1160 1161 !! The number of integrals that fit into the given MAX_MEMORY 1162 1163 private_lib = actual_x_data%lib 1164 1165 !!!!!!!! Missing part on the density matrix 1166 1167 !! Initialize schwarz screening matrices for near field estimates and boxing screening matrices 1168 !! for far field estimates. The update is only performed if the geomtry of the system changed. 1169 !! If the system is periodic, then the corresponding routines are called and some variables 1170 !! are initialized 1171 1172 IF (.NOT. shm_master_x_data%screen_funct_is_initialized) THEN 1173 CALL calc_pair_dist_radii(qs_env, basis_parameter, & 1174 shm_master_x_data%pair_dist_radii_pgf, max_set, max_pgf, eps_schwarz, & 1175 n_threads, i_thread) 1176 CALL calc_screening_functions(qs_env, basis_parameter, private_lib, shm_master_x_data%potential_parameter, & 1177 shm_master_x_data%screen_funct_coeffs_set, & 1178 shm_master_x_data%screen_funct_coeffs_kind, & 1179 shm_master_x_data%screen_funct_coeffs_pgf, & 1180 shm_master_x_data%pair_dist_radii_pgf, & 1181 max_set, max_pgf, n_threads, i_thread, p_work) 1182 1183 shm_master_x_data%screen_funct_is_initialized = .TRUE. 1184 END IF 1185 1186 screen_coeffs_set => shm_master_x_data%screen_funct_coeffs_set 1187 screen_coeffs_kind => shm_master_x_data%screen_funct_coeffs_kind 1188 screen_coeffs_pgf => shm_master_x_data%screen_funct_coeffs_pgf 1189 radii_pgf => shm_master_x_data%pair_dist_radii_pgf 1190 1191 CALL timestop(handle) 1192 1193 END SUBROUTINE prepare_integral_calc 1194 1195! ************************************************************************************************** 1196!> \brief ... 1197!> \param mp2_env ... 1198!> \param l_max ... 1199!> \param para_env ... 1200!> \param mp2_potential_parameter ... 1201! ************************************************************************************************** 1202 SUBROUTINE get_potential_parameters(mp2_env, l_max, para_env, mp2_potential_parameter) 1203 1204 TYPE(mp2_type), POINTER :: mp2_env 1205 INTEGER, INTENT(IN) :: l_max 1206 TYPE(cp_para_env_type), POINTER :: para_env 1207 TYPE(hfx_potential_type), INTENT(OUT) :: mp2_potential_parameter 1208 1209 INTEGER :: unit_id 1210 1211 IF (mp2_env%potential_parameter%potential_type == do_potential_TShPSC) THEN 1212 IF (l_max > init_TShPSC_lmax) THEN 1213 IF (para_env%ionode) THEN 1214 CALL open_file(unit_number=unit_id, file_name=mp2_env%potential_parameter%filename) 1215 END IF 1216 CALL init(l_max, unit_id, para_env%mepos, para_env%group) 1217 IF (para_env%ionode) THEN 1218 CALL close_file(unit_id) 1219 END IF 1220 init_TShPSC_lmax = l_max 1221 END IF 1222 mp2_potential_parameter%cutoff_radius = mp2_env%potential_parameter%cutoff_radius/2.0_dp 1223 ELSE IF (mp2_env%potential_parameter%potential_type == do_potential_long) THEN 1224 mp2_potential_parameter%omega = mp2_env%potential_parameter%omega 1225 END IF 1226 1227 mp2_potential_parameter%potential_type = mp2_env%potential_parameter%potential_type 1228 1229 END SUBROUTINE get_potential_parameters 1230 1231! ************************************************************************************************** 1232!> \brief ... 1233!> \param natom ... 1234!> \param kind_of ... 1235!> \param basis_parameter ... 1236!> \param ncos_max ... 1237!> \param nsgf_max ... 1238! ************************************************************************************************** 1239 SUBROUTINE get_ncos_and_ncsgf(natom, kind_of, basis_parameter, ncos_max, nsgf_max) 1240 1241 INTEGER, INTENT(IN) :: natom 1242 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of 1243 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter 1244 INTEGER, INTENT(INOUT) :: ncos_max, nsgf_max 1245 1246 INTEGER :: iatom, ikind, iset, nseta 1247 INTEGER, DIMENSION(:), POINTER :: la_max, npgfa, nsgfa 1248 1249 DO iatom = 1, natom 1250 ikind = kind_of(iatom) 1251 nseta = basis_parameter(ikind)%nset 1252 npgfa => basis_parameter(ikind)%npgf 1253 la_max => basis_parameter(ikind)%lmax 1254 nsgfa => basis_parameter(ikind)%nsgf 1255 DO iset = 1, nseta 1256 ncos_max = MAX(ncos_max, ncoset(la_max(iset))) 1257 nsgf_max = MAX(nsgf_max, nsgfa(iset)) 1258 ENDDO 1259 ENDDO 1260 1261 END SUBROUTINE get_ncos_and_ncsgf 1262 1263! ************************************************************************************************** 1264!> \brief ... 1265!> \param max_contraction ... 1266!> \param max_set ... 1267!> \param natom ... 1268!> \param max_pgf ... 1269!> \param kind_of ... 1270!> \param basis_parameter ... 1271! ************************************************************************************************** 1272 SUBROUTINE get_max_contraction(max_contraction, max_set, natom, max_pgf, kind_of, basis_parameter) 1273 1274 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :), & 1275 INTENT(OUT) :: max_contraction 1276 INTEGER, INTENT(IN) :: max_set, natom 1277 INTEGER, INTENT(OUT) :: max_pgf 1278 INTEGER, DIMENSION(:), INTENT(IN) :: kind_of 1279 TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter 1280 1281 INTEGER :: i, jatom, jkind, jset, ncob, nsetb, sgfb 1282 INTEGER, DIMENSION(:), POINTER :: lb_max, npgfb, nsgfb 1283 INTEGER, DIMENSION(:, :), POINTER :: first_sgfb 1284 REAL(KIND=dp), DIMENSION(:, :), POINTER :: sphi_b 1285 1286 ALLOCATE (max_contraction(max_set, natom)) 1287 1288 max_contraction = 0.0_dp 1289 max_pgf = 0 1290 DO jatom = 1, natom 1291 jkind = kind_of(jatom) 1292 lb_max => basis_parameter(jkind)%lmax 1293 nsetb = basis_parameter(jkind)%nset 1294 npgfb => basis_parameter(jkind)%npgf 1295 first_sgfb => basis_parameter(jkind)%first_sgf 1296 sphi_b => basis_parameter(jkind)%sphi 1297 nsgfb => basis_parameter(jkind)%nsgf 1298 DO jset = 1, nsetb 1299 ! takes the primitive to contracted transformation into account 1300 ncob = npgfb(jset)*ncoset(lb_max(jset)) 1301 sgfb = first_sgfb(1, jset) 1302 ! if the primitives are assumed to be all of max_val2, max_val2*p2s_b becomes 1303 ! the maximum value after multiplication with sphi_b 1304 max_contraction(jset, jatom) = MAXVAL((/(SUM(ABS(sphi_b(1:ncob, i))), i=sgfb, sgfb + nsgfb(jset) - 1)/)) 1305 max_pgf = MAX(max_pgf, npgfb(jset)) 1306 ENDDO 1307 ENDDO 1308 1309 END SUBROUTINE get_max_contraction 1310 1311! ************************************************************************************************** 1312!> \brief ... 1313!> \param matrix ... 1314!> \param size_matrix ... 1315!> \return ... 1316! ************************************************************************************************** 1317 FUNCTION calc_cond_num(matrix, size_matrix) RESULT(cond_num) 1318 TYPE(cp_fm_type), INTENT(IN), POINTER :: matrix 1319 INTEGER, INTENT(IN) :: size_matrix 1320 REAL(KIND=dp) :: cond_num 1321 1322 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals 1323 TYPE(cp_fm_type), POINTER :: matrix_diag, matrix_work 1324 1325 ! create a work matrix that later contains the eigenvectors (not needed here) 1326 NULLIFY (matrix_diag) 1327 NULLIFY (matrix_work) 1328 1329 CALL cp_fm_create(matrix_diag, matrix%matrix_struct, name="fm_matrix_diag") 1330 CALL cp_fm_create(matrix_work, matrix%matrix_struct, name="fm_matrix_work") 1331 1332 CALL cp_fm_set_all(matrix=matrix_diag, alpha=0.0_dp) 1333 CALL cp_fm_set_all(matrix=matrix_work, alpha=0.0_dp) 1334 1335 CALL cp_fm_to_fm(source=matrix, destination=matrix_diag) 1336 1337 ALLOCATE (evals(size_matrix)) 1338 1339 evals = 0.0_dp 1340 CALL choose_eigv_solver(matrix=matrix_diag, eigenvectors=matrix_work, eigenvalues=evals) 1341 1342 cond_num = MAXVAL(ABS(evals))/MINVAL(ABS(evals)) 1343 1344 CALL cp_fm_release(matrix_diag) 1345 CALL cp_fm_release(matrix_work) 1346 1347 DEALLOCATE (evals) 1348 END FUNCTION calc_cond_num 1349 1350END MODULE mp2_ri_libint 1351