1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Initialization for solid harmonic Gaussian (SHG) integral scheme. Scheme for calculation 8!> of contracted, spherical Gaussian integrals using the solid harmonics. Initialization of 9!> the contraction matrices 10!> \par History 11!> created [08.2016] 12!> \author Dorothea Golze 13! ************************************************************************************************** 14MODULE generic_shg_integrals_init 15 USE basis_set_types, ONLY: gto_basis_set_type 16 USE kinds, ONLY: dp 17 USE mathconstants, ONLY: fac,& 18 ifac,& 19 pi 20 USE memory_utilities, ONLY: reallocate 21 USE orbital_pointers, ONLY: indso,& 22 nsoset 23 USE spherical_harmonics, ONLY: clebsch_gordon,& 24 clebsch_gordon_deallocate,& 25 clebsch_gordon_init 26#include "../base/base_uses.f90" 27 28 IMPLICIT NONE 29 30 PRIVATE 31 32 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_shg_integrals_init' 33 34 PUBLIC :: contraction_matrix_shg, contraction_matrix_shg_mix, contraction_matrix_shg_rx2m, & 35 get_clebsch_gordon_coefficients 36 37! ************************************************************************************************** 38 39CONTAINS 40 41! ************************************************************************************************** 42!> \brief contraction matrix for SHG integrals 43!> \param basis ... 44!> \param scon_shg contraction matrix 45! ************************************************************************************************** 46 SUBROUTINE contraction_matrix_shg(basis, scon_shg) 47 48 TYPE(gto_basis_set_type), POINTER :: basis 49 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: scon_shg 50 51 CHARACTER(len=*), PARAMETER :: routineN = 'contraction_matrix_shg', & 52 routineP = moduleN//':'//routineN 53 54 INTEGER :: ipgf, iset, ishell, l, maxpgf, maxshell, & 55 nset 56 INTEGER, DIMENSION(:), POINTER :: npgf, nshell 57 REAL(KIND=dp) :: aif, gcc 58 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: norm 59 REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet 60 61 nset = basis%nset 62 npgf => basis%npgf 63 nshell => basis%nshell 64 zet => basis%zet 65 66 maxpgf = SIZE(basis%gcc, 1) 67 maxshell = SIZE(basis%gcc, 2) 68 ALLOCATE (norm(basis%nset, maxshell)) 69 ALLOCATE (scon_shg(maxpgf, maxshell, nset)) 70 scon_shg = 0.0_dp 71 72 CALL basis_norm_shg(basis, norm) 73 74 DO iset = 1, nset 75 DO ishell = 1, nshell(iset) 76 l = basis%l(ishell, iset) 77 DO ipgf = 1, npgf(iset) 78 aif = 1.0_dp/((2._dp*zet(ipgf, iset))**l) 79 gcc = basis%gcc(ipgf, ishell, iset) 80 scon_shg(ipgf, ishell, iset) = norm(iset, ishell)*gcc*aif 81 END DO 82 END DO 83 END DO 84 85 DEALLOCATE (norm) 86 87 END SUBROUTINE contraction_matrix_shg 88 89!*************************************************************************************************** 90!> \brief normalization solid harmonic Gaussians (SHG) 91!> \param basis ... 92!> \param norm ... 93! ************************************************************************************************** 94 SUBROUTINE basis_norm_shg(basis, norm) 95 96 TYPE(gto_basis_set_type), POINTER :: basis 97 REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: norm 98 99 CHARACTER(len=*), PARAMETER :: routineN = 'basis_norm_shg', routineP = moduleN//':'//routineN 100 101 INTEGER :: ipgf, iset, ishell, jpgf, l 102 REAL(KIND=dp) :: aai, aaj, cci, ccj, expa, ppl 103 104 norm = 0._dp 105 106 DO iset = 1, basis%nset 107 DO ishell = 1, basis%nshell(iset) 108 l = basis%l(ishell, iset) 109 expa = 0.5_dp*REAL(2*l + 3, dp) 110 ppl = fac(2*l + 2)*pi**(1.5_dp)/fac(l + 1) 111 ppl = ppl/(2._dp**REAL(2*l + 1, dp)) 112 ppl = ppl/REAL(2*l + 1, dp) 113 DO ipgf = 1, basis%npgf(iset) 114 cci = basis%gcc(ipgf, ishell, iset) 115 aai = basis%zet(ipgf, iset) 116 DO jpgf = 1, basis%npgf(iset) 117 ccj = basis%gcc(jpgf, ishell, iset) 118 aaj = basis%zet(jpgf, iset) 119 norm(iset, ishell) = norm(iset, ishell) + cci*ccj*ppl/(aai + aaj)**expa 120 END DO 121 END DO 122 norm(iset, ishell) = 1.0_dp/SQRT(norm(iset, ishell)) 123 END DO 124 END DO 125 126 END SUBROUTINE basis_norm_shg 127 128! ************************************************************************************************** 129!> \brief mixed contraction matrix for SHG integrals [aba] and [abb] for orbital and ri basis 130!> at the same atom 131!> \param orb_basis orbital basis 132!> \param ri_basis ... 133!> \param orb_index index for orbital basis 134!> \param ri_index index for ri basis 135!> \param scon_mix mixed contraction matrix 136! ************************************************************************************************** 137 SUBROUTINE contraction_matrix_shg_mix(orb_basis, ri_basis, orb_index, ri_index, scon_mix) 138 139 TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis 140 INTEGER, DIMENSION(:, :, :), POINTER :: orb_index, ri_index 141 REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: scon_mix 142 143 CHARACTER(len=*), PARAMETER :: routineN = 'contraction_matrix_shg_mix', & 144 routineP = moduleN//':'//routineN 145 146 INTEGER :: forb, fri, iil, il, ipgf, iset, ishell, jpgf, jset, jshell, l, l1, l2, lmax_orb, & 147 lmax_ri, maxpgf_orb, maxpgf_ri, maxshell_orb, maxshell_ri, nf_orb, nf_ri, nl, nl_max, & 148 nset_orb, nset_ri 149 INTEGER, DIMENSION(:), POINTER :: npgf_orb, npgf_ri, nshell_orb, nshell_ri 150 REAL(KIND=dp) :: cjf, const, const1, const2, gcc_orb, & 151 gcc_ri, prefac, scon1, scon2, zet 152 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: shg_fac 153 REAL(KIND=dp), DIMENSION(:, :), POINTER :: norm_orb, norm_ri, zet_orb, zet_ri 154 155 nset_orb = orb_basis%nset 156 npgf_orb => orb_basis%npgf 157 nshell_orb => orb_basis%nshell 158 zet_orb => orb_basis%zet 159 nset_ri = ri_basis%nset 160 npgf_ri => ri_basis%npgf 161 nshell_ri => ri_basis%nshell 162 zet_ri => ri_basis%zet 163 164 maxpgf_orb = SIZE(orb_basis%gcc, 1) 165 maxshell_orb = SIZE(orb_basis%gcc, 2) 166 ALLOCATE (norm_orb(nset_orb, maxshell_orb)) 167 maxpgf_ri = SIZE(ri_basis%gcc, 1) 168 maxshell_ri = SIZE(ri_basis%gcc, 2) 169 ALLOCATE (norm_ri(nset_ri, maxshell_ri)) 170 171 CALL basis_norm_shg(orb_basis, norm_orb) 172 CALL basis_norm_shg(ri_basis, norm_ri) 173 174 ALLOCATE (orb_index(maxpgf_orb, maxshell_orb, nset_orb)) 175 ALLOCATE (ri_index(maxpgf_ri, maxshell_ri, nset_ri)) 176 177 !** index orbital basis set 178 nf_orb = 0 179 DO iset = 1, nset_orb 180 DO ishell = 1, nshell_orb(iset) 181 DO ipgf = 1, npgf_orb(iset) 182 nf_orb = nf_orb + 1 183 orb_index(ipgf, ishell, iset) = nf_orb 184 END DO 185 END DO 186 END DO 187 188 !** index ri basis set 189 nf_ri = 0 190 DO iset = 1, nset_ri 191 DO ishell = 1, nshell_ri(iset) 192 DO ipgf = 1, npgf_ri(iset) 193 nf_ri = nf_ri + 1 194 ri_index(ipgf, ishell, iset) = nf_ri 195 END DO 196 END DO 197 END DO 198 199 lmax_orb = MAXVAL(orb_basis%lmax) 200 lmax_ri = MAXVAL(ri_basis%lmax) 201 nl_max = INT((lmax_orb + lmax_ri)/2) + 1 202 ALLOCATE (scon_mix(nl_max, nf_ri, nf_orb, nl_max)) 203 scon_mix = 0.0_dp 204 205 ALLOCATE (shg_fac(0:nl_max - 1)) 206 shg_fac(0) = 1.0_dp 207 208 DO iset = 1, nset_orb 209 DO ishell = 1, nshell_orb(iset) 210 l1 = orb_basis%l(ishell, iset) 211 const1 = SQRT(1.0_dp/REAL(2*l1 + 1, dp)) 212 DO jset = 1, nset_ri 213 DO jshell = 1, nshell_ri(jset) 214 l2 = ri_basis%l(jshell, jset) 215 const2 = SQRT(1.0_dp/REAL(2*l2 + 1, dp)) 216 nl = INT((l1 + l2)/2) 217 IF (l1 == 0 .OR. l2 == 0) nl = 0 218 DO il = 0, nl 219 l = l1 + l2 - 2*il 220 const = const1*const2*2.0_dp*SQRT(pi*REAL(2*l + 1, dp)) 221 DO iil = 1, il 222 shg_fac(iil) = fac(l + iil - 1)*ifac(l)*REAL(l, dp) & 223 *fac(il)/fac(il - iil)/fac(iil) 224 ENDDO 225 DO ipgf = 1, npgf_orb(iset) 226 forb = orb_index(ipgf, ishell, iset) 227 gcc_orb = orb_basis%gcc(ipgf, ishell, iset) 228 scon1 = norm_orb(iset, ishell)*gcc_orb 229 DO jpgf = 1, npgf_ri(jset) 230 fri = ri_index(jpgf, jshell, jset) 231 gcc_ri = ri_basis%gcc(jpgf, jshell, jset) 232 scon2 = norm_ri(jset, jshell)*gcc_ri 233 zet = zet_orb(ipgf, iset) + zet_ri(jpgf, jset) 234 cjf = 1.0_dp/((2._dp*zet)**l) 235 prefac = const*cjf*scon1*scon2 236 DO iil = 0, il 237 scon_mix(iil + 1, fri, forb, il + 1) = prefac*shg_fac(iil)/zet**iil 238 ENDDO 239 ENDDO 240 ENDDO 241 ENDDO 242 ENDDO 243 ENDDO 244 ENDDO 245 ENDDO 246 247 DEALLOCATE (norm_orb, norm_ri, shg_fac) 248 249 END SUBROUTINE contraction_matrix_shg_mix 250 251! ************************************************************************************************** 252!> \brief ... 253!> \param basis ... 254!> \param m ... 255!> \param scon_shg ... 256!> \param scon_rx2m ... 257! ************************************************************************************************** 258 SUBROUTINE contraction_matrix_shg_rx2m(basis, m, scon_shg, scon_rx2m) 259 260 TYPE(gto_basis_set_type), POINTER :: basis 261 INTEGER, INTENT(IN) :: m 262 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: scon_shg 263 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: scon_rx2m 264 265 CHARACTER(len=*), PARAMETER :: routineN = 'contraction_matrix_shg_rx2m', & 266 routineP = moduleN//':'//routineN 267 268 INTEGER :: ipgf, iset, ishell, j, l, maxpgf, & 269 maxshell, nset 270 INTEGER, DIMENSION(:), POINTER :: npgf, nshell 271 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: shg_fac 272 REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet 273 274 npgf => basis%npgf 275 nshell => basis%nshell 276 zet => basis%zet 277 nset = basis%nset 278 279 maxpgf = SIZE(basis%gcc, 1) 280 maxshell = SIZE(basis%gcc, 2) 281 ALLOCATE (scon_rx2m(maxpgf, m + 1, maxshell, nset)) 282 scon_rx2m = 0.0_dp 283 ALLOCATE (shg_fac(0:m)) 284 shg_fac(0) = 1.0_dp 285 286 DO iset = 1, nset 287 DO ishell = 1, nshell(iset) 288 l = basis%l(ishell, iset) 289 DO j = 1, m 290 shg_fac(j) = fac(l + j - 1)*ifac(l)*REAL(l, dp) & 291 *fac(m)/fac(m - j)/fac(j) 292 ENDDO 293 DO ipgf = 1, npgf(iset) 294 DO j = 0, m 295 scon_rx2m(ipgf, j + 1, ishell, iset) = scon_shg(ipgf, ishell, iset) & 296 *shg_fac(j)/zet(ipgf, iset)**j 297 ENDDO 298 ENDDO 299 ENDDO 300 ENDDO 301 302 DEALLOCATE (shg_fac) 303 304 END SUBROUTINE contraction_matrix_shg_rx2m 305 306! ************************************************************************************************** 307!> \brief calculate the Clebsch-Gordon (CG) coefficients for expansion of the 308!> product of two spherical harmonic Gaussians 309!> \param my_cg matrix storing CG coefficients 310!> \param cg_none0_list list of none-zero CG coefficients 311!> \param ncg_none0 number of none-zero CG coefficients 312!> \param maxl1 maximal l quantum number of 1st spherical function 313!> \param maxl2 maximal l quantum number of 2nd spherical function 314! ************************************************************************************************** 315 SUBROUTINE get_clebsch_gordon_coefficients(my_cg, cg_none0_list, ncg_none0, & 316 maxl1, maxl2) 317 318 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_cg 319 INTEGER, DIMENSION(:, :, :), POINTER :: cg_none0_list 320 INTEGER, DIMENSION(:, :), POINTER :: ncg_none0 321 INTEGER, INTENT(IN) :: maxl1, maxl2 322 323 CHARACTER(len=*), PARAMETER :: routineN = 'get_clebsch_gordon_coefficients', & 324 routineP = moduleN//':'//routineN 325 326 INTEGER :: il, ilist, iso, iso1, iso2, l1, l1l2, & 327 l2, lc1, lc2, lp, m1, m2, maxl, mm, & 328 mp, nlist, nlist_max, nsfunc, nsfunc1, & 329 nsfunc2 330 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rga 331 332 nlist_max = 6 333 nsfunc1 = nsoset(maxl1) 334 nsfunc2 = nsoset(maxl2) 335 maxl = maxl1 + maxl2 336 nsfunc = nsoset(maxl) 337 338 CALL clebsch_gordon_init(maxl) 339 340 ALLOCATE (my_cg(nsfunc1, nsfunc2, nsfunc)) 341 my_cg = 0.0_dp 342 ALLOCATE (ncg_none0(nsfunc1, nsfunc2)) 343 ncg_none0 = 0 344 ALLOCATE (cg_none0_list(nsfunc1, nsfunc2, nlist_max)) 345 cg_none0_list = 0 346 347 ALLOCATE (rga(maxl, 2)) 348 rga = 0.0_dp 349 DO lc1 = 0, maxl1 350 DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1) 351 l1 = indso(1, iso1) 352 m1 = indso(2, iso1) 353 DO lc2 = 0, maxl2 354 DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2) 355 l2 = indso(1, iso2) 356 m2 = indso(2, iso2) 357 CALL clebsch_gordon(l1, m1, l2, m2, rga) 358 l1l2 = l1 + l2 359 mp = m1 + m2 360 mm = m1 - m2 361 IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN 362 mp = -ABS(mp) 363 mm = -ABS(mm) 364 ELSE 365 mp = ABS(mp) 366 mm = ABS(mm) 367 END IF 368 DO lp = MOD(l1 + l2, 2), l1l2, 2 369 il = lp/2 + 1 370 IF (ABS(mp) <= lp) THEN 371 IF (mp >= 0) THEN 372 iso = nsoset(lp - 1) + lp + 1 + mp 373 ELSE 374 iso = nsoset(lp - 1) + lp + 1 - ABS(mp) 375 END IF 376 my_cg(iso1, iso2, iso) = rga(il, 1) 377 ENDIF 378 IF (mp /= mm .AND. ABS(mm) <= lp) THEN 379 IF (mm >= 0) THEN 380 iso = nsoset(lp - 1) + lp + 1 + mm 381 ELSE 382 iso = nsoset(lp - 1) + lp + 1 - ABS(mm) 383 END IF 384 my_cg(iso1, iso2, iso) = rga(il, 2) 385 ENDIF 386 END DO 387 nlist = 0 388 DO ilist = 1, nsfunc 389 IF (ABS(my_cg(iso1, iso2, ilist)) > 1.E-8_dp) THEN 390 nlist = nlist + 1 391 IF (nlist > nlist_max) THEN 392 CALL reallocate(cg_none0_list, 1, nsfunc1, 1, nsfunc2, 1, nlist) 393 nlist_max = nlist 394 ENDIF 395 cg_none0_list(iso1, iso2, nlist) = ilist 396 ENDIF 397 ENDDO 398 ncg_none0(iso1, iso2) = nlist 399 ENDDO ! iso2 400 ENDDO ! lc2 401 ENDDO ! iso1 402 ENDDO ! lc1 403 404 DEALLOCATE (rga) 405 CALL clebsch_gordon_deallocate() 406 407 END SUBROUTINE get_clebsch_gordon_coefficients 408 409END MODULE generic_shg_integrals_init 410