1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5MODULE grid_collocate 6 USE cell_types, ONLY: cell_type 7 USE cube_utils, ONLY: compute_cube_center,& 8 cube_info_type,& 9 return_cube,& 10 return_cube_nonortho 11 USE d3_poly, ONLY: poly_cp2k2d3 12 USE gauss_colloc, ONLY: collocGauss 13 USE grid_modify_pab_block, ONLY: prepare_adb_m_dab,& 14 prepare_ardb_m_darb,& 15 prepare_dab_p_adb,& 16 prepare_dadb,& 17 prepare_diadib,& 18 prepare_diiadiib,& 19 prepare_dijadijb 20 USE kinds, ONLY: dp,& 21 int_8 22 USE mathconstants, ONLY: fac 23 USE orbital_pointers, ONLY: coset,& 24 ncoset 25 USE realspace_grid_types, ONLY: realspace_grid_type 26#include "../base/base_uses.f90" 27 28 IMPLICIT NONE 29 30 PRIVATE 31 32 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'grid_collocate' 33 34 PUBLIC :: collocate_pgf_product 35 36 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_AB = 100 37 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DADB = 200 38 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ADBmDAB_X = 301 39 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ADBmDAB_Y = 302 40 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ADBmDAB_Z = 303 41 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ARDBmDARB_XX = 411 42 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ARDBmDARB_XY = 412 43 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ARDBmDARB_XZ = 413 44 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ARDBmDARB_YX = 421 45 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ARDBmDARB_YY = 422 46 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ARDBmDARB_YZ = 423 47 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ARDBmDARB_ZX = 431 48 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ARDBmDARB_ZY = 432 49 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_ARDBmDARB_ZZ = 433 50 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DABpADB_X = 501 51 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DABpADB_Y = 502 52 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DABpADB_Z = 503 53 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DX = 601 54 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DY = 602 55 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DZ = 603 56 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DXDY = 701 57 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DYDZ = 702 58 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DZDX = 703 59 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DXDX = 801 60 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DYDY = 802 61 INTEGER, PARAMETER, PUBLIC :: GRID_FUNC_DZDZ = 803 62 63CONTAINS 64 65! ************************************************************************************************** 66!> \brief low level collocation of primitive gaussian functions 67!> \param la_max ... 68!> \param zeta ... 69!> \param la_min ... 70!> \param lb_max ... 71!> \param zetb ... 72!> \param lb_min ... 73!> \param ra ... 74!> \param rab ... 75!> \param scale ... 76!> \param pab ... 77!> \param o1 ... 78!> \param o2 ... 79!> \param rsgrid ... 80!> \param cell ... 81!> \param cube_info ... 82!> \param ga_gb_function ... 83!> \param radius ... 84!> \param use_subpatch ... 85!> \param subpatch_pattern ... 86! ************************************************************************************************** 87 SUBROUTINE collocate_pgf_product(la_max, zeta, la_min, & 88 lb_max, zetb, lb_min, & 89 ra, rab, scale, pab, o1, o2, & 90 rsgrid, cell, cube_info, & 91 ga_gb_function, radius, & 92 use_subpatch, subpatch_pattern) 93 94 INTEGER, INTENT(IN) :: la_max 95 REAL(KIND=dp), INTENT(IN) :: zeta 96 INTEGER, INTENT(IN) :: la_min, lb_max 97 REAL(KIND=dp), INTENT(IN) :: zetb 98 INTEGER, INTENT(IN) :: lb_min 99 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rab 100 REAL(KIND=dp), INTENT(IN) :: scale 101 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab 102 INTEGER, INTENT(IN) :: o1, o2 103 TYPE(realspace_grid_type) :: rsgrid 104 TYPE(cell_type), POINTER :: cell 105 TYPE(cube_info_type), INTENT(IN) :: cube_info 106 INTEGER, INTENT(IN) :: ga_gb_function 107 REAL(KIND=dp), INTENT(IN) :: radius 108 LOGICAL, OPTIONAL :: use_subpatch 109 INTEGER(KIND=int_8), INTENT(IN), OPTIONAL :: subpatch_pattern 110 111 CHARACTER(len=*), PARAMETER :: routineN = 'collocate_pgf_product', & 112 routineP = moduleN//':'//routineN 113 114 INTEGER :: ider1, ider2, idir, ir, la_max_local, la_min_local, lb_max_local, lb_min_local, & 115 lxa, lxb, lya, lyb, lza, lzb, o1_local, o2_local 116 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab_local 117 118 ! it's a choice to compute lX_min/max, pab here, 119 ! this way we get the same radius as we use for the corresponding density 120 SELECT CASE (ga_gb_function) 121 CASE (GRID_FUNC_DADB) 122 la_max_local = la_max + 1 123 la_min_local = MAX(la_min - 1, 0) 124 lb_max_local = lb_max + 1 125 lb_min_local = MAX(lb_min - 1, 0) 126 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 127 ! is equivalent to mapping pab with 0.5 * (nabla pgf_a) . (nabla pgf_b) 128 ! (ddx pgf_a ) (ddx pgf_b) = (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x})*(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) 129 ! cleaner would possibly be to touch pzyx directly (avoiding the following allocate) 130 ALLOCATE (pab_local(ncoset(la_max_local), ncoset(lb_max_local))) 131 pab_local = 0.0_dp 132 DO lxa = 0, la_max 133 DO lxb = 0, lb_max 134 DO lya = 0, la_max - lxa 135 DO lyb = 0, lb_max - lxb 136 DO lza = MAX(la_min - lxa - lya, 0), la_max - lxa - lya 137 DO lzb = MAX(lb_min - lxb - lyb, 0), lb_max - lxb - lyb 138 139 ! this element of pab results in 12 elements of pab_local 140 CALL prepare_dadb(pab_local, pab, lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 141 142 ENDDO 143 ENDDO 144 ENDDO 145 ENDDO 146 ENDDO 147 ENDDO 148 o1_local = 0 149 o2_local = 0 150 pab_local = pab_local*0.5_dp 151 CASE (GRID_FUNC_ADBmDAB_X, GRID_FUNC_ADBmDAB_Y, GRID_FUNC_ADBmDAB_Z) 152 idir = MODULO(ga_gb_function, 10) 153 la_max_local = la_max + 1 154 la_min_local = MAX(la_min - 1, 0) 155 lb_max_local = lb_max + 1 156 lb_min_local = MAX(lb_min - 1, 0) 157 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 158 ! is equivalent to mapping pab with 159 ! pgf_a (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) pgf_b 160 ! ( pgf_a ) (ddx pgf_b) - (ddx pgf_a)( pgf_b ) = 161 ! pgf_a *(lbx pgf_{b-1x} - 2*zetb*pgf_{b+1x}) - 162 ! (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_b 163 164 ALLOCATE (pab_local(ncoset(la_max_local), ncoset(lb_max_local))) 165 pab_local = 0.0_dp 166 DO lxa = 0, la_max 167 DO lxb = 0, lb_max 168 DO lya = 0, la_max - lxa 169 DO lyb = 0, lb_max - lxb 170 DO lza = MAX(la_min - lxa - lya, 0), la_max - lxa - lya 171 DO lzb = MAX(lb_min - lxb - lyb, 0), lb_max - lxb - lyb 172 ! this element of pab results in 4 elements of pab_local 173 CALL prepare_adb_m_dab(pab_local, pab, idir, & 174 lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 175 END DO 176 END DO 177 END DO 178 END DO 179 END DO 180 END DO 181 o1_local = 0 182 o2_local = 0 183 CASE (GRID_FUNC_DABpADB_X, GRID_FUNC_DABpADB_Y, GRID_FUNC_DABpADB_Z) 184 idir = MODULO(ga_gb_function, 10) 185 la_max_local = la_max + 1 186 la_min_local = MAX(la_min - 1, 0) 187 lb_max_local = lb_max + 1 188 lb_min_local = MAX(lb_min - 1, 0) 189 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 190 ! is equivalent to mapping pab with 191 ! pgf_a (nabla_{idir} pgf_b) + (nabla_{idir} pgf_a) pgf_b 192 ! ( pgf_a ) (ddx pgf_b) + (ddx pgf_a)( pgf_b ) = 193 ! pgf_a *(lbx pgf_{b-1x} + 2*zetb*pgf_{b+1x}) + 194 ! (lax pgf_{a-1x} + 2*zeta*pgf_{a+1x}) pgf_b 195 196 ALLOCATE (pab_local(ncoset(la_max_local), ncoset(lb_max_local))) 197 pab_local = 0.0_dp 198 DO lxa = 0, la_max 199 DO lxb = 0, lb_max 200 DO lya = 0, la_max - lxa 201 DO lyb = 0, lb_max - lxb 202 DO lza = MAX(la_min - lxa - lya, 0), la_max - lxa - lya 203 DO lzb = MAX(lb_min - lxb - lyb, 0), lb_max - lxb - lyb 204 ! this element of pab results in 4 elements of pab_local 205 CALL prepare_dab_p_adb(pab_local, pab, idir, & 206 lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 207 END DO 208 END DO 209 END DO 210 END DO 211 END DO 212 END DO 213 o1_local = 0 214 o2_local = 0 215 CASE (GRID_FUNC_DX, GRID_FUNC_DY, GRID_FUNC_DZ) 216 ider1 = MODULO(ga_gb_function, 10) 217 la_max_local = la_max + 1 218 la_min_local = MAX(la_min - 1, 0) 219 lb_max_local = lb_max + 1 220 lb_min_local = MAX(lb_min - 1, 0) 221 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 222 ! is equivalent to mapping pab with 223 ! d_{ider1} pgf_a d_{ider1} pgf_b 224 ! dx pgf_a dx pgf_b = 225 ! (lax pgf_{a-1x})*(lbx pgf_{b-1x}) - 2*zetb*lax*pgf_{a-1x}*pgf_{b+1x} - 226 ! lbx pgf_{b-1x}*2*zeta*pgf_{a+1x}+ 4*zeta*zetab*pgf_{a+1x}pgf_{b+1x} 227 228 ALLOCATE (pab_local(ncoset(la_max_local), ncoset(lb_max_local))) 229 pab_local = 0.0_dp 230 DO lxa = 0, la_max 231 DO lxb = 0, lb_max 232 DO lya = 0, la_max - lxa 233 DO lyb = 0, lb_max - lxb 234 DO lza = MAX(la_min - lxa - lya, 0), la_max - lxa - lya 235 DO lzb = MAX(lb_min - lxb - lyb, 0), lb_max - lxb - lyb 236 ! this element of pab results in 4 elements of pab_local 237 CALL prepare_dIadIb(pab_local, pab, ider1, & 238 lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 239 END DO 240 END DO 241 END DO 242 END DO 243 END DO 244 END DO 245 o1_local = 0 246 o2_local = 0 247 CASE (GRID_FUNC_DXDY, GRID_FUNC_DYDZ, GRID_FUNC_DZDX) 248 ider1 = MODULO(ga_gb_function, 10) 249 ider2 = ider1 + 1 250 IF (ider2 > 3) ider2 = ider1 - 2 251 la_max_local = la_max + 2 252 la_min_local = MAX(la_min - 2, 0) 253 lb_max_local = lb_max + 2 254 lb_min_local = MAX(lb_min - 2, 0) 255 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 256 ! is equivalent to mapping pab with 257 ! d_{ider1} pgf_a d_{ider1} pgf_b 258 ALLOCATE (pab_local(ncoset(la_max_local), ncoset(lb_max_local))) 259 pab_local = 0.0_dp 260 DO lxa = 0, la_max 261 DO lxb = 0, lb_max 262 DO lya = 0, la_max - lxa 263 DO lyb = 0, lb_max - lxb 264 DO lza = MAX(la_min - lxa - lya, 0), la_max - lxa - lya 265 DO lzb = MAX(lb_min - lxb - lyb, 0), lb_max - lxb - lyb 266 ! this element of pab results in 16 elements of pab_local 267 CALL prepare_dijadijb(pab_local, pab, ider1, ider2, & 268 lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 269 END DO 270 END DO 271 END DO 272 END DO 273 END DO 274 END DO 275 o1_local = 0 276 o2_local = 0 277 CASE (GRID_FUNC_DXDX, GRID_FUNC_DYDY, GRID_FUNC_DZDZ) 278 ider1 = MODULO(ga_gb_function, 10) 279 la_max_local = la_max + 2 280 la_min_local = MAX(la_min - 2, 0) 281 lb_max_local = lb_max + 2 282 lb_min_local = MAX(lb_min - 2, 0) 283 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 284 ! is equivalent to mapping pab with 285 ! dd_{ider1} pgf_a dd_{ider1} pgf_b 286 287 ALLOCATE (pab_local(ncoset(la_max_local), ncoset(lb_max_local))) 288 pab_local = 0.0_dp 289 DO lxa = 0, la_max 290 DO lxb = 0, lb_max 291 DO lya = 0, la_max - lxa 292 DO lyb = 0, lb_max - lxb 293 DO lza = MAX(la_min - lxa - lya, 0), la_max - lxa - lya 294 DO lzb = MAX(lb_min - lxb - lyb, 0), lb_max - lxb - lyb 295 ! this element of pab results in 9 elements of pab_local 296 CALL prepare_diiadiib(pab_local, pab, ider1, & 297 lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 298 END DO 299 END DO 300 END DO 301 END DO 302 END DO 303 END DO 304 o1_local = 0 305 o2_local = 0 306 CASE (GRID_FUNC_ARDBmDARB_XX, GRID_FUNC_ARDBmDARB_XY, GRID_FUNC_ARDBmDARB_XZ, & 307 GRID_FUNC_ARDBmDARB_YX, GRID_FUNC_ARDBmDARB_YY, GRID_FUNC_ARDBmDARB_YZ, & 308 GRID_FUNC_ARDBmDARB_ZX, GRID_FUNC_ARDBmDARB_ZY, GRID_FUNC_ARDBmDARB_ZZ) 309 ir = MODULO(ga_gb_function, 10) 310 idir = MODULO(ga_gb_function - ir, 100)/10 311 la_max_local = la_max + 1 312 la_min_local = MAX(la_min - 1, 0) 313 lb_max_local = lb_max + 2 314 lb_min_local = MAX(lb_min - 1, 0) 315 ! create a new pab_local so that mapping pab_local with pgf_a pgf_b 316 ! is equivalent to mapping pab with 317 ! pgf_a (r-Rb)_{ir} (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) (r-Rb)_{ir} pgf_b 318 ! ( pgf_a )(r-Rb)_{ir} (ddx pgf_b) - (ddx pgf_a) (r-Rb)_{ir} ( pgf_b ) = 319 ! pgf_a *(lbx pgf_{b-1x+1ir} - 2*zetb*pgf_{b+1x+1ir}) - 320 ! (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_{b+1ir} 321 322 ALLOCATE (pab_local(ncoset(la_max_local), ncoset(lb_max_local))) 323 pab_local = 0.0_dp 324 DO lxa = 0, la_max 325 DO lxb = 0, lb_max 326 DO lya = 0, la_max - lxa 327 DO lyb = 0, lb_max - lxb 328 DO lza = MAX(la_min - lxa - lya, 0), la_max - lxa - lya 329 DO lzb = MAX(lb_min - lxb - lyb, 0), lb_max - lxb - lyb 330 331 ! this element of pab results in 4 elements of pab_local 332 CALL prepare_ardb_m_darb(pab_local, pab, idir, ir, & 333 lxa, lya, lza, lxb, lyb, lzb, o1, o2, zeta, zetb) 334 END DO 335 END DO 336 END DO 337 END DO 338 END DO 339 END DO 340 o1_local = 0 341 o2_local = 0 342 CASE (GRID_FUNC_AB) 343 la_max_local = la_max 344 la_min_local = la_min 345 lb_max_local = lb_max 346 lb_min_local = lb_min 347 pab_local => pab 348 o1_local = o1 349 o2_local = o2 350 CASE DEFAULT 351 CPASSERT(.FALSE.) 352 END SELECT 353 354 CALL collocate_pgf_product_part2(la_max_local, zeta, la_min_local, & 355 lb_max_local, zetb, lb_min_local, & 356 ra, rab, scale, pab_local, & 357 o1_local, o2_local, & 358 rsgrid, cell, cube_info, & 359 radius, & 360 use_subpatch, subpatch_pattern, & 361 lp=la_max_local + lb_max_local, & 362 lmax=MAX(la_max_local, lb_max_local)) 363 364 IF (ga_gb_function /= GRID_FUNC_AB) THEN 365 DEALLOCATE (pab_local) 366 ENDIF 367 368 END SUBROUTINE collocate_pgf_product 369 370! ************************************************************************************************** 371!> \brief After lp and lmax are determined they can be used to allocate arrays on the stack. 372!> \param la_max_local ... 373!> \param zeta ... 374!> \param la_min_local ... 375!> \param lb_max_local ... 376!> \param zetb ... 377!> \param lb_min_local ... 378!> \param ra ... 379!> \param rab ... 380!> \param scale ... 381!> \param pab_local ... 382!> \param o1_local ... 383!> \param o2_local ... 384!> \param rsgrid ... 385!> \param cell ... 386!> \param cube_info ... 387!> \param radius ... 388!> \param use_subpatch ... 389!> \param subpatch_pattern ... 390!> \param lp ... 391!> \param lmax ... 392! ************************************************************************************************** 393 SUBROUTINE collocate_pgf_product_part2(la_max_local, zeta, la_min_local, & 394 lb_max_local, zetb, lb_min_local, & 395 ra, rab, scale, pab_local, & 396 o1_local, o2_local, & 397 rsgrid, cell, cube_info, & 398 radius, & 399 use_subpatch, subpatch_pattern, & 400 lp, lmax) 401 402 INTEGER, INTENT(IN) :: la_max_local 403 REAL(KIND=dp), INTENT(IN) :: zeta 404 INTEGER, INTENT(IN) :: la_min_local, lb_max_local 405 REAL(KIND=dp), INTENT(IN) :: zetb 406 INTEGER, INTENT(IN) :: lb_min_local 407 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rab 408 REAL(KIND=dp), INTENT(IN) :: scale 409 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab_local 410 INTEGER, INTENT(IN) :: o1_local, o2_local 411 TYPE(realspace_grid_type) :: rsgrid 412 TYPE(cell_type), POINTER :: cell 413 TYPE(cube_info_type), INTENT(IN) :: cube_info 414 REAL(KIND=dp), INTENT(IN) :: radius 415 LOGICAL, OPTIONAL :: use_subpatch 416 INTEGER(KIND=int_8), OPTIONAL, INTENT(IN):: subpatch_pattern 417 INTEGER, INTENT(IN) :: lp, lmax 418 419 CHARACTER(len=*), PARAMETER :: routineN = 'collocate_pgf_product_part2', & 420 routineP = moduleN//':'//routineN 421 422 INTEGER :: cmax, gridbounds(2, 3), i, ico, icoef, ig, jco, k, l, lxp, lyp, lzp, lxpm, iaxis, & 423 length, lxa, lxb, lxy, lxyz, lya, lyb, lza, lzb, offset, start 424 INTEGER, DIMENSION(3) :: cubecenter, lb_cube, ng, ub_cube 425 INTEGER, DIMENSION(:), POINTER :: sphere_bounds 426 LOGICAL :: subpatch_collocate 427 REAL(KIND=dp) :: a, b, binomial_k_lxa, binomial_l_lxb, pg, prefactor, rpg, zetp, f, p_ele, & 428 t_exp_1, t_exp_2, t_exp_min_1, t_exp_min_2, t_exp_plus_1, t_exp_plus_2 429 REAL(KIND=dp), DIMENSION(3) :: dr, roffset, rb, rp 430 REAL(KIND=dp), DIMENSION(:, :, :), & 431 POINTER :: grid 432 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: map 433 REAL(kind=dp), DIMENSION(0:lp, 0:lmax, 0:lmax, 3) :: alpha 434 REAL(kind=dp), DIMENSION(((lp + 1)*(lp + 2)*(lp + 3))/6) :: coef_xyz 435 REAL(kind=dp), DIMENSION(((lp + 1)*(lp + 2))/2) :: coef_xyt 436 REAL(kind=dp), DIMENSION(0:lp) :: coef_xtt 437 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: pol_z 438 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: pol_y 439 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: pol_x 440 441 subpatch_collocate = .FALSE. 442 443 IF (PRESENT(use_subpatch)) THEN 444 IF (use_subpatch) THEN 445 subpatch_collocate = .TRUE. 446 CPASSERT(PRESENT(subpatch_pattern)) 447 ENDIF 448 ENDIF 449 450 ! check to avoid overflows 451 a = MAXVAL(ABS(rsgrid%desc%dh)) 452 a = 300._dp/(a*a) 453 ! CPASSERT(zetp < a) 454 455 a = MAXVAL(ABS(rsgrid%desc%dh)) 456 IF (radius .LT. a/2.0_dp) THEN 457 RETURN 458 END IF 459 460 zetp = zeta + zetb 461 f = zetb/zetp 462 rp(:) = ra(:) + f*rab(:) 463 rb(:) = ra(:) + rab(:) 464 prefactor = scale*EXP(-zeta*f*DOT_PRODUCT(rab, rab)) 465 466 ng(:) = rsgrid%desc%npts(:) 467 grid => rsgrid%r(:, :, :) 468 gridbounds(1, 1) = LBOUND(GRID, 1) 469 gridbounds(2, 1) = UBOUND(GRID, 1) 470 gridbounds(1, 2) = LBOUND(GRID, 2) 471 gridbounds(2, 2) = UBOUND(GRID, 2) 472 gridbounds(1, 3) = LBOUND(GRID, 3) 473 gridbounds(2, 3) = UBOUND(GRID, 3) 474 475! *** initialise the coefficient matrix, we transform the sum 476! 477! sum_{lxa,lya,lza,lxb,lyb,lzb} P_{lxa,lya,lza,lxb,lyb,lzb} * 478! (x-a_x)**lxa (y-a_y)**lya (z-a_z)**lza (x-b_x)**lxb (y-a_y)**lya (z-a_z)**lza 479! 480! into 481! 482! sum_{lxp,lyp,lzp} P_{lxp,lyp,lzp} (x-p_x)**lxp (y-p_y)**lyp (z-p_z)**lzp 483! 484! where p is center of the product gaussian, and lp = la_max + lb_max 485! (current implementation is l**7) 486! 487! compute polynomial expansion coefs -> (x-a)**lxa (x-b)**lxb -> sum_{ls} alpha(ls,lxa,lxb,1)*(x-p)**ls 488! 489! *** make the alpha matrix *** 490 alpha(:, :, :, :) = 0.0_dp 491 DO iaxis = 1, 3 492 DO lxa = 0, la_max_local 493 DO lxb = 0, lb_max_local 494 binomial_k_lxa = 1.0_dp 495 a = 1.0_dp 496 DO k = 0, lxa 497 binomial_l_lxb = 1.0_dp 498 b = 1.0_dp 499 DO l = 0, lxb 500 alpha(lxa - l + lxb - k, lxa, lxb, iaxis) = alpha(lxa - l + lxb - k, lxa, lxb, iaxis) + & 501 binomial_k_lxa*binomial_l_lxb*a*b 502 binomial_l_lxb = binomial_l_lxb*REAL(lxb - l, dp)/REAL(l + 1, dp) 503 b = b*(rp(iaxis) - (ra(iaxis) + rab(iaxis))) 504 ENDDO 505 binomial_k_lxa = binomial_k_lxa*REAL(lxa - k, dp)/REAL(k + 1, dp) 506 a = a*(-ra(iaxis) + rp(iaxis)) 507 ENDDO 508 ENDDO 509 ENDDO 510 ENDDO 511 512! 513! compute P_{lxp,lyp,lzp} given P_{lxa,lya,lza,lxb,lyb,lzb} and alpha(ls,lxa,lxb,1) 514! use a three step procedure 515! we don't store zeros, so counting is done using lxyz,lxy in order to have contiguous memory access in collocate_fast.F 516! 517 lxyz = 0 518 DO lzp = 0, lp 519 DO lyp = 0, lp - lzp 520 DO lxp = 0, lp - lzp - lyp 521 lxyz = lxyz + 1 522 coef_xyz(lxyz) = 0.0_dp 523 ENDDO 524 ENDDO 525 ENDDO 526 DO lzb = 0, lb_max_local 527 DO lza = 0, la_max_local 528 lxy = 0 529 DO lyp = 0, lp - lza - lzb 530 DO lxp = 0, lp - lza - lzb - lyp 531 lxy = lxy + 1 532 coef_xyt(lxy) = 0.0_dp 533 ENDDO 534 lxy = lxy + lza + lzb 535 ENDDO 536 DO lyb = 0, lb_max_local - lzb 537 DO lya = 0, la_max_local - lza 538 lxpm = (lb_max_local - lzb - lyb) + (la_max_local - lza - lya) 539 coef_xtt(0:lxpm) = 0.0_dp 540 DO lxb = MAX(lb_min_local - lzb - lyb, 0), lb_max_local - lzb - lyb 541 DO lxa = MAX(la_min_local - lza - lya, 0), la_max_local - lza - lya 542 ico = coset(lxa, lya, lza) 543 jco = coset(lxb, lyb, lzb) 544 p_ele = prefactor*pab_local(o1_local + ico, o2_local + jco) 545 DO lxp = 0, lxa + lxb 546 coef_xtt(lxp) = coef_xtt(lxp) + p_ele*alpha(lxp, lxa, lxb, 1) 547 ENDDO 548 ENDDO 549 ENDDO 550 lxy = 0 551 DO lyp = 0, lya + lyb 552 DO lxp = 0, lp - lza - lzb - lya - lyb 553 lxy = lxy + 1 554 coef_xyt(lxy) = coef_xyt(lxy) + alpha(lyp, lya, lyb, 2)*coef_xtt(lxp) 555 ENDDO 556 lxy = lxy + lza + lzb + lya + lyb - lyp 557 ENDDO 558 ENDDO 559 ENDDO 560 lxyz = 0 561 DO lzp = 0, lza + lzb 562 lxy = 0 563 DO lyp = 0, lp - lza - lzb 564 DO lxp = 0, lp - lza - lzb - lyp 565 lxy = lxy + 1; lxyz = lxyz + 1 566 coef_xyz(lxyz) = coef_xyz(lxyz) + alpha(lzp, lza, lzb, 3)*coef_xyt(lxy) 567 ENDDO 568 lxy = lxy + lza + lzb; lxyz = lxyz + lza + lzb - lzp 569 ENDDO 570 DO lyp = lp - lza - lzb + 1, lp - lzp 571 DO lxp = 0, lp - lyp - lzp 572 lxyz = lxyz + 1 573 ENDDO 574 ENDDO 575 ENDDO 576 ENDDO 577 ENDDO 578 579 IF (subpatch_collocate) THEN 580 CALL collocate_general_subpatch() 581 ELSE 582 IF (rsgrid%desc%orthorhombic) THEN 583 CALL collocate_ortho() 584 ! CALL collocate_general() 585 ELSE 586 CALL collocate_general_wings() 587 !CALL collocate_general_opt() 588 END IF 589 END IF 590 591 CONTAINS 592 593! ************************************************************************************************** 594!> \brief this treats efficiently the orthogonal case 595! ************************************************************************************************** 596 SUBROUTINE collocate_ortho() 597 598! *** properties of the grid *** 599 600 ! notice we're in the ortho case 601 dr(1) = rsgrid%desc%dh(1, 1) 602 dr(2) = rsgrid%desc%dh(2, 2) 603 dr(3) = rsgrid%desc%dh(3, 3) 604 605! *** get the sub grid properties for the given radius *** 606 CALL return_cube(cube_info, radius, lb_cube, ub_cube, sphere_bounds) 607 cmax = MAXVAL(ub_cube) 608 609! *** position of the gaussian product 610! 611! this is the actual definition of the position on the grid 612! i.e. a point rp(:) gets here grid coordinates 613! MODULO(rp(:)/dr(:),ng(:))+1 614! hence (0.0,0.0,0.0) in real space is rsgrid%lb on the rsgrid ((1,1,1) on grid) 615! 616 617 ALLOCATE (map(-cmax:cmax, 3)) 618 CALL compute_cube_center(cubecenter, rsgrid%desc, zeta, zetb, ra, rab) 619 roffset(:) = rp(:) - REAL(cubecenter(:), dp)*dr(:) 620! *** a mapping so that the ig corresponds to the right grid point 621 DO i = 1, 3 622 IF (rsgrid%desc%perd(i) == 1) THEN 623 start = lb_cube(i) 624 DO 625 offset = MODULO(cubecenter(i) + start, ng(i)) + 1 - start 626 length = MIN(ub_cube(i), ng(i) - offset) - start 627 DO ig = start, start + length 628 map(ig, i) = ig + offset 629 END DO 630 IF (start + length .GE. ub_cube(i)) EXIT 631 start = start + length + 1 632 END DO 633 ELSE 634 ! this takes partial grid + border regions into account 635 offset = MODULO(cubecenter(i) + lb_cube(i) + rsgrid%desc%lb(i) - rsgrid%lb_local(i), ng(i)) + 1 - lb_cube(i) 636 ! check for out of bounds 637 IF (ub_cube(i) + offset > UBOUND(grid, i) .OR. lb_cube(i) + offset < LBOUND(grid, i)) THEN 638 CPASSERT(.FALSE.) 639 ENDIF 640 DO ig = lb_cube(i), ub_cube(i) 641 map(ig, i) = ig + offset 642 END DO 643 END IF 644 ENDDO 645 ALLOCATE (pol_z(1:2, 0:lp, -cmax:0)) 646 ALLOCATE (pol_y(1:2, 0:lp, -cmax:0)) 647 ALLOCATE (pol_x(0:lp, -cmax:cmax)) 648 649#include "prep.f90" 650#include "call_collocate.f90" 651 652 ! deallocation needed to pass around a pgi bug.. 653 DEALLOCATE (pol_z) 654 DEALLOCATE (pol_y) 655 DEALLOCATE (pol_x) 656 DEALLOCATE (map) 657 658 END SUBROUTINE collocate_ortho 659 660! ************************************************************************************************** 661!> \brief this is a general 'optimized' routine to do the collocation 662! ************************************************************************************************** 663 SUBROUTINE collocate_general_opt() 664 665 INTEGER :: i, i_index, il, ilx, ily, ilz, index_max(3), index_min(3), ismax, ismin, j, & 666 j_index, jl, jlx, jly, jlz, k, k_index, kl, klx, kly, klz, lpx, lpy, lpz, lx, ly, lz, & 667 offset(3) 668 INTEGER, ALLOCATABLE, DIMENSION(:) :: grid_map 669 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: coef_map 670 REAL(KIND=dp) :: a, b, c, d, di, dip, dj, djp, dk, dkp, & 671 exp0i, exp1i, exp2i, gp(3), & 672 hmatgrid(3, 3), pointj(3), pointk(3), & 673 res, v(3) 674 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coef_ijk 675 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: hmatgridp 676 677! 678! transform P_{lxp,lyp,lzp} into a P_{lip,ljp,lkp} such that 679! sum_{lxp,lyp,lzp} P_{lxp,lyp,lzp} (x-x_p)**lxp (y-y_p)**lyp (z-z_p)**lzp = 680! sum_{lip,ljp,lkp} P_{lip,ljp,lkp} (i-i_p)**lip (j-j_p)**ljp (k-k_p)**lkp 681! 682 683 ALLOCATE (coef_ijk(((lp + 1)*(lp + 2)*(lp + 3))/6)) 684 685 ! aux mapping array to simplify life 686 ALLOCATE (coef_map(0:lp, 0:lp, 0:lp)) 687 coef_map = HUGE(coef_map) 688 lxyz = 0 689 DO lzp = 0, lp 690 DO lyp = 0, lp - lzp 691 DO lxp = 0, lp - lzp - lyp 692 lxyz = lxyz + 1 693 coef_ijk(lxyz) = 0.0_dp 694 coef_map(lxp, lyp, lzp) = lxyz 695 ENDDO 696 ENDDO 697 ENDDO 698 699 ! cell hmat in grid points 700 hmatgrid = rsgrid%desc%dh 701 702 ! center in grid coords 703 gp = MATMUL(rsgrid%desc%dh_inv, rp) 704 cubecenter(:) = FLOOR(gp) 705 706 ! transform using multinomials 707 ALLOCATE (hmatgridp(3, 3, 0:lp)) 708 hmatgridp(:, :, 0) = 1.0_dp 709 DO k = 1, lp 710 hmatgridp(:, :, k) = hmatgridp(:, :, k - 1)*hmatgrid(:, :) 711 ENDDO 712 713 lpx = lp 714 DO klx = 0, lpx 715 DO jlx = 0, lpx - klx 716 DO ilx = 0, lpx - klx - jlx 717 lx = ilx + jlx + klx 718 lpy = lp - lx 719 DO kly = 0, lpy 720 DO jly = 0, lpy - kly 721 DO ily = 0, lpy - kly - jly 722 ly = ily + jly + kly 723 lpz = lp - lx - ly 724 DO klz = 0, lpz 725 DO jlz = 0, lpz - klz 726 DO ilz = 0, lpz - klz - jlz 727 lz = ilz + jlz + klz 728 729 il = ilx + ily + ilz 730 jl = jlx + jly + jlz 731 kl = klx + kly + klz 732 coef_ijk(coef_map(il, jl, kl)) = & 733 coef_ijk(coef_map(il, jl, kl)) + coef_xyz(coef_map(lx, ly, lz))* & 734 hmatgridp(1, 1, ilx)*hmatgridp(1, 2, jlx)*hmatgridp(1, 3, klx)* & 735 hmatgridp(2, 1, ily)*hmatgridp(2, 2, jly)*hmatgridp(2, 3, kly)* & 736 hmatgridp(3, 1, ilz)*hmatgridp(3, 2, jlz)*hmatgridp(3, 3, klz)* & 737 fac(lx)*fac(ly)*fac(lz)/ & 738 (fac(ilx)*fac(ily)*fac(ilz)*fac(jlx)*fac(jly)*fac(jlz)*fac(klx)*fac(kly)*fac(klz)) 739 ENDDO 740 ENDDO 741 ENDDO 742 ENDDO 743 ENDDO 744 ENDDO 745 ENDDO 746 ENDDO 747 ENDDO 748 749 CALL return_cube_nonortho(cube_info, radius, index_min, index_max, rp) 750 751 offset(:) = MODULO(index_min(:) + rsgrid%desc%lb(:) - rsgrid%lb_local(:), ng(:)) + 1 752 753 ALLOCATE (grid_map(index_min(1):index_max(1))) 754 DO i = index_min(1), index_max(1) 755 grid_map(i) = MODULO(i, ng(1)) + 1 756 IF (rsgrid%desc%perd(1) == 1) THEN 757 grid_map(i) = MODULO(i, ng(1)) + 1 758 ELSE 759 grid_map(i) = i - index_min(1) + offset(1) 760 ENDIF 761 ENDDO 762 763 ! go over the grid, but cycle if the point is not within the radius 764 DO k = index_min(3), index_max(3) 765 dk = k - gp(3) 766 pointk = hmatgrid(:, 3)*dk 767 768 IF (rsgrid%desc%perd(3) == 1) THEN 769 k_index = MODULO(k, ng(3)) + 1 770 ELSE 771 k_index = k - index_min(3) + offset(3) 772 ENDIF 773 774 coef_xyt = 0.0_dp 775 lxyz = 0 776 dkp = 1.0_dp 777 DO kl = 0, lp 778 lxy = 0 779 DO jl = 0, lp - kl 780 DO il = 0, lp - kl - jl 781 lxyz = lxyz + 1; lxy = lxy + 1 782 coef_xyt(lxy) = coef_xyt(lxy) + coef_ijk(lxyz)*dkp 783 ENDDO 784 lxy = lxy + kl 785 ENDDO 786 dkp = dkp*dk 787 ENDDO 788 789 DO j = index_min(2), index_max(2) 790 dj = j - gp(2) 791 pointj = pointk + hmatgrid(:, 2)*dj 792 IF (rsgrid%desc%perd(2) == 1) THEN 793 j_index = MODULO(j, ng(2)) + 1 794 ELSE 795 j_index = j - index_min(2) + offset(2) 796 ENDIF 797 798 coef_xtt = 0.0_dp 799 lxy = 0 800 djp = 1.0_dp 801 DO jl = 0, lp 802 DO il = 0, lp - jl 803 lxy = lxy + 1 804 coef_xtt(il) = coef_xtt(il) + coef_xyt(lxy)*djp 805 ENDDO 806 djp = djp*dj 807 ENDDO 808 809 ! find bounds for the inner loop 810 ! based on a quadratic equation in i 811 ! a*i**2+b*i+c=radius**2 812 v = pointj - gp(1)*hmatgrid(:, 1) 813 a = DOT_PRODUCT(hmatgrid(:, 1), hmatgrid(:, 1)) 814 b = 2*DOT_PRODUCT(v, hmatgrid(:, 1)) 815 c = DOT_PRODUCT(v, v) 816 d = b*b - 4*a*(c - radius**2) 817 818 IF (d < 0) THEN 819 CYCLE 820 ELSE 821 d = SQRT(d) 822 ismin = CEILING((-b - d)/(2*a)) 823 ismax = FLOOR((-b + d)/(2*a)) 824 ENDIF 825 ! prepare for computing -zetp*rsq 826 a = -zetp*a 827 b = -zetp*b 828 c = -zetp*c 829 i = ismin - 1 830 831 ! the recursion relation might have to be done 832 ! from the center of the gaussian (in both directions) 833 ! instead as the current implementation from an edge 834 exp2i = EXP((a*i + b)*i + c) 835 exp1i = EXP(2*a*i + a + b) 836 exp0i = EXP(2*a) 837 838 DO i = ismin, ismax 839 di = i - gp(1) 840 841 ! polynomial terms 842 res = 0.0_dp 843 dip = 1.0_dp 844 DO il = 0, lp 845 res = res + coef_xtt(il)*dip 846 dip = dip*di 847 ENDDO 848 849 ! the exponential recursion 850 exp2i = exp2i*exp1i 851 exp1i = exp1i*exp0i 852 res = res*exp2i 853 854 i_index = grid_map(i) 855 grid(i_index, j_index, k_index) = grid(i_index, j_index, k_index) + res 856 ENDDO 857 ENDDO 858 ENDDO 859 !t2=nanotime_ia32() 860 !write(*,*) t2-t1 861 ! deallocation needed to pass around a pgi bug.. 862 DEALLOCATE (coef_ijk) 863 DEALLOCATE (coef_map) 864 DEALLOCATE (hmatgridp) 865 DEALLOCATE (grid_map) 866 867 END SUBROUTINE collocate_general_opt 868 869! ************************************************************************************************** 870!> \brief ... 871! ************************************************************************************************** 872 SUBROUTINE collocate_general_subpatch() 873 INTEGER, DIMENSION(2, 3) :: local_b 874 INTEGER, DIMENSION(3) :: local_s, periodic 875 REAL(dp), DIMENSION((lp+1)*(lp+2)*(lp+3)/6) :: poly_d3 876 877 periodic = 1 ! cell%perd 878 CALL poly_cp2k2d3(coef_xyz, lp, poly_d3) 879 local_b(1, :) = rsgrid%lb_real - rsgrid%desc%lb 880 local_b(2, :) = rsgrid%ub_real - rsgrid%desc%lb 881 local_s = rsgrid%lb_real - rsgrid%lb_local 882 IF (BTEST(subpatch_pattern, 0)) local_b(1, 1) = local_b(1, 1) - rsgrid%desc%border 883 IF (BTEST(subpatch_pattern, 1)) local_b(2, 1) = local_b(2, 1) + rsgrid%desc%border 884 IF (BTEST(subpatch_pattern, 2)) local_b(1, 2) = local_b(1, 2) - rsgrid%desc%border 885 IF (BTEST(subpatch_pattern, 3)) local_b(2, 2) = local_b(2, 2) + rsgrid%desc%border 886 IF (BTEST(subpatch_pattern, 4)) local_b(1, 3) = local_b(1, 3) - rsgrid%desc%border 887 IF (BTEST(subpatch_pattern, 5)) local_b(2, 3) = local_b(2, 3) + rsgrid%desc%border 888 IF (BTEST(subpatch_pattern, 0)) local_s(1) = local_s(1) - rsgrid%desc%border 889 IF (BTEST(subpatch_pattern, 2)) local_s(2) = local_s(2) - rsgrid%desc%border 890 IF (BTEST(subpatch_pattern, 4)) local_s(3) = local_s(3) - rsgrid%desc%border 891 CALL collocGauss(h=cell%hmat, h_inv=cell%h_inv, & 892 grid=grid, poly=poly_d3, alphai=zetp, posi=rp, max_r2=radius*radius, & 893 periodic=periodic, gdim=ng, local_bounds=local_b, local_shift=local_s) 894 ! defaults: local_shift=(/0,0,0/),poly_shift=(/0.0_dp,0.0_dp,0.0_dp/),scale=1.0_dp, 895 896 END SUBROUTINE 897 898! ************************************************************************************************** 899!> \brief ... 900! ************************************************************************************************** 901 SUBROUTINE collocate_general_wings() 902 INTEGER, DIMENSION(2, 3) :: local_b 903 INTEGER, DIMENSION(3) :: periodic 904 REAL(dp), DIMENSION((lp+1)*(lp+2)*(lp+3)/6) :: poly_d3 905 REAL(dp), DIMENSION(3) :: local_shift, rShifted 906 907 periodic = 1 ! cell%perd 908 CALL poly_cp2k2d3(coef_xyz, lp, poly_d3) 909 local_b(1, :) = 0 910 local_b(2, :) = MIN(rsgrid%desc%npts - 1, rsgrid%ub_local - rsgrid%lb_local) 911 local_shift = REAL(rsgrid%desc%lb - rsgrid%lb_local, dp)/REAL(rsgrid%desc%npts, dp) 912 rShifted(1) = rp(1) + cell%hmat(1, 1)*local_shift(1) & 913 + cell%hmat(1, 2)*local_shift(2) & 914 + cell%hmat(1, 3)*local_shift(3) 915 rShifted(2) = rp(2) + cell%hmat(2, 1)*local_shift(1) & 916 + cell%hmat(2, 2)*local_shift(2) & 917 + cell%hmat(2, 3)*local_shift(3) 918 rShifted(3) = rp(3) + cell%hmat(3, 1)*local_shift(1) & 919 + cell%hmat(3, 2)*local_shift(2) & 920 + cell%hmat(3, 3)*local_shift(3) 921 CALL collocGauss(h=cell%hmat, h_inv=cell%h_inv, & 922 grid=grid, poly=poly_d3, alphai=zetp, posi=rShifted, max_r2=radius*radius, & 923 periodic=periodic, gdim=ng, local_bounds=local_b) 924 ! defaults: local_shift=(/0,0,0/),poly_shift=(/0.0_dp,0.0_dp,0.0_dp/),scale=1.0_dp, 925 926 END SUBROUTINE 927 928! ************************************************************************************************** 929!> \brief this is a general 'reference' routine to do the collocation 930! ************************************************************************************************** 931 SUBROUTINE collocate_general() 932 933 INTEGER :: i, index_max(3), index_min(3), & 934 ipoint(3), j, k 935 REAL(KIND=dp) :: inv_ng(3), point(3), primpt 936 937! still hard-wired (see MODULO) 938 939 CPASSERT(ALL(rsgrid%desc%perd == 1)) 940 941 CALL return_cube_nonortho(cube_info, radius, index_min, index_max, rp) 942 943 inv_ng = 1.0_dp/ng 944 945 ! go over the grid, but cycle if the point is not within the radius 946 DO k = index_min(3), index_max(3) 947 DO j = index_min(2), index_max(2) 948 DO i = index_min(1), index_max(1) 949 ! point in real space 950 point = MATMUL(cell%hmat, REAL((/i, j, k/), KIND=dp)*inv_ng) 951 ! primitive_value of point 952 primpt = primitive_value(point) 953 ! skip if outside of the sphere 954 IF (SUM((point - rp)**2) > radius**2) CYCLE 955 ! point on the grid (including pbc) 956 ipoint = MODULO((/i, j, k/), ng) + 1 957 ! add to grid 958 grid(ipoint(1), ipoint(2), ipoint(3)) = grid(ipoint(1), ipoint(2), ipoint(3)) + primpt 959 ENDDO 960 ENDDO 961 ENDDO 962 963 END SUBROUTINE collocate_general 964 965! ************************************************************************************************** 966!> \brief ... 967!> \param point ... 968!> \return ... 969! ************************************************************************************************** 970 FUNCTION primitive_value(point) RESULT(res) 971 REAL(KIND=dp) :: point(3), res 972 973 REAL(KIND=dp) :: dra(3), drap(3), drb(3), drbp(3), myexp, & 974 pdrap 975 976 res = 0.0_dp 977 myexp = EXP(-zetp*SUM((point - rp)**2))*prefactor 978 dra = point - ra 979 drb = point - rb 980 drap(1) = 1.0_dp 981 DO lxa = 0, la_max_local 982 drbp(1) = 1.0_dp 983 DO lxb = 0, lb_max_local 984 drap(2) = 1.0_dp 985 DO lya = 0, la_max_local - lxa 986 drbp(2) = 1.0_dp 987 DO lyb = 0, lb_max_local - lxb 988 drap(3) = 1.0_dp 989 DO lza = 1, MAX(la_min_local - lxa - lya, 0) 990 drap(3) = drap(3)*dra(3) 991 ENDDO 992 DO lza = MAX(la_min_local - lxa - lya, 0), la_max_local - lxa - lya 993 drbp(3) = 1.0_dp 994 DO lzb = 1, MAX(lb_min_local - lxb - lyb, 0) 995 drbp(3) = drbp(3)*drb(3) 996 ENDDO 997 ico = coset(lxa, lya, lza) 998 pdrap = PRODUCT(drap) 999 DO lzb = MAX(lb_min_local - lxb - lyb, 0), lb_max_local - lxb - lyb 1000 jco = coset(lxb, lyb, lzb) 1001 res = res + pab_local(ico + o1_local, jco + o2_local)*myexp*pdrap*PRODUCT(drbp) 1002 drbp(3) = drbp(3)*drb(3) 1003 ENDDO 1004 drap(3) = drap(3)*dra(3) 1005 ENDDO 1006 drbp(2) = drbp(2)*drb(2) 1007 ENDDO 1008 drap(2) = drap(2)*dra(2) 1009 ENDDO 1010 drbp(1) = drbp(1)*drb(1) 1011 ENDDO 1012 drap(1) = drap(1)*dra(1) 1013 ENDDO 1014 1015 END FUNCTION primitive_value 1016 1017 END SUBROUTINE collocate_pgf_product_part2 1018 1019END MODULE grid_collocate 1020