1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Ewald sums to represent integrals in direct and reciprocal lattice. 8!> \par History 9!> 2015 09 created 10!> \author Patrick Seewald 11! ************************************************************************************************** 12 13MODULE eri_mme_lattice_summation 14#:include "eri_mme_unroll.fypp" 15 16 USE ao_util, ONLY: exp_radius 17 USE eri_mme_gaussian, ONLY: create_gaussian_overlap_dist_to_hermite, & 18 create_hermite_to_cartesian, & 19 eri_mme_coulomb, eri_mme_yukawa, & 20 eri_mme_longrange 21 USE kinds, ONLY: dp, & 22 int_8 23 USE mathconstants, ONLY: gaussi, & 24 pi, & 25 twopi 26 USE orbital_pointers, ONLY: coset, & 27 indco, & 28 ncoset 29#include "../base/base_uses.f90" 30 31 IMPLICIT NONE 32 33 PRIVATE 34 35 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE. 36 37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'eri_mme_lattice_summation' 38 39 PUBLIC :: & 40 ellipsoid_bounds, & 41 eri_mme_2c_get_bounds, & 42 eri_mme_2c_get_rads, & 43 eri_mme_3c_get_bounds, & 44 eri_mme_3c_get_rads, & 45 get_l, & 46 pgf_sum_2c_gspace_1d, & 47 pgf_sum_2c_gspace_1d_deltal, & 48 pgf_sum_2c_gspace_3d, & 49 pgf_sum_2c_rspace_1d, & 50 pgf_sum_2c_rspace_3d, & 51 pgf_sum_3c_1d, & 52 pgf_sum_3c_3d 53 54 INTEGER, PARAMETER, PRIVATE :: exp_w = 50, div_w = 10 55 56CONTAINS 57 58! ************************************************************************************************** 59!> \brief Get summation radii for 2c integrals 60!> \param la_max ... 61!> \param lb_max ... 62!> \param zeta ... 63!> \param zetb ... 64!> \param a_mm ... 65!> \param G_min ... 66!> \param R_min ... 67!> \param sum_precision ... 68!> \param G_rad ... 69!> \param R_rad ... 70! ************************************************************************************************** 71 SUBROUTINE eri_mme_2c_get_rads(la_max, lb_max, zeta, zetb, a_mm, G_min, R_min, sum_precision, G_rad, R_rad) 72 INTEGER, INTENT(IN) :: la_max, lb_max 73 REAL(KIND=dp), INTENT(IN) :: zeta, zetb, a_mm, G_min, R_min, & 74 sum_precision 75 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: G_rad, R_rad 76 77 INTEGER :: l_max 78 REAL(KIND=dp) :: alpha_G, alpha_R, G_res, R_res 79 80 l_max = la_max + lb_max 81 alpha_G = a_mm + 0.25_dp/zeta + 0.25_dp/zetb 82 alpha_R = 0.25_dp/alpha_G 83 84 G_res = 0.5_dp*G_min 85 R_res = 0.5_dp*R_min 86 87 IF (PRESENT(G_rad)) G_rad = exp_radius(l_max, alpha_G, sum_precision, 1.0_dp, epsabs=G_res) 88 IF (PRESENT(R_rad)) R_rad = exp_radius(l_max, alpha_R, sum_precision, 1.0_dp, epsabs=R_res) 89 90 END SUBROUTINE 91 92! ************************************************************************************************** 93!> \brief Get summation radii for 3c integrals 94!> \param la_max ... 95!> \param lb_max ... 96!> \param lc_max ... 97!> \param zeta ... 98!> \param zetb ... 99!> \param zetc ... 100!> \param a_mm ... 101!> \param G_min ... 102!> \param R_min ... 103!> \param sum_precision ... 104!> \param G_rads_1 ... 105!> \param R_rads_2 ... 106!> \param R_rads_3 ... 107! ************************************************************************************************** 108 SUBROUTINE eri_mme_3c_get_rads(la_max, lb_max, lc_max, zeta, zetb, zetc, a_mm, G_min, R_min, & 109 sum_precision, G_rads_1, R_rads_2, R_rads_3) 110 INTEGER, INTENT(IN) :: la_max, lb_max, lc_max 111 REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc, a_mm 112 REAL(KIND=dp), INTENT(IN), OPTIONAL :: G_min, R_min 113 REAL(KIND=dp), INTENT(IN) :: sum_precision 114 REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: G_rads_1, R_rads_2 115 REAL(KIND=dp), DIMENSION(2), INTENT(OUT), OPTIONAL :: R_rads_3 116 117 CHARACTER(LEN=*), PARAMETER :: routineN = 'eri_mme_3c_get_rads', & 118 routineP = moduleN//':'//routineN 119 120 REAL(KIND=dp) :: alpha, alpha_R, beta, G_res, gamma, R_res 121 122 ! exponents in G space 123 alpha = 0.25_dp/zeta 124 beta = 0.25_dp/zetb 125 gamma = 0.25_dp/zetc + a_mm 126 127 ! Summation radii and number of summands for all lattice summation methods 128 ! sum method 1 129 IF (PRESENT(G_rads_1)) THEN 130 G_res = 0.5_dp*G_min 131 G_rads_1(1) = exp_radius(la_max, alpha, sum_precision, 1.0_dp, G_res) 132 G_rads_1(2) = exp_radius(lb_max, beta, sum_precision, 1.0_dp, G_res) 133 G_rads_1(3) = exp_radius(lc_max, gamma, sum_precision, 1.0_dp, G_res) 134 ENDIF 135 136 ! sum method 2 137 IF (PRESENT(R_rads_2)) THEN 138 R_res = 0.5_dp*R_min 139 R_rads_2(1) = exp_radius(lb_max + lc_max, 0.25_dp/(beta + gamma), sum_precision, 1.0_dp, epsabs=R_res) 140 R_rads_2(2) = exp_radius(lc_max + la_max, 0.25_dp/(alpha + gamma), sum_precision, 1.0_dp, epsabs=R_res) 141 R_rads_2(3) = exp_radius(lb_max + la_max, 0.25_dp/(alpha + beta), sum_precision, 1.0_dp, epsabs=R_res) 142 ENDIF 143 144 ! sum method 3 145 146 IF (PRESENT(R_rads_3)) THEN 147 R_res = 0.5_dp*R_min 148 alpha_R = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm) 149 R_rads_3(1) = exp_radius(la_max + lb_max, zeta*zetb/(zeta + zetb), sum_precision, 1.0_dp, R_res) 150 R_rads_3(2) = exp_radius(la_max + lb_max + lc_max, alpha_R, sum_precision, 1.0_dp, R_res) 151 ENDIF 152 153 END SUBROUTINE 154 155! ************************************************************************************************** 156!> \brief Get summation bounds for 2c integrals 157!> \param hmat ... 158!> \param h_inv ... 159!> \param vol ... 160!> \param is_ortho ... 161!> \param G_min ... 162!> \param R_min ... 163!> \param la_max ... 164!> \param lb_max ... 165!> \param zeta ... 166!> \param zetb ... 167!> \param a_mm ... 168!> \param sum_precision ... 169!> \param n_sum_1d ... 170!> \param n_sum_3d ... 171!> \param G_bounds ... 172!> \param G_rad ... 173!> \param R_bounds ... 174!> \param R_rad ... 175! ************************************************************************************************** 176 SUBROUTINE eri_mme_2c_get_bounds(hmat, h_inv, vol, is_ortho, G_min, R_min, la_max, lb_max, & 177 zeta, zetb, a_mm, sum_precision, n_sum_1d, n_sum_3d, & 178 G_bounds, G_rad, R_bounds, R_rad) 179 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv 180 REAL(KIND=dp), INTENT(IN) :: vol 181 LOGICAL, INTENT(IN) :: is_ortho 182 REAL(KIND=dp), INTENT(IN) :: G_min, R_min 183 INTEGER, INTENT(IN) :: la_max, lb_max 184 REAL(KIND=dp), INTENT(IN) :: zeta, zetb, a_mm, sum_precision 185 INTEGER(KIND=int_8), DIMENSION(2, 3), INTENT(OUT) :: n_sum_1d 186 INTEGER(KIND=int_8), DIMENSION(2), INTENT(OUT) :: n_sum_3d 187 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: G_bounds 188 REAL(KIND=dp), INTENT(OUT) :: G_rad 189 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: R_bounds 190 REAL(KIND=dp), INTENT(OUT) :: R_rad 191 192 INTEGER :: i_xyz 193 REAL(KIND=dp) :: ns_G, ns_R, vol_G 194 195 CALL eri_mme_2c_get_rads(la_max, lb_max, zeta, zetb, a_mm, G_min, R_min, sum_precision, G_rad, R_rad) 196 197 G_bounds = ellipsoid_bounds(G_rad, TRANSPOSE(hmat)/(2.0_dp*pi)) 198 R_bounds = ellipsoid_bounds(R_rad, h_inv) 199 200 vol_G = twopi**3/vol 201 202 IF (is_ortho) THEN 203 DO i_xyz = 1, 3 204 ns_G = 2.0_dp*G_bounds(i_xyz) 205 ns_R = 2.0_dp*R_bounds(i_xyz) 206 n_sum_1d(1, i_xyz) = nsum_2c_gspace_1d(ns_G, la_max, lb_max) 207 n_sum_1d(2, i_xyz) = nsum_2c_rspace_1d(ns_R, la_max, lb_max) 208 ENDDO 209 ELSE 210 ns_G = 4._dp/3._dp*pi*G_rad**3/vol_G 211 ns_R = 4._dp/3._dp*pi*R_rad**3/vol 212 n_sum_3d(1) = nsum_2c_gspace_3d(ns_G, la_max, lb_max) 213 n_sum_3d(2) = nsum_2c_rspace_3d(ns_R, la_max, lb_max) 214 ENDIF 215 216 END SUBROUTINE 217 218! ************************************************************************************************** 219!> \brief Get summation bounds for 3c integrals 220!> \param hmat ... 221!> \param h_inv ... 222!> \param vol ... 223!> \param is_ortho ... 224!> \param G_min ... 225!> \param R_min ... 226!> \param la_max ... 227!> \param lb_max ... 228!> \param lc_max ... 229!> \param zeta ... 230!> \param zetb ... 231!> \param zetc ... 232!> \param a_mm ... 233!> \param sum_precision ... 234!> \param n_sum_1d ... 235!> \param n_sum_3d ... 236!> \param G_bounds_1 ... 237!> \param G_rads_1 ... 238!> \param R_bounds_2 ... 239!> \param R_rads_2 ... 240!> \param R_bounds_3 ... 241!> \param R_rads_3 ... 242! ************************************************************************************************** 243 SUBROUTINE eri_mme_3c_get_bounds(hmat, h_inv, vol, is_ortho, G_min, R_min, la_max, lb_max, lc_max, & 244 zeta, zetb, zetc, a_mm, sum_precision, n_sum_1d, n_sum_3d, & 245 G_bounds_1, G_rads_1, R_bounds_2, R_rads_2, R_bounds_3, R_rads_3) 246 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv 247 REAL(KIND=dp), INTENT(IN) :: vol 248 LOGICAL, INTENT(IN) :: is_ortho 249 REAL(KIND=dp), INTENT(IN) :: G_min, R_min 250 INTEGER, INTENT(IN) :: la_max, lb_max, lc_max 251 REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc, a_mm, sum_precision 252 INTEGER(KIND=int_8), DIMENSION(3, 3), INTENT(OUT) :: n_sum_1d 253 INTEGER(KIND=int_8), DIMENSION(3), INTENT(OUT) :: n_sum_3d 254 REAL(KIND=dp), DIMENSION(3, 3) :: G_bounds_1 255 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: G_rads_1 256 REAL(KIND=dp), DIMENSION(3, 3) :: R_bounds_2 257 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: R_rads_2 258 REAL(KIND=dp), DIMENSION(2, 3) :: R_bounds_3 259 REAL(KIND=dp), DIMENSION(2), INTENT(OUT) :: R_rads_3 260 261 INTEGER :: i, i_xyz, order_1, order_2 262 REAL(KIND=dp) :: ns1_G1, ns1_G2, ns2_G, ns2_R, ns3_R1, & 263 ns3_R2, vol_G 264 CALL eri_mme_3c_get_rads(la_max, lb_max, lc_max, zeta, zetb, zetc, a_mm, G_min, R_min, sum_precision, & 265 G_rads_1=G_rads_1, R_rads_2=R_rads_2, R_rads_3=R_rads_3) 266 267 vol_G = twopi**3/vol 268 269 order_1 = MAXLOC(G_rads_1, DIM=1) 270 order_2 = MINLOC(G_rads_1, DIM=1) 271 272 DO i = 1, 3 273 G_bounds_1(i, :) = ellipsoid_bounds(G_rads_1(i), TRANSPOSE(hmat)/(2.0_dp*pi)) 274 ENDDO 275 276 DO i = 1, 3 277 R_bounds_2(i, :) = ellipsoid_bounds(R_rads_2(i), h_inv) 278 ENDDO 279 280 DO i = 1, 2 281 R_bounds_3(i, :) = ellipsoid_bounds(R_rads_3(i), h_inv) 282 ENDDO 283 284 IF (is_ortho) THEN 285 DO i_xyz = 1, 3 286 287 ns3_R1 = 2.0_dp*R_bounds_3(1, i_xyz) 288 ns3_R2 = 2.0_dp*R_bounds_3(2, i_xyz) 289 290 n_sum_1d(3, i_xyz) = nsum_3c_rspace_1d(ns3_R1, ns3_R2) 291 ENDDO 292 293 ELSE 294 295 order_1 = MAXLOC(G_rads_1, DIM=1) 296 order_2 = MINLOC(G_rads_1, DIM=1) 297 298 SELECT CASE (order_1) 299 CASE (1) 300 ns1_G1 = 4._dp/3._dp*pi*G_rads_1(2)**3/vol_G 301 ns1_G2 = 4._dp/3._dp*pi*G_rads_1(3)**3/vol_G 302 CASE (2) 303 ns1_G1 = 4._dp/3._dp*pi*G_rads_1(1)**3/vol_G 304 ns1_G2 = 4._dp/3._dp*pi*G_rads_1(3)**3/vol_G 305 CASE (3) 306 ns1_G1 = 4._dp/3._dp*pi*G_rads_1(1)**3/vol_G 307 ns1_G2 = 4._dp/3._dp*pi*G_rads_1(2)**3/vol_G 308 END SELECT 309 310 n_sum_3d(1) = nsum_3c_gspace_3d(ns1_G1, ns1_G2, la_max, lb_max, lc_max) 311 312 ns2_G = 4._dp/3._dp*pi*G_rads_1(order_2)**3/vol_G 313 ns2_R = 4._dp/3._dp*pi*R_rads_2(order_2)**3/vol 314 315 ns3_R1 = 4._dp/3._dp*pi*R_rads_3(1)**3/vol 316 ns3_R2 = 4._dp/3._dp*pi*R_rads_3(2)**3/vol 317 318 SELECT CASE (order_2) 319 CASE (1) 320 n_sum_3d(2) = nsum_product_3c_gspace_3d(ns2_G, ns2_R, la_max, lb_max, lc_max) 321 CASE (2) 322 n_sum_3d(2) = nsum_product_3c_gspace_3d(ns2_G, ns2_R, lb_max, la_max, lc_max) 323 CASE (3) 324 n_sum_3d(2) = nsum_product_3c_gspace_3d(ns2_G, ns2_R, lc_max, lb_max, la_max) 325 END SELECT 326 327 n_sum_3d(3) = nsum_3c_rspace_3d(ns3_R1, ns3_R2, la_max, lb_max, lc_max) 328 ENDIF 329 330 END SUBROUTINE 331 332! ************************************************************************************************** 333!> \brief Roughly estimated number of floating point operations 334!> \param ns_G ... 335!> \param l ... 336!> \param m ... 337!> \return ... 338! ************************************************************************************************** 339 PURE FUNCTION nsum_2c_gspace_1d(ns_G, l, m) 340 REAL(KIND=dp), INTENT(IN) :: ns_G 341 INTEGER, INTENT(IN) :: l, m 342 INTEGER(KIND=int_8) :: nsum_2c_gspace_1d 343 344 nsum_2c_gspace_1d = NINT(ns_G*(2*exp_w + (l + m + 1)*5), KIND=int_8) 345 END FUNCTION 346 347! ************************************************************************************************** 348!> \brief Compute Ewald-like sum for 2-center ERIs in G space in 1 dimension 349!> S_G(l, alpha) = (-i)^l*inv_lgth*sum_G( C(l, alpha, G) exp(iGR) ), with 350!> C(l, alpha, r) = r^l exp(-alpha*r^2), 351!> dG = inv_lgth*twopi and G = -G_bound*dG, (-G_bound + 1)*dG, ..., G_bound*dG 352!> for all l < = l_max. 353!> \param S_G ... 354!> \param R ... 355!> \param alpha ... 356!> \param inv_lgth ... 357!> \param G_c ... 358!> \note S_G is real. 359! ************************************************************************************************** 360 PURE SUBROUTINE pgf_sum_2c_gspace_1d(S_G, R, alpha, inv_lgth, G_c) 361 REAL(KIND=dp), DIMENSION(0:), INTENT(OUT) :: S_G 362 REAL(KIND=dp), INTENT(IN) :: R, alpha, inv_lgth, G_c 363 364 CHARACTER(LEN=*), PARAMETER :: routineN = 'pgf_sum_2c_gspace_1d', & 365 routineP = moduleN//':'//routineN 366 367 COMPLEX(KIND=dp) :: exp_tot 368 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: S_G_c 369 INTEGER :: gg, l, l_max 370 REAL(KIND=dp) :: dG, G, G_pow_l 371 372 dG = inv_lgth*twopi 373 l_max = UBOUND(S_G, 1) 374 375 ALLOCATE (S_G_c(0:l_max)) 376 S_G_c(:) = 0.0_dp 377 DO gg = -FLOOR(G_c), FLOOR(G_c) 378 G = gg*dG 379 exp_tot = EXP(-alpha*G**2)*EXP(gaussi*G*R) ! cost: 2 exp_w flops 380 G_pow_l = 1.0_dp 381 DO l = 0, l_max 382 S_G_c(l) = S_G_c(l) + G_pow_l*(-1.0_dp)**l*exp_tot ! cost: 4 flops 383 G_pow_l = G_pow_l*G ! cost: 1 flop 384 ENDDO 385 ENDDO 386 387 S_G(:) = REAL(S_G_c(0:l_max)*i_pow((/(l, l=0, l_max)/)))*inv_lgth 388 END SUBROUTINE pgf_sum_2c_gspace_1d 389 390! ************************************************************************************************** 391!> \brief Roughly estimated number of floating point operations 392!> \param ns_G ... 393!> \param l ... 394!> \param m ... 395!> \return ... 396! ************************************************************************************************** 397 PURE FUNCTION nsum_2c_gspace_3d(ns_G, l, m) 398 REAL(KIND=dp), INTENT(IN) :: ns_G 399 INTEGER, INTENT(IN) :: l, m 400 INTEGER(KIND=int_8) :: nsum_2c_gspace_3d 401 402 nsum_2c_gspace_3d = NINT(ns_G*(2*exp_w + ncoset(l + m)*7), KIND=int_8) 403 END FUNCTION 404 405! ************************************************************************************************** 406!> \brief As pgf_sum_2c_gspace_1d but 3d sum required for non-orthorhombic cells 407!> \param S_G ... 408!> \param l_max ... 409!> \param R ... 410!> \param alpha ... 411!> \param h_inv ... 412!> \param G_c ... 413!> \param G_rad ... 414!> \param vol ... 415!> \param coulomb ... 416!> \note MMME Method is not very efficient for non-orthorhombic cells 417! ************************************************************************************************** 418 PURE SUBROUTINE pgf_sum_2c_gspace_3d(S_G, l_max, R, alpha, h_inv, G_c, G_rad, vol, potential, pot_par) 419 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: S_G 420 INTEGER, INTENT(IN) :: l_max 421 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R 422 REAL(KIND=dp), INTENT(IN) :: alpha 423 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: h_inv 424 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G_c 425 REAL(KIND=dp), INTENT(IN) :: G_rad, vol 426 INTEGER, INTENT(IN), OPTIONAL :: potential 427 REAL(KIND=dp), INTENT(IN), OPTIONAL :: pot_par 428 429 COMPLEX(KIND=dp) :: exp_tot 430 INTEGER, DIMENSION(3) :: l_xyz 431 INTEGER :: gx, gy, gz, k, l, lco, lx, ly, lz 432 COMPLEX(KIND=dp), DIMENSION(ncoset(l_max)) :: Ig 433 REAL(KIND=dp) :: G_rads_sq, G_sq 434 REAL(KIND=dp), DIMENSION(3) :: G, G_x, G_y, G_z 435 REAL(KIND=dp), DIMENSION(3, 0:l_max) :: G_pow_l 436 REAL(KIND=dp), DIMENSION(3, 3) :: ht 437 438 ht = twopi*TRANSPOSE(h_inv) 439 Ig(:) = 0.0_dp 440 441 G_rads_sq = G_rad**2 442 443 DO gx = -FLOOR(G_c(1)), FLOOR(G_c(1)) 444 G_x = ht(:, 1)*gx 445 DO gy = -FLOOR(G_c(2)), FLOOR(G_c(2)) 446 G_y = ht(:, 2)*gy 447 DO gz = -FLOOR(G_c(3)), FLOOR(G_c(3)) 448 G_z = ht(:, 3)*gz 449 450 G = G_x + G_y + G_z 451 G_sq = G(1)**2 + G(2)**2 + G(3)**2 452 IF (G_sq .GT. G_rads_sq) CYCLE 453 454 IF (PRESENT(potential)) THEN 455 IF (gx .EQ. 0 .AND. gy .EQ. 0 .AND. gz .EQ. 0) CYCLE 456 ENDIF 457 458 exp_tot = EXP(-alpha*G_sq)*EXP(gaussi*DOT_PRODUCT(G, R)) ! cost: 2 exp_w flops 459 IF (PRESENT(potential)) THEN 460 SELECT CASE (potential) 461 CASE (eri_mme_coulomb) 462 exp_tot = exp_tot/G_sq 463 CASE (eri_mme_yukawa) 464 exp_tot = exp_tot/(G_sq + pot_par**2) 465 !exp_tot = exp_tot/G_sq 466 CASE (eri_mme_longrange) 467 exp_tot = exp_tot/G_sq*EXP(-G_sq/pot_par**2) 468 END SELECT 469 ENDIF 470 DO k = 1, 3 471 G_pow_l(k, 0) = 1.0_dp 472 DO l = 1, l_max 473 G_pow_l(k, l) = G_pow_l(k, l - 1)*G(k) 474 ENDDO 475 ENDDO 476 DO lco = 1, ncoset(l_max) 477 CALL get_l(lco, l, lx, ly, lz) 478 l_xyz = [lx, ly, lz] 479 Ig(coset(lx, ly, lz)) = Ig(coset(lx, ly, lz)) + & ! cost: 7 flops 480 G_pow_l(1, lx)*G_pow_l(2, ly)*G_pow_l(3, lz)* & 481 exp_tot*(-1.0_dp)**l*i_pow(l) 482 ENDDO 483 ENDDO 484 ENDDO 485 ENDDO 486 S_G(:) = REAL(Ig(:), KIND=dp)/vol 487 END SUBROUTINE pgf_sum_2c_gspace_3d 488 489! ************************************************************************************************** 490!> \brief Roughly estimated number of floating point operations 491!> \param ns_R ... 492!> \param l ... 493!> \param m ... 494!> \return ... 495! ************************************************************************************************** 496 PURE FUNCTION nsum_2c_rspace_1d(ns_R, l, m) 497 REAL(KIND=dp), INTENT(IN) :: ns_R 498 INTEGER, INTENT(IN) :: l, m 499 INTEGER(KIND=int_8) :: nsum_2c_rspace_1d 500 501 nsum_2c_rspace_1d = NINT(ns_R*(exp_w + (l + m + 1)*3), KIND=int_8) 502 END FUNCTION 503 504! ************************************************************************************************** 505!> \brief Compute Ewald-like sum for 2-center ERIs in R space in 1 dimension 506!> S_R(l, alpha) = SQRT(alpha/pi) sum_R'( H(l, alpha, R-R') ), 507!> with H(l, alpha, R) = (-d/dR)^l exp(-alpha*R^2), 508!> dR = lgth and R' = -R_min*dR, (-R_min + 1)*dR, ..., R_max*dR, 509!> for all l < = l_max. 510!> \param S_R ... 511!> \param R ... 512!> \param alpha ... 513!> \param lgth ... 514!> \param R_c ... 515!> \note result is equivalent to pgf_sum_2c_gspace_1d with 516!> S_R(l, alpha) = S_G(l, 1/(4*alpha)) 517! ************************************************************************************************** 518 PURE SUBROUTINE pgf_sum_2c_rspace_1d(S_R, R, alpha, lgth, R_c) 519 REAL(KIND=dp), DIMENSION(0:), INTENT(OUT) :: S_R 520 REAL(KIND=dp), INTENT(IN) :: R, alpha, lgth, R_c 521 522 CHARACTER(LEN=*), PARAMETER :: routineN = 'pgf_sum_2c_rspace_1d', & 523 routineP = moduleN//':'//routineN 524 525 INTEGER :: l, l_max, rr 526 REAL(KIND=dp) :: dR, exp_tot, R_pow_l, Rp 527 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_to_c 528 529 dR = lgth 530 l_max = UBOUND(S_R, 1) 531 532 ! 1) compute sum over C(l, alpha, R - R') instead of H(l, alpha, R - R') 533 S_R(:) = 0.0_dp 534 Rp = R - R_c*dR 535 DO rr = CEILING(-R_c - R/dR), FLOOR(R_c - R/dR) 536 Rp = R + rr*dR 537 exp_tot = EXP(-alpha*Rp**2) ! cost: exp_w flops 538 R_pow_l = 1.0_dp 539 DO l = 0, l_max 540 S_R(l) = S_R(l) + R_pow_l*exp_tot ! cost: 2 flops 541 R_pow_l = R_pow_l*Rp ! cost: 1 flop 542 ENDDO 543 ENDDO 544 545 ! 2) C --> H 546 CALL create_hermite_to_cartesian(alpha, l_max, h_to_c) 547 S_R = MATMUL(TRANSPOSE(h_to_c(0:l_max, 0:l_max)), S_R)*SQRT(alpha/pi) 548 END SUBROUTINE pgf_sum_2c_rspace_1d 549 550! ************************************************************************************************** 551!> \brief Roughly estimated number of floating point operations 552!> \param ns_R ... 553!> \param l ... 554!> \param m ... 555!> \return ... 556! ************************************************************************************************** 557 PURE FUNCTION nsum_2c_rspace_3d(ns_R, l, m) 558 REAL(KIND=dp), INTENT(IN) :: ns_R 559 INTEGER, INTENT(IN) :: l, m 560 INTEGER(KIND=int_8) :: nsum_2c_rspace_3d 561 562 nsum_2c_rspace_3d = NINT(ns_R*(exp_w + ncoset(l + m)*(4 + ncoset(l + m)*4)), KIND=int_8) 563 END FUNCTION 564 565! ************************************************************************************************** 566!> \brief As pgf_sum_2c_rspace_1d but 3d sum required for non-orthorhombic cells 567!> \param S_R ... 568!> \param l_max ... 569!> \param R ... 570!> \param alpha ... 571!> \param hmat ... 572!> \param h_inv ... 573!> \param R_c ... 574!> \param R_rad ... 575!> \note MMME Method is not very efficient for non-orthorhombic cells 576! ************************************************************************************************** 577 PURE SUBROUTINE pgf_sum_2c_rspace_3d(S_R, l_max, R, alpha, hmat, h_inv, R_c, R_rad) 578 REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: S_R 579 INTEGER, INTENT(IN) :: l_max 580 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R 581 REAL(KIND=dp), INTENT(IN) :: alpha 582 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv 583 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R_c 584 REAL(KIND=dp), INTENT(IN) :: R_rad 585 586 INTEGER :: k, l, lco, llx, lly, llz, lx, ly, lz, & 587 sx, sy, sz 588 REAL(KIND=dp) :: exp_tot, R_rad_sq, R_sq 589 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_to_c 590 REAL(KIND=dp), DIMENSION(3) :: R_l, R_r, Rp, Rx, Ry, Rz, s_shift 591 REAL(KIND=dp), DIMENSION(3, 0:l_max) :: R_pow_l 592 REAL(KIND=dp), DIMENSION(ncoset(l_max)) :: S_R_C 593 594 S_R(:) = 0.0_dp 595 S_R_C(:) = 0.0_dp 596 597 s_shift = MATMUL(h_inv, -R) 598 R_l = -R_c + s_shift 599 R_r = R_c + s_shift 600 601 R_rad_sq = R_rad**2 602 603 DO sx = CEILING(R_l(1)), FLOOR(R_r(1)) 604 Rx = hmat(:, 1)*sx 605 DO sy = CEILING(R_l(2)), FLOOR(R_r(2)) 606 Ry = hmat(:, 2)*sy 607 DO sz = CEILING(R_l(3)), FLOOR(R_r(3)) 608 Rz = hmat(:, 3)*sz 609 Rp = Rx + Ry + Rz 610 R_sq = (Rp(1) + R(1))**2 + (Rp(2) + R(2))**2 + (Rp(3) + R(3))**2 611 IF (R_sq .GT. R_rad_sq) CYCLE 612 exp_tot = EXP(-alpha*R_sq) ! cost: exp_w flops 613 DO k = 1, 3 614 R_pow_l(k, 0) = 1.0_dp 615 DO l = 1, l_max 616 R_pow_l(k, l) = R_pow_l(k, l - 1)*(Rp(k) + R(k)) 617 ENDDO 618 ENDDO 619 DO lco = 1, ncoset(l_max) 620 CALL get_l(lco, l, lx, ly, lz) 621 S_R_C(coset(lx, ly, lz)) = S_R_C(coset(lx, ly, lz)) + R_pow_l(1, lx)*R_pow_l(2, ly)*R_pow_l(3, lz)*exp_tot ! cost: 4 flops 622 ENDDO 623 ENDDO 624 ENDDO 625 ENDDO 626 627 CALL create_hermite_to_cartesian(alpha, l_max, h_to_c) 628 629 DO lco = 1, ncoset(l_max) 630 CALL get_l(lco, l, lx, ly, lz) 631 632 DO llx = 0, lx 633 DO lly = 0, ly 634 DO llz = 0, lz 635 S_R(coset(lx, ly, lz)) = S_R(coset(lx, ly, lz)) + & ! cost: 4 flops 636 h_to_c(llx, lx)*h_to_c(lly, ly)*h_to_c(llz, lz)* & 637 S_R_C(coset(llx, lly, llz)) 638 ENDDO 639 ENDDO 640 ENDDO 641 ENDDO 642 S_R(:) = S_R(:)*(alpha/pi)**1.5_dp 643 644 END SUBROUTINE pgf_sum_2c_rspace_3d 645 646! ************************************************************************************************** 647!> \brief Compute 1d sum 648!> S_G(l, alpha) = inv_lgth*sum_G( C(l, alpha, delta_l, G) ) with 649!> C(l, alpha, delta_l, G) = prefactor*|G|^(l-delta_l) exp(-alpha*G^2) 650!> if G not equal 0 651!> C(l = 0, alpha, delta_l, 0) = 1, C(l>0, alpha, delta_l, 0) = 0 652!> dG = inv_lgth*twopi and G = -G_bound*dG, (-G_bound + 1)*dG, ..., G_bound*dG 653!> for all l < = l_max. 654!> \param S_G ... 655!> \param alpha ... 656!> \param inv_lgth ... 657!> \param G_min ... 658!> \param G_c ... 659!> \param delta_l ... 660!> \param prefactor ... 661!> \note needed for cutoff error estimate 662! ************************************************************************************************** 663 PURE SUBROUTINE pgf_sum_2c_gspace_1d_deltal(S_G, alpha, inv_lgth, G_min, G_c, delta_l, prefactor) 664 REAL(KIND=dp), DIMENSION(0:), INTENT(OUT) :: S_G 665 REAL(KIND=dp), INTENT(IN) :: alpha, inv_lgth 666 INTEGER, INTENT(IN) :: G_min, G_c 667 REAL(KIND=dp), INTENT(IN) :: delta_l, prefactor 668 669 CHARACTER(LEN=*), PARAMETER :: routineN = 'pgf_sum_2c_gspace_1d_deltal', & 670 routineP = moduleN//':'//routineN 671 672 INTEGER :: k, k0, k1, l, l_max 673 REAL(KIND=dp) :: dG, exp_tot, G, prefac 674 675 prefac = prefactor*inv_lgth 676 dG = inv_lgth*twopi ! positive 677 l_max = UBOUND(S_G, 1) 678 679 S_G(:) = 0.0_dp 680 IF (0 .LT. G_min) THEN 681 k0 = G_min; k1 = 0 682 ELSE IF (G_c .LT. 0) THEN 683 k0 = 0; k1 = G_c 684 ELSE ! Gmin <= 0 /\ 0 <= Gc 685 S_G(0) = prefac 686 k0 = 1; k1 = -1 687 ENDIF 688 ! positive range 689 IF (0 .LT. k0) THEN 690 DO k = k0, G_c 691 G = k*dG; exp_tot = EXP(-alpha*G**2)*prefac 692 DO l = 0, l_max 693 S_G(l) = S_G(l) + G**(l - delta_l)*exp_tot 694 ENDDO 695 ENDDO 696 ENDIF 697 ! negative range 698 IF (k1 .LT. 0) THEN 699 DO k = G_min, k1 700 G = k*dG; exp_tot = EXP(-alpha*G**2)*prefac 701 DO l = 0, l_max 702 S_G(l) = S_G(l) + ABS(G)**(l - delta_l)*exp_tot 703 ENDDO 704 ENDDO 705 ENDIF 706 END SUBROUTINE pgf_sum_2c_gspace_1d_deltal 707 708! ************************************************************************************************** 709!> \brief Compute Ewald-like sum for 3-center integrals in 1 dimension 710!> S_G(l, m, n, alpha, beta, gamma) = i^(l+m+n)*(-1)^(l+m)*inv_lgth^2* 711!> sum_G sum_G'( exp(i G R1) 712!> C(l,alpha,G) C(m,beta,G'-G) C(n,gamma,G') exp(i G' R2) ) 713!> for all l < = l_max, m <= m_max, n <= n_max. 714!> a_mm is the minimax exponent. 715!> alpha = 1/(4 zeta), beta = 1/(4 zetb), gamma = 1/(4 zetc) + a_mm 716!> R1 = RB-RA; R2 = RC-RB 717!> Note on method / order arguments: 718!> Three equivalent methods (Poisson summation) to compute this sum over 719!> Cartesian Gaussians C or Hermite Gaussians H and 720!> reciprocal lattice vectors G or direct lattice vectors R: 721!> - method 1: sum_G sum_G' C(G) C(G,G') C(G') 722!> - method 2: sum_G sum_R C(G) C(R) 723!> - method 3: sum_R sum_R' H(R, R') 724!> The order parameter selects the Gaussian functions over which the sum is performed 725!> method 1: order = 1, 2, 3 726!> method 2: order = 1, 2, 3 727!> method 3: order = 1 728!> If method and order are not present, the method / order that converges fastest is 729!> automatically chosen. 730!> \param S_G ... 731!> \param RA ... 732!> \param RB ... 733!> \param RC ... 734!> \param zeta ... 735!> \param zetb ... 736!> \param zetc ... 737!> \param a_mm ... 738!> \param lgth ... 739!> \param G_bounds_1 ... 740!> \param R_bounds_2 ... 741!> \param R_bounds_3 ... 742!> \param method ... 743!> \param method_out ... 744!> \param order ... 745! ************************************************************************************************** 746 SUBROUTINE pgf_sum_3c_1d(S_G, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_bounds_3) 747 REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(OUT) :: S_G 748 REAL(KIND=dp), INTENT(IN) :: RA, RB, RC, zeta, zetb, zetc, a_mm, lgth 749 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_bounds_3 750 751 CHARACTER(LEN=*), PARAMETER :: routineN = 'pgf_sum_3c_1d', routineP = moduleN//':'//routineN 752 753 INTEGER :: l_max, m_max, n_max 754 INTEGER :: nR1, nR2 755 INTEGER :: prop_exp 756 757 l_max = UBOUND(S_G, 1) 758 m_max = UBOUND(S_G, 2) 759 n_max = UBOUND(S_G, 3) 760 761 nR1 = 2*FLOOR(R_bounds_3(1)) + 1 762 nR2 = 2*FLOOR(R_bounds_3(2)) + 1 763 764 IF (nR1*nR2 > 1 + nR1*2) THEN 765 prop_exp = 1 766 ELSE 767 prop_exp = 0 768 ENDIF 769 770 IF (MAXVAL([l_max, m_max, n_max]) .GT. ${lmax_unroll}$) THEN 771 CALL pgf_sum_3c_rspace_1d_generic(S_G, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_bounds_3) 772 ELSE 773#:for l_max in range(0,lmax_unroll+1) 774 IF (l_max == ${l_max}$) THEN 775#:for m_max in range(0,lmax_unroll+1) 776 IF (m_max == ${m_max}$) THEN 777#:for n_max in range(0, lmax_unroll+1) 778 IF (n_max == ${n_max}$) THEN 779#:for prop_exp in range(0,2) 780 IF (prop_exp == ${prop_exp}$) THEN 781 CALL pgf_sum_3c_rspace_1d_${l_max}$_${m_max}$_${n_max}$_exp_${prop_exp}$ (S_G, RA, RB, RC, & 782 zeta, zetb, zetc, a_mm, lgth, R_bounds_3) 783 RETURN 784 ENDIF 785#:endfor 786 ENDIF 787#:endfor 788 ENDIF 789#:endfor 790 ENDIF 791#:endfor 792 ENDIF 793 794 END SUBROUTINE pgf_sum_3c_1d 795 796! ************************************************************************************************** 797!> \brief Roughly estimated number of floating point operations 798!> \param ns_G1 ... 799!> \param ns_G2 ... 800!> \param l ... 801!> \param m ... 802!> \param n ... 803!> \return ... 804! ************************************************************************************************** 805 PURE FUNCTION nsum_3c_gspace_1d() 806 INTEGER(KIND=int_8) :: nsum_3c_gspace_1d 807 808 nsum_3c_gspace_1d = 15 809 END FUNCTION 810 811! ************************************************************************************************** 812!> \brief Roughly estimated number of floating point operations 813!> \param ns_G ... 814!> \param ns_R ... 815!> \param l ... 816!> \param m ... 817!> \param n ... 818!> \return ... 819! ************************************************************************************************** 820 PURE FUNCTION nsum_product_3c_gspace_1d(ns_G, ns_R) 821 REAL(KIND=dp), INTENT(IN) :: ns_G, ns_R 822 INTEGER(KIND=int_8) :: nsum_product_3c_gspace_1d 823 824 nsum_product_3c_gspace_1d = MIN(19, NINT(ns_G*(3 + ns_R*2))) 825 END FUNCTION 826 827! ************************************************************************************************** 828!> \brief Roughly estimated number of floating point operations 829!> \param ns_R1 ... 830!> \param ns_R2 ... 831!> \param l ... 832!> \param m ... 833!> \param n ... 834!> \return ... 835! ************************************************************************************************** 836 PURE FUNCTION nsum_3c_rspace_1d(ns_R1, ns_R2) 837 REAL(KIND=dp), INTENT(IN) :: ns_R1, ns_R2 838 INTEGER(KIND=int_8) :: nsum_3c_rspace_1d 839 840 nsum_3c_rspace_1d = NINT(MIN((4 + ns_R1*2), ns_R1*(ns_R2 + 1)), KIND=int_8) 841 END FUNCTION 842 843! ************************************************************************************************** 844!> \brief Helper routine: compute SQRT(alpha/pi) (-1)^n sum_(R, R') sum_{t=0}^{l+m} E(t,l,m) H(RC - P(R) - R', t + n, alpha) 845!> with alpha = 1.0_dp/((a + b + c)/((a + b)*c) + 4.0_dp*a_mm), 846!> P(R) = (a*(RA + R) + b*RB)/(a + b) 847!> \param S_R ... 848!> \param RA ... 849!> \param RB ... 850!> \param RC ... 851!> \param zeta ... 852!> \param zetb ... 853!> \param zetc ... 854!> \param a_mm ... 855!> \param lgth ... 856!> \param R_c ... 857! ************************************************************************************************** 858 PURE SUBROUTINE pgf_sum_3c_rspace_1d_generic(S_R, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_c) 859 REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(OUT) :: S_R 860 REAL(KIND=dp), INTENT(IN) :: RA, RB, RC, zeta, zetb, zetc, a_mm, lgth 861 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_c 862 863 INTEGER :: ll, mm, l, k, l_max, m, m_max, n, n_max, rr1, rr2, t, l_tot_max 864 REAL(KIND=dp) :: alpha, dR, exp_tot, R, R1, R2, R_offset, & 865 R_pow_t, R_tmp, rr1_delta, rr2_delta, c1, c2, c3 866 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: S_R_t 867 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_to_c 868 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: E 869 870 dR = lgth 871 alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm) 872 l_max = UBOUND(S_R, 1) 873 m_max = UBOUND(S_R, 2) 874 n_max = UBOUND(S_R, 3) 875 l_tot_max = l_max + m_max + n_max 876 877 ALLOCATE (S_R_t(0:l_max + m_max + n_max)) 878 ALLOCATE (E(-1:l_max + m_max + 1, -1:l_max, -1:m_max)) 879 880 S_R(:, :, :) = 0.0_dp 881 882 R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb) 883 884 ! inline CALL create_hermite_to_cartesian(alpha, l_tot_max, h_to_c) 885 ALLOCATE (h_to_c(-1:l_tot_max + 1, 0:l_tot_max)) 886 h_to_c(:, :) = 0.0_dp 887 h_to_c(0, 0) = 1.0_dp 888 DO l = 0, l_tot_max - 1 889 DO k = 0, l + 1 890 h_to_c(k, l + 1) = -(k + 1)*h_to_c(k + 1, l) + 2.0_dp*alpha*h_to_c(k - 1, l) 891 ENDDO 892 ENDDO 893 894 rr1_delta = (RA - RB)/dR 895 DO rr1 = CEILING(-R_c(1) + rr1_delta), FLOOR(R_c(1) + rr1_delta) 896 S_R_t(:) = 0.0_dp 897 R1 = rr1*dR 898 R_tmp = R_offset + R1*zeta/(zeta + zetb) 899 rr2_delta = -R_tmp/dR 900 DO rr2 = CEILING(-R_c(2) + rr2_delta), FLOOR(R_c(2) + rr2_delta) 901 R2 = rr2*dR 902 R = R_tmp + R2 903 exp_tot = EXP(-alpha*R**2) ! cost: exp_w flops 904 R_pow_t = 1.0_dp 905 DO t = 0, l_max + m_max + n_max 906 S_R_t(t) = S_R_t(t) + R_pow_t*exp_tot ! cost: 2 flops 907 R_pow_t = R_pow_t*R ! cost: 1 flop 908 ENDDO 909 ENDDO 910 911 ! C --> H 912 S_R_t(:) = MATMUL(TRANSPOSE(h_to_c(0:l_max + m_max + n_max, 0:l_max + m_max + n_max)), S_R_t)*SQRT(alpha/pi) 913 914 ! H --> HH 915 !inline CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA-R1, RB, 2, E) 916 917 E(:, :, :) = 0.0_dp 918 E(0, 0, 0) = EXP(-zeta*zetb/(zeta + zetb)*(RA - R1 - RB)**2) 919 920 c1 = 1.0_dp/(zeta + zetb) 921 c2 = 2.0_dp*(zetb/(zeta + zetb))*(RB - (RA - R1)) 922 c3 = 2.0_dp*(zeta/(zeta + zetb))*(RA - R1 - RB) 923 924 DO mm = 0, m_max - 1 925 DO ll = 0, l_max - 1 926 DO t = 0, ll + mm + 1 927 E(t, ll + 1, mm) = zeta*(c1*E(t - 1, ll, mm) + & 928 c2*E(t, ll, mm) + & 929 2*(t + 1)*E(t + 1, ll, mm) - & 930 2*ll*E(t, ll - 1, mm)) 931 E(t, ll, mm + 1) = zetb*(c1*E(t - 1, ll, mm) + & 932 c3*E(t, ll, mm) + & 933 2*(t + 1)*E(t + 1, ll, mm) - & 934 2*mm*E(t, ll, mm - 1)) 935 ENDDO 936 ENDDO 937 ENDDO 938 939 DO ll = 0, l_max - 1 940 DO t = 0, ll + m_max + 1 941 E(t, ll + 1, m_max) = zeta*(c1*E(t - 1, ll, m_max) + & 942 c2*E(t, ll, m_max) + & 943 2*(t + 1)*E(t + 1, ll, m_max) - & 944 2*ll*E(t, ll - 1, m_max)) 945 ENDDO 946 ENDDO 947 948 DO mm = 0, m_max - 1 949 DO t = 0, l_max + mm + 1 950 E(t, l_max, mm + 1) = zetb*(c1*E(t - 1, l_max, mm) + & 951 c3*E(t, l_max, mm) + & 952 2*(t + 1)*E(t + 1, l_max, mm) - & 953 2*mm*E(t, l_max, mm - 1)) 954 ENDDO 955 ENDDO 956 957 DO n = 0, n_max 958 DO m = 0, m_max 959 DO l = 0, l_max 960 DO t = 0, l + m 961 S_R(l, m, n) = S_R(l, m, n) + E(t, l, m)*(-1)**n*S_R_t(t + n) ! cost: 5 flops 962 ENDDO 963 ENDDO 964 ENDDO 965 ENDDO 966 ENDDO 967 968 S_R = S_R*pi**(-0.5_dp)*((zeta + zetb)/(zeta*zetb))**(-0.5_dp) 969 END SUBROUTINE 970 971! ************************************************************************************************** 972!> \brief Helper routine: compute SQRT(alpha/pi) (-1)^n sum_(R, R') sum_{t=0}^{l+m} E(t,l,m) H(RC - P(R) - R', t + n, alpha) 973!> with alpha = 1.0_dp/((a + b + c)/((a + b)*c) + 4.0_dp*a_mm), 974!> P(R) = (a*(RA + R) + b*RB)/(a + b) 975!> \param S_R ... 976!> \param RA ... 977!> \param RB ... 978!> \param RC ... 979!> \param zeta ... 980!> \param zetb ... 981!> \param zetc ... 982!> \param a_mm ... 983!> \param lgth ... 984!> \param R_c ... 985! ************************************************************************************************** 986#:for prop_exp in range(0,2) 987#:for l_max in range(0, lmax_unroll+1) 988#:for m_max in range(0, lmax_unroll+1) 989#:for n_max in range(0, lmax_unroll+1) 990#:set l_tot_max = l_max + m_max + n_max 991 PURE SUBROUTINE pgf_sum_3c_rspace_1d_${l_max}$_${m_max}$_${n_max}$_exp_${prop_exp}$ ( & 992 S_R, RA, RB, RC, zeta, zetb, zetc, a_mm, lgth, R_c) 993 REAL(KIND=dp), DIMENSION(0:, 0:, 0:), INTENT(OUT) :: S_R 994 REAL(KIND=dp), INTENT(IN) :: RA, RB, RC, zeta, zetb, zetc, a_mm, lgth 995 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_c 996 INTEGER :: rr1, rr2, rr2_l, rr2_r, rr1_l, rr1_r 997 REAL(KIND=dp) :: alpha, alpha_E, dR, exp2_Rsq, R, R1, R_offset, & 998 R_pow_t, R_tmp, rr1_delta, rr2_delta 999 1000#:if l_tot_max > 0 1001 REAL(KIND=dp) :: c1, c2, c3 1002#:endif 1003 REAL(KIND=dp) :: ${", ".join(["S_R_t_{}".format(t) for t in range(0,l_tot_max+1)])}$ 1004 REAL(KIND=dp) :: ${", ".join(["S_R_t2_{}".format(t) for t in range(0,l_tot_max+1)])}$ 1005 REAL(KIND=dp) :: ${", ".join([", ".join(["h_to_c_{}_{}".format(l1,l2) for l1 in range(0,l2+1)]) for l2 in range(0,l_tot_max+1)])}$ 1006 REAL(KIND=dp) :: ${", ".join([", ".join([", ".join(["E_{}_{}_{}".format(t,l,m) for t in range(0,l+m+1)]) for l in range(0,l_max+1)]) for m in range(0,m_max+1)])}$ 1007 1008#:if prop_exp 1009 REAL(KIND=dp) :: exp2_2RdR, exp_dRsq, exp_2dRsq !exp_E_dRsq, exp_E_2dRsq, exp_E_2RdR, exp_E_Rsq 1010#:endif 1011 1012 dR = lgth 1013 alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm) 1014 1015 S_R(:, :, :) = 0.0_dp 1016 1017 R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb) 1018 1019 h_to_c_0_0 = SQRT(alpha/pi) 1020 1021#:for l in range(0, l_tot_max) 1022#:for k in range(0, l+2) 1023#:if k<l or k>0 1024 h_to_c_${k}$_${l+1}$ = #{if k<l}#-${k+1}$*h_to_c_${k+1}$_${l}$#{endif}# #{if k > 0}#+2*alpha*h_to_c_${k-1}$_${l}$#{endif}# 1025#:else 1026 h_to_c_${k}$_${l+1}$ = 0.0_dp 1027#:endif 1028#:endfor 1029#:endfor 1030 1031#:if prop_exp 1032 exp_dRsq = exp(-alpha*dR*dR) 1033 exp_2dRsq = exp_dRsq*exp_dRsq 1034#:endif 1035 1036 rr1_delta = (RA - RB)/dR 1037 1038 rr1_l = CEILING(-R_c(1) + rr1_delta) 1039 rr1_r = FLOOR(R_c(1) + rr1_delta) 1040 1041 R1 = rr1_l*dR 1042 1043 alpha_E = zeta*zetb/(zeta + zetb) 1044 1045 DO rr1 = rr1_l, rr1_r 1046#:for t in range(0, l_tot_max + 1) 1047 S_R_t_${t}$ = 0.0_dp 1048 S_R_t2_${t}$ = 0.0_dp 1049#:endfor 1050 R_tmp = R_offset + R1*zeta/(zeta + zetb) 1051 rr2_delta = -R_tmp/dR 1052 1053 rr2_l = CEILING(-R_c(2) + rr2_delta) 1054 rr2_r = FLOOR(R_c(2) + rr2_delta) 1055 1056 R = R_tmp + (rr2_l)*dR 1057 1058#:if prop_exp 1059 exp2_2RdR = exp(-2*alpha*R*dR) 1060 exp2_Rsq = exp(-alpha*R*R) 1061#:endif 1062 1063 DO rr2 = rr2_l, rr2_r 1064 R_pow_t = 1.0_dp 1065#:if not prop_exp 1066 exp2_Rsq = exp(-alpha*R*R) 1067#:endif 1068#:for t in range(0, l_tot_max + 1) 1069 S_R_t_${t}$ = S_R_t_${t}$+R_pow_t*exp2_Rsq 1070#:if t<l_tot_max 1071 R_pow_t = R_pow_t*R 1072#:endif 1073#:endfor 1074 1075#:if prop_exp 1076 exp2_Rsq = exp2_Rsq*exp_dRsq*exp2_2RdR 1077 exp2_2RdR = exp2_2RdR*exp_2dRsq 1078#:endif 1079 R = R + dR 1080 ENDDO 1081 1082 ! C --> H 1083#:for l in range(0, l_tot_max+1) 1084#:for k in range(0, l+1) 1085 S_R_t2_${l}$ = S_R_t2_${l}$+h_to_c_${k}$_${l}$*S_R_t_${k}$ 1086#:endfor 1087#:endfor 1088 1089 ! H --> HH 1090 E_0_0_0 = exp(-alpha_E*(RA - RB - R1)*(RA - RB - R1)) 1091 1092#:if l_tot_max > 0 1093 c1 = 1.0_dp/(zeta + zetb) 1094 c2 = 2.0_dp*(zetb/(zeta + zetb))*(RB - (RA - R1)) 1095 c3 = 2.0_dp*(zeta/(zeta + zetb))*(RA - R1 - RB) 1096#:endif 1097 1098#:for m in range(0,m_max+1) 1099#:for l in range(0,l_max+1) 1100#:for t in range(0,l+m+2) 1101#:if l < l_max 1102 E_${t}$_${l+1}$_${m}$ = zeta*(#{if t>0}# c1*E_${t-1}$_${l}$_${m}$#{endif}# & 1103 #{if t<=l+m}# +c2*E_${t}$_${l}$_${m}$&#{endif}# 1104 #{if t<l+m}# +${2*(t+1)}$*E_${t+1}$_${l}$_${m}$ &#{endif}# 1105 #{if l>0 and t<=l-1+m}#-${2*l}$*E_${t}$_${l-1}$_${m}$#{endif}#) 1106#:endif 1107#:if m < m_max 1108 E_${t}$_${l}$_${m+1}$ = zetb*(#{if t>0}# c1*E_${t-1}$_${l}$_${m}$#{endif}# & 1109 #{if t<=l+m}#+c3*E_${t}$_${l}$_${m}$&#{endif}# 1110 #{if t<l+m}# +${2*(t+1)}$*E_${t+1}$_${l}$_${m}$ &#{endif}# 1111 #{if m>0 and t<=m-1+l}#-${2*m}$*E_${t}$_${l}$_${m-1}$#{endif}#) 1112#:endif 1113#:endfor 1114#:endfor 1115#:endfor 1116 1117#:for n in range(0, n_max+1) 1118#:for m in range(0, m_max+1) 1119#:for l in range(0, l_max+1) 1120#:for t in range(0, l+m+1) 1121 S_R(${l}$, ${m}$, ${n}$) = S_R(${l}$, ${m}$, ${n}$) + E_${t}$_${l}$_${m}$*(${(-1)**n}$)*S_R_t2_${t+n}$ ! cost: 5 flops 1122#:endfor 1123#:endfor 1124#:endfor 1125#:endfor 1126 R1 = R1 + dR 1127 ENDDO 1128 1129 S_R = S_R*pi**(-0.5_dp)*((zeta + zetb)/(zeta*zetb))**(-0.5_dp) 1130 END SUBROUTINE 1131#:endfor 1132#:endfor 1133#:endfor 1134#:endfor 1135 1136! ************************************************************************************************** 1137!> \brief As pgf_sum_3c_1d but 3d sum required for non-orthorhombic cells 1138!> \param S_G ... 1139!> \param la_max ... 1140!> \param lb_max ... 1141!> \param lc_max ... 1142!> \param RA ... 1143!> \param RB ... 1144!> \param RC ... 1145!> \param zeta ... 1146!> \param zetb ... 1147!> \param zetc ... 1148!> \param a_mm ... 1149!> \param hmat ... 1150!> \param h_inv ... 1151!> \param vol ... 1152!> \param G_bounds_1 ... 1153!> \param R_bounds_2 ... 1154!> \param R_bounds_3 ... 1155!> \param G_rads_1 ... 1156!> \param R_rads_2 ... 1157!> \param R_rads_3 ... 1158!> \param method ... 1159!> \param method_out ... 1160!> \param order ... 1161! ************************************************************************************************** 1162 SUBROUTINE pgf_sum_3c_3d(S_G, la_max, lb_max, lc_max, RA, RB, RC, & 1163 zeta, zetb, zetc, a_mm, hmat, h_inv, vol, & 1164 G_bounds_1, R_bounds_2, R_bounds_3, & 1165 G_rads_1, R_rads_2, R_rads_3, & 1166 method, method_out, order) 1167 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_G 1168 INTEGER, INTENT(IN) :: la_max, lb_max, lc_max 1169 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: RA, RB, RC 1170 REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc, a_mm 1171 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv 1172 REAL(KIND=dp), INTENT(IN) :: vol 1173 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: G_bounds_1, R_bounds_2 1174 REAL(KIND=dp), DIMENSION(2, 3), INTENT(IN) :: R_bounds_3 1175 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G_rads_1, R_rads_2 1176 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_rads_3 1177 INTEGER, INTENT(IN) :: method 1178 INTEGER, INTENT(OUT), OPTIONAL :: method_out 1179 INTEGER, INTENT(IN), OPTIONAL :: order 1180 1181 INTEGER :: l_max, m_max, n_max, sum_method, & 1182 sum_order 1183 LOGICAL :: assert_stm 1184 REAL(KIND=dp) :: alpha, beta, G_rad, gamma, R_rad 1185 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: S_G_tmp 1186 REAL(KIND=dp), DIMENSION(3) :: G_bound, R1, R2, R_bound 1187 REAL(KIND=dp), DIMENSION(3, 3) :: ht 1188 1189 IF (PRESENT(order)) THEN 1190 sum_order = order 1191 ELSE 1192 sum_order = 0 1193 ENDIF 1194 1195 sum_method = method 1196 1197 alpha = 0.25_dp/zeta 1198 beta = 0.25_dp/zetb 1199 gamma = 0.25_dp/zetc + a_mm 1200 1201 l_max = la_max 1202 m_max = lb_max 1203 n_max = lc_max 1204 1205 R1 = RB - RA 1206 R2 = RC - RB 1207 1208 ht = twopi*TRANSPOSE(h_inv) 1209 1210 SELECT CASE (sum_method) 1211 CASE (1) ! sum_G sum_G' C(G) C(G,G') C(G') 1212 1213 IF (sum_order .EQ. 0) THEN 1214 sum_order = MAXLOC(G_bounds_1(:, 1), DIM=1) 1215 assert_stm = MINLOC(G_bounds_1(:, 1), DIM=1) == & 1216 MINLOC(G_bounds_1(:, 2), DIM=1) .AND. & 1217 MINLOC(G_bounds_1(:, 1), DIM=1) == & 1218 MINLOC(G_bounds_1(:, 3), DIM=1) 1219 CPASSERT(assert_stm) 1220 ENDIF 1221 1222 CALL pgf_sum_3c_gspace_3d(S_G, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, ht, vol, G_bounds_1, G_rads_1, sum_order) 1223 1224 CASE (2) ! sum_G sum_R C(G) C(R) 1225 IF (sum_order .EQ. 0) THEN 1226 sum_order = MINLOC(G_bounds_1(:, 1), DIM=1) 1227 assert_stm = MINLOC(G_bounds_1(:, 1), DIM=1) == & 1228 MINLOC(G_bounds_1(:, 2), DIM=1) .AND. & 1229 MINLOC(G_bounds_1(:, 1), DIM=1) == & 1230 MINLOC(G_bounds_1(:, 3), DIM=1) 1231 CPASSERT(assert_stm) 1232 ENDIF 1233 R_rad = R_rads_2(sum_order) 1234 G_rad = G_rads_1(sum_order) 1235 R_bound = R_bounds_2(sum_order, :) 1236 G_bound = G_bounds_1(sum_order, :) 1237 SELECT CASE (sum_order) 1238 CASE (1) 1239 ALLOCATE (S_G_tmp(ncoset(l_max), ncoset(m_max), ncoset(n_max))) 1240 CALL pgf_sum_product_3c_gspace_3d(S_G_tmp, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, hmat, h_inv, vol, & 1241 R_bound, G_bound, R_rad, G_rad) 1242 S_G = RESHAPE(S_G_tmp, SHAPE(S_G), order=[1, 2, 3]) 1243 CASE (2) 1244 ALLOCATE (S_G_tmp(ncoset(m_max), ncoset(l_max), ncoset(n_max))) 1245 CALL pgf_sum_product_3c_gspace_3d(S_G_tmp, m_max, l_max, n_max, -R1, R1 + R2, beta, alpha, gamma, hmat, h_inv, vol, & 1246 R_bound, G_bound, R_rad, G_rad) 1247 S_G = RESHAPE(S_G_tmp, SHAPE(S_G), order=[2, 1, 3]) 1248 CASE (3) 1249 ALLOCATE (S_G_tmp(ncoset(n_max), ncoset(m_max), ncoset(l_max))) 1250 CALL pgf_sum_product_3c_gspace_3d(S_G_tmp, n_max, m_max, l_max, -R2, -R1, gamma, beta, alpha, hmat, h_inv, vol, & 1251 R_bound, G_bound, R_rad, G_rad) 1252 S_G = RESHAPE(S_G_tmp, SHAPE(S_G), order=[3, 2, 1]) 1253 END SELECT 1254 CASE (3) ! sum_R sum_R' H(R, R') 1255 CALL pgf_sum_3c_rspace_3d(S_G, l_max, m_max, n_max, RA, RB, RC, zeta, zetb, zetc, a_mm, hmat, h_inv, & 1256 R_bounds_3, R_rads_3) 1257 S_G = S_G*pi**(-1.5_dp)*((zeta + zetb)/(zeta*zetb))**(-1.5_dp) 1258 END SELECT 1259 1260 IF (PRESENT(method_out)) method_out = sum_method 1261 1262 END SUBROUTINE pgf_sum_3c_3d 1263 1264! ************************************************************************************************** 1265!> \brief Roughly estimated number of floating point operations 1266!> \param ns_G1 ... 1267!> \param ns_G2 ... 1268!> \param l ... 1269!> \param m ... 1270!> \param n ... 1271!> \return ... 1272! ************************************************************************************************** 1273 PURE FUNCTION nsum_3c_gspace_3d(ns_G1, ns_G2, l, m, n) 1274 REAL(KIND=dp), INTENT(IN) :: ns_G1, ns_G2 1275 INTEGER, INTENT(IN) :: l, m, n 1276 INTEGER(KIND=int_8) :: nsum_3c_gspace_3d 1277 1278 nsum_3c_gspace_3d = NINT(ns_G1*ns_G2*(5*exp_w + ncoset(l)*ncoset(m)*ncoset(n)*4), KIND=int_8) 1279 1280 END FUNCTION 1281 1282! ************************************************************************************************** 1283!> \brief ... 1284!> \param S_G ... 1285!> \param l_max ... 1286!> \param m_max ... 1287!> \param n_max ... 1288!> \param R1 ... 1289!> \param R2 ... 1290!> \param alpha ... 1291!> \param beta ... 1292!> \param gamma ... 1293!> \param ht ... 1294!> \param vol ... 1295!> \param G_c ... 1296!> \param G_rad ... 1297!> \param sum_order ... 1298!> \param coulomb ... 1299! ************************************************************************************************** 1300 PURE SUBROUTINE pgf_sum_3c_gspace_3d(S_G, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, & 1301 ht, vol, G_c, G_rad, sum_order, coulomb) 1302 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_G 1303 INTEGER, INTENT(IN) :: l_max, m_max, n_max 1304 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R1, R2 1305 REAL(KIND=dp), INTENT(IN) :: alpha, beta, gamma 1306 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: ht 1307 REAL(KIND=dp), INTENT(IN) :: vol 1308 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: G_c 1309 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G_rad 1310 INTEGER, INTENT(IN) :: sum_order 1311 LOGICAL, INTENT(IN), OPTIONAL :: coulomb 1312 1313 INTEGER, DIMENSION(3) :: G1c, G2c, G3c 1314 INTEGER :: g1x, g1y, g1z, g2x, g2y, g2z, g3x, g3y, & 1315 g3z 1316 COMPLEX(KIND=dp), DIMENSION(ncoset(l_max), ncoset( & 1317 m_max), ncoset(n_max)) :: S_G_c 1318 LOGICAL :: use_coulomb 1319 REAL(KIND=dp) :: G1_sq, G2_sq, G3_sq 1320 REAL(KIND=dp), DIMENSION(3) :: G1, G1_x, G1_y, G1_z, G2, G2_x, G2_y, & 1321 G2_z, G3, G3_x, G3_y, G3_z, G_rads_sq 1322 1323 S_G_c(:, :, :) = 0.0_dp 1324 1325 G1c = FLOOR(G_c(1, :)) 1326 G2c = FLOOR(G_c(2, :)) 1327 G3c = FLOOR(G_c(3, :)) 1328 1329 ! we can not precompute exponentials for general cell 1330 ! Perform double G sum 1331 G_rads_sq = G_rad**2 1332 1333 IF (PRESENT(coulomb)) THEN 1334 use_coulomb = coulomb 1335 ELSE 1336 use_coulomb = .FALSE. 1337 ENDIF 1338 1339 SELECT CASE (sum_order) 1340 CASE (1) 1341 DO g2x = -G2c(1), G2c(1) 1342 G2_x = ht(:, 1)*g2x 1343 DO g2y = -G2c(2), G2c(2) 1344 G2_y = ht(:, 2)*g2y 1345 DO g2z = -G2c(3), G2c(3) 1346 G2_z = ht(:, 3)*g2z 1347 G2 = G2_x + G2_y + G2_z 1348 G2_sq = G2(1)**2 + G2(2)**2 + G2(3)**2 1349 IF (G2_sq .GT. G_rads_sq(2)) CYCLE 1350 DO g3x = -G3c(1), G3c(1) 1351 G3_x = ht(:, 1)*g3x 1352 DO g3y = -G3c(2), G3c(2) 1353 G3_y = ht(:, 2)*g3y 1354 DO g3z = -G3c(3), G3c(3) 1355 G3_z = ht(:, 3)*g3z 1356 G3 = G3_x + G3_y + G3_z 1357 G3_sq = G3(1)**2 + G3(2)**2 + G3(3)**2 1358 IF (G3_sq .GT. G_rads_sq(3)) CYCLE 1359 G1 = G3 - G2 1360 G1_sq = G1(1)**2 + G1(2)**2 + G1(3)**2 1361 IF (G1_sq .GT. G_rads_sq(1)) CYCLE 1362 IF (use_coulomb) THEN 1363 IF (g3x == 0 .AND. g3y == 0 .AND. g3z == 0) CYCLE 1364 ENDIF 1365 CALL pgf_product_3c_gspace_3d(S_G_c, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, & 1366 alpha, beta, gamma, R1, R2, use_coulomb) 1367 ENDDO 1368 ENDDO 1369 ENDDO 1370 ENDDO 1371 ENDDO 1372 ENDDO 1373 CASE (2) 1374 DO g1x = -G1c(1), G1c(1) 1375 G1_x = ht(:, 1)*g1x 1376 DO g1y = -G1c(2), G1c(2) 1377 G1_y = ht(:, 2)*g1y 1378 DO g1z = -G1c(3), G1c(3) 1379 G1_z = ht(:, 3)*g1z 1380 G1 = G1_x + G1_y + G1_z 1381 G1_sq = G1(1)**2 + G1(2)**2 + G1(3)**2 1382 IF (G1_sq .GT. G_rads_sq(1)) CYCLE 1383 DO g3x = -G3c(1), G3c(1) 1384 G3_x = ht(:, 1)*g3x 1385 DO g3y = -G3c(2), G3c(2) 1386 G3_y = ht(:, 2)*g3y 1387 DO g3z = -G3c(3), G3c(3) 1388 G3_z = ht(:, 3)*g3z 1389 G3 = G3_x + G3_y + G3_z 1390 G3_sq = G3(1)**2 + G3(2)**2 + G3(3)**2 1391 IF (G3_sq .GT. G_rads_sq(3)) CYCLE 1392 G2 = G3 - G1 1393 G2_sq = G2(1)**2 + G2(2)**2 + G2(3)**2 1394 IF (G2_sq .GT. G_rads_sq(2)) CYCLE 1395 IF (use_coulomb) THEN 1396 IF (g3x == 0 .AND. g3y == 0 .AND. g3z == 0) CYCLE 1397 ENDIF 1398 CALL pgf_product_3c_gspace_3d(S_G_c, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, & 1399 alpha, beta, gamma, R1, R2, use_coulomb) 1400 ENDDO 1401 ENDDO 1402 ENDDO 1403 ENDDO 1404 ENDDO 1405 ENDDO 1406 CASE (3) 1407 DO g1x = -G1c(1), G1c(1) 1408 G1_x = ht(:, 1)*g1x 1409 DO g1y = -G1c(2), G1c(2) 1410 G1_y = ht(:, 2)*g1y 1411 DO g1z = -G1c(3), G1c(3) 1412 G1_z = ht(:, 3)*g1z 1413 G1 = G1_x + G1_y + G1_z 1414 G1_sq = G1(1)**2 + G1(2)**2 + G1(3)**2 1415 IF (G1_sq .GT. G_rads_sq(1)) CYCLE 1416 DO g2x = -G2c(1), G2c(1) 1417 G2_x = ht(:, 1)*g2x 1418 DO g2y = -G2c(2), G2c(2) 1419 G2_y = ht(:, 2)*g2y 1420 DO g2z = -G2c(3), G2c(3) 1421 G2_z = ht(:, 3)*g2z 1422 G2 = G2_x + G2_y + G2_z 1423 G2_sq = G2(1)**2 + G2(2)**2 + G2(3)**2 1424 IF (G2_sq .GT. G_rads_sq(2)) CYCLE 1425 G3 = G1 + G2 1426 G3_sq = G3(1)**2 + G3(2)**2 + G3(3)**2 1427 IF (G3_sq .GT. G_rads_sq(3)) CYCLE 1428 IF (use_coulomb) THEN 1429 IF (g1x + g2x == 0 .AND. g1y + g2y == 0 .AND. g1z + g2z == 0) CYCLE 1430 ENDIF 1431 CALL pgf_product_3c_gspace_3d(S_G_c, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, & 1432 alpha, beta, gamma, R1, R2, use_coulomb) 1433 ENDDO 1434 ENDDO 1435 ENDDO 1436 ENDDO 1437 ENDDO 1438 ENDDO 1439 END SELECT 1440 S_G = REAL(S_G_c, KIND=dp)/vol**2 1441 1442 END SUBROUTINE 1443 1444! ************************************************************************************************** 1445!> \brief ... 1446!> \param S_G ... 1447!> \param G1 ... 1448!> \param G1_sq ... 1449!> \param G2 ... 1450!> \param G2_sq ... 1451!> \param G3 ... 1452!> \param G3_sq ... 1453!> \param l_max ... 1454!> \param m_max ... 1455!> \param n_max ... 1456!> \param alpha ... 1457!> \param beta ... 1458!> \param gamma ... 1459!> \param R1 ... 1460!> \param R2 ... 1461!> \param coulomb ... 1462! ************************************************************************************************** 1463 PURE SUBROUTINE pgf_product_3c_gspace_3d(S_G, G1, G1_sq, G2, G2_sq, G3, G3_sq, l_max, m_max, n_max, & 1464 alpha, beta, gamma, R1, R2, coulomb) 1465 1466 COMPLEX(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_G 1467 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G1 1468 REAL(KIND=dp), INTENT(IN) :: G1_sq 1469 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G2 1470 REAL(KIND=dp), INTENT(IN) :: G2_sq 1471 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G3 1472 REAL(KIND=dp), INTENT(IN) :: G3_sq 1473 INTEGER, INTENT(IN) :: l_max, m_max, n_max 1474 REAL(KIND=dp), INTENT(IN) :: alpha, beta, gamma 1475 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R1, R2 1476 LOGICAL, INTENT(IN) :: coulomb 1477 1478 COMPLEX(KIND=dp) :: exp_tot 1479 INTEGER :: k, l, lco, lx, ly, lz, m, mco, mx, my, & 1480 mz, n, nco, nx, ny, nz 1481 COMPLEX(KIND=dp), DIMENSION(ncoset(n_max)) :: S_G_n 1482 COMPLEX(KIND=dp), DIMENSION(ncoset(m_max)) :: S_G_m 1483 COMPLEX(KIND=dp), DIMENSION(ncoset(l_max)) :: S_G_l 1484 REAL(KIND=dp), DIMENSION(3, 0:l_max) :: G1_pow_l 1485 REAL(KIND=dp), DIMENSION(3, 0:m_max) :: G2_pow_m 1486 REAL(KIND=dp), DIMENSION(3, 0:n_max) :: G3_pow_n 1487 1488 exp_tot = EXP(gaussi*DOT_PRODUCT(G1, R1))*EXP(-alpha*G1_sq)* & ! cost: 5 exp_w flops 1489 EXP(-beta*G2_sq)* & 1490 EXP(-gamma*G3_sq)*EXP(gaussi*DOT_PRODUCT(G3, R2)) 1491 1492 IF (coulomb) exp_tot = exp_tot/G3_sq 1493 1494 DO k = 1, 3 1495 G1_pow_l(k, 0) = 1.0_dp 1496 DO l = 1, l_max 1497 G1_pow_l(k, l) = G1_pow_l(k, l - 1)*G1(k) 1498 ENDDO 1499 G2_pow_m(k, 0) = 1.0_dp 1500 DO m = 1, m_max 1501 G2_pow_m(k, m) = G2_pow_m(k, m - 1)*G2(k) 1502 ENDDO 1503 G3_pow_n(k, 0) = 1.0_dp 1504 DO n = 1, n_max 1505 G3_pow_n(k, n) = G3_pow_n(k, n - 1)*G3(k) 1506 ENDDO 1507 ENDDO 1508 1509 DO lco = 1, ncoset(l_max) 1510 CALL get_l(lco, l, lx, ly, lz) 1511 S_G_l(lco) = G1_pow_l(1, lx)*G1_pow_l(2, ly)*G1_pow_l(3, lz)*(-1.0_dp)**l*i_pow(l) 1512 ENDDO 1513 1514 DO mco = 1, ncoset(m_max) 1515 CALL get_l(mco, m, mx, my, mz) 1516 S_G_m(mco) = G2_pow_m(1, mx)*G2_pow_m(2, my)*G2_pow_m(3, mz)*(-1.0_dp)**m*i_pow(m) 1517 ENDDO 1518 1519 DO nco = 1, ncoset(n_max) 1520 CALL get_l(nco, n, nx, ny, nz) 1521 S_G_n(nco) = G3_pow_n(1, nx)*G3_pow_n(2, ny)*G3_pow_n(3, nz)*i_pow(n) 1522 ENDDO 1523 1524 DO nco = 1, ncoset(n_max) 1525 DO mco = 1, ncoset(m_max) 1526 DO lco = 1, ncoset(l_max) 1527 S_G(lco, mco, nco) = S_G(lco, mco, nco) + & ! cost: 4 flops 1528 S_G_l(lco)*S_G_m(mco)*S_G_n(nco)* & 1529 exp_tot 1530 ENDDO 1531 ENDDO 1532 ENDDO 1533 1534 END SUBROUTINE pgf_product_3c_gspace_3d 1535 1536! ************************************************************************************************** 1537!> \brief Roughly estimated number of floating point operations 1538!> \param ns_G ... 1539!> \param ns_R ... 1540!> \param l ... 1541!> \param m ... 1542!> \param n ... 1543!> \return ... 1544! ************************************************************************************************** 1545 PURE FUNCTION nsum_product_3c_gspace_3d(ns_G, ns_R, l, m, n) 1546 REAL(KIND=dp), INTENT(IN) :: ns_G, ns_R 1547 INTEGER, INTENT(IN) :: l, m, n 1548 INTEGER(KIND=int_8) :: nsum_product_3c_gspace_3d 1549 1550 nsum_product_3c_gspace_3d = & 1551 NINT( & 1552 ns_G*( & 1553 (exp_w*2) + & 1554 ns_R*(exp_w*2 + ncoset(l + m)*7) + & 1555 3*nsum_gaussian_overlap(l, m, 1) + & 1556 ncoset(l)*ncoset(m)*(ncoset(l + m)*4 + ncoset(n)*8)), & 1557 KIND=int_8) 1558 END FUNCTION 1559 1560! ************************************************************************************************** 1561!> \brief ... 1562!> \param S_G ... 1563!> \param l_max ... 1564!> \param m_max ... 1565!> \param n_max ... 1566!> \param R1 ... 1567!> \param R2 ... 1568!> \param alpha ... 1569!> \param beta ... 1570!> \param gamma ... 1571!> \param hmat ... 1572!> \param h_inv ... 1573!> \param vol ... 1574!> \param R_c ... 1575!> \param G_c ... 1576!> \param R_rad ... 1577!> \param G_rad ... 1578! ************************************************************************************************** 1579 PURE SUBROUTINE pgf_sum_product_3c_gspace_3d(S_G, l_max, m_max, n_max, R1, R2, alpha, beta, gamma, & 1580 hmat, h_inv, vol, R_c, G_c, R_rad, G_rad) 1581 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_G 1582 INTEGER, INTENT(IN) :: l_max, m_max, n_max 1583 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R1, R2 1584 REAL(KIND=dp), INTENT(IN) :: alpha, beta, gamma 1585 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv 1586 REAL(KIND=dp), INTENT(IN) :: vol 1587 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R_c, G_c 1588 REAL(KIND=dp), INTENT(IN) :: R_rad, G_rad 1589 1590 COMPLEX(KIND=dp) :: exp_tot 1591 INTEGER :: gx, gy, gz, k, l, lco, lnco, lx, ly, lz, & 1592 m, mco, n, nco 1593 COMPLEX(KIND=dp), & 1594 DIMENSION(ncoset(m_max), ncoset(n_max)) :: S_R 1595 COMPLEX(KIND=dp), DIMENSION(ncoset(l_max), ncoset( & 1596 m_max), ncoset(n_max)) :: S_G_c 1597 REAL(KIND=dp) :: G_rad_sq, G_sq, R_rad_sq 1598 REAL(KIND=dp), DIMENSION(3) :: G, G_x, G_y, G_z 1599 REAL(KIND=dp), DIMENSION(3, 0:l_max) :: G_pow_l 1600 REAL(KIND=dp), DIMENSION(3, 3) :: ht 1601 REAL(KIND=dp), DIMENSION(ncoset(l_max)) :: S_G_c_l 1602 1603 S_G_c(:, :, :) = 0.0_dp 1604 1605 G_rad_sq = G_rad**2 1606 R_rad_sq = R_rad**2 1607 1608 lnco = ncoset(l_max) 1609 1610 ht = twopi*TRANSPOSE(h_inv) 1611 DO gx = -FLOOR(G_c(1)), FLOOR(G_c(1)) 1612 G_x = ht(:, 1)*gx 1613 DO gy = -FLOOR(G_c(2)), FLOOR(G_c(2)) 1614 G_y = ht(:, 2)*gy 1615 DO gz = -FLOOR(G_c(3)), FLOOR(G_c(3)) 1616 G_z = ht(:, 3)*gz 1617 G = G_x + G_y + G_z 1618 G_sq = G(1)**2 + G(2)**2 + G(3)**2 1619 IF (G_sq .GT. G_rad_sq) CYCLE 1620 1621 exp_tot = EXP(-alpha*G_sq)*EXP(gaussi*DOT_PRODUCT(G, R1)) ! cost: exp_w*2 flops 1622 1623 DO k = 1, 3 1624 G_pow_l(k, 0) = 1.0_dp 1625 DO l = 1, l_max 1626 G_pow_l(k, l) = G_pow_l(k, l - 1)*G(k) 1627 ENDDO 1628 ENDDO 1629 1630 CALL pgf_sum_product_3c_rspace_3d(S_R, m_max, n_max, G, R2, beta, gamma, hmat, h_inv, vol, R_c, R_rad_sq) 1631 1632 DO lco = 1, ncoset(l_max) 1633 CALL get_l(lco, l, lx, ly, lz) 1634 S_G_c_l(lco) = G_pow_l(1, lx)*G_pow_l(2, ly)*G_pow_l(3, lz)*(-1.0_dp)**l 1635 ENDDO 1636 1637 DO nco = 1, ncoset(n_max) 1638 CALL get_l(nco, n) 1639 DO mco = 1, ncoset(m_max) 1640 CALL get_l(mco, m) 1641 DO lco = 1, ncoset(l_max) 1642 CALL get_l(lco, l) 1643 S_G_c(lco, mco, nco) = & ! cost: 8 flops 1644 S_G_c(lco, mco, nco) + & 1645 S_G_c_l(lco)* & 1646 exp_tot*i_pow(l + m + n)*(-1.0_dp)**m*S_R(mco, nco) 1647 ENDDO 1648 ENDDO 1649 ENDDO 1650 ENDDO 1651 ENDDO 1652 ENDDO 1653 S_G = REAL(S_G_c, KIND=dp)/vol**2 1654 1655 END SUBROUTINE pgf_sum_product_3c_gspace_3d 1656 1657! ************************************************************************************************** 1658!> \brief ... 1659!> \param S_R ... 1660!> \param l_max ... 1661!> \param m_max ... 1662!> \param G ... 1663!> \param R ... 1664!> \param alpha ... 1665!> \param beta ... 1666!> \param hmat ... 1667!> \param h_inv ... 1668!> \param vol ... 1669!> \param R_c ... 1670!> \param R_rad_sq ... 1671! ************************************************************************************************** 1672 PURE SUBROUTINE pgf_sum_product_3c_rspace_3d(S_R, l_max, m_max, G, R, alpha, beta, hmat, h_inv, vol, R_c, R_rad_sq) 1673 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: S_R 1674 INTEGER, INTENT(IN) :: l_max, m_max 1675 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: G, R 1676 REAL(KIND=dp), INTENT(IN) :: alpha, beta 1677 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv 1678 REAL(KIND=dp), INTENT(IN) :: vol 1679 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: R_c 1680 REAL(KIND=dp), INTENT(IN) :: R_rad_sq 1681 1682 COMPLEX(KIND=dp) :: exp_tot 1683 INTEGER :: k, l, lco, lx, ly, lz, m, mco, mx, my, & 1684 mz, sx, sy, sz, t, tco, tx, ty, tz 1685 COMPLEX(KIND=dp), DIMENSION(ncoset(l_max + m_max)) :: S_R_t 1686 REAL(KIND=dp) :: c1, c2, Rp_sq 1687 REAL(KIND=dp), & 1688 DIMENSION(-1:l_max + m_max + 1, -1:l_max, -1:m_max) :: E1, E2, E3 1689 REAL(KIND=dp), DIMENSION(3) :: R_l, R_r, Rp, Rx, Ry, Rz, s_shift 1690 REAL(KIND=dp), DIMENSION(3, 0:l_max + m_max) :: R_pow_t 1691 1692 c1 = 0.25_dp/(alpha + beta) 1693 c2 = alpha/(alpha + beta) 1694 1695 S_R_t(:) = 0.0_dp 1696 S_R(:, :) = 0.0_dp 1697 1698 s_shift = MATMUL(h_inv, R) 1699 R_l = -R_c + s_shift 1700 R_r = R_c + s_shift 1701 1702 DO sx = CEILING(R_l(1)), FLOOR(R_r(1)) 1703 Rx = hmat(:, 1)*sx 1704 DO sy = CEILING(R_l(2)), FLOOR(R_r(2)) 1705 Ry = hmat(:, 2)*sy 1706 DO sz = CEILING(R_l(3)), FLOOR(R_r(3)) 1707 Rz = hmat(:, 3)*sz 1708 Rp = Rx + Ry + Rz - R 1709 Rp_sq = Rp(1)**2 + Rp(2)**2 + Rp(3)**2 1710 IF (Rp_sq .GT. R_rad_sq) CYCLE 1711 1712 exp_tot = EXP(-c1*Rp_sq)*EXP(-gaussi*c2*DOT_PRODUCT(Rp, G)) ! cost: exp_w*2 flops 1713 DO k = 1, 3 1714 R_pow_t(k, 0) = 1.0_dp 1715 DO t = 1, l_max + m_max 1716 R_pow_t(k, t) = R_pow_t(k, t - 1)*Rp(k) 1717 ENDDO 1718 ENDDO 1719 1720 DO tco = 1, ncoset(l_max + m_max) 1721 CALL get_l(tco, t, tx, ty, tz) 1722 S_R_t(tco) = S_R_t(tco) + & ! cost: 7 flops 1723 R_pow_t(1, tx)*R_pow_t(2, ty)*R_pow_t(3, tz)* & 1724 (-1.0_dp)**t*i_pow(t)*exp_tot 1725 ENDDO 1726 1727 ENDDO 1728 ENDDO 1729 ENDDO 1730 1731 CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, alpha, beta, G(1), 0.0_dp, 1, E1) 1732 CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, alpha, beta, G(2), 0.0_dp, 1, E2) 1733 CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, alpha, beta, G(3), 0.0_dp, 1, E3) 1734 1735 DO mco = 1, ncoset(m_max) 1736 CALL get_l(mco, m, mx, my, mz) 1737 DO lco = 1, ncoset(l_max) 1738 CALL get_l(lco, l, lx, ly, lz) 1739 DO tx = 0, lx + mx 1740 DO ty = 0, ly + my 1741 DO tz = 0, lz + mz 1742 tco = coset(tx, ty, tz) 1743 S_R(lco, mco) = S_R(lco, mco) + & ! cost: 4 flops 1744 E1(tx, lx, mx)*E2(ty, ly, my)*E3(tz, lz, mz)*S_R_t(tco) 1745 1746 ENDDO 1747 ENDDO 1748 ENDDO 1749 ENDDO 1750 ENDDO 1751 1752 S_R(:, :) = S_R(:, :)*vol/(twopi)**3*(pi/(alpha + beta))**1.5_dp 1753 1754 END SUBROUTINE pgf_sum_product_3c_rspace_3d 1755 1756! ************************************************************************************************** 1757!> \brief Roughly estimated number of floating point operations 1758!> \param ns_R1 ... 1759!> \param ns_R2 ... 1760!> \param l ... 1761!> \param m ... 1762!> \param n ... 1763!> \return ... 1764! ************************************************************************************************** 1765 PURE FUNCTION nsum_3c_rspace_3d(ns_R1, ns_R2, l, m, n) 1766 REAL(KIND=dp), INTENT(IN) :: ns_R1, ns_R2 1767 INTEGER, INTENT(IN) :: l, m, n 1768 INTEGER(KIND=int_8) :: nsum_3c_rspace_3d 1769 1770 nsum_3c_rspace_3d = & 1771 NINT( & 1772 ns_R1*( & 1773 ns_R2*(exp_w + ncoset(l + m + n)*4) + & 1774 3*nsum_gaussian_overlap(l, m, 2) + & 1775 ncoset(m)*ncoset(l)*( & 1776 ncoset(l + m)*2 + ncoset(n)*ncoset(l + m)*4)), & 1777 KIND=int_8) 1778 1779 END FUNCTION 1780 1781! ************************************************************************************************** 1782!> \brief ... 1783!> \param S_R ... 1784!> \param l_max ... 1785!> \param m_max ... 1786!> \param n_max ... 1787!> \param RA ... 1788!> \param RB ... 1789!> \param RC ... 1790!> \param zeta ... 1791!> \param zetb ... 1792!> \param zetc ... 1793!> \param a_mm ... 1794!> \param hmat ... 1795!> \param h_inv ... 1796!> \param R_c ... 1797!> \param R_rad ... 1798! ************************************************************************************************** 1799 SUBROUTINE pgf_sum_3c_rspace_3d(S_R, l_max, m_max, n_max, RA, RB, RC, zeta, zetb, zetc, a_mm, & 1800 hmat, h_inv, R_c, R_rad) 1801 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: S_R 1802 INTEGER, INTENT(IN) :: l_max, m_max, n_max 1803 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: RA, RB, RC 1804 REAL(KIND=dp), INTENT(IN) :: zeta, zetb, zetc, a_mm 1805 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat, h_inv 1806 REAL(KIND=dp), DIMENSION(2, 3), INTENT(IN) :: R_c 1807 REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: R_rad 1808 1809 INTEGER :: k, l, lco, lx, ly, lz, m, mco, mx, my, & 1810 mz, n, nco, nx, ny, nz, s1x, s1y, s1z, & 1811 s2x, s2y, s2z, t, tco, tnco, ttco, & 1812 ttx, tty, ttz, tx, ty, tz 1813 REAL(KIND=dp) :: alpha, exp_tot, R1_sq, R_sq 1814 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_to_c 1815 REAL(KIND=dp), & 1816 DIMENSION(-1:l_max + m_max + 1, -1:l_max, -1:m_max) :: E1, E2, E3 1817 REAL(KIND=dp), DIMENSION(2) :: R_rad_sq 1818 REAL(KIND=dp), DIMENSION(3) :: R, R1, R1_l, R1_r, R1_tmp, R1x, R1y, & 1819 R1z, R2_l, R2_r, R2x, R2y, R2z, & 1820 R_offset, R_tmp, s1_shift, s2_shift 1821 REAL(KIND=dp), DIMENSION(3, 0:l_max + m_max + n_max) :: R_pow_t 1822 REAL(KIND=dp), DIMENSION(ncoset(l_max + m_max), & 1823 ncoset(l_max), ncoset(m_max)) :: Eco 1824 REAL(KIND=dp), & 1825 DIMENSION(ncoset(l_max + m_max + n_max)) :: S_R_t 1826 REAL(KIND=dp), DIMENSION(ncoset(l_max + m_max + n_max) & 1827 , ncoset(l_max + m_max + n_max)) :: h_to_c_tco 1828 1829 alpha = 1.0_dp/((zeta + zetb + zetc)/((zeta + zetb)*zetc) + 4.0_dp*a_mm) 1830 1831 Eco(:, :, :) = 0.0_dp 1832 S_R(:, :, :) = 0.0_dp 1833 h_to_c_tco(:, :) = 0.0_dp 1834 1835 R_offset = RC - (zeta*RA + zetb*RB)/(zeta + zetb) 1836 1837 CALL create_hermite_to_cartesian(alpha, l_max + m_max + n_max, h_to_c) 1838 1839 DO tco = 1, ncoset(l_max + m_max + n_max) 1840 CALL get_l(tco, t, tx, ty, tz) 1841 DO ttx = 0, tx 1842 DO tty = 0, ty 1843 DO ttz = 0, tz 1844 ttco = coset(ttx, tty, ttz) 1845 h_to_c_tco(ttco, tco) = h_to_c(ttx, tx)*h_to_c(tty, ty)*h_to_c(ttz, tz) 1846 ENDDO 1847 ENDDO 1848 ENDDO 1849 ENDDO 1850 1851 s1_shift = MATMUL(h_inv, RA - RB) 1852 R1_l = -R_c(1, :) + s1_shift 1853 R1_r = R_c(1, :) + s1_shift 1854 1855 R_rad_sq = R_rad**2 1856 1857 DO s1x = CEILING(R1_l(1)), FLOOR(R1_r(1)) 1858 R1x = hmat(:, 1)*s1x 1859 DO s1y = CEILING(R1_l(2)), FLOOR(R1_r(2)) 1860 R1y = hmat(:, 2)*s1y 1861 DO s1z = CEILING(R1_l(3)), FLOOR(R1_r(3)) 1862 R1z = hmat(:, 3)*s1z 1863 R1 = R1x + R1y + R1z 1864 S_R_t(:) = 0.0_dp 1865 R1_tmp = R1 - (RA - RB) 1866 R1_sq = R1_tmp(1)**2 + R1_tmp(2)**2 + R1_tmp(3)**2 1867 1868 IF (R1_sq .GT. R_rad_sq(1)) CYCLE 1869 1870 R_tmp = R_offset + R1*zeta/(zeta + zetb) 1871 s2_shift = -MATMUL(h_inv, R_tmp) 1872 R2_l = -R_c(2, :) + s2_shift 1873 R2_r = R_c(2, :) + s2_shift 1874 DO s2x = CEILING(R2_l(1)), FLOOR(R2_r(1)) 1875 R2x = hmat(:, 1)*s2x 1876 DO s2y = CEILING(R2_l(2)), FLOOR(R2_r(2)) 1877 R2y = hmat(:, 2)*s2y 1878 DO s2z = CEILING(R2_l(3)), FLOOR(R2_r(3)) 1879 R2z = hmat(:, 3)*s2z 1880 R = R_tmp + R2x + R2y + R2z 1881 R_sq = R(1)**2 + R(2)**2 + R(3)**2 1882 1883 IF (R_sq .GT. R_rad_sq(2)) CYCLE 1884 1885 exp_tot = EXP(-alpha*R_sq) ! cost: exp_w flops 1886 DO k = 1, 3 1887 R_pow_t(k, 0) = 1.0_dp 1888 DO t = 1, l_max + m_max + n_max 1889 R_pow_t(k, t) = R_pow_t(k, t - 1)*R(k) 1890 ENDDO 1891 ENDDO 1892 DO tco = 1, ncoset(l_max + m_max + n_max) 1893 CALL get_l(tco, t, tx, ty, tz) 1894 S_R_t(tco) = S_R_t(tco) + R_pow_t(1, tx)*R_pow_t(2, ty)*R_pow_t(3, tz)*exp_tot ! cost: 4 flops 1895 ENDDO 1896 ENDDO 1897 ENDDO 1898 ENDDO 1899 1900 S_R_t(:) = MATMUL(TRANSPOSE(h_to_c_tco), S_R_t)*(alpha/pi)**1.5_dp 1901 1902 CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA(1) - R1(1), RB(1), 2, E1) 1903 CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA(2) - R1(2), RB(2), 2, E2) 1904 CALL create_gaussian_overlap_dist_to_hermite(l_max, m_max, zeta, zetb, RA(3) - R1(3), RB(3), 2, E3) 1905 1906 DO mco = 1, ncoset(m_max) 1907 CALL get_l(mco, m, mx, my, mz) 1908 DO lco = 1, ncoset(l_max) 1909 CALL get_l(lco, l, lx, ly, lz) 1910 DO tx = 0, lx + mx 1911 DO ty = 0, ly + my 1912 DO tz = 0, lz + mz 1913 tco = coset(tx, ty, tz) 1914 Eco(tco, lco, mco) = E1(tx, lx, mx)*E2(ty, ly, my)*E3(tz, lz, mz) ! cost: 2 flops 1915 ENDDO 1916 ENDDO 1917 ENDDO 1918 ENDDO 1919 ENDDO 1920 1921 DO nco = 1, ncoset(n_max) 1922 CALL get_l(nco, n, nx, ny, nz) 1923 DO tco = 1, ncoset(l_max + m_max) 1924 CALL get_l(tco, t, tx, ty, tz) 1925 tnco = coset(tx + nx, ty + ny, tz + nz) 1926 S_R(:, :, nco) = S_R(:, :, nco) + & ! cost: 4 flops 1927 Eco(tco, :, :)* & 1928 (-1)**n*S_R_t(tnco) 1929 1930 ENDDO 1931 ENDDO 1932 ENDDO 1933 ENDDO 1934 ENDDO 1935 END SUBROUTINE pgf_sum_3c_rspace_3d 1936 1937! ************************************************************************************************** 1938!> \brief Compute bounding box for ellipsoid. This is needed in order to find summation bounds for 1939!> sphere for sums over non-orthogonal lattice vectors. 1940!> \param s_rad sphere radius 1941!> \param s_to_e sphere to ellipsoid trafo 1942!> \return ... 1943! ************************************************************************************************** 1944 PURE FUNCTION ellipsoid_bounds(s_rad, s_to_e) 1945 REAL(KIND=dp), INTENT(IN) :: s_rad 1946 REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: s_to_e 1947 REAL(KIND=dp), DIMENSION(3) :: ellipsoid_bounds 1948 1949 INTEGER :: i_xyz 1950 1951 DO i_xyz = 1, 3 1952 ellipsoid_bounds(i_xyz) = SQRT(s_to_e(i_xyz, 1)**2 + s_to_e(i_xyz, 2)**2 + s_to_e(i_xyz, 3)**2)*s_rad 1953 ENDDO 1954 1955 END FUNCTION ellipsoid_bounds 1956 1957! ************************************************************************************************** 1958!> \brief Roughly estimated number of floating point operations 1959!> \param l ... 1960!> \param m ... 1961!> \param H_or_C_product ... 1962!> \return ... 1963! ************************************************************************************************** 1964 PURE FUNCTION nsum_gaussian_overlap(l, m, H_or_C_product) 1965 INTEGER, INTENT(IN) :: l, m, H_or_C_product 1966 INTEGER :: nsum_gaussian_overlap 1967 1968 INTEGER :: loop 1969 1970 nsum_gaussian_overlap = exp_w 1971 loop = (m + 1)*(l + 1)*(m + l + 2) 1972 IF (H_or_C_product .EQ. 1) THEN 1973 nsum_gaussian_overlap = nsum_gaussian_overlap + loop*16 1974 ELSE 1975 nsum_gaussian_overlap = nsum_gaussian_overlap + loop*32 1976 ENDIF 1977 END FUNCTION 1978 1979! ************************************************************************************************** 1980!> \brief ... 1981!> \param lco ... 1982!> \param l ... 1983!> \param lx ... 1984!> \param ly ... 1985!> \param lz ... 1986! ************************************************************************************************** 1987 PURE ELEMENTAL SUBROUTINE get_l(lco, l, lx, ly, lz) 1988 INTEGER, INTENT(IN) :: lco 1989 INTEGER, INTENT(OUT) :: l 1990 INTEGER, INTENT(OUT), OPTIONAL :: lx, ly, lz 1991 1992 l = SUM(indco(:, lco)) 1993 IF (PRESENT(lx)) lx = indco(1, lco) 1994 IF (PRESENT(ly)) ly = indco(2, lco) 1995 IF (PRESENT(lz)) lz = indco(3, lco) 1996 END SUBROUTINE 1997 1998! ************************************************************************************************** 1999!> \brief ... 2000!> \param i ... 2001!> \return ... 2002! ************************************************************************************************** 2003 PURE ELEMENTAL FUNCTION i_pow(i) 2004 INTEGER, INTENT(IN) :: i 2005 COMPLEX(KIND=dp) :: i_pow 2006 2007 COMPLEX(KIND=dp), DIMENSION(0:3), PARAMETER :: & 2008 ip = (/(1.0_dp, 0.0_dp), (0.0_dp, 1.0_dp), (-1.0_dp, 0.0_dp), (0.0_dp, -1.0_dp)/) 2009 2010 i_pow = ip(MOD(i, 4)) 2011 2012 END FUNCTION 2013 2014END MODULE eri_mme_lattice_summation 2015