1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief initializes the environment for lri 8!> lri : local resolution of the identity 9!> \par History 10!> created [06.2015] 11!> \author Dorothea Golze 12! ************************************************************************************************** 13MODULE lri_environment_init 14 USE ao_util, ONLY: exp_radius 15 USE atomic_kind_types, ONLY: atomic_kind_type,& 16 get_atomic_kind_set 17 USE basis_set_types, ONLY: copy_gto_basis_set,& 18 gto_basis_set_type 19 USE bibliography, ONLY: Golze2017a,& 20 Golze2017b,& 21 cite_reference 22 USE cp_control_types, ONLY: dft_control_type 23 USE generic_shg_integrals, ONLY: int_overlap_aba_shg 24 USE generic_shg_integrals_init, ONLY: contraction_matrix_shg,& 25 contraction_matrix_shg_mix,& 26 get_clebsch_gordon_coefficients 27 USE input_section_types, ONLY: section_vals_type,& 28 section_vals_val_get 29 USE kinds, ONLY: dp 30 USE lri_environment_types, ONLY: deallocate_bas_properties,& 31 lri_env_create,& 32 lri_environment_type 33 USE mathconstants, ONLY: fac,& 34 pi 35 USE mathlib, ONLY: invert_matrix 36 USE qs_environment_types, ONLY: get_qs_env,& 37 qs_environment_type 38 USE qs_kind_types, ONLY: get_qs_kind,& 39 qs_kind_type 40#include "./base/base_uses.f90" 41 42 IMPLICIT NONE 43 44 PRIVATE 45 46! ************************************************************************************************** 47 48 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_environment_init' 49 50 PUBLIC :: lri_env_init, lri_env_basis, lri_basis_init 51 52! ************************************************************************************************** 53 54CONTAINS 55 56! ************************************************************************************************** 57!> \brief initializes the lri env 58!> \param lri_env ... 59!> \param lri_section ... 60! ************************************************************************************************** 61 SUBROUTINE lri_env_init(lri_env, lri_section) 62 63 TYPE(lri_environment_type), POINTER :: lri_env 64 TYPE(section_vals_type), POINTER :: lri_section 65 66 CHARACTER(len=*), PARAMETER :: routineN = 'lri_env_init', routineP = moduleN//':'//routineN 67 68 REAL(KIND=dp), DIMENSION(:), POINTER :: radii 69 70 NULLIFY (lri_env) 71 CALL lri_env_create(lri_env) 72 73 ! init keywords 74 ! use RI for local pp terms 75 CALL section_vals_val_get(lri_section, "RI_STATISTIC", & 76 l_val=lri_env%statistics) 77 ! use exact one-center terms 78 CALL section_vals_val_get(lri_section, "EXACT_1C_TERMS", & 79 l_val=lri_env%exact_1c_terms) 80 ! use RI for local pp terms 81 CALL section_vals_val_get(lri_section, "PPL_RI", & 82 l_val=lri_env%ppl_ri) 83 ! check for debug (OS scheme) 84 CALL section_vals_val_get(lri_section, "DEBUG_LRI_INTEGRALS", & 85 l_val=lri_env%debug) 86 ! integrals based on solid harmonic Gaussians 87 CALL section_vals_val_get(lri_section, "SHG_LRI_INTEGRALS", & 88 l_val=lri_env%use_shg_integrals) 89 ! how to calculate inverse/pseuodinverse of overlap 90 CALL section_vals_val_get(lri_section, "LRI_OVERLAP_MATRIX", & 91 i_val=lri_env%lri_overlap_inv) 92 CALL section_vals_val_get(lri_section, "MAX_CONDITION_NUM", & 93 r_val=lri_env%cond_max) 94 ! integrals threshold (aba, abb) 95 CALL section_vals_val_get(lri_section, "EPS_O3_INT", & 96 r_val=lri_env%eps_o3_int) 97 ! RI SINV 98 CALL section_vals_val_get(lri_section, "RI_SINV", & 99 c_val=lri_env%ri_sinv_app) 100 ! Distant Pair Approximation 101 CALL section_vals_val_get(lri_section, "DISTANT_PAIR_APPROXIMATION", & 102 l_val=lri_env%distant_pair_approximation) 103 CALL section_vals_val_get(lri_section, "DISTANT_PAIR_METHOD", & 104 c_val=lri_env%distant_pair_method) 105 CALL section_vals_val_get(lri_section, "DISTANT_PAIR_RADII", r_vals=radii) 106 CPASSERT(SIZE(radii) == 2) 107 CPASSERT(radii(2) > radii(1)) 108 CPASSERT(radii(1) > 0.0_dp) 109 lri_env%r_in = radii(1) 110 lri_env%r_out = radii(2) 111 112 CALL cite_reference(Golze2017b) 113 IF (lri_env%use_shg_integrals) CALL cite_reference(Golze2017a) 114 115 END SUBROUTINE lri_env_init 116! ************************************************************************************************** 117!> \brief initializes the lri env 118!> \param ri_type ... 119!> \param qs_env ... 120!> \param lri_env ... 121!> \param qs_kind_set ... 122! ************************************************************************************************** 123 SUBROUTINE lri_env_basis(ri_type, qs_env, lri_env, qs_kind_set) 124 125 CHARACTER(len=*), INTENT(IN) :: ri_type 126 TYPE(qs_environment_type), POINTER :: qs_env 127 TYPE(lri_environment_type), POINTER :: lri_env 128 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 129 130 CHARACTER(len=*), PARAMETER :: routineN = 'lri_env_basis', routineP = moduleN//':'//routineN 131 132 INTEGER :: i, i1, i2, iat, ikind, ip, ipgf, iset, ishell, jp, l, lmax_ikind_orb, & 133 lmax_ikind_ri, maxl_orb, maxl_ri, n1, n2, natom, nbas, nkind, nribas, nspin 134 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of 135 REAL(KIND=dp) :: gcc, rad, rai, raj, xradius, zeta 136 REAL(KIND=dp), DIMENSION(:), POINTER :: int_aux, norm 137 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 138 TYPE(dft_control_type), POINTER :: dft_control 139 TYPE(gto_basis_set_type), POINTER :: orb_basis_set, ri_basis_set 140 141 ! initialize the basic basis sets (orb and ri) 142 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set) 143 nkind = SIZE(atomic_kind_set) 144 ALLOCATE (lri_env%orb_basis(nkind), lri_env%ri_basis(nkind)) 145 maxl_orb = 0 146 maxl_ri = 0 147 DO ikind = 1, nkind 148 NULLIFY (orb_basis_set, ri_basis_set) 149 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB") 150 IF (ri_type == "LRI") THEN 151 CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="LRI_AUX") 152 ELSE IF (ri_type == "RI") THEN 153 CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="RI_HXC") 154 ELSE 155 CPABORT('ri_type') 156 END IF 157 NULLIFY (lri_env%orb_basis(ikind)%gto_basis_set) 158 NULLIFY (lri_env%ri_basis(ikind)%gto_basis_set) 159 IF (ASSOCIATED(orb_basis_set)) THEN 160 CALL copy_gto_basis_set(orb_basis_set, lri_env%orb_basis(ikind)%gto_basis_set) 161 CALL copy_gto_basis_set(ri_basis_set, lri_env%ri_basis(ikind)%gto_basis_set) 162 END IF 163 lmax_ikind_orb = MAXVAL(lri_env%orb_basis(ikind)%gto_basis_set%lmax) 164 lmax_ikind_ri = MAXVAL(lri_env%ri_basis(ikind)%gto_basis_set%lmax) 165 maxl_orb = MAX(maxl_orb, lmax_ikind_orb) 166 maxl_ri = MAX(maxl_ri, lmax_ikind_ri) 167 END DO 168 169 IF (ri_type == "LRI") THEN 170 ! CG coefficients needed for lri integrals 171 IF (ASSOCIATED(lri_env%cg_shg)) THEN 172 CALL get_clebsch_gordon_coefficients(lri_env%cg_shg%cg_coeff, & 173 lri_env%cg_shg%cg_none0_list, & 174 lri_env%cg_shg%ncg_none0, & 175 maxl_orb, maxl_ri) 176 ENDIF 177 CALL lri_basis_init(lri_env) 178 ! distant pair approximation 179 IF (lri_env%distant_pair_approximation) THEN 180 ! 181 SELECT CASE (lri_env%distant_pair_method) 182 CASE ("EW") 183 ! equal weight of 0.5 184 CASE ("AW") 185 ALLOCATE (lri_env%aradius(nkind)) 186 DO i = 1, nkind 187 orb_basis_set => lri_env%orb_basis(i)%gto_basis_set 188 lri_env%aradius(i) = orb_basis_set%kind_radius 189 END DO 190 CASE ("SW") 191 ALLOCATE (lri_env%wbas(nkind)) 192 DO i = 1, nkind 193 orb_basis_set => lri_env%orb_basis(i)%gto_basis_set 194 n1 = orb_basis_set%nsgf 195 ALLOCATE (lri_env%wbas(i)%vec(n1)) 196 DO iset = 1, orb_basis_set%nset 197 i1 = orb_basis_set%first_sgf(1, iset) 198 n2 = orb_basis_set%nshell(iset) 199 i2 = orb_basis_set%last_sgf(n2, iset) 200 lri_env%wbas(i)%vec(i1:i2) = orb_basis_set%set_radius(iset) 201 END DO 202 END DO 203 CASE ("LW") 204 ALLOCATE (lri_env%wbas(nkind)) 205 DO i = 1, nkind 206 orb_basis_set => lri_env%orb_basis(i)%gto_basis_set 207 n1 = orb_basis_set%nsgf 208 ALLOCATE (lri_env%wbas(i)%vec(n1)) 209 DO iset = 1, orb_basis_set%nset 210 DO ishell = 1, orb_basis_set%nshell(iset) 211 i1 = orb_basis_set%first_sgf(ishell, iset) 212 i2 = orb_basis_set%last_sgf(ishell, iset) 213 l = orb_basis_set%l(ishell, iset) 214 xradius = 0.0_dp 215 DO ipgf = 1, orb_basis_set%npgf(iset) 216 gcc = orb_basis_set%gcc(ipgf, ishell, iset) 217 zeta = orb_basis_set%zet(ipgf, iset) 218 rad = exp_radius(l, zeta, 1.e-5_dp, gcc) 219 xradius = MAX(xradius, rad) 220 END DO 221 lri_env%wbas(i)%vec(i1:i2) = xradius 222 END DO 223 END DO 224 END DO 225 CASE DEFAULT 226 CPABORT("Unknown DISTANT_PAIR_METHOD in LRI") 227 END SELECT 228 ! 229 ALLOCATE (lri_env%wmat(nkind, nkind)) 230 SELECT CASE (lri_env%distant_pair_method) 231 CASE ("EW") 232 ! equal weight of 0.5 233 DO i1 = 1, nkind 234 n1 = lri_env%orb_basis(i1)%gto_basis_set%nsgf 235 DO i2 = 1, nkind 236 n2 = lri_env%orb_basis(i2)%gto_basis_set%nsgf 237 ALLOCATE (lri_env%wmat(i1, i2)%mat(n1, n2)) 238 lri_env%wmat(i1, i2)%mat(:, :) = 0.5_dp 239 END DO 240 END DO 241 CASE ("AW") 242 DO i1 = 1, nkind 243 n1 = lri_env%orb_basis(i1)%gto_basis_set%nsgf 244 DO i2 = 1, nkind 245 n2 = lri_env%orb_basis(i2)%gto_basis_set%nsgf 246 ALLOCATE (lri_env%wmat(i1, i2)%mat(n1, n2)) 247 rai = lri_env%aradius(i1)**2 248 raj = lri_env%aradius(i2)**2 249 IF (raj > rai) THEN 250 lri_env%wmat(i1, i2)%mat(:, :) = 1.0_dp 251 ELSE 252 lri_env%wmat(i1, i2)%mat(:, :) = 0.0_dp 253 END IF 254 END DO 255 END DO 256 CASE ("SW", "LW") 257 DO i1 = 1, nkind 258 n1 = lri_env%orb_basis(i1)%gto_basis_set%nsgf 259 DO i2 = 1, nkind 260 n2 = lri_env%orb_basis(i2)%gto_basis_set%nsgf 261 ALLOCATE (lri_env%wmat(i1, i2)%mat(n1, n2)) 262 DO ip = 1, SIZE(lri_env%wbas(i1)%vec) 263 rai = lri_env%wbas(i1)%vec(ip)**2 264 DO jp = 1, SIZE(lri_env%wbas(i2)%vec) 265 raj = lri_env%wbas(i2)%vec(jp)**2 266 IF (raj > rai) THEN 267 lri_env%wmat(i1, i2)%mat(ip, jp) = 1.0_dp 268 ELSE 269 lri_env%wmat(i1, i2)%mat(ip, jp) = 0.0_dp 270 END IF 271 END DO 272 END DO 273 END DO 274 END DO 275 END SELECT 276 END IF 277 ELSE IF (ri_type == "RI") THEN 278 ALLOCATE (lri_env%ri_fit) 279 NULLIFY (lri_env%ri_fit%nvec) 280 NULLIFY (lri_env%ri_fit%bas_ptr) 281 CALL get_qs_env(qs_env=qs_env, natom=natom) 282 ! initialize pointers to RI basis vector 283 ALLOCATE (lri_env%ri_fit%bas_ptr(2, natom)) 284 ALLOCATE (kind_of(natom)) 285 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of) 286 nbas = 0 287 DO iat = 1, natom 288 ikind = kind_of(iat) 289 nribas = lri_env%ri_basis(ikind)%gto_basis_set%nsgf 290 lri_env%ri_fit%bas_ptr(1, iat) = nbas + 1 291 lri_env%ri_fit%bas_ptr(2, iat) = nbas + nribas 292 nbas = nbas + nribas 293 END DO 294 ! initialize vector t 295 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control) 296 nspin = dft_control%nspins 297 ALLOCATE (lri_env%ri_fit%tvec(nbas, nspin), lri_env%ri_fit%rm1t(nbas, nspin)) 298 ! initialize vector a, expansion of density 299 ALLOCATE (lri_env%ri_fit%avec(nbas, nspin)) 300 ! initialize vector fout, R^(-1)*(f-p*n) 301 ALLOCATE (lri_env%ri_fit%fout(nbas, nspin)) 302 ! initialize vector with RI basis integrated 303 NULLIFY (norm, int_aux) 304 nbas = lri_env%ri_fit%bas_ptr(2, natom) 305 ALLOCATE (lri_env%ri_fit%nvec(nbas), lri_env%ri_fit%rm1n(nbas)) 306 ikind = 0 307 DO iat = 1, natom 308 IF (ikind /= kind_of(iat)) THEN 309 ikind = kind_of(iat) 310 ri_basis_set => lri_env%ri_basis(ikind)%gto_basis_set 311 IF (ASSOCIATED(norm)) DEALLOCATE (norm) 312 IF (ASSOCIATED(int_aux)) DEALLOCATE (int_aux) 313 CALL basis_norm_s_func(ri_basis_set, norm) 314 CALL basis_int(ri_basis_set, int_aux, norm) 315 END IF 316 nbas = SIZE(int_aux) 317 i1 = lri_env%ri_fit%bas_ptr(1, iat) 318 i2 = lri_env%ri_fit%bas_ptr(2, iat) 319 lri_env%ri_fit%nvec(i1:i2) = int_aux(1:nbas) 320 END DO 321 IF (ASSOCIATED(norm)) DEALLOCATE (norm) 322 IF (ASSOCIATED(int_aux)) DEALLOCATE (int_aux) 323 DEALLOCATE (kind_of) 324 ELSE 325 CPABORT('ri_type') 326 END IF 327 328 END SUBROUTINE lri_env_basis 329 330! ************************************************************************************************** 331!> \brief initializes the lri basis: calculates the norm, self-overlap 332!> and integral of the ri basis 333!> \param lri_env ... 334! ************************************************************************************************** 335 SUBROUTINE lri_basis_init(lri_env) 336 TYPE(lri_environment_type), POINTER :: lri_env 337 338 CHARACTER(len=*), PARAMETER :: routineN = 'lri_basis_init', routineP = moduleN//':'//routineN 339 340 INTEGER :: ikind, nkind 341 INTEGER, DIMENSION(:, :, :), POINTER :: orb_index, ri_index 342 REAL(KIND=dp) :: delta 343 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: dovlp3 344 REAL(KIND=dp), DIMENSION(:), POINTER :: orb_norm_r, ri_int_fbas, ri_norm_r, & 345 ri_norm_s 346 REAL(KIND=dp), DIMENSION(:, :), POINTER :: orb_ovlp, ri_ovlp, ri_ovlp_inv 347 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_orb, scon_ri 348 REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scon_mix 349 TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis 350 351 IF (ASSOCIATED(lri_env)) THEN 352 IF (ASSOCIATED(lri_env%orb_basis)) THEN 353 CPASSERT(ASSOCIATED(lri_env%ri_basis)) 354 nkind = SIZE(lri_env%orb_basis) 355 CALL deallocate_bas_properties(lri_env) 356 ALLOCATE (lri_env%bas_prop(nkind)) 357 DO ikind = 1, nkind 358 NULLIFY (orb_basis, ri_basis) 359 orb_basis => lri_env%orb_basis(ikind)%gto_basis_set 360 IF (ASSOCIATED(orb_basis)) THEN 361 ri_basis => lri_env%ri_basis(ikind)%gto_basis_set 362 CPASSERT(ASSOCIATED(ri_basis)) 363 NULLIFY (ri_norm_r) 364 CALL basis_norm_radial(ri_basis, ri_norm_r) 365 NULLIFY (orb_norm_r) 366 CALL basis_norm_radial(orb_basis, orb_norm_r) 367 NULLIFY (ri_norm_s) 368 CALL basis_norm_s_func(ri_basis, ri_norm_s) 369 NULLIFY (ri_int_fbas) 370 CALL basis_int(ri_basis, ri_int_fbas, ri_norm_s) 371 lri_env%bas_prop(ikind)%int_fbas => ri_int_fbas 372 NULLIFY (ri_ovlp) 373 CALL basis_ovlp(ri_basis, ri_ovlp, ri_norm_r) 374 lri_env%bas_prop(ikind)%ri_ovlp => ri_ovlp 375 NULLIFY (orb_ovlp) 376 CALL basis_ovlp(orb_basis, orb_ovlp, orb_norm_r) 377 lri_env%bas_prop(ikind)%orb_ovlp => orb_ovlp 378 NULLIFY (scon_ri) 379 CALL contraction_matrix_shg(ri_basis, scon_ri) 380 lri_env%bas_prop(ikind)%scon_ri => scon_ri 381 NULLIFY (scon_orb) 382 CALL contraction_matrix_shg(orb_basis, scon_orb) 383 lri_env%bas_prop(ikind)%scon_orb => scon_orb 384 NULLIFY (scon_mix) 385 CALL contraction_matrix_shg_mix(orb_basis, ri_basis, & 386 orb_index, ri_index, scon_mix) 387 lri_env%bas_prop(ikind)%scon_mix => scon_mix 388 lri_env%bas_prop(ikind)%orb_index => orb_index 389 lri_env%bas_prop(ikind)%ri_index => ri_index 390 ALLOCATE (lri_env%bas_prop(ikind)%ovlp3(orb_basis%nsgf, orb_basis%nsgf, ri_basis%nsgf)) 391 ALLOCATE (dovlp3(orb_basis%nsgf, orb_basis%nsgf, ri_basis%nsgf, 3)) 392 CALL int_overlap_aba_shg(lri_env%bas_prop(ikind)%ovlp3, dovlp3, (/0.0_dp, 0.0_dp, 0.0_dp/), & 393 orb_basis, orb_basis, ri_basis, scon_orb, & 394 scon_mix, orb_index, ri_index, & 395 lri_env%cg_shg%cg_coeff, & 396 lri_env%cg_shg%cg_none0_list, & 397 lri_env%cg_shg%ncg_none0, & 398 calculate_forces=.FALSE.) 399 DEALLOCATE (orb_norm_r, ri_norm_r, ri_norm_s) 400 DEALLOCATE (dovlp3) 401 ALLOCATE (ri_ovlp_inv(ri_basis%nsgf, ri_basis%nsgf)) 402 CALL invert_matrix(ri_ovlp, ri_ovlp_inv, delta, improve=.TRUE.) 403 lri_env%bas_prop(ikind)%ri_ovlp_inv => ri_ovlp_inv 404 END IF 405 END DO 406 END IF 407 END IF 408 409 END SUBROUTINE lri_basis_init 410 411! ************************************************************************************************** 412!> \brief normalization for a contracted Gaussian s-function, 413!> spherical = cartesian Gaussian for s-functions 414!> \param basis ... 415!> \param norm ... 416! ************************************************************************************************** 417 SUBROUTINE basis_norm_s_func(basis, norm) 418 419 TYPE(gto_basis_set_type), POINTER :: basis 420 REAL(dp), DIMENSION(:), POINTER :: norm 421 422 CHARACTER(len=*), PARAMETER :: routineN = 'basis_norm_s_func', & 423 routineP = moduleN//':'//routineN 424 425 INTEGER :: ipgf, iset, isgf, ishell, jpgf, l, nbas 426 REAL(KIND=dp) :: aai, aaj, cci, ccj, expa, ppl 427 428 NULLIFY (norm) 429 430 nbas = basis%nsgf 431 ALLOCATE (norm(nbas)) 432 norm = 0._dp 433 434 DO iset = 1, basis%nset 435 DO ishell = 1, basis%nshell(iset) 436 l = basis%l(ishell, iset) 437 IF (l /= 0) CYCLE 438 expa = 0.5_dp*REAL(2*l + 3, dp) 439 ppl = pi**(3._dp/2._dp) 440 DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset) 441 DO ipgf = 1, basis%npgf(iset) 442 cci = basis%gcc(ipgf, ishell, iset) 443 aai = basis%zet(ipgf, iset) 444 DO jpgf = 1, basis%npgf(iset) 445 ccj = basis%gcc(jpgf, ishell, iset) 446 aaj = basis%zet(jpgf, iset) 447 norm(isgf) = norm(isgf) + cci*ccj*ppl/(aai + aaj)**expa 448 END DO 449 END DO 450 norm(isgf) = 1.0_dp/SQRT(norm(isgf)) 451 END DO 452 END DO 453 END DO 454 455 END SUBROUTINE basis_norm_s_func 456 457! ************************************************************************************************** 458!> \brief normalization for radial part of contracted spherical Gaussian 459!> functions 460!> \param basis ... 461!> \param norm ... 462! ************************************************************************************************** 463 SUBROUTINE basis_norm_radial(basis, norm) 464 465 TYPE(gto_basis_set_type), POINTER :: basis 466 REAL(dp), DIMENSION(:), POINTER :: norm 467 468 CHARACTER(len=*), PARAMETER :: routineN = 'basis_norm_radial', & 469 routineP = moduleN//':'//routineN 470 471 INTEGER :: ipgf, iset, isgf, ishell, jpgf, l, nbas 472 REAL(KIND=dp) :: aai, aaj, cci, ccj, expa, ppl 473 474 NULLIFY (norm) 475 476 nbas = basis%nsgf 477 ALLOCATE (norm(nbas)) 478 norm = 0._dp 479 480 DO iset = 1, basis%nset 481 DO ishell = 1, basis%nshell(iset) 482 l = basis%l(ishell, iset) 483 expa = 0.5_dp*REAL(2*l + 3, dp) 484 ppl = fac(2*l + 2)*SQRT(pi)/2._dp**REAL(2*l + 3, dp)/fac(l + 1) 485 DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset) 486 DO ipgf = 1, basis%npgf(iset) 487 cci = basis%gcc(ipgf, ishell, iset) 488 aai = basis%zet(ipgf, iset) 489 DO jpgf = 1, basis%npgf(iset) 490 ccj = basis%gcc(jpgf, ishell, iset) 491 aaj = basis%zet(jpgf, iset) 492 norm(isgf) = norm(isgf) + cci*ccj*ppl/(aai + aaj)**expa 493 END DO 494 END DO 495 norm(isgf) = 1.0_dp/SQRT(norm(isgf)) 496 END DO 497 END DO 498 END DO 499 500 END SUBROUTINE basis_norm_radial 501 502!***************************************************************************** 503!> \brief integral over a single (contracted) lri auxiliary basis function, 504!> integral is zero for all but s-functions 505!> \param basis ... 506!> \param int_aux ... 507!> \param norm ... 508! ************************************************************************************************** 509 SUBROUTINE basis_int(basis, int_aux, norm) 510 511 TYPE(gto_basis_set_type), POINTER :: basis 512 REAL(dp), DIMENSION(:), POINTER :: int_aux, norm 513 514 CHARACTER(len=*), PARAMETER :: routineN = 'basis_int', routineP = moduleN//':'//routineN 515 516 INTEGER :: ipgf, iset, isgf, ishell, l, nbas 517 REAL(KIND=dp) :: aa, cc, pp 518 519 nbas = basis%nsgf 520 ALLOCATE (int_aux(nbas)) 521 int_aux = 0._dp 522 523 DO iset = 1, basis%nset 524 DO ishell = 1, basis%nshell(iset) 525 l = basis%l(ishell, iset) 526 IF (l /= 0) CYCLE 527 DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset) 528 DO ipgf = 1, basis%npgf(iset) 529 cc = basis%gcc(ipgf, ishell, iset) 530 aa = basis%zet(ipgf, iset) 531 pp = (pi/aa)**(3._dp/2._dp) 532 int_aux(isgf) = int_aux(isgf) + norm(isgf)*cc*pp 533 END DO 534 END DO 535 END DO 536 END DO 537 538 END SUBROUTINE basis_int 539 540!***************************************************************************** 541!> \brief self-overlap of lri basis for contracted spherical Gaussians. 542!> Overlap of radial part. Norm contains only normalization of radial 543!> part. Norm and overlap of spherical harmonics not explicitly 544!> calculated since this cancels for the self-overlap anyway. 545!> \param basis ... 546!> \param ovlp ... 547!> \param norm ... 548! ************************************************************************************************** 549 SUBROUTINE basis_ovlp(basis, ovlp, norm) 550 551 TYPE(gto_basis_set_type), POINTER :: basis 552 REAL(dp), DIMENSION(:, :), POINTER :: ovlp 553 REAL(dp), DIMENSION(:), POINTER :: norm 554 555 CHARACTER(len=*), PARAMETER :: routineN = 'basis_ovlp', routineP = moduleN//':'//routineN 556 557 INTEGER :: ipgf, iset, isgf, ishell, jpgf, jset, & 558 jsgf, jshell, l, li, lj, m_i, m_j, nbas 559 REAL(KIND=dp) :: aai, aaj, cci, ccj, expa, norm_i, & 560 norm_j, oo, ppl 561 562 nbas = basis%nsgf 563 ALLOCATE (ovlp(nbas, nbas)) 564 ovlp = 0._dp 565 566 DO iset = 1, basis%nset 567 DO ishell = 1, basis%nshell(iset) 568 li = basis%l(ishell, iset) 569 DO jset = 1, basis%nset 570 DO jshell = 1, basis%nshell(jset) 571 lj = basis%l(jshell, jset) 572 IF (li == lj) THEN 573 l = li 574 expa = 0.5_dp*REAL(2*l + 3, dp) 575 ppl = fac(2*l + 2)*SQRT(pi)/2._dp**REAL(2*l + 3, dp)/fac(l + 1) 576 DO isgf = basis%first_sgf(ishell, iset), basis%last_sgf(ishell, iset) 577 m_i = basis%m(isgf) 578 DO jsgf = basis%first_sgf(jshell, jset), basis%last_sgf(jshell, jset) 579 m_j = basis%m(jsgf) 580 IF (m_i == m_j) THEN 581 DO ipgf = 1, basis%npgf(iset) 582 cci = basis%gcc(ipgf, ishell, iset) 583 aai = basis%zet(ipgf, iset) 584 norm_i = norm(isgf) 585 DO jpgf = 1, basis%npgf(jset) 586 ccj = basis%gcc(jpgf, jshell, jset) 587 aaj = basis%zet(jpgf, jset) 588 oo = 1._dp/(aai + aaj)**expa 589 norm_j = norm(jsgf) 590 ovlp(isgf, jsgf) = ovlp(isgf, jsgf) + norm_i*norm_j*ppl*cci*ccj*oo 591 END DO 592 END DO 593 ENDIF 594 END DO 595 END DO 596 END IF 597 END DO 598 END DO 599 END DO 600 END DO 601 602 END SUBROUTINE basis_ovlp 603 604END MODULE lri_environment_init 605