|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00010 MODULE semi_empirical_int_utils 00011 00012 USE f77_blas 00013 USE input_constants, ONLY: do_method_pchg,& 00014 do_se_IS_kdso_d 00015 USE kinds, ONLY: dp 00016 USE semi_empirical_int3_utils, ONLY: charg_int_3,& 00017 dcharg_int_3,& 00018 ijkl_low_3 00019 USE semi_empirical_int_arrays, ONLY: & 00020 CLMp, CLMxx, CLMxy, CLMyy, CLMz, CLMzp, CLMzz, clm_d, clm_sp, & 00021 ijkl_ind, indexa, indexb, int2c_type 00022 USE semi_empirical_types, ONLY: rotmat_type,& 00023 se_int_control_type,& 00024 se_int_screen_type,& 00025 se_taper_type,& 00026 semi_empirical_type 00027 #include "cp_common_uses.h" 00028 00029 IMPLICIT NONE 00030 #include "semi_empirical_int_debug.h" 00031 00032 PRIVATE 00033 LOGICAL, PARAMETER, PRIVATE :: debug_this_module=.FALSE. 00034 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'semi_empirical_int_utils' 00035 00036 PUBLIC :: ijkl_sp, ijkl_d, rotmat, rot_2el_2c_first, store_2el_2c_diag,& 00037 d_ijkl_sp, d_ijkl_d 00038 00039 CONTAINS 00040 00041 ! ***************************************************************************** 00051 FUNCTION ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control,& 00052 se_int_screen, itype, error) RESULT(res) 00053 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00054 INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic 00055 REAL(KIND=dp), INTENT(IN) :: r 00056 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00057 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 00058 INTEGER, INTENT(IN) :: itype 00059 TYPE(cp_error_type), INTENT(inout) :: error 00060 REAL(KIND=dp) :: res 00061 00062 CHARACTER(len=*), PARAMETER :: routineN = 'ijkl_sp', 00063 routineP = moduleN//':'//routineN 00064 00065 res = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,& 00066 se_int_control%integral_screening, se_int_control%shortrange,& 00067 se_int_control%pc_coulomb_int, se_int_control%max_multipole,& 00068 itype, charg_int_nri, error) 00069 00070 ! If only the shortrange component is requested we can skip the rest 00071 IF ((.NOT.se_int_control%pc_coulomb_int).AND.(itype/=do_method_pchg)) THEN 00072 ! Handle the 1/r^3 term, this term is ALWAYS false for KDSO-D integrals 00073 IF (se_int_control%shortrange.AND.se_int_control%do_ewald_r3) THEN 00074 res = res - ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r,& 00075 itype, charg_int_3, error) 00076 END IF 00077 END IF 00078 END FUNCTION ijkl_sp 00079 00080 ! ***************************************************************************** 00090 FUNCTION d_ijkl_sp(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control,& 00091 se_int_screen, itype, error) RESULT(res) 00092 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00093 INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic 00094 REAL(KIND=dp), INTENT(IN) :: r 00095 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00096 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 00097 INTEGER, INTENT(IN) :: itype 00098 TYPE(cp_error_type), INTENT(inout) :: error 00099 REAL(KIND=dp) :: res 00100 00101 CHARACTER(len=*), PARAMETER :: routineN = 'd_ijkl_sp', 00102 routineP = moduleN//':'//routineN 00103 00104 REAL(KIND=dp) :: dfs, srd 00105 00106 IF (se_int_control%integral_screening==do_se_IS_kdso_d) THEN 00107 res = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,& 00108 se_int_control%integral_screening, .FALSE.,& 00109 se_int_control%pc_coulomb_int, se_int_control%max_multipole,& 00110 itype, dcharg_int_nri, error) 00111 00112 IF (.NOT.se_int_control%pc_coulomb_int) THEN 00113 dfs = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,& 00114 se_int_control%integral_screening, .FALSE., .FALSE., & 00115 se_int_control%max_multipole, itype, dcharg_int_nri_fs, error) 00116 res = res + dfs*se_int_screen%dft 00117 00118 ! In case we need the shortrange part we have to evaluate an additional derivative 00119 ! to handle the derivative of the Tapering term 00120 IF (se_int_control%shortrange) THEN 00121 srd = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,& 00122 se_int_control%integral_screening, .FALSE., .TRUE., & 00123 se_int_control%max_multipole, itype, dcharg_int_nri, error) 00124 res = res - srd 00125 END IF 00126 END IF 00127 ELSE 00128 res = ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,& 00129 se_int_control%integral_screening, se_int_control%shortrange,& 00130 se_int_control%pc_coulomb_int, se_int_control%max_multipole,& 00131 itype, dcharg_int_nri, error) 00132 END IF 00133 00134 ! If only the shortrange component is requested we can skip the rest 00135 IF ((.NOT.se_int_control%pc_coulomb_int).AND.(itype/=do_method_pchg)) THEN 00136 ! Handle the 1/r^3 term, this term is ALWAYS false for KDSO-D integrals 00137 IF (se_int_control%shortrange.AND.se_int_control%do_ewald_r3) THEN 00138 res = res - ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r,& 00139 itype, dcharg_int_3, error) 00140 END IF 00141 END IF 00142 00143 END FUNCTION d_ijkl_sp 00144 00145 ! ***************************************************************************** 00156 FUNCTION ijkl_sp_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,& 00157 iscreen, shortrange, pc_coulomb_int, max_multipole, itype, eval, error) RESULT(res) 00158 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00159 INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic 00160 REAL(KIND=dp), INTENT(IN) :: r 00161 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 00162 INTEGER, INTENT(IN) :: iscreen 00163 LOGICAL, INTENT(IN) :: shortrange, pc_coulomb_int 00164 INTEGER, INTENT(IN) :: max_multipole, itype 00165 REAL(KIND=dp) :: eval 00166 TYPE(cp_error_type), INTENT(inout) :: error 00167 REAL(KIND=dp) :: res 00168 00169 CHARACTER(len=*), PARAMETER :: routineN = 'ijkl_sp_low', 00170 routineP = moduleN//':'//routineN 00171 00172 INTEGER :: ccc, l1, l1max, l1min, l2, 00173 l2max, l2min, lij, lkl, lmin, 00174 m 00175 REAL(KIND=dp) :: add, chrg, dij, dkl, fact_ij, 00176 fact_kl, fact_screen, pij, 00177 pkl, s1, sum 00178 00179 l1min = ABS (li-lj) 00180 l1max = li + lj 00181 lij = indexb(li+1, lj+1) 00182 l2min = ABS (lk-ll) 00183 l2max = lk + ll 00184 lkl = indexb(lk+1, ll+1) 00185 00186 l1max = MIN (l1max, 2) 00187 l1min = MIN (l1min, 2) 00188 l2max = MIN (l2max, 2) 00189 l2min = MIN (l2min, 2) 00190 sum = 0.0_dp 00191 dij = 0.0_dp 00192 dkl = 0.0_dp 00193 fact_ij = 1.0_dp 00194 fact_kl = 1.0_dp 00195 fact_screen = 1.0_dp 00196 IF (lij==3) fact_ij = SQRT(2.0_dp) 00197 IF (lkl==3) fact_kl = SQRT(2.0_dp) 00198 IF (.NOT.pc_coulomb_int) THEN 00199 IF (iscreen==do_se_IS_kdso_d) fact_screen = se_int_screen%ft 00200 ! Standard value of the integral 00201 DO l1 = l1min, l1max 00202 IF (l1 == 0) THEN 00203 IF (lij == 1) THEN 00204 pij = sepi%ko(1) 00205 IF (ic == -1 .OR. ic == 1) THEN 00206 pij = sepi%ko(9) 00207 END IF 00208 ELSE IF (lij == 3) THEN 00209 pij = sepi%ko(7) 00210 END IF 00211 ELSE 00212 dij = sepi%cs(lij)*fact_ij 00213 pij = sepi%ko(lij) 00214 END IF 00215 ! 00216 DO l2 = l2min, l2max 00217 IF (l2 == 0) THEN 00218 IF (lkl == 1) THEN 00219 pkl = sepj%ko(1) 00220 IF (ic == -1 .OR. ic == 2) THEN 00221 pkl = sepj%ko(9) 00222 END IF 00223 ELSE IF (lkl == 3) THEN 00224 pkl = sepj%ko(7) 00225 END IF 00226 ELSE 00227 dkl = sepj%cs(lkl)*fact_kl 00228 pkl = sepj%ko(lkl) 00229 END IF 00230 IF (itype==do_method_pchg) THEN 00231 add = 0.0_dp 00232 ELSE 00233 add = (pij+pkl) ** 2 00234 END IF 00235 lmin = MAX(l1, l2) 00236 s1 = 0.0_dp 00237 DO m = -lmin, lmin 00238 ccc = clm_sp(ij, l1, m) * clm_sp(kl, l2, m) 00239 IF (ABS(ccc) > EPSILON(0.0_dp)) THEN 00240 chrg = eval(r, l1, l2, clm_sp(ij, l1, m), clm_sp(kl, l2, m), dij, dkl, add, fact_screen, error) 00241 s1 = s1 + chrg 00242 END IF 00243 END DO 00244 sum = sum + s1 00245 END DO 00246 END DO 00247 res = sum 00248 END IF 00249 ! Shortrange: Possibly computes pure Coulomb and subtract from the original integral value 00250 IF (shortrange.OR.pc_coulomb_int) THEN 00251 sum = 0.0_dp 00252 dij = 0.0_dp 00253 dkl = 0.0_dp 00254 add = 0.0_dp 00255 fact_screen = 0.0_dp 00256 DO l1 = l1min, l1max 00257 IF (l1 > max_multipole) CYCLE 00258 IF (l1 /= 0) THEN 00259 dij = sepi%cs(lij)*fact_ij 00260 END IF 00261 ! 00262 DO l2 = l2min, l2max 00263 IF (l2 > max_multipole) CYCLE 00264 IF (l2 /= 0) THEN 00265 dkl = sepj%cs(lkl)*fact_kl 00266 END IF 00267 lmin = MAX(l1, l2) 00268 s1 = 0.0_dp 00269 DO m = -lmin, lmin 00270 ccc = clm_sp(ij, l1, m) * clm_sp(kl, l2, m) 00271 IF (ABS(ccc) > EPSILON(0.0_dp)) THEN 00272 chrg = eval(r, l1, l2, clm_sp(ij, l1, m), clm_sp(kl, l2, m), dij, dkl, add, fact_screen, error) 00273 s1 = s1 + chrg 00274 END IF 00275 END DO 00276 sum = sum + s1 00277 END DO 00278 END DO 00279 IF (pc_coulomb_int) res = sum 00280 IF (shortrange) res = res - sum 00281 END IF 00282 END FUNCTION ijkl_sp_low 00283 00284 ! ***************************************************************************** 00297 FUNCTION charg_int_nri(r, l1_i, l2_i, m1_i, m2_i, da_i, db_i, add0, fact_screen, error) RESULT(charg) 00298 REAL(KIND=dp), INTENT(in) :: r 00299 INTEGER, INTENT(in) :: l1_i, l2_i, m1_i, m2_i 00300 REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen 00301 TYPE(cp_error_type), INTENT(inout) :: error 00302 REAL(KIND=dp) :: charg 00303 00304 CHARACTER(len=*), PARAMETER :: routineN = 'charg_int_nri', 00305 routineP = moduleN//':'//routineN 00306 00307 INTEGER :: l1, l2, m1, m2 00308 LOGICAL :: failure 00309 REAL(KIND=dp) :: add, da, db, dxdx, dxqxz, 00310 dzdz, dzqxx, dzqzz, fact, 00311 qqxx, qqzz, qxxqxx, qxxqyy, 00312 qxzqxz, xyxy, zzzz 00313 00314 failure = .FALSE. 00315 ! Computing only Integral Values 00316 IF ( l1_i < l2_i) THEN 00317 l1 = l1_i 00318 l2 = l2_i 00319 m1 = m1_i 00320 m2 = m2_i 00321 da = da_i 00322 db = db_i 00323 fact = 1.0_dp 00324 ELSE IF ( l1_i > l2_i ) THEN 00325 l1 = l2_i 00326 l2 = l1_i 00327 m1 = m2_i 00328 m2 = m1_i 00329 da = db_i 00330 db = da_i 00331 fact = (-1.0_dp)**(l1+l2) 00332 ELSE IF ( l1_i == l2_i ) THEN 00333 l1 = l1_i 00334 l2 = l2_i 00335 IF ( m1_i <= m2_i ) THEN 00336 m1 = m1_i 00337 m2 = m2_i 00338 da = da_i 00339 db = db_i 00340 ELSE 00341 m1 = m2_i 00342 m2 = m1_i 00343 da = db_i 00344 db = da_i 00345 END IF 00346 fact = 1.0_dp 00347 END IF 00348 add = add0*fact_screen 00349 charg = 0.0_dp 00350 ! Q - Q. 00351 IF (l1 == 0 .AND. l2 == 0) THEN 00352 charg = fact/SQRT(r**2+add) 00353 RETURN 00354 END IF 00355 ! Q - Z. 00356 IF (l1 == 0 .AND. l2 == 1 .AND. m2 == CLMz) THEN 00357 charg = 1.0_dp/SQRT((r+db)**2+add) - 1.0_dp/SQRT((r-db)**2+add) 00358 charg = charg*0.5_dp*fact 00359 RETURN 00360 END IF 00361 ! Z - Z. 00362 IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMz .AND. m2 == CLMz ) THEN 00363 dzdz = & 00364 +1.0_dp/SQRT((r+da-db)**2+add) + 1.0_dp/SQRT((r-da+db)**2+add) & 00365 -1.0_dp/SQRT((r-da-db)**2+add) - 1.0_dp/SQRT((r+da+db)**2+add) 00366 charg = dzdz*0.25_dp*fact 00367 RETURN 00368 END IF 00369 ! X - X 00370 IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMp .AND. m2 == CLMp ) THEN 00371 dxdx = 2.0_dp/SQRT(r**2+(da-db)**2+add) - 2.0_dp/SQRT(r**2+(da+db)**2+add) 00372 charg = dxdx*0.25_dp*fact 00373 RETURN 00374 END IF 00375 ! Q - ZZ 00376 IF (l1 == 0 .AND. l2 == 2 .AND. m2 == CLMzz) THEN 00377 qqzz = 1.0_dp/SQRT((r-db)**2+add) - 2.0_dp/SQRT(r**2+add) + 1.0_dp/SQRT((r+db)**2+add) 00378 charg = qqzz*0.25_dp*fact 00379 RETURN 00380 END IF 00381 ! Q - XX 00382 IF (l1 == 0 .AND. l2 == 2 .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN 00383 qqxx = - 1.0_dp/SQRT(r**2+add) + 1.0_dp/SQRT(r**2+add+db**2) 00384 charg = qqxx*0.5_dp*fact 00385 RETURN 00386 END IF 00387 ! Z - ZZ 00388 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. m2 == CLMzz ) THEN 00389 dzqzz = & 00390 +1.0_dp/SQRT((r-da-db)**2+add) - 2.0_dp/SQRT((r-da)**2+add) & 00391 +1.0_dp/SQRT((r-da+db)**2+add) - 1.0_dp/SQRT((r+da-db)**2+add) & 00392 +2.0_dp/SQRT((r+da)**2+add) - 1.0_dp/SQRT((r+da+db)**2+add) 00393 charg = dzqzz*0.125_dp*fact 00394 RETURN 00395 END IF 00396 ! Z - XX 00397 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. ( m2 == CLMyy .OR. m2 == CLMxx )) THEN 00398 dzqxx = & 00399 +1.0_dp/SQRT((r+da)**2+add) - 1.0_dp/SQRT((r+da)**2+add+db**2) & 00400 -1.0_dp/SQRT((r-da)**2+add) + 1.0_dp/SQRT((r-da)**2+add+db**2) 00401 charg = dzqxx*0.25_dp*fact 00402 RETURN 00403 END IF 00404 ! ZZ - ZZ 00405 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. m2 == CLMzz ) THEN 00406 zzzz = & 00407 +1.0_dp/SQRT((r-da-db)**2+add) + 1.0_dp/SQRT((r+da+db)**2+add) & 00408 +1.0_dp/SQRT((r-da+db)**2+add) + 1.0_dp/SQRT((r+da-db)**2+add) 00409 xyxy = & 00410 +1.0_dp/SQRT((r-da)**2+add) + 1.0_dp/SQRT((r+da)**2+add) & 00411 +1.0_dp/SQRT((r-db)**2+add) + 1.0_dp/SQRT((r+db)**2+add) & 00412 -2.0_dp/SQRT(r**2+add) 00413 charg = (zzzz*0.0625_dp - xyxy*0.125_dp)*fact 00414 RETURN 00415 END IF 00416 ! ZZ - XX 00417 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. ( m2 == CLMxx .OR. m2 == CLMyy )) THEN 00418 zzzz = & 00419 -1.0_dp/SQRT((r+da)**2+add) + 1.0_dp/SQRT((r+da)**2+db**2+add) & 00420 -1.0_dp/SQRT((r-da)**2+add) + 1.0_dp/SQRT((r-da)**2+db**2+add) 00421 xyxy = & 00422 +1.0_dp/SQRT(r**2+db**2+add) - 1.0_dp/SQRT(r**2+add) 00423 charg = (zzzz*0.125_dp - xyxy*0.25_dp)*fact 00424 RETURN 00425 END IF 00426 ! X - ZX 00427 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMp .AND. m2 == CLMzp ) THEN 00428 db = db/2.0_dp 00429 dxqxz = & 00430 -1.0_dp/SQRT((r-db)**2+(da-db)**2+add) + 1.0_dp/SQRT((r+db)**2+(da-db)**2+add) & 00431 +1.0_dp/SQRT((r-db)**2+(da+db)**2+add) - 1.0_dp/SQRT((r+db)**2+(da+db)**2+add) 00432 charg = dxqxz*0.25_dp*fact 00433 RETURN 00434 END IF 00435 ! ZX - ZX 00436 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzp .AND. m2 == CLMzp ) THEN 00437 da = da / 2.0_dp 00438 db = db / 2.0_dp 00439 qxzqxz = & 00440 +1.0_dp/SQRT((r+da-db)**2+(da-db)**2+add) - 1.0_dp/SQRT((r+da+db)**2+(da-db)**2+add) & 00441 -1.0_dp/SQRT((r-da-db)**2+(da-db)**2+add) + 1.0_dp/SQRT((r-da+db)**2+(da-db)**2+add) & 00442 -1.0_dp/SQRT((r+da-db)**2+(da+db)**2+add) + 1.0_dp/SQRT((r+da+db)**2+(da+db)**2+add) & 00443 +1.0_dp/SQRT((r-da-db)**2+(da+db)**2+add) - 1.0_dp/SQRT((r-da+db)**2+(da+db)**2+add) 00444 charg = qxzqxz*0.125_dp*fact 00445 RETURN 00446 END IF 00447 ! XX - XX 00448 IF (l1 == 2 .AND. l2 == 2 .AND. (((m1 == CLMyy) .AND. (m2 == CLMyy)).OR.((m1 == CLMxx) .AND. (m2 == CLMxx)))) THEN 00449 qxxqxx = & 00450 + 2.0_dp/SQRT(r**2+add) + 1.0_dp/SQRT(r**2+(da-db)**2+add) & 00451 + 1.0_dp/SQRT(r**2+(da+db)**2+add) - 2.0_dp/SQRT(r**2+da**2+add) & 00452 - 2.0_dp/SQRT(r**2+db**2+add) 00453 charg = qxxqxx*0.125_dp*fact 00454 RETURN 00455 END IF 00456 ! XX - YY 00457 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMyy .AND. m2 == CLMxx) THEN 00458 qxxqyy = & 00459 + 1.0_dp/SQRT(r**2+add) - 1.0_dp/SQRT(r**2+da**2+add) & 00460 - 1.0_dp/SQRT(r**2+db**2+add) + 1.0_dp/SQRT(r**2+da**2+db**2+add) 00461 charg = qxxqyy*0.25_dp*fact 00462 RETURN 00463 END IF 00464 ! XY - XY 00465 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMxy .AND. m2 == CLMxy) THEN 00466 qxxqxx = & 00467 + 2.0_dp/SQRT(r**2+add) + 1.0_dp/SQRT(r**2+(da-db)**2+add) & 00468 + 1.0_dp/SQRT(r**2+(da+db)**2+add) - 2.0_dp/SQRT(r**2+da**2+add) & 00469 - 2.0_dp/SQRT(r**2+db**2+add) 00470 qxxqyy = & 00471 + 1.0_dp/SQRT(r**2+add) - 1.0_dp/SQRT(r**2+da**2+add) & 00472 - 1.0_dp/SQRT(r**2+db**2+add) + 1.0_dp/SQRT(r**2+da**2+db**2+add) 00473 charg = 0.5_dp * (qxxqxx*0.125_dp - qxxqyy*0.25_dp) * fact 00474 RETURN 00475 END IF 00476 ! We should NEVER reach this point 00477 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 00478 END FUNCTION charg_int_nri 00479 00480 ! ***************************************************************************** 00495 FUNCTION dcharg_int_nri(r, l1_i, l2_i, m1_i, m2_i, da_i, db_i, add0, fact_screen, error) RESULT(charg) 00496 REAL(KIND=dp), INTENT(in) :: r 00497 INTEGER, INTENT(in) :: l1_i, l2_i, m1_i, m2_i 00498 REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen 00499 TYPE(cp_error_type), INTENT(inout) :: error 00500 REAL(KIND=dp) :: charg 00501 00502 CHARACTER(len=*), PARAMETER :: routineN = 'dcharg_int_nri', 00503 routineP = moduleN//':'//routineN 00504 00505 INTEGER :: l1, l2, m1, m2 00506 LOGICAL :: failure 00507 REAL(KIND=dp) :: add, da, db, dxdx, dxqxz, 00508 dzdz, dzqxx, dzqzz, fact, 00509 qqxx, qqzz, qxxqxx, qxxqyy, 00510 qxzqxz, xyxy, zzzz 00511 00512 failure = .FALSE. 00513 ! Computing only Integral Derivatives 00514 IF ( l1_i < l2_i) THEN 00515 l1 = l1_i 00516 l2 = l2_i 00517 m1 = m1_i 00518 m2 = m2_i 00519 da = da_i 00520 db = db_i 00521 fact = 1.0_dp 00522 ELSE IF ( l1_i > l2_i ) THEN 00523 l1 = l2_i 00524 l2 = l1_i 00525 m1 = m2_i 00526 m2 = m1_i 00527 da = db_i 00528 db = da_i 00529 fact = (-1.0_dp)**(l1+l2) 00530 ELSE IF ( l1_i == l2_i ) THEN 00531 l1 = l1_i 00532 l2 = l2_i 00533 IF ( m1_i <= m2_i ) THEN 00534 m1 = m1_i 00535 m2 = m2_i 00536 da = da_i 00537 db = db_i 00538 ELSE 00539 m1 = m2_i 00540 m2 = m1_i 00541 da = db_i 00542 db = da_i 00543 END IF 00544 fact = 1.0_dp 00545 END IF 00546 charg = 0.0_dp 00547 add = add0 * fact_screen 00548 ! Q - Q. 00549 IF (l1 == 0 .AND. l2 == 0) THEN 00550 charg = r/SQRT(r**2+add)**3 00551 charg = -charg*fact 00552 RETURN 00553 END IF 00554 ! Q - Z. 00555 IF (l1 == 0 .AND. l2 == 1 .AND. m2 == CLMz) THEN 00556 charg = (r+db)/SQRT((r+db)**2+add)**3 - (r-db)/SQRT((r-db)**2+add)**3 00557 charg = -charg*0.5_dp*fact 00558 RETURN 00559 END IF 00560 ! Z - Z. 00561 IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMz .AND. m2 == CLMz ) THEN 00562 dzdz = & 00563 +(r+da-db)/SQRT((r+da-db)**2+add)**3 + (r-da+db)/SQRT((r-da+db)**2+add)**3 & 00564 -(r-da-db)/SQRT((r-da-db)**2+add)**3 - (r+da+db)/SQRT((r+da+db)**2+add)**3 00565 charg = -dzdz*0.25_dp*fact 00566 RETURN 00567 END IF 00568 ! X - X 00569 IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMp .AND. m2 == CLMp ) THEN 00570 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 00571 charg = -dxdx*0.25_dp*fact 00572 RETURN 00573 END IF 00574 ! Q - ZZ 00575 IF (l1 == 0 .AND. l2 == 2 .AND. m2 == CLMzz) THEN 00576 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 00577 charg = -qqzz*0.25_dp*fact 00578 RETURN 00579 END IF 00580 ! Q - XX 00581 IF (l1 == 0 .AND. l2 == 2 .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN 00582 qqxx = -r/SQRT(r**2+add)**3 + r/SQRT(r**2+add+db**2)**3 00583 charg = -qqxx*0.5_dp*fact 00584 RETURN 00585 END IF 00586 ! Z - ZZ 00587 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. m2 == CLMzz ) THEN 00588 dzqzz = & 00589 +(r-da-db)/SQRT((r-da-db)**2+add)**3 - 2.0_dp*(r-da)/SQRT((r-da)**2+add)**3 & 00590 +(r-da+db)/SQRT((r-da+db)**2+add)**3 - (r+da-db)/SQRT((r+da-db)**2+add)**3 & 00591 +2.0_dp*(r+da)/SQRT((r+da)**2+add)**3 - (r+da+db)/SQRT((r+da+db)**2+add)**3 00592 charg = -dzqzz*0.125_dp*fact 00593 RETURN 00594 END IF 00595 ! Z - XX 00596 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. ( m2 == CLMyy .OR. m2 == CLMxx )) THEN 00597 dzqxx = & 00598 +(r+da)/SQRT((r+da)**2+add)**3 - (r+da)/SQRT((r+da)**2+add+db**2)**3 & 00599 -(r-da)/SQRT((r-da)**2+add)**3 + (r-da)/SQRT((r-da)**2+add+db**2)**3 00600 charg = -dzqxx*0.25_dp*fact 00601 RETURN 00602 END IF 00603 ! ZZ - ZZ 00604 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. m2 == CLMzz ) THEN 00605 zzzz = & 00606 +(r-da-db)/SQRT((r-da-db)**2+add)**3 + (r+da+db)/SQRT((r+da+db)**2+add)**3 & 00607 +(r-da+db)/SQRT((r-da+db)**2+add)**3 + (r+da-db)/SQRT((r+da-db)**2+add)**3 00608 xyxy = & 00609 +(r-da)/SQRT((r-da)**2+add)**3 + (r+da)/SQRT((r+da)**2+add)**3 & 00610 +(r-db)/SQRT((r-db)**2+add)**3 + (r+db)/SQRT((r+db)**2+add)**3 & 00611 -2.0_dp*r/SQRT(r**2+add)**3 00612 charg = -(zzzz*0.0625_dp - xyxy*0.125_dp)*fact 00613 RETURN 00614 END IF 00615 ! ZZ - XX 00616 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. ( m2 == CLMxx .OR. m2 == CLMyy )) THEN 00617 zzzz = & 00618 -(r+da)/SQRT((r+da)**2+add)**3 + (r+da)/SQRT((r+da)**2+db**2+add)**3 & 00619 -(r-da)/SQRT((r-da)**2+add)**3 + (r-da)/SQRT((r-da)**2+db**2+add)**3 00620 xyxy = r/SQRT(r**2+db**2+add)**3 - r/SQRT(r**2+add)**3 00621 charg = -(zzzz*0.125_dp - xyxy*0.25_dp)*fact 00622 RETURN 00623 END IF 00624 ! X - ZX 00625 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMp .AND. m2 == CLMzp ) THEN 00626 db = db/2.0_dp 00627 dxqxz = & 00628 -(r-db)/SQRT((r-db)**2+(da-db)**2+add)**3 + (r+db)/SQRT((r+db)**2+(da-db)**2+add)**3 & 00629 +(r-db)/SQRT((r-db)**2+(da+db)**2+add)**3 - (r+db)/SQRT((r+db)**2+(da+db)**2+add)**3 00630 charg = -dxqxz*0.25_dp*fact 00631 RETURN 00632 END IF 00633 ! ZX - ZX 00634 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzp .AND. m2 == CLMzp ) THEN 00635 da = da / 2.0_dp 00636 db = db / 2.0_dp 00637 qxzqxz = & 00638 +(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 & 00639 -(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 & 00640 -(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 & 00641 +(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 00642 charg = -qxzqxz*0.125_dp*fact 00643 RETURN 00644 END IF 00645 ! XX - XX 00646 IF (l1 == 2 .AND. l2 == 2 .AND. (((m1 == CLMyy) .AND. (m2 == CLMyy)).OR.((m1 == CLMxx) .AND. (m2 == CLMxx)))) THEN 00647 qxxqxx = & 00648 + 2.0_dp*r/SQRT(r**2+add)**3 + r/SQRT(r**2+(da-db)**2+add)**3 & 00649 + r/SQRT(r**2+(da+db)**2+add)**3 - 2.0_dp*r/SQRT(r**2+da**2+add)**3 & 00650 - 2.0_dp*r/SQRT(r**2+db**2+add)**3 00651 charg = -qxxqxx*0.125_dp*fact 00652 RETURN 00653 END IF 00654 ! XX - YY 00655 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMyy .AND. m2 == CLMxx) THEN 00656 qxxqyy = & 00657 + r/SQRT(r**2+add)**3 - r/SQRT(r**2+da**2+add)**3 & 00658 - r/SQRT(r**2+db**2+add)**3 + r/SQRT(r**2+da**2+db**2+add)**3 00659 charg = -qxxqyy*0.25_dp*fact 00660 RETURN 00661 END IF 00662 ! XY - XY 00663 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMxy .AND. m2 == CLMxy) THEN 00664 qxxqxx = & 00665 + 2.0_dp*r/SQRT(r**2+add)**3 + r/SQRT(r**2+(da-db)**2+add)**3 & 00666 + r/SQRT(r**2+(da+db)**2+add)**3 - 2.0_dp*r/SQRT(r**2+da**2+add)**3 & 00667 - 2.0_dp*r/SQRT(r**2+db**2+add)**3 00668 qxxqyy = & 00669 + r/SQRT(r**2+add)**3 - r/SQRT(r**2+da**2+add)**3 & 00670 - r/SQRT(r**2+db**2+add)**3 + r/SQRT(r**2+da**2+db**2+add)**3 00671 charg = -0.5_dp * (qxxqxx*0.125_dp - qxxqyy*0.25_dp) * fact 00672 RETURN 00673 END IF 00674 ! We should NEVER reach this point 00675 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 00676 END FUNCTION dcharg_int_nri 00677 00678 ! ***************************************************************************** 00694 FUNCTION dcharg_int_nri_fs(r, l1_i, l2_i, m1_i, m2_i, da_i, db_i, add0, fact_screen, error) RESULT(charg) 00695 REAL(KIND=dp), INTENT(in) :: r 00696 INTEGER, INTENT(in) :: l1_i, l2_i, m1_i, m2_i 00697 REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen 00698 TYPE(cp_error_type), INTENT(inout) :: error 00699 REAL(KIND=dp) :: charg 00700 00701 CHARACTER(len=*), PARAMETER :: routineN = 'dcharg_int_nri_fs', 00702 routineP = moduleN//':'//routineN 00703 00704 INTEGER :: l1, l2, m1, m2 00705 LOGICAL :: failure 00706 REAL(KIND=dp) :: add, da, db, dxdx, dxqxz, 00707 dzdz, dzqxx, dzqzz, fact, 00708 qqxx, qqzz, qxxqxx, qxxqyy, 00709 qxzqxz, xyxy, zzzz 00710 00711 failure = .FALSE. 00712 ! Computing only Integral Derivatives 00713 IF ( l1_i < l2_i) THEN 00714 l1 = l1_i 00715 l2 = l2_i 00716 m1 = m1_i 00717 m2 = m2_i 00718 da = da_i 00719 db = db_i 00720 fact = 1.0_dp 00721 ELSE IF ( l1_i > l2_i ) THEN 00722 l1 = l2_i 00723 l2 = l1_i 00724 m1 = m2_i 00725 m2 = m1_i 00726 da = db_i 00727 db = da_i 00728 fact = (-1.0_dp)**(l1+l2) 00729 ELSE IF ( l1_i == l2_i ) THEN 00730 l1 = l1_i 00731 l2 = l2_i 00732 IF ( m1_i <= m2_i ) THEN 00733 m1 = m1_i 00734 m2 = m2_i 00735 da = da_i 00736 db = db_i 00737 ELSE 00738 m1 = m2_i 00739 m2 = m1_i 00740 da = db_i 00741 db = da_i 00742 END IF 00743 fact = 1.0_dp 00744 END IF 00745 charg = 0.0_dp 00746 add = add0 * fact_screen 00747 ! The 0.5 factor handles the derivative of the SQRT 00748 fact = fact * 0.5_dp 00749 ! Q - Q. 00750 IF (l1 == 0 .AND. l2 == 0) THEN 00751 charg = add0/SQRT(r**2+add)**3 00752 charg = -charg*fact 00753 RETURN 00754 END IF 00755 ! Q - Z. 00756 IF (l1 == 0 .AND. l2 == 1 .AND. m2 == CLMz) THEN 00757 charg = add0/SQRT((r+db)**2+add)**3 - add0/SQRT((r-db)**2+add)**3 00758 charg = -charg*0.5_dp*fact 00759 RETURN 00760 END IF 00761 ! Z - Z. 00762 IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMz .AND. m2 == CLMz ) THEN 00763 dzdz = & 00764 +add0/SQRT((r+da-db)**2+add)**3 + add0/SQRT((r-da+db)**2+add)**3 & 00765 -add0/SQRT((r-da-db)**2+add)**3 - add0/SQRT((r+da+db)**2+add)**3 00766 charg = -dzdz*0.25_dp*fact 00767 RETURN 00768 END IF 00769 ! X - X 00770 IF (l1 == 1 .AND. l2 == 1 .AND. m1 == CLMp .AND. m2 == CLMp ) THEN 00771 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 00772 charg = -dxdx*0.25_dp*fact 00773 RETURN 00774 END IF 00775 ! Q - ZZ 00776 IF (l1 == 0 .AND. l2 == 2 .AND. m2 == CLMzz) THEN 00777 qqzz = add0/SQRT((r-db)**2+add)**3 - 2.0_dp*add0/SQRT(r**2+add)**3 + add0/SQRT((r+db)**2+add)**3 00778 charg = -qqzz*0.25_dp*fact 00779 RETURN 00780 END IF 00781 ! Q - XX 00782 IF (l1 == 0 .AND. l2 == 2 .AND. (m2 == CLMyy .OR. m2 == CLMxx)) THEN 00783 qqxx = -add0/SQRT(r**2+add)**3 + add0/SQRT(r**2+add+db**2)**3 00784 charg = -qqxx*0.5_dp*fact 00785 RETURN 00786 END IF 00787 ! Z - ZZ 00788 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. m2 == CLMzz ) THEN 00789 dzqzz = & 00790 +add0/SQRT((r-da-db)**2+add)**3 - 2.0_dp*add0/SQRT((r-da)**2+add)**3 & 00791 +add0/SQRT((r-da+db)**2+add)**3 - add0/SQRT((r+da-db)**2+add)**3 & 00792 +2.0_dp*add0/SQRT((r+da)**2+add)**3 - add0/SQRT((r+da+db)**2+add)**3 00793 charg = -dzqzz*0.125_dp*fact 00794 RETURN 00795 END IF 00796 ! Z - XX 00797 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMz .AND. ( m2 == CLMyy .OR. m2 == CLMxx )) THEN 00798 dzqxx = & 00799 +add0/SQRT((r+da)**2+add)**3 - add0/SQRT((r+da)**2+add+db**2)**3 & 00800 -add0/SQRT((r-da)**2+add)**3 + add0/SQRT((r-da)**2+add+db**2)**3 00801 charg = -dzqxx*0.25_dp*fact 00802 RETURN 00803 END IF 00804 ! ZZ - ZZ 00805 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. m2 == CLMzz ) THEN 00806 zzzz = & 00807 +add0/SQRT((r-da-db)**2+add)**3 + add0/SQRT((r+da+db)**2+add)**3 & 00808 +add0/SQRT((r-da+db)**2+add)**3 + add0/SQRT((r+da-db)**2+add)**3 00809 xyxy = & 00810 +add0/SQRT((r-da)**2+add)**3 + add0/SQRT((r+da)**2+add)**3 & 00811 +add0/SQRT((r-db)**2+add)**3 + add0/SQRT((r+db)**2+add)**3 & 00812 -2.0_dp*add0/SQRT(r**2+add)**3 00813 charg = -(zzzz*0.0625_dp - xyxy*0.125_dp)*fact 00814 RETURN 00815 END IF 00816 ! ZZ - XX 00817 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzz .AND. ( m2 == CLMxx .OR. m2 == CLMyy )) THEN 00818 zzzz = & 00819 -add0/SQRT((r+da)**2+add)**3 + add0/SQRT((r+da)**2+db**2+add)**3 & 00820 -add0/SQRT((r-da)**2+add)**3 + add0/SQRT((r-da)**2+db**2+add)**3 00821 xyxy =add0/SQRT(r**2+db**2+add)**3 - add0/SQRT(r**2+add)**3 00822 charg = -(zzzz*0.125_dp - xyxy*0.25_dp)*fact 00823 RETURN 00824 END IF 00825 ! X - ZX 00826 IF (l1 == 1 .AND. l2 == 2 .AND. m1 == CLMp .AND. m2 == CLMzp ) THEN 00827 db = db/2.0_dp 00828 dxqxz = & 00829 -add0/SQRT((r-db)**2+(da-db)**2+add)**3 + add0/SQRT((r+db)**2+(da-db)**2+add)**3 & 00830 +add0/SQRT((r-db)**2+(da+db)**2+add)**3 - add0/SQRT((r+db)**2+(da+db)**2+add)**3 00831 charg = -dxqxz*0.25_dp*fact 00832 RETURN 00833 END IF 00834 ! ZX - ZX 00835 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMzp .AND. m2 == CLMzp ) THEN 00836 da = da / 2.0_dp 00837 db = db / 2.0_dp 00838 qxzqxz = & 00839 +add0/SQRT((r+da-db)**2+(da-db)**2+add)**3 - add0/SQRT((r+da+db)**2+(da-db)**2+add)**3 & 00840 -add0/SQRT((r-da-db)**2+(da-db)**2+add)**3 + add0/SQRT((r-da+db)**2+(da-db)**2+add)**3 & 00841 -add0/SQRT((r+da-db)**2+(da+db)**2+add)**3 + add0/SQRT((r+da+db)**2+(da+db)**2+add)**3 & 00842 +add0/SQRT((r-da-db)**2+(da+db)**2+add)**3 - add0/SQRT((r-da+db)**2+(da+db)**2+add)**3 00843 charg = -qxzqxz*0.125_dp*fact 00844 RETURN 00845 END IF 00846 ! XX - XX 00847 IF (l1 == 2 .AND. l2 == 2 .AND. (((m1 == CLMyy) .AND. (m2 == CLMyy)).OR.((m1 == CLMxx) .AND. (m2 == CLMxx)))) THEN 00848 qxxqxx = & 00849 + 2.0_dp*add0/SQRT(r**2+add)**3 + add0/SQRT(r**2+(da-db)**2+add)**3 & 00850 + add0/SQRT(r**2+(da+db)**2+add)**3 - 2.0_dp*add0/SQRT(r**2+da**2+add)**3 & 00851 - 2.0_dp*add0/SQRT(r**2+db**2+add)**3 00852 charg = -qxxqxx*0.125_dp*fact 00853 RETURN 00854 END IF 00855 ! XX - YY 00856 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMyy .AND. m2 == CLMxx) THEN 00857 qxxqyy = & 00858 + add0/SQRT(r**2+add)**3 - add0/SQRT(r**2+da**2+add)**3 & 00859 - add0/SQRT(r**2+db**2+add)**3 + add0/SQRT(r**2+da**2+db**2+add)**3 00860 charg = -qxxqyy*0.25_dp*fact 00861 RETURN 00862 END IF 00863 ! XY - XY 00864 IF (l1 == 2 .AND. l2 == 2 .AND. m1 == CLMxy .AND. m2 == CLMxy) THEN 00865 qxxqxx = & 00866 + 2.0_dp*add0/SQRT(r**2+add)**3 + add0/SQRT(r**2+(da-db)**2+add)**3 & 00867 + add0/SQRT(r**2+(da+db)**2+add)**3 - 2.0_dp*add0/SQRT(r**2+da**2+add)**3 & 00868 - 2.0_dp*add0/SQRT(r**2+db**2+add)**3 00869 qxxqyy = & 00870 + add0/SQRT(r**2+add)**3 - add0/SQRT(r**2+da**2+add)**3 & 00871 - add0/SQRT(r**2+db**2+add)**3 + add0/SQRT(r**2+da**2+db**2+add)**3 00872 charg = -0.5_dp * (qxxqxx*0.125_dp - qxxqyy*0.25_dp) * fact 00873 RETURN 00874 END IF 00875 ! We should NEVER reach this point 00876 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 00877 END FUNCTION dcharg_int_nri_fs 00878 00879 ! ***************************************************************************** 00895 FUNCTION ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control, & 00896 se_int_screen, itype, error) RESULT(res) 00897 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00898 INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic 00899 REAL(KIND=dp), INTENT(IN) :: r 00900 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00901 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 00902 INTEGER, INTENT(IN) :: itype 00903 TYPE(cp_error_type), INTENT(inout) :: error 00904 REAL(KIND=dp) :: res 00905 00906 CHARACTER(len=*), PARAMETER :: routineN = 'ijkl_d', 00907 routineP = moduleN//':'//routineN 00908 00909 res = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,& 00910 se_int_control%integral_screening, se_int_control%shortrange,& 00911 se_int_control%pc_coulomb_int, se_int_control%max_multipole,& 00912 itype, charg_int_ri, error) 00913 00914 ! If only the shortrange component is requested we can skip the rest 00915 IF ((.NOT.se_int_control%pc_coulomb_int).AND.(itype/=do_method_pchg)) THEN 00916 ! Handle the 1/r^3 term, this term is ALWAYS false for KDSO-D integrals 00917 IF (se_int_control%shortrange.AND.se_int_control%do_ewald_r3) THEN 00918 res = res - ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r,& 00919 itype, charg_int_3, error) 00920 END IF 00921 END IF 00922 END FUNCTION ijkl_d 00923 00924 ! ***************************************************************************** 00940 FUNCTION d_ijkl_d(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_control,& 00941 se_int_screen, itype, error) RESULT(res) 00942 TYPE(semi_empirical_type), POINTER :: sepi, sepj 00943 INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic 00944 REAL(KIND=dp), INTENT(IN) :: r 00945 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 00946 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 00947 INTEGER, INTENT(IN) :: itype 00948 TYPE(cp_error_type), INTENT(inout) :: error 00949 REAL(KIND=dp) :: res 00950 00951 CHARACTER(len=*), PARAMETER :: routineN = 'd_ijkl_d', 00952 routineP = moduleN//':'//routineN 00953 00954 REAL(KIND=dp) :: dfs, srd 00955 00956 IF (se_int_control%integral_screening==do_se_IS_kdso_d) THEN 00957 res = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,& 00958 se_int_control%integral_screening, .FALSE.,& 00959 se_int_control%pc_coulomb_int, se_int_control%max_multipole,& 00960 itype, dcharg_int_ri, error) 00961 00962 IF (.NOT.se_int_control%pc_coulomb_int) THEN 00963 dfs = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,& 00964 se_int_control%integral_screening, .FALSE., .FALSE.,& 00965 se_int_control%max_multipole, itype, dcharg_int_ri_fs, error) 00966 res = res + dfs*se_int_screen%dft 00967 00968 ! In case we need the shortrange part we have to evaluate an additional derivative 00969 ! to handle the derivative of the Tapering term 00970 IF (se_int_control%shortrange) THEN 00971 srd = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,& 00972 se_int_control%integral_screening, .FALSE., .TRUE.,& 00973 se_int_control%max_multipole, itype, dcharg_int_ri, error) 00974 res = res - srd 00975 END IF 00976 END IF 00977 ELSE 00978 res = ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,& 00979 se_int_control%integral_screening, se_int_control%shortrange,& 00980 se_int_control%pc_coulomb_int, se_int_control%max_multipole,& 00981 itype, dcharg_int_ri, error) 00982 END IF 00983 00984 ! If only the shortrange component is requested we can skip the rest 00985 IF ((.NOT.se_int_control%pc_coulomb_int).AND.(itype/=do_method_pchg)) THEN 00986 ! Handle the 1/r^3 term, this term is ALWAYS false for KDSO-D integrals 00987 IF (se_int_control%shortrange.AND.se_int_control%do_ewald_r3) THEN 00988 res = res - ijkl_low_3(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r,& 00989 itype, dcharg_int_3, error) 00990 END IF 00991 END IF 00992 00993 END FUNCTION d_ijkl_d 00994 00995 ! ***************************************************************************** 01011 FUNCTION ijkl_d_low(sepi, sepj, ij, kl, li, lj, lk, ll, ic, r, se_int_screen,& 01012 iscreen, shortrange, pc_coulomb_int, max_multipole, itype, eval, error) RESULT(res) 01013 TYPE(semi_empirical_type), POINTER :: sepi, sepj 01014 INTEGER, INTENT(IN) :: ij, kl, li, lj, lk, ll, ic 01015 REAL(KIND=dp), INTENT(IN) :: r 01016 TYPE(se_int_screen_type), INTENT(IN) :: se_int_screen 01017 INTEGER, INTENT(IN) :: iscreen 01018 LOGICAL, INTENT(IN) :: shortrange, pc_coulomb_int 01019 INTEGER, INTENT(IN) :: max_multipole, itype 01020 REAL(KIND=dp) :: eval 01021 TYPE(cp_error_type), INTENT(inout) :: error 01022 REAL(KIND=dp) :: res 01023 01024 CHARACTER(len=*), PARAMETER :: routineN = 'ijkl_d_low', 01025 routineP = moduleN//':'//routineN 01026 01027 INTEGER :: l1, l1max, l1min, l2, l2max, 01028 l2min, lij, lkl, lmin, m, mm 01029 REAL(KIND=dp) :: add, ccc, chrg, dij, dkl, 01030 fact_screen, pij, pkl, s1, sum 01031 01032 l1min = ABS (li-lj) 01033 l1max = li + lj 01034 lij = indexb(li+1, lj+1) 01035 l2min = ABS (lk-ll) 01036 l2max = lk + ll 01037 lkl = indexb(lk+1, ll+1) 01038 01039 l1max = MIN (l1max, 2) 01040 l1min = MIN (l1min, 2) 01041 l2max = MIN (l2max, 2) 01042 l2min = MIN (l2min, 2) 01043 sum = 0.0_dp 01044 dij = 0.0_dp 01045 dkl = 0.0_dp 01046 fact_screen = 1.0_dp 01047 IF (.NOT.pc_coulomb_int) THEN 01048 IF (iscreen==do_se_IS_kdso_d) fact_screen = se_int_screen%ft 01049 ! Standard value of the integral 01050 DO l1 = l1min, l1max 01051 IF (l1 == 0) THEN 01052 IF (lij == 1) THEN 01053 pij = sepi%ko(1) 01054 IF (ic == 1) THEN 01055 pij = sepi%ko(9) 01056 END IF 01057 ELSE IF (lij == 3) THEN 01058 pij = sepi%ko(7) 01059 ELSE IF (lij == 6) THEN 01060 pij = sepi%ko(8) 01061 END IF 01062 ELSE 01063 dij = sepi%cs(lij) 01064 pij = sepi%ko(lij) 01065 END IF 01066 ! 01067 DO l2 = l2min, l2max 01068 IF (l2 == 0) THEN 01069 IF (lkl == 1) THEN 01070 pkl = sepj%ko(1) 01071 IF (ic == 2) THEN 01072 pkl = sepj%ko(9) 01073 END IF 01074 ELSE IF (lkl == 3) THEN 01075 pkl = sepj%ko(7) 01076 ELSE IF (lkl == 6) THEN 01077 pkl = sepj%ko(8) 01078 END IF 01079 ELSE 01080 dkl = sepj%cs(lkl) 01081 pkl = sepj%ko(lkl) 01082 END IF 01083 IF (itype==do_method_pchg) THEN 01084 add = 0.0_dp 01085 ELSE 01086 add = (pij+pkl) ** 2 01087 END IF 01088 lmin = MIN(l1, l2) 01089 s1 = 0.0_dp 01090 DO m = -lmin, lmin 01091 ccc = clm_d(ij, l1, m) * clm_d(kl, l2, m) 01092 IF (ABS(ccc) > EPSILON(0.0_dp)) THEN 01093 mm = ABS (m) 01094 chrg = eval(r, l1, l2, mm, dij, dkl, add, fact_screen, error) 01095 s1 = s1 + chrg * ccc 01096 END IF 01097 END DO 01098 sum = sum + s1 01099 END DO 01100 END DO 01101 res = sum 01102 END IF 01103 ! Shortrange: Possibly computes pure Coulomb and subtract from the original integral valeu 01104 IF (shortrange.OR.pc_coulomb_int) THEN 01105 sum = 0.0_dp 01106 dij = 0.0_dp 01107 dkl = 0.0_dp 01108 add = 0.0_dp 01109 fact_screen = 0.0_dp 01110 DO l1 = l1min, l1max 01111 IF (l1 > max_multipole) CYCLE 01112 IF (l1 /= 0) THEN 01113 dij = sepi%cs(lij) 01114 END IF 01115 ! 01116 DO l2 = l2min, l2max 01117 IF (l2 > max_multipole) CYCLE 01118 IF (l2 /= 0) THEN 01119 dkl = sepj%cs(lkl) 01120 END IF 01121 lmin = MIN(l1, l2) 01122 s1 = 0.0_dp 01123 DO m = -lmin, lmin 01124 ccc = clm_d(ij, l1, m) * clm_d(kl, l2, m) 01125 IF (ABS(ccc) > EPSILON(0.0_dp)) THEN 01126 mm = ABS (m) 01127 chrg = eval(r, l1, l2, mm, dij, dkl, add, fact_screen, error) 01128 s1 = s1 + chrg * ccc 01129 END IF 01130 END DO 01131 sum = sum + s1 01132 END DO 01133 END DO 01134 IF (pc_coulomb_int) res = sum 01135 IF (shortrange) res = res - sum 01136 END IF 01137 END FUNCTION ijkl_d_low 01138 01139 ! ***************************************************************************** 01152 FUNCTION charg_int_ri(r, l1_i, l2_i, m, da_i, db_i, add0, fact_screen, error) RESULT(charg) 01153 REAL(KIND=dp), INTENT(in) :: r 01154 INTEGER, INTENT(in) :: l1_i, l2_i, m 01155 REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen 01156 TYPE(cp_error_type), INTENT(inout) :: error 01157 REAL(KIND=dp) :: charg 01158 01159 CHARACTER(len=*), PARAMETER :: routineN = 'charg_int_ri', 01160 routineP = moduleN//':'//routineN 01161 01162 INTEGER :: l1, l2 01163 LOGICAL :: failure 01164 REAL(KIND=dp) :: aa, ab, add, da, db, dxdx, 01165 dxqxz, dzdz, dzqzz, fact, 01166 qqzz, qxzqxz, xyxy, zzzz 01167 01168 failure = .FALSE. 01169 IF ( l1_i <= l2_i) THEN 01170 l1 = l1_i 01171 l2 = l2_i 01172 da = da_i 01173 db = db_i 01174 fact = 1.0_dp 01175 ELSE IF ( l1_i > l2_i ) THEN 01176 l1 = l2_i 01177 l2 = l1_i 01178 da = db_i 01179 db = da_i 01180 fact = (-1.0_dp)**(l1+l2) 01181 END IF 01182 charg = 0.0_dp 01183 add = add0*fact_screen 01184 ! Q - Q. 01185 IF (l1 == 0 .AND. l2 == 0) THEN 01186 charg = fact/SQRT(r**2+add) 01187 RETURN 01188 END IF 01189 ! Q - Z. 01190 IF (l1 == 0 .AND. l2 == 1) THEN 01191 charg = 1.0_dp/SQRT((r+db)**2+add) - 1.0_dp/SQRT((r-db)**2+add) 01192 charg = charg*0.5_dp*fact 01193 RETURN 01194 END IF 01195 ! Z - Z. 01196 IF (l1 == 1 .AND. l2 == 1 .AND. m == 0) THEN 01197 dzdz = & 01198 +1.0_dp/SQRT((r+da-db)**2+add) + 1.0_dp/SQRT((r-da+db)**2+add) & 01199 -1.0_dp/SQRT((r-da-db)**2+add) - 1.0_dp/SQRT((r+da+db)**2+add) 01200 charg = dzdz*0.25_dp*fact 01201 RETURN 01202 END IF 01203 ! X - X 01204 IF (l1 == 1 .AND. l2 == 1 .AND. m == 1) THEN 01205 dxdx = 2.0_dp/SQRT(r**2+(da-db)**2+add) - 2.0_dp/SQRT(r**2+(da+db)**2+add) 01206 charg = dxdx*0.25_dp*fact 01207 RETURN 01208 END IF 01209 ! Q - ZZ 01210 IF (l1 == 0 .AND. l2 == 2) THEN 01211 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) 01212 charg = qqzz*0.25_dp*fact 01213 RETURN 01214 END IF 01215 ! Z - ZZ 01216 IF (l1 == 1 .AND. l2 == 2 .AND. m == 0) THEN 01217 dzqzz = & 01218 +1.0_dp/SQRT((r-da-db)**2+add) - 2.0_dp/SQRT((r-da)**2+db**2+add) & 01219 +1.0_dp/SQRT((r+db-da)**2+add) - 1.0_dp/SQRT((r-db+da)**2+add) & 01220 +2.0_dp/SQRT((r+da)**2+db**2+add) - 1.0_dp/SQRT((r+da+db)**2+add) 01221 charg = dzqzz*0.125_dp*fact 01222 RETURN 01223 END IF 01224 ! ZZ - ZZ 01225 IF (l1 == 2 .AND. l2 == 2 .AND. m == 0) THEN 01226 zzzz = & 01227 +1.0_dp/SQRT((r-da-db)**2+add) + 1.0_dp/SQRT((r+da+db)**2+add) & 01228 +1.0_dp/SQRT((r-da+db)**2+add) + 1.0_dp/SQRT((r+da-db)**2+add) & 01229 -2.0_dp/SQRT((r-da)**2+db**2+add) - 2.0_dp/SQRT((r-db)**2+da**2+add) & 01230 -2.0_dp/SQRT((r+da)**2+db**2+add) - 2.0_dp/SQRT((r+db)**2+da**2+add) & 01231 +2.0_dp/SQRT(r**2+(da-db)**2+add) + 2.0_dp/SQRT(r**2+(da+db)**2+add) 01232 xyxy = & 01233 +4.0_dp/SQRT(r**2+(da-db)**2+add) + 4.0_dp/SQRT(r**2+(da+db)**2+add) & 01234 -8.0_dp/SQRT(r**2+da**2+db**2+add) 01235 charg = (zzzz*0.0625_dp - xyxy*0.015625_dp)*fact 01236 RETURN 01237 END IF 01238 ! X - ZX 01239 IF (l1 == 1 .AND. l2 == 2 .AND. m == 1) THEN 01240 ab = db/SQRT(2.0_dp) 01241 dxqxz = & 01242 -2.0_dp/SQRT((r-ab)**2+(da-ab)**2+add) + 2.0_dp/SQRT((r+ab)**2+(da-ab)**2+add) & 01243 +2.0_dp/SQRT((r-ab)**2+(da+ab)**2+add) - 2.0_dp/SQRT((r+ab)**2+(da+ab)**2+add) 01244 charg = dxqxz*0.125_dp*fact 01245 RETURN 01246 END IF 01247 ! ZX - ZX 01248 IF (l1 == 2 .AND. l2 == 2 .AND. m == 1) THEN 01249 aa = da/SQRT (2.0_dp) 01250 ab = db/SQRT (2.0_dp) 01251 qxzqxz = & 01252 +2.0_dp/SQRT((r+aa-ab)**2+(aa-ab)**2+add) - 2.0_dp/SQRT((r+aa+ab)**2+(aa-ab)**2+add) & 01253 -2.0_dp/SQRT((r-aa-ab)**2+(aa-ab)**2+add) + 2.0_dp/SQRT((r-aa+ab)**2+(aa-ab)**2+add) & 01254 -2.0_dp/SQRT((r+aa-ab)**2+(aa+ab)**2+add) + 2.0_dp/SQRT((r+aa+ab)**2+(aa+ab)**2+add) & 01255 +2.0_dp/SQRT((r-aa-ab)**2+(aa+ab)**2+add) - 2.0_dp/SQRT((r-aa+ab)**2+(aa+ab)**2+add) 01256 charg = qxzqxz*0.0625_dp*fact 01257 RETURN 01258 END IF 01259 ! XX - XX 01260 IF (l1 == 2 .AND. l2 == 2 .AND. m == 2) THEN 01261 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) 01262 charg = xyxy*0.0625_dp*fact 01263 RETURN 01264 END IF 01265 ! We should NEVER reach this point 01266 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 01267 END FUNCTION charg_int_ri 01268 01269 ! ***************************************************************************** 01283 FUNCTION dcharg_int_ri(r, l1_i, l2_i, m, da_i, db_i, add0, fact_screen, error) RESULT(charg) 01284 REAL(KIND=dp), INTENT(in) :: r 01285 INTEGER, INTENT(in) :: l1_i, l2_i, m 01286 REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen 01287 TYPE(cp_error_type), INTENT(inout) :: error 01288 REAL(KIND=dp) :: charg 01289 01290 CHARACTER(len=*), PARAMETER :: routineN = 'dcharg_int_ri', 01291 routineP = moduleN//':'//routineN 01292 01293 INTEGER :: l1, l2 01294 LOGICAL :: failure 01295 REAL(KIND=dp) :: aa, ab, add, da, db, dxdx, 01296 dxqxz, dzdz, dzqzz, fact, 01297 qqzz, qxzqxz, xyxy, zzzz 01298 01299 failure = .FALSE. 01300 IF ( l1_i <= l2_i) THEN 01301 l1 = l1_i 01302 l2 = l2_i 01303 da = da_i 01304 db = db_i 01305 fact = 1.0_dp 01306 ELSE IF ( l1_i > l2_i ) THEN 01307 l1 = l2_i 01308 l2 = l1_i 01309 da = db_i 01310 db = da_i 01311 fact = (-1.0_dp)**(l1+l2) 01312 END IF 01313 charg = 0.0_dp 01314 add = add0*fact_screen 01315 ! Q - Q. 01316 IF (l1 == 0 .AND. l2 == 0) THEN 01317 charg = r/SQRT(r**2+add)**3 01318 charg = -fact*charg 01319 RETURN 01320 END IF 01321 ! Q - Z. 01322 IF (l1 == 0 .AND. l2 == 1) THEN 01323 charg = (r+db)/SQRT((r+db)**2+add)**3 - (r-db)/SQRT((r-db)**2+add)**3 01324 charg = -charg*0.5_dp*fact 01325 RETURN 01326 END IF 01327 ! Z - Z. 01328 IF (l1 == 1 .AND. l2 == 1 .AND. m == 0) THEN 01329 dzdz = & 01330 +(r+da-db)/SQRT((r+da-db)**2+add)**3 + (r-da+db)/SQRT((r-da+db)**2+add)**3 & 01331 -(r-da-db)/SQRT((r-da-db)**2+add)**3 - (r+da+db)/SQRT((r+da+db)**2+add)**3 01332 charg = -dzdz*0.25_dp*fact 01333 RETURN 01334 END IF 01335 ! X - X 01336 IF (l1 == 1 .AND. l2 == 1 .AND. m == 1) THEN 01337 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 01338 charg = -dxdx*0.25_dp*fact 01339 RETURN 01340 END IF 01341 ! Q - ZZ 01342 IF (l1 == 0 .AND. l2 == 2) THEN 01343 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 01344 charg = -qqzz*0.25_dp*fact 01345 RETURN 01346 END IF 01347 ! Z - ZZ 01348 IF (l1 == 1 .AND. l2 == 2 .AND. m == 0) THEN 01349 dzqzz = & 01350 +(r-da-db)/SQRT((r-da-db)**2+add)**3 - 2.0_dp*(r-da)/SQRT((r-da)**2+db**2+add)**3 & 01351 +(r+db-da)/SQRT((r+db-da)**2+add)**3 - (r-db+da)/SQRT((r-db+da)**2+add)**3 & 01352 +2.0_dp*(r+da)/SQRT((r+da)**2+db**2+add)**3 - (r+da+db)/SQRT((r+da+db)**2+add)**3 01353 charg = -dzqzz*0.125_dp*fact 01354 RETURN 01355 END IF 01356 ! ZZ - ZZ 01357 IF (l1 == 2 .AND. l2 == 2 .AND. m == 0) THEN 01358 zzzz = & 01359 +(r-da-db)/SQRT((r-da-db)**2+add)**3 + (r+da+db)/SQRT((r+da+db)**2+add)**3 & 01360 +(r-da+db)/SQRT((r-da+db)**2+add)**3 + (r+da-db)/SQRT((r+da-db)**2+add)**3 & 01361 -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 & 01362 -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 & 01363 +2.0_dp*r/SQRT(r**2+(da-db)**2+add)**3 + 2.0_dp*r/SQRT(r**2+(da+db)**2+add)**3 01364 xyxy = & 01365 +4.0_dp*r/SQRT(r**2+(da-db)**2+add)**3 + 4.0_dp*r/SQRT(r**2+(da+db)**2+add)**3 & 01366 -8.0_dp*r/SQRT(r**2+da**2+db**2+add)**3 01367 charg = -(zzzz*0.0625_dp - xyxy*0.015625_dp)*fact 01368 RETURN 01369 END IF 01370 ! X - ZX 01371 IF (l1 == 1 .AND. l2 == 2 .AND. m == 1) THEN 01372 ab = db/SQRT(2.0_dp) 01373 dxqxz = & 01374 -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 & 01375 +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 01376 charg = -dxqxz*0.125_dp*fact 01377 RETURN 01378 END IF 01379 ! ZX - ZX 01380 IF (l1 == 2 .AND. l2 == 2 .AND. m == 1) THEN 01381 aa = da/SQRT (2.0_dp) 01382 ab = db/SQRT (2.0_dp) 01383 qxzqxz = & 01384 +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 & 01385 -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 & 01386 -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 & 01387 +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 01388 charg = -qxzqxz*0.0625_dp*fact 01389 RETURN 01390 END IF 01391 ! XX - XX 01392 IF (l1 == 2 .AND. l2 == 2 .AND. m == 2) THEN 01393 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 01394 charg = -xyxy*0.0625_dp*fact 01395 RETURN 01396 END IF 01397 ! We should NEVER reach this point 01398 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 01399 END FUNCTION dcharg_int_ri 01400 01401 ! ***************************************************************************** 01416 FUNCTION dcharg_int_ri_fs(r, l1_i, l2_i, m, da_i, db_i, add0, fact_screen, error) RESULT(charg) 01417 REAL(KIND=dp), INTENT(in) :: r 01418 INTEGER, INTENT(in) :: l1_i, l2_i, m 01419 REAL(KIND=dp), INTENT(in) :: da_i, db_i, add0, fact_screen 01420 TYPE(cp_error_type), INTENT(inout) :: error 01421 REAL(KIND=dp) :: charg 01422 01423 CHARACTER(len=*), PARAMETER :: routineN = 'dcharg_int_ri_fs', 01424 routineP = moduleN//':'//routineN 01425 01426 INTEGER :: l1, l2 01427 LOGICAL :: failure 01428 REAL(KIND=dp) :: aa, ab, add, da, db, dxdx, 01429 dxqxz, dzdz, dzqzz, fact, 01430 qqzz, qxzqxz, xyxy, zzzz 01431 01432 failure = .FALSE. 01433 IF ( l1_i <= l2_i) THEN 01434 l1 = l1_i 01435 l2 = l2_i 01436 da = da_i 01437 db = db_i 01438 fact = 1.0_dp 01439 ELSE IF ( l1_i > l2_i ) THEN 01440 l1 = l2_i 01441 l2 = l1_i 01442 da = db_i 01443 db = da_i 01444 fact = (-1.0_dp)**(l1+l2) 01445 END IF 01446 charg = 0.0_dp 01447 add = add0*fact_screen 01448 ! The 0.5 factor handles the derivative of the SQRT 01449 fact = fact * 0.5_dp 01450 ! Q - Q. 01451 IF (l1 == 0 .AND. l2 == 0) THEN 01452 charg = add0/SQRT(r**2+add)**3 01453 charg = -fact*charg 01454 RETURN 01455 END IF 01456 ! Q - Z. 01457 IF (l1 == 0 .AND. l2 == 1) THEN 01458 charg = add0/SQRT((r+db)**2+add)**3 - add0/SQRT((r-db)**2+add)**3 01459 charg = -charg*0.5_dp*fact 01460 RETURN 01461 END IF 01462 ! Z - Z. 01463 IF (l1 == 1 .AND. l2 == 1 .AND. m == 0) THEN 01464 dzdz = & 01465 +add0/SQRT((r+da-db)**2+add)**3 + add0/SQRT((r-da+db)**2+add)**3 & 01466 -add0/SQRT((r-da-db)**2+add)**3 - add0/SQRT((r+da+db)**2+add)**3 01467 charg = -dzdz*0.25_dp*fact 01468 RETURN 01469 END IF 01470 ! X - X 01471 IF (l1 == 1 .AND. l2 == 1 .AND. m == 1) THEN 01472 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 01473 charg = -dxdx*0.25_dp*fact 01474 RETURN 01475 END IF 01476 ! Q - ZZ 01477 IF (l1 == 0 .AND. l2 == 2) THEN 01478 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 01479 charg = -qqzz*0.25_dp*fact 01480 RETURN 01481 END IF 01482 ! Z - ZZ 01483 IF (l1 == 1 .AND. l2 == 2 .AND. m == 0) THEN 01484 dzqzz = & 01485 +add0/SQRT((r-da-db)**2+add)**3 - 2.0_dp*add0/SQRT((r-da)**2+db**2+add)**3 & 01486 +add0/SQRT((r+db-da)**2+add)**3 - add0/SQRT((r-db+da)**2+add)**3 & 01487 +2.0_dp*add0/SQRT((r+da)**2+db**2+add)**3 - add0/SQRT((r+da+db)**2+add)**3 01488 charg = -dzqzz*0.125_dp*fact 01489 RETURN 01490 END IF 01491 ! ZZ - ZZ 01492 IF (l1 == 2 .AND. l2 == 2 .AND. m == 0) THEN 01493 zzzz = & 01494 +add0/SQRT((r-da-db)**2+add)**3 + add0/SQRT((r+da+db)**2+add)**3 & 01495 +add0/SQRT((r-da+db)**2+add)**3 + add0/SQRT((r+da-db)**2+add)**3 & 01496 -2.0_dp*add0/SQRT((r-da)**2+db**2+add)**3 - 2.0_dp*add0/SQRT((r-db)**2+da**2+add)**3 & 01497 -2.0_dp*add0/SQRT((r+da)**2+db**2+add)**3 - 2.0_dp*add0/SQRT((r+db)**2+da**2+add)**3 & 01498 +2.0_dp*add0/SQRT(r**2+(da-db)**2+add)**3 + 2.0_dp*add0/SQRT(r**2+(da+db)**2+add)**3 01499 xyxy = & 01500 +4.0_dp*add0/SQRT(r**2+(da-db)**2+add)**3 + 4.0_dp*add0/SQRT(r**2+(da+db)**2+add)**3 & 01501 -8.0_dp*add0/SQRT(r**2+da**2+db**2+add)**3 01502 charg = -(zzzz*0.0625_dp - xyxy*0.015625_dp)*fact 01503 RETURN 01504 END IF 01505 ! X - ZX 01506 IF (l1 == 1 .AND. l2 == 2 .AND. m == 1) THEN 01507 ab = db/SQRT(2.0_dp) 01508 dxqxz = & 01509 -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 & 01510 +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 01511 charg = -dxqxz*0.125_dp*fact 01512 RETURN 01513 END IF 01514 ! ZX - ZX 01515 IF (l1 == 2 .AND. l2 == 2 .AND. m == 1) THEN 01516 aa = da/SQRT (2.0_dp) 01517 ab = db/SQRT (2.0_dp) 01518 qxzqxz = & 01519 +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 & 01520 -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 & 01521 -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 & 01522 +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 01523 charg = -qxzqxz*0.0625_dp*fact 01524 RETURN 01525 END IF 01526 ! XX - XX 01527 IF (l1 == 2 .AND. l2 == 2 .AND. m == 2) THEN 01528 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 -& 01529 8.0_dp*add0/SQRT(r**2+da**2+db**2+add)**3 01530 charg = -xyxy*0.0625_dp*fact 01531 RETURN 01532 END IF 01533 ! We should NEVER reach this point 01534 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 01535 END FUNCTION dcharg_int_ri_fs 01536 01537 ! ***************************************************************************** 01544 RECURSIVE SUBROUTINE rotmat (sepi, sepj, rjiv, r, ij_matrix, do_derivatives,& 01545 do_invert, debug_invert, error) 01546 TYPE(semi_empirical_type), POINTER :: sepi, sepj 01547 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rjiv 01548 REAL(KIND=dp), INTENT(IN) :: r 01549 TYPE(rotmat_type), POINTER :: ij_matrix 01550 LOGICAL, INTENT(IN) :: do_derivatives 01551 LOGICAL, INTENT(OUT), OPTIONAL :: do_invert 01552 LOGICAL, INTENT(IN), OPTIONAL :: debug_invert 01553 TYPE(cp_error_type), INTENT(inout) :: error 01554 01555 CHARACTER(len=*), PARAMETER :: routineN = 'rotmat', 01556 routineP = moduleN//':'//routineN 01557 01558 INTEGER :: imap(3), k, l 01559 LOGICAL :: dbg_inv, eval, failure 01560 REAL(KIND=dp) :: b, c2a, c2b, ca, ca2, cb, 01561 cb2, check, pt5sq3, r2i, s2a, 01562 s2b, sa, sb, sb2, sqb, sqb2i, 01563 x11, x22, x33 01564 REAL(KIND=dp), DIMENSION(3) :: b_d, c2a_d, c2b_d, ca2_d, ca_d, cb2_d, 01565 cb_d, r_d, s2a_d, s2b_d, sa_d, sb2_d, sb_d, sqb_d, x11_d, x22_d, x33_d 01566 REAL(KIND=dp), DIMENSION(3, 3) :: p 01567 REAL(KIND=dp), DIMENSION(3, 3, 3) :: p_d 01568 REAL(KIND=dp), DIMENSION(3, 5, 5) :: d_d 01569 REAL(KIND=dp), DIMENSION(5, 5) :: d 01570 01571 failure = .FALSE. 01572 CPPostcondition(ASSOCIATED(ij_matrix),cp_failure_level,routineP,error,failure) 01573 IF (PRESENT(do_invert)) do_invert = .FALSE. 01574 IF ((sepi%natorb>1).OR.(sepj%natorb>1)) THEN 01575 ! Compute Geomtric data and interatomic distance 01576 ! CA = COS(PHI) , SA = SIN(PHI) 01577 ! CB = COS(THETA) , SB = SIN(THETA) 01578 ! C2A = COS(2*PHI) , S2A = SIN(2*PHI) 01579 ! C2B = COS(2*THETA), S2B = SIN(2*PHI) 01580 eval = .TRUE. 01581 dbg_inv = .FALSE. 01582 pt5sq3 = 0.5_dp*SQRT(3.0_dp) 01583 imap(1) = 1 01584 imap(2) = 2 01585 imap(3) = 3 01586 check = rjiv(3)/r 01587 IF (PRESENT(debug_invert)) dbg_inv = debug_invert 01588 ! Check for special case: both atoms lie on the Z Axis 01589 IF (ABS(check) > 0.99999999_dp) THEN 01590 IF (PRESENT(do_invert)) THEN 01591 imap(1) = 3 01592 imap(2) = 2 01593 imap(3) = 1 01594 eval = .TRUE. 01595 do_invert = .TRUE. 01596 ELSE 01597 sa = 0.0_dp 01598 sb = 0.0_dp 01599 IF (check < 0.0_dp) THEN 01600 ca = -1.0_dp 01601 cb = -1.0_dp 01602 ELSE IF (check > 0.0_dp) THEN 01603 ca = 1.0_dp 01604 cb = 1.0_dp 01605 ELSE 01606 ca = 0.0_dp 01607 cb = 0.0_dp 01608 END IF 01609 eval = .FALSE. 01610 END IF 01611 END IF 01612 IF (dbg_inv) THEN 01613 ! When debugging the derivatives of the rotation matrix we must possibly force 01614 ! the inversion of the reference frame 01615 CPPostcondition(.NOT.do_derivatives,cp_failure_level,routineP,error,failure) 01616 imap(1) = 3 01617 imap(2) = 2 01618 imap(3) = 1 01619 eval = .TRUE. 01620 END IF 01621 IF (eval) THEN 01622 x11 = rjiv(imap(1)) 01623 x22 = rjiv(imap(2)) 01624 x33 = rjiv(imap(3)) 01625 cb = x33/r 01626 b = x11**2 + x22**2 01627 sqb = SQRT(b) 01628 ca = x11/sqb 01629 sa = x22/sqb 01630 sb = sqb/r 01631 END IF 01632 ! Calculate rotation matrix elements 01633 p(1, 1) = ca * sb 01634 p(2, 1) = ca * cb 01635 p(3, 1) = -sa 01636 p(1, 2) = sa * sb 01637 p(2, 2) = sa * cb 01638 p(3, 2) = ca 01639 p(1, 3) = cb 01640 p(2, 3) = -sb 01641 p(3, 3) = 0.0_dp 01642 IF (sepi%dorb.OR.sepj%dorb) THEN 01643 ca2 = ca**2 01644 cb2 = cb**2 01645 sb2 = sb**2 01646 c2a = 2.0_dp * ca2 - 1.0_dp 01647 c2b = 2.0_dp * cb2 - 1.0_dp 01648 s2a = 2.0_dp * sa * ca 01649 s2b = 2.0_dp * sb * cb 01650 d(1, 1) = pt5sq3 * c2a * sb2 01651 d(2, 1) = 0.5_dp * c2a * s2b 01652 d(3, 1) = -s2a * sb 01653 d(4, 1) = c2a * (cb2+0.5_dp*sb2) 01654 d(5, 1) = -s2a * cb 01655 d(1, 2) = pt5sq3 * ca * s2b 01656 d(2, 2) = ca * c2b 01657 d(3, 2) = -sa * cb 01658 d(4, 2) = -0.5_dp * ca * s2b 01659 d(5, 2) = sa * sb 01660 d(1, 3) = cb2 - 0.5_dp * sb2 01661 d(2, 3) = -pt5sq3 * s2b 01662 d(3, 3) = 0.0_dp 01663 d(4, 3) = pt5sq3 * sb2 01664 d(5, 3) = 0.0_dp 01665 d(1, 4) = pt5sq3 * sa * s2b 01666 d(2, 4) = sa * c2b 01667 d(3, 4) = ca * cb 01668 d(4, 4) = -0.5_dp * sa * s2b 01669 d(5, 4) = -ca * sb 01670 d(1, 5) = pt5sq3 * s2a * sb2 01671 d(2, 5) = 0.5_dp * s2a * s2b 01672 d(3, 5) = c2a * sb 01673 d(4, 5) = s2a * (cb2+0.5_dp*sb2) 01674 d(5, 5) = c2a * cb 01675 END IF 01676 ! Rotation Elements for : S-P 01677 DO k = 1, 3 01678 DO l = 1, 3 01679 ij_matrix%sp(k, l) = p(k, l) 01680 END DO 01681 END DO 01682 ! Rotation Elements for : P-P 01683 DO k = 1, 3 01684 ij_matrix%pp(1, k, k) = p(k, 1) * p(k, 1) 01685 ij_matrix%pp(2, k, k) = p(k, 2) * p(k, 2) 01686 ij_matrix%pp(3, k, k) = p(k, 3) * p(k, 3) 01687 ij_matrix%pp(4, k, k) = p(k, 1) * p(k, 2) 01688 ij_matrix%pp(5, k, k) = p(k, 1) * p(k, 3) 01689 ij_matrix%pp(6, k, k) = p(k, 2) * p(k, 3) 01690 IF (k /= 1) THEN 01691 DO l = 1, k - 1 01692 ij_matrix%pp(1, k, l) = 2.0_dp * p(k, 1) * p(l, 1) 01693 ij_matrix%pp(2, k, l) = 2.0_dp * p(k, 2) * p(l, 2) 01694 ij_matrix%pp(3, k, l) = 2.0_dp * p(k, 3) * p(l, 3) 01695 ij_matrix%pp(4, k, l) = p(k, 1) * p(l, 2) + p(k, 2) * p(l, 1) 01696 ij_matrix%pp(5, k, l) = p(k, 1) * p(l, 3) + p(k, 3) * p(l, 1) 01697 ij_matrix%pp(6, k, l) = p(k, 2) * p(l, 3) + p(k, 3) * p(l, 2) 01698 END DO 01699 END IF 01700 END DO 01701 IF (sepi%dorb.OR.sepj%dorb) THEN 01702 ! Rotation Elements for : S-D 01703 DO k = 1, 5 01704 DO l = 1, 5 01705 ij_matrix%sd(k, l) = d(k, l) 01706 END DO 01707 END DO 01708 ! Rotation Elements for : D-P 01709 DO k = 1, 5 01710 DO l = 1, 3 01711 ij_matrix%pd( 1, k, l) = d(k, 1) * p(l, 1) 01712 ij_matrix%pd( 2, k, l) = d(k, 1) * p(l, 2) 01713 ij_matrix%pd( 3, k, l) = d(k, 1) * p(l, 3) 01714 ij_matrix%pd( 4, k, l) = d(k, 2) * p(l, 1) 01715 ij_matrix%pd( 5, k, l) = d(k, 2) * p(l, 2) 01716 ij_matrix%pd( 6, k, l) = d(k, 2) * p(l, 3) 01717 ij_matrix%pd( 7, k, l) = d(k, 3) * p(l, 1) 01718 ij_matrix%pd( 8, k, l) = d(k, 3) * p(l, 2) 01719 ij_matrix%pd( 9, k, l) = d(k, 3) * p(l, 3) 01720 ij_matrix%pd(10, k, l) = d(k, 4) * p(l, 1) 01721 ij_matrix%pd(11, k, l) = d(k, 4) * p(l, 2) 01722 ij_matrix%pd(12, k, l) = d(k, 4) * p(l, 3) 01723 ij_matrix%pd(13, k, l) = d(k, 5) * p(l, 1) 01724 ij_matrix%pd(14, k, l) = d(k, 5) * p(l, 2) 01725 ij_matrix%pd(15, k, l) = d(k, 5) * p(l, 3) 01726 END DO 01727 END DO 01728 ! Rotation Elements for : D-D 01729 DO k = 1, 5 01730 ij_matrix%dd( 1, k, k) = d(k, 1) * d(k, 1) 01731 ij_matrix%dd( 2, k, k) = d(k, 2) * d(k, 2) 01732 ij_matrix%dd( 3, k, k) = d(k, 3) * d(k, 3) 01733 ij_matrix%dd( 4, k, k) = d(k, 4) * d(k, 4) 01734 ij_matrix%dd( 5, k, k) = d(k, 5) * d(k, 5) 01735 ij_matrix%dd( 6, k, k) = d(k, 1) * d(k, 2) 01736 ij_matrix%dd( 7, k, k) = d(k, 1) * d(k, 3) 01737 ij_matrix%dd( 8, k, k) = d(k, 2) * d(k, 3) 01738 ij_matrix%dd( 9, k, k) = d(k, 1) * d(k, 4) 01739 ij_matrix%dd(10, k, k) = d(k, 2) * d(k, 4) 01740 ij_matrix%dd(11, k, k) = d(k, 3) * d(k, 4) 01741 ij_matrix%dd(12, k, k) = d(k, 1) * d(k, 5) 01742 ij_matrix%dd(13, k, k) = d(k, 2) * d(k, 5) 01743 ij_matrix%dd(14, k, k) = d(k, 3) * d(k, 5) 01744 ij_matrix%dd(15, k, k) = d(k, 4) * d(k, 5) 01745 IF (k /= 1) THEN 01746 DO l = 1, k - 1 01747 ij_matrix%dd( 1, k, l) = 2.0_dp * d(k, 1) * d(l, 1) 01748 ij_matrix%dd( 2, k, l) = 2.0_dp * d(k, 2) * d(l, 2) 01749 ij_matrix%dd( 3, k, l) = 2.0_dp * d(k, 3) * d(l, 3) 01750 ij_matrix%dd( 4, k, l) = 2.0_dp * d(k, 4) * d(l, 4) 01751 ij_matrix%dd( 5, k, l) = 2.0_dp * d(k, 5) * d(l, 5) 01752 ij_matrix%dd( 6, k, l) = d(k, 1) * d(l, 2) + d(k, 2) * d(l, 1) 01753 ij_matrix%dd( 7, k, l) = d(k, 1) * d(l, 3) + d(k, 3) * d(l, 1) 01754 ij_matrix%dd( 8, k, l) = d(k, 2) * d(l, 3) + d(k, 3) * d(l, 2) 01755 ij_matrix%dd( 9, k, l) = d(k, 1) * d(l, 4) + d(k, 4) * d(l, 1) 01756 ij_matrix%dd(10, k, l) = d(k, 2) * d(l, 4) + d(k, 4) * d(l, 2) 01757 ij_matrix%dd(11, k, l) = d(k, 3) * d(l, 4) + d(k, 4) * d(l, 3) 01758 ij_matrix%dd(12, k, l) = d(k, 1) * d(l, 5) + d(k, 5) * d(l, 1) 01759 ij_matrix%dd(13, k, l) = d(k, 2) * d(l, 5) + d(k, 5) * d(l, 2) 01760 ij_matrix%dd(14, k, l) = d(k, 3) * d(l, 5) + d(k, 5) * d(l, 3) 01761 ij_matrix%dd(15, k, l) = d(k, 4) * d(l, 5) + d(k, 5) * d(l, 4) 01762 END DO 01763 END IF 01764 END DO 01765 END IF 01766 ! Evaluate Derivatives if required 01767 IF (do_derivatives) THEN 01768 ! This condition is necessary because derivatives uses the invertion of the 01769 ! axis for treating the divergence point along z-axis 01770 CPPostcondition(eval,cp_failure_level,routineP,error,failure) 01771 x11_d = 0.0_dp; x11_d(1) = 1.0_dp 01772 x22_d = 0.0_dp; x22_d(2) = 1.0_dp 01773 x33_d = 0.0_dp; x33_d(3) = 1.0_dp 01774 r_d = (/x11,x22,x33/)/r 01775 b_d = 2.0_dp*(x11*x11_d + x22*x22_d) 01776 sqb_d = (0.5_dp/sqb)*b_d 01777 r2i = 1.0_dp/r**2 01778 sqb2i = 1.0_dp/sqb**2 01779 cb_d = (x33_d*r-x33*r_d)*r2i 01780 ca_d = (x11_d*sqb-x11*sqb_d)*sqb2i 01781 sa_d = (x22_d*sqb-x22*sqb_d)*sqb2i 01782 sb_d = (sqb_d*r-sqb*r_d)*r2i 01783 ! Calculate derivatives of rotation matrix elements 01784 p_d(:, 1, 1) = ca_d * sb + ca * sb_d 01785 p_d(:, 2, 1) = ca_d * cb + ca * cb_d 01786 p_d(:, 3, 1) = -sa_d 01787 p_d(:, 1, 2) = sa_d * sb + sa * sb_d 01788 p_d(:, 2, 2) = sa_d * cb + sa * cb_d 01789 p_d(:, 3, 2) = ca_d 01790 p_d(:, 1, 3) = cb_d 01791 p_d(:, 2, 3) = -sb_d 01792 p_d(:, 3, 3) = 0.0_dp 01793 IF (sepi%dorb.OR.sepj%dorb) THEN 01794 ca2_d = 2.0_dp * ca * ca_d 01795 cb2_d = 2.0_dp * cb * cb_d 01796 sb2_d = 2.0_dp * sb * sb_d 01797 c2a_d = 2.0_dp * ca2_d 01798 c2b_d = 2.0_dp * cb2_d 01799 s2a_d = 2.0_dp * ( sa_d * ca + sa * ca_d ) 01800 s2b_d = 2.0_dp * ( sb_d * cb + sb * cb_d ) 01801 d_d(:, 1, 1) = pt5sq3 * ( c2a_d * sb2 + c2a * sb2_d ) 01802 d_d(:, 2, 1) = 0.5_dp * ( c2a_d * s2b + c2a * s2b_d ) 01803 d_d(:, 3, 1) = -s2a_d * sb - s2a * sb_d 01804 d_d(:, 4, 1) = c2a_d * (cb2+0.5_dp*sb2) + c2a * (cb2_d+0.5_dp*sb2_d) 01805 d_d(:, 5, 1) = -s2a_d * cb - s2a * cb_d 01806 d_d(:, 1, 2) = pt5sq3 * ( ca_d * s2b + ca * s2b_d ) 01807 d_d(:, 2, 2) = ca_d * c2b + ca * c2b_d 01808 d_d(:, 3, 2) = -sa_d * cb - sa * cb_d 01809 d_d(:, 4, 2) = -0.5_dp * ( ca_d * s2b + ca * s2b_d ) 01810 d_d(:, 5, 2) = sa_d * sb + sa * sb_d 01811 d_d(:, 1, 3) = cb2_d - 0.5_dp * sb2_d 01812 d_d(:, 2, 3) = -pt5sq3 * s2b_d 01813 d_d(:, 3, 3) = 0.0_dp 01814 d_d(:, 4, 3) = pt5sq3 * sb2_d 01815 d_d(:, 5, 3) = 0.0_dp 01816 d_d(:, 1, 4) = pt5sq3 * ( sa_d * s2b + sa * s2b_d ) 01817 d_d(:, 2, 4) = sa_d * c2b + sa * c2b_d 01818 d_d(:, 3, 4) = ca_d * cb + ca * cb_d 01819 d_d(:, 4, 4) = -0.5_dp * ( sa_d * s2b + sa * s2b_d ) 01820 d_d(:, 5, 4) = -ca_d * sb -ca * sb_d 01821 d_d(:, 1, 5) = pt5sq3 * ( s2a_d * sb2 + s2a * sb2_d ) 01822 d_d(:, 2, 5) = 0.5_dp * ( s2a_d * s2b + s2a * s2b_d ) 01823 d_d(:, 3, 5) = c2a_d * sb + c2a * sb_d 01824 d_d(:, 4, 5) = s2a_d * (cb2+0.5_dp*sb2) + s2a * (cb2_d+0.5_dp*sb2_d) 01825 d_d(:, 5, 5) = c2a_d * cb + c2a * cb_d 01826 END IF 01827 ! Derivative for Rotation Elements for : S-P 01828 DO k = 1, 3 01829 DO l = 1, 3 01830 ij_matrix%sp_d(:, k, l) = p_d(:, k, l) 01831 END DO 01832 END DO 01833 ! Derivative for Rotation Elements for : P-P 01834 DO k = 1, 3 01835 ij_matrix%pp_d(:, 1, k, k) = p_d(:, k, 1) * p(k, 1) + p(k, 1) * p_d(:, k, 1) 01836 ij_matrix%pp_d(:, 2, k, k) = p_d(:, k, 2) * p(k, 2) + p(k, 2) * p_d(:, k, 2) 01837 ij_matrix%pp_d(:, 3, k, k) = p_d(:, k, 3) * p(k, 3) + p(k, 3) * p_d(:, k, 3) 01838 ij_matrix%pp_d(:, 4, k, k) = p_d(:, k, 1) * p(k, 2) + p(k, 1) * p_d(:, k, 2) 01839 ij_matrix%pp_d(:, 5, k, k) = p_d(:, k, 1) * p(k, 3) + p(k, 1) * p_d(:, k, 3) 01840 ij_matrix%pp_d(:, 6, k, k) = p_d(:, k, 2) * p(k, 3) + p(k, 2) * p_d(:, k, 3) 01841 IF (k /= 1) THEN 01842 DO l = 1, k - 1 01843 ij_matrix%pp_d(:, 1, k, l) = 2.0_dp * ( p_d(:, k, 1) * p(l, 1) + p(k, 1) * p_d(:, l, 1) ) 01844 ij_matrix%pp_d(:, 2, k, l) = 2.0_dp * ( p_d(:, k, 2) * p(l, 2) + p(k, 2) * p_d(:, l, 2) ) 01845 ij_matrix%pp_d(:, 3, k, l) = 2.0_dp * ( p_d(:, k, 3) * p(l, 3) + p(k, 3) * p_d(:, l, 3) ) 01846 ij_matrix%pp_d(:, 4, k, l) = ( p_d(:, k, 1) * p(l, 2) + p(k, 1) * p_d(:, l, 2) ) + & 01847 ( p_d(:, k, 2) * p(l, 1) + p(k, 2) * p_d(:, l, 1) ) 01848 ij_matrix%pp_d(:, 5, k, l) = ( p_d(:, k, 1) * p(l, 3) + p(k, 1) * p_d(:, l, 3) ) + & 01849 ( p_d(:, k, 3) * p(l, 1) + p(k, 3) * p_d(:, l, 1) ) 01850 ij_matrix%pp_d(:, 6, k, l) = ( p_d(:, k, 2) * p(l, 3) + p(k, 2) * p_d(:, l, 3) ) + & 01851 ( p_d(:, k, 3) * p(l, 2) + p(k, 3) * p_d(:, l, 2) ) 01852 END DO 01853 END IF 01854 END DO 01855 IF (sepi%dorb.OR.sepj%dorb) THEN 01856 ! Rotation Elements for : S-D 01857 DO k = 1, 5 01858 DO l = 1, 5 01859 ij_matrix%sd_d(:, k, l) = d_d(:, k, l) 01860 END DO 01861 END DO 01862 ! Rotation Elements for : D-P 01863 DO k = 1, 5 01864 DO l = 1, 3 01865 ij_matrix%pd_d(:, 1, k, l) = ( d_d(:, k, 1) * p(l, 1) + d(k, 1) * p_d(:, l, 1) ) 01866 ij_matrix%pd_d(:, 2, k, l) = ( d_d(:, k, 1) * p(l, 2) + d(k, 1) * p_d(:, l, 2) ) 01867 ij_matrix%pd_d(:, 3, k, l) = ( d_d(:, k, 1) * p(l, 3) + d(k, 1) * p_d(:, l, 3) ) 01868 ij_matrix%pd_d(:, 4, k, l) = ( d_d(:, k, 2) * p(l, 1) + d(k, 2) * p_d(:, l, 1) ) 01869 ij_matrix%pd_d(:, 5, k, l) = ( d_d(:, k, 2) * p(l, 2) + d(k, 2) * p_d(:, l, 2) ) 01870 ij_matrix%pd_d(:, 6, k, l) = ( d_d(:, k, 2) * p(l, 3) + d(k, 2) * p_d(:, l, 3) ) 01871 ij_matrix%pd_d(:, 7, k, l) = ( d_d(:, k, 3) * p(l, 1) + d(k, 3) * p_d(:, l, 1) ) 01872 ij_matrix%pd_d(:, 8, k, l) = ( d_d(:, k, 3) * p(l, 2) + d(k, 3) * p_d(:, l, 2) ) 01873 ij_matrix%pd_d(:, 9, k, l) = ( d_d(:, k, 3) * p(l, 3) + d(k, 3) * p_d(:, l, 3) ) 01874 ij_matrix%pd_d(:, 10, k, l) = ( d_d(:, k, 4) * p(l, 1) + d(k, 4) * p_d(:, l, 1) ) 01875 ij_matrix%pd_d(:, 11, k, l) = ( d_d(:, k, 4) * p(l, 2) + d(k, 4) * p_d(:, l, 2) ) 01876 ij_matrix%pd_d(:, 12, k, l) = ( d_d(:, k, 4) * p(l, 3) + d(k, 4) * p_d(:, l, 3) ) 01877 ij_matrix%pd_d(:, 13, k, l) = ( d_d(:, k, 5) * p(l, 1) + d(k, 5) * p_d(:, l, 1) ) 01878 ij_matrix%pd_d(:, 14, k, l) = ( d_d(:, k, 5) * p(l, 2) + d(k, 5) * p_d(:, l, 2) ) 01879 ij_matrix%pd_d(:, 15, k, l) = ( d_d(:, k, 5) * p(l, 3) + d(k, 5) * p_d(:, l, 3) ) 01880 END DO 01881 END DO 01882 ! Rotation Elements for : D-D 01883 DO k = 1, 5 01884 ij_matrix%dd_d(:, 1, k, k) = ( d_d(:, k, 1) * d(k, 1) + d(k, 1) * d_d(:, k, 1) ) 01885 ij_matrix%dd_d(:, 2, k, k) = ( d_d(:, k, 2) * d(k, 2) + d(k, 2) * d_d(:, k, 2) ) 01886 ij_matrix%dd_d(:, 3, k, k) = ( d_d(:, k, 3) * d(k, 3) + d(k, 3) * d_d(:, k, 3) ) 01887 ij_matrix%dd_d(:, 4, k, k) = ( d_d(:, k, 4) * d(k, 4) + d(k, 4) * d_d(:, k, 4) ) 01888 ij_matrix%dd_d(:, 5, k, k) = ( d_d(:, k, 5) * d(k, 5) + d(k, 5) * d_d(:, k, 5) ) 01889 ij_matrix%dd_d(:, 6, k, k) = ( d_d(:, k, 1) * d(k, 2) + d(k, 1) * d_d(:, k, 2) ) 01890 ij_matrix%dd_d(:, 7, k, k) = ( d_d(:, k, 1) * d(k, 3) + d(k, 1) * d_d(:, k, 3) ) 01891 ij_matrix%dd_d(:, 8, k, k) = ( d_d(:, k, 2) * d(k, 3) + d(k, 2) * d_d(:, k, 3) ) 01892 ij_matrix%dd_d(:, 9, k, k) = ( d_d(:, k, 1) * d(k, 4) + d(k, 1) * d_d(:, k, 4) ) 01893 ij_matrix%dd_d(:, 10, k, k) = ( d_d(:, k, 2) * d(k, 4) + d(k, 2) * d_d(:, k, 4) ) 01894 ij_matrix%dd_d(:, 11, k, k) = ( d_d(:, k, 3) * d(k, 4) + d(k, 3) * d_d(:, k, 4) ) 01895 ij_matrix%dd_d(:, 12, k, k) = ( d_d(:, k, 1) * d(k, 5) + d(k, 1) * d_d(:, k, 5) ) 01896 ij_matrix%dd_d(:, 13, k, k) = ( d_d(:, k, 2) * d(k, 5) + d(k, 2) * d_d(:, k, 5) ) 01897 ij_matrix%dd_d(:, 14, k, k) = ( d_d(:, k, 3) * d(k, 5) + d(k, 3) * d_d(:, k, 5) ) 01898 ij_matrix%dd_d(:, 15, k, k) = ( d_d(:, k, 4) * d(k, 5) + d(k, 4) * d_d(:, k, 5) ) 01899 IF (k /= 1) THEN 01900 DO l = 1, k - 1 01901 ij_matrix%dd_d(:, 1, k, l) = 2.0_dp * ( d_d(:, k, 1) * d(l, 1) + d(k, 1) * d_d(:, l, 1) ) 01902 ij_matrix%dd_d(:, 2, k, l) = 2.0_dp * ( d_d(:, k, 2) * d(l, 2) + d(k, 2) * d_d(:, l, 2) ) 01903 ij_matrix%dd_d(:, 3, k, l) = 2.0_dp * ( d_d(:, k, 3) * d(l, 3) + d(k, 3) * d_d(:, l, 3) ) 01904 ij_matrix%dd_d(:, 4, k, l) = 2.0_dp * ( d_d(:, k, 4) * d(l, 4) + d(k, 4) * d_d(:, l, 4) ) 01905 ij_matrix%dd_d(:, 5, k, l) = 2.0_dp * ( d_d(:, k, 5) * d(l, 5) + d(k, 5) * d_d(:, l, 5) ) 01906 ij_matrix%dd_d(:, 6, k, l) = ( d_d(:, k, 1) * d(l, 2) + d(k, 1) * d_d(:, l, 2) ) + & 01907 ( d_d(:, k, 2) * d(l, 1) + d(k, 2) * d_d(:, l, 1) ) 01908 ij_matrix%dd_d(:, 7, k, l) = ( d_d(:, k, 1) * d(l, 3) + d(k, 1) * d_d(:, l, 3) ) + & 01909 ( d_d(:, k, 3) * d(l, 1) + d(k, 3) * d_d(:, l, 1) ) 01910 ij_matrix%dd_d(:, 8, k, l) = ( d_d(:, k, 2) * d(l, 3) + d(k, 2) * d_d(:, l, 3) ) + & 01911 ( d_d(:, k, 3) * d(l, 2) + d(k, 3) * d_d(:, l, 2) ) 01912 ij_matrix%dd_d(:, 9, k, l) = ( d_d(:, k, 1) * d(l, 4) + d(k, 1) * d_d(:, l, 4) ) + & 01913 ( d_d(:, k, 4) * d(l, 1) + d(k, 4) * d_d(:, l, 1) ) 01914 ij_matrix%dd_d(:, 10, k, l) = ( d_d(:, k, 2) * d(l, 4) + d(k, 2) * d_d(:, l, 4) ) + & 01915 ( d_d(:, k, 4) * d(l, 2) + d(k, 4) * d_d(:, l, 2) ) 01916 ij_matrix%dd_d(:, 11, k, l) = ( d_d(:, k, 3) * d(l, 4) + d(k, 3) * d_d(:, l, 4) ) + & 01917 ( d_d(:, k, 4) * d(l, 3) + d(k, 4) * d_d(:, l, 3) ) 01918 ij_matrix%dd_d(:, 12, k, l) = ( d_d(:, k, 1) * d(l, 5) + d(k, 1) * d_d(:, l, 5) ) + & 01919 ( d_d(:, k, 5) * d(l, 1) + d(k, 5) * d_d(:, l, 1) ) 01920 ij_matrix%dd_d(:, 13, k, l) = ( d_d(:, k, 2) * d(l, 5) + d(k, 2) * d_d(:, l, 5) ) + & 01921 ( d_d(:, k, 5) * d(l, 2) + d(k, 5) * d_d(:, l, 2) ) 01922 ij_matrix%dd_d(:, 14, k, l) = ( d_d(:, k, 3) * d(l, 5) + d(k, 3) * d_d(:, l, 5) ) + & 01923 ( d_d(:, k, 5) * d(l, 3) + d(k, 5) * d_d(:, l, 3) ) 01924 ij_matrix%dd_d(:, 15, k, l) = ( d_d(:, k, 4) * d(l, 5) + d(k, 4) * d_d(:, l, 5) ) + & 01925 ( d_d(:, k, 5) * d(l, 4) + d(k, 5) * d_d(:, l, 4) ) 01926 END DO 01927 END IF 01928 END DO 01929 END IF 01930 IF (debug_this_module) THEN 01931 CALL check_rotmat_der(sepi, sepj, rjiv, ij_matrix, do_invert=do_invert, error=error) 01932 END IF 01933 END IF 01934 END IF 01935 END SUBROUTINE rotmat 01936 01937 ! ***************************************************************************** 01944 RECURSIVE SUBROUTINE rot_2el_2c_first (sepi, sepj, rijv, se_int_control, se_taper,& 01945 invert, ii, kk, rep, logv, ij_matrix, v, lgrad, rep_d, v_d, logv_d, drij,& 01946 error) 01947 TYPE(semi_empirical_type), POINTER :: sepi, sepj 01948 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: rijv 01949 TYPE(se_int_control_type), INTENT(IN) :: se_int_control 01950 TYPE(se_taper_type), POINTER :: se_taper 01951 LOGICAL, INTENT(IN) :: invert 01952 INTEGER, INTENT(IN) :: ii, kk 01953 REAL(KIND=dp), DIMENSION(491), 01954 INTENT(IN) :: rep 01955 LOGICAL, DIMENSION(45, 45), INTENT(OUT) :: logv 01956 TYPE(rotmat_type), POINTER :: ij_matrix 01957 REAL(KIND=dp), DIMENSION(45, 45), 01958 INTENT(OUT) :: v 01959 LOGICAL, INTENT(IN) :: lgrad 01960 REAL(KIND=dp), DIMENSION(491), 01961 INTENT(IN), OPTIONAL :: rep_d 01962 REAL(KIND=dp), DIMENSION(3, 45, 45), 01963 INTENT(OUT), OPTIONAL :: v_d 01964 LOGICAL, DIMENSION(45, 45), 01965 INTENT(OUT), OPTIONAL :: logv_d 01966 REAL(KIND=dp), DIMENSION(3), 01967 INTENT(IN), OPTIONAL :: drij 01968 TYPE(cp_error_type), INTENT(inout) :: error 01969 01970 CHARACTER(len=*), PARAMETER :: routineN = 'rot_2el_2c_first', 01971 routineP = moduleN//':'//routineN 01972 01973 INTEGER :: i, i1, ij, j1, k, k1, kl, l, 01974 l1, limkl, ll, mm, nd 01975 LOGICAL :: failure 01976 REAL(KIND=dp) :: wrepp, wrepp_d(3) 01977 01978 failure = .FALSE. 01979 IF (lgrad) THEN 01980 CPPostcondition(PRESENT( rep_d),cp_failure_level,routineP,error,failure) 01981 CPPostcondition(PRESENT( v_d),cp_failure_level,routineP,error,failure) 01982 CPPostcondition(PRESENT(logv_d),cp_failure_level,routineP,error,failure) 01983 CPPostcondition(PRESENT( drij),cp_failure_level,routineP,error,failure) 01984 END IF 01985 limkl = indexb(kk, kk) 01986 DO k = 1, limkl 01987 DO i = 1, 45 01988 logv(i, k) = .FALSE. 01989 v(i, k) = 0.0_dp 01990 END DO 01991 END DO 01992 ! 01993 DO i1 = 1, ii 01994 DO j1 = 1, i1 01995 ij = indexa(i1, j1) 01996 ! 01997 DO k1 = 1, kk 01998 ! 01999 DO l1 = 1, k1 02000 kl = indexa(k1, l1) 02001 nd = ijkl_ind(ij, kl) 02002 IF (nd /= 0) THEN 02003 ! 02004 wrepp = rep(nd) 02005 ll = indexb(k1, l1) 02006 mm = int2c_type(ll) 02007 ! 02008 IF (mm == 1) THEN 02009 v(ij, 1) = wrepp 02010 ELSE IF (mm == 2) THEN 02011 k = k1 - 1 02012 v(ij, 2) = v(ij, 2) + ij_matrix%sp(k, 1) * wrepp 02013 v(ij, 4) = v(ij, 4) + ij_matrix%sp(k, 2) * wrepp 02014 v(ij, 7) = v(ij, 7) + ij_matrix%sp(k, 3) * wrepp 02015 ELSE IF (mm == 3) THEN 02016 k = k1 - 1 02017 l = l1 - 1 02018 v(ij, 3) = v(ij, 3) + ij_matrix%pp(1, k, l) * wrepp 02019 v(ij, 6) = v(ij, 6) + ij_matrix%pp(2, k, l) * wrepp 02020 v(ij, 10) = v(ij, 10) + ij_matrix%pp(3, k, l) * wrepp 02021 v(ij, 5) = v(ij, 5) + ij_matrix%pp(4, k, l) * wrepp 02022 v(ij, 8) = v(ij, 8) + ij_matrix%pp(5, k, l) * wrepp 02023 v(ij, 9) = v(ij, 9) + ij_matrix%pp(6, k, l) * wrepp 02024 ELSE IF (mm == 4) THEN 02025 k = k1 - 4 02026 v(ij, 11) = v(ij, 11) + ij_matrix%sd(k, 1) * wrepp 02027 v(ij, 16) = v(ij, 16) + ij_matrix%sd(k, 2) * wrepp 02028 v(ij, 22) = v(ij, 22) + ij_matrix%sd(k, 3) * wrepp 02029 v(ij, 29) = v(ij, 29) + ij_matrix%sd(k, 4) * wrepp 02030 v(ij, 37) = v(ij, 37) + ij_matrix%sd(k, 5) * wrepp 02031 ELSE IF (mm == 5) THEN 02032 k = k1 - 4 02033 l = l1 - 1 02034 v(ij, 12) = v(ij, 12) + ij_matrix%pd( 1, k, l) * wrepp 02035 v(ij, 13) = v(ij, 13) + ij_matrix%pd( 2, k, l) * wrepp 02036 v(ij, 14) = v(ij, 14) + ij_matrix%pd( 3, k, l) * wrepp 02037 v(ij, 17) = v(ij, 17) + ij_matrix%pd( 4, k, l) * wrepp 02038 v(ij, 18) = v(ij, 18) + ij_matrix%pd( 5, k, l) * wrepp 02039 v(ij, 19) = v(ij, 19) + ij_matrix%pd( 6, k, l) * wrepp 02040 v(ij, 23) = v(ij, 23) + ij_matrix%pd( 7, k, l) * wrepp 02041 v(ij, 24) = v(ij, 24) + ij_matrix%pd( 8, k, l) * wrepp 02042 v(ij, 25) = v(ij, 25) + ij_matrix%pd( 9, k, l) * wrepp 02043 v(ij, 30) = v(ij, 30) + ij_matrix%pd(10, k, l) * wrepp 02044 v(ij, 31) = v(ij, 31) + ij_matrix%pd(11, k, l) * wrepp 02045 v(ij, 32) = v(ij, 32) + ij_matrix%pd(12, k, l) * wrepp 02046 v(ij, 38) = v(ij, 38) + ij_matrix%pd(13, k, l) * wrepp 02047 v(ij, 39) = v(ij, 39) + ij_matrix%pd(14, k, l) * wrepp 02048 v(ij, 40) = v(ij, 40) + ij_matrix%pd(15, k, l) * wrepp 02049 ELSE IF (mm == 6) THEN 02050 k = k1 - 4 02051 l = l1 - 4 02052 v(ij, 15) = v(ij, 15) + ij_matrix%dd( 1, k, l) * wrepp 02053 v(ij, 21) = v(ij, 21) + ij_matrix%dd( 2, k, l) * wrepp 02054 v(ij, 28) = v(ij, 28) + ij_matrix%dd( 3, k, l) * wrepp 02055 v(ij, 36) = v(ij, 36) + ij_matrix%dd( 4, k, l) * wrepp 02056 v(ij, 45) = v(ij, 45) + ij_matrix%dd( 5, k, l) * wrepp 02057 v(ij, 20) = v(ij, 20) + ij_matrix%dd( 6, k, l) * wrepp 02058 v(ij, 26) = v(ij, 26) + ij_matrix%dd( 7, k, l) * wrepp 02059 v(ij, 27) = v(ij, 27) + ij_matrix%dd( 8, k, l) * wrepp 02060 v(ij, 33) = v(ij, 33) + ij_matrix%dd( 9, k, l) * wrepp 02061 v(ij, 34) = v(ij, 34) + ij_matrix%dd(10, k, l) * wrepp 02062 v(ij, 35) = v(ij, 35) + ij_matrix%dd(11, k, l) * wrepp 02063 v(ij, 41) = v(ij, 41) + ij_matrix%dd(12, k, l) * wrepp 02064 v(ij, 42) = v(ij, 42) + ij_matrix%dd(13, k, l) * wrepp 02065 v(ij, 43) = v(ij, 43) + ij_matrix%dd(14, k, l) * wrepp 02066 v(ij, 44) = v(ij, 44) + ij_matrix%dd(15, k, l) * wrepp 02067 END IF 02068 END IF 02069 END DO 02070 END DO 02071 DO kl = 1, limkl 02072 logv(ij, kl) = (ABS (v(ij, kl)) > 0.0_dp) 02073 END DO 02074 END DO 02075 END DO 02076 ! Gradients 02077 IF (lgrad) THEN 02078 DO k = 1, limkl 02079 DO i = 1, 45 02080 logv_d(i, k) = .FALSE. 02081 v_d(1, i, k) = 0.0_dp 02082 v_d(2, i, k) = 0.0_dp 02083 v_d(3, i, k) = 0.0_dp 02084 END DO 02085 END DO 02086 DO i1 = 1, ii 02087 DO j1 = 1, i1 02088 ij = indexa(i1, j1) 02089 ! 02090 DO k1 = 1, kk 02091 ! 02092 DO l1 = 1, k1 02093 kl = indexa(k1, l1) 02094 nd = ijkl_ind(ij, kl) 02095 IF (nd /= 0) THEN 02096 ! 02097 wrepp = rep(nd) 02098 wrepp_d = rep_d(nd)*drij 02099 ll = indexb(k1, l1) 02100 mm = int2c_type(ll) 02101 ! 02102 IF (mm == 1) THEN 02103 v_d(1, ij, 1) = wrepp_d(1) 02104 v_d(2, ij, 1) = wrepp_d(2) 02105 v_d(3, ij, 1) = wrepp_d(3) 02106 ELSE IF (mm == 2) THEN 02107 k = k1 - 1 02108 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) 02109 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) 02110 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) 02111 02112 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) 02113 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) 02114 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) 02115 02116 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) 02117 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) 02118 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) 02119 ELSE IF (mm == 3) THEN 02120 k = k1 - 1 02121 l = l1 - 1 02122 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) 02123 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) 02124 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) 02125 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) 02126 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) 02127 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) 02128 02129 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) 02130 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) 02131 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) 02132 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) 02133 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) 02134 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) 02135 02136 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) 02137 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) 02138 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) 02139 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) 02140 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) 02141 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) 02142 ELSE IF (mm == 4) THEN 02143 k = k1 - 4 02144 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) 02145 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) 02146 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) 02147 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) 02148 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) 02149 02150 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) 02151 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) 02152 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) 02153 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) 02154 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) 02155 02156 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) 02157 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) 02158 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) 02159 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) 02160 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) 02161 ELSE IF (mm == 5) THEN 02162 k = k1 - 4 02163 l = l1 - 1 02164 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) 02165 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) 02166 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) 02167 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) 02168 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) 02169 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) 02170 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) 02171 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) 02172 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) 02173 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) 02174 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) 02175 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) 02176 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) 02177 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) 02178 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) 02179 02180 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) 02181 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) 02182 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) 02183 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) 02184 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) 02185 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) 02186 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) 02187 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) 02188 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) 02189 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) 02190 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) 02191 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) 02192 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) 02193 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) 02194 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) 02195 02196 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) 02197 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) 02198 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) 02199 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) 02200 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) 02201 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) 02202 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) 02203 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) 02204 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) 02205 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) 02206 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) 02207 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) 02208 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) 02209 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) 02210 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) 02211 ELSE IF (mm == 6) THEN 02212 k = k1 - 4 02213 l = l1 - 4 02214 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) 02215 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) 02216 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) 02217 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) 02218 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) 02219 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) 02220 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) 02221 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) 02222 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) 02223 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) 02224 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) 02225 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) 02226 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) 02227 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) 02228 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) 02229 02230 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) 02231 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) 02232 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) 02233 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) 02234 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) 02235 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) 02236 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) 02237 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) 02238 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) 02239 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) 02240 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) 02241 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) 02242 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) 02243 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) 02244 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) 02245 02246 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) 02247 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) 02248 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) 02249 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) 02250 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) 02251 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) 02252 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) 02253 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) 02254 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) 02255 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) 02256 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) 02257 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) 02258 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) 02259 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) 02260 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) 02261 END IF 02262 END IF 02263 END DO 02264 END DO 02265 DO kl = 1, limkl 02266 logv_d(ij, kl) = (ANY(ABS(v_d(1:3, ij, kl)) > 0.0_dp)) 02267 END DO 02268 END DO 02269 END DO 02270 IF (debug_this_module) THEN 02271 CALL rot_2el_2c_first_debug(sepi, sepj, rijv, se_int_control, se_taper, invert, ii, kk, v_d, error) 02272 END IF 02273 END IF 02274 END SUBROUTINE rot_2el_2c_first 02275 02276 ! ***************************************************************************** 02282 SUBROUTINE store_2el_2c_diag (limij, limkl, ww, w, ww_dx, ww_dy, ww_dz, dw, error) 02283 02284 INTEGER, INTENT(IN) :: limij, limkl 02285 REAL(KIND=dp), DIMENSION(limkl, limij), 02286 INTENT(IN), OPTIONAL :: ww 02287 REAL(KIND=dp), DIMENSION(:), 02288 INTENT(INOUT), OPTIONAL :: w 02289 REAL(KIND=dp), DIMENSION(limkl, limij), 02290 INTENT(IN), OPTIONAL :: ww_dx, ww_dy, ww_dz 02291 REAL(KIND=dp), DIMENSION(:, :), 02292 INTENT(INOUT), OPTIONAL :: dw 02293 TYPE(cp_error_type), INTENT(inout) :: error 02294 02295 CHARACTER(len=*), PARAMETER :: routineN = 'store_2el_2c_diag', 02296 routineP = moduleN//':'//routineN 02297 02298 INTEGER :: ij, kl, l 02299 LOGICAL :: failure 02300 02301 failure = .FALSE. 02302 l = 0 02303 IF (PRESENT(ww).AND.PRESENT(w)) THEN 02304 DO ij = 1, limij 02305 DO kl = 1, limkl 02306 l = l + 1 02307 w(l) = ww(kl, ij) 02308 END DO 02309 END DO 02310 ELSE IF (PRESENT(ww_dx).AND.PRESENT(ww_dy).AND.PRESENT(ww_dz).AND.PRESENT(dw)) THEN 02311 DO ij = 1, limij 02312 DO kl = 1, limkl 02313 l = l + 1 02314 dw(1,l) = ww_dx(kl, ij) 02315 dw(2,l) = ww_dy(kl, ij) 02316 dw(3,l) = ww_dz(kl, ij) 02317 END DO 02318 END DO 02319 ELSE 02320 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 02321 END IF 02322 02323 END SUBROUTINE store_2el_2c_diag 02324 02325 END MODULE semi_empirical_int_utils
1.7.3