1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Utilities for Integrals for semi-empiric methods 8!> \author Teodoro Laino (03.2008) [tlaino] - University of Zurich 9! ************************************************************************************************** 10MODULE semi_empirical_int_utils 11 12 USE input_constants, ONLY: do_method_pchg,& 13 do_se_IS_kdso_d 14 USE kinds, ONLY: dp 15 USE semi_empirical_int3_utils, ONLY: charg_int_3,& 16 dcharg_int_3,& 17 ijkl_low_3 18 USE semi_empirical_int_arrays, ONLY: & 19 CLMp, CLMxx, CLMxy, CLMyy, CLMz, CLMzp, CLMzz, clm_d, clm_sp, ijkl_ind, indexa, indexb, & 20 int2c_type 21 USE semi_empirical_types, ONLY: rotmat_type,& 22 se_int_control_type,& 23 se_int_screen_type,& 24 se_taper_type,& 25 semi_empirical_type 26#include "./base/base_uses.f90" 27 28 IMPLICIT NONE 29#include "semi_empirical_int_debug.h" 30 31 PRIVATE 32 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE. 33 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_utils' 34 35 PUBLIC :: ijkl_sp, ijkl_d, rotmat, rot_2el_2c_first, store_2el_2c_diag, & 36 d_ijkl_sp, d_ijkl_d 37 38CONTAINS 39 40! ************************************************************************************************** 41!> \brief General driver for computing semi-empirical integrals <ij|kl> with 42!> sp basis set. This code uses the old definitions of quadrupoles and 43!> therefore cannot be used for integrals involving d-orbitals (which 44!> require a definition of quadrupoles based on the rotational invariant 45!> property) 46!> 47!> \param sepi ... 48!> \param sepj ... 49!> \param ij ... 50!> \param kl ... 51!> \param li ... 52!> \param lj ... 53!> \param lk ... 54!> \param ll ... 55!> \param ic ... 56!> \param r ... 57!> \param se_int_control ... 58!> \param se_int_screen ... 59!> \param itype ... 60!> \return ... 61!> \date 04.2008 [tlaino] 62!> \author Teodoro Laino [tlaino] - University of Zurich 63! ************************************************************************************************** 64 FUNCTION ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control, & 65 se_int_screen, itype) RESULT(res) 66 TYPE(semi_empirical_type), POINTER :: sepi, sepj 67 INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic 68 REAL(KIND=dp), INTENT(IN) :: r 69 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 70 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 71 INTEGER, INTENT(IN) :: itype 72 REAL(KIND=dp) :: res 73 74 CHARACTER(len=*), PARAMETER :: routineN = 'ijkl_sp', routineP = moduleN//':'//routineN 75 76 res = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, & 77 se_int_control%integral_screening, se_int_control%shortrange, & 78 se_int_control%pc_coulomb_int, se_int_control%max_multipole, & 79 itype, charg_int_nri) 80 81 ! If only the shortrange component is requested we can skip the rest 82 IF ((.NOT. se_int_control%pc_coulomb_int) .AND. (itype /= do_method_pchg)) THEN 83 ! Handle the 1/r^3 term, this term is ALWAYS false for KDSO-D integrals 84 IF (se_int_control%shortrange .AND. se_int_control%do_ewald_r3) THEN 85 res = res - ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, & 86 itype, charg_int_3) 87 END IF 88 END IF 89 END FUNCTION ijkl_sp 90 91! ************************************************************************************************** 92!> \brief General driver for computing derivatives of semi-empirical integrals 93!> <ij|kl> with sp basis set. 94!> This code uses the old definitions of quadrupoles and therefore 95!> cannot be used for integrals involving d-orbitals (which requires a 96!> definition of quadrupoles based on the rotational invariant property) 97!> 98!> \param sepi ... 99!> \param sepj ... 100!> \param ij ... 101!> \param kl ... 102!> \param li ... 103!> \param lj ... 104!> \param lk ... 105!> \param ll ... 106!> \param ic ... 107!> \param r ... 108!> \param se_int_control ... 109!> \param se_int_screen ... 110!> \param itype ... 111!> \return ... 112!> \date 05.2008 [tlaino] 113!> \author Teodoro Laino [tlaino] - University of Zurich 114! ************************************************************************************************** 115 FUNCTION d_ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control, & 116 se_int_screen, itype) RESULT(res) 117 TYPE(semi_empirical_type), POINTER :: sepi, sepj 118 INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic 119 REAL(KIND=dp), INTENT(IN) :: r 120 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 121 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 122 INTEGER, INTENT(IN) :: itype 123 REAL(KIND=dp) :: res 124 125 CHARACTER(len=*), PARAMETER :: routineN = 'd_ijkl_sp', routineP = moduleN//':'//routineN 126 127 REAL(KIND=dp) :: dfs, srd 128 129 IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN 130 res = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, & 131 se_int_control%integral_screening, .FALSE., & 132 se_int_control%pc_coulomb_int, se_int_control%max_multipole, & 133 itype, dcharg_int_nri) 134 135 IF (.NOT. se_int_control%pc_coulomb_int) THEN 136 dfs = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, & 137 se_int_control%integral_screening, .FALSE., .FALSE., & 138 se_int_control%max_multipole, itype, dcharg_int_nri_fs) 139 res = res + dfs*se_int_screen%dft 140 141 ! In case we need the shortrange part we have to evaluate an additional derivative 142 ! to handle the derivative of the Tapering term 143 IF (se_int_control%shortrange) THEN 144 srd = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, & 145 se_int_control%integral_screening, .FALSE., .TRUE., & 146 se_int_control%max_multipole, itype, dcharg_int_nri) 147 res = res - srd 148 END IF 149 END IF 150 ELSE 151 res = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, & 152 se_int_control%integral_screening, se_int_control%shortrange, & 153 se_int_control%pc_coulomb_int, se_int_control%max_multipole, & 154 itype, dcharg_int_nri) 155 END IF 156 157 ! If only the shortrange component is requested we can skip the rest 158 IF ((.NOT. se_int_control%pc_coulomb_int) .AND. (itype /= do_method_pchg)) THEN 159 ! Handle the 1/r^3 term, this term is ALWAYS false for KDSO-D integrals 160 IF (se_int_control%shortrange .AND. se_int_control%do_ewald_r3) THEN 161 res = res - ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, & 162 itype, dcharg_int_3) 163 END IF 164 END IF 165 166 END FUNCTION d_ijkl_sp 167 168! ************************************************************************************************** 169!> \brief Low level general driver for computing semi-empirical integrals 170!> <ij|kl> and their derivatives with sp basis set only. 171!> This code uses the old definitions of quadrupoles and 172!> therefore cannot be used for integrals involving d-orbitals (which 173!> require a definition of quadrupoles based on the rotational invariant 174!> property) 175!> 176!> \param sepi ... 177!> \param sepj ... 178!> \param ij ... 179!> \param kl ... 180!> \param li ... 181!> \param lj ... 182!> \param lk ... 183!> \param ll ... 184!> \param ic ... 185!> \param r ... 186!> \param se_int_screen ... 187!> \param iscreen ... 188!> \param shortrange ... 189!> \param pc_coulomb_int ... 190!> \param max_multipole ... 191!> \param itype ... 192!> \param eval a function without explicit interface 193!> \return ... 194!> \date 05.2008 [tlaino] 195!> \author Teodoro Laino [tlaino] - University of Zurich 196! ************************************************************************************************** 197 FUNCTION ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, & 198 iscreen, shortrange, pc_coulomb_int, max_multipole, itype, eval) RESULT(res) 199 TYPE(semi_empirical_type), POINTER :: sepi, sepj 200 INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic 201 REAL(KIND=dp), INTENT(IN) :: r 202 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 203 INTEGER, INTENT(IN) :: iscreen 204 LOGICAL, INTENT(IN) :: shortrange, pc_coulomb_int 205 INTEGER, INTENT(IN) :: max_multipole, itype 206 REAL(KIND=dp) :: eval, res 207 208 CHARACTER(len=*), PARAMETER :: routineN = 'ijkl_sp_low', routineP = moduleN//':'//routineN 209 210 INTEGER :: ccc, l1, l1max, l1min, l2, l2max, l2min, & 211 lij, lkl, lmin, m 212 REAL(KIND=dp) :: add, chrg, dij, dkl, fact_ij, fact_kl, & 213 fact_screen, pij, pkl, s1, sum 214 215 l1min = ABS(li - lj) 216 l1max = li + lj 217 lij = indexb(li + 1, lj + 1) 218 l2min = ABS(lk - ll) 219 l2max = lk + ll 220 lkl = indexb(lk + 1, ll + 1) 221 222 l1max = MIN(l1max, 2) 223 l1min = MIN(l1min, 2) 224 l2max = MIN(l2max, 2) 225 l2min = MIN(l2min, 2) 226 sum = 0.0_dp 227 dij = 0.0_dp 228 dkl = 0.0_dp 229 fact_ij = 1.0_dp 230 fact_kl = 1.0_dp 231 fact_screen = 1.0_dp 232 IF (lij == 3) fact_ij = SQRT(2.0_dp) 233 IF (lkl == 3) fact_kl = SQRT(2.0_dp) 234 IF (.NOT. pc_coulomb_int) THEN 235 IF (iscreen == do_se_IS_kdso_d) fact_screen = se_int_screen%ft 236 ! Standard value of the integral 237 DO l1 = l1min, l1max 238 IF (l1 == 0) THEN 239 IF (lij == 1) THEN 240 pij = sepi%ko(1) 241 IF (ic == -1 .OR. ic == 1) THEN 242 pij = sepi%ko(9) 243 END IF 244 ELSE IF (lij == 3) THEN 245 pij = sepi%ko(7) 246 END IF 247 ELSE 248 dij = sepi%cs(lij)*fact_ij 249 pij = sepi%ko(lij) 250 END IF 251 ! 252 DO l2 = l2min, l2max 253 IF (l2 == 0) THEN 254 IF (lkl == 1) THEN 255 pkl = sepj%ko(1) 256 IF (ic == -1 .OR. ic == 2) THEN 257 pkl = sepj%ko(9) 258 END IF 259 ELSE IF (lkl == 3) THEN 260 pkl = sepj%ko(7) 261 END IF 262 ELSE 263 dkl = sepj%cs(lkl)*fact_kl 264 pkl = sepj%ko(lkl) 265 END IF 266 IF (itype == do_method_pchg) THEN 267 add = 0.0_dp 268 ELSE 269 add = (pij + pkl)**2 270 END IF 271 lmin = MAX(l1, l2) 272 s1 = 0.0_dp 273 DO m = -lmin, lmin 274 ccc = clm_sp(ij, l1, m)*clm_sp(kl, l2, m) 275 IF (ABS(ccc) > EPSILON(0.0_dp)) THEN 276 chrg = eval(r, l1, l2, clm_sp(ij, l1, m), clm_sp(kl, l2, m), dij, dkl, add, fact_screen) 277 s1 = s1 + chrg 278 END IF 279 END DO 280 sum = sum + s1 281 END DO 282 END DO 283 res = sum 284 END IF 285 ! Shortrange: Possibly computes pure Coulomb and subtract from the original integral value 286 IF (shortrange .OR. pc_coulomb_int) THEN 287 sum = 0.0_dp 288 dij = 0.0_dp 289 dkl = 0.0_dp 290 add = 0.0_dp 291 fact_screen = 0.0_dp 292 DO l1 = l1min, l1max 293 IF (l1 > max_multipole) CYCLE 294 IF (l1 /= 0) THEN 295 dij = sepi%cs(lij)*fact_ij 296 END IF 297 ! 298 DO l2 = l2min, l2max 299 IF (l2 > max_multipole) CYCLE 300 IF (l2 /= 0) THEN 301 dkl = sepj%cs(lkl)*fact_kl 302 END IF 303 lmin = MAX(l1, l2) 304 s1 = 0.0_dp 305 DO m = -lmin, lmin 306 ccc = clm_sp(ij, l1, m)*clm_sp(kl, l2, m) 307 IF (ABS(ccc) > EPSILON(0.0_dp)) THEN 308 chrg = eval(r, l1, l2, clm_sp(ij, l1, m), clm_sp(kl, l2, m), dij, dkl, add, fact_screen) 309 s1 = s1 + chrg 310 END IF 311 END DO 312 sum = sum + s1 313 END DO 314 END DO 315 IF (pc_coulomb_int) res = sum 316 IF (shortrange) res = res - sum 317 END IF 318 END FUNCTION ijkl_sp_low 319 320! ************************************************************************************************** 321!> \brief Interaction function between two point-charge configurations NDDO sp-code 322!> Non-Rotational Invariant definition of quadrupoles 323!> r - Distance r12 324!> l1,m - Quantum numbers for multipole of configuration 1 325!> l2,m - Quantum numbers for multipole of configuration 2 326!> da - charge separation of configuration 1 327!> db - charge separation of configuration 2 328!> add - additive term 329!> 330!> \param r ... 331!> \param l1_i ... 332!> \param l2_i ... 333!> \param m1_i ... 334!> \param m2_i ... 335!> \param da_i ... 336!> \param db_i ... 337!> \param add0 ... 338!> \param fact_screen ... 339!> \return ... 340!> \date 04.2008 [tlaino] 341!> \author Teodoro Laino [tlaino] - University of Zurich 342! ************************************************************************************************** 343 FUNCTION charg_int_nri(r, l1_i, l2_i, m1_i, m2_i, da_i, db_i, add0, fact_screen) RESULT(charg) 344 REAL(KIND=dp), INTENT(in) :: r 345 INTEGER, INTENT(in) :: l1_i, l2_i, m1_i, m2_i 346 REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen 347 REAL(KIND=dp) :: charg 348 349 CHARACTER(len=*), PARAMETER :: routineN = 'charg_int_nri', routineP = moduleN//':'//routineN 350 351 INTEGER :: l1, l2, m1, m2 352 REAL(KIND=dp) :: add, da, db, dxdx, dxqxz, dzdz, dzqxx, & 353 dzqzz, fact, qqxx, qqzz, qxxqxx, & 354 qxxqyy, qxzqxz, xyxy, zzzz 355 356! Computing only Integral Values 357 358 IF (l1_i < l2_i) THEN 359 l1 = l1_i 360 l2 = l2_i 361 m1 = m1_i 362 m2 = m2_i 363 da = da_i 364 db = db_i 365 fact = 1.0_dp 366 ELSE IF (l1_i > l2_i) THEN 367 l1 = l2_i 368 l2 = l1_i 369 m1 = m2_i 370 m2 = m1_i 371 da = db_i 372 db = da_i 373 fact = (-1.0_dp)**(l1 + l2) 374 ELSE IF (l1_i == l2_i) THEN 375 l1 = l1_i 376 l2 = l2_i 377 IF (m1_i <= m2_i) THEN 378 m1 = m1_i 379 m2 = m2_i 380 da = da_i 381 db = db_i 382 ELSE 383 m1 = m2_i 384 m2 = m1_i 385 da = db_i 386 db = da_i 387 END IF 388 fact = 1.0_dp 389 END IF 390 add = add0*fact_screen 391 charg = 0.0_dp 392 ! Q - Q. 393 IF (l1 == 0 .AND. l2 == 0) THEN 394 charg = fact/SQRT(r**2 + add) 395 RETURN 396 END IF 397 ! Q - Z. 398 IF (l1 == 0 .AND. l2 == 1 .AND. m2 == CLMz) THEN 399 charg = 1.0_dp/SQRT((r + db)**2 + add) - 1.0_dp/SQRT((r - db)**2 + add) 400 charg = charg*0.5_dp*fact 401 RETURN 402 END IF 403 ! Z - Z. 404 IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMz .AND. m2 == CLMz) THEN 405 dzdz = & 406 +1.0_dp/SQRT((r + da - db)**2 + add) + 1.0_dp/SQRT((r - da + db)**2 + add) & 407 - 1.0_dp/SQRT((r - da - db)**2 + add) - 1.0_dp/SQRT((r + da + db)**2 + add) 408 charg = dzdz*0.25_dp*fact 409 RETURN 410 END IF 411 ! X - X 412 IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMp .AND. m2 == CLMp) THEN 413 dxdx = 2.0_dp/SQRT(r**2 + (da - db)**2 + add) - 2.0_dp/SQRT(r**2 + (da + db)**2 + add) 414 charg = dxdx*0.25_dp*fact 415 RETURN 416 END IF 417 ! Q - ZZ 418 IF (l1 == 0 .AND. l2 == 2 .AND. m2 == CLMzz) THEN 419 qqzz = 1.0_dp/SQRT((r - db)**2 + add) - 2.0_dp/SQRT(r**2 + add) + 1.0_dp/SQRT((r + db)**2 + add) 420 charg = qqzz*0.25_dp*fact 421 RETURN 422 END IF 423 ! Q - XX 424 IF (l1 == 0 .AND. l2 == 2 .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN 425 qqxx = -1.0_dp/SQRT(r**2 + add) + 1.0_dp/SQRT(r**2 + add + db**2) 426 charg = qqxx*0.5_dp*fact 427 RETURN 428 END IF 429 ! Z - ZZ 430 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. m2 == CLMzz) THEN 431 dzqzz = & 432 +1.0_dp/SQRT((r - da - db)**2 + add) - 2.0_dp/SQRT((r - da)**2 + add) & 433 + 1.0_dp/SQRT((r - da + db)**2 + add) - 1.0_dp/SQRT((r + da - db)**2 + add) & 434 + 2.0_dp/SQRT((r + da)**2 + add) - 1.0_dp/SQRT((r + da + db)**2 + add) 435 charg = dzqzz*0.125_dp*fact 436 RETURN 437 END IF 438 ! Z - XX 439 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN 440 dzqxx = & 441 +1.0_dp/SQRT((r + da)**2 + add) - 1.0_dp/SQRT((r + da)**2 + add + db**2) & 442 - 1.0_dp/SQRT((r - da)**2 + add) + 1.0_dp/SQRT((r - da)**2 + add + db**2) 443 charg = dzqxx*0.25_dp*fact 444 RETURN 445 END IF 446 ! ZZ - ZZ 447 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. m2 == CLMzz) THEN 448 zzzz = & 449 +1.0_dp/SQRT((r - da - db)**2 + add) + 1.0_dp/SQRT((r + da + db)**2 + add) & 450 + 1.0_dp/SQRT((r - da + db)**2 + add) + 1.0_dp/SQRT((r + da - db)**2 + add) 451 xyxy = & 452 +1.0_dp/SQRT((r - da)**2 + add) + 1.0_dp/SQRT((r + da)**2 + add) & 453 + 1.0_dp/SQRT((r - db)**2 + add) + 1.0_dp/SQRT((r + db)**2 + add) & 454 - 2.0_dp/SQRT(r**2 + add) 455 charg = (zzzz*0.0625_dp - xyxy*0.125_dp)*fact 456 RETURN 457 END IF 458 ! ZZ - XX 459 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. (m2 == CLMxx .OR. m2 == CLMyy)) THEN 460 zzzz = & 461 -1.0_dp/SQRT((r + da)**2 + add) + 1.0_dp/SQRT((r + da)**2 + db**2 + add) & 462 - 1.0_dp/SQRT((r - da)**2 + add) + 1.0_dp/SQRT((r - da)**2 + db**2 + add) 463 xyxy = & 464 +1.0_dp/SQRT(r**2 + db**2 + add) - 1.0_dp/SQRT(r**2 + add) 465 charg = (zzzz*0.125_dp - xyxy*0.25_dp)*fact 466 RETURN 467 END IF 468 ! X - ZX 469 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMp .AND. m2 == CLMzp) THEN 470 db = db/2.0_dp 471 dxqxz = & 472 -1.0_dp/SQRT((r - db)**2 + (da - db)**2 + add) + 1.0_dp/SQRT((r + db)**2 + (da - db)**2 + add) & 473 + 1.0_dp/SQRT((r - db)**2 + (da + db)**2 + add) - 1.0_dp/SQRT((r + db)**2 + (da + db)**2 + add) 474 charg = dxqxz*0.25_dp*fact 475 RETURN 476 END IF 477 ! ZX - ZX 478 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzp .AND. m2 == CLMzp) THEN 479 da = da/2.0_dp 480 db = db/2.0_dp 481 qxzqxz = & 482 +1.0_dp/SQRT((r + da - db)**2 + (da - db)**2 + add) - 1.0_dp/SQRT((r + da + db)**2 + (da - db)**2 + add) & 483 - 1.0_dp/SQRT((r - da - db)**2 + (da - db)**2 + add) + 1.0_dp/SQRT((r - da + db)**2 + (da - db)**2 + add) & 484 - 1.0_dp/SQRT((r + da - db)**2 + (da + db)**2 + add) + 1.0_dp/SQRT((r + da + db)**2 + (da + db)**2 + add) & 485 + 1.0_dp/SQRT((r - da - db)**2 + (da + db)**2 + add) - 1.0_dp/SQRT((r - da + db)**2 + (da + db)**2 + add) 486 charg = qxzqxz*0.125_dp*fact 487 RETURN 488 END IF 489 ! XX - XX 490 IF (l1 == 2 .AND. l2 == 2 .AND. (((m1 == CLMyy) .AND. (m2 == CLMyy)) .OR. ((m1 == CLMxx) .AND. (m2 == CLMxx)))) THEN 491 qxxqxx = & 492 +2.0_dp/SQRT(r**2 + add) + 1.0_dp/SQRT(r**2 + (da - db)**2 + add) & 493 + 1.0_dp/SQRT(r**2 + (da + db)**2 + add) - 2.0_dp/SQRT(r**2 + da**2 + add) & 494 - 2.0_dp/SQRT(r**2 + db**2 + add) 495 charg = qxxqxx*0.125_dp*fact 496 RETURN 497 END IF 498 ! XX - YY 499 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMyy .AND. m2 == CLMxx) THEN 500 qxxqyy = & 501 +1.0_dp/SQRT(r**2 + add) - 1.0_dp/SQRT(r**2 + da**2 + add) & 502 - 1.0_dp/SQRT(r**2 + db**2 + add) + 1.0_dp/SQRT(r**2 + da**2 + db**2 + add) 503 charg = qxxqyy*0.25_dp*fact 504 RETURN 505 END IF 506 ! XY - XY 507 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMxy .AND. m2 == CLMxy) THEN 508 qxxqxx = & 509 +2.0_dp/SQRT(r**2 + add) + 1.0_dp/SQRT(r**2 + (da - db)**2 + add) & 510 + 1.0_dp/SQRT(r**2 + (da + db)**2 + add) - 2.0_dp/SQRT(r**2 + da**2 + add) & 511 - 2.0_dp/SQRT(r**2 + db**2 + add) 512 qxxqyy = & 513 +1.0_dp/SQRT(r**2 + add) - 1.0_dp/SQRT(r**2 + da**2 + add) & 514 - 1.0_dp/SQRT(r**2 + db**2 + add) + 1.0_dp/SQRT(r**2 + da**2 + db**2 + add) 515 charg = 0.5_dp*(qxxqxx*0.125_dp - qxxqyy*0.25_dp)*fact 516 RETURN 517 END IF 518 ! We should NEVER reach this point 519 CPABORT("") 520 END FUNCTION charg_int_nri 521 522! ************************************************************************************************** 523!> \brief Derivatives of interaction function between two point-charge 524!> configurations NDDO sp-code. 525!> Non-Rotational Invariant definition of quadrupoles 526!> 527!> r - Distance r12 528!> l1,m - Quantum numbers for multipole of configuration 1 529!> l2,m - Quantum numbers for multipole of configuration 2 530!> da - charge separation of configuration 1 531!> db - charge separation of configuration 2 532!> add - additive term 533!> 534!> \param r ... 535!> \param l1_i ... 536!> \param l2_i ... 537!> \param m1_i ... 538!> \param m2_i ... 539!> \param da_i ... 540!> \param db_i ... 541!> \param add0 ... 542!> \param fact_screen ... 543!> \return ... 544!> \date 04.2008 [tlaino] 545!> \author Teodoro Laino [tlaino] - University of Zurich 546! ************************************************************************************************** 547 FUNCTION dcharg_int_nri(r, l1_i, l2_i, m1_i, m2_i, da_i, db_i, add0, fact_screen) RESULT(charg) 548 REAL(KIND=dp), INTENT(in) :: r 549 INTEGER, INTENT(in) :: l1_i, l2_i, m1_i, m2_i 550 REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen 551 REAL(KIND=dp) :: charg 552 553 CHARACTER(len=*), PARAMETER :: routineN = 'dcharg_int_nri', routineP = moduleN//':'//routineN 554 555 INTEGER :: l1, l2, m1, m2 556 REAL(KIND=dp) :: add, da, db, dxdx, dxqxz, dzdz, dzqxx, & 557 dzqzz, fact, qqxx, qqzz, qxxqxx, & 558 qxxqyy, qxzqxz, xyxy, zzzz 559 560! Computing only Integral Derivatives 561 562 IF (l1_i < l2_i) THEN 563 l1 = l1_i 564 l2 = l2_i 565 m1 = m1_i 566 m2 = m2_i 567 da = da_i 568 db = db_i 569 fact = 1.0_dp 570 ELSE IF (l1_i > l2_i) THEN 571 l1 = l2_i 572 l2 = l1_i 573 m1 = m2_i 574 m2 = m1_i 575 da = db_i 576 db = da_i 577 fact = (-1.0_dp)**(l1 + l2) 578 ELSE IF (l1_i == l2_i) THEN 579 l1 = l1_i 580 l2 = l2_i 581 IF (m1_i <= m2_i) THEN 582 m1 = m1_i 583 m2 = m2_i 584 da = da_i 585 db = db_i 586 ELSE 587 m1 = m2_i 588 m2 = m1_i 589 da = db_i 590 db = da_i 591 END IF 592 fact = 1.0_dp 593 END IF 594 charg = 0.0_dp 595 add = add0*fact_screen 596 ! Q - Q. 597 IF (l1 == 0 .AND. l2 == 0) THEN 598 charg = r/SQRT(r**2 + add)**3 599 charg = -charg*fact 600 RETURN 601 END IF 602 ! Q - Z. 603 IF (l1 == 0 .AND. l2 == 1 .AND. m2 == CLMz) THEN 604 charg = (r + db)/SQRT((r + db)**2 + add)**3 - (r - db)/SQRT((r - db)**2 + add)**3 605 charg = -charg*0.5_dp*fact 606 RETURN 607 END IF 608 ! Z - Z. 609 IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMz .AND. m2 == CLMz) THEN 610 dzdz = & 611 +(r + da - db)/SQRT((r + da - db)**2 + add)**3 + (r - da + db)/SQRT((r - da + db)**2 + add)**3 & 612 - (r - da - db)/SQRT((r - da - db)**2 + add)**3 - (r + da + db)/SQRT((r + da + db)**2 + add)**3 613 charg = -dzdz*0.25_dp*fact 614 RETURN 615 END IF 616 ! X - X 617 IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMp .AND. m2 == CLMp) THEN 618 dxdx = 2.0_dp*r/SQRT(r**2 + (da - db)**2 + add)**3 - 2.0_dp*r/SQRT(r**2 + (da + db)**2 + add)**3 619 charg = -dxdx*0.25_dp*fact 620 RETURN 621 END IF 622 ! Q - ZZ 623 IF (l1 == 0 .AND. l2 == 2 .AND. m2 == CLMzz) THEN 624 qqzz = (r - db)/SQRT((r - db)**2 + add)**3 - 2.0_dp*r/SQRT(r**2 + add)**3 + (r + db)/SQRT((r + db)**2 + add)**3 625 charg = -qqzz*0.25_dp*fact 626 RETURN 627 END IF 628 ! Q - XX 629 IF (l1 == 0 .AND. l2 == 2 .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN 630 qqxx = -r/SQRT(r**2 + add)**3 + r/SQRT(r**2 + add + db**2)**3 631 charg = -qqxx*0.5_dp*fact 632 RETURN 633 END IF 634 ! Z - ZZ 635 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. m2 == CLMzz) THEN 636 dzqzz = & 637 +(r - da - db)/SQRT((r - da - db)**2 + add)**3 - 2.0_dp*(r - da)/SQRT((r - da)**2 + add)**3 & 638 + (r - da + db)/SQRT((r - da + db)**2 + add)**3 - (r + da - db)/SQRT((r + da - db)**2 + add)**3 & 639 + 2.0_dp*(r + da)/SQRT((r + da)**2 + add)**3 - (r + da + db)/SQRT((r + da + db)**2 + add)**3 640 charg = -dzqzz*0.125_dp*fact 641 RETURN 642 END IF 643 ! Z - XX 644 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN 645 dzqxx = & 646 +(r + da)/SQRT((r + da)**2 + add)**3 - (r + da)/SQRT((r + da)**2 + add + db**2)**3 & 647 - (r - da)/SQRT((r - da)**2 + add)**3 + (r - da)/SQRT((r - da)**2 + add + db**2)**3 648 charg = -dzqxx*0.25_dp*fact 649 RETURN 650 END IF 651 ! ZZ - ZZ 652 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. m2 == CLMzz) THEN 653 zzzz = & 654 +(r - da - db)/SQRT((r - da - db)**2 + add)**3 + (r + da + db)/SQRT((r + da + db)**2 + add)**3 & 655 + (r - da + db)/SQRT((r - da + db)**2 + add)**3 + (r + da - db)/SQRT((r + da - db)**2 + add)**3 656 xyxy = & 657 +(r - da)/SQRT((r - da)**2 + add)**3 + (r + da)/SQRT((r + da)**2 + add)**3 & 658 + (r - db)/SQRT((r - db)**2 + add)**3 + (r + db)/SQRT((r + db)**2 + add)**3 & 659 - 2.0_dp*r/SQRT(r**2 + add)**3 660 charg = -(zzzz*0.0625_dp - xyxy*0.125_dp)*fact 661 RETURN 662 END IF 663 ! ZZ - XX 664 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. (m2 == CLMxx .OR. m2 == CLMyy)) THEN 665 zzzz = & 666 -(r + da)/SQRT((r + da)**2 + add)**3 + (r + da)/SQRT((r + da)**2 + db**2 + add)**3 & 667 - (r - da)/SQRT((r - da)**2 + add)**3 + (r - da)/SQRT((r - da)**2 + db**2 + add)**3 668 xyxy = r/SQRT(r**2 + db**2 + add)**3 - r/SQRT(r**2 + add)**3 669 charg = -(zzzz*0.125_dp - xyxy*0.25_dp)*fact 670 RETURN 671 END IF 672 ! X - ZX 673 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMp .AND. m2 == CLMzp) THEN 674 db = db/2.0_dp 675 dxqxz = & 676 -(r - db)/SQRT((r - db)**2 + (da - db)**2 + add)**3 + (r + db)/SQRT((r + db)**2 + (da - db)**2 + add)**3 & 677 + (r - db)/SQRT((r - db)**2 + (da + db)**2 + add)**3 - (r + db)/SQRT((r + db)**2 + (da + db)**2 + add)**3 678 charg = -dxqxz*0.25_dp*fact 679 RETURN 680 END IF 681 ! ZX - ZX 682 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzp .AND. m2 == CLMzp) THEN 683 da = da/2.0_dp 684 db = db/2.0_dp 685 qxzqxz = & 686 +(r + da - db)/SQRT((r + da - db)**2 + (da - db)**2 + add)**3 - (r + da + db)/SQRT((r + da + db)**2 + (da - db)**2 + add)**3 & 687 - (r - da - db)/SQRT((r - da - db)**2 + (da - db)**2 + add)**3 + (r - da + db)/SQRT((r - da + db)**2 + (da - db)**2 + add)**3 & 688 - (r + da - db)/SQRT((r + da - db)**2 + (da + db)**2 + add)**3 + (r + da + db)/SQRT((r + da + db)**2 + (da + db)**2 + add)**3 & 689 + (r - da - db)/SQRT((r - da - db)**2 + (da + db)**2 + add)**3 - (r - da + db)/SQRT((r - da + db)**2 + (da + db)**2 + add)**3 690 charg = -qxzqxz*0.125_dp*fact 691 RETURN 692 END IF 693 ! XX - XX 694 IF (l1 == 2 .AND. l2 == 2 .AND. (((m1 == CLMyy) .AND. (m2 == CLMyy)) .OR. ((m1 == CLMxx) .AND. (m2 == CLMxx)))) THEN 695 qxxqxx = & 696 +2.0_dp*r/SQRT(r**2 + add)**3 + r/SQRT(r**2 + (da - db)**2 + add)**3 & 697 + r/SQRT(r**2 + (da + db)**2 + add)**3 - 2.0_dp*r/SQRT(r**2 + da**2 + add)**3 & 698 - 2.0_dp*r/SQRT(r**2 + db**2 + add)**3 699 charg = -qxxqxx*0.125_dp*fact 700 RETURN 701 END IF 702 ! XX - YY 703 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMyy .AND. m2 == CLMxx) THEN 704 qxxqyy = & 705 +r/SQRT(r**2 + add)**3 - r/SQRT(r**2 + da**2 + add)**3 & 706 - r/SQRT(r**2 + db**2 + add)**3 + r/SQRT(r**2 + da**2 + db**2 + add)**3 707 charg = -qxxqyy*0.25_dp*fact 708 RETURN 709 END IF 710 ! XY - XY 711 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMxy .AND. m2 == CLMxy) THEN 712 qxxqxx = & 713 +2.0_dp*r/SQRT(r**2 + add)**3 + r/SQRT(r**2 + (da - db)**2 + add)**3 & 714 + r/SQRT(r**2 + (da + db)**2 + add)**3 - 2.0_dp*r/SQRT(r**2 + da**2 + add)**3 & 715 - 2.0_dp*r/SQRT(r**2 + db**2 + add)**3 716 qxxqyy = & 717 +r/SQRT(r**2 + add)**3 - r/SQRT(r**2 + da**2 + add)**3 & 718 - r/SQRT(r**2 + db**2 + add)**3 + r/SQRT(r**2 + da**2 + db**2 + add)**3 719 charg = -0.5_dp*(qxxqxx*0.125_dp - qxxqyy*0.25_dp)*fact 720 RETURN 721 END IF 722 ! We should NEVER reach this point 723 CPABORT("") 724 END FUNCTION dcharg_int_nri 725 726! ************************************************************************************************** 727!> \brief Derivatives of interaction function between two point-charge 728!> configurations NDDO sp-code. The derivative takes care of the screening 729!> term only. 730!> Non-Rotational Invariant definition of quadrupoles 731!> 732!> r - Distance r12 733!> l1,m - Quantum numbers for multipole of configuration 1 734!> l2,m - Quantum numbers for multipole of configuration 2 735!> da - charge separation of configuration 1 736!> db - charge separation of configuration 2 737!> add - additive term 738!> 739!> \param r ... 740!> \param l1_i ... 741!> \param l2_i ... 742!> \param m1_i ... 743!> \param m2_i ... 744!> \param da_i ... 745!> \param db_i ... 746!> \param add0 ... 747!> \param fact_screen ... 748!> \return ... 749!> \date 04.2008 [tlaino] 750!> \author Teodoro Laino [tlaino] - University of Zurich 751! ************************************************************************************************** 752 FUNCTION dcharg_int_nri_fs(r, l1_i, l2_i, m1_i, m2_i, da_i, db_i, add0, fact_screen) RESULT(charg) 753 REAL(KIND=dp), INTENT(in) :: r 754 INTEGER, INTENT(in) :: l1_i, l2_i, m1_i, m2_i 755 REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen 756 REAL(KIND=dp) :: charg 757 758 CHARACTER(len=*), PARAMETER :: routineN = 'dcharg_int_nri_fs', & 759 routineP = moduleN//':'//routineN 760 761 INTEGER :: l1, l2, m1, m2 762 REAL(KIND=dp) :: add, da, db, dxdx, dxqxz, dzdz, dzqxx, & 763 dzqzz, fact, qqxx, qqzz, qxxqxx, & 764 qxxqyy, qxzqxz, xyxy, zzzz 765 766! Computing only Integral Derivatives 767 768 IF (l1_i < l2_i) THEN 769 l1 = l1_i 770 l2 = l2_i 771 m1 = m1_i 772 m2 = m2_i 773 da = da_i 774 db = db_i 775 fact = 1.0_dp 776 ELSE IF (l1_i > l2_i) THEN 777 l1 = l2_i 778 l2 = l1_i 779 m1 = m2_i 780 m2 = m1_i 781 da = db_i 782 db = da_i 783 fact = (-1.0_dp)**(l1 + l2) 784 ELSE IF (l1_i == l2_i) THEN 785 l1 = l1_i 786 l2 = l2_i 787 IF (m1_i <= m2_i) THEN 788 m1 = m1_i 789 m2 = m2_i 790 da = da_i 791 db = db_i 792 ELSE 793 m1 = m2_i 794 m2 = m1_i 795 da = db_i 796 db = da_i 797 END IF 798 fact = 1.0_dp 799 END IF 800 charg = 0.0_dp 801 add = add0*fact_screen 802 ! The 0.5 factor handles the derivative of the SQRT 803 fact = fact*0.5_dp 804 ! Q - Q. 805 IF (l1 == 0 .AND. l2 == 0) THEN 806 charg = add0/SQRT(r**2 + add)**3 807 charg = -charg*fact 808 RETURN 809 END IF 810 ! Q - Z. 811 IF (l1 == 0 .AND. l2 == 1 .AND. m2 == CLMz) THEN 812 charg = add0/SQRT((r + db)**2 + add)**3 - add0/SQRT((r - db)**2 + add)**3 813 charg = -charg*0.5_dp*fact 814 RETURN 815 END IF 816 ! Z - Z. 817 IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMz .AND. m2 == CLMz) THEN 818 dzdz = & 819 +add0/SQRT((r + da - db)**2 + add)**3 + add0/SQRT((r - da + db)**2 + add)**3 & 820 - add0/SQRT((r - da - db)**2 + add)**3 - add0/SQRT((r + da + db)**2 + add)**3 821 charg = -dzdz*0.25_dp*fact 822 RETURN 823 END IF 824 ! X - X 825 IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMp .AND. m2 == CLMp) THEN 826 dxdx = 2.0_dp*add0/SQRT(r**2 + (da - db)**2 + add)**3 - 2.0_dp*add0/SQRT(r**2 + (da + db)**2 + add)**3 827 charg = -dxdx*0.25_dp*fact 828 RETURN 829 END IF 830 ! Q - ZZ 831 IF (l1 == 0 .AND. l2 == 2 .AND. m2 == CLMzz) THEN 832 qqzz = add0/SQRT((r - db)**2 + add)**3 - 2.0_dp*add0/SQRT(r**2 + add)**3 + add0/SQRT((r + db)**2 + add)**3 833 charg = -qqzz*0.25_dp*fact 834 RETURN 835 END IF 836 ! Q - XX 837 IF (l1 == 0 .AND. l2 == 2 .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN 838 qqxx = -add0/SQRT(r**2 + add)**3 + add0/SQRT(r**2 + add + db**2)**3 839 charg = -qqxx*0.5_dp*fact 840 RETURN 841 END IF 842 ! Z - ZZ 843 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. m2 == CLMzz) THEN 844 dzqzz = & 845 +add0/SQRT((r - da - db)**2 + add)**3 - 2.0_dp*add0/SQRT((r - da)**2 + add)**3 & 846 + add0/SQRT((r - da + db)**2 + add)**3 - add0/SQRT((r + da - db)**2 + add)**3 & 847 + 2.0_dp*add0/SQRT((r + da)**2 + add)**3 - add0/SQRT((r + da + db)**2 + add)**3 848 charg = -dzqzz*0.125_dp*fact 849 RETURN 850 END IF 851 ! Z - XX 852 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN 853 dzqxx = & 854 +add0/SQRT((r + da)**2 + add)**3 - add0/SQRT((r + da)**2 + add + db**2)**3 & 855 - add0/SQRT((r - da)**2 + add)**3 + add0/SQRT((r - da)**2 + add + db**2)**3 856 charg = -dzqxx*0.25_dp*fact 857 RETURN 858 END IF 859 ! ZZ - ZZ 860 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. m2 == CLMzz) THEN 861 zzzz = & 862 +add0/SQRT((r - da - db)**2 + add)**3 + add0/SQRT((r + da + db)**2 + add)**3 & 863 + add0/SQRT((r - da + db)**2 + add)**3 + add0/SQRT((r + da - db)**2 + add)**3 864 xyxy = & 865 +add0/SQRT((r - da)**2 + add)**3 + add0/SQRT((r + da)**2 + add)**3 & 866 + add0/SQRT((r - db)**2 + add)**3 + add0/SQRT((r + db)**2 + add)**3 & 867 - 2.0_dp*add0/SQRT(r**2 + add)**3 868 charg = -(zzzz*0.0625_dp - xyxy*0.125_dp)*fact 869 RETURN 870 END IF 871 ! ZZ - XX 872 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. (m2 == CLMxx .OR. m2 == CLMyy)) THEN 873 zzzz = & 874 -add0/SQRT((r + da)**2 + add)**3 + add0/SQRT((r + da)**2 + db**2 + add)**3 & 875 - add0/SQRT((r - da)**2 + add)**3 + add0/SQRT((r - da)**2 + db**2 + add)**3 876 xyxy = add0/SQRT(r**2 + db**2 + add)**3 - add0/SQRT(r**2 + add)**3 877 charg = -(zzzz*0.125_dp - xyxy*0.25_dp)*fact 878 RETURN 879 END IF 880 ! X - ZX 881 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMp .AND. m2 == CLMzp) THEN 882 db = db/2.0_dp 883 dxqxz = & 884 -add0/SQRT((r - db)**2 + (da - db)**2 + add)**3 + add0/SQRT((r + db)**2 + (da - db)**2 + add)**3 & 885 + add0/SQRT((r - db)**2 + (da + db)**2 + add)**3 - add0/SQRT((r + db)**2 + (da + db)**2 + add)**3 886 charg = -dxqxz*0.25_dp*fact 887 RETURN 888 END IF 889 ! ZX - ZX 890 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzp .AND. m2 == CLMzp) THEN 891 da = da/2.0_dp 892 db = db/2.0_dp 893 qxzqxz = & 894 +add0/SQRT((r + da - db)**2 + (da - db)**2 + add)**3 - add0/SQRT((r + da + db)**2 + (da - db)**2 + add)**3 & 895 - add0/SQRT((r - da - db)**2 + (da - db)**2 + add)**3 + add0/SQRT((r - da + db)**2 + (da - db)**2 + add)**3 & 896 - add0/SQRT((r + da - db)**2 + (da + db)**2 + add)**3 + add0/SQRT((r + da + db)**2 + (da + db)**2 + add)**3 & 897 + add0/SQRT((r - da - db)**2 + (da + db)**2 + add)**3 - add0/SQRT((r - da + db)**2 + (da + db)**2 + add)**3 898 charg = -qxzqxz*0.125_dp*fact 899 RETURN 900 END IF 901 ! XX - XX 902 IF (l1 == 2 .AND. l2 == 2 .AND. (((m1 == CLMyy) .AND. (m2 == CLMyy)) .OR. ((m1 == CLMxx) .AND. (m2 == CLMxx)))) THEN 903 qxxqxx = & 904 +2.0_dp*add0/SQRT(r**2 + add)**3 + add0/SQRT(r**2 + (da - db)**2 + add)**3 & 905 + add0/SQRT(r**2 + (da + db)**2 + add)**3 - 2.0_dp*add0/SQRT(r**2 + da**2 + add)**3 & 906 - 2.0_dp*add0/SQRT(r**2 + db**2 + add)**3 907 charg = -qxxqxx*0.125_dp*fact 908 RETURN 909 END IF 910 ! XX - YY 911 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMyy .AND. m2 == CLMxx) THEN 912 qxxqyy = & 913 +add0/SQRT(r**2 + add)**3 - add0/SQRT(r**2 + da**2 + add)**3 & 914 - add0/SQRT(r**2 + db**2 + add)**3 + add0/SQRT(r**2 + da**2 + db**2 + add)**3 915 charg = -qxxqyy*0.25_dp*fact 916 RETURN 917 END IF 918 ! XY - XY 919 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMxy .AND. m2 == CLMxy) THEN 920 qxxqxx = & 921 +2.0_dp*add0/SQRT(r**2 + add)**3 + add0/SQRT(r**2 + (da - db)**2 + add)**3 & 922 + add0/SQRT(r**2 + (da + db)**2 + add)**3 - 2.0_dp*add0/SQRT(r**2 + da**2 + add)**3 & 923 - 2.0_dp*add0/SQRT(r**2 + db**2 + add)**3 924 qxxqyy = & 925 +add0/SQRT(r**2 + add)**3 - add0/SQRT(r**2 + da**2 + add)**3 & 926 - add0/SQRT(r**2 + db**2 + add)**3 + add0/SQRT(r**2 + da**2 + db**2 + add)**3 927 charg = -0.5_dp*(qxxqxx*0.125_dp - qxxqyy*0.25_dp)*fact 928 RETURN 929 END IF 930 ! We should NEVER reach this point 931 CPABORT("") 932 END FUNCTION dcharg_int_nri_fs 933 934! ************************************************************************************************** 935!> \brief General driver for computing semi-empirical integrals <ij|kl> 936!> involving d-orbitals. 937!> The choice of the linear quadrupole was REALLY unhappy 938!> in the first development of the NDDO codes. That choice makes 939!> impossible the merging of the integral code with or without d-orbitals 940!> unless a reparametrization of all NDDO codes for s and p orbitals will 941!> be performed.. more over the choice of the linear quadrupole does not make 942!> calculations rotational invariants (of course the rotational invariant 943!> can be forced). The definitions of quadrupoles for d-orbitals is the 944!> correct one in order to have the rotational invariant property by 945!> construction.. 946!> 947!> \param sepi ... 948!> \param sepj ... 949!> \param ij ... 950!> \param kl ... 951!> \param li ... 952!> \param lj ... 953!> \param lk ... 954!> \param ll ... 955!> \param ic ... 956!> \param r ... 957!> \param se_int_control ... 958!> \param se_int_screen ... 959!> \param itype ... 960!> \return ... 961!> \date 03.2008 [tlaino] 962!> \author Teodoro Laino [tlaino] - University of Zurich 963! ************************************************************************************************** 964 FUNCTION ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control, & 965 se_int_screen, itype) RESULT(res) 966 TYPE(semi_empirical_type), POINTER :: sepi, sepj 967 INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic 968 REAL(KIND=dp), INTENT(IN) :: r 969 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 970 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 971 INTEGER, INTENT(IN) :: itype 972 REAL(KIND=dp) :: res 973 974 CHARACTER(len=*), PARAMETER :: routineN = 'ijkl_d', routineP = moduleN//':'//routineN 975 976 res = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, & 977 se_int_control%integral_screening, se_int_control%shortrange, & 978 se_int_control%pc_coulomb_int, se_int_control%max_multipole, & 979 itype, charg_int_ri) 980 981 ! If only the shortrange component is requested we can skip the rest 982 IF ((.NOT. se_int_control%pc_coulomb_int) .AND. (itype /= do_method_pchg)) THEN 983 ! Handle the 1/r^3 term, this term is ALWAYS false for KDSO-D integrals 984 IF (se_int_control%shortrange .AND. se_int_control%do_ewald_r3) THEN 985 res = res - ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, & 986 itype, charg_int_3) 987 END IF 988 END IF 989 END FUNCTION ijkl_d 990 991! ************************************************************************************************** 992!> \brief General driver for computing the derivatives of semi-empirical integrals <ij|kl> 993!> involving d-orbitals. 994!> The choice of the linear quadrupole was REALLY unhappy 995!> in the first development of the NDDO codes. That choice makes 996!> impossible the merging of the integral code with or without d-orbitals 997!> unless a reparametrization of all NDDO codes for s and p orbitals will 998!> be performed.. more over the choice of the linear quadrupole does not make 999!> calculations rotational invariants (of course the rotational invariant 1000!> can be forced). The definitions of quadrupoles for d-orbitals is the 1001!> correct one in order to have the rotational invariant property by 1002!> construction.. 1003!> 1004!> \param sepi ... 1005!> \param sepj ... 1006!> \param ij ... 1007!> \param kl ... 1008!> \param li ... 1009!> \param lj ... 1010!> \param lk ... 1011!> \param ll ... 1012!> \param ic ... 1013!> \param r ... 1014!> \param se_int_control ... 1015!> \param se_int_screen ... 1016!> \param itype ... 1017!> \return ... 1018!> \date 03.2008 [tlaino] 1019!> \author Teodoro Laino [tlaino] - University of Zurich 1020! ************************************************************************************************** 1021 FUNCTION d_ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control, & 1022 se_int_screen, itype) RESULT(res) 1023 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1024 INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic 1025 REAL(KIND=dp), INTENT(IN) :: r 1026 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 1027 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 1028 INTEGER, INTENT(IN) :: itype 1029 REAL(KIND=dp) :: res 1030 1031 CHARACTER(len=*), PARAMETER :: routineN = 'd_ijkl_d', routineP = moduleN//':'//routineN 1032 1033 REAL(KIND=dp) :: dfs, srd 1034 1035 IF (se_int_control%integral_screening == do_se_IS_kdso_d) THEN 1036 res = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, & 1037 se_int_control%integral_screening, .FALSE., & 1038 se_int_control%pc_coulomb_int, se_int_control%max_multipole, & 1039 itype, dcharg_int_ri) 1040 1041 IF (.NOT. se_int_control%pc_coulomb_int) THEN 1042 dfs = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, & 1043 se_int_control%integral_screening, .FALSE., .FALSE., & 1044 se_int_control%max_multipole, itype, dcharg_int_ri_fs) 1045 res = res + dfs*se_int_screen%dft 1046 1047 ! In case we need the shortrange part we have to evaluate an additional derivative 1048 ! to handle the derivative of the Tapering term 1049 IF (se_int_control%shortrange) THEN 1050 srd = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, & 1051 se_int_control%integral_screening, .FALSE., .TRUE., & 1052 se_int_control%max_multipole, itype, dcharg_int_ri) 1053 res = res - srd 1054 END IF 1055 END IF 1056 ELSE 1057 res = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, & 1058 se_int_control%integral_screening, se_int_control%shortrange, & 1059 se_int_control%pc_coulomb_int, se_int_control%max_multipole, & 1060 itype, dcharg_int_ri) 1061 END IF 1062 1063 ! If only the shortrange component is requested we can skip the rest 1064 IF ((.NOT. se_int_control%pc_coulomb_int) .AND. (itype /= do_method_pchg)) THEN 1065 ! Handle the 1/r^3 term, this term is ALWAYS false for KDSO-D integrals 1066 IF (se_int_control%shortrange .AND. se_int_control%do_ewald_r3) THEN 1067 res = res - ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, & 1068 itype, dcharg_int_3) 1069 END IF 1070 END IF 1071 1072 END FUNCTION d_ijkl_d 1073 1074! ************************************************************************************************** 1075!> \brief Low level general driver for computing semi-empirical integrals <ij|kl> 1076!> and their derivatives involving d-orbitals. 1077!> The choice of the linear quadrupole was REALLY unhappy 1078!> in the first development of the NDDO codes. That choice makes 1079!> impossible the merging of the integral code with or without d-orbitals 1080!> unless a reparametrization of all NDDO codes for s and p orbitals will 1081!> be performed.. more over the choice of the linear quadrupole does not make 1082!> calculations rotational invariants (of course the rotational invariant 1083!> can be forced). The definitions of quadrupoles for d-orbitals is the 1084!> correct one in order to have the rotational invariant property by 1085!> construction.. 1086!> 1087!> \param sepi ... 1088!> \param sepj ... 1089!> \param ij ... 1090!> \param kl ... 1091!> \param li ... 1092!> \param lj ... 1093!> \param lk ... 1094!> \param ll ... 1095!> \param ic ... 1096!> \param r ... 1097!> \param se_int_screen ... 1098!> \param iscreen ... 1099!> \param shortrange ... 1100!> \param pc_coulomb_int ... 1101!> \param max_multipole ... 1102!> \param itype ... 1103!> \param eval a function without explicit interface 1104!> \return ... 1105!> \date 03.2008 [tlaino] 1106!> \author Teodoro Laino [tlaino] - University of Zurich 1107! ************************************************************************************************** 1108 FUNCTION ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen, & 1109 iscreen, shortrange, pc_coulomb_int, max_multipole, itype, eval) RESULT(res) 1110 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1111 INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic 1112 REAL(KIND=dp), INTENT(IN) :: r 1113 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 1114 INTEGER, INTENT(IN) :: iscreen 1115 LOGICAL, INTENT(IN) :: shortrange, pc_coulomb_int 1116 INTEGER, INTENT(IN) :: max_multipole, itype 1117 REAL(KIND=dp) :: eval, res 1118 1119 CHARACTER(len=*), PARAMETER :: routineN = 'ijkl_d_low', routineP = moduleN//':'//routineN 1120 1121 INTEGER :: l1, l1max, l1min, l2, l2max, l2min, lij, & 1122 lkl, lmin, m, mm 1123 REAL(KIND=dp) :: add, ccc, chrg, dij, dkl, fact_screen, & 1124 pij, pkl, s1, sum 1125 1126 l1min = ABS(li - lj) 1127 l1max = li + lj 1128 lij = indexb(li + 1, lj + 1) 1129 l2min = ABS(lk - ll) 1130 l2max = lk + ll 1131 lkl = indexb(lk + 1, ll + 1) 1132 1133 l1max = MIN(l1max, 2) 1134 l1min = MIN(l1min, 2) 1135 l2max = MIN(l2max, 2) 1136 l2min = MIN(l2min, 2) 1137 sum = 0.0_dp 1138 dij = 0.0_dp 1139 dkl = 0.0_dp 1140 fact_screen = 1.0_dp 1141 IF (.NOT. pc_coulomb_int) THEN 1142 IF (iscreen == do_se_IS_kdso_d) fact_screen = se_int_screen%ft 1143 ! Standard value of the integral 1144 DO l1 = l1min, l1max 1145 IF (l1 == 0) THEN 1146 IF (lij == 1) THEN 1147 pij = sepi%ko(1) 1148 IF (ic == 1) THEN 1149 pij = sepi%ko(9) 1150 END IF 1151 ELSE IF (lij == 3) THEN 1152 pij = sepi%ko(7) 1153 ELSE IF (lij == 6) THEN 1154 pij = sepi%ko(8) 1155 END IF 1156 ELSE 1157 dij = sepi%cs(lij) 1158 pij = sepi%ko(lij) 1159 END IF 1160 ! 1161 DO l2 = l2min, l2max 1162 IF (l2 == 0) THEN 1163 IF (lkl == 1) THEN 1164 pkl = sepj%ko(1) 1165 IF (ic == 2) THEN 1166 pkl = sepj%ko(9) 1167 END IF 1168 ELSE IF (lkl == 3) THEN 1169 pkl = sepj%ko(7) 1170 ELSE IF (lkl == 6) THEN 1171 pkl = sepj%ko(8) 1172 END IF 1173 ELSE 1174 dkl = sepj%cs(lkl) 1175 pkl = sepj%ko(lkl) 1176 END IF 1177 IF (itype == do_method_pchg) THEN 1178 add = 0.0_dp 1179 ELSE 1180 add = (pij + pkl)**2 1181 END IF 1182 lmin = MIN(l1, l2) 1183 s1 = 0.0_dp 1184 DO m = -lmin, lmin 1185 ccc = clm_d(ij, l1, m)*clm_d(kl, l2, m) 1186 IF (ABS(ccc) > EPSILON(0.0_dp)) THEN 1187 mm = ABS(m) 1188 chrg = eval(r, l1, l2, mm, dij, dkl, add, fact_screen) 1189 s1 = s1 + chrg*ccc 1190 END IF 1191 END DO 1192 sum = sum + s1 1193 END DO 1194 END DO 1195 res = sum 1196 END IF 1197 ! Shortrange: Possibly computes pure Coulomb and subtract from the original integral valeu 1198 IF (shortrange .OR. pc_coulomb_int) THEN 1199 sum = 0.0_dp 1200 dij = 0.0_dp 1201 dkl = 0.0_dp 1202 add = 0.0_dp 1203 fact_screen = 0.0_dp 1204 DO l1 = l1min, l1max 1205 IF (l1 > max_multipole) CYCLE 1206 IF (l1 /= 0) THEN 1207 dij = sepi%cs(lij) 1208 END IF 1209 ! 1210 DO l2 = l2min, l2max 1211 IF (l2 > max_multipole) CYCLE 1212 IF (l2 /= 0) THEN 1213 dkl = sepj%cs(lkl) 1214 END IF 1215 lmin = MIN(l1, l2) 1216 s1 = 0.0_dp 1217 DO m = -lmin, lmin 1218 ccc = clm_d(ij, l1, m)*clm_d(kl, l2, m) 1219 IF (ABS(ccc) > EPSILON(0.0_dp)) THEN 1220 mm = ABS(m) 1221 chrg = eval(r, l1, l2, mm, dij, dkl, add, fact_screen) 1222 s1 = s1 + chrg*ccc 1223 END IF 1224 END DO 1225 sum = sum + s1 1226 END DO 1227 END DO 1228 IF (pc_coulomb_int) res = sum 1229 IF (shortrange) res = res - sum 1230 END IF 1231 END FUNCTION ijkl_d_low 1232 1233! ************************************************************************************************** 1234!> \brief Interaction function between two point-charge configurations (MNDO/d) 1235!> Rotational invariant property built-in in the quadrupole definition 1236!> r - Distance r12 1237!> l1,m - Quantum numbers for multipole of configuration 1 1238!> l2,m - Quantum numbers for multipole of configuration 2 1239!> da - charge separation of configuration 1 1240!> db - charge separation of configuration 2 1241!> add - additive term 1242!> 1243!> \param r ... 1244!> \param l1_i ... 1245!> \param l2_i ... 1246!> \param m ... 1247!> \param da_i ... 1248!> \param db_i ... 1249!> \param add0 ... 1250!> \param fact_screen ... 1251!> \return ... 1252!> \date 03.2008 [tlaino] 1253!> \author Teodoro Laino [tlaino] - University of Zurich 1254! ************************************************************************************************** 1255 FUNCTION charg_int_ri(r, l1_i, l2_i, m, da_i, db_i, add0, fact_screen) RESULT(charg) 1256 REAL(KIND=dp), INTENT(in) :: r 1257 INTEGER, INTENT(in) :: l1_i, l2_i, m 1258 REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen 1259 REAL(KIND=dp) :: charg 1260 1261 CHARACTER(len=*), PARAMETER :: routineN = 'charg_int_ri', routineP = moduleN//':'//routineN 1262 1263 INTEGER :: l1, l2 1264 REAL(KIND=dp) :: aa, ab, add, da, db, dxdx, dxqxz, dzdz, & 1265 dzqzz, fact, qqzz, qxzqxz, xyxy, zzzz 1266 1267 IF (l1_i <= l2_i) THEN 1268 l1 = l1_i 1269 l2 = l2_i 1270 da = da_i 1271 db = db_i 1272 fact = 1.0_dp 1273 ELSE IF (l1_i > l2_i) THEN 1274 l1 = l2_i 1275 l2 = l1_i 1276 da = db_i 1277 db = da_i 1278 fact = (-1.0_dp)**(l1 + l2) 1279 END IF 1280 charg = 0.0_dp 1281 add = add0*fact_screen 1282 ! Q - Q. 1283 IF (l1 == 0 .AND. l2 == 0) THEN 1284 charg = fact/SQRT(r**2 + add) 1285 RETURN 1286 END IF 1287 ! Q - Z. 1288 IF (l1 == 0 .AND. l2 == 1) THEN 1289 charg = 1.0_dp/SQRT((r + db)**2 + add) - 1.0_dp/SQRT((r - db)**2 + add) 1290 charg = charg*0.5_dp*fact 1291 RETURN 1292 END IF 1293 ! Z - Z. 1294 IF (l1 == 1 .AND. l2 == 1 .AND. m == 0) THEN 1295 dzdz = & 1296 +1.0_dp/SQRT((r + da - db)**2 + add) + 1.0_dp/SQRT((r - da + db)**2 + add) & 1297 - 1.0_dp/SQRT((r - da - db)**2 + add) - 1.0_dp/SQRT((r + da + db)**2 + add) 1298 charg = dzdz*0.25_dp*fact 1299 RETURN 1300 END IF 1301 ! X - X 1302 IF (l1 == 1 .AND. l2 == 1 .AND. m == 1) THEN 1303 dxdx = 2.0_dp/SQRT(r**2 + (da - db)**2 + add) - 2.0_dp/SQRT(r**2 + (da + db)**2 + add) 1304 charg = dxdx*0.25_dp*fact 1305 RETURN 1306 END IF 1307 ! Q - ZZ 1308 IF (l1 == 0 .AND. l2 == 2) THEN 1309 qqzz = 1.0_dp/SQRT((r - db)**2 + add) - 2.0_dp/SQRT(r**2 + db**2 + add) + 1.0_dp/SQRT((r + db)**2 + add) 1310 charg = qqzz*0.25_dp*fact 1311 RETURN 1312 END IF 1313 ! Z - ZZ 1314 IF (l1 == 1 .AND. l2 == 2 .AND. m == 0) THEN 1315 dzqzz = & 1316 +1.0_dp/SQRT((r - da - db)**2 + add) - 2.0_dp/SQRT((r - da)**2 + db**2 + add) & 1317 + 1.0_dp/SQRT((r + db - da)**2 + add) - 1.0_dp/SQRT((r - db + da)**2 + add) & 1318 + 2.0_dp/SQRT((r + da)**2 + db**2 + add) - 1.0_dp/SQRT((r + da + db)**2 + add) 1319 charg = dzqzz*0.125_dp*fact 1320 RETURN 1321 END IF 1322 ! ZZ - ZZ 1323 IF (l1 == 2 .AND. l2 == 2 .AND. m == 0) THEN 1324 zzzz = & 1325 +1.0_dp/SQRT((r - da - db)**2 + add) + 1.0_dp/SQRT((r + da + db)**2 + add) & 1326 + 1.0_dp/SQRT((r - da + db)**2 + add) + 1.0_dp/SQRT((r + da - db)**2 + add) & 1327 - 2.0_dp/SQRT((r - da)**2 + db**2 + add) - 2.0_dp/SQRT((r - db)**2 + da**2 + add) & 1328 - 2.0_dp/SQRT((r + da)**2 + db**2 + add) - 2.0_dp/SQRT((r + db)**2 + da**2 + add) & 1329 + 2.0_dp/SQRT(r**2 + (da - db)**2 + add) + 2.0_dp/SQRT(r**2 + (da + db)**2 + add) 1330 xyxy = & 1331 +4.0_dp/SQRT(r**2 + (da - db)**2 + add) + 4.0_dp/SQRT(r**2 + (da + db)**2 + add) & 1332 - 8.0_dp/SQRT(r**2 + da**2 + db**2 + add) 1333 charg = (zzzz*0.0625_dp - xyxy*0.015625_dp)*fact 1334 RETURN 1335 END IF 1336 ! X - ZX 1337 IF (l1 == 1 .AND. l2 == 2 .AND. m == 1) THEN 1338 ab = db/SQRT(2.0_dp) 1339 dxqxz = & 1340 -2.0_dp/SQRT((r - ab)**2 + (da - ab)**2 + add) + 2.0_dp/SQRT((r + ab)**2 + (da - ab)**2 + add) & 1341 + 2.0_dp/SQRT((r - ab)**2 + (da + ab)**2 + add) - 2.0_dp/SQRT((r + ab)**2 + (da + ab)**2 + add) 1342 charg = dxqxz*0.125_dp*fact 1343 RETURN 1344 END IF 1345 ! ZX - ZX 1346 IF (l1 == 2 .AND. l2 == 2 .AND. m == 1) THEN 1347 aa = da/SQRT(2.0_dp) 1348 ab = db/SQRT(2.0_dp) 1349 qxzqxz = & 1350 +2.0_dp/SQRT((r + aa - ab)**2 + (aa - ab)**2 + add) - 2.0_dp/SQRT((r + aa + ab)**2 + (aa - ab)**2 + add) & 1351 - 2.0_dp/SQRT((r - aa - ab)**2 + (aa - ab)**2 + add) + 2.0_dp/SQRT((r - aa + ab)**2 + (aa - ab)**2 + add) & 1352 - 2.0_dp/SQRT((r + aa - ab)**2 + (aa + ab)**2 + add) + 2.0_dp/SQRT((r + aa + ab)**2 + (aa + ab)**2 + add) & 1353 + 2.0_dp/SQRT((r - aa - ab)**2 + (aa + ab)**2 + add) - 2.0_dp/SQRT((r - aa + ab)**2 + (aa + ab)**2 + add) 1354 charg = qxzqxz*0.0625_dp*fact 1355 RETURN 1356 END IF 1357 ! XX - XX 1358 IF (l1 == 2 .AND. l2 == 2 .AND. m == 2) THEN 1359 xyxy = 4.0_dp/SQRT(r**2 + (da - db)**2 + add) + 4.0_dp/SQRT(r**2 + (da + db)**2 + add) - 8.0_dp/SQRT(r**2 + da**2 + db**2 + add) 1360 charg = xyxy*0.0625_dp*fact 1361 RETURN 1362 END IF 1363 ! We should NEVER reach this point 1364 CPABORT("") 1365 END FUNCTION charg_int_ri 1366 1367! ************************************************************************************************** 1368!> \brief Derivatives of the interaction function between two point-charge 1369!> configurations (MNDO/d) 1370!> Rotational invariant property built-in in the quadrupole definition 1371!> r - Distance r12 1372!> l1,m - Quantum numbers for multipole of configuration 1 1373!> l2,m - Quantum numbers for multipole of configuration 2 1374!> da - charge separation of configuration 1 1375!> db - charge separation of configuration 2 1376!> add - additive term 1377!> 1378!> \param r ... 1379!> \param l1_i ... 1380!> \param l2_i ... 1381!> \param m ... 1382!> \param da_i ... 1383!> \param db_i ... 1384!> \param add0 ... 1385!> \param fact_screen ... 1386!> \return ... 1387!> \date 03.2008 [tlaino] 1388!> \author Teodoro Laino [tlaino] - University of Zurich 1389! ************************************************************************************************** 1390 FUNCTION dcharg_int_ri(r, l1_i, l2_i, m, da_i, db_i, add0, fact_screen) RESULT(charg) 1391 REAL(KIND=dp), INTENT(in) :: r 1392 INTEGER, INTENT(in) :: l1_i, l2_i, m 1393 REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen 1394 REAL(KIND=dp) :: charg 1395 1396 CHARACTER(len=*), PARAMETER :: routineN = 'dcharg_int_ri', routineP = moduleN//':'//routineN 1397 1398 INTEGER :: l1, l2 1399 REAL(KIND=dp) :: aa, ab, add, da, db, dxdx, dxqxz, dzdz, & 1400 dzqzz, fact, qqzz, qxzqxz, xyxy, zzzz 1401 1402 IF (l1_i <= l2_i) THEN 1403 l1 = l1_i 1404 l2 = l2_i 1405 da = da_i 1406 db = db_i 1407 fact = 1.0_dp 1408 ELSE IF (l1_i > l2_i) THEN 1409 l1 = l2_i 1410 l2 = l1_i 1411 da = db_i 1412 db = da_i 1413 fact = (-1.0_dp)**(l1 + l2) 1414 END IF 1415 charg = 0.0_dp 1416 add = add0*fact_screen 1417 ! Q - Q. 1418 IF (l1 == 0 .AND. l2 == 0) THEN 1419 charg = r/SQRT(r**2 + add)**3 1420 charg = -fact*charg 1421 RETURN 1422 END IF 1423 ! Q - Z. 1424 IF (l1 == 0 .AND. l2 == 1) THEN 1425 charg = (r + db)/SQRT((r + db)**2 + add)**3 - (r - db)/SQRT((r - db)**2 + add)**3 1426 charg = -charg*0.5_dp*fact 1427 RETURN 1428 END IF 1429 ! Z - Z. 1430 IF (l1 == 1 .AND. l2 == 1 .AND. m == 0) THEN 1431 dzdz = & 1432 +(r + da - db)/SQRT((r + da - db)**2 + add)**3 + (r - da + db)/SQRT((r - da + db)**2 + add)**3 & 1433 - (r - da - db)/SQRT((r - da - db)**2 + add)**3 - (r + da + db)/SQRT((r + da + db)**2 + add)**3 1434 charg = -dzdz*0.25_dp*fact 1435 RETURN 1436 END IF 1437 ! X - X 1438 IF (l1 == 1 .AND. l2 == 1 .AND. m == 1) THEN 1439 dxdx = 2.0_dp*r/SQRT(r**2 + (da - db)**2 + add)**3 - 2.0_dp*r/SQRT(r**2 + (da + db)**2 + add)**3 1440 charg = -dxdx*0.25_dp*fact 1441 RETURN 1442 END IF 1443 ! Q - ZZ 1444 IF (l1 == 0 .AND. l2 == 2) THEN 1445 qqzz = (r - db)/SQRT((r - db)**2 + add)**3 - 2.0_dp*r/SQRT(r**2 + db**2 + add)**3 + (r + db)/SQRT((r + db)**2 + add)**3 1446 charg = -qqzz*0.25_dp*fact 1447 RETURN 1448 END IF 1449 ! Z - ZZ 1450 IF (l1 == 1 .AND. l2 == 2 .AND. m == 0) THEN 1451 dzqzz = & 1452 +(r - da - db)/SQRT((r - da - db)**2 + add)**3 - 2.0_dp*(r - da)/SQRT((r - da)**2 + db**2 + add)**3 & 1453 + (r + db - da)/SQRT((r + db - da)**2 + add)**3 - (r - db + da)/SQRT((r - db + da)**2 + add)**3 & 1454 + 2.0_dp*(r + da)/SQRT((r + da)**2 + db**2 + add)**3 - (r + da + db)/SQRT((r + da + db)**2 + add)**3 1455 charg = -dzqzz*0.125_dp*fact 1456 RETURN 1457 END IF 1458 ! ZZ - ZZ 1459 IF (l1 == 2 .AND. l2 == 2 .AND. m == 0) THEN 1460 zzzz = & 1461 +(r - da - db)/SQRT((r - da - db)**2 + add)**3 + (r + da + db)/SQRT((r + da + db)**2 + add)**3 & 1462 + (r - da + db)/SQRT((r - da + db)**2 + add)**3 + (r + da - db)/SQRT((r + da - db)**2 + add)**3 & 1463 - 2.0_dp*(r - da)/SQRT((r - da)**2 + db**2 + add)**3 - 2.0_dp*(r - db)/SQRT((r - db)**2 + da**2 + add)**3 & 1464 - 2.0_dp*(r + da)/SQRT((r + da)**2 + db**2 + add)**3 - 2.0_dp*(r + db)/SQRT((r + db)**2 + da**2 + add)**3 & 1465 + 2.0_dp*r/SQRT(r**2 + (da - db)**2 + add)**3 + 2.0_dp*r/SQRT(r**2 + (da + db)**2 + add)**3 1466 xyxy = & 1467 +4.0_dp*r/SQRT(r**2 + (da - db)**2 + add)**3 + 4.0_dp*r/SQRT(r**2 + (da + db)**2 + add)**3 & 1468 - 8.0_dp*r/SQRT(r**2 + da**2 + db**2 + add)**3 1469 charg = -(zzzz*0.0625_dp - xyxy*0.015625_dp)*fact 1470 RETURN 1471 END IF 1472 ! X - ZX 1473 IF (l1 == 1 .AND. l2 == 2 .AND. m == 1) THEN 1474 ab = db/SQRT(2.0_dp) 1475 dxqxz = & 1476 -2.0_dp*(r - ab)/SQRT((r - ab)**2 + (da - ab)**2 + add)**3 + 2.0_dp*(r + ab)/SQRT((r + ab)**2 + (da - ab)**2 + add)**3 & 1477 + 2.0_dp*(r - ab)/SQRT((r - ab)**2 + (da + ab)**2 + add)**3 - 2.0_dp*(r + ab)/SQRT((r + ab)**2 + (da + ab)**2 + add)**3 1478 charg = -dxqxz*0.125_dp*fact 1479 RETURN 1480 END IF 1481 ! ZX - ZX 1482 IF (l1 == 2 .AND. l2 == 2 .AND. m == 1) THEN 1483 aa = da/SQRT(2.0_dp) 1484 ab = db/SQRT(2.0_dp) 1485 qxzqxz = & 1486 +2.0_dp*(r+aa-ab)/SQRT((r+aa-ab)**2+(aa-ab)**2+add)**3-2.0_dp*(r+aa+ab)/SQRT((r+aa+ab)**2+(aa-ab)**2+add)**3 & 1487 -2.0_dp*(r-aa-ab)/SQRT((r-aa-ab)**2+(aa-ab)**2+add)**3+2.0_dp*(r-aa+ab)/SQRT((r-aa+ab)**2+(aa-ab)**2+add)**3 & 1488 -2.0_dp*(r+aa-ab)/SQRT((r+aa-ab)**2+(aa+ab)**2+add)**3+2.0_dp*(r+aa+ab)/SQRT((r+aa+ab)**2+(aa+ab)**2+add)**3 & 1489 +2.0_dp*(r-aa-ab)/SQRT((r-aa-ab)**2+(aa+ab)**2+add)**3-2.0_dp*(r-aa+ab)/SQRT((r-aa+ab)**2+(aa+ab)**2+add)**3 1490 charg = -qxzqxz*0.0625_dp*fact 1491 RETURN 1492 END IF 1493 ! XX - XX 1494 IF (l1 == 2 .AND. l2 == 2 .AND. m == 2) THEN 1495 xyxy = 4.0_dp*r/SQRT(r**2+(da-db)**2+add)**3+4.0_dp*r/SQRT(r**2+(da+db)**2+add)**3-8.0_dp*r/SQRT(r**2+da**2+db**2+add)**3 1496 charg = -xyxy*0.0625_dp*fact 1497 RETURN 1498 END IF 1499 ! We should NEVER reach this point 1500 CPABORT("") 1501 END FUNCTION dcharg_int_ri 1502 1503! ************************************************************************************************** 1504!> \brief Derivatives of the interaction function between two point-charge 1505!> configurations (MNDO/d). This derivative takes into account only the 1506!> tapering term 1507!> Rotational invariant property built-in in the quadrupole definition 1508!> r - Distance r12 1509!> l1,m - Quantum numbers for multipole of configuration 1 1510!> l2,m - Quantum numbers for multipole of configuration 2 1511!> da - charge separation of configuration 1 1512!> db - charge separation of configuration 2 1513!> add - additive term 1514!> 1515!> \param r ... 1516!> \param l1_i ... 1517!> \param l2_i ... 1518!> \param m ... 1519!> \param da_i ... 1520!> \param db_i ... 1521!> \param add0 ... 1522!> \param fact_screen ... 1523!> \return ... 1524!> \date 03.2008 [tlaino] 1525!> \author Teodoro Laino [tlaino] - University of Zurich 1526! ************************************************************************************************** 1527 FUNCTION dcharg_int_ri_fs(r, l1_i, l2_i, m, da_i, db_i, add0, fact_screen) RESULT(charg) 1528 REAL(KIND=dp), INTENT(in) :: r 1529 INTEGER, INTENT(in) :: l1_i, l2_i, m 1530 REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen 1531 REAL(KIND=dp) :: charg 1532 1533 CHARACTER(len=*), PARAMETER :: routineN = 'dcharg_int_ri_fs', & 1534 routineP = moduleN//':'//routineN 1535 1536 INTEGER :: l1, l2 1537 REAL(KIND=dp) :: aa, ab, add, da, db, dxdx, dxqxz, dzdz, & 1538 dzqzz, fact, qqzz, qxzqxz, xyxy, zzzz 1539 1540 IF (l1_i <= l2_i) THEN 1541 l1 = l1_i 1542 l2 = l2_i 1543 da = da_i 1544 db = db_i 1545 fact = 1.0_dp 1546 ELSE IF (l1_i > l2_i) THEN 1547 l1 = l2_i 1548 l2 = l1_i 1549 da = db_i 1550 db = da_i 1551 fact = (-1.0_dp)**(l1 + l2) 1552 END IF 1553 charg = 0.0_dp 1554 add = add0*fact_screen 1555 ! The 0.5 factor handles the derivative of the SQRT 1556 fact = fact*0.5_dp 1557 ! Q - Q. 1558 IF (l1 == 0 .AND. l2 == 0) THEN 1559 charg = add0/SQRT(r**2 + add)**3 1560 charg = -fact*charg 1561 RETURN 1562 END IF 1563 ! Q - Z. 1564 IF (l1 == 0 .AND. l2 == 1) THEN 1565 charg = add0/SQRT((r + db)**2 + add)**3 - add0/SQRT((r - db)**2 + add)**3 1566 charg = -charg*0.5_dp*fact 1567 RETURN 1568 END IF 1569 ! Z - Z. 1570 IF (l1 == 1 .AND. l2 == 1 .AND. m == 0) THEN 1571 dzdz = & 1572 +add0/SQRT((r + da - db)**2 + add)**3 + add0/SQRT((r - da + db)**2 + add)**3 & 1573 - add0/SQRT((r - da - db)**2 + add)**3 - add0/SQRT((r + da + db)**2 + add)**3 1574 charg = -dzdz*0.25_dp*fact 1575 RETURN 1576 END IF 1577 ! X - X 1578 IF (l1 == 1 .AND. l2 == 1 .AND. m == 1) THEN 1579 dxdx = 2.0_dp*add0/SQRT(r**2 + (da - db)**2 + add)**3 - 2.0_dp*add0/SQRT(r**2 + (da + db)**2 + add)**3 1580 charg = -dxdx*0.25_dp*fact 1581 RETURN 1582 END IF 1583 ! Q - ZZ 1584 IF (l1 == 0 .AND. l2 == 2) THEN 1585 qqzz = add0/SQRT((r - db)**2 + add)**3 - 2.0_dp*add0/SQRT(r**2 + db**2 + add)**3 + add0/SQRT((r + db)**2 + add)**3 1586 charg = -qqzz*0.25_dp*fact 1587 RETURN 1588 END IF 1589 ! Z - ZZ 1590 IF (l1 == 1 .AND. l2 == 2 .AND. m == 0) THEN 1591 dzqzz = & 1592 +add0/SQRT((r - da - db)**2 + add)**3 - 2.0_dp*add0/SQRT((r - da)**2 + db**2 + add)**3 & 1593 + add0/SQRT((r + db - da)**2 + add)**3 - add0/SQRT((r - db + da)**2 + add)**3 & 1594 + 2.0_dp*add0/SQRT((r + da)**2 + db**2 + add)**3 - add0/SQRT((r + da + db)**2 + add)**3 1595 charg = -dzqzz*0.125_dp*fact 1596 RETURN 1597 END IF 1598 ! ZZ - ZZ 1599 IF (l1 == 2 .AND. l2 == 2 .AND. m == 0) THEN 1600 zzzz = & 1601 +add0/SQRT((r - da - db)**2 + add)**3 + add0/SQRT((r + da + db)**2 + add)**3 & 1602 + add0/SQRT((r - da + db)**2 + add)**3 + add0/SQRT((r + da - db)**2 + add)**3 & 1603 - 2.0_dp*add0/SQRT((r - da)**2 + db**2 + add)**3 - 2.0_dp*add0/SQRT((r - db)**2 + da**2 + add)**3 & 1604 - 2.0_dp*add0/SQRT((r + da)**2 + db**2 + add)**3 - 2.0_dp*add0/SQRT((r + db)**2 + da**2 + add)**3 & 1605 + 2.0_dp*add0/SQRT(r**2 + (da - db)**2 + add)**3 + 2.0_dp*add0/SQRT(r**2 + (da + db)**2 + add)**3 1606 xyxy = & 1607 +4.0_dp*add0/SQRT(r**2 + (da - db)**2 + add)**3 + 4.0_dp*add0/SQRT(r**2 + (da + db)**2 + add)**3 & 1608 - 8.0_dp*add0/SQRT(r**2 + da**2 + db**2 + add)**3 1609 charg = -(zzzz*0.0625_dp - xyxy*0.015625_dp)*fact 1610 RETURN 1611 END IF 1612 ! X - ZX 1613 IF (l1 == 1 .AND. l2 == 2 .AND. m == 1) THEN 1614 ab = db/SQRT(2.0_dp) 1615 dxqxz = & 1616 -2.0_dp*add0/SQRT((r - ab)**2 + (da - ab)**2 + add)**3 + 2.0_dp*add0/SQRT((r + ab)**2 + (da - ab)**2 + add)**3 & 1617 + 2.0_dp*add0/SQRT((r - ab)**2 + (da + ab)**2 + add)**3 - 2.0_dp*add0/SQRT((r + ab)**2 + (da + ab)**2 + add)**3 1618 charg = -dxqxz*0.125_dp*fact 1619 RETURN 1620 END IF 1621 ! ZX - ZX 1622 IF (l1 == 2 .AND. l2 == 2 .AND. m == 1) THEN 1623 aa = da/SQRT(2.0_dp) 1624 ab = db/SQRT(2.0_dp) 1625 qxzqxz = & 1626 +2.0_dp*add0/SQRT((r + aa - ab)**2 + (aa - ab)**2 + add)**3 - 2.0_dp*add0/SQRT((r + aa + ab)**2 + (aa - ab)**2 + add)**3 & 1627 - 2.0_dp*add0/SQRT((r - aa - ab)**2 + (aa - ab)**2 + add)**3 + 2.0_dp*add0/SQRT((r - aa + ab)**2 + (aa - ab)**2 + add)**3 & 1628 - 2.0_dp*add0/SQRT((r + aa - ab)**2 + (aa + ab)**2 + add)**3 + 2.0_dp*add0/SQRT((r + aa + ab)**2 + (aa + ab)**2 + add)**3 & 1629 + 2.0_dp*add0/SQRT((r - aa - ab)**2 + (aa + ab)**2 + add)**3 - 2.0_dp*add0/SQRT((r - aa + ab)**2 + (aa + ab)**2 + add)**3 1630 charg = -qxzqxz*0.0625_dp*fact 1631 RETURN 1632 END IF 1633 ! XX - XX 1634 IF (l1 == 2 .AND. l2 == 2 .AND. m == 2) THEN 1635 xyxy = 4.0_dp*add0/SQRT(r**2 + (da - db)**2 + add)**3 + 4.0_dp*add0/SQRT(r**2 + (da + db)**2 + add)**3 - & 1636 8.0_dp*add0/SQRT(r**2 + da**2 + db**2 + add)**3 1637 charg = -xyxy*0.0625_dp*fact 1638 RETURN 1639 END IF 1640 ! We should NEVER reach this point 1641 CPABORT("") 1642 END FUNCTION dcharg_int_ri_fs 1643 1644! ************************************************************************************************** 1645!> \brief Computes the general rotation matrices for the integrals 1646!> between I and J (J>=I) 1647!> 1648!> \param sepi ... 1649!> \param sepj ... 1650!> \param rjiv ... 1651!> \param r ... 1652!> \param ij_matrix ... 1653!> \param do_derivatives ... 1654!> \param do_invert ... 1655!> \param debug_invert ... 1656!> \date 04.2008 [tlaino] 1657!> \author Teodoro Laino [tlaino] - University of Zurich 1658! ************************************************************************************************** 1659 RECURSIVE SUBROUTINE rotmat(sepi, sepj, rjiv, r, ij_matrix, do_derivatives, & 1660 do_invert, debug_invert) 1661 TYPE(semi_empirical_type), POINTER :: sepi, sepj 1662 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rjiv 1663 REAL(KIND=dp), INTENT(IN) :: r 1664 TYPE(rotmat_type), POINTER :: ij_matrix 1665 LOGICAL, INTENT(IN) :: do_derivatives 1666 LOGICAL, INTENT(OUT), OPTIONAL :: do_invert 1667 LOGICAL, INTENT(IN), OPTIONAL :: debug_invert 1668 1669 CHARACTER(len=*), PARAMETER :: routineN = 'rotmat', routineP = moduleN//':'//routineN 1670 1671 INTEGER :: imap(3), k, l 1672 LOGICAL :: dbg_inv, eval 1673 REAL(KIND=dp) :: b, c2a, c2b, ca, ca2, cb, cb2, check, & 1674 pt5sq3, r2i, s2a, s2b, sa, sb, sb2, & 1675 sqb, sqb2i, x11, x22, x33 1676 REAL(KIND=dp), DIMENSION(3) :: b_d, c2a_d, c2b_d, ca2_d, ca_d, cb2_d, & 1677 cb_d, r_d, s2a_d, s2b_d, sa_d, sb2_d, & 1678 sb_d, sqb_d, x11_d, x22_d, x33_d 1679 REAL(KIND=dp), DIMENSION(3, 3) :: p 1680 REAL(KIND=dp), DIMENSION(3, 3, 3) :: p_d 1681 REAL(KIND=dp), DIMENSION(3, 5, 5) :: d_d 1682 REAL(KIND=dp), DIMENSION(5, 5) :: d 1683 1684 CPASSERT(ASSOCIATED(ij_matrix)) 1685 IF (PRESENT(do_invert)) do_invert = .FALSE. 1686 IF ((sepi%natorb > 1) .OR. (sepj%natorb > 1)) THEN 1687 ! Compute Geomtric data and interatomic distance 1688 ! CA = COS(PHI) , SA = SIN(PHI) 1689 ! CB = COS(THETA) , SB = SIN(THETA) 1690 ! C2A = COS(2*PHI) , S2A = SIN(2*PHI) 1691 ! C2B = COS(2*THETA), S2B = SIN(2*PHI) 1692 eval = .TRUE. 1693 dbg_inv = .FALSE. 1694 pt5sq3 = 0.5_dp*SQRT(3.0_dp) 1695 imap(1) = 1 1696 imap(2) = 2 1697 imap(3) = 3 1698 check = rjiv(3)/r 1699 IF (PRESENT(debug_invert)) dbg_inv = debug_invert 1700 ! Check for special case: both atoms lie on the Z Axis 1701 IF (ABS(check) > 0.99999999_dp) THEN 1702 IF (PRESENT(do_invert)) THEN 1703 imap(1) = 3 1704 imap(2) = 2 1705 imap(3) = 1 1706 eval = .TRUE. 1707 do_invert = .TRUE. 1708 ELSE 1709 sa = 0.0_dp 1710 sb = 0.0_dp 1711 IF (check < 0.0_dp) THEN 1712 ca = -1.0_dp 1713 cb = -1.0_dp 1714 ELSE IF (check > 0.0_dp) THEN 1715 ca = 1.0_dp 1716 cb = 1.0_dp 1717 ELSE 1718 ca = 0.0_dp 1719 cb = 0.0_dp 1720 END IF 1721 eval = .FALSE. 1722 END IF 1723 END IF 1724 IF (dbg_inv) THEN 1725 ! When debugging the derivatives of the rotation matrix we must possibly force 1726 ! the inversion of the reference frame 1727 CPASSERT(.NOT. do_derivatives) 1728 imap(1) = 3 1729 imap(2) = 2 1730 imap(3) = 1 1731 eval = .TRUE. 1732 END IF 1733 IF (eval) THEN 1734 x11 = rjiv(imap(1)) 1735 x22 = rjiv(imap(2)) 1736 x33 = rjiv(imap(3)) 1737 cb = x33/r 1738 b = x11**2 + x22**2 1739 sqb = SQRT(b) 1740 ca = x11/sqb 1741 sa = x22/sqb 1742 sb = sqb/r 1743 END IF 1744 ! Calculate rotation matrix elements 1745 p(1, 1) = ca*sb 1746 p(2, 1) = ca*cb 1747 p(3, 1) = -sa 1748 p(1, 2) = sa*sb 1749 p(2, 2) = sa*cb 1750 p(3, 2) = ca 1751 p(1, 3) = cb 1752 p(2, 3) = -sb 1753 p(3, 3) = 0.0_dp 1754 IF (sepi%dorb .OR. sepj%dorb) THEN 1755 ca2 = ca**2 1756 cb2 = cb**2 1757 sb2 = sb**2 1758 c2a = 2.0_dp*ca2 - 1.0_dp 1759 c2b = 2.0_dp*cb2 - 1.0_dp 1760 s2a = 2.0_dp*sa*ca 1761 s2b = 2.0_dp*sb*cb 1762 d(1, 1) = pt5sq3*c2a*sb2 1763 d(2, 1) = 0.5_dp*c2a*s2b 1764 d(3, 1) = -s2a*sb 1765 d(4, 1) = c2a*(cb2 + 0.5_dp*sb2) 1766 d(5, 1) = -s2a*cb 1767 d(1, 2) = pt5sq3*ca*s2b 1768 d(2, 2) = ca*c2b 1769 d(3, 2) = -sa*cb 1770 d(4, 2) = -0.5_dp*ca*s2b 1771 d(5, 2) = sa*sb 1772 d(1, 3) = cb2 - 0.5_dp*sb2 1773 d(2, 3) = -pt5sq3*s2b 1774 d(3, 3) = 0.0_dp 1775 d(4, 3) = pt5sq3*sb2 1776 d(5, 3) = 0.0_dp 1777 d(1, 4) = pt5sq3*sa*s2b 1778 d(2, 4) = sa*c2b 1779 d(3, 4) = ca*cb 1780 d(4, 4) = -0.5_dp*sa*s2b 1781 d(5, 4) = -ca*sb 1782 d(1, 5) = pt5sq3*s2a*sb2 1783 d(2, 5) = 0.5_dp*s2a*s2b 1784 d(3, 5) = c2a*sb 1785 d(4, 5) = s2a*(cb2 + 0.5_dp*sb2) 1786 d(5, 5) = c2a*cb 1787 END IF 1788 ! Rotation Elements for : S-P 1789 DO k = 1, 3 1790 DO l = 1, 3 1791 ij_matrix%sp(k, l) = p(k, l) 1792 END DO 1793 END DO 1794 ! Rotation Elements for : P-P 1795 DO k = 1, 3 1796 ij_matrix%pp(1, k, k) = p(k, 1)*p(k, 1) 1797 ij_matrix%pp(2, k, k) = p(k, 2)*p(k, 2) 1798 ij_matrix%pp(3, k, k) = p(k, 3)*p(k, 3) 1799 ij_matrix%pp(4, k, k) = p(k, 1)*p(k, 2) 1800 ij_matrix%pp(5, k, k) = p(k, 1)*p(k, 3) 1801 ij_matrix%pp(6, k, k) = p(k, 2)*p(k, 3) 1802 IF (k /= 1) THEN 1803 DO l = 1, k - 1 1804 ij_matrix%pp(1, k, l) = 2.0_dp*p(k, 1)*p(l, 1) 1805 ij_matrix%pp(2, k, l) = 2.0_dp*p(k, 2)*p(l, 2) 1806 ij_matrix%pp(3, k, l) = 2.0_dp*p(k, 3)*p(l, 3) 1807 ij_matrix%pp(4, k, l) = p(k, 1)*p(l, 2) + p(k, 2)*p(l, 1) 1808 ij_matrix%pp(5, k, l) = p(k, 1)*p(l, 3) + p(k, 3)*p(l, 1) 1809 ij_matrix%pp(6, k, l) = p(k, 2)*p(l, 3) + p(k, 3)*p(l, 2) 1810 END DO 1811 END IF 1812 END DO 1813 IF (sepi%dorb .OR. sepj%dorb) THEN 1814 ! Rotation Elements for : S-D 1815 DO k = 1, 5 1816 DO l = 1, 5 1817 ij_matrix%sd(k, l) = d(k, l) 1818 END DO 1819 END DO 1820 ! Rotation Elements for : D-P 1821 DO k = 1, 5 1822 DO l = 1, 3 1823 ij_matrix%pd(1, k, l) = d(k, 1)*p(l, 1) 1824 ij_matrix%pd(2, k, l) = d(k, 1)*p(l, 2) 1825 ij_matrix%pd(3, k, l) = d(k, 1)*p(l, 3) 1826 ij_matrix%pd(4, k, l) = d(k, 2)*p(l, 1) 1827 ij_matrix%pd(5, k, l) = d(k, 2)*p(l, 2) 1828 ij_matrix%pd(6, k, l) = d(k, 2)*p(l, 3) 1829 ij_matrix%pd(7, k, l) = d(k, 3)*p(l, 1) 1830 ij_matrix%pd(8, k, l) = d(k, 3)*p(l, 2) 1831 ij_matrix%pd(9, k, l) = d(k, 3)*p(l, 3) 1832 ij_matrix%pd(10, k, l) = d(k, 4)*p(l, 1) 1833 ij_matrix%pd(11, k, l) = d(k, 4)*p(l, 2) 1834 ij_matrix%pd(12, k, l) = d(k, 4)*p(l, 3) 1835 ij_matrix%pd(13, k, l) = d(k, 5)*p(l, 1) 1836 ij_matrix%pd(14, k, l) = d(k, 5)*p(l, 2) 1837 ij_matrix%pd(15, k, l) = d(k, 5)*p(l, 3) 1838 END DO 1839 END DO 1840 ! Rotation Elements for : D-D 1841 DO k = 1, 5 1842 ij_matrix%dd(1, k, k) = d(k, 1)*d(k, 1) 1843 ij_matrix%dd(2, k, k) = d(k, 2)*d(k, 2) 1844 ij_matrix%dd(3, k, k) = d(k, 3)*d(k, 3) 1845 ij_matrix%dd(4, k, k) = d(k, 4)*d(k, 4) 1846 ij_matrix%dd(5, k, k) = d(k, 5)*d(k, 5) 1847 ij_matrix%dd(6, k, k) = d(k, 1)*d(k, 2) 1848 ij_matrix%dd(7, k, k) = d(k, 1)*d(k, 3) 1849 ij_matrix%dd(8, k, k) = d(k, 2)*d(k, 3) 1850 ij_matrix%dd(9, k, k) = d(k, 1)*d(k, 4) 1851 ij_matrix%dd(10, k, k) = d(k, 2)*d(k, 4) 1852 ij_matrix%dd(11, k, k) = d(k, 3)*d(k, 4) 1853 ij_matrix%dd(12, k, k) = d(k, 1)*d(k, 5) 1854 ij_matrix%dd(13, k, k) = d(k, 2)*d(k, 5) 1855 ij_matrix%dd(14, k, k) = d(k, 3)*d(k, 5) 1856 ij_matrix%dd(15, k, k) = d(k, 4)*d(k, 5) 1857 IF (k /= 1) THEN 1858 DO l = 1, k - 1 1859 ij_matrix%dd(1, k, l) = 2.0_dp*d(k, 1)*d(l, 1) 1860 ij_matrix%dd(2, k, l) = 2.0_dp*d(k, 2)*d(l, 2) 1861 ij_matrix%dd(3, k, l) = 2.0_dp*d(k, 3)*d(l, 3) 1862 ij_matrix%dd(4, k, l) = 2.0_dp*d(k, 4)*d(l, 4) 1863 ij_matrix%dd(5, k, l) = 2.0_dp*d(k, 5)*d(l, 5) 1864 ij_matrix%dd(6, k, l) = d(k, 1)*d(l, 2) + d(k, 2)*d(l, 1) 1865 ij_matrix%dd(7, k, l) = d(k, 1)*d(l, 3) + d(k, 3)*d(l, 1) 1866 ij_matrix%dd(8, k, l) = d(k, 2)*d(l, 3) + d(k, 3)*d(l, 2) 1867 ij_matrix%dd(9, k, l) = d(k, 1)*d(l, 4) + d(k, 4)*d(l, 1) 1868 ij_matrix%dd(10, k, l) = d(k, 2)*d(l, 4) + d(k, 4)*d(l, 2) 1869 ij_matrix%dd(11, k, l) = d(k, 3)*d(l, 4) + d(k, 4)*d(l, 3) 1870 ij_matrix%dd(12, k, l) = d(k, 1)*d(l, 5) + d(k, 5)*d(l, 1) 1871 ij_matrix%dd(13, k, l) = d(k, 2)*d(l, 5) + d(k, 5)*d(l, 2) 1872 ij_matrix%dd(14, k, l) = d(k, 3)*d(l, 5) + d(k, 5)*d(l, 3) 1873 ij_matrix%dd(15, k, l) = d(k, 4)*d(l, 5) + d(k, 5)*d(l, 4) 1874 END DO 1875 END IF 1876 END DO 1877 END IF 1878 ! Evaluate Derivatives if required 1879 IF (do_derivatives) THEN 1880 ! This condition is necessary because derivatives uses the invertion of the 1881 ! axis for treating the divergence point along z-axis 1882 CPASSERT(eval) 1883 x11_d = 0.0_dp; x11_d(1) = 1.0_dp 1884 x22_d = 0.0_dp; x22_d(2) = 1.0_dp 1885 x33_d = 0.0_dp; x33_d(3) = 1.0_dp 1886 r_d = (/x11, x22, x33/)/r 1887 b_d = 2.0_dp*(x11*x11_d + x22*x22_d) 1888 sqb_d = (0.5_dp/sqb)*b_d 1889 r2i = 1.0_dp/r**2 1890 sqb2i = 1.0_dp/sqb**2 1891 cb_d = (x33_d*r - x33*r_d)*r2i 1892 ca_d = (x11_d*sqb - x11*sqb_d)*sqb2i 1893 sa_d = (x22_d*sqb - x22*sqb_d)*sqb2i 1894 sb_d = (sqb_d*r - sqb*r_d)*r2i 1895 ! Calculate derivatives of rotation matrix elements 1896 p_d(:, 1, 1) = ca_d*sb + ca*sb_d 1897 p_d(:, 2, 1) = ca_d*cb + ca*cb_d 1898 p_d(:, 3, 1) = -sa_d 1899 p_d(:, 1, 2) = sa_d*sb + sa*sb_d 1900 p_d(:, 2, 2) = sa_d*cb + sa*cb_d 1901 p_d(:, 3, 2) = ca_d 1902 p_d(:, 1, 3) = cb_d 1903 p_d(:, 2, 3) = -sb_d 1904 p_d(:, 3, 3) = 0.0_dp 1905 IF (sepi%dorb .OR. sepj%dorb) THEN 1906 ca2_d = 2.0_dp*ca*ca_d 1907 cb2_d = 2.0_dp*cb*cb_d 1908 sb2_d = 2.0_dp*sb*sb_d 1909 c2a_d = 2.0_dp*ca2_d 1910 c2b_d = 2.0_dp*cb2_d 1911 s2a_d = 2.0_dp*(sa_d*ca + sa*ca_d) 1912 s2b_d = 2.0_dp*(sb_d*cb + sb*cb_d) 1913 d_d(:, 1, 1) = pt5sq3*(c2a_d*sb2 + c2a*sb2_d) 1914 d_d(:, 2, 1) = 0.5_dp*(c2a_d*s2b + c2a*s2b_d) 1915 d_d(:, 3, 1) = -s2a_d*sb - s2a*sb_d 1916 d_d(:, 4, 1) = c2a_d*(cb2 + 0.5_dp*sb2) + c2a*(cb2_d + 0.5_dp*sb2_d) 1917 d_d(:, 5, 1) = -s2a_d*cb - s2a*cb_d 1918 d_d(:, 1, 2) = pt5sq3*(ca_d*s2b + ca*s2b_d) 1919 d_d(:, 2, 2) = ca_d*c2b + ca*c2b_d 1920 d_d(:, 3, 2) = -sa_d*cb - sa*cb_d 1921 d_d(:, 4, 2) = -0.5_dp*(ca_d*s2b + ca*s2b_d) 1922 d_d(:, 5, 2) = sa_d*sb + sa*sb_d 1923 d_d(:, 1, 3) = cb2_d - 0.5_dp*sb2_d 1924 d_d(:, 2, 3) = -pt5sq3*s2b_d 1925 d_d(:, 3, 3) = 0.0_dp 1926 d_d(:, 4, 3) = pt5sq3*sb2_d 1927 d_d(:, 5, 3) = 0.0_dp 1928 d_d(:, 1, 4) = pt5sq3*(sa_d*s2b + sa*s2b_d) 1929 d_d(:, 2, 4) = sa_d*c2b + sa*c2b_d 1930 d_d(:, 3, 4) = ca_d*cb + ca*cb_d 1931 d_d(:, 4, 4) = -0.5_dp*(sa_d*s2b + sa*s2b_d) 1932 d_d(:, 5, 4) = -ca_d*sb - ca*sb_d 1933 d_d(:, 1, 5) = pt5sq3*(s2a_d*sb2 + s2a*sb2_d) 1934 d_d(:, 2, 5) = 0.5_dp*(s2a_d*s2b + s2a*s2b_d) 1935 d_d(:, 3, 5) = c2a_d*sb + c2a*sb_d 1936 d_d(:, 4, 5) = s2a_d*(cb2 + 0.5_dp*sb2) + s2a*(cb2_d + 0.5_dp*sb2_d) 1937 d_d(:, 5, 5) = c2a_d*cb + c2a*cb_d 1938 END IF 1939 ! Derivative for Rotation Elements for : S-P 1940 DO k = 1, 3 1941 DO l = 1, 3 1942 ij_matrix%sp_d(:, k, l) = p_d(:, k, l) 1943 END DO 1944 END DO 1945 ! Derivative for Rotation Elements for : P-P 1946 DO k = 1, 3 1947 ij_matrix%pp_d(:, 1, k, k) = p_d(:, k, 1)*p(k, 1) + p(k, 1)*p_d(:, k, 1) 1948 ij_matrix%pp_d(:, 2, k, k) = p_d(:, k, 2)*p(k, 2) + p(k, 2)*p_d(:, k, 2) 1949 ij_matrix%pp_d(:, 3, k, k) = p_d(:, k, 3)*p(k, 3) + p(k, 3)*p_d(:, k, 3) 1950 ij_matrix%pp_d(:, 4, k, k) = p_d(:, k, 1)*p(k, 2) + p(k, 1)*p_d(:, k, 2) 1951 ij_matrix%pp_d(:, 5, k, k) = p_d(:, k, 1)*p(k, 3) + p(k, 1)*p_d(:, k, 3) 1952 ij_matrix%pp_d(:, 6, k, k) = p_d(:, k, 2)*p(k, 3) + p(k, 2)*p_d(:, k, 3) 1953 IF (k /= 1) THEN 1954 DO l = 1, k - 1 1955 ij_matrix%pp_d(:, 1, k, l) = 2.0_dp*(p_d(:, k, 1)*p(l, 1) + p(k, 1)*p_d(:, l, 1)) 1956 ij_matrix%pp_d(:, 2, k, l) = 2.0_dp*(p_d(:, k, 2)*p(l, 2) + p(k, 2)*p_d(:, l, 2)) 1957 ij_matrix%pp_d(:, 3, k, l) = 2.0_dp*(p_d(:, k, 3)*p(l, 3) + p(k, 3)*p_d(:, l, 3)) 1958 ij_matrix%pp_d(:, 4, k, l) = (p_d(:, k, 1)*p(l, 2) + p(k, 1)*p_d(:, l, 2)) + & 1959 (p_d(:, k, 2)*p(l, 1) + p(k, 2)*p_d(:, l, 1)) 1960 ij_matrix%pp_d(:, 5, k, l) = (p_d(:, k, 1)*p(l, 3) + p(k, 1)*p_d(:, l, 3)) + & 1961 (p_d(:, k, 3)*p(l, 1) + p(k, 3)*p_d(:, l, 1)) 1962 ij_matrix%pp_d(:, 6, k, l) = (p_d(:, k, 2)*p(l, 3) + p(k, 2)*p_d(:, l, 3)) + & 1963 (p_d(:, k, 3)*p(l, 2) + p(k, 3)*p_d(:, l, 2)) 1964 END DO 1965 END IF 1966 END DO 1967 IF (sepi%dorb .OR. sepj%dorb) THEN 1968 ! Rotation Elements for : S-D 1969 DO k = 1, 5 1970 DO l = 1, 5 1971 ij_matrix%sd_d(:, k, l) = d_d(:, k, l) 1972 END DO 1973 END DO 1974 ! Rotation Elements for : D-P 1975 DO k = 1, 5 1976 DO l = 1, 3 1977 ij_matrix%pd_d(:, 1, k, l) = (d_d(:, k, 1)*p(l, 1) + d(k, 1)*p_d(:, l, 1)) 1978 ij_matrix%pd_d(:, 2, k, l) = (d_d(:, k, 1)*p(l, 2) + d(k, 1)*p_d(:, l, 2)) 1979 ij_matrix%pd_d(:, 3, k, l) = (d_d(:, k, 1)*p(l, 3) + d(k, 1)*p_d(:, l, 3)) 1980 ij_matrix%pd_d(:, 4, k, l) = (d_d(:, k, 2)*p(l, 1) + d(k, 2)*p_d(:, l, 1)) 1981 ij_matrix%pd_d(:, 5, k, l) = (d_d(:, k, 2)*p(l, 2) + d(k, 2)*p_d(:, l, 2)) 1982 ij_matrix%pd_d(:, 6, k, l) = (d_d(:, k, 2)*p(l, 3) + d(k, 2)*p_d(:, l, 3)) 1983 ij_matrix%pd_d(:, 7, k, l) = (d_d(:, k, 3)*p(l, 1) + d(k, 3)*p_d(:, l, 1)) 1984 ij_matrix%pd_d(:, 8, k, l) = (d_d(:, k, 3)*p(l, 2) + d(k, 3)*p_d(:, l, 2)) 1985 ij_matrix%pd_d(:, 9, k, l) = (d_d(:, k, 3)*p(l, 3) + d(k, 3)*p_d(:, l, 3)) 1986 ij_matrix%pd_d(:, 10, k, l) = (d_d(:, k, 4)*p(l, 1) + d(k, 4)*p_d(:, l, 1)) 1987 ij_matrix%pd_d(:, 11, k, l) = (d_d(:, k, 4)*p(l, 2) + d(k, 4)*p_d(:, l, 2)) 1988 ij_matrix%pd_d(:, 12, k, l) = (d_d(:, k, 4)*p(l, 3) + d(k, 4)*p_d(:, l, 3)) 1989 ij_matrix%pd_d(:, 13, k, l) = (d_d(:, k, 5)*p(l, 1) + d(k, 5)*p_d(:, l, 1)) 1990 ij_matrix%pd_d(:, 14, k, l) = (d_d(:, k, 5)*p(l, 2) + d(k, 5)*p_d(:, l, 2)) 1991 ij_matrix%pd_d(:, 15, k, l) = (d_d(:, k, 5)*p(l, 3) + d(k, 5)*p_d(:, l, 3)) 1992 END DO 1993 END DO 1994 ! Rotation Elements for : D-D 1995 DO k = 1, 5 1996 ij_matrix%dd_d(:, 1, k, k) = (d_d(:, k, 1)*d(k, 1) + d(k, 1)*d_d(:, k, 1)) 1997 ij_matrix%dd_d(:, 2, k, k) = (d_d(:, k, 2)*d(k, 2) + d(k, 2)*d_d(:, k, 2)) 1998 ij_matrix%dd_d(:, 3, k, k) = (d_d(:, k, 3)*d(k, 3) + d(k, 3)*d_d(:, k, 3)) 1999 ij_matrix%dd_d(:, 4, k, k) = (d_d(:, k, 4)*d(k, 4) + d(k, 4)*d_d(:, k, 4)) 2000 ij_matrix%dd_d(:, 5, k, k) = (d_d(:, k, 5)*d(k, 5) + d(k, 5)*d_d(:, k, 5)) 2001 ij_matrix%dd_d(:, 6, k, k) = (d_d(:, k, 1)*d(k, 2) + d(k, 1)*d_d(:, k, 2)) 2002 ij_matrix%dd_d(:, 7, k, k) = (d_d(:, k, 1)*d(k, 3) + d(k, 1)*d_d(:, k, 3)) 2003 ij_matrix%dd_d(:, 8, k, k) = (d_d(:, k, 2)*d(k, 3) + d(k, 2)*d_d(:, k, 3)) 2004 ij_matrix%dd_d(:, 9, k, k) = (d_d(:, k, 1)*d(k, 4) + d(k, 1)*d_d(:, k, 4)) 2005 ij_matrix%dd_d(:, 10, k, k) = (d_d(:, k, 2)*d(k, 4) + d(k, 2)*d_d(:, k, 4)) 2006 ij_matrix%dd_d(:, 11, k, k) = (d_d(:, k, 3)*d(k, 4) + d(k, 3)*d_d(:, k, 4)) 2007 ij_matrix%dd_d(:, 12, k, k) = (d_d(:, k, 1)*d(k, 5) + d(k, 1)*d_d(:, k, 5)) 2008 ij_matrix%dd_d(:, 13, k, k) = (d_d(:, k, 2)*d(k, 5) + d(k, 2)*d_d(:, k, 5)) 2009 ij_matrix%dd_d(:, 14, k, k) = (d_d(:, k, 3)*d(k, 5) + d(k, 3)*d_d(:, k, 5)) 2010 ij_matrix%dd_d(:, 15, k, k) = (d_d(:, k, 4)*d(k, 5) + d(k, 4)*d_d(:, k, 5)) 2011 IF (k /= 1) THEN 2012 DO l = 1, k - 1 2013 ij_matrix%dd_d(:, 1, k, l) = 2.0_dp*(d_d(:, k, 1)*d(l, 1) + d(k, 1)*d_d(:, l, 1)) 2014 ij_matrix%dd_d(:, 2, k, l) = 2.0_dp*(d_d(:, k, 2)*d(l, 2) + d(k, 2)*d_d(:, l, 2)) 2015 ij_matrix%dd_d(:, 3, k, l) = 2.0_dp*(d_d(:, k, 3)*d(l, 3) + d(k, 3)*d_d(:, l, 3)) 2016 ij_matrix%dd_d(:, 4, k, l) = 2.0_dp*(d_d(:, k, 4)*d(l, 4) + d(k, 4)*d_d(:, l, 4)) 2017 ij_matrix%dd_d(:, 5, k, l) = 2.0_dp*(d_d(:, k, 5)*d(l, 5) + d(k, 5)*d_d(:, l, 5)) 2018 ij_matrix%dd_d(:, 6, k, l) = (d_d(:, k, 1)*d(l, 2) + d(k, 1)*d_d(:, l, 2)) + & 2019 (d_d(:, k, 2)*d(l, 1) + d(k, 2)*d_d(:, l, 1)) 2020 ij_matrix%dd_d(:, 7, k, l) = (d_d(:, k, 1)*d(l, 3) + d(k, 1)*d_d(:, l, 3)) + & 2021 (d_d(:, k, 3)*d(l, 1) + d(k, 3)*d_d(:, l, 1)) 2022 ij_matrix%dd_d(:, 8, k, l) = (d_d(:, k, 2)*d(l, 3) + d(k, 2)*d_d(:, l, 3)) + & 2023 (d_d(:, k, 3)*d(l, 2) + d(k, 3)*d_d(:, l, 2)) 2024 ij_matrix%dd_d(:, 9, k, l) = (d_d(:, k, 1)*d(l, 4) + d(k, 1)*d_d(:, l, 4)) + & 2025 (d_d(:, k, 4)*d(l, 1) + d(k, 4)*d_d(:, l, 1)) 2026 ij_matrix%dd_d(:, 10, k, l) = (d_d(:, k, 2)*d(l, 4) + d(k, 2)*d_d(:, l, 4)) + & 2027 (d_d(:, k, 4)*d(l, 2) + d(k, 4)*d_d(:, l, 2)) 2028 ij_matrix%dd_d(:, 11, k, l) = (d_d(:, k, 3)*d(l, 4) + d(k, 3)*d_d(:, l, 4)) + & 2029 (d_d(:, k, 4)*d(l, 3) + d(k, 4)*d_d(:, l, 3)) 2030 ij_matrix%dd_d(:, 12, k, l) = (d_d(:, k, 1)*d(l, 5) + d(k, 1)*d_d(:, l, 5)) + & 2031 (d_d(:, k, 5)*d(l, 1) + d(k, 5)*d_d(:, l, 1)) 2032 ij_matrix%dd_d(:, 13, k, l) = (d_d(:, k, 2)*d(l, 5) + d(k, 2)*d_d(:, l, 5)) + & 2033 (d_d(:, k, 5)*d(l, 2) + d(k, 5)*d_d(:, l, 2)) 2034 ij_matrix%dd_d(:, 14, k, l) = (d_d(:, k, 3)*d(l, 5) + d(k, 3)*d_d(:, l, 5)) + & 2035 (d_d(:, k, 5)*d(l, 3) + d(k, 5)*d_d(:, l, 3)) 2036 ij_matrix%dd_d(:, 15, k, l) = (d_d(:, k, 4)*d(l, 5) + d(k, 4)*d_d(:, l, 5)) + & 2037 (d_d(:, k, 5)*d(l, 4) + d(k, 5)*d_d(:, l, 4)) 2038 END DO 2039 END IF 2040 END DO 2041 END IF 2042 IF (debug_this_module) THEN 2043 CALL check_rotmat_der(sepi, sepj, rjiv, ij_matrix, do_invert=do_invert) 2044 END IF 2045 END IF 2046 END IF 2047 END SUBROUTINE rotmat 2048 2049! ************************************************************************************************** 2050!> \brief First Step of the rotation of the two-electron two-center integrals in 2051!> SPD basis 2052!> 2053!> \param sepi ... 2054!> \param sepj ... 2055!> \param rijv ... 2056!> \param se_int_control ... 2057!> \param se_taper ... 2058!> \param invert ... 2059!> \param ii ... 2060!> \param kk ... 2061!> \param rep ... 2062!> \param logv ... 2063!> \param ij_matrix ... 2064!> \param v ... 2065!> \param lgrad ... 2066!> \param rep_d ... 2067!> \param v_d ... 2068!> \param logv_d ... 2069!> \param drij ... 2070!> \date 04.2008 [tlaino] 2071!> \author Teodoro Laino [tlaino] - University of Zurich 2072! ************************************************************************************************** 2073 RECURSIVE SUBROUTINE rot_2el_2c_first(sepi, sepj, rijv, se_int_control, se_taper, & 2074 invert, ii, kk, rep, logv, ij_matrix, v, lgrad, rep_d, v_d, logv_d, drij) 2075 TYPE(semi_empirical_type), POINTER :: sepi, sepj 2076 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rijv 2077 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 2078 TYPE(se_taper_type), POINTER :: se_taper 2079 LOGICAL, INTENT(IN) :: invert 2080 INTEGER, INTENT(IN) :: ii, kk 2081 REAL(KIND=dp), DIMENSION(491), INTENT(IN) :: rep 2082 LOGICAL, DIMENSION(45, 45), INTENT(OUT) :: logv 2083 TYPE(rotmat_type), POINTER :: ij_matrix 2084 REAL(KIND=dp), DIMENSION(45, 45), INTENT(OUT) :: v 2085 LOGICAL, INTENT(IN) :: lgrad 2086 REAL(KIND=dp), DIMENSION(491), INTENT(IN), & 2087 OPTIONAL :: rep_d 2088 REAL(KIND=dp), DIMENSION(3, 45, 45), INTENT(OUT), & 2089 OPTIONAL :: v_d 2090 LOGICAL, DIMENSION(45, 45), INTENT(OUT), OPTIONAL :: logv_d 2091 REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: drij 2092 2093 CHARACTER(len=*), PARAMETER :: routineN = 'rot_2el_2c_first', & 2094 routineP = moduleN//':'//routineN 2095 2096 INTEGER :: i, i1, ij, j1, k, k1, kl, l, l1, limkl, & 2097 ll, mm, nd 2098 REAL(KIND=dp) :: wrepp, wrepp_d(3) 2099 2100 IF (lgrad) THEN 2101 CPASSERT(PRESENT(rep_d)) 2102 CPASSERT(PRESENT(v_d)) 2103 CPASSERT(PRESENT(logv_d)) 2104 CPASSERT(PRESENT(drij)) 2105 END IF 2106 limkl = indexb(kk, kk) 2107 DO k = 1, limkl 2108 DO i = 1, 45 2109 logv(i, k) = .FALSE. 2110 v(i, k) = 0.0_dp 2111 END DO 2112 END DO 2113 ! 2114 DO i1 = 1, ii 2115 DO j1 = 1, i1 2116 ij = indexa(i1, j1) 2117 ! 2118 DO k1 = 1, kk 2119 ! 2120 DO l1 = 1, k1 2121 kl = indexa(k1, l1) 2122 nd = ijkl_ind(ij, kl) 2123 IF (nd /= 0) THEN 2124 ! 2125 wrepp = rep(nd) 2126 ll = indexb(k1, l1) 2127 mm = int2c_type(ll) 2128 ! 2129 IF (mm == 1) THEN 2130 v(ij, 1) = wrepp 2131 ELSE IF (mm == 2) THEN 2132 k = k1 - 1 2133 v(ij, 2) = v(ij, 2) + ij_matrix%sp(k, 1)*wrepp 2134 v(ij, 4) = v(ij, 4) + ij_matrix%sp(k, 2)*wrepp 2135 v(ij, 7) = v(ij, 7) + ij_matrix%sp(k, 3)*wrepp 2136 ELSE IF (mm == 3) THEN 2137 k = k1 - 1 2138 l = l1 - 1 2139 v(ij, 3) = v(ij, 3) + ij_matrix%pp(1, k, l)*wrepp 2140 v(ij, 6) = v(ij, 6) + ij_matrix%pp(2, k, l)*wrepp 2141 v(ij, 10) = v(ij, 10) + ij_matrix%pp(3, k, l)*wrepp 2142 v(ij, 5) = v(ij, 5) + ij_matrix%pp(4, k, l)*wrepp 2143 v(ij, 8) = v(ij, 8) + ij_matrix%pp(5, k, l)*wrepp 2144 v(ij, 9) = v(ij, 9) + ij_matrix%pp(6, k, l)*wrepp 2145 ELSE IF (mm == 4) THEN 2146 k = k1 - 4 2147 v(ij, 11) = v(ij, 11) + ij_matrix%sd(k, 1)*wrepp 2148 v(ij, 16) = v(ij, 16) + ij_matrix%sd(k, 2)*wrepp 2149 v(ij, 22) = v(ij, 22) + ij_matrix%sd(k, 3)*wrepp 2150 v(ij, 29) = v(ij, 29) + ij_matrix%sd(k, 4)*wrepp 2151 v(ij, 37) = v(ij, 37) + ij_matrix%sd(k, 5)*wrepp 2152 ELSE IF (mm == 5) THEN 2153 k = k1 - 4 2154 l = l1 - 1 2155 v(ij, 12) = v(ij, 12) + ij_matrix%pd(1, k, l)*wrepp 2156 v(ij, 13) = v(ij, 13) + ij_matrix%pd(2, k, l)*wrepp 2157 v(ij, 14) = v(ij, 14) + ij_matrix%pd(3, k, l)*wrepp 2158 v(ij, 17) = v(ij, 17) + ij_matrix%pd(4, k, l)*wrepp 2159 v(ij, 18) = v(ij, 18) + ij_matrix%pd(5, k, l)*wrepp 2160 v(ij, 19) = v(ij, 19) + ij_matrix%pd(6, k, l)*wrepp 2161 v(ij, 23) = v(ij, 23) + ij_matrix%pd(7, k, l)*wrepp 2162 v(ij, 24) = v(ij, 24) + ij_matrix%pd(8, k, l)*wrepp 2163 v(ij, 25) = v(ij, 25) + ij_matrix%pd(9, k, l)*wrepp 2164 v(ij, 30) = v(ij, 30) + ij_matrix%pd(10, k, l)*wrepp 2165 v(ij, 31) = v(ij, 31) + ij_matrix%pd(11, k, l)*wrepp 2166 v(ij, 32) = v(ij, 32) + ij_matrix%pd(12, k, l)*wrepp 2167 v(ij, 38) = v(ij, 38) + ij_matrix%pd(13, k, l)*wrepp 2168 v(ij, 39) = v(ij, 39) + ij_matrix%pd(14, k, l)*wrepp 2169 v(ij, 40) = v(ij, 40) + ij_matrix%pd(15, k, l)*wrepp 2170 ELSE IF (mm == 6) THEN 2171 k = k1 - 4 2172 l = l1 - 4 2173 v(ij, 15) = v(ij, 15) + ij_matrix%dd(1, k, l)*wrepp 2174 v(ij, 21) = v(ij, 21) + ij_matrix%dd(2, k, l)*wrepp 2175 v(ij, 28) = v(ij, 28) + ij_matrix%dd(3, k, l)*wrepp 2176 v(ij, 36) = v(ij, 36) + ij_matrix%dd(4, k, l)*wrepp 2177 v(ij, 45) = v(ij, 45) + ij_matrix%dd(5, k, l)*wrepp 2178 v(ij, 20) = v(ij, 20) + ij_matrix%dd(6, k, l)*wrepp 2179 v(ij, 26) = v(ij, 26) + ij_matrix%dd(7, k, l)*wrepp 2180 v(ij, 27) = v(ij, 27) + ij_matrix%dd(8, k, l)*wrepp 2181 v(ij, 33) = v(ij, 33) + ij_matrix%dd(9, k, l)*wrepp 2182 v(ij, 34) = v(ij, 34) + ij_matrix%dd(10, k, l)*wrepp 2183 v(ij, 35) = v(ij, 35) + ij_matrix%dd(11, k, l)*wrepp 2184 v(ij, 41) = v(ij, 41) + ij_matrix%dd(12, k, l)*wrepp 2185 v(ij, 42) = v(ij, 42) + ij_matrix%dd(13, k, l)*wrepp 2186 v(ij, 43) = v(ij, 43) + ij_matrix%dd(14, k, l)*wrepp 2187 v(ij, 44) = v(ij, 44) + ij_matrix%dd(15, k, l)*wrepp 2188 END IF 2189 END IF 2190 END DO 2191 END DO 2192 DO kl = 1, limkl 2193 logv(ij, kl) = (ABS(v(ij, kl)) > 0.0_dp) 2194 END DO 2195 END DO 2196 END DO 2197 ! Gradients 2198 IF (lgrad) THEN 2199 DO k = 1, limkl 2200 DO i = 1, 45 2201 logv_d(i, k) = .FALSE. 2202 v_d(1, i, k) = 0.0_dp 2203 v_d(2, i, k) = 0.0_dp 2204 v_d(3, i, k) = 0.0_dp 2205 END DO 2206 END DO 2207 DO i1 = 1, ii 2208 DO j1 = 1, i1 2209 ij = indexa(i1, j1) 2210 ! 2211 DO k1 = 1, kk 2212 ! 2213 DO l1 = 1, k1 2214 kl = indexa(k1, l1) 2215 nd = ijkl_ind(ij, kl) 2216 IF (nd /= 0) THEN 2217 ! 2218 wrepp = rep(nd) 2219 wrepp_d = rep_d(nd)*drij 2220 ll = indexb(k1, l1) 2221 mm = int2c_type(ll) 2222 ! 2223 IF (mm == 1) THEN 2224 v_d(1, ij, 1) = wrepp_d(1) 2225 v_d(2, ij, 1) = wrepp_d(2) 2226 v_d(3, ij, 1) = wrepp_d(3) 2227 ELSE IF (mm == 2) THEN 2228 k = k1 - 1 2229 v_d(1, ij, 2) = v_d(1, ij, 2) + ij_matrix%sp_d(1, k, 1)*wrepp + ij_matrix%sp(k, 1)*wrepp_d(1) 2230 v_d(1, ij, 4) = v_d(1, ij, 4) + ij_matrix%sp_d(1, k, 2)*wrepp + ij_matrix%sp(k, 2)*wrepp_d(1) 2231 v_d(1, ij, 7) = v_d(1, ij, 7) + ij_matrix%sp_d(1, k, 3)*wrepp + ij_matrix%sp(k, 3)*wrepp_d(1) 2232 2233 v_d(2, ij, 2) = v_d(2, ij, 2) + ij_matrix%sp_d(2, k, 1)*wrepp + ij_matrix%sp(k, 1)*wrepp_d(2) 2234 v_d(2, ij, 4) = v_d(2, ij, 4) + ij_matrix%sp_d(2, k, 2)*wrepp + ij_matrix%sp(k, 2)*wrepp_d(2) 2235 v_d(2, ij, 7) = v_d(2, ij, 7) + ij_matrix%sp_d(2, k, 3)*wrepp + ij_matrix%sp(k, 3)*wrepp_d(2) 2236 2237 v_d(3, ij, 2) = v_d(3, ij, 2) + ij_matrix%sp_d(3, k, 1)*wrepp + ij_matrix%sp(k, 1)*wrepp_d(3) 2238 v_d(3, ij, 4) = v_d(3, ij, 4) + ij_matrix%sp_d(3, k, 2)*wrepp + ij_matrix%sp(k, 2)*wrepp_d(3) 2239 v_d(3, ij, 7) = v_d(3, ij, 7) + ij_matrix%sp_d(3, k, 3)*wrepp + ij_matrix%sp(k, 3)*wrepp_d(3) 2240 ELSE IF (mm == 3) THEN 2241 k = k1 - 1 2242 l = l1 - 1 2243 v_d(1, ij, 3) = v_d(1, ij, 3) + ij_matrix%pp_d(1, 1, k, l)*wrepp + ij_matrix%pp(1, k, l)*wrepp_d(1) 2244 v_d(1, ij, 6) = v_d(1, ij, 6) + ij_matrix%pp_d(1, 2, k, l)*wrepp + ij_matrix%pp(2, k, l)*wrepp_d(1) 2245 v_d(1, ij, 10) = v_d(1, ij, 10) + ij_matrix%pp_d(1, 3, k, l)*wrepp + ij_matrix%pp(3, k, l)*wrepp_d(1) 2246 v_d(1, ij, 5) = v_d(1, ij, 5) + ij_matrix%pp_d(1, 4, k, l)*wrepp + ij_matrix%pp(4, k, l)*wrepp_d(1) 2247 v_d(1, ij, 8) = v_d(1, ij, 8) + ij_matrix%pp_d(1, 5, k, l)*wrepp + ij_matrix%pp(5, k, l)*wrepp_d(1) 2248 v_d(1, ij, 9) = v_d(1, ij, 9) + ij_matrix%pp_d(1, 6, k, l)*wrepp + ij_matrix%pp(6, k, l)*wrepp_d(1) 2249 2250 v_d(2, ij, 3) = v_d(2, ij, 3) + ij_matrix%pp_d(2, 1, k, l)*wrepp + ij_matrix%pp(1, k, l)*wrepp_d(2) 2251 v_d(2, ij, 6) = v_d(2, ij, 6) + ij_matrix%pp_d(2, 2, k, l)*wrepp + ij_matrix%pp(2, k, l)*wrepp_d(2) 2252 v_d(2, ij, 10) = v_d(2, ij, 10) + ij_matrix%pp_d(2, 3, k, l)*wrepp + ij_matrix%pp(3, k, l)*wrepp_d(2) 2253 v_d(2, ij, 5) = v_d(2, ij, 5) + ij_matrix%pp_d(2, 4, k, l)*wrepp + ij_matrix%pp(4, k, l)*wrepp_d(2) 2254 v_d(2, ij, 8) = v_d(2, ij, 8) + ij_matrix%pp_d(2, 5, k, l)*wrepp + ij_matrix%pp(5, k, l)*wrepp_d(2) 2255 v_d(2, ij, 9) = v_d(2, ij, 9) + ij_matrix%pp_d(2, 6, k, l)*wrepp + ij_matrix%pp(6, k, l)*wrepp_d(2) 2256 2257 v_d(3, ij, 3) = v_d(3, ij, 3) + ij_matrix%pp_d(3, 1, k, l)*wrepp + ij_matrix%pp(1, k, l)*wrepp_d(3) 2258 v_d(3, ij, 6) = v_d(3, ij, 6) + ij_matrix%pp_d(3, 2, k, l)*wrepp + ij_matrix%pp(2, k, l)*wrepp_d(3) 2259 v_d(3, ij, 10) = v_d(3, ij, 10) + ij_matrix%pp_d(3, 3, k, l)*wrepp + ij_matrix%pp(3, k, l)*wrepp_d(3) 2260 v_d(3, ij, 5) = v_d(3, ij, 5) + ij_matrix%pp_d(3, 4, k, l)*wrepp + ij_matrix%pp(4, k, l)*wrepp_d(3) 2261 v_d(3, ij, 8) = v_d(3, ij, 8) + ij_matrix%pp_d(3, 5, k, l)*wrepp + ij_matrix%pp(5, k, l)*wrepp_d(3) 2262 v_d(3, ij, 9) = v_d(3, ij, 9) + ij_matrix%pp_d(3, 6, k, l)*wrepp + ij_matrix%pp(6, k, l)*wrepp_d(3) 2263 ELSE IF (mm == 4) THEN 2264 k = k1 - 4 2265 v_d(1, ij, 11) = v_d(1, ij, 11) + ij_matrix%sd_d(1, k, 1)*wrepp + ij_matrix%sd(k, 1)*wrepp_d(1) 2266 v_d(1, ij, 16) = v_d(1, ij, 16) + ij_matrix%sd_d(1, k, 2)*wrepp + ij_matrix%sd(k, 2)*wrepp_d(1) 2267 v_d(1, ij, 22) = v_d(1, ij, 22) + ij_matrix%sd_d(1, k, 3)*wrepp + ij_matrix%sd(k, 3)*wrepp_d(1) 2268 v_d(1, ij, 29) = v_d(1, ij, 29) + ij_matrix%sd_d(1, k, 4)*wrepp + ij_matrix%sd(k, 4)*wrepp_d(1) 2269 v_d(1, ij, 37) = v_d(1, ij, 37) + ij_matrix%sd_d(1, k, 5)*wrepp + ij_matrix%sd(k, 5)*wrepp_d(1) 2270 2271 v_d(2, ij, 11) = v_d(2, ij, 11) + ij_matrix%sd_d(2, k, 1)*wrepp + ij_matrix%sd(k, 1)*wrepp_d(2) 2272 v_d(2, ij, 16) = v_d(2, ij, 16) + ij_matrix%sd_d(2, k, 2)*wrepp + ij_matrix%sd(k, 2)*wrepp_d(2) 2273 v_d(2, ij, 22) = v_d(2, ij, 22) + ij_matrix%sd_d(2, k, 3)*wrepp + ij_matrix%sd(k, 3)*wrepp_d(2) 2274 v_d(2, ij, 29) = v_d(2, ij, 29) + ij_matrix%sd_d(2, k, 4)*wrepp + ij_matrix%sd(k, 4)*wrepp_d(2) 2275 v_d(2, ij, 37) = v_d(2, ij, 37) + ij_matrix%sd_d(2, k, 5)*wrepp + ij_matrix%sd(k, 5)*wrepp_d(2) 2276 2277 v_d(3, ij, 11) = v_d(3, ij, 11) + ij_matrix%sd_d(3, k, 1)*wrepp + ij_matrix%sd(k, 1)*wrepp_d(3) 2278 v_d(3, ij, 16) = v_d(3, ij, 16) + ij_matrix%sd_d(3, k, 2)*wrepp + ij_matrix%sd(k, 2)*wrepp_d(3) 2279 v_d(3, ij, 22) = v_d(3, ij, 22) + ij_matrix%sd_d(3, k, 3)*wrepp + ij_matrix%sd(k, 3)*wrepp_d(3) 2280 v_d(3, ij, 29) = v_d(3, ij, 29) + ij_matrix%sd_d(3, k, 4)*wrepp + ij_matrix%sd(k, 4)*wrepp_d(3) 2281 v_d(3, ij, 37) = v_d(3, ij, 37) + ij_matrix%sd_d(3, k, 5)*wrepp + ij_matrix%sd(k, 5)*wrepp_d(3) 2282 ELSE IF (mm == 5) THEN 2283 k = k1 - 4 2284 l = l1 - 1 2285 v_d(1, ij, 12) = v_d(1, ij, 12) + ij_matrix%pd_d(1, 1, k, l)*wrepp + ij_matrix%pd(1, k, l)*wrepp_d(1) 2286 v_d(1, ij, 13) = v_d(1, ij, 13) + ij_matrix%pd_d(1, 2, k, l)*wrepp + ij_matrix%pd(2, k, l)*wrepp_d(1) 2287 v_d(1, ij, 14) = v_d(1, ij, 14) + ij_matrix%pd_d(1, 3, k, l)*wrepp + ij_matrix%pd(3, k, l)*wrepp_d(1) 2288 v_d(1, ij, 17) = v_d(1, ij, 17) + ij_matrix%pd_d(1, 4, k, l)*wrepp + ij_matrix%pd(4, k, l)*wrepp_d(1) 2289 v_d(1, ij, 18) = v_d(1, ij, 18) + ij_matrix%pd_d(1, 5, k, l)*wrepp + ij_matrix%pd(5, k, l)*wrepp_d(1) 2290 v_d(1, ij, 19) = v_d(1, ij, 19) + ij_matrix%pd_d(1, 6, k, l)*wrepp + ij_matrix%pd(6, k, l)*wrepp_d(1) 2291 v_d(1, ij, 23) = v_d(1, ij, 23) + ij_matrix%pd_d(1, 7, k, l)*wrepp + ij_matrix%pd(7, k, l)*wrepp_d(1) 2292 v_d(1, ij, 24) = v_d(1, ij, 24) + ij_matrix%pd_d(1, 8, k, l)*wrepp + ij_matrix%pd(8, k, l)*wrepp_d(1) 2293 v_d(1, ij, 25) = v_d(1, ij, 25) + ij_matrix%pd_d(1, 9, k, l)*wrepp + ij_matrix%pd(9, k, l)*wrepp_d(1) 2294 v_d(1, ij, 30) = v_d(1, ij, 30) + ij_matrix%pd_d(1, 10, k, l)*wrepp + ij_matrix%pd(10, k, l)*wrepp_d(1) 2295 v_d(1, ij, 31) = v_d(1, ij, 31) + ij_matrix%pd_d(1, 11, k, l)*wrepp + ij_matrix%pd(11, k, l)*wrepp_d(1) 2296 v_d(1, ij, 32) = v_d(1, ij, 32) + ij_matrix%pd_d(1, 12, k, l)*wrepp + ij_matrix%pd(12, k, l)*wrepp_d(1) 2297 v_d(1, ij, 38) = v_d(1, ij, 38) + ij_matrix%pd_d(1, 13, k, l)*wrepp + ij_matrix%pd(13, k, l)*wrepp_d(1) 2298 v_d(1, ij, 39) = v_d(1, ij, 39) + ij_matrix%pd_d(1, 14, k, l)*wrepp + ij_matrix%pd(14, k, l)*wrepp_d(1) 2299 v_d(1, ij, 40) = v_d(1, ij, 40) + ij_matrix%pd_d(1, 15, k, l)*wrepp + ij_matrix%pd(15, k, l)*wrepp_d(1) 2300 2301 v_d(2, ij, 12) = v_d(2, ij, 12) + ij_matrix%pd_d(2, 1, k, l)*wrepp + ij_matrix%pd(1, k, l)*wrepp_d(2) 2302 v_d(2, ij, 13) = v_d(2, ij, 13) + ij_matrix%pd_d(2, 2, k, l)*wrepp + ij_matrix%pd(2, k, l)*wrepp_d(2) 2303 v_d(2, ij, 14) = v_d(2, ij, 14) + ij_matrix%pd_d(2, 3, k, l)*wrepp + ij_matrix%pd(3, k, l)*wrepp_d(2) 2304 v_d(2, ij, 17) = v_d(2, ij, 17) + ij_matrix%pd_d(2, 4, k, l)*wrepp + ij_matrix%pd(4, k, l)*wrepp_d(2) 2305 v_d(2, ij, 18) = v_d(2, ij, 18) + ij_matrix%pd_d(2, 5, k, l)*wrepp + ij_matrix%pd(5, k, l)*wrepp_d(2) 2306 v_d(2, ij, 19) = v_d(2, ij, 19) + ij_matrix%pd_d(2, 6, k, l)*wrepp + ij_matrix%pd(6, k, l)*wrepp_d(2) 2307 v_d(2, ij, 23) = v_d(2, ij, 23) + ij_matrix%pd_d(2, 7, k, l)*wrepp + ij_matrix%pd(7, k, l)*wrepp_d(2) 2308 v_d(2, ij, 24) = v_d(2, ij, 24) + ij_matrix%pd_d(2, 8, k, l)*wrepp + ij_matrix%pd(8, k, l)*wrepp_d(2) 2309 v_d(2, ij, 25) = v_d(2, ij, 25) + ij_matrix%pd_d(2, 9, k, l)*wrepp + ij_matrix%pd(9, k, l)*wrepp_d(2) 2310 v_d(2, ij, 30) = v_d(2, ij, 30) + ij_matrix%pd_d(2, 10, k, l)*wrepp + ij_matrix%pd(10, k, l)*wrepp_d(2) 2311 v_d(2, ij, 31) = v_d(2, ij, 31) + ij_matrix%pd_d(2, 11, k, l)*wrepp + ij_matrix%pd(11, k, l)*wrepp_d(2) 2312 v_d(2, ij, 32) = v_d(2, ij, 32) + ij_matrix%pd_d(2, 12, k, l)*wrepp + ij_matrix%pd(12, k, l)*wrepp_d(2) 2313 v_d(2, ij, 38) = v_d(2, ij, 38) + ij_matrix%pd_d(2, 13, k, l)*wrepp + ij_matrix%pd(13, k, l)*wrepp_d(2) 2314 v_d(2, ij, 39) = v_d(2, ij, 39) + ij_matrix%pd_d(2, 14, k, l)*wrepp + ij_matrix%pd(14, k, l)*wrepp_d(2) 2315 v_d(2, ij, 40) = v_d(2, ij, 40) + ij_matrix%pd_d(2, 15, k, l)*wrepp + ij_matrix%pd(15, k, l)*wrepp_d(2) 2316 2317 v_d(3, ij, 12) = v_d(3, ij, 12) + ij_matrix%pd_d(3, 1, k, l)*wrepp + ij_matrix%pd(1, k, l)*wrepp_d(3) 2318 v_d(3, ij, 13) = v_d(3, ij, 13) + ij_matrix%pd_d(3, 2, k, l)*wrepp + ij_matrix%pd(2, k, l)*wrepp_d(3) 2319 v_d(3, ij, 14) = v_d(3, ij, 14) + ij_matrix%pd_d(3, 3, k, l)*wrepp + ij_matrix%pd(3, k, l)*wrepp_d(3) 2320 v_d(3, ij, 17) = v_d(3, ij, 17) + ij_matrix%pd_d(3, 4, k, l)*wrepp + ij_matrix%pd(4, k, l)*wrepp_d(3) 2321 v_d(3, ij, 18) = v_d(3, ij, 18) + ij_matrix%pd_d(3, 5, k, l)*wrepp + ij_matrix%pd(5, k, l)*wrepp_d(3) 2322 v_d(3, ij, 19) = v_d(3, ij, 19) + ij_matrix%pd_d(3, 6, k, l)*wrepp + ij_matrix%pd(6, k, l)*wrepp_d(3) 2323 v_d(3, ij, 23) = v_d(3, ij, 23) + ij_matrix%pd_d(3, 7, k, l)*wrepp + ij_matrix%pd(7, k, l)*wrepp_d(3) 2324 v_d(3, ij, 24) = v_d(3, ij, 24) + ij_matrix%pd_d(3, 8, k, l)*wrepp + ij_matrix%pd(8, k, l)*wrepp_d(3) 2325 v_d(3, ij, 25) = v_d(3, ij, 25) + ij_matrix%pd_d(3, 9, k, l)*wrepp + ij_matrix%pd(9, k, l)*wrepp_d(3) 2326 v_d(3, ij, 30) = v_d(3, ij, 30) + ij_matrix%pd_d(3, 10, k, l)*wrepp + ij_matrix%pd(10, k, l)*wrepp_d(3) 2327 v_d(3, ij, 31) = v_d(3, ij, 31) + ij_matrix%pd_d(3, 11, k, l)*wrepp + ij_matrix%pd(11, k, l)*wrepp_d(3) 2328 v_d(3, ij, 32) = v_d(3, ij, 32) + ij_matrix%pd_d(3, 12, k, l)*wrepp + ij_matrix%pd(12, k, l)*wrepp_d(3) 2329 v_d(3, ij, 38) = v_d(3, ij, 38) + ij_matrix%pd_d(3, 13, k, l)*wrepp + ij_matrix%pd(13, k, l)*wrepp_d(3) 2330 v_d(3, ij, 39) = v_d(3, ij, 39) + ij_matrix%pd_d(3, 14, k, l)*wrepp + ij_matrix%pd(14, k, l)*wrepp_d(3) 2331 v_d(3, ij, 40) = v_d(3, ij, 40) + ij_matrix%pd_d(3, 15, k, l)*wrepp + ij_matrix%pd(15, k, l)*wrepp_d(3) 2332 ELSE IF (mm == 6) THEN 2333 k = k1 - 4 2334 l = l1 - 4 2335 v_d(1, ij, 15) = v_d(1, ij, 15) + ij_matrix%dd_d(1, 1, k, l)*wrepp + ij_matrix%dd(1, k, l)*wrepp_d(1) 2336 v_d(1, ij, 21) = v_d(1, ij, 21) + ij_matrix%dd_d(1, 2, k, l)*wrepp + ij_matrix%dd(2, k, l)*wrepp_d(1) 2337 v_d(1, ij, 28) = v_d(1, ij, 28) + ij_matrix%dd_d(1, 3, k, l)*wrepp + ij_matrix%dd(3, k, l)*wrepp_d(1) 2338 v_d(1, ij, 36) = v_d(1, ij, 36) + ij_matrix%dd_d(1, 4, k, l)*wrepp + ij_matrix%dd(4, k, l)*wrepp_d(1) 2339 v_d(1, ij, 45) = v_d(1, ij, 45) + ij_matrix%dd_d(1, 5, k, l)*wrepp + ij_matrix%dd(5, k, l)*wrepp_d(1) 2340 v_d(1, ij, 20) = v_d(1, ij, 20) + ij_matrix%dd_d(1, 6, k, l)*wrepp + ij_matrix%dd(6, k, l)*wrepp_d(1) 2341 v_d(1, ij, 26) = v_d(1, ij, 26) + ij_matrix%dd_d(1, 7, k, l)*wrepp + ij_matrix%dd(7, k, l)*wrepp_d(1) 2342 v_d(1, ij, 27) = v_d(1, ij, 27) + ij_matrix%dd_d(1, 8, k, l)*wrepp + ij_matrix%dd(8, k, l)*wrepp_d(1) 2343 v_d(1, ij, 33) = v_d(1, ij, 33) + ij_matrix%dd_d(1, 9, k, l)*wrepp + ij_matrix%dd(9, k, l)*wrepp_d(1) 2344 v_d(1, ij, 34) = v_d(1, ij, 34) + ij_matrix%dd_d(1, 10, k, l)*wrepp + ij_matrix%dd(10, k, l)*wrepp_d(1) 2345 v_d(1, ij, 35) = v_d(1, ij, 35) + ij_matrix%dd_d(1, 11, k, l)*wrepp + ij_matrix%dd(11, k, l)*wrepp_d(1) 2346 v_d(1, ij, 41) = v_d(1, ij, 41) + ij_matrix%dd_d(1, 12, k, l)*wrepp + ij_matrix%dd(12, k, l)*wrepp_d(1) 2347 v_d(1, ij, 42) = v_d(1, ij, 42) + ij_matrix%dd_d(1, 13, k, l)*wrepp + ij_matrix%dd(13, k, l)*wrepp_d(1) 2348 v_d(1, ij, 43) = v_d(1, ij, 43) + ij_matrix%dd_d(1, 14, k, l)*wrepp + ij_matrix%dd(14, k, l)*wrepp_d(1) 2349 v_d(1, ij, 44) = v_d(1, ij, 44) + ij_matrix%dd_d(1, 15, k, l)*wrepp + ij_matrix%dd(15, k, l)*wrepp_d(1) 2350 2351 v_d(2, ij, 15) = v_d(2, ij, 15) + ij_matrix%dd_d(2, 1, k, l)*wrepp + ij_matrix%dd(1, k, l)*wrepp_d(2) 2352 v_d(2, ij, 21) = v_d(2, ij, 21) + ij_matrix%dd_d(2, 2, k, l)*wrepp + ij_matrix%dd(2, k, l)*wrepp_d(2) 2353 v_d(2, ij, 28) = v_d(2, ij, 28) + ij_matrix%dd_d(2, 3, k, l)*wrepp + ij_matrix%dd(3, k, l)*wrepp_d(2) 2354 v_d(2, ij, 36) = v_d(2, ij, 36) + ij_matrix%dd_d(2, 4, k, l)*wrepp + ij_matrix%dd(4, k, l)*wrepp_d(2) 2355 v_d(2, ij, 45) = v_d(2, ij, 45) + ij_matrix%dd_d(2, 5, k, l)*wrepp + ij_matrix%dd(5, k, l)*wrepp_d(2) 2356 v_d(2, ij, 20) = v_d(2, ij, 20) + ij_matrix%dd_d(2, 6, k, l)*wrepp + ij_matrix%dd(6, k, l)*wrepp_d(2) 2357 v_d(2, ij, 26) = v_d(2, ij, 26) + ij_matrix%dd_d(2, 7, k, l)*wrepp + ij_matrix%dd(7, k, l)*wrepp_d(2) 2358 v_d(2, ij, 27) = v_d(2, ij, 27) + ij_matrix%dd_d(2, 8, k, l)*wrepp + ij_matrix%dd(8, k, l)*wrepp_d(2) 2359 v_d(2, ij, 33) = v_d(2, ij, 33) + ij_matrix%dd_d(2, 9, k, l)*wrepp + ij_matrix%dd(9, k, l)*wrepp_d(2) 2360 v_d(2, ij, 34) = v_d(2, ij, 34) + ij_matrix%dd_d(2, 10, k, l)*wrepp + ij_matrix%dd(10, k, l)*wrepp_d(2) 2361 v_d(2, ij, 35) = v_d(2, ij, 35) + ij_matrix%dd_d(2, 11, k, l)*wrepp + ij_matrix%dd(11, k, l)*wrepp_d(2) 2362 v_d(2, ij, 41) = v_d(2, ij, 41) + ij_matrix%dd_d(2, 12, k, l)*wrepp + ij_matrix%dd(12, k, l)*wrepp_d(2) 2363 v_d(2, ij, 42) = v_d(2, ij, 42) + ij_matrix%dd_d(2, 13, k, l)*wrepp + ij_matrix%dd(13, k, l)*wrepp_d(2) 2364 v_d(2, ij, 43) = v_d(2, ij, 43) + ij_matrix%dd_d(2, 14, k, l)*wrepp + ij_matrix%dd(14, k, l)*wrepp_d(2) 2365 v_d(2, ij, 44) = v_d(2, ij, 44) + ij_matrix%dd_d(2, 15, k, l)*wrepp + ij_matrix%dd(15, k, l)*wrepp_d(2) 2366 2367 v_d(3, ij, 15) = v_d(3, ij, 15) + ij_matrix%dd_d(3, 1, k, l)*wrepp + ij_matrix%dd(1, k, l)*wrepp_d(3) 2368 v_d(3, ij, 21) = v_d(3, ij, 21) + ij_matrix%dd_d(3, 2, k, l)*wrepp + ij_matrix%dd(2, k, l)*wrepp_d(3) 2369 v_d(3, ij, 28) = v_d(3, ij, 28) + ij_matrix%dd_d(3, 3, k, l)*wrepp + ij_matrix%dd(3, k, l)*wrepp_d(3) 2370 v_d(3, ij, 36) = v_d(3, ij, 36) + ij_matrix%dd_d(3, 4, k, l)*wrepp + ij_matrix%dd(4, k, l)*wrepp_d(3) 2371 v_d(3, ij, 45) = v_d(3, ij, 45) + ij_matrix%dd_d(3, 5, k, l)*wrepp + ij_matrix%dd(5, k, l)*wrepp_d(3) 2372 v_d(3, ij, 20) = v_d(3, ij, 20) + ij_matrix%dd_d(3, 6, k, l)*wrepp + ij_matrix%dd(6, k, l)*wrepp_d(3) 2373 v_d(3, ij, 26) = v_d(3, ij, 26) + ij_matrix%dd_d(3, 7, k, l)*wrepp + ij_matrix%dd(7, k, l)*wrepp_d(3) 2374 v_d(3, ij, 27) = v_d(3, ij, 27) + ij_matrix%dd_d(3, 8, k, l)*wrepp + ij_matrix%dd(8, k, l)*wrepp_d(3) 2375 v_d(3, ij, 33) = v_d(3, ij, 33) + ij_matrix%dd_d(3, 9, k, l)*wrepp + ij_matrix%dd(9, k, l)*wrepp_d(3) 2376 v_d(3, ij, 34) = v_d(3, ij, 34) + ij_matrix%dd_d(3, 10, k, l)*wrepp + ij_matrix%dd(10, k, l)*wrepp_d(3) 2377 v_d(3, ij, 35) = v_d(3, ij, 35) + ij_matrix%dd_d(3, 11, k, l)*wrepp + ij_matrix%dd(11, k, l)*wrepp_d(3) 2378 v_d(3, ij, 41) = v_d(3, ij, 41) + ij_matrix%dd_d(3, 12, k, l)*wrepp + ij_matrix%dd(12, k, l)*wrepp_d(3) 2379 v_d(3, ij, 42) = v_d(3, ij, 42) + ij_matrix%dd_d(3, 13, k, l)*wrepp + ij_matrix%dd(13, k, l)*wrepp_d(3) 2380 v_d(3, ij, 43) = v_d(3, ij, 43) + ij_matrix%dd_d(3, 14, k, l)*wrepp + ij_matrix%dd(14, k, l)*wrepp_d(3) 2381 v_d(3, ij, 44) = v_d(3, ij, 44) + ij_matrix%dd_d(3, 15, k, l)*wrepp + ij_matrix%dd(15, k, l)*wrepp_d(3) 2382 END IF 2383 END IF 2384 END DO 2385 END DO 2386 DO kl = 1, limkl 2387 logv_d(ij, kl) = (ANY(ABS(v_d(1:3, ij, kl)) > 0.0_dp)) 2388 END DO 2389 END DO 2390 END DO 2391 IF (debug_this_module) THEN 2392 CALL rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, v_d) 2393 END IF 2394 END IF 2395 END SUBROUTINE rot_2el_2c_first 2396 2397! ************************************************************************************************** 2398!> \brief Store the two-electron two-center integrals in diagonal form 2399!> 2400!> \param limij ... 2401!> \param limkl ... 2402!> \param ww ... 2403!> \param w ... 2404!> \param ww_dx ... 2405!> \param ww_dy ... 2406!> \param ww_dz ... 2407!> \param dw ... 2408!> \date 04.2008 [tlaino] 2409!> \author Teodoro Laino [tlaino] - University of Zurich 2410! ************************************************************************************************** 2411 SUBROUTINE store_2el_2c_diag(limij, limkl, ww, w, ww_dx, ww_dy, ww_dz, dw) 2412 2413 INTEGER, INTENT(IN) :: limij, limkl 2414 REAL(KIND=dp), DIMENSION(limkl, limij), & 2415 INTENT(IN), OPTIONAL :: ww 2416 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), & 2417 OPTIONAL :: w 2418 REAL(KIND=dp), DIMENSION(limkl, limij), & 2419 INTENT(IN), OPTIONAL :: ww_dx, ww_dy, ww_dz 2420 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), & 2421 OPTIONAL :: dw 2422 2423 CHARACTER(len=*), PARAMETER :: routineN = 'store_2el_2c_diag', & 2424 routineP = moduleN//':'//routineN 2425 2426 INTEGER :: ij, kl, l 2427 2428 l = 0 2429 IF (PRESENT(ww) .AND. PRESENT(w)) THEN 2430 DO ij = 1, limij 2431 DO kl = 1, limkl 2432 l = l + 1 2433 w(l) = ww(kl, ij) 2434 END DO 2435 END DO 2436 ELSE IF (PRESENT(ww_dx) .AND. PRESENT(ww_dy) .AND. PRESENT(ww_dz) .AND. PRESENT(dw)) THEN 2437 DO ij = 1, limij 2438 DO kl = 1, limkl 2439 l = l + 1 2440 dw(1, l) = ww_dx(kl, ij) 2441 dw(2, l) = ww_dy(kl, ij) 2442 dw(3, l) = ww_dz(kl, ij) 2443 END DO 2444 END DO 2445 ELSE 2446 CPABORT("") 2447 END IF 2448 2449 END SUBROUTINE store_2el_2c_diag 2450 2451END MODULE semi_empirical_int_utils 2452