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