1#!-------------------------------------------------------------------------------------------------! 2#! CP2K: A general program to perform molecular dynamics simulations ! 3#! Copyright (C) 2000 - 2019 CP2K developers group ! 4#!-------------------------------------------------------------------------------------------------! 5#:mute 6 7#:def ewalds_multipole_sr_macro(mode, damping=False, store_energy=False, store_forces=False) 8 ! Compute the Short Range constribution according the task 9 IF (debug_this_module) THEN 10 f = HUGE(0.0_dp) 11 tij = HUGE(0.0_dp) 12 tij_a = HUGE(0.0_dp) 13 tij_ab = HUGE(0.0_dp) 14 tij_abc = HUGE(0.0_dp) 15 tij_abcd = HUGE(0.0_dp) 16 tij_abcde = HUGE(0.0_dp) 17 END IF 18 r = SQRT(rab2) 19 irab2 = 1.0_dp/rab2 20 ir = 1.0_dp/r 21 22 ! Compute the radial function 23#:if mode=="SCREENED_COULOMB_ERFC" 24 ! code for point multipole with screening 25 IF (debug_this_module.AND.debug_r_space.AND.(.NOT.debug_g_space)) THEN 26 f(0) = ir 27 tmp = 0.0_dp 28 ELSE 29 f(0) = erfc(alpha*r)*ir 30 tmp = EXP(-alpha**2*rab2)*oorootpi 31 END IF 32 fac = 1.0_dp 33 DO i = 1, 5 34 fac = fac*REAL(2*i-1,KIND=dp) 35 f(i) = irab2*(f(i-1)+ tmp*((2.0_dp*alpha**2)**i)/(fac*alpha)) 36 END DO 37 38#:elif mode=="SCREENED_COULOMB_GAUSS" 39 ! code for gaussian multipole with screening 40 IF (debug_this_module.AND.debug_r_space.AND.(.NOT.debug_g_space)) THEN 41 f(0) = ir 42 tmp1 = 0.0_dp 43 tmp2 = 0.0_dp 44 ELSE 45 f(0) = erf(beta*r)*ir - erf(alpha*r)*ir 46 tmp1 = EXP(-alpha**2*rab2)*oorootpi 47 tmp2 = EXP(-beta**2*rab2)*oorootpi 48 END IF 49 fac = 1.0_dp 50 DO i = 1, 5 51 fac = fac*REAL(2*i-1,KIND=dp) 52 f(i) = irab2*(f(i-1) + tmp1*((2.0_dp*alpha**2)**i)/(fac*alpha) - tmp2*((2.0_dp*beta**2)**i)/(fac*beta)) 53 END DO 54 55#:elif mode=="SCREENED_COULOMB_ERF" 56 IF (debug_this_module.AND.debug_r_space.AND.(.NOT.debug_g_space)) THEN 57 f(0) = ir 58 tmp = 0.0_dp 59 ELSE 60 f(0) = erf(alpha*r)*ir 61 tmp = EXP(-alpha**2*rab2)*oorootpi 62 END IF 63 fac = 1.0_dp 64 DO i = 1, 5 65 fac = fac*REAL(2*i-1,KIND=dp) 66 f(i) = irab2*(f(i-1) - tmp*((2.0_dp*alpha**2)**i)/(fac*alpha)) 67 END DO 68 69#:elif mode=="PURE_COULOMB" 70 ! code for point multipole without screening 71 f(0) = ir 72 DO i = 1, 5 73 f(i) = irab2*f(i-1) 74 END DO 75 76#:elif mode=="PURE_COULOMB_GAUSS" 77 ! code for gaussian multipole without screening 78 IF (debug_this_module.AND.debug_r_space.AND.(.NOT.debug_g_space)) THEN 79 f(0) = ir 80 tmp = 0.0_dp 81 ELSE 82 f(0) = erf(beta*r)*ir 83 tmp = EXP(-beta**2*rab2)*oorootpi 84 END IF 85 fac = 1.0_dp 86 PRINT *, "CHECK", f(0) 87 DO i = 1, 5 88 fac = fac*REAL(2*i-1,KIND=dp) 89 f(i) = irab2*(f(i-1) - tmp*((2.0_dp*beta**2)**i)/(fac*beta)) 90 END DO 91#:else 92#:stop "Unkown mode: "+mode 93#:endif 94 95 ! Compute the Tensor components 96 force_eval = do_stress 97 IF (task(1,1)) THEN 98 tij = f(0)*fac_ij 99 force_eval = do_forces .OR.do_efield1 100 END IF 101 IF (task(2,2)) force_eval = force_eval.OR.do_efield0 102 IF (task(1,2).OR.force_eval) THEN 103 force_eval = do_stress 104 tij_a = - rab*f(1)*fac_ij 105 IF (task(1,2)) force_eval = force_eval.OR. do_forces 106 END IF 107 IF (task(1,1)) force_eval = force_eval.OR.do_efield2 108 IF (task(3,3)) force_eval = force_eval.OR.do_efield0 109 IF (task(2,2).OR.task(3,1).OR.force_eval) THEN 110 force_eval = do_stress 111 DO b = 1,3 112 DO a = 1,3 113 tmp = rab(a)*rab(b)*fac_ij 114 tij_ab(a,b) = 3.0_dp*tmp*f(2) 115 IF (a==b) tij_ab(a,b) = tij_ab(a,b) - f(1)*fac_ij 116 END DO 117 END DO 118 IF (task(2,2).OR.task(3,1)) force_eval = force_eval.OR. do_forces 119 END IF 120 IF (task(2,2)) force_eval = force_eval.OR.do_efield2 121 IF (task(3,3)) force_eval = force_eval.OR.do_efield1 122 IF (task(3,2).OR.force_eval) THEN 123 force_eval = do_stress 124 DO c = 1, 3 125 DO b = 1, 3 126 DO a = 1, 3 127 tmp = rab(a)*rab(b)*rab(c)*fac_ij 128 tij_abc(a,b,c) = - 15.0_dp*tmp*f(3) 129 tmp = 3.0_dp*f(2)*fac_ij 130 IF (a==b) tij_abc(a,b,c) = tij_abc(a,b,c) + tmp*rab(c) 131 IF (a==c) tij_abc(a,b,c) = tij_abc(a,b,c) + tmp*rab(b) 132 IF (b==c) tij_abc(a,b,c) = tij_abc(a,b,c) + tmp*rab(a) 133 END DO 134 END DO 135 END DO 136 IF (task(3,2)) force_eval = force_eval.OR. do_forces 137 END IF 138 IF (task(3,3).OR.force_eval) THEN 139 force_eval = do_stress 140 DO d = 1, 3 141 DO c = 1, 3 142 DO b = 1, 3 143 DO a = 1, 3 144 tmp = rab(a)*rab(b)*rab(c)*rab(d)*fac_ij 145 tij_abcd(a,b,c,d) = 105.0_dp*tmp*f(4) 146 tmp1 = 15.0_dp*f(3)*fac_ij 147 tmp2 = 3.0_dp*f(2)*fac_ij 148 IF (a==b) THEN 149 tij_abcd(a,b,c,d) = tij_abcd(a,b,c,d) - tmp1*rab(c)*rab(d) 150 IF (c==d) tij_abcd(a,b,c,d) = tij_abcd(a,b,c,d) + tmp2 151 END IF 152 IF (a==c) THEN 153 tij_abcd(a,b,c,d) = tij_abcd(a,b,c,d) - tmp1*rab(b)*rab(d) 154 IF (b==d) tij_abcd(a,b,c,d) = tij_abcd(a,b,c,d) + tmp2 155 END IF 156 IF (a==d) tij_abcd(a,b,c,d) = tij_abcd(a,b,c,d) - tmp1*rab(b)*rab(c) 157 IF (b==c) THEN 158 tij_abcd(a,b,c,d) = tij_abcd(a,b,c,d) - tmp1*rab(a)*rab(d) 159 IF (a==d) tij_abcd(a,b,c,d) = tij_abcd(a,b,c,d) + tmp2 160 END IF 161 IF (b==d) tij_abcd(a,b,c,d) = tij_abcd(a,b,c,d) - tmp1*rab(a)*rab(c) 162 IF (c==d) tij_abcd(a,b,c,d) = tij_abcd(a,b,c,d) - tmp1*rab(a)*rab(b) 163 END DO 164 END DO 165 END DO 166 END DO 167 IF (task(3,3)) force_eval = force_eval.OR. do_forces 168 END IF 169 IF (force_eval) THEN 170 force_eval = do_stress 171 DO e = 1, 3 172 DO d = 1, 3 173 DO c = 1, 3 174 DO b = 1, 3 175 DO a = 1, 3 176 tmp = rab(a)*rab(b)*rab(c)*rab(d)*rab(e)*fac_ij 177 tij_abcde(a,b,c,d,e) = -945.0_dp*tmp*f(5) 178 tmp1 = 105.0_dp*f(4)*fac_ij 179 tmp2 = 15.0_dp*f(3)*fac_ij 180 IF (a==b) THEN 181 tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) + tmp1*rab(c)*rab(d)*rab(e) 182 IF (c==d) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(e) 183 IF (c==e) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(d) 184 IF (d==e) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(c) 185 END IF 186 IF (a==c) THEN 187 tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) + tmp1*rab(b)*rab(d)*rab(e) 188 IF (b==d) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(e) 189 IF (b==e) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(d) 190 IF (d==e) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(b) 191 END IF 192 IF (a==d) THEN 193 tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) + tmp1*rab(b)*rab(c)*rab(e) 194 IF (b==c) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(e) 195 IF (b==e) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(c) 196 IF (c==e) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(b) 197 END IF 198 IF (a==e) THEN 199 tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) + tmp1*rab(b)*rab(c)*rab(d) 200 IF (b==c) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(d) 201 IF (b==d) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(c) 202 IF (c==d) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(b) 203 END IF 204 IF (b==c) THEN 205 tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) + tmp1*rab(a)*rab(d)*rab(e) 206 IF (d==e) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(a) 207 END IF 208 IF (b==d) THEN 209 tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) + tmp1*rab(a)*rab(c)*rab(e) 210 IF (c==e) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(a) 211 END IF 212 IF (b==e) THEN 213 tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) + tmp1*rab(a)*rab(c)*rab(d) 214 IF (c==d) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) - tmp2*rab(a) 215 END IF 216 IF (c==d) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) + tmp1*rab(a)*rab(b)*rab(e) 217 IF (c==e) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) + tmp1*rab(a)*rab(b)*rab(d) 218 IF (d==e) tij_abcde(a,b,c,d,e) = tij_abcde(a,b,c,d,e) + tmp1*rab(a)*rab(b)*rab(c) 219 END DO 220 END DO 221 END DO 222 END DO 223 END DO 224 END IF 225 eloc = 0.0_dp 226 fr = 0.0_dp 227 ef0_i = 0.0_dp 228 ef0_j = 0.0_dp 229 ef1_j = 0.0_dp 230 ef1_i = 0.0_dp 231 ef2_j = 0.0_dp 232 ef2_i = 0.0_dp 233 234#:if damping 235 236 ! Initialize the damping function. 237 IF (kind_a==ikind) THEN 238 ! for atom i 239 SELECT CASE (itype_ij) 240 CASE (tang_toennies) 241 dampsumfi = 1.0_dp 242 xf = 1.0_dp 243 factorial = 1.0_dp 244 DO kk = 1, nkdamp_ij 245 xf = xf*dampa_ij*r 246 factorial = factorial * REAL(kk, KIND=dp) 247 dampsumfi = dampsumfi + (xf/factorial) 248 END DO 249 dampaexpi = dexp(-dampa_ij * r) 250 dampfunci = dampsumfi * dampaexpi * dampfac_ij 251 dampfuncdiffi = -dampa_ij * dampaexpi * & 252 dampfac_ij * (((dampa_ij * r) ** nkdamp_ij) / & 253 factorial) 254 CASE DEFAULT 255 dampfunci=0.0_dp 256 dampfuncdiffi=0.0_dp 257 END SELECT 258 259 ! for atom j 260 SELECT CASE (itype_ji) 261 CASE (tang_toennies) 262 dampsumfj = 1.0_dp 263 xf = 1.0_dp 264 factorial = 1.0_dp 265 DO kk = 1, nkdamp_ji 266 xf = xf*dampa_ji*r 267 factorial = factorial * REAL(kk, KIND=dp) 268 dampsumfj = dampsumfj + (xf/factorial) 269 END DO 270 dampaexpj = dexp(-dampa_ji * r) 271 dampfuncj = dampsumfj * dampaexpj * dampfac_ji 272 dampfuncdiffj = -dampa_ji * dampaexpj * & 273 dampfac_ji * (((dampa_ji * r) ** nkdamp_ji) / & 274 factorial) 275 CASE DEFAULT 276 dampfuncj = 0.0_dp 277 dampfuncdiffj = 0.0_dp 278 END SELECT 279 ELSE 280 SELECT CASE (itype_ij) 281 CASE(tang_toennies) 282 dampsumfj = 1.0_dp 283 xf = 1.0_dp 284 factorial = 1.0_dp 285 DO kk = 1, nkdamp_ij 286 xf = xf*dampa_ij*r 287 factorial = factorial * REAL(kk, KIND=dp) 288 dampsumfj = dampsumfj + (xf/factorial) 289 END DO 290 dampaexpj = dexp(-dampa_ij * r) 291 dampfuncj = dampsumfj * dampaexpj * dampfac_ij 292 dampfuncdiffj = -dampa_ij * dampaexpj * & 293 dampfac_ij * (((dampa_ij * r) ** nkdamp_ij) / & 294 factorial) 295 CASE DEFAULT 296 dampfuncj=0.0_dp 297 dampfuncdiffj=0.0_dp 298 END SELECT 299 300 !for j 301 SELECT CASE (itype_ji) 302 CASE (tang_toennies) 303 dampsumfi = 1.0_dp 304 xf = 1.0_dp 305 factorial = 1.0_dp 306 DO kk = 1, nkdamp_ji 307 xf = xf*dampa_ji*r 308 factorial = factorial * REAL(kk, KIND=dp) 309 dampsumfi = dampsumfi + (xf/factorial) 310 END DO 311 dampaexpi = dexp(-dampa_ji * r) 312 dampfunci = dampsumfi * dampaexpi * dampfac_ji 313 dampfuncdiffi = -dampa_ji * dampaexpi * & 314 dampfac_ji * (((dampa_ji * r) ** nkdamp_ji) / & 315 factorial) 316 CASE DEFAULT 317 dampfunci = 0.0_dp 318 dampfuncdiffi = 0.0_dp 319 END SELECT 320 END IF 321 322 damptij_a = -rab*dampfunci*fac_ij*irab2*ir 323 damptji_a = -rab*dampfuncj*fac_ij*irab2*ir 324 DO b = 1,3 325 DO a = 1,3 326 tmp = rab(a)*rab(b)*fac_ij 327 damptij_ab(a,b) = tmp*(-dampfuncdiffi*irab2*irab2+3.0_dp*dampfunci*irab2*irab2*ir) 328 damptji_ab(a,b) = tmp*(-dampfuncdiffj*irab2*irab2+3.0_dp*dampfuncj*irab2*irab2*ir) 329 IF (a==b) damptij_ab(a,b) = damptij_ab(a,b) - dampfunci*fac_ij*irab2*ir 330 IF (a==b) damptji_ab(a,b) = damptji_ab(a,b) - dampfuncj*fac_ij*irab2*ir 331 END DO 332 END DO 333 334#:endif 335 336 ! Initialize the charge, dipole and quadrupole for atom A and B 337 IF (debug_this_module) THEN 338 ch_j = HUGE(0.0_dp) 339 ch_i = HUGE(0.0_dp) 340 dp_j = HUGE(0.0_dp) 341 dp_i = HUGE(0.0_dp) 342 qp_j = HUGE(0.0_dp) 343 qp_i = HUGE(0.0_dp) 344 END IF 345 IF (ANY(task(1,:))) THEN 346 ch_j = charges(atom_a) 347 ch_i = charges(atom_b) 348 END IF 349 IF (ANY(task(2,:))) THEN 350 dp_j = dipoles(:,atom_a) 351 dp_i = dipoles(:,atom_b) 352 END IF 353 IF (ANY(task(3,:))) THEN 354 qp_j = quadrupoles(:,:,atom_a) 355 qp_i = quadrupoles(:,:,atom_b) 356 END IF 357 IF (task(1,1)) THEN 358 ! Charge - Charge 359 eloc = eloc + ch_i*tij*ch_j 360 ! Forces on particle i (locally b) 361 IF (do_forces.OR.do_stress) THEN 362 fr(1) = fr(1) - ch_j * tij_a(1) * ch_i 363 fr(2) = fr(2) - ch_j * tij_a(2) * ch_i 364 fr(3) = fr(3) - ch_j * tij_a(3) * ch_i 365 END IF 366 ! Electric fields 367 IF (do_efield) THEN 368 ! Potential 369 IF (do_efield0) THEN 370 ef0_i = ef0_i + tij * ch_j 371 372 ef0_j = ef0_j + tij * ch_i 373 END IF 374 ! Electric field 375 IF (do_efield1) THEN 376 ef1_i(1) = ef1_i(1) - tij_a(1) * ch_j 377 ef1_i(2) = ef1_i(2) - tij_a(2) * ch_j 378 ef1_i(3) = ef1_i(3) - tij_a(3) * ch_j 379 380 ef1_j(1) = ef1_j(1) + tij_a(1) * ch_i 381 ef1_j(2) = ef1_j(2) + tij_a(2) * ch_i 382 ef1_j(3) = ef1_j(3) + tij_a(3) * ch_i 383 384#:if damping 385 ef1_i(1) = ef1_i(1) + damptij_a(1) * ch_j 386 ef1_i(2) = ef1_i(2) + damptij_a(2) * ch_j 387 ef1_i(3) = ef1_i(3) + damptij_a(3) * ch_j 388 389 ef1_j(1) = ef1_j(1) - damptji_a(1) * ch_i 390 ef1_j(2) = ef1_j(2) - damptji_a(2) * ch_i 391 ef1_j(3) = ef1_j(3) - damptji_a(3) * ch_i 392#:endif 393 394 END IF 395 ! Electric field gradient 396 IF (do_efield2) THEN 397 ef2_i(1,1) = ef2_i(1,1) - tij_ab(1,1) * ch_j 398 ef2_i(2,1) = ef2_i(2,1) - tij_ab(2,1) * ch_j 399 ef2_i(3,1) = ef2_i(3,1) - tij_ab(3,1) * ch_j 400 ef2_i(1,2) = ef2_i(1,2) - tij_ab(1,2) * ch_j 401 ef2_i(2,2) = ef2_i(2,2) - tij_ab(2,2) * ch_j 402 ef2_i(3,2) = ef2_i(3,2) - tij_ab(3,2) * ch_j 403 ef2_i(1,3) = ef2_i(1,3) - tij_ab(1,3) * ch_j 404 ef2_i(2,3) = ef2_i(2,3) - tij_ab(2,3) * ch_j 405 ef2_i(3,3) = ef2_i(3,3) - tij_ab(3,3) * ch_j 406 407 ef2_j(1,1) = ef2_j(1,1) - tij_ab(1,1) * ch_i 408 ef2_j(2,1) = ef2_j(2,1) - tij_ab(2,1) * ch_i 409 ef2_j(3,1) = ef2_j(3,1) - tij_ab(3,1) * ch_i 410 ef2_j(1,2) = ef2_j(1,2) - tij_ab(1,2) * ch_i 411 ef2_j(2,2) = ef2_j(2,2) - tij_ab(2,2) * ch_i 412 ef2_j(3,2) = ef2_j(3,2) - tij_ab(3,2) * ch_i 413 ef2_j(1,3) = ef2_j(1,3) - tij_ab(1,3) * ch_i 414 ef2_j(2,3) = ef2_j(2,3) - tij_ab(2,3) * ch_i 415 ef2_j(3,3) = ef2_j(3,3) - tij_ab(3,3) * ch_i 416 END IF 417 END IF 418 END IF 419 IF (task(2,2)) THEN 420 ! Dipole - Dipole 421 tmp= - (dp_i(1)*(tij_ab(1,1)*dp_j(1)+& 422 tij_ab(2,1)*dp_j(2)+& 423 tij_ab(3,1)*dp_j(3))+& 424 dp_i(2)*(tij_ab(1,2)*dp_j(1)+& 425 tij_ab(2,2)*dp_j(2)+& 426 tij_ab(3,2)*dp_j(3))+& 427 dp_i(3)*(tij_ab(1,3)*dp_j(1)+& 428 tij_ab(2,3)*dp_j(2)+& 429 tij_ab(3,3)*dp_j(3))) 430 eloc = eloc + tmp 431 ! Forces on particle i (locally b) 432 IF (do_forces.OR.do_stress) THEN 433 DO k = 1, 3 434 fr(k) = fr(k) + dp_i(1)*(tij_abc(1,1,k)*dp_j(1)+& 435 tij_abc(2,1,k)*dp_j(2)+& 436 tij_abc(3,1,k)*dp_j(3))& 437 + dp_i(2)*(tij_abc(1,2,k)*dp_j(1)+& 438 tij_abc(2,2,k)*dp_j(2)+& 439 tij_abc(3,2,k)*dp_j(3))& 440 + dp_i(3)*(tij_abc(1,3,k)*dp_j(1)+& 441 tij_abc(2,3,k)*dp_j(2)+& 442 tij_abc(3,3,k)*dp_j(3)) 443 END DO 444 END IF 445 ! Electric fields 446 IF (do_efield) THEN 447 ! Potential 448 IF (do_efield0) THEN 449 ef0_i = ef0_i - (tij_a(1)*dp_j(1)+& 450 tij_a(2)*dp_j(2)+& 451 tij_a(3)*dp_j(3)) 452 453 ef0_j = ef0_j + (tij_a(1)*dp_i(1)+& 454 tij_a(2)*dp_i(2)+& 455 tij_a(3)*dp_i(3)) 456 END IF 457 ! Electric field 458 IF (do_efield1) THEN 459 ef1_i(1) = ef1_i(1) + (tij_ab(1,1)*dp_j(1)+& 460 tij_ab(2,1)*dp_j(2)+& 461 tij_ab(3,1)*dp_j(3)) 462 ef1_i(2) = ef1_i(2) + (tij_ab(1,2)*dp_j(1)+& 463 tij_ab(2,2)*dp_j(2)+& 464 tij_ab(3,2)*dp_j(3)) 465 ef1_i(3) = ef1_i(3) + (tij_ab(1,3)*dp_j(1)+& 466 tij_ab(2,3)*dp_j(2)+& 467 tij_ab(3,3)*dp_j(3)) 468 469 ef1_j(1) = ef1_j(1) + (tij_ab(1,1)*dp_i(1)+& 470 tij_ab(2,1)*dp_i(2)+& 471 tij_ab(3,1)*dp_i(3)) 472 ef1_j(2) = ef1_j(2) + (tij_ab(1,2)*dp_i(1)+& 473 tij_ab(2,2)*dp_i(2)+& 474 tij_ab(3,2)*dp_i(3)) 475 ef1_j(3) = ef1_j(3) + (tij_ab(1,3)*dp_i(1)+& 476 tij_ab(2,3)*dp_i(2)+& 477 tij_ab(3,3)*dp_i(3)) 478 END IF 479 ! Electric field gradient 480 IF (do_efield2) THEN 481 ef2_i(1,1) = ef2_i(1,1) + (tij_abc(1,1,1)*dp_j(1)+& 482 tij_abc(2,1,1)*dp_j(2)+& 483 tij_abc(3,1,1)*dp_j(3)) 484 ef2_i(1,2) = ef2_i(1,2) + (tij_abc(1,1,2)*dp_j(1)+& 485 tij_abc(2,1,2)*dp_j(2)+& 486 tij_abc(3,1,2)*dp_j(3)) 487 ef2_i(1,3) = ef2_i(1,3) + (tij_abc(1,1,3)*dp_j(1)+& 488 tij_abc(2,1,3)*dp_j(2)+& 489 tij_abc(3,1,3)*dp_j(3)) 490 ef2_i(2,1) = ef2_i(2,1) + (tij_abc(1,2,1)*dp_j(1)+& 491 tij_abc(2,2,1)*dp_j(2)+& 492 tij_abc(3,2,1)*dp_j(3)) 493 ef2_i(2,2) = ef2_i(2,2) + (tij_abc(1,2,2)*dp_j(1)+& 494 tij_abc(2,2,2)*dp_j(2)+& 495 tij_abc(3,2,2)*dp_j(3)) 496 ef2_i(2,3) = ef2_i(2,3) + (tij_abc(1,2,3)*dp_j(1)+& 497 tij_abc(2,2,3)*dp_j(2)+& 498 tij_abc(3,2,3)*dp_j(3)) 499 ef2_i(3,1) = ef2_i(3,1) + (tij_abc(1,3,1)*dp_j(1)+& 500 tij_abc(2,3,1)*dp_j(2)+& 501 tij_abc(3,3,1)*dp_j(3)) 502 ef2_i(3,2) = ef2_i(3,2) + (tij_abc(1,3,2)*dp_j(1)+& 503 tij_abc(2,3,2)*dp_j(2)+& 504 tij_abc(3,3,2)*dp_j(3)) 505 ef2_i(3,3) = ef2_i(3,3) + (tij_abc(1,3,3)*dp_j(1)+& 506 tij_abc(2,3,3)*dp_j(2)+& 507 tij_abc(3,3,3)*dp_j(3)) 508 509 ef2_j(1,1) = ef2_j(1,1) - (tij_abc(1,1,1)*dp_i(1)+& 510 tij_abc(2,1,1)*dp_i(2)+& 511 tij_abc(3,1,1)*dp_i(3)) 512 ef2_j(1,2) = ef2_j(1,2) - (tij_abc(1,1,2)*dp_i(1)+& 513 tij_abc(2,1,2)*dp_i(2)+& 514 tij_abc(3,1,2)*dp_i(3)) 515 ef2_j(1,3) = ef2_j(1,3) - (tij_abc(1,1,3)*dp_i(1)+& 516 tij_abc(2,1,3)*dp_i(2)+& 517 tij_abc(3,1,3)*dp_i(3)) 518 ef2_j(2,1) = ef2_j(2,1) - (tij_abc(1,2,1)*dp_i(1)+& 519 tij_abc(2,2,1)*dp_i(2)+& 520 tij_abc(3,2,1)*dp_i(3)) 521 ef2_j(2,2) = ef2_j(2,2) - (tij_abc(1,2,2)*dp_i(1)+& 522 tij_abc(2,2,2)*dp_i(2)+& 523 tij_abc(3,2,2)*dp_i(3)) 524 ef2_j(2,3) = ef2_j(2,3) - (tij_abc(1,2,3)*dp_i(1)+& 525 tij_abc(2,2,3)*dp_i(2)+& 526 tij_abc(3,2,3)*dp_i(3)) 527 ef2_j(3,1) = ef2_j(3,1) - (tij_abc(1,3,1)*dp_i(1)+& 528 tij_abc(2,3,1)*dp_i(2)+& 529 tij_abc(3,3,1)*dp_i(3)) 530 ef2_j(3,2) = ef2_j(3,2) - (tij_abc(1,3,2)*dp_i(1)+& 531 tij_abc(2,3,2)*dp_i(2)+& 532 tij_abc(3,3,2)*dp_i(3)) 533 ef2_j(3,3) = ef2_j(3,3) - (tij_abc(1,3,3)*dp_i(1)+& 534 tij_abc(2,3,3)*dp_i(2)+& 535 tij_abc(3,3,3)*dp_i(3)) 536 END IF 537 END IF 538 END IF 539 IF (task(2,1)) THEN 540 ! Dipole - Charge 541 tmp= ch_j*(tij_a(1)*dp_i(1)+& 542 tij_a(2)*dp_i(2)+& 543 tij_a(3)*dp_i(3))& 544 - ch_i*(tij_a(1)*dp_j(1)+& 545 tij_a(2)*dp_j(2)+& 546 tij_a(3)*dp_j(3)) 547#:if damping 548 tmp= tmp- ch_j*(damptij_a(1)*dp_i(1)+& 549 damptij_a(2)*dp_i(2)+& 550 damptij_a(3)*dp_i(3))& 551 + ch_i*(damptji_a(1)*dp_j(1)+& 552 damptji_a(2)*dp_j(2)+& 553 damptji_a(3)*dp_j(3)) 554#:endif 555 eloc = eloc + tmp 556 ! Forces on particle i (locally b) 557 IF (do_forces.OR.do_stress) THEN 558 DO k = 1, 3 559 fr(k) = fr(k) - ch_j *(tij_ab(1,k)*dp_i(1)+& 560 tij_ab(2,k)*dp_i(2)+& 561 tij_ab(3,k)*dp_i(3))& 562 + ch_i *(tij_ab(1,k)*dp_j(1)+& 563 tij_ab(2,k)*dp_j(2)+& 564 tij_ab(3,k)*dp_j(3)) 565#:if damping 566 fr(k) = fr(k) + ch_j *(damptij_ab(1,k)*dp_i(1)+& 567 damptij_ab(2,k)*dp_i(2)+& 568 damptij_ab(3,k)*dp_i(3))& 569 - ch_i *(damptji_ab(1,k)*dp_j(1)+& 570 damptji_ab(2,k)*dp_j(2)+& 571 damptji_ab(3,k)*dp_j(3)) 572#:endif 573 END DO 574 END IF 575 END IF 576 IF (task(3,3)) THEN 577 ! Quadrupole - Quadrupole 578 fac = 1.0_dp/9.0_dp 579 tmp11 = qp_i(1,1)*(tij_abcd(1,1,1,1)*qp_j(1,1)+& 580 tij_abcd(2,1,1,1)*qp_j(2,1)+& 581 tij_abcd(3,1,1,1)*qp_j(3,1)+& 582 tij_abcd(1,2,1,1)*qp_j(1,2)+& 583 tij_abcd(2,2,1,1)*qp_j(2,2)+& 584 tij_abcd(3,2,1,1)*qp_j(3,2)+& 585 tij_abcd(1,3,1,1)*qp_j(1,3)+& 586 tij_abcd(2,3,1,1)*qp_j(2,3)+& 587 tij_abcd(3,3,1,1)*qp_j(3,3)) 588 tmp21 = qp_i(2,1)*(tij_abcd(1,1,1,2)*qp_j(1,1)+& 589 tij_abcd(2,1,1,2)*qp_j(2,1)+& 590 tij_abcd(3,1,1,2)*qp_j(3,1)+& 591 tij_abcd(1,2,1,2)*qp_j(1,2)+& 592 tij_abcd(2,2,1,2)*qp_j(2,2)+& 593 tij_abcd(3,2,1,2)*qp_j(3,2)+& 594 tij_abcd(1,3,1,2)*qp_j(1,3)+& 595 tij_abcd(2,3,1,2)*qp_j(2,3)+& 596 tij_abcd(3,3,1,2)*qp_j(3,3)) 597 tmp31 = qp_i(3,1)*(tij_abcd(1,1,1,3)*qp_j(1,1)+& 598 tij_abcd(2,1,1,3)*qp_j(2,1)+& 599 tij_abcd(3,1,1,3)*qp_j(3,1)+& 600 tij_abcd(1,2,1,3)*qp_j(1,2)+& 601 tij_abcd(2,2,1,3)*qp_j(2,2)+& 602 tij_abcd(3,2,1,3)*qp_j(3,2)+& 603 tij_abcd(1,3,1,3)*qp_j(1,3)+& 604 tij_abcd(2,3,1,3)*qp_j(2,3)+& 605 tij_abcd(3,3,1,3)*qp_j(3,3)) 606 tmp22 = qp_i(2,2)*(tij_abcd(1,1,2,2)*qp_j(1,1)+& 607 tij_abcd(2,1,2,2)*qp_j(2,1)+& 608 tij_abcd(3,1,2,2)*qp_j(3,1)+& 609 tij_abcd(1,2,2,2)*qp_j(1,2)+& 610 tij_abcd(2,2,2,2)*qp_j(2,2)+& 611 tij_abcd(3,2,2,2)*qp_j(3,2)+& 612 tij_abcd(1,3,2,2)*qp_j(1,3)+& 613 tij_abcd(2,3,2,2)*qp_j(2,3)+& 614 tij_abcd(3,3,2,2)*qp_j(3,3)) 615 tmp32 = qp_i(3,2)*(tij_abcd(1,1,2,3)*qp_j(1,1)+& 616 tij_abcd(2,1,2,3)*qp_j(2,1)+& 617 tij_abcd(3,1,2,3)*qp_j(3,1)+& 618 tij_abcd(1,2,2,3)*qp_j(1,2)+& 619 tij_abcd(2,2,2,3)*qp_j(2,2)+& 620 tij_abcd(3,2,2,3)*qp_j(3,2)+& 621 tij_abcd(1,3,2,3)*qp_j(1,3)+& 622 tij_abcd(2,3,2,3)*qp_j(2,3)+& 623 tij_abcd(3,3,2,3)*qp_j(3,3)) 624 tmp33 = qp_i(3,3)*(tij_abcd(1,1,3,3)*qp_j(1,1)+& 625 tij_abcd(2,1,3,3)*qp_j(2,1)+& 626 tij_abcd(3,1,3,3)*qp_j(3,1)+& 627 tij_abcd(1,2,3,3)*qp_j(1,2)+& 628 tij_abcd(2,2,3,3)*qp_j(2,2)+& 629 tij_abcd(3,2,3,3)*qp_j(3,2)+& 630 tij_abcd(1,3,3,3)*qp_j(1,3)+& 631 tij_abcd(2,3,3,3)*qp_j(2,3)+& 632 tij_abcd(3,3,3,3)*qp_j(3,3)) 633 tmp12 = tmp21 634 tmp13 = tmp31 635 tmp23 = tmp32 636 tmp = tmp11 + tmp12 + tmp13 + & 637 tmp21 + tmp22 + tmp23 + & 638 tmp31 + tmp32 + tmp33 639 640 eloc = eloc + fac*tmp 641 ! Forces on particle i (locally b) 642 IF (do_forces.OR.do_stress) THEN 643 DO k = 1, 3 644 tmp11 = qp_i(1,1)*(tij_abcde(1,1,1,1,k)*qp_j(1,1)+& 645 tij_abcde(2,1,1,1,k)*qp_j(2,1)+& 646 tij_abcde(3,1,1,1,k)*qp_j(3,1)+& 647 tij_abcde(1,2,1,1,k)*qp_j(1,2)+& 648 tij_abcde(2,2,1,1,k)*qp_j(2,2)+& 649 tij_abcde(3,2,1,1,k)*qp_j(3,2)+& 650 tij_abcde(1,3,1,1,k)*qp_j(1,3)+& 651 tij_abcde(2,3,1,1,k)*qp_j(2,3)+& 652 tij_abcde(3,3,1,1,k)*qp_j(3,3)) 653 tmp21 = qp_i(2,1)*(tij_abcde(1,1,2,1,k)*qp_j(1,1)+& 654 tij_abcde(2,1,2,1,k)*qp_j(2,1)+& 655 tij_abcde(3,1,2,1,k)*qp_j(3,1)+& 656 tij_abcde(1,2,2,1,k)*qp_j(1,2)+& 657 tij_abcde(2,2,2,1,k)*qp_j(2,2)+& 658 tij_abcde(3,2,2,1,k)*qp_j(3,2)+& 659 tij_abcde(1,3,2,1,k)*qp_j(1,3)+& 660 tij_abcde(2,3,2,1,k)*qp_j(2,3)+& 661 tij_abcde(3,3,2,1,k)*qp_j(3,3)) 662 tmp31 = qp_i(3,1)*(tij_abcde(1,1,3,1,k)*qp_j(1,1)+& 663 tij_abcde(2,1,3,1,k)*qp_j(2,1)+& 664 tij_abcde(3,1,3,1,k)*qp_j(3,1)+& 665 tij_abcde(1,2,3,1,k)*qp_j(1,2)+& 666 tij_abcde(2,2,3,1,k)*qp_j(2,2)+& 667 tij_abcde(3,2,3,1,k)*qp_j(3,2)+& 668 tij_abcde(1,3,3,1,k)*qp_j(1,3)+& 669 tij_abcde(2,3,3,1,k)*qp_j(2,3)+& 670 tij_abcde(3,3,3,1,k)*qp_j(3,3)) 671 tmp22 = qp_i(2,2)*(tij_abcde(1,1,2,2,k)*qp_j(1,1)+& 672 tij_abcde(2,1,2,2,k)*qp_j(2,1)+& 673 tij_abcde(3,1,2,2,k)*qp_j(3,1)+& 674 tij_abcde(1,2,2,2,k)*qp_j(1,2)+& 675 tij_abcde(2,2,2,2,k)*qp_j(2,2)+& 676 tij_abcde(3,2,2,2,k)*qp_j(3,2)+& 677 tij_abcde(1,3,2,2,k)*qp_j(1,3)+& 678 tij_abcde(2,3,2,2,k)*qp_j(2,3)+& 679 tij_abcde(3,3,2,2,k)*qp_j(3,3)) 680 tmp32 = qp_i(3,2)*(tij_abcde(1,1,3,2,k)*qp_j(1,1)+& 681 tij_abcde(2,1,3,2,k)*qp_j(2,1)+& 682 tij_abcde(3,1,3,2,k)*qp_j(3,1)+& 683 tij_abcde(1,2,3,2,k)*qp_j(1,2)+& 684 tij_abcde(2,2,3,2,k)*qp_j(2,2)+& 685 tij_abcde(3,2,3,2,k)*qp_j(3,2)+& 686 tij_abcde(1,3,3,2,k)*qp_j(1,3)+& 687 tij_abcde(2,3,3,2,k)*qp_j(2,3)+& 688 tij_abcde(3,3,3,2,k)*qp_j(3,3)) 689 tmp33 = qp_i(3,3)*(tij_abcde(1,1,3,3,k)*qp_j(1,1)+& 690 tij_abcde(2,1,3,3,k)*qp_j(2,1)+& 691 tij_abcde(3,1,3,3,k)*qp_j(3,1)+& 692 tij_abcde(1,2,3,3,k)*qp_j(1,2)+& 693 tij_abcde(2,2,3,3,k)*qp_j(2,2)+& 694 tij_abcde(3,2,3,3,k)*qp_j(3,2)+& 695 tij_abcde(1,3,3,3,k)*qp_j(1,3)+& 696 tij_abcde(2,3,3,3,k)*qp_j(2,3)+& 697 tij_abcde(3,3,3,3,k)*qp_j(3,3)) 698 tmp12 = tmp21 699 tmp13 = tmp31 700 tmp23 = tmp32 701 fr(k) = fr(k) - fac * ( tmp11 + tmp12 + tmp13 +& 702 tmp21 + tmp22 + tmp23 +& 703 tmp31 + tmp32 + tmp33 ) 704 END DO 705 END IF 706 ! Electric field 707 IF (do_efield) THEN 708 fac = 1.0_dp/3.0_dp 709 ! Potential 710 IF (do_efield0) THEN 711 ef0_i = ef0_i + fac*(tij_ab(1,1)*qp_j(1,1)+& 712 tij_ab(2,1)*qp_j(2,1)+& 713 tij_ab(3,1)*qp_j(3,1)+& 714 tij_ab(1,2)*qp_j(1,2)+& 715 tij_ab(2,2)*qp_j(2,2)+& 716 tij_ab(3,2)*qp_j(3,2)+& 717 tij_ab(1,3)*qp_j(1,3)+& 718 tij_ab(2,3)*qp_j(2,3)+& 719 tij_ab(3,3)*qp_j(3,3)) 720 721 ef0_j = ef0_j + fac*(tij_ab(1,1)*qp_i(1,1)+& 722 tij_ab(2,1)*qp_i(2,1)+& 723 tij_ab(3,1)*qp_i(3,1)+& 724 tij_ab(1,2)*qp_i(1,2)+& 725 tij_ab(2,2)*qp_i(2,2)+& 726 tij_ab(3,2)*qp_i(3,2)+& 727 tij_ab(1,3)*qp_i(1,3)+& 728 tij_ab(2,3)*qp_i(2,3)+& 729 tij_ab(3,3)*qp_i(3,3)) 730 END IF 731 ! Electric field 732 IF (do_efield1) THEN 733 ef1_i(1) = ef1_i(1) - fac*(tij_abc(1,1,1)*qp_j(1,1)+& 734 tij_abc(2,1,1)*qp_j(2,1)+& 735 tij_abc(3,1,1)*qp_j(3,1)+& 736 tij_abc(1,2,1)*qp_j(1,2)+& 737 tij_abc(2,2,1)*qp_j(2,2)+& 738 tij_abc(3,2,1)*qp_j(3,2)+& 739 tij_abc(1,3,1)*qp_j(1,3)+& 740 tij_abc(2,3,1)*qp_j(2,3)+& 741 tij_abc(3,3,1)*qp_j(3,3)) 742 ef1_i(2) = ef1_i(2) - fac*(tij_abc(1,1,2)*qp_j(1,1)+& 743 tij_abc(2,1,2)*qp_j(2,1)+& 744 tij_abc(3,1,2)*qp_j(3,1)+& 745 tij_abc(1,2,2)*qp_j(1,2)+& 746 tij_abc(2,2,2)*qp_j(2,2)+& 747 tij_abc(3,2,2)*qp_j(3,2)+& 748 tij_abc(1,3,2)*qp_j(1,3)+& 749 tij_abc(2,3,2)*qp_j(2,3)+& 750 tij_abc(3,3,2)*qp_j(3,3)) 751 ef1_i(3) = ef1_i(3) - fac*(tij_abc(1,1,3)*qp_j(1,1)+& 752 tij_abc(2,1,3)*qp_j(2,1)+& 753 tij_abc(3,1,3)*qp_j(3,1)+& 754 tij_abc(1,2,3)*qp_j(1,2)+& 755 tij_abc(2,2,3)*qp_j(2,2)+& 756 tij_abc(3,2,3)*qp_j(3,2)+& 757 tij_abc(1,3,3)*qp_j(1,3)+& 758 tij_abc(2,3,3)*qp_j(2,3)+& 759 tij_abc(3,3,3)*qp_j(3,3)) 760 761 ef1_j(1) = ef1_j(1) + fac*(tij_abc(1,1,1)*qp_i(1,1)+& 762 tij_abc(2,1,1)*qp_i(2,1)+& 763 tij_abc(3,1,1)*qp_i(3,1)+& 764 tij_abc(1,2,1)*qp_i(1,2)+& 765 tij_abc(2,2,1)*qp_i(2,2)+& 766 tij_abc(3,2,1)*qp_i(3,2)+& 767 tij_abc(1,3,1)*qp_i(1,3)+& 768 tij_abc(2,3,1)*qp_i(2,3)+& 769 tij_abc(3,3,1)*qp_i(3,3)) 770 ef1_j(2) = ef1_j(2) + fac*(tij_abc(1,1,2)*qp_i(1,1)+& 771 tij_abc(2,1,2)*qp_i(2,1)+& 772 tij_abc(3,1,2)*qp_i(3,1)+& 773 tij_abc(1,2,2)*qp_i(1,2)+& 774 tij_abc(2,2,2)*qp_i(2,2)+& 775 tij_abc(3,2,2)*qp_i(3,2)+& 776 tij_abc(1,3,2)*qp_i(1,3)+& 777 tij_abc(2,3,2)*qp_i(2,3)+& 778 tij_abc(3,3,2)*qp_i(3,3)) 779 ef1_j(3) = ef1_j(3) + fac*(tij_abc(1,1,3)*qp_i(1,1)+& 780 tij_abc(2,1,3)*qp_i(2,1)+& 781 tij_abc(3,1,3)*qp_i(3,1)+& 782 tij_abc(1,2,3)*qp_i(1,2)+& 783 tij_abc(2,2,3)*qp_i(2,2)+& 784 tij_abc(3,2,3)*qp_i(3,2)+& 785 tij_abc(1,3,3)*qp_i(1,3)+& 786 tij_abc(2,3,3)*qp_i(2,3)+& 787 tij_abc(3,3,3)*qp_i(3,3)) 788 END IF 789 ! Electric field gradient 790 IF (do_efield2) THEN 791 tmp11 = fac *(tij_abcd(1,1,1,1)*qp_j(1,1)+& 792 tij_abcd(2,1,1,1)*qp_j(2,1)+& 793 tij_abcd(3,1,1,1)*qp_j(3,1)+& 794 tij_abcd(1,2,1,1)*qp_j(1,2)+& 795 tij_abcd(2,2,1,1)*qp_j(2,2)+& 796 tij_abcd(3,2,1,1)*qp_j(3,2)+& 797 tij_abcd(1,3,1,1)*qp_j(1,3)+& 798 tij_abcd(2,3,1,1)*qp_j(2,3)+& 799 tij_abcd(3,3,1,1)*qp_j(3,3)) 800 tmp12 = fac *(tij_abcd(1,1,1,2)*qp_j(1,1)+& 801 tij_abcd(2,1,1,2)*qp_j(2,1)+& 802 tij_abcd(3,1,1,2)*qp_j(3,1)+& 803 tij_abcd(1,2,1,2)*qp_j(1,2)+& 804 tij_abcd(2,2,1,2)*qp_j(2,2)+& 805 tij_abcd(3,2,1,2)*qp_j(3,2)+& 806 tij_abcd(1,3,1,2)*qp_j(1,3)+& 807 tij_abcd(2,3,1,2)*qp_j(2,3)+& 808 tij_abcd(3,3,1,2)*qp_j(3,3)) 809 tmp13 = fac *(tij_abcd(1,1,1,3)*qp_j(1,1)+& 810 tij_abcd(2,1,1,3)*qp_j(2,1)+& 811 tij_abcd(3,1,1,3)*qp_j(3,1)+& 812 tij_abcd(1,2,1,3)*qp_j(1,2)+& 813 tij_abcd(2,2,1,3)*qp_j(2,2)+& 814 tij_abcd(3,2,1,3)*qp_j(3,2)+& 815 tij_abcd(1,3,1,3)*qp_j(1,3)+& 816 tij_abcd(2,3,1,3)*qp_j(2,3)+& 817 tij_abcd(3,3,1,3)*qp_j(3,3)) 818 tmp22 = fac *(tij_abcd(1,1,2,2)*qp_j(1,1)+& 819 tij_abcd(2,1,2,2)*qp_j(2,1)+& 820 tij_abcd(3,1,2,2)*qp_j(3,1)+& 821 tij_abcd(1,2,2,2)*qp_j(1,2)+& 822 tij_abcd(2,2,2,2)*qp_j(2,2)+& 823 tij_abcd(3,2,2,2)*qp_j(3,2)+& 824 tij_abcd(1,3,2,2)*qp_j(1,3)+& 825 tij_abcd(2,3,2,2)*qp_j(2,3)+& 826 tij_abcd(3,3,2,2)*qp_j(3,3)) 827 tmp23 = fac *(tij_abcd(1,1,2,3)*qp_j(1,1)+& 828 tij_abcd(2,1,2,3)*qp_j(2,1)+& 829 tij_abcd(3,1,2,3)*qp_j(3,1)+& 830 tij_abcd(1,2,2,3)*qp_j(1,2)+& 831 tij_abcd(2,2,2,3)*qp_j(2,2)+& 832 tij_abcd(3,2,2,3)*qp_j(3,2)+& 833 tij_abcd(1,3,2,3)*qp_j(1,3)+& 834 tij_abcd(2,3,2,3)*qp_j(2,3)+& 835 tij_abcd(3,3,2,3)*qp_j(3,3)) 836 tmp33 = fac *(tij_abcd(1,1,3,3)*qp_j(1,1)+& 837 tij_abcd(2,1,3,3)*qp_j(2,1)+& 838 tij_abcd(3,1,3,3)*qp_j(3,1)+& 839 tij_abcd(1,2,3,3)*qp_j(1,2)+& 840 tij_abcd(2,2,3,3)*qp_j(2,2)+& 841 tij_abcd(3,2,3,3)*qp_j(3,2)+& 842 tij_abcd(1,3,3,3)*qp_j(1,3)+& 843 tij_abcd(2,3,3,3)*qp_j(2,3)+& 844 tij_abcd(3,3,3,3)*qp_j(3,3)) 845 846 ef2_i(1,1) = ef2_i(1,1) - tmp11 847 ef2_i(1,2) = ef2_i(1,2) - tmp12 848 ef2_i(1,3) = ef2_i(1,3) - tmp13 849 ef2_i(2,1) = ef2_i(2,1) - tmp12 850 ef2_i(2,2) = ef2_i(2,2) - tmp22 851 ef2_i(2,3) = ef2_i(2,3) - tmp23 852 ef2_i(3,1) = ef2_i(3,1) - tmp13 853 ef2_i(3,2) = ef2_i(3,2) - tmp23 854 ef2_i(3,3) = ef2_i(3,3) - tmp33 855 856 tmp11 = fac *(tij_abcd(1,1,1,1)*qp_i(1,1)+& 857 tij_abcd(2,1,1,1)*qp_i(2,1)+& 858 tij_abcd(3,1,1,1)*qp_i(3,1)+& 859 tij_abcd(1,2,1,1)*qp_i(1,2)+& 860 tij_abcd(2,2,1,1)*qp_i(2,2)+& 861 tij_abcd(3,2,1,1)*qp_i(3,2)+& 862 tij_abcd(1,3,1,1)*qp_i(1,3)+& 863 tij_abcd(2,3,1,1)*qp_i(2,3)+& 864 tij_abcd(3,3,1,1)*qp_i(3,3)) 865 tmp12 = fac *(tij_abcd(1,1,1,2)*qp_i(1,1)+& 866 tij_abcd(2,1,1,2)*qp_i(2,1)+& 867 tij_abcd(3,1,1,2)*qp_i(3,1)+& 868 tij_abcd(1,2,1,2)*qp_i(1,2)+& 869 tij_abcd(2,2,1,2)*qp_i(2,2)+& 870 tij_abcd(3,2,1,2)*qp_i(3,2)+& 871 tij_abcd(1,3,1,2)*qp_i(1,3)+& 872 tij_abcd(2,3,1,2)*qp_i(2,3)+& 873 tij_abcd(3,3,1,2)*qp_i(3,3)) 874 tmp13 = fac *(tij_abcd(1,1,1,3)*qp_i(1,1)+& 875 tij_abcd(2,1,1,3)*qp_i(2,1)+& 876 tij_abcd(3,1,1,3)*qp_i(3,1)+& 877 tij_abcd(1,2,1,3)*qp_i(1,2)+& 878 tij_abcd(2,2,1,3)*qp_i(2,2)+& 879 tij_abcd(3,2,1,3)*qp_i(3,2)+& 880 tij_abcd(1,3,1,3)*qp_i(1,3)+& 881 tij_abcd(2,3,1,3)*qp_i(2,3)+& 882 tij_abcd(3,3,1,3)*qp_i(3,3)) 883 tmp22 = fac *(tij_abcd(1,1,2,2)*qp_i(1,1)+& 884 tij_abcd(2,1,2,2)*qp_i(2,1)+& 885 tij_abcd(3,1,2,2)*qp_i(3,1)+& 886 tij_abcd(1,2,2,2)*qp_i(1,2)+& 887 tij_abcd(2,2,2,2)*qp_i(2,2)+& 888 tij_abcd(3,2,2,2)*qp_i(3,2)+& 889 tij_abcd(1,3,2,2)*qp_i(1,3)+& 890 tij_abcd(2,3,2,2)*qp_i(2,3)+& 891 tij_abcd(3,3,2,2)*qp_i(3,3)) 892 tmp23 = fac *(tij_abcd(1,1,2,3)*qp_i(1,1)+& 893 tij_abcd(2,1,2,3)*qp_i(2,1)+& 894 tij_abcd(3,1,2,3)*qp_i(3,1)+& 895 tij_abcd(1,2,2,3)*qp_i(1,2)+& 896 tij_abcd(2,2,2,3)*qp_i(2,2)+& 897 tij_abcd(3,2,2,3)*qp_i(3,2)+& 898 tij_abcd(1,3,2,3)*qp_i(1,3)+& 899 tij_abcd(2,3,2,3)*qp_i(2,3)+& 900 tij_abcd(3,3,2,3)*qp_i(3,3)) 901 tmp33 = fac *(tij_abcd(1,1,3,3)*qp_i(1,1)+& 902 tij_abcd(2,1,3,3)*qp_i(2,1)+& 903 tij_abcd(3,1,3,3)*qp_i(3,1)+& 904 tij_abcd(1,2,3,3)*qp_i(1,2)+& 905 tij_abcd(2,2,3,3)*qp_i(2,2)+& 906 tij_abcd(3,2,3,3)*qp_i(3,2)+& 907 tij_abcd(1,3,3,3)*qp_i(1,3)+& 908 tij_abcd(2,3,3,3)*qp_i(2,3)+& 909 tij_abcd(3,3,3,3)*qp_i(3,3)) 910 911 ef2_j(1,1) = ef2_j(1,1) - tmp11 912 ef2_j(1,2) = ef2_j(1,2) - tmp12 913 ef2_j(1,3) = ef2_j(1,3) - tmp13 914 ef2_j(2,1) = ef2_j(2,1) - tmp12 915 ef2_j(2,2) = ef2_j(2,2) - tmp22 916 ef2_j(2,3) = ef2_j(2,3) - tmp23 917 ef2_j(3,1) = ef2_j(3,1) - tmp13 918 ef2_j(3,2) = ef2_j(3,2) - tmp23 919 ef2_j(3,3) = ef2_j(3,3) - tmp33 920 END IF 921 END IF 922 END IF 923 IF (task(3,2)) THEN 924 ! Quadrupole - Dipole 925 fac = 1.0_dp/3.0_dp 926 ! Dipole i (locally B) - Quadrupole j (locally A) 927 tmp_ij = dp_i(1)*(tij_abc(1,1,1)*qp_j(1,1)+& 928 tij_abc(2,1,1)*qp_j(2,1)+& 929 tij_abc(3,1,1)*qp_j(3,1)+& 930 tij_abc(1,2,1)*qp_j(1,2)+& 931 tij_abc(2,2,1)*qp_j(2,2)+& 932 tij_abc(3,2,1)*qp_j(3,2)+& 933 tij_abc(1,3,1)*qp_j(1,3)+& 934 tij_abc(2,3,1)*qp_j(2,3)+& 935 tij_abc(3,3,1)*qp_j(3,3))+& 936 dp_i(2)*(tij_abc(1,1,2)*qp_j(1,1)+& 937 tij_abc(2,1,2)*qp_j(2,1)+& 938 tij_abc(3,1,2)*qp_j(3,1)+& 939 tij_abc(1,2,2)*qp_j(1,2)+& 940 tij_abc(2,2,2)*qp_j(2,2)+& 941 tij_abc(3,2,2)*qp_j(3,2)+& 942 tij_abc(1,3,2)*qp_j(1,3)+& 943 tij_abc(2,3,2)*qp_j(2,3)+& 944 tij_abc(3,3,2)*qp_j(3,3))+& 945 dp_i(3)*(tij_abc(1,1,3)*qp_j(1,1)+& 946 tij_abc(2,1,3)*qp_j(2,1)+& 947 tij_abc(3,1,3)*qp_j(3,1)+& 948 tij_abc(1,2,3)*qp_j(1,2)+& 949 tij_abc(2,2,3)*qp_j(2,2)+& 950 tij_abc(3,2,3)*qp_j(3,2)+& 951 tij_abc(1,3,3)*qp_j(1,3)+& 952 tij_abc(2,3,3)*qp_j(2,3)+& 953 tij_abc(3,3,3)*qp_j(3,3)) 954 955 ! Dipole j (locally A) - Quadrupole i (locally B) 956 tmp_ji = dp_j(1)*(tij_abc(1,1,1)*qp_i(1,1)+& 957 tij_abc(2,1,1)*qp_i(2,1)+& 958 tij_abc(3,1,1)*qp_i(3,1)+& 959 tij_abc(1,2,1)*qp_i(1,2)+& 960 tij_abc(2,2,1)*qp_i(2,2)+& 961 tij_abc(3,2,1)*qp_i(3,2)+& 962 tij_abc(1,3,1)*qp_i(1,3)+& 963 tij_abc(2,3,1)*qp_i(2,3)+& 964 tij_abc(3,3,1)*qp_i(3,3))+& 965 dp_j(2)*(tij_abc(1,1,2)*qp_i(1,1)+& 966 tij_abc(2,1,2)*qp_i(2,1)+& 967 tij_abc(3,1,2)*qp_i(3,1)+& 968 tij_abc(1,2,2)*qp_i(1,2)+& 969 tij_abc(2,2,2)*qp_i(2,2)+& 970 tij_abc(3,2,2)*qp_i(3,2)+& 971 tij_abc(1,3,2)*qp_i(1,3)+& 972 tij_abc(2,3,2)*qp_i(2,3)+& 973 tij_abc(3,3,2)*qp_i(3,3))+& 974 dp_j(3)*(tij_abc(1,1,3)*qp_i(1,1)+& 975 tij_abc(2,1,3)*qp_i(2,1)+& 976 tij_abc(3,1,3)*qp_i(3,1)+& 977 tij_abc(1,2,3)*qp_i(1,2)+& 978 tij_abc(2,2,3)*qp_i(2,2)+& 979 tij_abc(3,2,3)*qp_i(3,2)+& 980 tij_abc(1,3,3)*qp_i(1,3)+& 981 tij_abc(2,3,3)*qp_i(2,3)+& 982 tij_abc(3,3,3)*qp_i(3,3)) 983 984 tmp= fac * (tmp_ij - tmp_ji) 985 eloc = eloc + tmp 986 IF (do_forces.OR.do_stress) THEN 987 DO k = 1, 3 988 ! Dipole i (locally B) - Quadrupole j (locally A) 989 tmp_ij = dp_i(1)*(tij_abcd(1,1,1,k)*qp_j(1,1)+& 990 tij_abcd(2,1,1,k)*qp_j(2,1)+& 991 tij_abcd(3,1,1,k)*qp_j(3,1)+& 992 tij_abcd(1,2,1,k)*qp_j(1,2)+& 993 tij_abcd(2,2,1,k)*qp_j(2,2)+& 994 tij_abcd(3,2,1,k)*qp_j(3,2)+& 995 tij_abcd(1,3,1,k)*qp_j(1,3)+& 996 tij_abcd(2,3,1,k)*qp_j(2,3)+& 997 tij_abcd(3,3,1,k)*qp_j(3,3))+& 998 dp_i(2)*(tij_abcd(1,1,2,k)*qp_j(1,1)+& 999 tij_abcd(2,1,2,k)*qp_j(2,1)+& 1000 tij_abcd(3,1,2,k)*qp_j(3,1)+& 1001 tij_abcd(1,2,2,k)*qp_j(1,2)+& 1002 tij_abcd(2,2,2,k)*qp_j(2,2)+& 1003 tij_abcd(3,2,2,k)*qp_j(3,2)+& 1004 tij_abcd(1,3,2,k)*qp_j(1,3)+& 1005 tij_abcd(2,3,2,k)*qp_j(2,3)+& 1006 tij_abcd(3,3,2,k)*qp_j(3,3))+& 1007 dp_i(3)*(tij_abcd(1,1,3,k)*qp_j(1,1)+& 1008 tij_abcd(2,1,3,k)*qp_j(2,1)+& 1009 tij_abcd(3,1,3,k)*qp_j(3,1)+& 1010 tij_abcd(1,2,3,k)*qp_j(1,2)+& 1011 tij_abcd(2,2,3,k)*qp_j(2,2)+& 1012 tij_abcd(3,2,3,k)*qp_j(3,2)+& 1013 tij_abcd(1,3,3,k)*qp_j(1,3)+& 1014 tij_abcd(2,3,3,k)*qp_j(2,3)+& 1015 tij_abcd(3,3,3,k)*qp_j(3,3)) 1016 1017 ! Dipole j (locally A) - Quadrupole i (locally B) 1018 tmp_ji = dp_j(1)*(tij_abcd(1,1,1,k)*qp_i(1,1)+& 1019 tij_abcd(2,1,1,k)*qp_i(2,1)+& 1020 tij_abcd(3,1,1,k)*qp_i(3,1)+& 1021 tij_abcd(1,2,1,k)*qp_i(1,2)+& 1022 tij_abcd(2,2,1,k)*qp_i(2,2)+& 1023 tij_abcd(3,2,1,k)*qp_i(3,2)+& 1024 tij_abcd(1,3,1,k)*qp_i(1,3)+& 1025 tij_abcd(2,3,1,k)*qp_i(2,3)+& 1026 tij_abcd(3,3,1,k)*qp_i(3,3))+& 1027 dp_j(2)*(tij_abcd(1,1,2,k)*qp_i(1,1)+& 1028 tij_abcd(2,1,2,k)*qp_i(2,1)+& 1029 tij_abcd(3,1,2,k)*qp_i(3,1)+& 1030 tij_abcd(1,2,2,k)*qp_i(1,2)+& 1031 tij_abcd(2,2,2,k)*qp_i(2,2)+& 1032 tij_abcd(3,2,2,k)*qp_i(3,2)+& 1033 tij_abcd(1,3,2,k)*qp_i(1,3)+& 1034 tij_abcd(2,3,2,k)*qp_i(2,3)+& 1035 tij_abcd(3,3,2,k)*qp_i(3,3))+& 1036 dp_j(3)*(tij_abcd(1,1,3,k)*qp_i(1,1)+& 1037 tij_abcd(2,1,3,k)*qp_i(2,1)+& 1038 tij_abcd(3,1,3,k)*qp_i(3,1)+& 1039 tij_abcd(1,2,3,k)*qp_i(1,2)+& 1040 tij_abcd(2,2,3,k)*qp_i(2,2)+& 1041 tij_abcd(3,2,3,k)*qp_i(3,2)+& 1042 tij_abcd(1,3,3,k)*qp_i(1,3)+& 1043 tij_abcd(2,3,3,k)*qp_i(2,3)+& 1044 tij_abcd(3,3,3,k)*qp_i(3,3)) 1045 1046 fr(k) = fr(k) - fac * (tmp_ij - tmp_ji) 1047 END DO 1048 END IF 1049 END IF 1050 IF (task(3,1)) THEN 1051 ! Quadrupole - Charge 1052 fac = 1.0_dp/3.0_dp 1053 1054 ! Quadrupole j (locally A) - Charge j (locally B) 1055 tmp_ij = ch_i * (tij_ab(1,1)*qp_j(1,1)+& 1056 tij_ab(2,1)*qp_j(2,1)+& 1057 tij_ab(3,1)*qp_j(3,1)+& 1058 tij_ab(1,2)*qp_j(1,2)+& 1059 tij_ab(2,2)*qp_j(2,2)+& 1060 tij_ab(3,2)*qp_j(3,2)+& 1061 tij_ab(1,3)*qp_j(1,3)+& 1062 tij_ab(2,3)*qp_j(2,3)+& 1063 tij_ab(3,3)*qp_j(3,3)) 1064 1065 ! Quadrupole i (locally B) - Charge j (locally A) 1066 tmp_ji = ch_j * (tij_ab(1,1)*qp_i(1,1)+& 1067 tij_ab(2,1)*qp_i(2,1)+& 1068 tij_ab(3,1)*qp_i(3,1)+& 1069 tij_ab(1,2)*qp_i(1,2)+& 1070 tij_ab(2,2)*qp_i(2,2)+& 1071 tij_ab(3,2)*qp_i(3,2)+& 1072 tij_ab(1,3)*qp_i(1,3)+& 1073 tij_ab(2,3)*qp_i(2,3)+& 1074 tij_ab(3,3)*qp_i(3,3)) 1075 1076 eloc = eloc + fac*(tmp_ij+tmp_ji) 1077 IF (do_forces.OR.do_stress) THEN 1078 DO k = 1, 3 1079 ! Quadrupole j (locally A) - Charge i (locally B) 1080 tmp_ij = ch_i * (tij_abc(1,1,k)*qp_j(1,1)+& 1081 tij_abc(2,1,k)*qp_j(2,1)+& 1082 tij_abc(3,1,k)*qp_j(3,1)+& 1083 tij_abc(1,2,k)*qp_j(1,2)+& 1084 tij_abc(2,2,k)*qp_j(2,2)+& 1085 tij_abc(3,2,k)*qp_j(3,2)+& 1086 tij_abc(1,3,k)*qp_j(1,3)+& 1087 tij_abc(2,3,k)*qp_j(2,3)+& 1088 tij_abc(3,3,k)*qp_j(3,3)) 1089 1090 ! Quadrupole i (locally B) - Charge j (locally A) 1091 tmp_ji = ch_j * (tij_abc(1,1,k)*qp_i(1,1)+& 1092 tij_abc(2,1,k)*qp_i(2,1)+& 1093 tij_abc(3,1,k)*qp_i(3,1)+& 1094 tij_abc(1,2,k)*qp_i(1,2)+& 1095 tij_abc(2,2,k)*qp_i(2,2)+& 1096 tij_abc(3,2,k)*qp_i(3,2)+& 1097 tij_abc(1,3,k)*qp_i(1,3)+& 1098 tij_abc(2,3,k)*qp_i(2,3)+& 1099 tij_abc(3,3,k)*qp_i(3,3)) 1100 1101 fr(k) = fr(k) - fac *(tmp_ij + tmp_ji) 1102 END DO 1103 END IF 1104 END IF 1105#:if store_energy 1106 energy = energy + eloc 1107#:endif 1108#:if store_forces 1109 IF (do_forces) THEN 1110 forces(1,atom_a) = forces(1,atom_a) - fr(1) 1111 forces(2,atom_a) = forces(2,atom_a) - fr(2) 1112 forces(3,atom_a) = forces(3,atom_a) - fr(3) 1113 forces(1,atom_b) = forces(1,atom_b) + fr(1) 1114 forces(2,atom_b) = forces(2,atom_b) + fr(2) 1115 forces(3,atom_b) = forces(3,atom_b) + fr(3) 1116 END IF 1117#:endif 1118 ! Electric fields 1119 IF (do_efield) THEN 1120 ! Potential 1121 IF (do_efield0) THEN 1122 efield0( atom_a) = efield0( atom_a) + ef0_j 1123 1124 efield0( atom_b) = efield0( atom_b) + ef0_i 1125 END IF 1126 ! Electric field 1127 IF (do_efield1) THEN 1128 efield1(1,atom_a) = efield1(1,atom_a) + ef1_j(1) 1129 efield1(2,atom_a) = efield1(2,atom_a) + ef1_j(2) 1130 efield1(3,atom_a) = efield1(3,atom_a) + ef1_j(3) 1131 1132 efield1(1,atom_b) = efield1(1,atom_b) + ef1_i(1) 1133 efield1(2,atom_b) = efield1(2,atom_b) + ef1_i(2) 1134 efield1(3,atom_b) = efield1(3,atom_b) + ef1_i(3) 1135 END IF 1136 ! Electric field gradient 1137 IF (do_efield2) THEN 1138 efield2(1,atom_a) = efield2(1,atom_a) + ef2_j(1,1) 1139 efield2(2,atom_a) = efield2(2,atom_a) + ef2_j(1,2) 1140 efield2(3,atom_a) = efield2(3,atom_a) + ef2_j(1,3) 1141 efield2(4,atom_a) = efield2(4,atom_a) + ef2_j(2,1) 1142 efield2(5,atom_a) = efield2(5,atom_a) + ef2_j(2,2) 1143 efield2(6,atom_a) = efield2(6,atom_a) + ef2_j(2,3) 1144 efield2(7,atom_a) = efield2(7,atom_a) + ef2_j(3,1) 1145 efield2(8,atom_a) = efield2(8,atom_a) + ef2_j(3,2) 1146 efield2(9,atom_a) = efield2(9,atom_a) + ef2_j(3,3) 1147 1148 efield2(1,atom_b) = efield2(1,atom_b) + ef2_i(1,1) 1149 efield2(2,atom_b) = efield2(2,atom_b) + ef2_i(1,2) 1150 efield2(3,atom_b) = efield2(3,atom_b) + ef2_i(1,3) 1151 efield2(4,atom_b) = efield2(4,atom_b) + ef2_i(2,1) 1152 efield2(5,atom_b) = efield2(5,atom_b) + ef2_i(2,2) 1153 efield2(6,atom_b) = efield2(6,atom_b) + ef2_i(2,3) 1154 efield2(7,atom_b) = efield2(7,atom_b) + ef2_i(3,1) 1155 efield2(8,atom_b) = efield2(8,atom_b) + ef2_i(3,2) 1156 efield2(9,atom_b) = efield2(9,atom_b) + ef2_i(3,3) 1157 END IF 1158 END IF 1159 IF (do_stress) THEN 1160 ptens11 = ptens11 + rab(1) * fr(1) 1161 ptens21 = ptens21 + rab(2) * fr(1) 1162 ptens31 = ptens31 + rab(3) * fr(1) 1163 ptens12 = ptens12 + rab(1) * fr(2) 1164 ptens22 = ptens22 + rab(2) * fr(2) 1165 ptens32 = ptens32 + rab(3) * fr(2) 1166 ptens13 = ptens13 + rab(1) * fr(3) 1167 ptens23 = ptens23 + rab(2) * fr(3) 1168 ptens33 = ptens33 + rab(3) * fr(3) 1169 END IF 1170 1171#:enddef ewalds_multipole_sr_macro 1172 1173#:endmute 1174