57#include "./base/base_uses.f90"
64 REAL(KIND=
dp),
DIMENSION(:), &
65 POINTER :: charge => null()
66 REAL(KIND=
dp),
DIMENSION(:, :), &
67 POINTER :: pos => null()
68 END TYPE charge_mono_type
69 TYPE multi_charge_type
70 TYPE(charge_mono_type),
DIMENSION(:), &
71 POINTER :: charge_typ => null()
72 END TYPE multi_charge_type
74 LOGICAL,
PRIVATE,
PARAMETER :: debug_this_module = .false.
75 LOGICAL,
PRIVATE,
PARAMETER :: debug_r_space = .false.
76 LOGICAL,
PRIVATE,
PARAMETER :: debug_g_space = .false.
77 LOGICAL,
PRIVATE,
PARAMETER :: debug_e_field = .false.
78 LOGICAL,
PRIVATE,
PARAMETER :: debug_e_field_en = .false.
79 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'ewalds_multipole'
129 cell, particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, &
130 task, do_correction_bonded, do_forces, do_stress, &
131 do_efield, radii, charges, dipoles, &
132 quadrupoles, forces_local, forces_glob, pv_local, pv_glob, efield0, efield1, &
133 efield2, iw, do_debug, atomic_kind_set, mm_section)
140 REAL(kind=
dp),
INTENT(INOUT) :: energy_local, energy_glob
141 REAL(kind=
dp),
INTENT(OUT) :: e_neut, e_self
142 LOGICAL,
DIMENSION(3),
INTENT(IN) :: task
143 LOGICAL,
INTENT(IN) :: do_correction_bonded, do_forces, &
145 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: radii, charges
146 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
147 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
148 POINTER :: quadrupoles
149 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT), &
150 OPTIONAL :: forces_local, forces_glob, pv_local, &
152 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT),
OPTIONAL :: efield0
153 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT), &
154 OPTIONAL :: efield1, efield2
155 INTEGER,
INTENT(IN) :: iw
156 LOGICAL,
INTENT(IN) :: do_debug
158 POINTER :: atomic_kind_set
161 CHARACTER(len=*),
PARAMETER :: routinen =
'ewald_multipole_evaluate'
163 INTEGER :: handle, i, j, size1, size2
164 LOGICAL :: check_debug, check_efield, check_forces, &
166 LOGICAL,
DIMENSION(3, 3) :: my_task
167 REAL(kind=
dp) :: e_bonded, e_bonded_t, e_rspace, &
168 e_rspace_t, energy_glob_t
169 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0_lr, efield0_sr
170 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: efield1_lr, efield1_sr, efield2_lr, &
176 CALL timeset(routinen, handle)
177 cpassert(
ASSOCIATED(nonbond_env))
178 check_debug = (debug_this_module .OR. debug_r_space .OR. debug_g_space .OR. debug_e_field .OR. debug_e_field_en) &
179 .EQV. debug_this_module
180 cpassert(check_debug)
181 check_forces = do_forces .EQV. (
PRESENT(forces_local) .AND.
PRESENT(forces_glob))
182 cpassert(check_forces)
183 check_efield = do_efield .EQV. (
PRESENT(efield0) .OR.
PRESENT(efield1) .OR.
PRESENT(efield2))
184 cpassert(check_efield)
186 IF (debug_this_module .AND. do_debug)
THEN
188 IF (debug_r_space)
THEN
189 CALL debug_ewald_multipoles(ewald_env, ewald_pw, nonbond_env, cell, &
190 particle_set, local_particles, iw, debug_r_space)
191 cpabort(
"Debug Multipole Requested: Real Part!")
194 IF (debug_e_field)
THEN
195 cpassert(
PRESENT(atomic_kind_set))
196 cpassert(
PRESENT(mm_section))
197 CALL debug_ewald_multipoles_fields(ewald_env, ewald_pw, nonbond_env, &
198 cell, particle_set, local_particles, radii, charges, dipoles, &
199 quadrupoles, task, iw, atomic_kind_set, mm_section)
200 cpabort(
"Debug Multipole Requested: POT+EFIELDS+GRAD!")
204 IF (debug_e_field_en)
THEN
205 CALL debug_ewald_multipoles_fields2(ewald_env, ewald_pw, nonbond_env, &
206 cell, particle_set, local_particles, radii, charges, dipoles, &
207 quadrupoles, task, iw)
208 cpabort(
"Debug Multipole Requested: POT+EFIELDS+GRAD to give the correct energy!")
218 do_task(1) = any(charges /= 0.0_dp)
220 do_task(2) = any(dipoles /= 0.0_dp)
222 do_task(3) = any(quadrupoles /= 0.0_dp)
228 my_task(j, i) = do_task(i) .AND. do_task(j)
229 my_task(i, j) = my_task(j, i)
234 NULLIFY (efield0_sr, efield0_lr, efield1_sr, efield1_lr, efield2_sr, efield2_lr)
236 IF (
PRESENT(efield0))
THEN
237 size1 =
SIZE(efield0)
238 ALLOCATE (efield0_sr(size1))
239 ALLOCATE (efield0_lr(size1))
243 IF (
PRESENT(efield1))
THEN
244 size1 =
SIZE(efield1, 1)
245 size2 =
SIZE(efield1, 2)
246 ALLOCATE (efield1_sr(size1, size2))
247 ALLOCATE (efield1_lr(size1, size2))
251 IF (
PRESENT(efield2))
THEN
252 size1 =
SIZE(efield2, 1)
253 size2 =
SIZE(efield2, 2)
254 ALLOCATE (efield2_sr(size1, size2))
255 ALLOCATE (efield2_lr(size1, size2))
263 IF ((.NOT. debug_g_space) .AND. (nonbond_env%do_nonbonded))
THEN
267 CALL ewald_multipole_sr(nonbond_env, ewald_env, atomic_kind_set, &
268 particle_set, cell, e_rspace, my_task, &
269 do_forces, do_efield, do_stress, radii, charges, dipoles, quadrupoles, &
270 forces_glob, pv_glob, efield0_sr, efield1_sr, efield2_sr)
271 energy_glob = energy_glob + e_rspace
273 IF (do_correction_bonded)
THEN
276 CALL ewald_multipole_bonded(nonbond_env, particle_set, ewald_env, &
277 cell, e_bonded, my_task, do_forces, do_efield, do_stress, &
278 charges, dipoles, quadrupoles, forces_glob, pv_glob, &
279 efield0_sr, efield1_sr, efield2_sr)
280 energy_glob = energy_glob + e_bonded
286 energy_local = 0.0_dp
287 IF (.NOT. debug_r_space)
THEN
289 CALL ewald_multipole_lr(ewald_env, ewald_pw, cell, particle_set, &
290 local_particles, energy_local, my_task, do_forces, do_efield, do_stress, &
291 charges, dipoles, quadrupoles, forces_local, pv_local, efield0_lr, efield1_lr, &
295 CALL ewald_multipole_self(ewald_env, cell, local_particles, e_self, &
296 e_neut, my_task, do_efield, radii, charges, dipoles, quadrupoles, &
297 efield0_lr, efield1_lr, efield2_lr)
302 energy_glob_t = energy_glob
303 e_rspace_t = e_rspace
304 e_bonded_t = e_bonded
305 CALL group%sum(energy_glob_t)
306 CALL group%sum(e_rspace_t)
307 CALL group%sum(e_bonded_t)
309 CALL ewald_multipole_print(iw, energy_local, e_rspace_t, e_bonded_t, e_self, e_neut)
313 IF (
PRESENT(efield0))
THEN
314 efield0 = efield0_sr + efield0_lr
315 CALL group%sum(efield0)
316 DEALLOCATE (efield0_sr)
317 DEALLOCATE (efield0_lr)
319 IF (
PRESENT(efield1))
THEN
320 efield1 = efield1_sr + efield1_lr
321 CALL group%sum(efield1)
322 DEALLOCATE (efield1_sr)
323 DEALLOCATE (efield1_lr)
325 IF (
PRESENT(efield2))
THEN
326 efield2 = efield2_sr + efield2_lr
327 CALL group%sum(efield2)
328 DEALLOCATE (efield2_sr)
329 DEALLOCATE (efield2_lr)
332 CALL timestop(handle)
359 SUBROUTINE ewald_multipole_sr(nonbond_env, ewald_env, atomic_kind_set, &
360 particle_set, cell, energy, task, &
361 do_forces, do_efield, do_stress, radii, charges, dipoles, quadrupoles, &
362 forces, pv, efield0, efield1, efield2)
366 POINTER :: atomic_kind_set
369 REAL(kind=
dp),
INTENT(INOUT) :: energy
370 LOGICAL,
DIMENSION(3, 3),
INTENT(IN) :: task
371 LOGICAL,
INTENT(IN) :: do_forces, do_efield, do_stress
372 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: radii, charges
373 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
374 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
375 POINTER :: quadrupoles
376 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT), &
377 OPTIONAL :: forces, pv
378 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0
379 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: efield1, efield2
381 CHARACTER(len=*),
PARAMETER :: routinen =
'ewald_multipole_SR'
383 INTEGER :: a, atom_a, atom_b, b, c, d, e, handle, i, iend, igrp, ikind, ilist, ipair, &
384 istart, itype_ij, itype_ji, jkind, k, kind_a, kind_b, kk, nkdamp_ij, nkdamp_ji, nkinds, &
386 INTEGER,
DIMENSION(:, :),
POINTER ::
list
387 LOGICAL :: do_efield0, do_efield1, do_efield2, &
389 REAL(kind=
dp) :: alpha, beta, ch_i, ch_j, dampa_ij, dampa_ji, dampaexpi, dampaexpj, &
390 dampfac_ij, dampfac_ji, dampfuncdiffi, dampfuncdiffj, dampfunci, dampfuncj, dampsumfi, &
391 dampsumfj, ef0_i, ef0_j, eloc,
fac, fac_ij, factorial, ir, irab2, ptens11, ptens12, &
392 ptens13, ptens21, ptens22, ptens23, ptens31, ptens32, ptens33, r, rab2, rab2_max, radius, &
393 rcut, tij, tmp, tmp1, tmp11, tmp12, tmp13, tmp2, tmp21, tmp22, tmp23, tmp31, tmp32, &
394 tmp33, tmp_ij, tmp_ji, xf
395 REAL(kind=
dp),
DIMENSION(0:5) :: f
396 REAL(kind=
dp),
DIMENSION(3) :: cell_v, cvi, damptij_a, damptji_a, dp_i, &
397 dp_j, ef1_i, ef1_j, fr, rab, tij_a
398 REAL(kind=
dp),
DIMENSION(3, 3) :: damptij_ab, damptji_ab, ef2_i, ef2_j, &
400 REAL(kind=
dp),
DIMENSION(3, 3, 3) :: tij_abc
401 REAL(kind=
dp),
DIMENSION(3, 3, 3, 3) :: tij_abcd
402 REAL(kind=
dp),
DIMENSION(3, 3, 3, 3, 3) :: tij_abcde
406 TYPE(
pos_type),
DIMENSION(:),
POINTER :: r_last_update, r_last_update_pbc
408 CALL timeset(routinen, handle)
409 NULLIFY (nonbonded, r_last_update, r_last_update_pbc)
410 do_efield0 = do_efield .AND.
ASSOCIATED(efield0)
411 do_efield1 = do_efield .AND.
ASSOCIATED(efield1)
412 do_efield2 = do_efield .AND.
ASSOCIATED(efield2)
414 ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp
415 ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp
416 ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp
420 r_last_update=r_last_update, r_last_update_pbc=r_last_update_pbc)
423 IF (debug_r_space)
THEN
424 rab2_max = huge(0.0_dp)
427 lists:
DO ilist = 1, nonbonded%nlists
428 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
429 npairs = neighbor_kind_pair%npairs
430 IF (npairs == 0) cycle lists
431 list => neighbor_kind_pair%list
432 cvi = neighbor_kind_pair%cell_vector
433 cell_v = matmul(cell%hmat, cvi)
434 kind_group_loop:
DO igrp = 1, neighbor_kind_pair%ngrp_kind
435 istart = neighbor_kind_pair%grp_kind_start(igrp)
436 iend = neighbor_kind_pair%grp_kind_end(igrp)
437 ikind = neighbor_kind_pair%ij_kind(1, igrp)
438 jkind = neighbor_kind_pair%ij_kind(2, igrp)
449 IF (
PRESENT(atomic_kind_set))
THEN
450 IF (
ASSOCIATED(atomic_kind_set(jkind)%damping))
THEN
451 damping_ij = atomic_kind_set(jkind)%damping%damp(ikind)
452 itype_ij = damping_ij%itype
453 nkdamp_ij = damping_ij%order
454 dampa_ij = damping_ij%bij
455 dampfac_ij = damping_ij%cij
458 IF (
ASSOCIATED(atomic_kind_set(ikind)%damping))
THEN
459 damping_ji = atomic_kind_set(ikind)%damping%damp(jkind)
460 itype_ji = damping_ji%itype
461 nkdamp_ji = damping_ji%order
462 dampa_ji = damping_ji%bij
463 dampfac_ji = damping_ji%cij
467 pairs:
DO ipair = istart, iend
468 IF (ipair <= neighbor_kind_pair%nscale)
THEN
471 fac_ij = neighbor_kind_pair%ei_scale(ipair)
472 IF (fac_ij <= 0) cycle pairs
476 atom_a =
list(1, ipair)
477 atom_b =
list(2, ipair)
478 kind_a = particle_set(atom_a)%atomic_kind%kind_number
479 kind_b = particle_set(atom_b)%atomic_kind%kind_number
480 IF (atom_a == atom_b) fac_ij = 0.5_dp
481 rab = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
483 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
484 IF (rab2 <= rab2_max)
THEN
485 IF (
PRESENT(radii))
THEN
486 radius = sqrt(radii(atom_a)*radii(atom_a) + radii(atom_b)*radii(atom_b))
490 IF (radius > 0.0_dp)
THEN
496 tij_ab = huge(0.0_dp)
497 tij_abc = huge(0.0_dp)
498 tij_abcd = huge(0.0_dp)
499 tij_abcde = huge(0.0_dp)
506 IF (debug_this_module .AND. debug_r_space .AND. (.NOT. debug_g_space))
THEN
511 f(0) = erf(beta*r)*ir - erf(alpha*r)*ir
518 f(i) = irab2*(f(i - 1) + tmp1*((2.0_dp*alpha**2)**i)/(
fac*alpha) - tmp2*((2.0_dp*beta**2)**i)/(
fac*beta))
523 force_eval = do_stress
526 force_eval = do_forces .OR. do_efield1
528 IF (task(2, 2)) force_eval = force_eval .OR. do_efield0
529 IF (task(1, 2) .OR. force_eval)
THEN
530 force_eval = do_stress
531 tij_a = -rab*f(1)*fac_ij
532 IF (task(1, 2)) force_eval = force_eval .OR. do_forces
534 IF (task(1, 1)) force_eval = force_eval .OR. do_efield2
535 IF (task(3, 3)) force_eval = force_eval .OR. do_efield0
536 IF (task(2, 2) .OR. task(3, 1) .OR. force_eval)
THEN
537 force_eval = do_stress
540 tmp = rab(a)*rab(b)*fac_ij
541 tij_ab(a, b) = 3.0_dp*tmp*f(2)
542 IF (a == b) tij_ab(a, b) = tij_ab(a, b) - f(1)*fac_ij
545 IF (task(2, 2) .OR. task(3, 1)) force_eval = force_eval .OR. do_forces
547 IF (task(2, 2)) force_eval = force_eval .OR. do_efield2
548 IF (task(3, 3)) force_eval = force_eval .OR. do_efield1
549 IF (task(3, 2) .OR. force_eval)
THEN
550 force_eval = do_stress
554 tmp = rab(a)*rab(b)*rab(c)*fac_ij
555 tij_abc(a, b, c) = -15.0_dp*tmp*f(3)
556 tmp = 3.0_dp*f(2)*fac_ij
557 IF (a == b) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(c)
558 IF (a == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(b)
559 IF (b == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(a)
563 IF (task(3, 2)) force_eval = force_eval .OR. do_forces
565 IF (task(3, 3) .OR. force_eval)
THEN
566 force_eval = do_stress
571 tmp = rab(a)*rab(b)*rab(c)*rab(d)*fac_ij
572 tij_abcd(a, b, c, d) = 105.0_dp*tmp*f(4)
573 tmp1 = 15.0_dp*f(3)*fac_ij
574 tmp2 = 3.0_dp*f(2)*fac_ij
576 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(c)*rab(d)
577 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
580 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(d)
581 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
583 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(c)
585 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(d)
586 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
588 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(c)
589 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(b)
594 IF (task(3, 3)) force_eval = force_eval .OR. do_forces
597 force_eval = do_stress
603 tmp = rab(a)*rab(b)*rab(c)*rab(d)*rab(e)*fac_ij
604 tij_abcde(a, b, c, d, e) = -945.0_dp*tmp*f(5)
605 tmp1 = 105.0_dp*f(4)*fac_ij
606 tmp2 = 15.0_dp*f(3)*fac_ij
608 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(c)*rab(d)*rab(e)
609 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
610 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
611 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
614 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(d)*rab(e)
615 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
616 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
617 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
620 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(e)
621 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
622 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
623 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
626 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(d)
627 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
628 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
629 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
632 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(d)*rab(e)
633 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
636 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(e)
637 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
640 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(d)
641 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
643 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(e)
644 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(d)
645 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(c)
663 IF (debug_this_module)
THEN
671 IF (any(task(1, :)))
THEN
672 ch_j = charges(atom_a)
673 ch_i = charges(atom_b)
675 IF (any(task(2, :)))
THEN
676 dp_j = dipoles(:, atom_a)
677 dp_i = dipoles(:, atom_b)
679 IF (any(task(3, :)))
THEN
680 qp_j = quadrupoles(:, :, atom_a)
681 qp_i = quadrupoles(:, :, atom_b)
685 eloc = eloc + ch_i*tij*ch_j
687 IF (do_forces .OR. do_stress)
THEN
688 fr(1) = fr(1) - ch_j*tij_a(1)*ch_i
689 fr(2) = fr(2) - ch_j*tij_a(2)*ch_i
690 fr(3) = fr(3) - ch_j*tij_a(3)*ch_i
696 ef0_i = ef0_i + tij*ch_j
698 ef0_j = ef0_j + tij*ch_i
702 ef1_i(1) = ef1_i(1) - tij_a(1)*ch_j
703 ef1_i(2) = ef1_i(2) - tij_a(2)*ch_j
704 ef1_i(3) = ef1_i(3) - tij_a(3)*ch_j
706 ef1_j(1) = ef1_j(1) + tij_a(1)*ch_i
707 ef1_j(2) = ef1_j(2) + tij_a(2)*ch_i
708 ef1_j(3) = ef1_j(3) + tij_a(3)*ch_i
714 ef2_i(1, 1) = ef2_i(1, 1) - tij_ab(1, 1)*ch_j
715 ef2_i(2, 1) = ef2_i(2, 1) - tij_ab(2, 1)*ch_j
716 ef2_i(3, 1) = ef2_i(3, 1) - tij_ab(3, 1)*ch_j
717 ef2_i(1, 2) = ef2_i(1, 2) - tij_ab(1, 2)*ch_j
718 ef2_i(2, 2) = ef2_i(2, 2) - tij_ab(2, 2)*ch_j
719 ef2_i(3, 2) = ef2_i(3, 2) - tij_ab(3, 2)*ch_j
720 ef2_i(1, 3) = ef2_i(1, 3) - tij_ab(1, 3)*ch_j
721 ef2_i(2, 3) = ef2_i(2, 3) - tij_ab(2, 3)*ch_j
722 ef2_i(3, 3) = ef2_i(3, 3) - tij_ab(3, 3)*ch_j
724 ef2_j(1, 1) = ef2_j(1, 1) - tij_ab(1, 1)*ch_i
725 ef2_j(2, 1) = ef2_j(2, 1) - tij_ab(2, 1)*ch_i
726 ef2_j(3, 1) = ef2_j(3, 1) - tij_ab(3, 1)*ch_i
727 ef2_j(1, 2) = ef2_j(1, 2) - tij_ab(1, 2)*ch_i
728 ef2_j(2, 2) = ef2_j(2, 2) - tij_ab(2, 2)*ch_i
729 ef2_j(3, 2) = ef2_j(3, 2) - tij_ab(3, 2)*ch_i
730 ef2_j(1, 3) = ef2_j(1, 3) - tij_ab(1, 3)*ch_i
731 ef2_j(2, 3) = ef2_j(2, 3) - tij_ab(2, 3)*ch_i
732 ef2_j(3, 3) = ef2_j(3, 3) - tij_ab(3, 3)*ch_i
738 tmp = -(dp_i(1)*(tij_ab(1, 1)*dp_j(1) + &
739 tij_ab(2, 1)*dp_j(2) + &
740 tij_ab(3, 1)*dp_j(3)) + &
741 dp_i(2)*(tij_ab(1, 2)*dp_j(1) + &
742 tij_ab(2, 2)*dp_j(2) + &
743 tij_ab(3, 2)*dp_j(3)) + &
744 dp_i(3)*(tij_ab(1, 3)*dp_j(1) + &
745 tij_ab(2, 3)*dp_j(2) + &
746 tij_ab(3, 3)*dp_j(3)))
749 IF (do_forces .OR. do_stress)
THEN
751 fr(k) = fr(k) + dp_i(1)*(tij_abc(1, 1, k)*dp_j(1) + &
752 tij_abc(2, 1, k)*dp_j(2) + &
753 tij_abc(3, 1, k)*dp_j(3)) &
754 + dp_i(2)*(tij_abc(1, 2, k)*dp_j(1) + &
755 tij_abc(2, 2, k)*dp_j(2) + &
756 tij_abc(3, 2, k)*dp_j(3)) &
757 + dp_i(3)*(tij_abc(1, 3, k)*dp_j(1) + &
758 tij_abc(2, 3, k)*dp_j(2) + &
759 tij_abc(3, 3, k)*dp_j(3))
766 ef0_i = ef0_i - (tij_a(1)*dp_j(1) + &
770 ef0_j = ef0_j + (tij_a(1)*dp_i(1) + &
776 ef1_i(1) = ef1_i(1) + (tij_ab(1, 1)*dp_j(1) + &
777 tij_ab(2, 1)*dp_j(2) + &
778 tij_ab(3, 1)*dp_j(3))
779 ef1_i(2) = ef1_i(2) + (tij_ab(1, 2)*dp_j(1) + &
780 tij_ab(2, 2)*dp_j(2) + &
781 tij_ab(3, 2)*dp_j(3))
782 ef1_i(3) = ef1_i(3) + (tij_ab(1, 3)*dp_j(1) + &
783 tij_ab(2, 3)*dp_j(2) + &
784 tij_ab(3, 3)*dp_j(3))
786 ef1_j(1) = ef1_j(1) + (tij_ab(1, 1)*dp_i(1) + &
787 tij_ab(2, 1)*dp_i(2) + &
788 tij_ab(3, 1)*dp_i(3))
789 ef1_j(2) = ef1_j(2) + (tij_ab(1, 2)*dp_i(1) + &
790 tij_ab(2, 2)*dp_i(2) + &
791 tij_ab(3, 2)*dp_i(3))
792 ef1_j(3) = ef1_j(3) + (tij_ab(1, 3)*dp_i(1) + &
793 tij_ab(2, 3)*dp_i(2) + &
794 tij_ab(3, 3)*dp_i(3))
798 ef2_i(1, 1) = ef2_i(1, 1) + (tij_abc(1, 1, 1)*dp_j(1) + &
799 tij_abc(2, 1, 1)*dp_j(2) + &
800 tij_abc(3, 1, 1)*dp_j(3))
801 ef2_i(1, 2) = ef2_i(1, 2) + (tij_abc(1, 1, 2)*dp_j(1) + &
802 tij_abc(2, 1, 2)*dp_j(2) + &
803 tij_abc(3, 1, 2)*dp_j(3))
804 ef2_i(1, 3) = ef2_i(1, 3) + (tij_abc(1, 1, 3)*dp_j(1) + &
805 tij_abc(2, 1, 3)*dp_j(2) + &
806 tij_abc(3, 1, 3)*dp_j(3))
807 ef2_i(2, 1) = ef2_i(2, 1) + (tij_abc(1, 2, 1)*dp_j(1) + &
808 tij_abc(2, 2, 1)*dp_j(2) + &
809 tij_abc(3, 2, 1)*dp_j(3))
810 ef2_i(2, 2) = ef2_i(2, 2) + (tij_abc(1, 2, 2)*dp_j(1) + &
811 tij_abc(2, 2, 2)*dp_j(2) + &
812 tij_abc(3, 2, 2)*dp_j(3))
813 ef2_i(2, 3) = ef2_i(2, 3) + (tij_abc(1, 2, 3)*dp_j(1) + &
814 tij_abc(2, 2, 3)*dp_j(2) + &
815 tij_abc(3, 2, 3)*dp_j(3))
816 ef2_i(3, 1) = ef2_i(3, 1) + (tij_abc(1, 3, 1)*dp_j(1) + &
817 tij_abc(2, 3, 1)*dp_j(2) + &
818 tij_abc(3, 3, 1)*dp_j(3))
819 ef2_i(3, 2) = ef2_i(3, 2) + (tij_abc(1, 3, 2)*dp_j(1) + &
820 tij_abc(2, 3, 2)*dp_j(2) + &
821 tij_abc(3, 3, 2)*dp_j(3))
822 ef2_i(3, 3) = ef2_i(3, 3) + (tij_abc(1, 3, 3)*dp_j(1) + &
823 tij_abc(2, 3, 3)*dp_j(2) + &
824 tij_abc(3, 3, 3)*dp_j(3))
826 ef2_j(1, 1) = ef2_j(1, 1) - (tij_abc(1, 1, 1)*dp_i(1) + &
827 tij_abc(2, 1, 1)*dp_i(2) + &
828 tij_abc(3, 1, 1)*dp_i(3))
829 ef2_j(1, 2) = ef2_j(1, 2) - (tij_abc(1, 1, 2)*dp_i(1) + &
830 tij_abc(2, 1, 2)*dp_i(2) + &
831 tij_abc(3, 1, 2)*dp_i(3))
832 ef2_j(1, 3) = ef2_j(1, 3) - (tij_abc(1, 1, 3)*dp_i(1) + &
833 tij_abc(2, 1, 3)*dp_i(2) + &
834 tij_abc(3, 1, 3)*dp_i(3))
835 ef2_j(2, 1) = ef2_j(2, 1) - (tij_abc(1, 2, 1)*dp_i(1) + &
836 tij_abc(2, 2, 1)*dp_i(2) + &
837 tij_abc(3, 2, 1)*dp_i(3))
838 ef2_j(2, 2) = ef2_j(2, 2) - (tij_abc(1, 2, 2)*dp_i(1) + &
839 tij_abc(2, 2, 2)*dp_i(2) + &
840 tij_abc(3, 2, 2)*dp_i(3))
841 ef2_j(2, 3) = ef2_j(2, 3) - (tij_abc(1, 2, 3)*dp_i(1) + &
842 tij_abc(2, 2, 3)*dp_i(2) + &
843 tij_abc(3, 2, 3)*dp_i(3))
844 ef2_j(3, 1) = ef2_j(3, 1) - (tij_abc(1, 3, 1)*dp_i(1) + &
845 tij_abc(2, 3, 1)*dp_i(2) + &
846 tij_abc(3, 3, 1)*dp_i(3))
847 ef2_j(3, 2) = ef2_j(3, 2) - (tij_abc(1, 3, 2)*dp_i(1) + &
848 tij_abc(2, 3, 2)*dp_i(2) + &
849 tij_abc(3, 3, 2)*dp_i(3))
850 ef2_j(3, 3) = ef2_j(3, 3) - (tij_abc(1, 3, 3)*dp_i(1) + &
851 tij_abc(2, 3, 3)*dp_i(2) + &
852 tij_abc(3, 3, 3)*dp_i(3))
858 tmp = ch_j*(tij_a(1)*dp_i(1) + &
861 - ch_i*(tij_a(1)*dp_j(1) + &
866 IF (do_forces .OR. do_stress)
THEN
868 fr(k) = fr(k) - ch_j*(tij_ab(1, k)*dp_i(1) + &
869 tij_ab(2, k)*dp_i(2) + &
870 tij_ab(3, k)*dp_i(3)) &
871 + ch_i*(tij_ab(1, k)*dp_j(1) + &
872 tij_ab(2, k)*dp_j(2) + &
873 tij_ab(3, k)*dp_j(3))
880 tmp11 = qp_i(1, 1)*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
881 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
882 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
883 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
884 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
885 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
886 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
887 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
888 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
889 tmp21 = qp_i(2, 1)*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
890 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
891 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
892 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
893 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
894 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
895 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
896 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
897 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
898 tmp31 = qp_i(3, 1)*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
899 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
900 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
901 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
902 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
903 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
904 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
905 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
906 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
907 tmp22 = qp_i(2, 2)*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
908 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
909 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
910 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
911 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
912 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
913 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
914 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
915 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
916 tmp32 = qp_i(3, 2)*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
917 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
918 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
919 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
920 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
921 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
922 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
923 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
924 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
925 tmp33 = qp_i(3, 3)*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
926 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
927 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
928 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
929 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
930 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
931 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
932 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
933 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
937 tmp = tmp11 + tmp12 + tmp13 + &
938 tmp21 + tmp22 + tmp23 + &
939 tmp31 + tmp32 + tmp33
941 eloc = eloc +
fac*tmp
943 IF (do_forces .OR. do_stress)
THEN
945 tmp11 = qp_i(1, 1)*(tij_abcde(1, 1, 1, 1, k)*qp_j(1, 1) + &
946 tij_abcde(2, 1, 1, 1, k)*qp_j(2, 1) + &
947 tij_abcde(3, 1, 1, 1, k)*qp_j(3, 1) + &
948 tij_abcde(1, 2, 1, 1, k)*qp_j(1, 2) + &
949 tij_abcde(2, 2, 1, 1, k)*qp_j(2, 2) + &
950 tij_abcde(3, 2, 1, 1, k)*qp_j(3, 2) + &
951 tij_abcde(1, 3, 1, 1, k)*qp_j(1, 3) + &
952 tij_abcde(2, 3, 1, 1, k)*qp_j(2, 3) + &
953 tij_abcde(3, 3, 1, 1, k)*qp_j(3, 3))
954 tmp21 = qp_i(2, 1)*(tij_abcde(1, 1, 2, 1, k)*qp_j(1, 1) + &
955 tij_abcde(2, 1, 2, 1, k)*qp_j(2, 1) + &
956 tij_abcde(3, 1, 2, 1, k)*qp_j(3, 1) + &
957 tij_abcde(1, 2, 2, 1, k)*qp_j(1, 2) + &
958 tij_abcde(2, 2, 2, 1, k)*qp_j(2, 2) + &
959 tij_abcde(3, 2, 2, 1, k)*qp_j(3, 2) + &
960 tij_abcde(1, 3, 2, 1, k)*qp_j(1, 3) + &
961 tij_abcde(2, 3, 2, 1, k)*qp_j(2, 3) + &
962 tij_abcde(3, 3, 2, 1, k)*qp_j(3, 3))
963 tmp31 = qp_i(3, 1)*(tij_abcde(1, 1, 3, 1, k)*qp_j(1, 1) + &
964 tij_abcde(2, 1, 3, 1, k)*qp_j(2, 1) + &
965 tij_abcde(3, 1, 3, 1, k)*qp_j(3, 1) + &
966 tij_abcde(1, 2, 3, 1, k)*qp_j(1, 2) + &
967 tij_abcde(2, 2, 3, 1, k)*qp_j(2, 2) + &
968 tij_abcde(3, 2, 3, 1, k)*qp_j(3, 2) + &
969 tij_abcde(1, 3, 3, 1, k)*qp_j(1, 3) + &
970 tij_abcde(2, 3, 3, 1, k)*qp_j(2, 3) + &
971 tij_abcde(3, 3, 3, 1, k)*qp_j(3, 3))
972 tmp22 = qp_i(2, 2)*(tij_abcde(1, 1, 2, 2, k)*qp_j(1, 1) + &
973 tij_abcde(2, 1, 2, 2, k)*qp_j(2, 1) + &
974 tij_abcde(3, 1, 2, 2, k)*qp_j(3, 1) + &
975 tij_abcde(1, 2, 2, 2, k)*qp_j(1, 2) + &
976 tij_abcde(2, 2, 2, 2, k)*qp_j(2, 2) + &
977 tij_abcde(3, 2, 2, 2, k)*qp_j(3, 2) + &
978 tij_abcde(1, 3, 2, 2, k)*qp_j(1, 3) + &
979 tij_abcde(2, 3, 2, 2, k)*qp_j(2, 3) + &
980 tij_abcde(3, 3, 2, 2, k)*qp_j(3, 3))
981 tmp32 = qp_i(3, 2)*(tij_abcde(1, 1, 3, 2, k)*qp_j(1, 1) + &
982 tij_abcde(2, 1, 3, 2, k)*qp_j(2, 1) + &
983 tij_abcde(3, 1, 3, 2, k)*qp_j(3, 1) + &
984 tij_abcde(1, 2, 3, 2, k)*qp_j(1, 2) + &
985 tij_abcde(2, 2, 3, 2, k)*qp_j(2, 2) + &
986 tij_abcde(3, 2, 3, 2, k)*qp_j(3, 2) + &
987 tij_abcde(1, 3, 3, 2, k)*qp_j(1, 3) + &
988 tij_abcde(2, 3, 3, 2, k)*qp_j(2, 3) + &
989 tij_abcde(3, 3, 3, 2, k)*qp_j(3, 3))
990 tmp33 = qp_i(3, 3)*(tij_abcde(1, 1, 3, 3, k)*qp_j(1, 1) + &
991 tij_abcde(2, 1, 3, 3, k)*qp_j(2, 1) + &
992 tij_abcde(3, 1, 3, 3, k)*qp_j(3, 1) + &
993 tij_abcde(1, 2, 3, 3, k)*qp_j(1, 2) + &
994 tij_abcde(2, 2, 3, 3, k)*qp_j(2, 2) + &
995 tij_abcde(3, 2, 3, 3, k)*qp_j(3, 2) + &
996 tij_abcde(1, 3, 3, 3, k)*qp_j(1, 3) + &
997 tij_abcde(2, 3, 3, 3, k)*qp_j(2, 3) + &
998 tij_abcde(3, 3, 3, 3, k)*qp_j(3, 3))
1002 fr(k) = fr(k) -
fac*(tmp11 + tmp12 + tmp13 + &
1003 tmp21 + tmp22 + tmp23 + &
1004 tmp31 + tmp32 + tmp33)
1011 IF (do_efield0)
THEN
1012 ef0_i = ef0_i +
fac*(tij_ab(1, 1)*qp_j(1, 1) + &
1013 tij_ab(2, 1)*qp_j(2, 1) + &
1014 tij_ab(3, 1)*qp_j(3, 1) + &
1015 tij_ab(1, 2)*qp_j(1, 2) + &
1016 tij_ab(2, 2)*qp_j(2, 2) + &
1017 tij_ab(3, 2)*qp_j(3, 2) + &
1018 tij_ab(1, 3)*qp_j(1, 3) + &
1019 tij_ab(2, 3)*qp_j(2, 3) + &
1020 tij_ab(3, 3)*qp_j(3, 3))
1022 ef0_j = ef0_j +
fac*(tij_ab(1, 1)*qp_i(1, 1) + &
1023 tij_ab(2, 1)*qp_i(2, 1) + &
1024 tij_ab(3, 1)*qp_i(3, 1) + &
1025 tij_ab(1, 2)*qp_i(1, 2) + &
1026 tij_ab(2, 2)*qp_i(2, 2) + &
1027 tij_ab(3, 2)*qp_i(3, 2) + &
1028 tij_ab(1, 3)*qp_i(1, 3) + &
1029 tij_ab(2, 3)*qp_i(2, 3) + &
1030 tij_ab(3, 3)*qp_i(3, 3))
1033 IF (do_efield1)
THEN
1034 ef1_i(1) = ef1_i(1) -
fac*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
1035 tij_abc(2, 1, 1)*qp_j(2, 1) + &
1036 tij_abc(3, 1, 1)*qp_j(3, 1) + &
1037 tij_abc(1, 2, 1)*qp_j(1, 2) + &
1038 tij_abc(2, 2, 1)*qp_j(2, 2) + &
1039 tij_abc(3, 2, 1)*qp_j(3, 2) + &
1040 tij_abc(1, 3, 1)*qp_j(1, 3) + &
1041 tij_abc(2, 3, 1)*qp_j(2, 3) + &
1042 tij_abc(3, 3, 1)*qp_j(3, 3))
1043 ef1_i(2) = ef1_i(2) -
fac*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
1044 tij_abc(2, 1, 2)*qp_j(2, 1) + &
1045 tij_abc(3, 1, 2)*qp_j(3, 1) + &
1046 tij_abc(1, 2, 2)*qp_j(1, 2) + &
1047 tij_abc(2, 2, 2)*qp_j(2, 2) + &
1048 tij_abc(3, 2, 2)*qp_j(3, 2) + &
1049 tij_abc(1, 3, 2)*qp_j(1, 3) + &
1050 tij_abc(2, 3, 2)*qp_j(2, 3) + &
1051 tij_abc(3, 3, 2)*qp_j(3, 3))
1052 ef1_i(3) = ef1_i(3) -
fac*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
1053 tij_abc(2, 1, 3)*qp_j(2, 1) + &
1054 tij_abc(3, 1, 3)*qp_j(3, 1) + &
1055 tij_abc(1, 2, 3)*qp_j(1, 2) + &
1056 tij_abc(2, 2, 3)*qp_j(2, 2) + &
1057 tij_abc(3, 2, 3)*qp_j(3, 2) + &
1058 tij_abc(1, 3, 3)*qp_j(1, 3) + &
1059 tij_abc(2, 3, 3)*qp_j(2, 3) + &
1060 tij_abc(3, 3, 3)*qp_j(3, 3))
1062 ef1_j(1) = ef1_j(1) +
fac*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
1063 tij_abc(2, 1, 1)*qp_i(2, 1) + &
1064 tij_abc(3, 1, 1)*qp_i(3, 1) + &
1065 tij_abc(1, 2, 1)*qp_i(1, 2) + &
1066 tij_abc(2, 2, 1)*qp_i(2, 2) + &
1067 tij_abc(3, 2, 1)*qp_i(3, 2) + &
1068 tij_abc(1, 3, 1)*qp_i(1, 3) + &
1069 tij_abc(2, 3, 1)*qp_i(2, 3) + &
1070 tij_abc(3, 3, 1)*qp_i(3, 3))
1071 ef1_j(2) = ef1_j(2) +
fac*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
1072 tij_abc(2, 1, 2)*qp_i(2, 1) + &
1073 tij_abc(3, 1, 2)*qp_i(3, 1) + &
1074 tij_abc(1, 2, 2)*qp_i(1, 2) + &
1075 tij_abc(2, 2, 2)*qp_i(2, 2) + &
1076 tij_abc(3, 2, 2)*qp_i(3, 2) + &
1077 tij_abc(1, 3, 2)*qp_i(1, 3) + &
1078 tij_abc(2, 3, 2)*qp_i(2, 3) + &
1079 tij_abc(3, 3, 2)*qp_i(3, 3))
1080 ef1_j(3) = ef1_j(3) +
fac*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
1081 tij_abc(2, 1, 3)*qp_i(2, 1) + &
1082 tij_abc(3, 1, 3)*qp_i(3, 1) + &
1083 tij_abc(1, 2, 3)*qp_i(1, 2) + &
1084 tij_abc(2, 2, 3)*qp_i(2, 2) + &
1085 tij_abc(3, 2, 3)*qp_i(3, 2) + &
1086 tij_abc(1, 3, 3)*qp_i(1, 3) + &
1087 tij_abc(2, 3, 3)*qp_i(2, 3) + &
1088 tij_abc(3, 3, 3)*qp_i(3, 3))
1091 IF (do_efield2)
THEN
1092 tmp11 =
fac*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
1093 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
1094 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
1095 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
1096 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
1097 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
1098 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
1099 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
1100 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
1101 tmp12 =
fac*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
1102 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
1103 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
1104 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
1105 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
1106 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
1107 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
1108 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
1109 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
1110 tmp13 =
fac*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
1111 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
1112 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
1113 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
1114 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
1115 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
1116 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
1117 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
1118 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
1119 tmp22 =
fac*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
1120 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
1121 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
1122 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
1123 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
1124 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
1125 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
1126 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
1127 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
1128 tmp23 =
fac*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
1129 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
1130 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
1131 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
1132 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
1133 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
1134 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
1135 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
1136 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
1137 tmp33 =
fac*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
1138 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
1139 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
1140 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
1141 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
1142 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
1143 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
1144 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
1145 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
1147 ef2_i(1, 1) = ef2_i(1, 1) - tmp11
1148 ef2_i(1, 2) = ef2_i(1, 2) - tmp12
1149 ef2_i(1, 3) = ef2_i(1, 3) - tmp13
1150 ef2_i(2, 1) = ef2_i(2, 1) - tmp12
1151 ef2_i(2, 2) = ef2_i(2, 2) - tmp22
1152 ef2_i(2, 3) = ef2_i(2, 3) - tmp23
1153 ef2_i(3, 1) = ef2_i(3, 1) - tmp13
1154 ef2_i(3, 2) = ef2_i(3, 2) - tmp23
1155 ef2_i(3, 3) = ef2_i(3, 3) - tmp33
1157 tmp11 =
fac*(tij_abcd(1, 1, 1, 1)*qp_i(1, 1) + &
1158 tij_abcd(2, 1, 1, 1)*qp_i(2, 1) + &
1159 tij_abcd(3, 1, 1, 1)*qp_i(3, 1) + &
1160 tij_abcd(1, 2, 1, 1)*qp_i(1, 2) + &
1161 tij_abcd(2, 2, 1, 1)*qp_i(2, 2) + &
1162 tij_abcd(3, 2, 1, 1)*qp_i(3, 2) + &
1163 tij_abcd(1, 3, 1, 1)*qp_i(1, 3) + &
1164 tij_abcd(2, 3, 1, 1)*qp_i(2, 3) + &
1165 tij_abcd(3, 3, 1, 1)*qp_i(3, 3))
1166 tmp12 =
fac*(tij_abcd(1, 1, 1, 2)*qp_i(1, 1) + &
1167 tij_abcd(2, 1, 1, 2)*qp_i(2, 1) + &
1168 tij_abcd(3, 1, 1, 2)*qp_i(3, 1) + &
1169 tij_abcd(1, 2, 1, 2)*qp_i(1, 2) + &
1170 tij_abcd(2, 2, 1, 2)*qp_i(2, 2) + &
1171 tij_abcd(3, 2, 1, 2)*qp_i(3, 2) + &
1172 tij_abcd(1, 3, 1, 2)*qp_i(1, 3) + &
1173 tij_abcd(2, 3, 1, 2)*qp_i(2, 3) + &
1174 tij_abcd(3, 3, 1, 2)*qp_i(3, 3))
1175 tmp13 =
fac*(tij_abcd(1, 1, 1, 3)*qp_i(1, 1) + &
1176 tij_abcd(2, 1, 1, 3)*qp_i(2, 1) + &
1177 tij_abcd(3, 1, 1, 3)*qp_i(3, 1) + &
1178 tij_abcd(1, 2, 1, 3)*qp_i(1, 2) + &
1179 tij_abcd(2, 2, 1, 3)*qp_i(2, 2) + &
1180 tij_abcd(3, 2, 1, 3)*qp_i(3, 2) + &
1181 tij_abcd(1, 3, 1, 3)*qp_i(1, 3) + &
1182 tij_abcd(2, 3, 1, 3)*qp_i(2, 3) + &
1183 tij_abcd(3, 3, 1, 3)*qp_i(3, 3))
1184 tmp22 =
fac*(tij_abcd(1, 1, 2, 2)*qp_i(1, 1) + &
1185 tij_abcd(2, 1, 2, 2)*qp_i(2, 1) + &
1186 tij_abcd(3, 1, 2, 2)*qp_i(3, 1) + &
1187 tij_abcd(1, 2, 2, 2)*qp_i(1, 2) + &
1188 tij_abcd(2, 2, 2, 2)*qp_i(2, 2) + &
1189 tij_abcd(3, 2, 2, 2)*qp_i(3, 2) + &
1190 tij_abcd(1, 3, 2, 2)*qp_i(1, 3) + &
1191 tij_abcd(2, 3, 2, 2)*qp_i(2, 3) + &
1192 tij_abcd(3, 3, 2, 2)*qp_i(3, 3))
1193 tmp23 =
fac*(tij_abcd(1, 1, 2, 3)*qp_i(1, 1) + &
1194 tij_abcd(2, 1, 2, 3)*qp_i(2, 1) + &
1195 tij_abcd(3, 1, 2, 3)*qp_i(3, 1) + &
1196 tij_abcd(1, 2, 2, 3)*qp_i(1, 2) + &
1197 tij_abcd(2, 2, 2, 3)*qp_i(2, 2) + &
1198 tij_abcd(3, 2, 2, 3)*qp_i(3, 2) + &
1199 tij_abcd(1, 3, 2, 3)*qp_i(1, 3) + &
1200 tij_abcd(2, 3, 2, 3)*qp_i(2, 3) + &
1201 tij_abcd(3, 3, 2, 3)*qp_i(3, 3))
1202 tmp33 =
fac*(tij_abcd(1, 1, 3, 3)*qp_i(1, 1) + &
1203 tij_abcd(2, 1, 3, 3)*qp_i(2, 1) + &
1204 tij_abcd(3, 1, 3, 3)*qp_i(3, 1) + &
1205 tij_abcd(1, 2, 3, 3)*qp_i(1, 2) + &
1206 tij_abcd(2, 2, 3, 3)*qp_i(2, 2) + &
1207 tij_abcd(3, 2, 3, 3)*qp_i(3, 2) + &
1208 tij_abcd(1, 3, 3, 3)*qp_i(1, 3) + &
1209 tij_abcd(2, 3, 3, 3)*qp_i(2, 3) + &
1210 tij_abcd(3, 3, 3, 3)*qp_i(3, 3))
1212 ef2_j(1, 1) = ef2_j(1, 1) - tmp11
1213 ef2_j(1, 2) = ef2_j(1, 2) - tmp12
1214 ef2_j(1, 3) = ef2_j(1, 3) - tmp13
1215 ef2_j(2, 1) = ef2_j(2, 1) - tmp12
1216 ef2_j(2, 2) = ef2_j(2, 2) - tmp22
1217 ef2_j(2, 3) = ef2_j(2, 3) - tmp23
1218 ef2_j(3, 1) = ef2_j(3, 1) - tmp13
1219 ef2_j(3, 2) = ef2_j(3, 2) - tmp23
1220 ef2_j(3, 3) = ef2_j(3, 3) - tmp33
1224 IF (task(3, 2))
THEN
1228 tmp_ij = dp_i(1)*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
1229 tij_abc(2, 1, 1)*qp_j(2, 1) + &
1230 tij_abc(3, 1, 1)*qp_j(3, 1) + &
1231 tij_abc(1, 2, 1)*qp_j(1, 2) + &
1232 tij_abc(2, 2, 1)*qp_j(2, 2) + &
1233 tij_abc(3, 2, 1)*qp_j(3, 2) + &
1234 tij_abc(1, 3, 1)*qp_j(1, 3) + &
1235 tij_abc(2, 3, 1)*qp_j(2, 3) + &
1236 tij_abc(3, 3, 1)*qp_j(3, 3)) + &
1237 dp_i(2)*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
1238 tij_abc(2, 1, 2)*qp_j(2, 1) + &
1239 tij_abc(3, 1, 2)*qp_j(3, 1) + &
1240 tij_abc(1, 2, 2)*qp_j(1, 2) + &
1241 tij_abc(2, 2, 2)*qp_j(2, 2) + &
1242 tij_abc(3, 2, 2)*qp_j(3, 2) + &
1243 tij_abc(1, 3, 2)*qp_j(1, 3) + &
1244 tij_abc(2, 3, 2)*qp_j(2, 3) + &
1245 tij_abc(3, 3, 2)*qp_j(3, 3)) + &
1246 dp_i(3)*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
1247 tij_abc(2, 1, 3)*qp_j(2, 1) + &
1248 tij_abc(3, 1, 3)*qp_j(3, 1) + &
1249 tij_abc(1, 2, 3)*qp_j(1, 2) + &
1250 tij_abc(2, 2, 3)*qp_j(2, 2) + &
1251 tij_abc(3, 2, 3)*qp_j(3, 2) + &
1252 tij_abc(1, 3, 3)*qp_j(1, 3) + &
1253 tij_abc(2, 3, 3)*qp_j(2, 3) + &
1254 tij_abc(3, 3, 3)*qp_j(3, 3))
1257 tmp_ji = dp_j(1)*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
1258 tij_abc(2, 1, 1)*qp_i(2, 1) + &
1259 tij_abc(3, 1, 1)*qp_i(3, 1) + &
1260 tij_abc(1, 2, 1)*qp_i(1, 2) + &
1261 tij_abc(2, 2, 1)*qp_i(2, 2) + &
1262 tij_abc(3, 2, 1)*qp_i(3, 2) + &
1263 tij_abc(1, 3, 1)*qp_i(1, 3) + &
1264 tij_abc(2, 3, 1)*qp_i(2, 3) + &
1265 tij_abc(3, 3, 1)*qp_i(3, 3)) + &
1266 dp_j(2)*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
1267 tij_abc(2, 1, 2)*qp_i(2, 1) + &
1268 tij_abc(3, 1, 2)*qp_i(3, 1) + &
1269 tij_abc(1, 2, 2)*qp_i(1, 2) + &
1270 tij_abc(2, 2, 2)*qp_i(2, 2) + &
1271 tij_abc(3, 2, 2)*qp_i(3, 2) + &
1272 tij_abc(1, 3, 2)*qp_i(1, 3) + &
1273 tij_abc(2, 3, 2)*qp_i(2, 3) + &
1274 tij_abc(3, 3, 2)*qp_i(3, 3)) + &
1275 dp_j(3)*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
1276 tij_abc(2, 1, 3)*qp_i(2, 1) + &
1277 tij_abc(3, 1, 3)*qp_i(3, 1) + &
1278 tij_abc(1, 2, 3)*qp_i(1, 2) + &
1279 tij_abc(2, 2, 3)*qp_i(2, 2) + &
1280 tij_abc(3, 2, 3)*qp_i(3, 2) + &
1281 tij_abc(1, 3, 3)*qp_i(1, 3) + &
1282 tij_abc(2, 3, 3)*qp_i(2, 3) + &
1283 tij_abc(3, 3, 3)*qp_i(3, 3))
1285 tmp =
fac*(tmp_ij - tmp_ji)
1287 IF (do_forces .OR. do_stress)
THEN
1290 tmp_ij = dp_i(1)*(tij_abcd(1, 1, 1, k)*qp_j(1, 1) + &
1291 tij_abcd(2, 1, 1, k)*qp_j(2, 1) + &
1292 tij_abcd(3, 1, 1, k)*qp_j(3, 1) + &
1293 tij_abcd(1, 2, 1, k)*qp_j(1, 2) + &
1294 tij_abcd(2, 2, 1, k)*qp_j(2, 2) + &
1295 tij_abcd(3, 2, 1, k)*qp_j(3, 2) + &
1296 tij_abcd(1, 3, 1, k)*qp_j(1, 3) + &
1297 tij_abcd(2, 3, 1, k)*qp_j(2, 3) + &
1298 tij_abcd(3, 3, 1, k)*qp_j(3, 3)) + &
1299 dp_i(2)*(tij_abcd(1, 1, 2, k)*qp_j(1, 1) + &
1300 tij_abcd(2, 1, 2, k)*qp_j(2, 1) + &
1301 tij_abcd(3, 1, 2, k)*qp_j(3, 1) + &
1302 tij_abcd(1, 2, 2, k)*qp_j(1, 2) + &
1303 tij_abcd(2, 2, 2, k)*qp_j(2, 2) + &
1304 tij_abcd(3, 2, 2, k)*qp_j(3, 2) + &
1305 tij_abcd(1, 3, 2, k)*qp_j(1, 3) + &
1306 tij_abcd(2, 3, 2, k)*qp_j(2, 3) + &
1307 tij_abcd(3, 3, 2, k)*qp_j(3, 3)) + &
1308 dp_i(3)*(tij_abcd(1, 1, 3, k)*qp_j(1, 1) + &
1309 tij_abcd(2, 1, 3, k)*qp_j(2, 1) + &
1310 tij_abcd(3, 1, 3, k)*qp_j(3, 1) + &
1311 tij_abcd(1, 2, 3, k)*qp_j(1, 2) + &
1312 tij_abcd(2, 2, 3, k)*qp_j(2, 2) + &
1313 tij_abcd(3, 2, 3, k)*qp_j(3, 2) + &
1314 tij_abcd(1, 3, 3, k)*qp_j(1, 3) + &
1315 tij_abcd(2, 3, 3, k)*qp_j(2, 3) + &
1316 tij_abcd(3, 3, 3, k)*qp_j(3, 3))
1319 tmp_ji = dp_j(1)*(tij_abcd(1, 1, 1, k)*qp_i(1, 1) + &
1320 tij_abcd(2, 1, 1, k)*qp_i(2, 1) + &
1321 tij_abcd(3, 1, 1, k)*qp_i(3, 1) + &
1322 tij_abcd(1, 2, 1, k)*qp_i(1, 2) + &
1323 tij_abcd(2, 2, 1, k)*qp_i(2, 2) + &
1324 tij_abcd(3, 2, 1, k)*qp_i(3, 2) + &
1325 tij_abcd(1, 3, 1, k)*qp_i(1, 3) + &
1326 tij_abcd(2, 3, 1, k)*qp_i(2, 3) + &
1327 tij_abcd(3, 3, 1, k)*qp_i(3, 3)) + &
1328 dp_j(2)*(tij_abcd(1, 1, 2, k)*qp_i(1, 1) + &
1329 tij_abcd(2, 1, 2, k)*qp_i(2, 1) + &
1330 tij_abcd(3, 1, 2, k)*qp_i(3, 1) + &
1331 tij_abcd(1, 2, 2, k)*qp_i(1, 2) + &
1332 tij_abcd(2, 2, 2, k)*qp_i(2, 2) + &
1333 tij_abcd(3, 2, 2, k)*qp_i(3, 2) + &
1334 tij_abcd(1, 3, 2, k)*qp_i(1, 3) + &
1335 tij_abcd(2, 3, 2, k)*qp_i(2, 3) + &
1336 tij_abcd(3, 3, 2, k)*qp_i(3, 3)) + &
1337 dp_j(3)*(tij_abcd(1, 1, 3, k)*qp_i(1, 1) + &
1338 tij_abcd(2, 1, 3, k)*qp_i(2, 1) + &
1339 tij_abcd(3, 1, 3, k)*qp_i(3, 1) + &
1340 tij_abcd(1, 2, 3, k)*qp_i(1, 2) + &
1341 tij_abcd(2, 2, 3, k)*qp_i(2, 2) + &
1342 tij_abcd(3, 2, 3, k)*qp_i(3, 2) + &
1343 tij_abcd(1, 3, 3, k)*qp_i(1, 3) + &
1344 tij_abcd(2, 3, 3, k)*qp_i(2, 3) + &
1345 tij_abcd(3, 3, 3, k)*qp_i(3, 3))
1347 fr(k) = fr(k) -
fac*(tmp_ij - tmp_ji)
1351 IF (task(3, 1))
THEN
1356 tmp_ij = ch_i*(tij_ab(1, 1)*qp_j(1, 1) + &
1357 tij_ab(2, 1)*qp_j(2, 1) + &
1358 tij_ab(3, 1)*qp_j(3, 1) + &
1359 tij_ab(1, 2)*qp_j(1, 2) + &
1360 tij_ab(2, 2)*qp_j(2, 2) + &
1361 tij_ab(3, 2)*qp_j(3, 2) + &
1362 tij_ab(1, 3)*qp_j(1, 3) + &
1363 tij_ab(2, 3)*qp_j(2, 3) + &
1364 tij_ab(3, 3)*qp_j(3, 3))
1367 tmp_ji = ch_j*(tij_ab(1, 1)*qp_i(1, 1) + &
1368 tij_ab(2, 1)*qp_i(2, 1) + &
1369 tij_ab(3, 1)*qp_i(3, 1) + &
1370 tij_ab(1, 2)*qp_i(1, 2) + &
1371 tij_ab(2, 2)*qp_i(2, 2) + &
1372 tij_ab(3, 2)*qp_i(3, 2) + &
1373 tij_ab(1, 3)*qp_i(1, 3) + &
1374 tij_ab(2, 3)*qp_i(2, 3) + &
1375 tij_ab(3, 3)*qp_i(3, 3))
1377 eloc = eloc +
fac*(tmp_ij + tmp_ji)
1378 IF (do_forces .OR. do_stress)
THEN
1381 tmp_ij = ch_i*(tij_abc(1, 1, k)*qp_j(1, 1) + &
1382 tij_abc(2, 1, k)*qp_j(2, 1) + &
1383 tij_abc(3, 1, k)*qp_j(3, 1) + &
1384 tij_abc(1, 2, k)*qp_j(1, 2) + &
1385 tij_abc(2, 2, k)*qp_j(2, 2) + &
1386 tij_abc(3, 2, k)*qp_j(3, 2) + &
1387 tij_abc(1, 3, k)*qp_j(1, 3) + &
1388 tij_abc(2, 3, k)*qp_j(2, 3) + &
1389 tij_abc(3, 3, k)*qp_j(3, 3))
1392 tmp_ji = ch_j*(tij_abc(1, 1, k)*qp_i(1, 1) + &
1393 tij_abc(2, 1, k)*qp_i(2, 1) + &
1394 tij_abc(3, 1, k)*qp_i(3, 1) + &
1395 tij_abc(1, 2, k)*qp_i(1, 2) + &
1396 tij_abc(2, 2, k)*qp_i(2, 2) + &
1397 tij_abc(3, 2, k)*qp_i(3, 2) + &
1398 tij_abc(1, 3, k)*qp_i(1, 3) + &
1399 tij_abc(2, 3, k)*qp_i(2, 3) + &
1400 tij_abc(3, 3, k)*qp_i(3, 3))
1402 fr(k) = fr(k) -
fac*(tmp_ij + tmp_ji)
1406 energy = energy + eloc
1408 forces(1, atom_a) = forces(1, atom_a) - fr(1)
1409 forces(2, atom_a) = forces(2, atom_a) - fr(2)
1410 forces(3, atom_a) = forces(3, atom_a) - fr(3)
1411 forces(1, atom_b) = forces(1, atom_b) + fr(1)
1412 forces(2, atom_b) = forces(2, atom_b) + fr(2)
1413 forces(3, atom_b) = forces(3, atom_b) + fr(3)
1418 IF (do_efield0)
THEN
1419 efield0(atom_a) = efield0(atom_a) + ef0_j
1421 efield0(atom_b) = efield0(atom_b) + ef0_i
1424 IF (do_efield1)
THEN
1425 efield1(1, atom_a) = efield1(1, atom_a) + ef1_j(1)
1426 efield1(2, atom_a) = efield1(2, atom_a) + ef1_j(2)
1427 efield1(3, atom_a) = efield1(3, atom_a) + ef1_j(3)
1429 efield1(1, atom_b) = efield1(1, atom_b) + ef1_i(1)
1430 efield1(2, atom_b) = efield1(2, atom_b) + ef1_i(2)
1431 efield1(3, atom_b) = efield1(3, atom_b) + ef1_i(3)
1434 IF (do_efield2)
THEN
1435 efield2(1, atom_a) = efield2(1, atom_a) + ef2_j(1, 1)
1436 efield2(2, atom_a) = efield2(2, atom_a) + ef2_j(1, 2)
1437 efield2(3, atom_a) = efield2(3, atom_a) + ef2_j(1, 3)
1438 efield2(4, atom_a) = efield2(4, atom_a) + ef2_j(2, 1)
1439 efield2(5, atom_a) = efield2(5, atom_a) + ef2_j(2, 2)
1440 efield2(6, atom_a) = efield2(6, atom_a) + ef2_j(2, 3)
1441 efield2(7, atom_a) = efield2(7, atom_a) + ef2_j(3, 1)
1442 efield2(8, atom_a) = efield2(8, atom_a) + ef2_j(3, 2)
1443 efield2(9, atom_a) = efield2(9, atom_a) + ef2_j(3, 3)
1445 efield2(1, atom_b) = efield2(1, atom_b) + ef2_i(1, 1)
1446 efield2(2, atom_b) = efield2(2, atom_b) + ef2_i(1, 2)
1447 efield2(3, atom_b) = efield2(3, atom_b) + ef2_i(1, 3)
1448 efield2(4, atom_b) = efield2(4, atom_b) + ef2_i(2, 1)
1449 efield2(5, atom_b) = efield2(5, atom_b) + ef2_i(2, 2)
1450 efield2(6, atom_b) = efield2(6, atom_b) + ef2_i(2, 3)
1451 efield2(7, atom_b) = efield2(7, atom_b) + ef2_i(3, 1)
1452 efield2(8, atom_b) = efield2(8, atom_b) + ef2_i(3, 2)
1453 efield2(9, atom_b) = efield2(9, atom_b) + ef2_i(3, 3)
1457 ptens11 = ptens11 + rab(1)*fr(1)
1458 ptens21 = ptens21 + rab(2)*fr(1)
1459 ptens31 = ptens31 + rab(3)*fr(1)
1460 ptens12 = ptens12 + rab(1)*fr(2)
1461 ptens22 = ptens22 + rab(2)*fr(2)
1462 ptens32 = ptens32 + rab(3)*fr(2)
1463 ptens13 = ptens13 + rab(1)*fr(3)
1464 ptens23 = ptens23 + rab(2)*fr(3)
1465 ptens33 = ptens33 + rab(3)*fr(3)
1472 tij_a = huge(0.0_dp)
1473 tij_ab = huge(0.0_dp)
1474 tij_abc = huge(0.0_dp)
1475 tij_abcd = huge(0.0_dp)
1476 tij_abcde = huge(0.0_dp)
1483 IF (debug_this_module .AND. debug_r_space .AND. (.NOT. debug_g_space))
THEN
1487 f(0) = erfc(alpha*r)*ir
1493 f(i) = irab2*(f(i - 1) + tmp*((2.0_dp*alpha**2)**i)/(
fac*alpha))
1498 force_eval = do_stress
1499 IF (task(1, 1))
THEN
1501 force_eval = do_forces .OR. do_efield1
1503 IF (task(2, 2)) force_eval = force_eval .OR. do_efield0
1504 IF (task(1, 2) .OR. force_eval)
THEN
1505 force_eval = do_stress
1506 tij_a = -rab*f(1)*fac_ij
1507 IF (task(1, 2)) force_eval = force_eval .OR. do_forces
1509 IF (task(1, 1)) force_eval = force_eval .OR. do_efield2
1510 IF (task(3, 3)) force_eval = force_eval .OR. do_efield0
1511 IF (task(2, 2) .OR. task(3, 1) .OR. force_eval)
THEN
1512 force_eval = do_stress
1515 tmp = rab(a)*rab(b)*fac_ij
1516 tij_ab(a, b) = 3.0_dp*tmp*f(2)
1517 IF (a == b) tij_ab(a, b) = tij_ab(a, b) - f(1)*fac_ij
1520 IF (task(2, 2) .OR. task(3, 1)) force_eval = force_eval .OR. do_forces
1522 IF (task(2, 2)) force_eval = force_eval .OR. do_efield2
1523 IF (task(3, 3)) force_eval = force_eval .OR. do_efield1
1524 IF (task(3, 2) .OR. force_eval)
THEN
1525 force_eval = do_stress
1529 tmp = rab(a)*rab(b)*rab(c)*fac_ij
1530 tij_abc(a, b, c) = -15.0_dp*tmp*f(3)
1531 tmp = 3.0_dp*f(2)*fac_ij
1532 IF (a == b) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(c)
1533 IF (a == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(b)
1534 IF (b == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(a)
1538 IF (task(3, 2)) force_eval = force_eval .OR. do_forces
1540 IF (task(3, 3) .OR. force_eval)
THEN
1541 force_eval = do_stress
1546 tmp = rab(a)*rab(b)*rab(c)*rab(d)*fac_ij
1547 tij_abcd(a, b, c, d) = 105.0_dp*tmp*f(4)
1548 tmp1 = 15.0_dp*f(3)*fac_ij
1549 tmp2 = 3.0_dp*f(2)*fac_ij
1551 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(c)*rab(d)
1552 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
1555 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(d)
1556 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
1558 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(c)
1560 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(d)
1561 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
1563 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(c)
1564 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(b)
1569 IF (task(3, 3)) force_eval = force_eval .OR. do_forces
1571 IF (force_eval)
THEN
1572 force_eval = do_stress
1578 tmp = rab(a)*rab(b)*rab(c)*rab(d)*rab(e)*fac_ij
1579 tij_abcde(a, b, c, d, e) = -945.0_dp*tmp*f(5)
1580 tmp1 = 105.0_dp*f(4)*fac_ij
1581 tmp2 = 15.0_dp*f(3)*fac_ij
1583 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(c)*rab(d)*rab(e)
1584 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
1585 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
1586 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
1589 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(d)*rab(e)
1590 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
1591 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
1592 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
1595 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(e)
1596 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
1597 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
1598 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
1601 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(d)
1602 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
1603 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
1604 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
1607 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(d)*rab(e)
1608 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
1611 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(e)
1612 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
1615 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(d)
1616 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
1618 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(e)
1619 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(d)
1620 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(c)
1638 IF (kind_a == ikind)
THEN
1640 SELECT CASE (itype_ij)
1645 DO kk = 1, nkdamp_ij
1647 factorial = factorial*real(kk, kind=
dp)
1648 dampsumfi = dampsumfi + (xf/factorial)
1650 dampaexpi = exp(-dampa_ij*r)
1651 dampfunci = dampsumfi*dampaexpi*dampfac_ij
1652 dampfuncdiffi = -dampa_ij*dampaexpi* &
1653 dampfac_ij*(((dampa_ij*r)**nkdamp_ij)/ &
1657 dampfuncdiffi = 0.0_dp
1661 SELECT CASE (itype_ji)
1666 DO kk = 1, nkdamp_ji
1668 factorial = factorial*real(kk, kind=
dp)
1669 dampsumfj = dampsumfj + (xf/factorial)
1671 dampaexpj = exp(-dampa_ji*r)
1672 dampfuncj = dampsumfj*dampaexpj*dampfac_ji
1673 dampfuncdiffj = -dampa_ji*dampaexpj* &
1674 dampfac_ji*(((dampa_ji*r)**nkdamp_ji)/ &
1678 dampfuncdiffj = 0.0_dp
1681 SELECT CASE (itype_ij)
1686 DO kk = 1, nkdamp_ij
1688 factorial = factorial*real(kk, kind=
dp)
1689 dampsumfj = dampsumfj + (xf/factorial)
1691 dampaexpj = exp(-dampa_ij*r)
1692 dampfuncj = dampsumfj*dampaexpj*dampfac_ij
1693 dampfuncdiffj = -dampa_ij*dampaexpj* &
1694 dampfac_ij*(((dampa_ij*r)**nkdamp_ij)/ &
1698 dampfuncdiffj = 0.0_dp
1702 SELECT CASE (itype_ji)
1707 DO kk = 1, nkdamp_ji
1709 factorial = factorial*real(kk, kind=
dp)
1710 dampsumfi = dampsumfi + (xf/factorial)
1712 dampaexpi = exp(-dampa_ji*r)
1713 dampfunci = dampsumfi*dampaexpi*dampfac_ji
1714 dampfuncdiffi = -dampa_ji*dampaexpi* &
1715 dampfac_ji*(((dampa_ji*r)**nkdamp_ji)/ &
1719 dampfuncdiffi = 0.0_dp
1723 damptij_a = -rab*dampfunci*fac_ij*irab2*ir
1724 damptji_a = -rab*dampfuncj*fac_ij*irab2*ir
1727 tmp = rab(a)*rab(b)*fac_ij
1728 damptij_ab(a, b) = tmp*(-dampfuncdiffi*irab2*irab2 + 3.0_dp*dampfunci*irab2*irab2*ir)
1729 damptji_ab(a, b) = tmp*(-dampfuncdiffj*irab2*irab2 + 3.0_dp*dampfuncj*irab2*irab2*ir)
1730 IF (a == b) damptij_ab(a, b) = damptij_ab(a, b) - dampfunci*fac_ij*irab2*ir
1731 IF (a == b) damptji_ab(a, b) = damptji_ab(a, b) - dampfuncj*fac_ij*irab2*ir
1737 IF (debug_this_module)
THEN
1745 IF (any(task(1, :)))
THEN
1746 ch_j = charges(atom_a)
1747 ch_i = charges(atom_b)
1749 IF (any(task(2, :)))
THEN
1750 dp_j = dipoles(:, atom_a)
1751 dp_i = dipoles(:, atom_b)
1753 IF (any(task(3, :)))
THEN
1754 qp_j = quadrupoles(:, :, atom_a)
1755 qp_i = quadrupoles(:, :, atom_b)
1757 IF (task(1, 1))
THEN
1759 eloc = eloc + ch_i*tij*ch_j
1761 IF (do_forces .OR. do_stress)
THEN
1762 fr(1) = fr(1) - ch_j*tij_a(1)*ch_i
1763 fr(2) = fr(2) - ch_j*tij_a(2)*ch_i
1764 fr(3) = fr(3) - ch_j*tij_a(3)*ch_i
1769 IF (do_efield0)
THEN
1770 ef0_i = ef0_i + tij*ch_j
1772 ef0_j = ef0_j + tij*ch_i
1775 IF (do_efield1)
THEN
1776 ef1_i(1) = ef1_i(1) - tij_a(1)*ch_j
1777 ef1_i(2) = ef1_i(2) - tij_a(2)*ch_j
1778 ef1_i(3) = ef1_i(3) - tij_a(3)*ch_j
1780 ef1_j(1) = ef1_j(1) + tij_a(1)*ch_i
1781 ef1_j(2) = ef1_j(2) + tij_a(2)*ch_i
1782 ef1_j(3) = ef1_j(3) + tij_a(3)*ch_i
1784 ef1_i(1) = ef1_i(1) + damptij_a(1)*ch_j
1785 ef1_i(2) = ef1_i(2) + damptij_a(2)*ch_j
1786 ef1_i(3) = ef1_i(3) + damptij_a(3)*ch_j
1788 ef1_j(1) = ef1_j(1) - damptji_a(1)*ch_i
1789 ef1_j(2) = ef1_j(2) - damptji_a(2)*ch_i
1790 ef1_j(3) = ef1_j(3) - damptji_a(3)*ch_i
1794 IF (do_efield2)
THEN
1795 ef2_i(1, 1) = ef2_i(1, 1) - tij_ab(1, 1)*ch_j
1796 ef2_i(2, 1) = ef2_i(2, 1) - tij_ab(2, 1)*ch_j
1797 ef2_i(3, 1) = ef2_i(3, 1) - tij_ab(3, 1)*ch_j
1798 ef2_i(1, 2) = ef2_i(1, 2) - tij_ab(1, 2)*ch_j
1799 ef2_i(2, 2) = ef2_i(2, 2) - tij_ab(2, 2)*ch_j
1800 ef2_i(3, 2) = ef2_i(3, 2) - tij_ab(3, 2)*ch_j
1801 ef2_i(1, 3) = ef2_i(1, 3) - tij_ab(1, 3)*ch_j
1802 ef2_i(2, 3) = ef2_i(2, 3) - tij_ab(2, 3)*ch_j
1803 ef2_i(3, 3) = ef2_i(3, 3) - tij_ab(3, 3)*ch_j
1805 ef2_j(1, 1) = ef2_j(1, 1) - tij_ab(1, 1)*ch_i
1806 ef2_j(2, 1) = ef2_j(2, 1) - tij_ab(2, 1)*ch_i
1807 ef2_j(3, 1) = ef2_j(3, 1) - tij_ab(3, 1)*ch_i
1808 ef2_j(1, 2) = ef2_j(1, 2) - tij_ab(1, 2)*ch_i
1809 ef2_j(2, 2) = ef2_j(2, 2) - tij_ab(2, 2)*ch_i
1810 ef2_j(3, 2) = ef2_j(3, 2) - tij_ab(3, 2)*ch_i
1811 ef2_j(1, 3) = ef2_j(1, 3) - tij_ab(1, 3)*ch_i
1812 ef2_j(2, 3) = ef2_j(2, 3) - tij_ab(2, 3)*ch_i
1813 ef2_j(3, 3) = ef2_j(3, 3) - tij_ab(3, 3)*ch_i
1817 IF (task(2, 2))
THEN
1819 tmp = -(dp_i(1)*(tij_ab(1, 1)*dp_j(1) + &
1820 tij_ab(2, 1)*dp_j(2) + &
1821 tij_ab(3, 1)*dp_j(3)) + &
1822 dp_i(2)*(tij_ab(1, 2)*dp_j(1) + &
1823 tij_ab(2, 2)*dp_j(2) + &
1824 tij_ab(3, 2)*dp_j(3)) + &
1825 dp_i(3)*(tij_ab(1, 3)*dp_j(1) + &
1826 tij_ab(2, 3)*dp_j(2) + &
1827 tij_ab(3, 3)*dp_j(3)))
1830 IF (do_forces .OR. do_stress)
THEN
1832 fr(k) = fr(k) + dp_i(1)*(tij_abc(1, 1, k)*dp_j(1) + &
1833 tij_abc(2, 1, k)*dp_j(2) + &
1834 tij_abc(3, 1, k)*dp_j(3)) &
1835 + dp_i(2)*(tij_abc(1, 2, k)*dp_j(1) + &
1836 tij_abc(2, 2, k)*dp_j(2) + &
1837 tij_abc(3, 2, k)*dp_j(3)) &
1838 + dp_i(3)*(tij_abc(1, 3, k)*dp_j(1) + &
1839 tij_abc(2, 3, k)*dp_j(2) + &
1840 tij_abc(3, 3, k)*dp_j(3))
1846 IF (do_efield0)
THEN
1847 ef0_i = ef0_i - (tij_a(1)*dp_j(1) + &
1848 tij_a(2)*dp_j(2) + &
1851 ef0_j = ef0_j + (tij_a(1)*dp_i(1) + &
1852 tij_a(2)*dp_i(2) + &
1856 IF (do_efield1)
THEN
1857 ef1_i(1) = ef1_i(1) + (tij_ab(1, 1)*dp_j(1) + &
1858 tij_ab(2, 1)*dp_j(2) + &
1859 tij_ab(3, 1)*dp_j(3))
1860 ef1_i(2) = ef1_i(2) + (tij_ab(1, 2)*dp_j(1) + &
1861 tij_ab(2, 2)*dp_j(2) + &
1862 tij_ab(3, 2)*dp_j(3))
1863 ef1_i(3) = ef1_i(3) + (tij_ab(1, 3)*dp_j(1) + &
1864 tij_ab(2, 3)*dp_j(2) + &
1865 tij_ab(3, 3)*dp_j(3))
1867 ef1_j(1) = ef1_j(1) + (tij_ab(1, 1)*dp_i(1) + &
1868 tij_ab(2, 1)*dp_i(2) + &
1869 tij_ab(3, 1)*dp_i(3))
1870 ef1_j(2) = ef1_j(2) + (tij_ab(1, 2)*dp_i(1) + &
1871 tij_ab(2, 2)*dp_i(2) + &
1872 tij_ab(3, 2)*dp_i(3))
1873 ef1_j(3) = ef1_j(3) + (tij_ab(1, 3)*dp_i(1) + &
1874 tij_ab(2, 3)*dp_i(2) + &
1875 tij_ab(3, 3)*dp_i(3))
1878 IF (do_efield2)
THEN
1879 ef2_i(1, 1) = ef2_i(1, 1) + (tij_abc(1, 1, 1)*dp_j(1) + &
1880 tij_abc(2, 1, 1)*dp_j(2) + &
1881 tij_abc(3, 1, 1)*dp_j(3))
1882 ef2_i(1, 2) = ef2_i(1, 2) + (tij_abc(1, 1, 2)*dp_j(1) + &
1883 tij_abc(2, 1, 2)*dp_j(2) + &
1884 tij_abc(3, 1, 2)*dp_j(3))
1885 ef2_i(1, 3) = ef2_i(1, 3) + (tij_abc(1, 1, 3)*dp_j(1) + &
1886 tij_abc(2, 1, 3)*dp_j(2) + &
1887 tij_abc(3, 1, 3)*dp_j(3))
1888 ef2_i(2, 1) = ef2_i(2, 1) + (tij_abc(1, 2, 1)*dp_j(1) + &
1889 tij_abc(2, 2, 1)*dp_j(2) + &
1890 tij_abc(3, 2, 1)*dp_j(3))
1891 ef2_i(2, 2) = ef2_i(2, 2) + (tij_abc(1, 2, 2)*dp_j(1) + &
1892 tij_abc(2, 2, 2)*dp_j(2) + &
1893 tij_abc(3, 2, 2)*dp_j(3))
1894 ef2_i(2, 3) = ef2_i(2, 3) + (tij_abc(1, 2, 3)*dp_j(1) + &
1895 tij_abc(2, 2, 3)*dp_j(2) + &
1896 tij_abc(3, 2, 3)*dp_j(3))
1897 ef2_i(3, 1) = ef2_i(3, 1) + (tij_abc(1, 3, 1)*dp_j(1) + &
1898 tij_abc(2, 3, 1)*dp_j(2) + &
1899 tij_abc(3, 3, 1)*dp_j(3))
1900 ef2_i(3, 2) = ef2_i(3, 2) + (tij_abc(1, 3, 2)*dp_j(1) + &
1901 tij_abc(2, 3, 2)*dp_j(2) + &
1902 tij_abc(3, 3, 2)*dp_j(3))
1903 ef2_i(3, 3) = ef2_i(3, 3) + (tij_abc(1, 3, 3)*dp_j(1) + &
1904 tij_abc(2, 3, 3)*dp_j(2) + &
1905 tij_abc(3, 3, 3)*dp_j(3))
1907 ef2_j(1, 1) = ef2_j(1, 1) - (tij_abc(1, 1, 1)*dp_i(1) + &
1908 tij_abc(2, 1, 1)*dp_i(2) + &
1909 tij_abc(3, 1, 1)*dp_i(3))
1910 ef2_j(1, 2) = ef2_j(1, 2) - (tij_abc(1, 1, 2)*dp_i(1) + &
1911 tij_abc(2, 1, 2)*dp_i(2) + &
1912 tij_abc(3, 1, 2)*dp_i(3))
1913 ef2_j(1, 3) = ef2_j(1, 3) - (tij_abc(1, 1, 3)*dp_i(1) + &
1914 tij_abc(2, 1, 3)*dp_i(2) + &
1915 tij_abc(3, 1, 3)*dp_i(3))
1916 ef2_j(2, 1) = ef2_j(2, 1) - (tij_abc(1, 2, 1)*dp_i(1) + &
1917 tij_abc(2, 2, 1)*dp_i(2) + &
1918 tij_abc(3, 2, 1)*dp_i(3))
1919 ef2_j(2, 2) = ef2_j(2, 2) - (tij_abc(1, 2, 2)*dp_i(1) + &
1920 tij_abc(2, 2, 2)*dp_i(2) + &
1921 tij_abc(3, 2, 2)*dp_i(3))
1922 ef2_j(2, 3) = ef2_j(2, 3) - (tij_abc(1, 2, 3)*dp_i(1) + &
1923 tij_abc(2, 2, 3)*dp_i(2) + &
1924 tij_abc(3, 2, 3)*dp_i(3))
1925 ef2_j(3, 1) = ef2_j(3, 1) - (tij_abc(1, 3, 1)*dp_i(1) + &
1926 tij_abc(2, 3, 1)*dp_i(2) + &
1927 tij_abc(3, 3, 1)*dp_i(3))
1928 ef2_j(3, 2) = ef2_j(3, 2) - (tij_abc(1, 3, 2)*dp_i(1) + &
1929 tij_abc(2, 3, 2)*dp_i(2) + &
1930 tij_abc(3, 3, 2)*dp_i(3))
1931 ef2_j(3, 3) = ef2_j(3, 3) - (tij_abc(1, 3, 3)*dp_i(1) + &
1932 tij_abc(2, 3, 3)*dp_i(2) + &
1933 tij_abc(3, 3, 3)*dp_i(3))
1937 IF (task(2, 1))
THEN
1939 tmp = ch_j*(tij_a(1)*dp_i(1) + &
1940 tij_a(2)*dp_i(2) + &
1942 - ch_i*(tij_a(1)*dp_j(1) + &
1943 tij_a(2)*dp_j(2) + &
1945 tmp = tmp - ch_j*(damptij_a(1)*dp_i(1) + &
1946 damptij_a(2)*dp_i(2) + &
1947 damptij_a(3)*dp_i(3)) &
1948 + ch_i*(damptji_a(1)*dp_j(1) + &
1949 damptji_a(2)*dp_j(2) + &
1950 damptji_a(3)*dp_j(3))
1953 IF (do_forces .OR. do_stress)
THEN
1955 fr(k) = fr(k) - ch_j*(tij_ab(1, k)*dp_i(1) + &
1956 tij_ab(2, k)*dp_i(2) + &
1957 tij_ab(3, k)*dp_i(3)) &
1958 + ch_i*(tij_ab(1, k)*dp_j(1) + &
1959 tij_ab(2, k)*dp_j(2) + &
1960 tij_ab(3, k)*dp_j(3))
1961 fr(k) = fr(k) + ch_j*(damptij_ab(1, k)*dp_i(1) + &
1962 damptij_ab(2, k)*dp_i(2) + &
1963 damptij_ab(3, k)*dp_i(3)) &
1964 - ch_i*(damptji_ab(1, k)*dp_j(1) + &
1965 damptji_ab(2, k)*dp_j(2) + &
1966 damptji_ab(3, k)*dp_j(3))
1970 IF (task(3, 3))
THEN
1973 tmp11 = qp_i(1, 1)*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
1974 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
1975 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
1976 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
1977 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
1978 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
1979 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
1980 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
1981 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
1982 tmp21 = qp_i(2, 1)*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
1983 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
1984 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
1985 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
1986 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
1987 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
1988 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
1989 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
1990 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
1991 tmp31 = qp_i(3, 1)*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
1992 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
1993 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
1994 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
1995 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
1996 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
1997 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
1998 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
1999 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
2000 tmp22 = qp_i(2, 2)*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
2001 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
2002 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
2003 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
2004 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
2005 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
2006 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
2007 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
2008 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
2009 tmp32 = qp_i(3, 2)*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
2010 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
2011 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
2012 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
2013 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
2014 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
2015 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
2016 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
2017 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
2018 tmp33 = qp_i(3, 3)*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
2019 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
2020 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
2021 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
2022 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
2023 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
2024 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
2025 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
2026 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
2030 tmp = tmp11 + tmp12 + tmp13 + &
2031 tmp21 + tmp22 + tmp23 + &
2032 tmp31 + tmp32 + tmp33
2034 eloc = eloc +
fac*tmp
2036 IF (do_forces .OR. do_stress)
THEN
2038 tmp11 = qp_i(1, 1)*(tij_abcde(1, 1, 1, 1, k)*qp_j(1, 1) + &
2039 tij_abcde(2, 1, 1, 1, k)*qp_j(2, 1) + &
2040 tij_abcde(3, 1, 1, 1, k)*qp_j(3, 1) + &
2041 tij_abcde(1, 2, 1, 1, k)*qp_j(1, 2) + &
2042 tij_abcde(2, 2, 1, 1, k)*qp_j(2, 2) + &
2043 tij_abcde(3, 2, 1, 1, k)*qp_j(3, 2) + &
2044 tij_abcde(1, 3, 1, 1, k)*qp_j(1, 3) + &
2045 tij_abcde(2, 3, 1, 1, k)*qp_j(2, 3) + &
2046 tij_abcde(3, 3, 1, 1, k)*qp_j(3, 3))
2047 tmp21 = qp_i(2, 1)*(tij_abcde(1, 1, 2, 1, k)*qp_j(1, 1) + &
2048 tij_abcde(2, 1, 2, 1, k)*qp_j(2, 1) + &
2049 tij_abcde(3, 1, 2, 1, k)*qp_j(3, 1) + &
2050 tij_abcde(1, 2, 2, 1, k)*qp_j(1, 2) + &
2051 tij_abcde(2, 2, 2, 1, k)*qp_j(2, 2) + &
2052 tij_abcde(3, 2, 2, 1, k)*qp_j(3, 2) + &
2053 tij_abcde(1, 3, 2, 1, k)*qp_j(1, 3) + &
2054 tij_abcde(2, 3, 2, 1, k)*qp_j(2, 3) + &
2055 tij_abcde(3, 3, 2, 1, k)*qp_j(3, 3))
2056 tmp31 = qp_i(3, 1)*(tij_abcde(1, 1, 3, 1, k)*qp_j(1, 1) + &
2057 tij_abcde(2, 1, 3, 1, k)*qp_j(2, 1) + &
2058 tij_abcde(3, 1, 3, 1, k)*qp_j(3, 1) + &
2059 tij_abcde(1, 2, 3, 1, k)*qp_j(1, 2) + &
2060 tij_abcde(2, 2, 3, 1, k)*qp_j(2, 2) + &
2061 tij_abcde(3, 2, 3, 1, k)*qp_j(3, 2) + &
2062 tij_abcde(1, 3, 3, 1, k)*qp_j(1, 3) + &
2063 tij_abcde(2, 3, 3, 1, k)*qp_j(2, 3) + &
2064 tij_abcde(3, 3, 3, 1, k)*qp_j(3, 3))
2065 tmp22 = qp_i(2, 2)*(tij_abcde(1, 1, 2, 2, k)*qp_j(1, 1) + &
2066 tij_abcde(2, 1, 2, 2, k)*qp_j(2, 1) + &
2067 tij_abcde(3, 1, 2, 2, k)*qp_j(3, 1) + &
2068 tij_abcde(1, 2, 2, 2, k)*qp_j(1, 2) + &
2069 tij_abcde(2, 2, 2, 2, k)*qp_j(2, 2) + &
2070 tij_abcde(3, 2, 2, 2, k)*qp_j(3, 2) + &
2071 tij_abcde(1, 3, 2, 2, k)*qp_j(1, 3) + &
2072 tij_abcde(2, 3, 2, 2, k)*qp_j(2, 3) + &
2073 tij_abcde(3, 3, 2, 2, k)*qp_j(3, 3))
2074 tmp32 = qp_i(3, 2)*(tij_abcde(1, 1, 3, 2, k)*qp_j(1, 1) + &
2075 tij_abcde(2, 1, 3, 2, k)*qp_j(2, 1) + &
2076 tij_abcde(3, 1, 3, 2, k)*qp_j(3, 1) + &
2077 tij_abcde(1, 2, 3, 2, k)*qp_j(1, 2) + &
2078 tij_abcde(2, 2, 3, 2, k)*qp_j(2, 2) + &
2079 tij_abcde(3, 2, 3, 2, k)*qp_j(3, 2) + &
2080 tij_abcde(1, 3, 3, 2, k)*qp_j(1, 3) + &
2081 tij_abcde(2, 3, 3, 2, k)*qp_j(2, 3) + &
2082 tij_abcde(3, 3, 3, 2, k)*qp_j(3, 3))
2083 tmp33 = qp_i(3, 3)*(tij_abcde(1, 1, 3, 3, k)*qp_j(1, 1) + &
2084 tij_abcde(2, 1, 3, 3, k)*qp_j(2, 1) + &
2085 tij_abcde(3, 1, 3, 3, k)*qp_j(3, 1) + &
2086 tij_abcde(1, 2, 3, 3, k)*qp_j(1, 2) + &
2087 tij_abcde(2, 2, 3, 3, k)*qp_j(2, 2) + &
2088 tij_abcde(3, 2, 3, 3, k)*qp_j(3, 2) + &
2089 tij_abcde(1, 3, 3, 3, k)*qp_j(1, 3) + &
2090 tij_abcde(2, 3, 3, 3, k)*qp_j(2, 3) + &
2091 tij_abcde(3, 3, 3, 3, k)*qp_j(3, 3))
2095 fr(k) = fr(k) -
fac*(tmp11 + tmp12 + tmp13 + &
2096 tmp21 + tmp22 + tmp23 + &
2097 tmp31 + tmp32 + tmp33)
2104 IF (do_efield0)
THEN
2105 ef0_i = ef0_i +
fac*(tij_ab(1, 1)*qp_j(1, 1) + &
2106 tij_ab(2, 1)*qp_j(2, 1) + &
2107 tij_ab(3, 1)*qp_j(3, 1) + &
2108 tij_ab(1, 2)*qp_j(1, 2) + &
2109 tij_ab(2, 2)*qp_j(2, 2) + &
2110 tij_ab(3, 2)*qp_j(3, 2) + &
2111 tij_ab(1, 3)*qp_j(1, 3) + &
2112 tij_ab(2, 3)*qp_j(2, 3) + &
2113 tij_ab(3, 3)*qp_j(3, 3))
2115 ef0_j = ef0_j +
fac*(tij_ab(1, 1)*qp_i(1, 1) + &
2116 tij_ab(2, 1)*qp_i(2, 1) + &
2117 tij_ab(3, 1)*qp_i(3, 1) + &
2118 tij_ab(1, 2)*qp_i(1, 2) + &
2119 tij_ab(2, 2)*qp_i(2, 2) + &
2120 tij_ab(3, 2)*qp_i(3, 2) + &
2121 tij_ab(1, 3)*qp_i(1, 3) + &
2122 tij_ab(2, 3)*qp_i(2, 3) + &
2123 tij_ab(3, 3)*qp_i(3, 3))
2126 IF (do_efield1)
THEN
2127 ef1_i(1) = ef1_i(1) -
fac*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
2128 tij_abc(2, 1, 1)*qp_j(2, 1) + &
2129 tij_abc(3, 1, 1)*qp_j(3, 1) + &
2130 tij_abc(1, 2, 1)*qp_j(1, 2) + &
2131 tij_abc(2, 2, 1)*qp_j(2, 2) + &
2132 tij_abc(3, 2, 1)*qp_j(3, 2) + &
2133 tij_abc(1, 3, 1)*qp_j(1, 3) + &
2134 tij_abc(2, 3, 1)*qp_j(2, 3) + &
2135 tij_abc(3, 3, 1)*qp_j(3, 3))
2136 ef1_i(2) = ef1_i(2) -
fac*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
2137 tij_abc(2, 1, 2)*qp_j(2, 1) + &
2138 tij_abc(3, 1, 2)*qp_j(3, 1) + &
2139 tij_abc(1, 2, 2)*qp_j(1, 2) + &
2140 tij_abc(2, 2, 2)*qp_j(2, 2) + &
2141 tij_abc(3, 2, 2)*qp_j(3, 2) + &
2142 tij_abc(1, 3, 2)*qp_j(1, 3) + &
2143 tij_abc(2, 3, 2)*qp_j(2, 3) + &
2144 tij_abc(3, 3, 2)*qp_j(3, 3))
2145 ef1_i(3) = ef1_i(3) -
fac*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
2146 tij_abc(2, 1, 3)*qp_j(2, 1) + &
2147 tij_abc(3, 1, 3)*qp_j(3, 1) + &
2148 tij_abc(1, 2, 3)*qp_j(1, 2) + &
2149 tij_abc(2, 2, 3)*qp_j(2, 2) + &
2150 tij_abc(3, 2, 3)*qp_j(3, 2) + &
2151 tij_abc(1, 3, 3)*qp_j(1, 3) + &
2152 tij_abc(2, 3, 3)*qp_j(2, 3) + &
2153 tij_abc(3, 3, 3)*qp_j(3, 3))
2155 ef1_j(1) = ef1_j(1) +
fac*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
2156 tij_abc(2, 1, 1)*qp_i(2, 1) + &
2157 tij_abc(3, 1, 1)*qp_i(3, 1) + &
2158 tij_abc(1, 2, 1)*qp_i(1, 2) + &
2159 tij_abc(2, 2, 1)*qp_i(2, 2) + &
2160 tij_abc(3, 2, 1)*qp_i(3, 2) + &
2161 tij_abc(1, 3, 1)*qp_i(1, 3) + &
2162 tij_abc(2, 3, 1)*qp_i(2, 3) + &
2163 tij_abc(3, 3, 1)*qp_i(3, 3))
2164 ef1_j(2) = ef1_j(2) +
fac*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
2165 tij_abc(2, 1, 2)*qp_i(2, 1) + &
2166 tij_abc(3, 1, 2)*qp_i(3, 1) + &
2167 tij_abc(1, 2, 2)*qp_i(1, 2) + &
2168 tij_abc(2, 2, 2)*qp_i(2, 2) + &
2169 tij_abc(3, 2, 2)*qp_i(3, 2) + &
2170 tij_abc(1, 3, 2)*qp_i(1, 3) + &
2171 tij_abc(2, 3, 2)*qp_i(2, 3) + &
2172 tij_abc(3, 3, 2)*qp_i(3, 3))
2173 ef1_j(3) = ef1_j(3) +
fac*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
2174 tij_abc(2, 1, 3)*qp_i(2, 1) + &
2175 tij_abc(3, 1, 3)*qp_i(3, 1) + &
2176 tij_abc(1, 2, 3)*qp_i(1, 2) + &
2177 tij_abc(2, 2, 3)*qp_i(2, 2) + &
2178 tij_abc(3, 2, 3)*qp_i(3, 2) + &
2179 tij_abc(1, 3, 3)*qp_i(1, 3) + &
2180 tij_abc(2, 3, 3)*qp_i(2, 3) + &
2181 tij_abc(3, 3, 3)*qp_i(3, 3))
2184 IF (do_efield2)
THEN
2185 tmp11 =
fac*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
2186 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
2187 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
2188 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
2189 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
2190 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
2191 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
2192 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
2193 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
2194 tmp12 =
fac*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
2195 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
2196 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
2197 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
2198 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
2199 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
2200 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
2201 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
2202 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
2203 tmp13 =
fac*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
2204 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
2205 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
2206 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
2207 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
2208 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
2209 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
2210 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
2211 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
2212 tmp22 =
fac*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
2213 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
2214 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
2215 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
2216 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
2217 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
2218 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
2219 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
2220 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
2221 tmp23 =
fac*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
2222 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
2223 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
2224 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
2225 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
2226 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
2227 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
2228 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
2229 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
2230 tmp33 =
fac*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
2231 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
2232 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
2233 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
2234 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
2235 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
2236 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
2237 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
2238 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
2240 ef2_i(1, 1) = ef2_i(1, 1) - tmp11
2241 ef2_i(1, 2) = ef2_i(1, 2) - tmp12
2242 ef2_i(1, 3) = ef2_i(1, 3) - tmp13
2243 ef2_i(2, 1) = ef2_i(2, 1) - tmp12
2244 ef2_i(2, 2) = ef2_i(2, 2) - tmp22
2245 ef2_i(2, 3) = ef2_i(2, 3) - tmp23
2246 ef2_i(3, 1) = ef2_i(3, 1) - tmp13
2247 ef2_i(3, 2) = ef2_i(3, 2) - tmp23
2248 ef2_i(3, 3) = ef2_i(3, 3) - tmp33
2250 tmp11 =
fac*(tij_abcd(1, 1, 1, 1)*qp_i(1, 1) + &
2251 tij_abcd(2, 1, 1, 1)*qp_i(2, 1) + &
2252 tij_abcd(3, 1, 1, 1)*qp_i(3, 1) + &
2253 tij_abcd(1, 2, 1, 1)*qp_i(1, 2) + &
2254 tij_abcd(2, 2, 1, 1)*qp_i(2, 2) + &
2255 tij_abcd(3, 2, 1, 1)*qp_i(3, 2) + &
2256 tij_abcd(1, 3, 1, 1)*qp_i(1, 3) + &
2257 tij_abcd(2, 3, 1, 1)*qp_i(2, 3) + &
2258 tij_abcd(3, 3, 1, 1)*qp_i(3, 3))
2259 tmp12 =
fac*(tij_abcd(1, 1, 1, 2)*qp_i(1, 1) + &
2260 tij_abcd(2, 1, 1, 2)*qp_i(2, 1) + &
2261 tij_abcd(3, 1, 1, 2)*qp_i(3, 1) + &
2262 tij_abcd(1, 2, 1, 2)*qp_i(1, 2) + &
2263 tij_abcd(2, 2, 1, 2)*qp_i(2, 2) + &
2264 tij_abcd(3, 2, 1, 2)*qp_i(3, 2) + &
2265 tij_abcd(1, 3, 1, 2)*qp_i(1, 3) + &
2266 tij_abcd(2, 3, 1, 2)*qp_i(2, 3) + &
2267 tij_abcd(3, 3, 1, 2)*qp_i(3, 3))
2268 tmp13 =
fac*(tij_abcd(1, 1, 1, 3)*qp_i(1, 1) + &
2269 tij_abcd(2, 1, 1, 3)*qp_i(2, 1) + &
2270 tij_abcd(3, 1, 1, 3)*qp_i(3, 1) + &
2271 tij_abcd(1, 2, 1, 3)*qp_i(1, 2) + &
2272 tij_abcd(2, 2, 1, 3)*qp_i(2, 2) + &
2273 tij_abcd(3, 2, 1, 3)*qp_i(3, 2) + &
2274 tij_abcd(1, 3, 1, 3)*qp_i(1, 3) + &
2275 tij_abcd(2, 3, 1, 3)*qp_i(2, 3) + &
2276 tij_abcd(3, 3, 1, 3)*qp_i(3, 3))
2277 tmp22 =
fac*(tij_abcd(1, 1, 2, 2)*qp_i(1, 1) + &
2278 tij_abcd(2, 1, 2, 2)*qp_i(2, 1) + &
2279 tij_abcd(3, 1, 2, 2)*qp_i(3, 1) + &
2280 tij_abcd(1, 2, 2, 2)*qp_i(1, 2) + &
2281 tij_abcd(2, 2, 2, 2)*qp_i(2, 2) + &
2282 tij_abcd(3, 2, 2, 2)*qp_i(3, 2) + &
2283 tij_abcd(1, 3, 2, 2)*qp_i(1, 3) + &
2284 tij_abcd(2, 3, 2, 2)*qp_i(2, 3) + &
2285 tij_abcd(3, 3, 2, 2)*qp_i(3, 3))
2286 tmp23 =
fac*(tij_abcd(1, 1, 2, 3)*qp_i(1, 1) + &
2287 tij_abcd(2, 1, 2, 3)*qp_i(2, 1) + &
2288 tij_abcd(3, 1, 2, 3)*qp_i(3, 1) + &
2289 tij_abcd(1, 2, 2, 3)*qp_i(1, 2) + &
2290 tij_abcd(2, 2, 2, 3)*qp_i(2, 2) + &
2291 tij_abcd(3, 2, 2, 3)*qp_i(3, 2) + &
2292 tij_abcd(1, 3, 2, 3)*qp_i(1, 3) + &
2293 tij_abcd(2, 3, 2, 3)*qp_i(2, 3) + &
2294 tij_abcd(3, 3, 2, 3)*qp_i(3, 3))
2295 tmp33 =
fac*(tij_abcd(1, 1, 3, 3)*qp_i(1, 1) + &
2296 tij_abcd(2, 1, 3, 3)*qp_i(2, 1) + &
2297 tij_abcd(3, 1, 3, 3)*qp_i(3, 1) + &
2298 tij_abcd(1, 2, 3, 3)*qp_i(1, 2) + &
2299 tij_abcd(2, 2, 3, 3)*qp_i(2, 2) + &
2300 tij_abcd(3, 2, 3, 3)*qp_i(3, 2) + &
2301 tij_abcd(1, 3, 3, 3)*qp_i(1, 3) + &
2302 tij_abcd(2, 3, 3, 3)*qp_i(2, 3) + &
2303 tij_abcd(3, 3, 3, 3)*qp_i(3, 3))
2305 ef2_j(1, 1) = ef2_j(1, 1) - tmp11
2306 ef2_j(1, 2) = ef2_j(1, 2) - tmp12
2307 ef2_j(1, 3) = ef2_j(1, 3) - tmp13
2308 ef2_j(2, 1) = ef2_j(2, 1) - tmp12
2309 ef2_j(2, 2) = ef2_j(2, 2) - tmp22
2310 ef2_j(2, 3) = ef2_j(2, 3) - tmp23
2311 ef2_j(3, 1) = ef2_j(3, 1) - tmp13
2312 ef2_j(3, 2) = ef2_j(3, 2) - tmp23
2313 ef2_j(3, 3) = ef2_j(3, 3) - tmp33
2317 IF (task(3, 2))
THEN
2321 tmp_ij = dp_i(1)*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
2322 tij_abc(2, 1, 1)*qp_j(2, 1) + &
2323 tij_abc(3, 1, 1)*qp_j(3, 1) + &
2324 tij_abc(1, 2, 1)*qp_j(1, 2) + &
2325 tij_abc(2, 2, 1)*qp_j(2, 2) + &
2326 tij_abc(3, 2, 1)*qp_j(3, 2) + &
2327 tij_abc(1, 3, 1)*qp_j(1, 3) + &
2328 tij_abc(2, 3, 1)*qp_j(2, 3) + &
2329 tij_abc(3, 3, 1)*qp_j(3, 3)) + &
2330 dp_i(2)*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
2331 tij_abc(2, 1, 2)*qp_j(2, 1) + &
2332 tij_abc(3, 1, 2)*qp_j(3, 1) + &
2333 tij_abc(1, 2, 2)*qp_j(1, 2) + &
2334 tij_abc(2, 2, 2)*qp_j(2, 2) + &
2335 tij_abc(3, 2, 2)*qp_j(3, 2) + &
2336 tij_abc(1, 3, 2)*qp_j(1, 3) + &
2337 tij_abc(2, 3, 2)*qp_j(2, 3) + &
2338 tij_abc(3, 3, 2)*qp_j(3, 3)) + &
2339 dp_i(3)*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
2340 tij_abc(2, 1, 3)*qp_j(2, 1) + &
2341 tij_abc(3, 1, 3)*qp_j(3, 1) + &
2342 tij_abc(1, 2, 3)*qp_j(1, 2) + &
2343 tij_abc(2, 2, 3)*qp_j(2, 2) + &
2344 tij_abc(3, 2, 3)*qp_j(3, 2) + &
2345 tij_abc(1, 3, 3)*qp_j(1, 3) + &
2346 tij_abc(2, 3, 3)*qp_j(2, 3) + &
2347 tij_abc(3, 3, 3)*qp_j(3, 3))
2350 tmp_ji = dp_j(1)*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
2351 tij_abc(2, 1, 1)*qp_i(2, 1) + &
2352 tij_abc(3, 1, 1)*qp_i(3, 1) + &
2353 tij_abc(1, 2, 1)*qp_i(1, 2) + &
2354 tij_abc(2, 2, 1)*qp_i(2, 2) + &
2355 tij_abc(3, 2, 1)*qp_i(3, 2) + &
2356 tij_abc(1, 3, 1)*qp_i(1, 3) + &
2357 tij_abc(2, 3, 1)*qp_i(2, 3) + &
2358 tij_abc(3, 3, 1)*qp_i(3, 3)) + &
2359 dp_j(2)*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
2360 tij_abc(2, 1, 2)*qp_i(2, 1) + &
2361 tij_abc(3, 1, 2)*qp_i(3, 1) + &
2362 tij_abc(1, 2, 2)*qp_i(1, 2) + &
2363 tij_abc(2, 2, 2)*qp_i(2, 2) + &
2364 tij_abc(3, 2, 2)*qp_i(3, 2) + &
2365 tij_abc(1, 3, 2)*qp_i(1, 3) + &
2366 tij_abc(2, 3, 2)*qp_i(2, 3) + &
2367 tij_abc(3, 3, 2)*qp_i(3, 3)) + &
2368 dp_j(3)*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
2369 tij_abc(2, 1, 3)*qp_i(2, 1) + &
2370 tij_abc(3, 1, 3)*qp_i(3, 1) + &
2371 tij_abc(1, 2, 3)*qp_i(1, 2) + &
2372 tij_abc(2, 2, 3)*qp_i(2, 2) + &
2373 tij_abc(3, 2, 3)*qp_i(3, 2) + &
2374 tij_abc(1, 3, 3)*qp_i(1, 3) + &
2375 tij_abc(2, 3, 3)*qp_i(2, 3) + &
2376 tij_abc(3, 3, 3)*qp_i(3, 3))
2378 tmp =
fac*(tmp_ij - tmp_ji)
2380 IF (do_forces .OR. do_stress)
THEN
2383 tmp_ij = dp_i(1)*(tij_abcd(1, 1, 1, k)*qp_j(1, 1) + &
2384 tij_abcd(2, 1, 1, k)*qp_j(2, 1) + &
2385 tij_abcd(3, 1, 1, k)*qp_j(3, 1) + &
2386 tij_abcd(1, 2, 1, k)*qp_j(1, 2) + &
2387 tij_abcd(2, 2, 1, k)*qp_j(2, 2) + &
2388 tij_abcd(3, 2, 1, k)*qp_j(3, 2) + &
2389 tij_abcd(1, 3, 1, k)*qp_j(1, 3) + &
2390 tij_abcd(2, 3, 1, k)*qp_j(2, 3) + &
2391 tij_abcd(3, 3, 1, k)*qp_j(3, 3)) + &
2392 dp_i(2)*(tij_abcd(1, 1, 2, k)*qp_j(1, 1) + &
2393 tij_abcd(2, 1, 2, k)*qp_j(2, 1) + &
2394 tij_abcd(3, 1, 2, k)*qp_j(3, 1) + &
2395 tij_abcd(1, 2, 2, k)*qp_j(1, 2) + &
2396 tij_abcd(2, 2, 2, k)*qp_j(2, 2) + &
2397 tij_abcd(3, 2, 2, k)*qp_j(3, 2) + &
2398 tij_abcd(1, 3, 2, k)*qp_j(1, 3) + &
2399 tij_abcd(2, 3, 2, k)*qp_j(2, 3) + &
2400 tij_abcd(3, 3, 2, k)*qp_j(3, 3)) + &
2401 dp_i(3)*(tij_abcd(1, 1, 3, k)*qp_j(1, 1) + &
2402 tij_abcd(2, 1, 3, k)*qp_j(2, 1) + &
2403 tij_abcd(3, 1, 3, k)*qp_j(3, 1) + &
2404 tij_abcd(1, 2, 3, k)*qp_j(1, 2) + &
2405 tij_abcd(2, 2, 3, k)*qp_j(2, 2) + &
2406 tij_abcd(3, 2, 3, k)*qp_j(3, 2) + &
2407 tij_abcd(1, 3, 3, k)*qp_j(1, 3) + &
2408 tij_abcd(2, 3, 3, k)*qp_j(2, 3) + &
2409 tij_abcd(3, 3, 3, k)*qp_j(3, 3))
2412 tmp_ji = dp_j(1)*(tij_abcd(1, 1, 1, k)*qp_i(1, 1) + &
2413 tij_abcd(2, 1, 1, k)*qp_i(2, 1) + &
2414 tij_abcd(3, 1, 1, k)*qp_i(3, 1) + &
2415 tij_abcd(1, 2, 1, k)*qp_i(1, 2) + &
2416 tij_abcd(2, 2, 1, k)*qp_i(2, 2) + &
2417 tij_abcd(3, 2, 1, k)*qp_i(3, 2) + &
2418 tij_abcd(1, 3, 1, k)*qp_i(1, 3) + &
2419 tij_abcd(2, 3, 1, k)*qp_i(2, 3) + &
2420 tij_abcd(3, 3, 1, k)*qp_i(3, 3)) + &
2421 dp_j(2)*(tij_abcd(1, 1, 2, k)*qp_i(1, 1) + &
2422 tij_abcd(2, 1, 2, k)*qp_i(2, 1) + &
2423 tij_abcd(3, 1, 2, k)*qp_i(3, 1) + &
2424 tij_abcd(1, 2, 2, k)*qp_i(1, 2) + &
2425 tij_abcd(2, 2, 2, k)*qp_i(2, 2) + &
2426 tij_abcd(3, 2, 2, k)*qp_i(3, 2) + &
2427 tij_abcd(1, 3, 2, k)*qp_i(1, 3) + &
2428 tij_abcd(2, 3, 2, k)*qp_i(2, 3) + &
2429 tij_abcd(3, 3, 2, k)*qp_i(3, 3)) + &
2430 dp_j(3)*(tij_abcd(1, 1, 3, k)*qp_i(1, 1) + &
2431 tij_abcd(2, 1, 3, k)*qp_i(2, 1) + &
2432 tij_abcd(3, 1, 3, k)*qp_i(3, 1) + &
2433 tij_abcd(1, 2, 3, k)*qp_i(1, 2) + &
2434 tij_abcd(2, 2, 3, k)*qp_i(2, 2) + &
2435 tij_abcd(3, 2, 3, k)*qp_i(3, 2) + &
2436 tij_abcd(1, 3, 3, k)*qp_i(1, 3) + &
2437 tij_abcd(2, 3, 3, k)*qp_i(2, 3) + &
2438 tij_abcd(3, 3, 3, k)*qp_i(3, 3))
2440 fr(k) = fr(k) -
fac*(tmp_ij - tmp_ji)
2444 IF (task(3, 1))
THEN
2449 tmp_ij = ch_i*(tij_ab(1, 1)*qp_j(1, 1) + &
2450 tij_ab(2, 1)*qp_j(2, 1) + &
2451 tij_ab(3, 1)*qp_j(3, 1) + &
2452 tij_ab(1, 2)*qp_j(1, 2) + &
2453 tij_ab(2, 2)*qp_j(2, 2) + &
2454 tij_ab(3, 2)*qp_j(3, 2) + &
2455 tij_ab(1, 3)*qp_j(1, 3) + &
2456 tij_ab(2, 3)*qp_j(2, 3) + &
2457 tij_ab(3, 3)*qp_j(3, 3))
2460 tmp_ji = ch_j*(tij_ab(1, 1)*qp_i(1, 1) + &
2461 tij_ab(2, 1)*qp_i(2, 1) + &
2462 tij_ab(3, 1)*qp_i(3, 1) + &
2463 tij_ab(1, 2)*qp_i(1, 2) + &
2464 tij_ab(2, 2)*qp_i(2, 2) + &
2465 tij_ab(3, 2)*qp_i(3, 2) + &
2466 tij_ab(1, 3)*qp_i(1, 3) + &
2467 tij_ab(2, 3)*qp_i(2, 3) + &
2468 tij_ab(3, 3)*qp_i(3, 3))
2470 eloc = eloc +
fac*(tmp_ij + tmp_ji)
2471 IF (do_forces .OR. do_stress)
THEN
2474 tmp_ij = ch_i*(tij_abc(1, 1, k)*qp_j(1, 1) + &
2475 tij_abc(2, 1, k)*qp_j(2, 1) + &
2476 tij_abc(3, 1, k)*qp_j(3, 1) + &
2477 tij_abc(1, 2, k)*qp_j(1, 2) + &
2478 tij_abc(2, 2, k)*qp_j(2, 2) + &
2479 tij_abc(3, 2, k)*qp_j(3, 2) + &
2480 tij_abc(1, 3, k)*qp_j(1, 3) + &
2481 tij_abc(2, 3, k)*qp_j(2, 3) + &
2482 tij_abc(3, 3, k)*qp_j(3, 3))
2485 tmp_ji = ch_j*(tij_abc(1, 1, k)*qp_i(1, 1) + &
2486 tij_abc(2, 1, k)*qp_i(2, 1) + &
2487 tij_abc(3, 1, k)*qp_i(3, 1) + &
2488 tij_abc(1, 2, k)*qp_i(1, 2) + &
2489 tij_abc(2, 2, k)*qp_i(2, 2) + &
2490 tij_abc(3, 2, k)*qp_i(3, 2) + &
2491 tij_abc(1, 3, k)*qp_i(1, 3) + &
2492 tij_abc(2, 3, k)*qp_i(2, 3) + &
2493 tij_abc(3, 3, k)*qp_i(3, 3))
2495 fr(k) = fr(k) -
fac*(tmp_ij + tmp_ji)
2499 energy = energy + eloc
2501 forces(1, atom_a) = forces(1, atom_a) - fr(1)
2502 forces(2, atom_a) = forces(2, atom_a) - fr(2)
2503 forces(3, atom_a) = forces(3, atom_a) - fr(3)
2504 forces(1, atom_b) = forces(1, atom_b) + fr(1)
2505 forces(2, atom_b) = forces(2, atom_b) + fr(2)
2506 forces(3, atom_b) = forces(3, atom_b) + fr(3)
2511 IF (do_efield0)
THEN
2512 efield0(atom_a) = efield0(atom_a) + ef0_j
2514 efield0(atom_b) = efield0(atom_b) + ef0_i
2517 IF (do_efield1)
THEN
2518 efield1(1, atom_a) = efield1(1, atom_a) + ef1_j(1)
2519 efield1(2, atom_a) = efield1(2, atom_a) + ef1_j(2)
2520 efield1(3, atom_a) = efield1(3, atom_a) + ef1_j(3)
2522 efield1(1, atom_b) = efield1(1, atom_b) + ef1_i(1)
2523 efield1(2, atom_b) = efield1(2, atom_b) + ef1_i(2)
2524 efield1(3, atom_b) = efield1(3, atom_b) + ef1_i(3)
2527 IF (do_efield2)
THEN
2528 efield2(1, atom_a) = efield2(1, atom_a) + ef2_j(1, 1)
2529 efield2(2, atom_a) = efield2(2, atom_a) + ef2_j(1, 2)
2530 efield2(3, atom_a) = efield2(3, atom_a) + ef2_j(1, 3)
2531 efield2(4, atom_a) = efield2(4, atom_a) + ef2_j(2, 1)
2532 efield2(5, atom_a) = efield2(5, atom_a) + ef2_j(2, 2)
2533 efield2(6, atom_a) = efield2(6, atom_a) + ef2_j(2, 3)
2534 efield2(7, atom_a) = efield2(7, atom_a) + ef2_j(3, 1)
2535 efield2(8, atom_a) = efield2(8, atom_a) + ef2_j(3, 2)
2536 efield2(9, atom_a) = efield2(9, atom_a) + ef2_j(3, 3)
2538 efield2(1, atom_b) = efield2(1, atom_b) + ef2_i(1, 1)
2539 efield2(2, atom_b) = efield2(2, atom_b) + ef2_i(1, 2)
2540 efield2(3, atom_b) = efield2(3, atom_b) + ef2_i(1, 3)
2541 efield2(4, atom_b) = efield2(4, atom_b) + ef2_i(2, 1)
2542 efield2(5, atom_b) = efield2(5, atom_b) + ef2_i(2, 2)
2543 efield2(6, atom_b) = efield2(6, atom_b) + ef2_i(2, 3)
2544 efield2(7, atom_b) = efield2(7, atom_b) + ef2_i(3, 1)
2545 efield2(8, atom_b) = efield2(8, atom_b) + ef2_i(3, 2)
2546 efield2(9, atom_b) = efield2(9, atom_b) + ef2_i(3, 3)
2550 ptens11 = ptens11 + rab(1)*fr(1)
2551 ptens21 = ptens21 + rab(2)*fr(1)
2552 ptens31 = ptens31 + rab(3)*fr(1)
2553 ptens12 = ptens12 + rab(1)*fr(2)
2554 ptens22 = ptens22 + rab(2)*fr(2)
2555 ptens32 = ptens32 + rab(3)*fr(2)
2556 ptens13 = ptens13 + rab(1)*fr(3)
2557 ptens23 = ptens23 + rab(2)*fr(3)
2558 ptens33 = ptens33 + rab(3)*fr(3)
2564 END DO kind_group_loop
2567 pv(1, 1) = pv(1, 1) + ptens11
2568 pv(1, 2) = pv(1, 2) + (ptens12 + ptens21)*0.5_dp
2569 pv(1, 3) = pv(1, 3) + (ptens13 + ptens31)*0.5_dp
2571 pv(2, 2) = pv(2, 2) + ptens22
2572 pv(2, 3) = pv(2, 3) + (ptens23 + ptens32)*0.5_dp
2575 pv(3, 3) = pv(3, 3) + ptens33
2578 CALL timestop(handle)
2579 END SUBROUTINE ewald_multipole_sr
2603 SUBROUTINE ewald_multipole_bonded(nonbond_env, particle_set, ewald_env, &
2604 cell, energy, task, do_forces, do_efield, do_stress, charges, &
2605 dipoles, quadrupoles, forces, pv, efield0, efield1, efield2)
2611 REAL(kind=
dp),
INTENT(INOUT) :: energy
2612 LOGICAL,
DIMENSION(3, 3),
INTENT(IN) :: task
2613 LOGICAL,
INTENT(IN) :: do_forces, do_efield, do_stress
2614 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: charges
2615 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
2616 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
2617 POINTER :: quadrupoles
2618 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT), &
2619 OPTIONAL :: forces, pv
2620 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0
2621 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: efield1, efield2
2623 CHARACTER(len=*),
PARAMETER :: routinen =
'ewald_multipole_bonded'
2625 INTEGER :: a, atom_a, atom_b, b, c, d, e, handle, &
2626 i, iend, igrp, ilist, ipair, istart, &
2628 INTEGER,
DIMENSION(:, :),
POINTER ::
list
2629 LOGICAL :: do_efield0, do_efield1, do_efield2, &
2631 REAL(kind=
dp) :: alpha, ch_i, ch_j, ef0_i, ef0_j, eloc,
fac, fac_ij, ir, irab2, ptens11, &
2632 ptens12, ptens13, ptens21, ptens22, ptens23, ptens31, ptens32, ptens33, r, rab2, tij, &
2633 tmp, tmp1, tmp11, tmp12, tmp13, tmp2, tmp21, tmp22, tmp23, tmp31, tmp32, tmp33, tmp_ij, &
2635 REAL(kind=
dp),
DIMENSION(0:5) :: f
2636 REAL(kind=
dp),
DIMENSION(3) :: dp_i, dp_j, ef1_i, ef1_j, fr, rab, tij_a
2637 REAL(kind=
dp),
DIMENSION(3, 3) :: ef2_i, ef2_j, qp_i, qp_j, tij_ab
2638 REAL(kind=
dp),
DIMENSION(3, 3, 3) :: tij_abc
2639 REAL(kind=
dp),
DIMENSION(3, 3, 3, 3) :: tij_abcd
2640 REAL(kind=
dp),
DIMENSION(3, 3, 3, 3, 3) :: tij_abcde
2644 CALL timeset(routinen, handle)
2645 do_efield0 = do_efield .AND.
ASSOCIATED(efield0)
2646 do_efield1 = do_efield .AND.
ASSOCIATED(efield1)
2647 do_efield2 = do_efield .AND.
ASSOCIATED(efield2)
2649 ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp
2650 ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp
2651 ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp
2657 lists:
DO ilist = 1, nonbonded%nlists
2658 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
2659 nscale = neighbor_kind_pair%nscale
2660 IF (nscale == 0) cycle lists
2661 list => neighbor_kind_pair%list
2662 kind_group_loop:
DO igrp = 1, neighbor_kind_pair%ngrp_kind
2663 istart = neighbor_kind_pair%grp_kind_start(igrp)
2664 IF (istart > nscale) cycle kind_group_loop
2665 iend = min(neighbor_kind_pair%grp_kind_end(igrp), nscale)
2666 pairs:
DO ipair = istart, iend
2668 fac_ij = -1.0_dp + neighbor_kind_pair%ei_scale(ipair)
2669 IF (fac_ij >= 0) cycle pairs
2671 atom_a =
list(1, ipair)
2672 atom_b =
list(2, ipair)
2674 rab = particle_set(atom_b)%r - particle_set(atom_a)%r
2675 rab =
pbc(rab, cell)
2676 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
2680 tij_a = huge(0.0_dp)
2681 tij_ab = huge(0.0_dp)
2682 tij_abc = huge(0.0_dp)
2683 tij_abcd = huge(0.0_dp)
2684 tij_abcde = huge(0.0_dp)
2690 IF (debug_this_module .AND. debug_r_space .AND. (.NOT. debug_g_space))
THEN
2694 f(0) = erf(alpha*r)*ir
2700 f(i) = irab2*(f(i - 1) - tmp*((2.0_dp*alpha**2)**i)/(
fac*alpha))
2705 force_eval = do_stress
2706 IF (task(1, 1))
THEN
2708 force_eval = do_forces .OR. do_efield1
2710 IF (task(2, 2)) force_eval = force_eval .OR. do_efield0
2711 IF (task(1, 2) .OR. force_eval)
THEN
2712 force_eval = do_stress
2713 tij_a = -rab*f(1)*fac_ij
2714 IF (task(1, 2)) force_eval = force_eval .OR. do_forces
2716 IF (task(1, 1)) force_eval = force_eval .OR. do_efield2
2717 IF (task(3, 3)) force_eval = force_eval .OR. do_efield0
2718 IF (task(2, 2) .OR. task(3, 1) .OR. force_eval)
THEN
2719 force_eval = do_stress
2722 tmp = rab(a)*rab(b)*fac_ij
2723 tij_ab(a, b) = 3.0_dp*tmp*f(2)
2724 IF (a == b) tij_ab(a, b) = tij_ab(a, b) - f(1)*fac_ij
2727 IF (task(2, 2) .OR. task(3, 1)) force_eval = force_eval .OR. do_forces
2729 IF (task(2, 2)) force_eval = force_eval .OR. do_efield2
2730 IF (task(3, 3)) force_eval = force_eval .OR. do_efield1
2731 IF (task(3, 2) .OR. force_eval)
THEN
2732 force_eval = do_stress
2736 tmp = rab(a)*rab(b)*rab(c)*fac_ij
2737 tij_abc(a, b, c) = -15.0_dp*tmp*f(3)
2738 tmp = 3.0_dp*f(2)*fac_ij
2739 IF (a == b) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(c)
2740 IF (a == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(b)
2741 IF (b == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(a)
2745 IF (task(3, 2)) force_eval = force_eval .OR. do_forces
2747 IF (task(3, 3) .OR. force_eval)
THEN
2748 force_eval = do_stress
2753 tmp = rab(a)*rab(b)*rab(c)*rab(d)*fac_ij
2754 tij_abcd(a, b, c, d) = 105.0_dp*tmp*f(4)
2755 tmp1 = 15.0_dp*f(3)*fac_ij
2756 tmp2 = 3.0_dp*f(2)*fac_ij
2758 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(c)*rab(d)
2759 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
2762 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(d)
2763 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
2765 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(c)
2767 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(d)
2768 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
2770 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(c)
2771 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(b)
2776 IF (task(3, 3)) force_eval = force_eval .OR. do_forces
2778 IF (force_eval)
THEN
2779 force_eval = do_stress
2785 tmp = rab(a)*rab(b)*rab(c)*rab(d)*rab(e)*fac_ij
2786 tij_abcde(a, b, c, d, e) = -945.0_dp*tmp*f(5)
2787 tmp1 = 105.0_dp*f(4)*fac_ij
2788 tmp2 = 15.0_dp*f(3)*fac_ij
2790 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(c)*rab(d)*rab(e)
2791 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
2792 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
2793 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
2796 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(d)*rab(e)
2797 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
2798 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
2799 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
2802 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(e)
2803 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
2804 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
2805 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
2808 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(d)
2809 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
2810 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
2811 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
2814 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(d)*rab(e)
2815 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
2818 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(e)
2819 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
2822 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(d)
2823 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
2825 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(e)
2826 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(d)
2827 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(c)
2845 IF (debug_this_module)
THEN
2853 IF (any(task(1, :)))
THEN
2854 ch_j = charges(atom_a)
2855 ch_i = charges(atom_b)
2857 IF (any(task(2, :)))
THEN
2858 dp_j = dipoles(:, atom_a)
2859 dp_i = dipoles(:, atom_b)
2861 IF (any(task(3, :)))
THEN
2862 qp_j = quadrupoles(:, :, atom_a)
2863 qp_i = quadrupoles(:, :, atom_b)
2865 IF (task(1, 1))
THEN
2867 eloc = eloc + ch_i*tij*ch_j
2869 IF (do_forces .OR. do_stress)
THEN
2870 fr(1) = fr(1) - ch_j*tij_a(1)*ch_i
2871 fr(2) = fr(2) - ch_j*tij_a(2)*ch_i
2872 fr(3) = fr(3) - ch_j*tij_a(3)*ch_i
2877 IF (do_efield0)
THEN
2878 ef0_i = ef0_i + tij*ch_j
2880 ef0_j = ef0_j + tij*ch_i
2883 IF (do_efield1)
THEN
2884 ef1_i(1) = ef1_i(1) - tij_a(1)*ch_j
2885 ef1_i(2) = ef1_i(2) - tij_a(2)*ch_j
2886 ef1_i(3) = ef1_i(3) - tij_a(3)*ch_j
2888 ef1_j(1) = ef1_j(1) + tij_a(1)*ch_i
2889 ef1_j(2) = ef1_j(2) + tij_a(2)*ch_i
2890 ef1_j(3) = ef1_j(3) + tij_a(3)*ch_i
2895 IF (do_efield2)
THEN
2896 ef2_i(1, 1) = ef2_i(1, 1) - tij_ab(1, 1)*ch_j
2897 ef2_i(2, 1) = ef2_i(2, 1) - tij_ab(2, 1)*ch_j
2898 ef2_i(3, 1) = ef2_i(3, 1) - tij_ab(3, 1)*ch_j
2899 ef2_i(1, 2) = ef2_i(1, 2) - tij_ab(1, 2)*ch_j
2900 ef2_i(2, 2) = ef2_i(2, 2) - tij_ab(2, 2)*ch_j
2901 ef2_i(3, 2) = ef2_i(3, 2) - tij_ab(3, 2)*ch_j
2902 ef2_i(1, 3) = ef2_i(1, 3) - tij_ab(1, 3)*ch_j
2903 ef2_i(2, 3) = ef2_i(2, 3) - tij_ab(2, 3)*ch_j
2904 ef2_i(3, 3) = ef2_i(3, 3) - tij_ab(3, 3)*ch_j
2906 ef2_j(1, 1) = ef2_j(1, 1) - tij_ab(1, 1)*ch_i
2907 ef2_j(2, 1) = ef2_j(2, 1) - tij_ab(2, 1)*ch_i
2908 ef2_j(3, 1) = ef2_j(3, 1) - tij_ab(3, 1)*ch_i
2909 ef2_j(1, 2) = ef2_j(1, 2) - tij_ab(1, 2)*ch_i
2910 ef2_j(2, 2) = ef2_j(2, 2) - tij_ab(2, 2)*ch_i
2911 ef2_j(3, 2) = ef2_j(3, 2) - tij_ab(3, 2)*ch_i
2912 ef2_j(1, 3) = ef2_j(1, 3) - tij_ab(1, 3)*ch_i
2913 ef2_j(2, 3) = ef2_j(2, 3) - tij_ab(2, 3)*ch_i
2914 ef2_j(3, 3) = ef2_j(3, 3) - tij_ab(3, 3)*ch_i
2918 IF (task(2, 2))
THEN
2920 tmp = -(dp_i(1)*(tij_ab(1, 1)*dp_j(1) + &
2921 tij_ab(2, 1)*dp_j(2) + &
2922 tij_ab(3, 1)*dp_j(3)) + &
2923 dp_i(2)*(tij_ab(1, 2)*dp_j(1) + &
2924 tij_ab(2, 2)*dp_j(2) + &
2925 tij_ab(3, 2)*dp_j(3)) + &
2926 dp_i(3)*(tij_ab(1, 3)*dp_j(1) + &
2927 tij_ab(2, 3)*dp_j(2) + &
2928 tij_ab(3, 3)*dp_j(3)))
2931 IF (do_forces .OR. do_stress)
THEN
2933 fr(k) = fr(k) + dp_i(1)*(tij_abc(1, 1, k)*dp_j(1) + &
2934 tij_abc(2, 1, k)*dp_j(2) + &
2935 tij_abc(3, 1, k)*dp_j(3)) &
2936 + dp_i(2)*(tij_abc(1, 2, k)*dp_j(1) + &
2937 tij_abc(2, 2, k)*dp_j(2) + &
2938 tij_abc(3, 2, k)*dp_j(3)) &
2939 + dp_i(3)*(tij_abc(1, 3, k)*dp_j(1) + &
2940 tij_abc(2, 3, k)*dp_j(2) + &
2941 tij_abc(3, 3, k)*dp_j(3))
2947 IF (do_efield0)
THEN
2948 ef0_i = ef0_i - (tij_a(1)*dp_j(1) + &
2949 tij_a(2)*dp_j(2) + &
2952 ef0_j = ef0_j + (tij_a(1)*dp_i(1) + &
2953 tij_a(2)*dp_i(2) + &
2957 IF (do_efield1)
THEN
2958 ef1_i(1) = ef1_i(1) + (tij_ab(1, 1)*dp_j(1) + &
2959 tij_ab(2, 1)*dp_j(2) + &
2960 tij_ab(3, 1)*dp_j(3))
2961 ef1_i(2) = ef1_i(2) + (tij_ab(1, 2)*dp_j(1) + &
2962 tij_ab(2, 2)*dp_j(2) + &
2963 tij_ab(3, 2)*dp_j(3))
2964 ef1_i(3) = ef1_i(3) + (tij_ab(1, 3)*dp_j(1) + &
2965 tij_ab(2, 3)*dp_j(2) + &
2966 tij_ab(3, 3)*dp_j(3))
2968 ef1_j(1) = ef1_j(1) + (tij_ab(1, 1)*dp_i(1) + &
2969 tij_ab(2, 1)*dp_i(2) + &
2970 tij_ab(3, 1)*dp_i(3))
2971 ef1_j(2) = ef1_j(2) + (tij_ab(1, 2)*dp_i(1) + &
2972 tij_ab(2, 2)*dp_i(2) + &
2973 tij_ab(3, 2)*dp_i(3))
2974 ef1_j(3) = ef1_j(3) + (tij_ab(1, 3)*dp_i(1) + &
2975 tij_ab(2, 3)*dp_i(2) + &
2976 tij_ab(3, 3)*dp_i(3))
2979 IF (do_efield2)
THEN
2980 ef2_i(1, 1) = ef2_i(1, 1) + (tij_abc(1, 1, 1)*dp_j(1) + &
2981 tij_abc(2, 1, 1)*dp_j(2) + &
2982 tij_abc(3, 1, 1)*dp_j(3))
2983 ef2_i(1, 2) = ef2_i(1, 2) + (tij_abc(1, 1, 2)*dp_j(1) + &
2984 tij_abc(2, 1, 2)*dp_j(2) + &
2985 tij_abc(3, 1, 2)*dp_j(3))
2986 ef2_i(1, 3) = ef2_i(1, 3) + (tij_abc(1, 1, 3)*dp_j(1) + &
2987 tij_abc(2, 1, 3)*dp_j(2) + &
2988 tij_abc(3, 1, 3)*dp_j(3))
2989 ef2_i(2, 1) = ef2_i(2, 1) + (tij_abc(1, 2, 1)*dp_j(1) + &
2990 tij_abc(2, 2, 1)*dp_j(2) + &
2991 tij_abc(3, 2, 1)*dp_j(3))
2992 ef2_i(2, 2) = ef2_i(2, 2) + (tij_abc(1, 2, 2)*dp_j(1) + &
2993 tij_abc(2, 2, 2)*dp_j(2) + &
2994 tij_abc(3, 2, 2)*dp_j(3))
2995 ef2_i(2, 3) = ef2_i(2, 3) + (tij_abc(1, 2, 3)*dp_j(1) + &
2996 tij_abc(2, 2, 3)*dp_j(2) + &
2997 tij_abc(3, 2, 3)*dp_j(3))
2998 ef2_i(3, 1) = ef2_i(3, 1) + (tij_abc(1, 3, 1)*dp_j(1) + &
2999 tij_abc(2, 3, 1)*dp_j(2) + &
3000 tij_abc(3, 3, 1)*dp_j(3))
3001 ef2_i(3, 2) = ef2_i(3, 2) + (tij_abc(1, 3, 2)*dp_j(1) + &
3002 tij_abc(2, 3, 2)*dp_j(2) + &
3003 tij_abc(3, 3, 2)*dp_j(3))
3004 ef2_i(3, 3) = ef2_i(3, 3) + (tij_abc(1, 3, 3)*dp_j(1) + &
3005 tij_abc(2, 3, 3)*dp_j(2) + &
3006 tij_abc(3, 3, 3)*dp_j(3))
3008 ef2_j(1, 1) = ef2_j(1, 1) - (tij_abc(1, 1, 1)*dp_i(1) + &
3009 tij_abc(2, 1, 1)*dp_i(2) + &
3010 tij_abc(3, 1, 1)*dp_i(3))
3011 ef2_j(1, 2) = ef2_j(1, 2) - (tij_abc(1, 1, 2)*dp_i(1) + &
3012 tij_abc(2, 1, 2)*dp_i(2) + &
3013 tij_abc(3, 1, 2)*dp_i(3))
3014 ef2_j(1, 3) = ef2_j(1, 3) - (tij_abc(1, 1, 3)*dp_i(1) + &
3015 tij_abc(2, 1, 3)*dp_i(2) + &
3016 tij_abc(3, 1, 3)*dp_i(3))
3017 ef2_j(2, 1) = ef2_j(2, 1) - (tij_abc(1, 2, 1)*dp_i(1) + &
3018 tij_abc(2, 2, 1)*dp_i(2) + &
3019 tij_abc(3, 2, 1)*dp_i(3))
3020 ef2_j(2, 2) = ef2_j(2, 2) - (tij_abc(1, 2, 2)*dp_i(1) + &
3021 tij_abc(2, 2, 2)*dp_i(2) + &
3022 tij_abc(3, 2, 2)*dp_i(3))
3023 ef2_j(2, 3) = ef2_j(2, 3) - (tij_abc(1, 2, 3)*dp_i(1) + &
3024 tij_abc(2, 2, 3)*dp_i(2) + &
3025 tij_abc(3, 2, 3)*dp_i(3))
3026 ef2_j(3, 1) = ef2_j(3, 1) - (tij_abc(1, 3, 1)*dp_i(1) + &
3027 tij_abc(2, 3, 1)*dp_i(2) + &
3028 tij_abc(3, 3, 1)*dp_i(3))
3029 ef2_j(3, 2) = ef2_j(3, 2) - (tij_abc(1, 3, 2)*dp_i(1) + &
3030 tij_abc(2, 3, 2)*dp_i(2) + &
3031 tij_abc(3, 3, 2)*dp_i(3))
3032 ef2_j(3, 3) = ef2_j(3, 3) - (tij_abc(1, 3, 3)*dp_i(1) + &
3033 tij_abc(2, 3, 3)*dp_i(2) + &
3034 tij_abc(3, 3, 3)*dp_i(3))
3038 IF (task(2, 1))
THEN
3040 tmp = ch_j*(tij_a(1)*dp_i(1) + &
3041 tij_a(2)*dp_i(2) + &
3043 - ch_i*(tij_a(1)*dp_j(1) + &
3044 tij_a(2)*dp_j(2) + &
3048 IF (do_forces .OR. do_stress)
THEN
3050 fr(k) = fr(k) - ch_j*(tij_ab(1, k)*dp_i(1) + &
3051 tij_ab(2, k)*dp_i(2) + &
3052 tij_ab(3, k)*dp_i(3)) &
3053 + ch_i*(tij_ab(1, k)*dp_j(1) + &
3054 tij_ab(2, k)*dp_j(2) + &
3055 tij_ab(3, k)*dp_j(3))
3059 IF (task(3, 3))
THEN
3062 tmp11 = qp_i(1, 1)*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
3063 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
3064 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
3065 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
3066 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
3067 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
3068 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
3069 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
3070 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
3071 tmp21 = qp_i(2, 1)*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
3072 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
3073 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
3074 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
3075 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
3076 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
3077 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
3078 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
3079 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
3080 tmp31 = qp_i(3, 1)*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
3081 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
3082 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
3083 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
3084 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
3085 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
3086 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
3087 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
3088 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
3089 tmp22 = qp_i(2, 2)*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
3090 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
3091 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
3092 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
3093 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
3094 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
3095 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
3096 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
3097 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
3098 tmp32 = qp_i(3, 2)*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
3099 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
3100 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
3101 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
3102 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
3103 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
3104 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
3105 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
3106 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
3107 tmp33 = qp_i(3, 3)*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
3108 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
3109 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
3110 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
3111 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
3112 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
3113 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
3114 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
3115 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
3119 tmp = tmp11 + tmp12 + tmp13 + &
3120 tmp21 + tmp22 + tmp23 + &
3121 tmp31 + tmp32 + tmp33
3123 eloc = eloc +
fac*tmp
3125 IF (do_forces .OR. do_stress)
THEN
3127 tmp11 = qp_i(1, 1)*(tij_abcde(1, 1, 1, 1, k)*qp_j(1, 1) + &
3128 tij_abcde(2, 1, 1, 1, k)*qp_j(2, 1) + &
3129 tij_abcde(3, 1, 1, 1, k)*qp_j(3, 1) + &
3130 tij_abcde(1, 2, 1, 1, k)*qp_j(1, 2) + &
3131 tij_abcde(2, 2, 1, 1, k)*qp_j(2, 2) + &
3132 tij_abcde(3, 2, 1, 1, k)*qp_j(3, 2) + &
3133 tij_abcde(1, 3, 1, 1, k)*qp_j(1, 3) + &
3134 tij_abcde(2, 3, 1, 1, k)*qp_j(2, 3) + &
3135 tij_abcde(3, 3, 1, 1, k)*qp_j(3, 3))
3136 tmp21 = qp_i(2, 1)*(tij_abcde(1, 1, 2, 1, k)*qp_j(1, 1) + &
3137 tij_abcde(2, 1, 2, 1, k)*qp_j(2, 1) + &
3138 tij_abcde(3, 1, 2, 1, k)*qp_j(3, 1) + &
3139 tij_abcde(1, 2, 2, 1, k)*qp_j(1, 2) + &
3140 tij_abcde(2, 2, 2, 1, k)*qp_j(2, 2) + &
3141 tij_abcde(3, 2, 2, 1, k)*qp_j(3, 2) + &
3142 tij_abcde(1, 3, 2, 1, k)*qp_j(1, 3) + &
3143 tij_abcde(2, 3, 2, 1, k)*qp_j(2, 3) + &
3144 tij_abcde(3, 3, 2, 1, k)*qp_j(3, 3))
3145 tmp31 = qp_i(3, 1)*(tij_abcde(1, 1, 3, 1, k)*qp_j(1, 1) + &
3146 tij_abcde(2, 1, 3, 1, k)*qp_j(2, 1) + &
3147 tij_abcde(3, 1, 3, 1, k)*qp_j(3, 1) + &
3148 tij_abcde(1, 2, 3, 1, k)*qp_j(1, 2) + &
3149 tij_abcde(2, 2, 3, 1, k)*qp_j(2, 2) + &
3150 tij_abcde(3, 2, 3, 1, k)*qp_j(3, 2) + &
3151 tij_abcde(1, 3, 3, 1, k)*qp_j(1, 3) + &
3152 tij_abcde(2, 3, 3, 1, k)*qp_j(2, 3) + &
3153 tij_abcde(3, 3, 3, 1, k)*qp_j(3, 3))
3154 tmp22 = qp_i(2, 2)*(tij_abcde(1, 1, 2, 2, k)*qp_j(1, 1) + &
3155 tij_abcde(2, 1, 2, 2, k)*qp_j(2, 1) + &
3156 tij_abcde(3, 1, 2, 2, k)*qp_j(3, 1) + &
3157 tij_abcde(1, 2, 2, 2, k)*qp_j(1, 2) + &
3158 tij_abcde(2, 2, 2, 2, k)*qp_j(2, 2) + &
3159 tij_abcde(3, 2, 2, 2, k)*qp_j(3, 2) + &
3160 tij_abcde(1, 3, 2, 2, k)*qp_j(1, 3) + &
3161 tij_abcde(2, 3, 2, 2, k)*qp_j(2, 3) + &
3162 tij_abcde(3, 3, 2, 2, k)*qp_j(3, 3))
3163 tmp32 = qp_i(3, 2)*(tij_abcde(1, 1, 3, 2, k)*qp_j(1, 1) + &
3164 tij_abcde(2, 1, 3, 2, k)*qp_j(2, 1) + &
3165 tij_abcde(3, 1, 3, 2, k)*qp_j(3, 1) + &
3166 tij_abcde(1, 2, 3, 2, k)*qp_j(1, 2) + &
3167 tij_abcde(2, 2, 3, 2, k)*qp_j(2, 2) + &
3168 tij_abcde(3, 2, 3, 2, k)*qp_j(3, 2) + &
3169 tij_abcde(1, 3, 3, 2, k)*qp_j(1, 3) + &
3170 tij_abcde(2, 3, 3, 2, k)*qp_j(2, 3) + &
3171 tij_abcde(3, 3, 3, 2, k)*qp_j(3, 3))
3172 tmp33 = qp_i(3, 3)*(tij_abcde(1, 1, 3, 3, k)*qp_j(1, 1) + &
3173 tij_abcde(2, 1, 3, 3, k)*qp_j(2, 1) + &
3174 tij_abcde(3, 1, 3, 3, k)*qp_j(3, 1) + &
3175 tij_abcde(1, 2, 3, 3, k)*qp_j(1, 2) + &
3176 tij_abcde(2, 2, 3, 3, k)*qp_j(2, 2) + &
3177 tij_abcde(3, 2, 3, 3, k)*qp_j(3, 2) + &
3178 tij_abcde(1, 3, 3, 3, k)*qp_j(1, 3) + &
3179 tij_abcde(2, 3, 3, 3, k)*qp_j(2, 3) + &
3180 tij_abcde(3, 3, 3, 3, k)*qp_j(3, 3))
3184 fr(k) = fr(k) -
fac*(tmp11 + tmp12 + tmp13 + &
3185 tmp21 + tmp22 + tmp23 + &
3186 tmp31 + tmp32 + tmp33)
3193 IF (do_efield0)
THEN
3194 ef0_i = ef0_i +
fac*(tij_ab(1, 1)*qp_j(1, 1) + &
3195 tij_ab(2, 1)*qp_j(2, 1) + &
3196 tij_ab(3, 1)*qp_j(3, 1) + &
3197 tij_ab(1, 2)*qp_j(1, 2) + &
3198 tij_ab(2, 2)*qp_j(2, 2) + &
3199 tij_ab(3, 2)*qp_j(3, 2) + &
3200 tij_ab(1, 3)*qp_j(1, 3) + &
3201 tij_ab(2, 3)*qp_j(2, 3) + &
3202 tij_ab(3, 3)*qp_j(3, 3))
3204 ef0_j = ef0_j +
fac*(tij_ab(1, 1)*qp_i(1, 1) + &
3205 tij_ab(2, 1)*qp_i(2, 1) + &
3206 tij_ab(3, 1)*qp_i(3, 1) + &
3207 tij_ab(1, 2)*qp_i(1, 2) + &
3208 tij_ab(2, 2)*qp_i(2, 2) + &
3209 tij_ab(3, 2)*qp_i(3, 2) + &
3210 tij_ab(1, 3)*qp_i(1, 3) + &
3211 tij_ab(2, 3)*qp_i(2, 3) + &
3212 tij_ab(3, 3)*qp_i(3, 3))
3215 IF (do_efield1)
THEN
3216 ef1_i(1) = ef1_i(1) -
fac*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
3217 tij_abc(2, 1, 1)*qp_j(2, 1) + &
3218 tij_abc(3, 1, 1)*qp_j(3, 1) + &
3219 tij_abc(1, 2, 1)*qp_j(1, 2) + &
3220 tij_abc(2, 2, 1)*qp_j(2, 2) + &
3221 tij_abc(3, 2, 1)*qp_j(3, 2) + &
3222 tij_abc(1, 3, 1)*qp_j(1, 3) + &
3223 tij_abc(2, 3, 1)*qp_j(2, 3) + &
3224 tij_abc(3, 3, 1)*qp_j(3, 3))
3225 ef1_i(2) = ef1_i(2) -
fac*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
3226 tij_abc(2, 1, 2)*qp_j(2, 1) + &
3227 tij_abc(3, 1, 2)*qp_j(3, 1) + &
3228 tij_abc(1, 2, 2)*qp_j(1, 2) + &
3229 tij_abc(2, 2, 2)*qp_j(2, 2) + &
3230 tij_abc(3, 2, 2)*qp_j(3, 2) + &
3231 tij_abc(1, 3, 2)*qp_j(1, 3) + &
3232 tij_abc(2, 3, 2)*qp_j(2, 3) + &
3233 tij_abc(3, 3, 2)*qp_j(3, 3))
3234 ef1_i(3) = ef1_i(3) -
fac*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
3235 tij_abc(2, 1, 3)*qp_j(2, 1) + &
3236 tij_abc(3, 1, 3)*qp_j(3, 1) + &
3237 tij_abc(1, 2, 3)*qp_j(1, 2) + &
3238 tij_abc(2, 2, 3)*qp_j(2, 2) + &
3239 tij_abc(3, 2, 3)*qp_j(3, 2) + &
3240 tij_abc(1, 3, 3)*qp_j(1, 3) + &
3241 tij_abc(2, 3, 3)*qp_j(2, 3) + &
3242 tij_abc(3, 3, 3)*qp_j(3, 3))
3244 ef1_j(1) = ef1_j(1) +
fac*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
3245 tij_abc(2, 1, 1)*qp_i(2, 1) + &
3246 tij_abc(3, 1, 1)*qp_i(3, 1) + &
3247 tij_abc(1, 2, 1)*qp_i(1, 2) + &
3248 tij_abc(2, 2, 1)*qp_i(2, 2) + &
3249 tij_abc(3, 2, 1)*qp_i(3, 2) + &
3250 tij_abc(1, 3, 1)*qp_i(1, 3) + &
3251 tij_abc(2, 3, 1)*qp_i(2, 3) + &
3252 tij_abc(3, 3, 1)*qp_i(3, 3))
3253 ef1_j(2) = ef1_j(2) +
fac*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
3254 tij_abc(2, 1, 2)*qp_i(2, 1) + &
3255 tij_abc(3, 1, 2)*qp_i(3, 1) + &
3256 tij_abc(1, 2, 2)*qp_i(1, 2) + &
3257 tij_abc(2, 2, 2)*qp_i(2, 2) + &
3258 tij_abc(3, 2, 2)*qp_i(3, 2) + &
3259 tij_abc(1, 3, 2)*qp_i(1, 3) + &
3260 tij_abc(2, 3, 2)*qp_i(2, 3) + &
3261 tij_abc(3, 3, 2)*qp_i(3, 3))
3262 ef1_j(3) = ef1_j(3) +
fac*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
3263 tij_abc(2, 1, 3)*qp_i(2, 1) + &
3264 tij_abc(3, 1, 3)*qp_i(3, 1) + &
3265 tij_abc(1, 2, 3)*qp_i(1, 2) + &
3266 tij_abc(2, 2, 3)*qp_i(2, 2) + &
3267 tij_abc(3, 2, 3)*qp_i(3, 2) + &
3268 tij_abc(1, 3, 3)*qp_i(1, 3) + &
3269 tij_abc(2, 3, 3)*qp_i(2, 3) + &
3270 tij_abc(3, 3, 3)*qp_i(3, 3))
3273 IF (do_efield2)
THEN
3274 tmp11 =
fac*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
3275 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
3276 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
3277 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
3278 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
3279 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
3280 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
3281 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
3282 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
3283 tmp12 =
fac*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
3284 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
3285 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
3286 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
3287 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
3288 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
3289 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
3290 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
3291 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
3292 tmp13 =
fac*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
3293 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
3294 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
3295 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
3296 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
3297 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
3298 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
3299 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
3300 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
3301 tmp22 =
fac*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
3302 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
3303 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
3304 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
3305 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
3306 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
3307 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
3308 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
3309 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
3310 tmp23 =
fac*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
3311 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
3312 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
3313 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
3314 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
3315 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
3316 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
3317 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
3318 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
3319 tmp33 =
fac*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
3320 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
3321 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
3322 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
3323 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
3324 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
3325 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
3326 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
3327 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
3329 ef2_i(1, 1) = ef2_i(1, 1) - tmp11
3330 ef2_i(1, 2) = ef2_i(1, 2) - tmp12
3331 ef2_i(1, 3) = ef2_i(1, 3) - tmp13
3332 ef2_i(2, 1) = ef2_i(2, 1) - tmp12
3333 ef2_i(2, 2) = ef2_i(2, 2) - tmp22
3334 ef2_i(2, 3) = ef2_i(2, 3) - tmp23
3335 ef2_i(3, 1) = ef2_i(3, 1) - tmp13
3336 ef2_i(3, 2) = ef2_i(3, 2) - tmp23
3337 ef2_i(3, 3) = ef2_i(3, 3) - tmp33
3339 tmp11 =
fac*(tij_abcd(1, 1, 1, 1)*qp_i(1, 1) + &
3340 tij_abcd(2, 1, 1, 1)*qp_i(2, 1) + &
3341 tij_abcd(3, 1, 1, 1)*qp_i(3, 1) + &
3342 tij_abcd(1, 2, 1, 1)*qp_i(1, 2) + &
3343 tij_abcd(2, 2, 1, 1)*qp_i(2, 2) + &
3344 tij_abcd(3, 2, 1, 1)*qp_i(3, 2) + &
3345 tij_abcd(1, 3, 1, 1)*qp_i(1, 3) + &
3346 tij_abcd(2, 3, 1, 1)*qp_i(2, 3) + &
3347 tij_abcd(3, 3, 1, 1)*qp_i(3, 3))
3348 tmp12 =
fac*(tij_abcd(1, 1, 1, 2)*qp_i(1, 1) + &
3349 tij_abcd(2, 1, 1, 2)*qp_i(2, 1) + &
3350 tij_abcd(3, 1, 1, 2)*qp_i(3, 1) + &
3351 tij_abcd(1, 2, 1, 2)*qp_i(1, 2) + &
3352 tij_abcd(2, 2, 1, 2)*qp_i(2, 2) + &
3353 tij_abcd(3, 2, 1, 2)*qp_i(3, 2) + &
3354 tij_abcd(1, 3, 1, 2)*qp_i(1, 3) + &
3355 tij_abcd(2, 3, 1, 2)*qp_i(2, 3) + &
3356 tij_abcd(3, 3, 1, 2)*qp_i(3, 3))
3357 tmp13 =
fac*(tij_abcd(1, 1, 1, 3)*qp_i(1, 1) + &
3358 tij_abcd(2, 1, 1, 3)*qp_i(2, 1) + &
3359 tij_abcd(3, 1, 1, 3)*qp_i(3, 1) + &
3360 tij_abcd(1, 2, 1, 3)*qp_i(1, 2) + &
3361 tij_abcd(2, 2, 1, 3)*qp_i(2, 2) + &
3362 tij_abcd(3, 2, 1, 3)*qp_i(3, 2) + &
3363 tij_abcd(1, 3, 1, 3)*qp_i(1, 3) + &
3364 tij_abcd(2, 3, 1, 3)*qp_i(2, 3) + &
3365 tij_abcd(3, 3, 1, 3)*qp_i(3, 3))
3366 tmp22 =
fac*(tij_abcd(1, 1, 2, 2)*qp_i(1, 1) + &
3367 tij_abcd(2, 1, 2, 2)*qp_i(2, 1) + &
3368 tij_abcd(3, 1, 2, 2)*qp_i(3, 1) + &
3369 tij_abcd(1, 2, 2, 2)*qp_i(1, 2) + &
3370 tij_abcd(2, 2, 2, 2)*qp_i(2, 2) + &
3371 tij_abcd(3, 2, 2, 2)*qp_i(3, 2) + &
3372 tij_abcd(1, 3, 2, 2)*qp_i(1, 3) + &
3373 tij_abcd(2, 3, 2, 2)*qp_i(2, 3) + &
3374 tij_abcd(3, 3, 2, 2)*qp_i(3, 3))
3375 tmp23 =
fac*(tij_abcd(1, 1, 2, 3)*qp_i(1, 1) + &
3376 tij_abcd(2, 1, 2, 3)*qp_i(2, 1) + &
3377 tij_abcd(3, 1, 2, 3)*qp_i(3, 1) + &
3378 tij_abcd(1, 2, 2, 3)*qp_i(1, 2) + &
3379 tij_abcd(2, 2, 2, 3)*qp_i(2, 2) + &
3380 tij_abcd(3, 2, 2, 3)*qp_i(3, 2) + &
3381 tij_abcd(1, 3, 2, 3)*qp_i(1, 3) + &
3382 tij_abcd(2, 3, 2, 3)*qp_i(2, 3) + &
3383 tij_abcd(3, 3, 2, 3)*qp_i(3, 3))
3384 tmp33 =
fac*(tij_abcd(1, 1, 3, 3)*qp_i(1, 1) + &
3385 tij_abcd(2, 1, 3, 3)*qp_i(2, 1) + &
3386 tij_abcd(3, 1, 3, 3)*qp_i(3, 1) + &
3387 tij_abcd(1, 2, 3, 3)*qp_i(1, 2) + &
3388 tij_abcd(2, 2, 3, 3)*qp_i(2, 2) + &
3389 tij_abcd(3, 2, 3, 3)*qp_i(3, 2) + &
3390 tij_abcd(1, 3, 3, 3)*qp_i(1, 3) + &
3391 tij_abcd(2, 3, 3, 3)*qp_i(2, 3) + &
3392 tij_abcd(3, 3, 3, 3)*qp_i(3, 3))
3394 ef2_j(1, 1) = ef2_j(1, 1) - tmp11
3395 ef2_j(1, 2) = ef2_j(1, 2) - tmp12
3396 ef2_j(1, 3) = ef2_j(1, 3) - tmp13
3397 ef2_j(2, 1) = ef2_j(2, 1) - tmp12
3398 ef2_j(2, 2) = ef2_j(2, 2) - tmp22
3399 ef2_j(2, 3) = ef2_j(2, 3) - tmp23
3400 ef2_j(3, 1) = ef2_j(3, 1) - tmp13
3401 ef2_j(3, 2) = ef2_j(3, 2) - tmp23
3402 ef2_j(3, 3) = ef2_j(3, 3) - tmp33
3406 IF (task(3, 2))
THEN
3410 tmp_ij = dp_i(1)*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
3411 tij_abc(2, 1, 1)*qp_j(2, 1) + &
3412 tij_abc(3, 1, 1)*qp_j(3, 1) + &
3413 tij_abc(1, 2, 1)*qp_j(1, 2) + &
3414 tij_abc(2, 2, 1)*qp_j(2, 2) + &
3415 tij_abc(3, 2, 1)*qp_j(3, 2) + &
3416 tij_abc(1, 3, 1)*qp_j(1, 3) + &
3417 tij_abc(2, 3, 1)*qp_j(2, 3) + &
3418 tij_abc(3, 3, 1)*qp_j(3, 3)) + &
3419 dp_i(2)*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
3420 tij_abc(2, 1, 2)*qp_j(2, 1) + &
3421 tij_abc(3, 1, 2)*qp_j(3, 1) + &
3422 tij_abc(1, 2, 2)*qp_j(1, 2) + &
3423 tij_abc(2, 2, 2)*qp_j(2, 2) + &
3424 tij_abc(3, 2, 2)*qp_j(3, 2) + &
3425 tij_abc(1, 3, 2)*qp_j(1, 3) + &
3426 tij_abc(2, 3, 2)*qp_j(2, 3) + &
3427 tij_abc(3, 3, 2)*qp_j(3, 3)) + &
3428 dp_i(3)*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
3429 tij_abc(2, 1, 3)*qp_j(2, 1) + &
3430 tij_abc(3, 1, 3)*qp_j(3, 1) + &
3431 tij_abc(1, 2, 3)*qp_j(1, 2) + &
3432 tij_abc(2, 2, 3)*qp_j(2, 2) + &
3433 tij_abc(3, 2, 3)*qp_j(3, 2) + &
3434 tij_abc(1, 3, 3)*qp_j(1, 3) + &
3435 tij_abc(2, 3, 3)*qp_j(2, 3) + &
3436 tij_abc(3, 3, 3)*qp_j(3, 3))
3439 tmp_ji = dp_j(1)*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
3440 tij_abc(2, 1, 1)*qp_i(2, 1) + &
3441 tij_abc(3, 1, 1)*qp_i(3, 1) + &
3442 tij_abc(1, 2, 1)*qp_i(1, 2) + &
3443 tij_abc(2, 2, 1)*qp_i(2, 2) + &
3444 tij_abc(3, 2, 1)*qp_i(3, 2) + &
3445 tij_abc(1, 3, 1)*qp_i(1, 3) + &
3446 tij_abc(2, 3, 1)*qp_i(2, 3) + &
3447 tij_abc(3, 3, 1)*qp_i(3, 3)) + &
3448 dp_j(2)*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
3449 tij_abc(2, 1, 2)*qp_i(2, 1) + &
3450 tij_abc(3, 1, 2)*qp_i(3, 1) + &
3451 tij_abc(1, 2, 2)*qp_i(1, 2) + &
3452 tij_abc(2, 2, 2)*qp_i(2, 2) + &
3453 tij_abc(3, 2, 2)*qp_i(3, 2) + &
3454 tij_abc(1, 3, 2)*qp_i(1, 3) + &
3455 tij_abc(2, 3, 2)*qp_i(2, 3) + &
3456 tij_abc(3, 3, 2)*qp_i(3, 3)) + &
3457 dp_j(3)*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
3458 tij_abc(2, 1, 3)*qp_i(2, 1) + &
3459 tij_abc(3, 1, 3)*qp_i(3, 1) + &
3460 tij_abc(1, 2, 3)*qp_i(1, 2) + &
3461 tij_abc(2, 2, 3)*qp_i(2, 2) + &
3462 tij_abc(3, 2, 3)*qp_i(3, 2) + &
3463 tij_abc(1, 3, 3)*qp_i(1, 3) + &
3464 tij_abc(2, 3, 3)*qp_i(2, 3) + &
3465 tij_abc(3, 3, 3)*qp_i(3, 3))
3467 tmp =
fac*(tmp_ij - tmp_ji)
3469 IF (do_forces .OR. do_stress)
THEN
3472 tmp_ij = dp_i(1)*(tij_abcd(1, 1, 1, k)*qp_j(1, 1) + &
3473 tij_abcd(2, 1, 1, k)*qp_j(2, 1) + &
3474 tij_abcd(3, 1, 1, k)*qp_j(3, 1) + &
3475 tij_abcd(1, 2, 1, k)*qp_j(1, 2) + &
3476 tij_abcd(2, 2, 1, k)*qp_j(2, 2) + &
3477 tij_abcd(3, 2, 1, k)*qp_j(3, 2) + &
3478 tij_abcd(1, 3, 1, k)*qp_j(1, 3) + &
3479 tij_abcd(2, 3, 1, k)*qp_j(2, 3) + &
3480 tij_abcd(3, 3, 1, k)*qp_j(3, 3)) + &
3481 dp_i(2)*(tij_abcd(1, 1, 2, k)*qp_j(1, 1) + &
3482 tij_abcd(2, 1, 2, k)*qp_j(2, 1) + &
3483 tij_abcd(3, 1, 2, k)*qp_j(3, 1) + &
3484 tij_abcd(1, 2, 2, k)*qp_j(1, 2) + &
3485 tij_abcd(2, 2, 2, k)*qp_j(2, 2) + &
3486 tij_abcd(3, 2, 2, k)*qp_j(3, 2) + &
3487 tij_abcd(1, 3, 2, k)*qp_j(1, 3) + &
3488 tij_abcd(2, 3, 2, k)*qp_j(2, 3) + &
3489 tij_abcd(3, 3, 2, k)*qp_j(3, 3)) + &
3490 dp_i(3)*(tij_abcd(1, 1, 3, k)*qp_j(1, 1) + &
3491 tij_abcd(2, 1, 3, k)*qp_j(2, 1) + &
3492 tij_abcd(3, 1, 3, k)*qp_j(3, 1) + &
3493 tij_abcd(1, 2, 3, k)*qp_j(1, 2) + &
3494 tij_abcd(2, 2, 3, k)*qp_j(2, 2) + &
3495 tij_abcd(3, 2, 3, k)*qp_j(3, 2) + &
3496 tij_abcd(1, 3, 3, k)*qp_j(1, 3) + &
3497 tij_abcd(2, 3, 3, k)*qp_j(2, 3) + &
3498 tij_abcd(3, 3, 3, k)*qp_j(3, 3))
3501 tmp_ji = dp_j(1)*(tij_abcd(1, 1, 1, k)*qp_i(1, 1) + &
3502 tij_abcd(2, 1, 1, k)*qp_i(2, 1) + &
3503 tij_abcd(3, 1, 1, k)*qp_i(3, 1) + &
3504 tij_abcd(1, 2, 1, k)*qp_i(1, 2) + &
3505 tij_abcd(2, 2, 1, k)*qp_i(2, 2) + &
3506 tij_abcd(3, 2, 1, k)*qp_i(3, 2) + &
3507 tij_abcd(1, 3, 1, k)*qp_i(1, 3) + &
3508 tij_abcd(2, 3, 1, k)*qp_i(2, 3) + &
3509 tij_abcd(3, 3, 1, k)*qp_i(3, 3)) + &
3510 dp_j(2)*(tij_abcd(1, 1, 2, k)*qp_i(1, 1) + &
3511 tij_abcd(2, 1, 2, k)*qp_i(2, 1) + &
3512 tij_abcd(3, 1, 2, k)*qp_i(3, 1) + &
3513 tij_abcd(1, 2, 2, k)*qp_i(1, 2) + &
3514 tij_abcd(2, 2, 2, k)*qp_i(2, 2) + &
3515 tij_abcd(3, 2, 2, k)*qp_i(3, 2) + &
3516 tij_abcd(1, 3, 2, k)*qp_i(1, 3) + &
3517 tij_abcd(2, 3, 2, k)*qp_i(2, 3) + &
3518 tij_abcd(3, 3, 2, k)*qp_i(3, 3)) + &
3519 dp_j(3)*(tij_abcd(1, 1, 3, k)*qp_i(1, 1) + &
3520 tij_abcd(2, 1, 3, k)*qp_i(2, 1) + &
3521 tij_abcd(3, 1, 3, k)*qp_i(3, 1) + &
3522 tij_abcd(1, 2, 3, k)*qp_i(1, 2) + &
3523 tij_abcd(2, 2, 3, k)*qp_i(2, 2) + &
3524 tij_abcd(3, 2, 3, k)*qp_i(3, 2) + &
3525 tij_abcd(1, 3, 3, k)*qp_i(1, 3) + &
3526 tij_abcd(2, 3, 3, k)*qp_i(2, 3) + &
3527 tij_abcd(3, 3, 3, k)*qp_i(3, 3))
3529 fr(k) = fr(k) -
fac*(tmp_ij - tmp_ji)
3533 IF (task(3, 1))
THEN
3538 tmp_ij = ch_i*(tij_ab(1, 1)*qp_j(1, 1) + &
3539 tij_ab(2, 1)*qp_j(2, 1) + &
3540 tij_ab(3, 1)*qp_j(3, 1) + &
3541 tij_ab(1, 2)*qp_j(1, 2) + &
3542 tij_ab(2, 2)*qp_j(2, 2) + &
3543 tij_ab(3, 2)*qp_j(3, 2) + &
3544 tij_ab(1, 3)*qp_j(1, 3) + &
3545 tij_ab(2, 3)*qp_j(2, 3) + &
3546 tij_ab(3, 3)*qp_j(3, 3))
3549 tmp_ji = ch_j*(tij_ab(1, 1)*qp_i(1, 1) + &
3550 tij_ab(2, 1)*qp_i(2, 1) + &
3551 tij_ab(3, 1)*qp_i(3, 1) + &
3552 tij_ab(1, 2)*qp_i(1, 2) + &
3553 tij_ab(2, 2)*qp_i(2, 2) + &
3554 tij_ab(3, 2)*qp_i(3, 2) + &
3555 tij_ab(1, 3)*qp_i(1, 3) + &
3556 tij_ab(2, 3)*qp_i(2, 3) + &
3557 tij_ab(3, 3)*qp_i(3, 3))
3559 eloc = eloc +
fac*(tmp_ij + tmp_ji)
3560 IF (do_forces .OR. do_stress)
THEN
3563 tmp_ij = ch_i*(tij_abc(1, 1, k)*qp_j(1, 1) + &
3564 tij_abc(2, 1, k)*qp_j(2, 1) + &
3565 tij_abc(3, 1, k)*qp_j(3, 1) + &
3566 tij_abc(1, 2, k)*qp_j(1, 2) + &
3567 tij_abc(2, 2, k)*qp_j(2, 2) + &
3568 tij_abc(3, 2, k)*qp_j(3, 2) + &
3569 tij_abc(1, 3, k)*qp_j(1, 3) + &
3570 tij_abc(2, 3, k)*qp_j(2, 3) + &
3571 tij_abc(3, 3, k)*qp_j(3, 3))
3574 tmp_ji = ch_j*(tij_abc(1, 1, k)*qp_i(1, 1) + &
3575 tij_abc(2, 1, k)*qp_i(2, 1) + &
3576 tij_abc(3, 1, k)*qp_i(3, 1) + &
3577 tij_abc(1, 2, k)*qp_i(1, 2) + &
3578 tij_abc(2, 2, k)*qp_i(2, 2) + &
3579 tij_abc(3, 2, k)*qp_i(3, 2) + &
3580 tij_abc(1, 3, k)*qp_i(1, 3) + &
3581 tij_abc(2, 3, k)*qp_i(2, 3) + &
3582 tij_abc(3, 3, k)*qp_i(3, 3))
3584 fr(k) = fr(k) -
fac*(tmp_ij + tmp_ji)
3588 energy = energy + eloc
3590 forces(1, atom_a) = forces(1, atom_a) - fr(1)
3591 forces(2, atom_a) = forces(2, atom_a) - fr(2)
3592 forces(3, atom_a) = forces(3, atom_a) - fr(3)
3593 forces(1, atom_b) = forces(1, atom_b) + fr(1)
3594 forces(2, atom_b) = forces(2, atom_b) + fr(2)
3595 forces(3, atom_b) = forces(3, atom_b) + fr(3)
3600 IF (do_efield0)
THEN
3601 efield0(atom_a) = efield0(atom_a) + ef0_j
3603 efield0(atom_b) = efield0(atom_b) + ef0_i
3606 IF (do_efield1)
THEN
3607 efield1(1, atom_a) = efield1(1, atom_a) + ef1_j(1)
3608 efield1(2, atom_a) = efield1(2, atom_a) + ef1_j(2)
3609 efield1(3, atom_a) = efield1(3, atom_a) + ef1_j(3)
3611 efield1(1, atom_b) = efield1(1, atom_b) + ef1_i(1)
3612 efield1(2, atom_b) = efield1(2, atom_b) + ef1_i(2)
3613 efield1(3, atom_b) = efield1(3, atom_b) + ef1_i(3)
3616 IF (do_efield2)
THEN
3617 efield2(1, atom_a) = efield2(1, atom_a) + ef2_j(1, 1)
3618 efield2(2, atom_a) = efield2(2, atom_a) + ef2_j(1, 2)
3619 efield2(3, atom_a) = efield2(3, atom_a) + ef2_j(1, 3)
3620 efield2(4, atom_a) = efield2(4, atom_a) + ef2_j(2, 1)
3621 efield2(5, atom_a) = efield2(5, atom_a) + ef2_j(2, 2)
3622 efield2(6, atom_a) = efield2(6, atom_a) + ef2_j(2, 3)
3623 efield2(7, atom_a) = efield2(7, atom_a) + ef2_j(3, 1)
3624 efield2(8, atom_a) = efield2(8, atom_a) + ef2_j(3, 2)
3625 efield2(9, atom_a) = efield2(9, atom_a) + ef2_j(3, 3)
3627 efield2(1, atom_b) = efield2(1, atom_b) + ef2_i(1, 1)
3628 efield2(2, atom_b) = efield2(2, atom_b) + ef2_i(1, 2)
3629 efield2(3, atom_b) = efield2(3, atom_b) + ef2_i(1, 3)
3630 efield2(4, atom_b) = efield2(4, atom_b) + ef2_i(2, 1)
3631 efield2(5, atom_b) = efield2(5, atom_b) + ef2_i(2, 2)
3632 efield2(6, atom_b) = efield2(6, atom_b) + ef2_i(2, 3)
3633 efield2(7, atom_b) = efield2(7, atom_b) + ef2_i(3, 1)
3634 efield2(8, atom_b) = efield2(8, atom_b) + ef2_i(3, 2)
3635 efield2(9, atom_b) = efield2(9, atom_b) + ef2_i(3, 3)
3639 ptens11 = ptens11 + rab(1)*fr(1)
3640 ptens21 = ptens21 + rab(2)*fr(1)
3641 ptens31 = ptens31 + rab(3)*fr(1)
3642 ptens12 = ptens12 + rab(1)*fr(2)
3643 ptens22 = ptens22 + rab(2)*fr(2)
3644 ptens32 = ptens32 + rab(3)*fr(2)
3645 ptens13 = ptens13 + rab(1)*fr(3)
3646 ptens23 = ptens23 + rab(2)*fr(3)
3647 ptens33 = ptens33 + rab(3)*fr(3)
3651 END DO kind_group_loop
3654 pv(1, 1) = pv(1, 1) + ptens11
3655 pv(1, 2) = pv(1, 2) + (ptens12 + ptens21)*0.5_dp
3656 pv(1, 3) = pv(1, 3) + (ptens13 + ptens31)*0.5_dp
3658 pv(2, 2) = pv(2, 2) + ptens22
3659 pv(2, 3) = pv(2, 3) + (ptens23 + ptens32)*0.5_dp
3662 pv(3, 3) = pv(3, 3) + ptens33
3665 CALL timestop(handle)
3666 END SUBROUTINE ewald_multipole_bonded
3691 SUBROUTINE ewald_multipole_lr(ewald_env, ewald_pw, cell, particle_set, &
3692 local_particles, energy, task, do_forces, do_efield, do_stress, &
3693 charges, dipoles, quadrupoles, forces, pv, efield0, efield1, efield2)
3699 REAL(kind=
dp),
INTENT(INOUT) :: energy
3700 LOGICAL,
DIMENSION(3, 3),
INTENT(IN) :: task
3701 LOGICAL,
INTENT(IN) :: do_forces, do_efield, do_stress
3702 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: charges
3703 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
3704 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
3705 POINTER :: quadrupoles
3706 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT), &
3707 OPTIONAL :: forces, pv
3708 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0
3709 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: efield1, efield2
3711 CHARACTER(len=*),
PARAMETER :: routinen =
'ewald_multipole_LR'
3713 COMPLEX(KIND=dp) :: atm_factor, atm_factor_st(3), cnjg_fac, &
3715 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: summe_ef
3716 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:, :) :: summe_st
3717 INTEGER :: gpt, handle, iparticle, iparticle_kind, iparticle_local, lp, mp, nnodes, &
3718 node, np, nparticle_kind, nparticle_local
3719 INTEGER,
DIMENSION(:, :),
POINTER :: bds
3720 LOGICAL :: do_efield0, do_efield1, do_efield2
3721 REAL(kind=
dp) :: alpha, denom, dipole_t(3), f0, factor, &
3722 four_alpha_sq, gauss, pref, q_t, tmp, &
3724 REAL(kind=
dp),
DIMENSION(3) :: tmp_v, vec
3725 REAL(kind=
dp),
DIMENSION(3, 3) :: pv_tmp
3726 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: rho0
3734 CALL timeset(routinen, handle)
3735 do_efield0 = do_efield .AND.
ASSOCIATED(efield0)
3736 do_efield1 = do_efield .AND.
ASSOCIATED(efield1)
3737 do_efield2 = do_efield .AND.
ASSOCIATED(efield2)
3741 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, dg=dg)
3742 CALL dg_get(dg, dg_rho0=dg_rho0)
3743 rho0 => dg_rho0%density%array
3744 pw_grid => pw_pool%pw_grid
3745 bds => pw_grid%bounds
3748 nparticle_kind =
SIZE(local_particles%n_el)
3750 DO iparticle_kind = 1, nparticle_kind
3751 nnodes = nnodes + local_particles%n_el(iparticle_kind)
3755 ALLOCATE (summe_ef(1:pw_grid%ngpts_cut))
3760 ALLOCATE (summe_st(3, 1:pw_grid%ngpts_cut))
3765 four_alpha_sq = 4.0_dp*alpha**2
3771 DO iparticle_kind = 1, nparticle_kind
3772 nparticle_local = local_particles%n_el(iparticle_kind)
3773 DO iparticle_local = 1, nparticle_local
3775 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
3776 vec = matmul(cell%h_inv, particle_set(iparticle)%r)
3778 exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
3781 IF (any(task(1, :)))
THEN
3782 q_t = q_t + charges(iparticle)
3784 IF (any(task(2, :)))
THEN
3785 dipole_t = dipole_t + dipoles(:, iparticle)
3787 IF (any(task(3, :)))
THEN
3788 trq_t = trq_t + quadrupoles(1, 1, iparticle) + &
3789 quadrupoles(2, 2, iparticle) + &
3790 quadrupoles(3, 3, iparticle)
3796 DO gpt = 1, pw_grid%ngpts_cut_local
3797 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
3798 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
3799 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
3807 DO iparticle_kind = 1, nparticle_kind
3808 nparticle_local = local_particles%n_el(iparticle_kind)
3809 DO iparticle_local = 1, nparticle_local
3811 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
3813 CALL get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges, &
3814 dipoles, quadrupoles)
3815 summe_tmp = exp_igr%ex(lp, node)*exp_igr%ey(mp, node)*exp_igr%ez(np, node)
3816 summe_ef(gpt) = summe_ef(gpt) + atm_factor*summe_tmp
3820 CALL get_atom_factor_stress(atm_factor_st, pw_grid, gpt, iparticle, task, &
3821 dipoles, quadrupoles)
3822 summe_st(1:3, gpt) = summe_st(1:3, gpt) + atm_factor_st(1:3)*summe_tmp
3828 CALL group%sum(dipole_t)
3829 CALL group%sum(trq_t)
3830 CALL group%sum(summe_ef)
3831 IF (do_stress)
CALL group%sum(summe_st)
3834 DO gpt = 1, pw_grid%ngpts_cut_local
3836 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
3837 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
3838 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
3844 IF (pw_grid%gsq(gpt) == 0.0_dp)
THEN
3846 energy = energy + (1.0_dp/6.0_dp)*dot_product(dipole_t, dipole_t) &
3847 - (1.0_dp/9.0_dp)*q_t*trq_t
3850 pv_tmp(1, 1) = pv_tmp(1, 1) + (1.0_dp/6.0_dp)*dot_product(dipole_t, dipole_t)
3851 pv_tmp(2, 2) = pv_tmp(2, 2) + (1.0_dp/6.0_dp)*dot_product(dipole_t, dipole_t)
3852 pv_tmp(3, 3) = pv_tmp(3, 3) + (1.0_dp/6.0_dp)*dot_product(dipole_t, dipole_t)
3855 IF (do_efield .AND. (debug_e_field_en .OR. (.NOT. debug_this_module)))
THEN
3860 DO iparticle_kind = 1, nparticle_kind
3861 nparticle_local = local_particles%n_el(iparticle_kind)
3862 DO iparticle_local = 1, nparticle_local
3864 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
3867 IF (do_efield0)
THEN
3868 efield0(iparticle) = efield0(iparticle)
3871 IF (do_efield1)
THEN
3872 efield1(1:3, iparticle) = efield1(1:3, iparticle) - (1.0_dp/6.0_dp)*dipole_t(1:3)
3875 IF (do_efield2)
THEN
3876 efield2(1, iparticle) = efield2(1, iparticle) - (1.0_dp/(18.0_dp))*q_t
3877 efield2(5, iparticle) = efield2(5, iparticle) - (1.0_dp/(18.0_dp))*q_t
3878 efield2(9, iparticle) = efield2(9, iparticle) - (1.0_dp/(18.0_dp))*q_t
3885 gauss = (rho0(lp, mp, np)*pw_grid%vol)**2/pw_grid%gsq(gpt)
3886 factor = gauss*real(summe_ef(gpt)*conjg(summe_ef(gpt)), kind=
dp)
3887 energy = energy + factor
3889 IF (do_forces .OR. do_efield)
THEN
3891 DO iparticle_kind = 1, nparticle_kind
3892 nparticle_local = local_particles%n_el(iparticle_kind)
3893 DO iparticle_local = 1, nparticle_local
3895 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
3896 fac = exp_igr%ex(lp, node)*exp_igr%ey(mp, node)*exp_igr%ez(np, node)
3897 cnjg_fac = conjg(
fac)
3901 CALL get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges, &
3902 dipoles, quadrupoles)
3904 tmp = gauss*aimag(summe_ef(gpt)*(cnjg_fac*conjg(atm_factor)))
3905 forces(1, node) = forces(1, node) + tmp*pw_grid%g(1, gpt)
3906 forces(2, node) = forces(2, node) + tmp*pw_grid%g(2, gpt)
3907 forces(3, node) = forces(3, node) + tmp*pw_grid%g(3, gpt)
3913 IF (do_efield0)
THEN
3914 efield0(iparticle) = efield0(iparticle) + gauss*real(
fac*conjg(summe_ef(gpt)), kind=
dp)
3917 IF (do_efield1)
THEN
3918 tmp = aimag(
fac*conjg(summe_ef(gpt)))*gauss
3919 efield1(1, iparticle) = efield1(1, iparticle) - tmp*pw_grid%g(1, gpt)
3920 efield1(2, iparticle) = efield1(2, iparticle) - tmp*pw_grid%g(2, gpt)
3921 efield1(3, iparticle) = efield1(3, iparticle) - tmp*pw_grid%g(3, gpt)
3924 IF (do_efield2)
THEN
3925 tmp_v(1) = real(
fac*conjg(summe_ef(gpt)), kind=
dp)*pw_grid%g(1, gpt)*gauss
3926 tmp_v(2) = real(
fac*conjg(summe_ef(gpt)), kind=
dp)*pw_grid%g(2, gpt)*gauss
3927 tmp_v(3) = real(
fac*conjg(summe_ef(gpt)), kind=
dp)*pw_grid%g(3, gpt)*gauss
3929 efield2(1, iparticle) = efield2(1, iparticle) + tmp_v(1)*pw_grid%g(1, gpt)
3930 efield2(2, iparticle) = efield2(2, iparticle) + tmp_v(1)*pw_grid%g(2, gpt)
3931 efield2(3, iparticle) = efield2(3, iparticle) + tmp_v(1)*pw_grid%g(3, gpt)
3932 efield2(4, iparticle) = efield2(4, iparticle) + tmp_v(2)*pw_grid%g(1, gpt)
3933 efield2(5, iparticle) = efield2(5, iparticle) + tmp_v(2)*pw_grid%g(2, gpt)
3934 efield2(6, iparticle) = efield2(6, iparticle) + tmp_v(2)*pw_grid%g(3, gpt)
3935 efield2(7, iparticle) = efield2(7, iparticle) + tmp_v(3)*pw_grid%g(1, gpt)
3936 efield2(8, iparticle) = efield2(8, iparticle) + tmp_v(3)*pw_grid%g(2, gpt)
3937 efield2(9, iparticle) = efield2(9, iparticle) + tmp_v(3)*pw_grid%g(3, gpt)
3948 denom = 1.0_dp/four_alpha_sq + 1.0_dp/pw_grid%gsq(gpt)
3949 pv_tmp(1, 1) = pv_tmp(1, 1) + factor*(1.0_dp - 2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(1, gpt)*denom)
3950 pv_tmp(1, 2) = pv_tmp(1, 2) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(2, gpt)*denom)
3951 pv_tmp(1, 3) = pv_tmp(1, 3) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(3, gpt)*denom)
3952 pv_tmp(2, 1) = pv_tmp(2, 1) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(1, gpt)*denom)
3953 pv_tmp(2, 2) = pv_tmp(2, 2) + factor*(1.0_dp - 2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(2, gpt)*denom)
3954 pv_tmp(2, 3) = pv_tmp(2, 3) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(3, gpt)*denom)
3955 pv_tmp(3, 1) = pv_tmp(3, 1) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(1, gpt)*denom)
3956 pv_tmp(3, 2) = pv_tmp(3, 2) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(2, gpt)*denom)
3957 pv_tmp(3, 3) = pv_tmp(3, 3) + factor*(1.0_dp - 2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(3, gpt)*denom)
3960 pv_tmp(1, 1) = pv_tmp(1, 1) + f0*pw_grid%g(1, gpt)*real(summe_st(1, gpt)*conjg(summe_ef(gpt)), kind=
dp)
3961 pv_tmp(1, 2) = pv_tmp(1, 2) + f0*pw_grid%g(1, gpt)*real(summe_st(2, gpt)*conjg(summe_ef(gpt)), kind=
dp)
3962 pv_tmp(1, 3) = pv_tmp(1, 3) + f0*pw_grid%g(1, gpt)*real(summe_st(3, gpt)*conjg(summe_ef(gpt)), kind=
dp)
3963 pv_tmp(2, 1) = pv_tmp(2, 1) + f0*pw_grid%g(2, gpt)*real(summe_st(1, gpt)*conjg(summe_ef(gpt)), kind=
dp)
3964 pv_tmp(2, 2) = pv_tmp(2, 2) + f0*pw_grid%g(2, gpt)*real(summe_st(2, gpt)*conjg(summe_ef(gpt)), kind=
dp)
3965 pv_tmp(2, 3) = pv_tmp(2, 3) + f0*pw_grid%g(2, gpt)*real(summe_st(3, gpt)*conjg(summe_ef(gpt)), kind=
dp)
3966 pv_tmp(3, 1) = pv_tmp(3, 1) + f0*pw_grid%g(3, gpt)*real(summe_st(1, gpt)*conjg(summe_ef(gpt)), kind=
dp)
3967 pv_tmp(3, 2) = pv_tmp(3, 2) + f0*pw_grid%g(3, gpt)*real(summe_st(2, gpt)*conjg(summe_ef(gpt)), kind=
dp)
3968 pv_tmp(3, 3) = pv_tmp(3, 3) + f0*pw_grid%g(3, gpt)*real(summe_st(3, gpt)*conjg(summe_ef(gpt)), kind=
dp)
3971 pref =
fourpi/pw_grid%vol
3972 energy = energy*pref
3975 DEALLOCATE (summe_ef)
3977 pv_tmp = pv_tmp*pref
3979 pv(1, 1) = pv(1, 1) + pv_tmp(1, 1)
3980 pv(1, 2) = pv(1, 2) + (pv_tmp(1, 2) + pv_tmp(2, 1))*0.5_dp
3981 pv(1, 3) = pv(1, 3) + (pv_tmp(1, 3) + pv_tmp(3, 1))*0.5_dp
3983 pv(2, 2) = pv(2, 2) + pv_tmp(2, 2)
3984 pv(2, 3) = pv(2, 3) + (pv_tmp(2, 3) + pv_tmp(3, 2))*0.5_dp
3987 pv(3, 3) = pv(3, 3) + pv_tmp(3, 3)
3988 DEALLOCATE (summe_st)
3991 forces = 2.0_dp*forces*pref
3993 IF (do_efield0)
THEN
3994 efield0 = 2.0_dp*efield0*pref
3996 IF (do_efield1)
THEN
3997 efield1 = 2.0_dp*efield1*pref
3999 IF (do_efield2)
THEN
4000 efield2 = 2.0_dp*efield2*pref
4002 CALL timestop(handle)
4004 END SUBROUTINE ewald_multipole_lr
4020 SUBROUTINE get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges, &
4021 dipoles, quadrupoles)
4022 COMPLEX(KIND=dp),
INTENT(OUT) :: atm_factor
4024 INTEGER,
INTENT(IN) :: gpt
4025 INTEGER :: iparticle
4026 LOGICAL,
DIMENSION(3, 3),
INTENT(IN) :: task
4027 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: charges
4028 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
4029 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
4030 POINTER :: quadrupoles
4032 COMPLEX(KIND=dp) :: tmp
4036 IF (task(1, 1))
THEN
4038 atm_factor = atm_factor + charges(iparticle)
4040 IF (task(2, 2))
THEN
4044 tmp = tmp + dipoles(i, iparticle)*pw_grid%g(i, gpt)
4046 atm_factor = atm_factor + tmp*cmplx(0.0_dp, -1.0_dp, kind=
dp)
4048 IF (task(3, 3))
THEN
4053 tmp = tmp + quadrupoles(j, i, iparticle)*pw_grid%g(j, gpt)*pw_grid%g(i, gpt)
4056 atm_factor = atm_factor - 1.0_dp/3.0_dp*tmp
4059 END SUBROUTINE get_atom_factor
4074 SUBROUTINE get_atom_factor_stress(atm_factor, pw_grid, gpt, iparticle, task, &
4075 dipoles, quadrupoles)
4076 COMPLEX(KIND=dp),
INTENT(OUT) :: atm_factor(3)
4078 INTEGER,
INTENT(IN) :: gpt
4079 INTEGER :: iparticle
4080 LOGICAL,
DIMENSION(3, 3),
INTENT(IN) :: task
4081 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
4082 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
4083 POINTER :: quadrupoles
4088 IF (any(task(2, :)))
THEN
4090 atm_factor = dipoles(:, iparticle)*cmplx(0.0_dp, -1.0_dp, kind=
dp)
4092 IF (any(task(3, :)))
THEN
4095 atm_factor(1) = atm_factor(1) - 1.0_dp/3.0_dp* &
4096 (quadrupoles(1, i, iparticle)*pw_grid%g(i, gpt) + &
4097 quadrupoles(i, 1, iparticle)*pw_grid%g(i, gpt))
4098 atm_factor(2) = atm_factor(2) - 1.0_dp/3.0_dp* &
4099 (quadrupoles(2, i, iparticle)*pw_grid%g(i, gpt) + &
4100 quadrupoles(i, 2, iparticle)*pw_grid%g(i, gpt))
4101 atm_factor(3) = atm_factor(3) - 1.0_dp/3.0_dp* &
4102 (quadrupoles(3, i, iparticle)*pw_grid%g(i, gpt) + &
4103 quadrupoles(i, 3, iparticle)*pw_grid%g(i, gpt))
4107 END SUBROUTINE get_atom_factor_stress
4128 SUBROUTINE ewald_multipole_self(ewald_env, cell, local_particles, e_self, &
4129 e_neut, task, do_efield, radii, charges, dipoles, quadrupoles, efield0, &
4134 REAL(kind=
dp),
INTENT(OUT) :: e_self, e_neut
4135 LOGICAL,
DIMENSION(3, 3),
INTENT(IN) :: task
4136 LOGICAL,
INTENT(IN) :: do_efield
4137 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: radii, charges
4138 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
4139 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
4140 POINTER :: quadrupoles
4141 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0
4142 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: efield1, efield2
4144 REAL(kind=
dp),
PARAMETER :: f23 = 2.0_dp/3.0_dp, &
4145 f415 = 4.0_dp/15.0_dp
4147 INTEGER :: ewald_type, i, iparticle, &
4148 iparticle_kind, iparticle_local, j, &
4150 LOGICAL :: do_efield0, do_efield1, do_efield2, &
4152 REAL(kind=
dp) :: alpha, ch_qu_self, ch_qu_self_tmp, &
4153 dipole_self, fac1, fac2, fac3, fac4, &
4154 q, q_neutg, q_self, q_sum, qu_qu_self, &
4158 CALL ewald_env_get(ewald_env, ewald_type=ewald_type, alpha=alpha, &
4161 do_efield0 = do_efield .AND.
ASSOCIATED(efield0)
4162 do_efield1 = do_efield .AND.
ASSOCIATED(efield1)
4163 do_efield2 = do_efield .AND.
ASSOCIATED(efield2)
4166 dipole_self = 0.0_dp
4170 fac2 = 6.0_dp*(f23**2)*(alpha**3)*
oorootpi
4171 fac3 = (2.0_dp*
oorootpi)*f23*alpha**3
4172 fac4 = (4.0_dp*
oorootpi)*f415*alpha**5
4173 lradii =
PRESENT(radii)
4176 DO iparticle_kind = 1,
SIZE(local_particles%n_el)
4177 nparticle_local = local_particles%n_el(iparticle_kind)
4178 DO iparticle_local = 1, nparticle_local
4179 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
4180 IF (any(task(1, :)))
THEN
4182 q = charges(iparticle)
4183 IF (lradii) radius = radii(iparticle)
4184 IF (radius > 0)
THEN
4185 q_neutg = q_neutg + 2.0_dp*q*radius**2
4187 q_self = q_self + q*q
4190 IF (do_efield0)
THEN
4191 efield0(iparticle) = efield0(iparticle) - q*fac1
4194 IF (task(1, 3))
THEN
4196 ch_qu_self_tmp = 0.0_dp
4198 ch_qu_self_tmp = ch_qu_self_tmp + quadrupoles(i, i, iparticle)*q
4200 ch_qu_self = ch_qu_self + ch_qu_self_tmp
4202 IF (do_efield2)
THEN
4203 efield2(1, iparticle) = efield2(1, iparticle) + fac2*q
4204 efield2(5, iparticle) = efield2(5, iparticle) + fac2*q
4205 efield2(9, iparticle) = efield2(9, iparticle) + fac2*q
4209 IF (any(task(2, :)))
THEN
4212 dipole_self = dipole_self + dipoles(i, iparticle)**2
4215 IF (do_efield1)
THEN
4216 efield1(1, iparticle) = efield1(1, iparticle) + fac3*dipoles(1, iparticle)
4217 efield1(2, iparticle) = efield1(2, iparticle) + fac3*dipoles(2, iparticle)
4218 efield1(3, iparticle) = efield1(3, iparticle) + fac3*dipoles(3, iparticle)
4221 IF (any(task(3, :)))
THEN
4225 qu_qu_self = qu_qu_self + quadrupoles(j, i, iparticle)**2
4229 IF (do_efield2)
THEN
4230 efield2(1, iparticle) = efield2(1, iparticle) + fac4*quadrupoles(1, 1, iparticle)
4231 efield2(2, iparticle) = efield2(2, iparticle) + fac4*quadrupoles(2, 1, iparticle)
4232 efield2(3, iparticle) = efield2(3, iparticle) + fac4*quadrupoles(3, 1, iparticle)
4233 efield2(4, iparticle) = efield2(4, iparticle) + fac4*quadrupoles(1, 2, iparticle)
4234 efield2(5, iparticle) = efield2(5, iparticle) + fac4*quadrupoles(2, 2, iparticle)
4235 efield2(6, iparticle) = efield2(6, iparticle) + fac4*quadrupoles(3, 2, iparticle)
4236 efield2(7, iparticle) = efield2(7, iparticle) + fac4*quadrupoles(1, 3, iparticle)
4237 efield2(8, iparticle) = efield2(8, iparticle) + fac4*quadrupoles(2, 3, iparticle)
4238 efield2(9, iparticle) = efield2(9, iparticle) + fac4*quadrupoles(3, 3, iparticle)
4244 CALL group%sum(q_neutg)
4245 CALL group%sum(q_self)
4246 CALL group%sum(q_sum)
4247 CALL group%sum(dipole_self)
4248 CALL group%sum(ch_qu_self)
4249 CALL group%sum(qu_qu_self)
4251 e_self = -(q_self + f23*(dipole_self - f23*ch_qu_self + f415*qu_qu_self*alpha**2)*alpha**2)*alpha*
oorootpi
4252 fac1 =
pi/(2.0_dp*cell%deth)
4253 e_neut = -q_sum*fac1*(q_sum/alpha**2 - q_neutg)
4256 DO iparticle_kind = 1,
SIZE(local_particles%n_el)
4257 nparticle_local = local_particles%n_el(iparticle_kind)
4258 DO iparticle_local = 1, nparticle_local
4259 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
4260 IF (any(task(1, :)))
THEN
4262 IF (do_efield0)
THEN
4263 efield0(iparticle) = efield0(iparticle) - q_sum*2.0_dp*fac1/alpha**2
4264 IF (lradii) radius = radii(iparticle)
4265 IF (radius > 0)
THEN
4266 q = charges(iparticle)
4267 efield0(iparticle) = efield0(iparticle) + fac1*radius**2*(q_sum + q)
4274 END SUBROUTINE ewald_multipole_self
4286 SUBROUTINE ewald_multipole_print(iw, e_gspace, e_rspace, e_bonded, e_self, e_neut)
4288 INTEGER,
INTENT(IN) :: iw
4289 REAL(kind=
dp),
INTENT(IN) :: e_gspace, e_rspace, e_bonded, e_self, &
4293 WRITE (iw,
'( A, A )')
' *********************************', &
4294 '**********************************************'
4295 WRITE (iw,
'( A, A, T35, A, T56, E25.15 )')
' INITIAL GSPACE ENERGY', &
4296 '[hartree]',
'= ', e_gspace
4297 WRITE (iw,
'( A, A, T35, A, T56, E25.15 )')
' INITIAL RSPACE ENERGY', &
4298 '[hartree]',
'= ', e_rspace
4299 WRITE (iw,
'( A, A, T35, A, T56, E25.15 )')
' BONDED CORRECTION', &
4300 '[hartree]',
'= ', e_bonded
4301 WRITE (iw,
'( A, A, T35, A, T56, E25.15 )')
' SELF ENERGY CORRECTION', &
4302 '[hartree]',
'= ', e_self
4303 WRITE (iw,
'( A, A, T35, A, T56, E25.15 )')
' NEUTRALIZ. BCKGR. ENERGY', &
4304 '[hartree]',
'= ', e_neut
4305 WRITE (iw,
'( A, A, T35, A, T56, E25.15 )')
' TOTAL ELECTROSTATIC EN.', &
4306 '[hartree]',
'= ', e_rspace + e_bonded + e_gspace + e_self + e_neut
4307 WRITE (iw,
'( A, A )')
' *********************************', &
4308 '**********************************************'
4310 END SUBROUTINE ewald_multipole_print
4325 SUBROUTINE debug_ewald_multipoles(ewald_env, ewald_pw, nonbond_env, cell, &
4326 particle_set, local_particles, iw, debug_r_space)
4332 POINTER :: particle_set
4334 INTEGER,
INTENT(IN) :: iw
4335 LOGICAL,
INTENT(IN) :: debug_r_space
4337 INTEGER :: nparticles
4338 LOGICAL,
DIMENSION(3) :: task
4339 REAL(kind=
dp) :: e_neut, e_self, g_energy, &
4340 r_energy, debug_energy
4341 REAL(kind=
dp),
POINTER,
DIMENSION(:) :: charges
4342 REAL(kind=
dp),
POINTER, &
4343 DIMENSION(:, :) :: dipoles, g_forces, g_pv, &
4344 r_forces, r_pv, e_field1, &
4346 REAL(kind=
dp),
POINTER, &
4347 DIMENSION(:, :, :) :: quadrupoles
4349 TYPE(multi_charge_type),
DIMENSION(:), &
4350 POINTER :: multipoles
4352 NULLIFY (multipoles, charges, dipoles, g_forces, g_pv, &
4353 r_forces, r_pv, e_field1, e_field2)
4358 nparticles =
SIZE(particle_set)
4361 ALLOCATE (charges(nparticles))
4362 ALLOCATE (dipoles(3, nparticles))
4363 ALLOCATE (quadrupoles(3, 3, nparticles))
4366 ALLOCATE (r_forces(3, nparticles))
4367 ALLOCATE (g_forces(3, nparticles))
4368 ALLOCATE (e_field1(3, nparticles))
4369 ALLOCATE (e_field2(3, nparticles))
4370 ALLOCATE (g_pv(3, 3))
4371 ALLOCATE (r_pv(3, 3))
4377 quadrupoles = 0.0_dp
4389 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2,
"CHARGE", echarge=-1.0_dp, &
4390 random_stream=random_stream, charges=charges)
4391 CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles,
"CHARGE", echarge=1.0_dp, &
4392 random_stream=random_stream, charges=charges)
4393 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
4396 WRITE (iw, *)
"DEBUG ENERGY (CHARGE-CHARGE): ", debug_energy
4398 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
4399 task, do_correction_bonded=.false., do_forces=.true., do_stress=.true., do_efield=.false., &
4400 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
4401 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.false.)
4402 CALL release_multi_type(multipoles)
4409 quadrupoles = 0.0_dp
4421 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2,
"CHARGE", echarge=-1.0_dp, &
4422 random_stream=random_stream, charges=charges)
4423 CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles,
"DIPOLE", echarge=0.5_dp, &
4424 random_stream=random_stream, dipoles=dipoles)
4425 WRITE (iw,
'("CHARGES",F15.9)') charges
4426 WRITE (iw,
'("DIPOLES",3F15.9)') dipoles
4427 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
4430 WRITE (iw, *)
"DEBUG ENERGY (CHARGE-DIPOLE): ", debug_energy
4432 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
4433 task, do_correction_bonded=.false., do_forces=.true., do_stress=.true., do_efield=.false., &
4434 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
4435 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.false.)
4436 CALL release_multi_type(multipoles)
4442 quadrupoles = 0.0_dp
4454 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2,
"DIPOLE", echarge=10000.0_dp, &
4455 random_stream=random_stream, dipoles=dipoles)
4456 CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles,
"DIPOLE", echarge=20000._dp, &
4457 random_stream=random_stream, dipoles=dipoles)
4458 WRITE (iw,
'("DIPOLES",3F15.9)') dipoles
4459 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
4462 WRITE (iw, *)
"DEBUG ENERGY (DIPOLE-DIPOLE): ", debug_energy
4464 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
4465 task, do_correction_bonded=.false., do_forces=.true., do_stress=.true., do_efield=.false., &
4466 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
4467 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.false.)
4468 CALL release_multi_type(multipoles)
4475 quadrupoles = 0.0_dp
4487 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2,
"CHARGE", echarge=-1.0_dp, &
4488 random_stream=random_stream, charges=charges)
4489 CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles,
"QUADRUPOLE", echarge=10.0_dp, &
4490 random_stream=random_stream, quadrupoles=quadrupoles)
4491 WRITE (iw,
'("CHARGES",F15.9)') charges
4492 WRITE (iw,
'("QUADRUPOLES",9F15.9)') quadrupoles
4493 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
4496 WRITE (iw, *)
"DEBUG ENERGY (CHARGE-QUADRUPOLE): ", debug_energy
4498 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
4499 task, do_correction_bonded=.false., do_forces=.true., do_stress=.true., do_efield=.false., &
4500 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
4501 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.false.)
4502 CALL release_multi_type(multipoles)
4509 quadrupoles = 0.0_dp
4521 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2,
"DIPOLE", echarge=10000.0_dp, &
4522 random_stream=random_stream, dipoles=dipoles)
4523 CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles,
"QUADRUPOLE", echarge=10000.0_dp, &
4524 random_stream=random_stream, quadrupoles=quadrupoles)
4525 WRITE (iw,
'("DIPOLES",3F15.9)') dipoles
4526 WRITE (iw,
'("QUADRUPOLES",9F15.9)') quadrupoles
4527 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
4530 WRITE (iw, *)
"DEBUG ENERGY (DIPOLE-QUADRUPOLE): ", debug_energy
4532 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
4533 task, do_correction_bonded=.false., do_forces=.true., do_stress=.true., do_efield=.false., &
4534 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
4535 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.false.)
4536 CALL release_multi_type(multipoles)
4542 quadrupoles = 0.0_dp
4554 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2,
"QUADRUPOLE", echarge=-20000.0_dp, &
4555 random_stream=random_stream, quadrupoles=quadrupoles)
4556 CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles,
"QUADRUPOLE", echarge=10000.0_dp, &
4557 random_stream=random_stream, quadrupoles=quadrupoles)
4558 WRITE (iw,
'("QUADRUPOLES",9F15.9)') quadrupoles
4559 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
4562 WRITE (iw, *)
"DEBUG ENERGY (QUADRUPOLE-QUADRUPOLE): ", debug_energy
4564 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
4565 task, do_correction_bonded=.false., do_forces=.true., do_stress=.true., do_efield=.false., &
4566 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
4567 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.false.)
4568 CALL release_multi_type(multipoles)
4570 DEALLOCATE (charges)
4571 DEALLOCATE (dipoles)
4572 DEALLOCATE (quadrupoles)
4573 DEALLOCATE (r_forces)
4574 DEALLOCATE (g_forces)
4575 DEALLOCATE (e_field1)
4576 DEALLOCATE (e_field2)
4592 SUBROUTINE debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, &
4593 energy, debug_r_space)
4597 TYPE(multi_charge_type),
DIMENSION(:),
POINTER :: multipoles
4598 REAL(kind=
dp),
INTENT(OUT) :: energy
4599 LOGICAL,
INTENT(IN) :: debug_r_space
4601 INTEGER :: atom_a, atom_b, icell, iend, igrp, &
4602 ikind, ilist, ipair, istart, jcell, &
4603 jkind, k, k1, kcell, l, l1, ncells, &
4605 INTEGER,
DIMENSION(:, :),
POINTER ::
list
4606 REAL(kind=
dp) :: fac_ij, q, r, rab2, rab2_max
4607 REAL(kind=
dp),
DIMENSION(3) :: cell_v, cvi, rab, rab0, rm
4610 TYPE(
pos_type),
DIMENSION(:),
POINTER :: r_last_update, r_last_update_pbc
4614 r_last_update=r_last_update, r_last_update_pbc=r_last_update_pbc)
4615 rab2_max = huge(0.0_dp)
4616 IF (debug_r_space)
THEN
4619 lists:
DO ilist = 1, nonbonded%nlists
4620 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
4621 npairs = neighbor_kind_pair%npairs
4622 IF (npairs == 0) cycle lists
4623 list => neighbor_kind_pair%list
4624 cvi = neighbor_kind_pair%cell_vector
4625 cell_v = matmul(cell%hmat, cvi)
4626 kind_group_loop:
DO igrp = 1, neighbor_kind_pair%ngrp_kind
4627 istart = neighbor_kind_pair%grp_kind_start(igrp)
4628 iend = neighbor_kind_pair%grp_kind_end(igrp)
4629 ikind = neighbor_kind_pair%ij_kind(1, igrp)
4630 jkind = neighbor_kind_pair%ij_kind(2, igrp)
4631 pairs:
DO ipair = istart, iend
4633 atom_a =
list(1, ipair)
4634 atom_b =
list(2, ipair)
4635 IF (atom_a == atom_b) fac_ij = 0.5_dp
4636 rab = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
4638 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
4639 IF (rab2 <= rab2_max)
THEN
4641 DO k = 1,
SIZE(multipoles(atom_a)%charge_typ)
4642 DO k1 = 1,
SIZE(multipoles(atom_a)%charge_typ(k)%charge)
4644 DO l = 1,
SIZE(multipoles(atom_b)%charge_typ)
4645 DO l1 = 1,
SIZE(multipoles(atom_b)%charge_typ(l)%charge)
4647 rm = rab + multipoles(atom_b)%charge_typ(l)%pos(:, l1) - multipoles(atom_a)%charge_typ(k)%pos(:, k1)
4649 q = multipoles(atom_b)%charge_typ(l)%charge(l1)*multipoles(atom_a)%charge_typ(k)%charge(k1)
4650 energy = energy + q/r*fac_ij
4659 END DO kind_group_loop
4665 DO atom_a = 1,
SIZE(particle_set)
4666 DO atom_b = atom_a,
SIZE(particle_set)
4668 IF (atom_a == atom_b) fac_ij = 0.5_dp
4669 rab0 = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
4671 DO icell = -ncells, ncells
4672 DO jcell = -ncells, ncells
4673 DO kcell = -ncells, ncells
4674 cell_v = matmul(cell%hmat, real([icell, jcell, kcell], kind=
dp))
4675 IF (all(cell_v == 0.0_dp) .AND. (atom_a == atom_b)) cycle
4677 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
4678 IF (rab2 <= rab2_max)
THEN
4680 DO k = 1,
SIZE(multipoles(atom_a)%charge_typ)
4681 DO k1 = 1,
SIZE(multipoles(atom_a)%charge_typ(k)%charge)
4683 DO l = 1,
SIZE(multipoles(atom_b)%charge_typ)
4684 DO l1 = 1,
SIZE(multipoles(atom_b)%charge_typ(l)%charge)
4686 rm = rab + multipoles(atom_b)%charge_typ(l)%pos(:, l1) - multipoles(atom_a)%charge_typ(k)%pos(:, k1)
4688 q = multipoles(atom_b)%charge_typ(l)%charge(l1)*multipoles(atom_a)%charge_typ(k)%charge(k1)
4689 energy = energy + q/r*fac_ij
4703 END SUBROUTINE debug_ewald_multipole_low
4720 SUBROUTINE create_multi_type(multipoles, idim, istart, iend, label, echarge, &
4721 random_stream, charges, dipoles, quadrupoles)
4722 TYPE(multi_charge_type),
DIMENSION(:),
POINTER :: multipoles
4723 INTEGER,
INTENT(IN) :: idim, istart, iend
4724 CHARACTER(LEN=*),
INTENT(IN) :: label
4725 REAL(kind=
dp),
INTENT(IN) :: echarge
4727 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: charges
4728 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
4729 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
4730 POINTER :: quadrupoles
4732 INTEGER :: i, isize, k, l, m
4733 REAL(kind=
dp) :: dx, r2, rvec(3), rvec1(3), rvec2(3)
4735 IF (
ASSOCIATED(multipoles))
THEN
4736 cpassert(
SIZE(multipoles) == idim)
4738 ALLOCATE (multipoles(idim))
4740 NULLIFY (multipoles(i)%charge_typ)
4744 IF (
ASSOCIATED(multipoles(i)%charge_typ))
THEN
4746 isize =
SIZE(multipoles(i)%charge_typ) + 1
4750 CALL reallocate_charge_type(multipoles(i)%charge_typ, 1, isize)
4753 cpassert(
PRESENT(charges))
4754 cpassert(
ASSOCIATED(charges))
4755 ALLOCATE (multipoles(i)%charge_typ(isize)%charge(1))
4756 ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 1))
4758 multipoles(i)%charge_typ(isize)%charge(1) = echarge
4759 multipoles(i)%charge_typ(isize)%pos(1:3, 1) = 0.0_dp
4760 charges(i) = charges(i) + echarge
4763 cpassert(
PRESENT(dipoles))
4764 cpassert(
ASSOCIATED(dipoles))
4765 ALLOCATE (multipoles(i)%charge_typ(isize)%charge(2))
4766 ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 2))
4767 CALL random_stream%fill(rvec)
4768 rvec = rvec/(2.0_dp*norm2(rvec))*dx
4769 multipoles(i)%charge_typ(isize)%charge(1) = echarge
4770 multipoles(i)%charge_typ(isize)%pos(1:3, 1) = rvec
4771 multipoles(i)%charge_typ(isize)%charge(2) = -echarge
4772 multipoles(i)%charge_typ(isize)%pos(1:3, 2) = -rvec
4774 dipoles(:, i) = dipoles(:, i) + 2.0_dp*echarge*rvec
4777 cpassert(
PRESENT(quadrupoles))
4778 cpassert(
ASSOCIATED(quadrupoles))
4779 ALLOCATE (multipoles(i)%charge_typ(isize)%charge(4))
4780 ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 4))
4781 CALL random_stream%fill(rvec1)
4782 CALL random_stream%fill(rvec2)
4783 rvec1 = rvec1/norm2(rvec1)
4784 rvec2 = rvec2 - dot_product(rvec2, rvec1)*rvec1
4785 rvec2 = rvec2/norm2(rvec2)
4787 rvec1 = rvec1/2.0_dp*dx
4788 rvec2 = rvec2/2.0_dp*dx
4796 multipoles(i)%charge_typ(isize)%charge(1) = -echarge
4797 multipoles(i)%charge_typ(isize)%pos(1:3, 1) = rvec1 + rvec2
4798 multipoles(i)%charge_typ(isize)%charge(2) = echarge
4799 multipoles(i)%charge_typ(isize)%pos(1:3, 2) = rvec1 - rvec2
4800 multipoles(i)%charge_typ(isize)%charge(3) = -echarge
4801 multipoles(i)%charge_typ(isize)%pos(1:3, 3) = -rvec1 - rvec2
4802 multipoles(i)%charge_typ(isize)%charge(4) = echarge
4803 multipoles(i)%charge_typ(isize)%pos(1:3, 4) = -rvec1 + rvec2
4806 r2 = dot_product(multipoles(i)%charge_typ(isize)%pos(:, k), multipoles(i)%charge_typ(isize)%pos(:, k))
4809 quadrupoles(m, l, i) = quadrupoles(m, l, i) + 3.0_dp*0.5_dp*multipoles(i)%charge_typ(isize)%charge(k)* &
4810 multipoles(i)%charge_typ(isize)%pos(l, k)* &
4811 multipoles(i)%charge_typ(isize)%pos(m, k)
4812 IF (m == l) quadrupoles(m, l, i) = quadrupoles(m, l, i) - 0.5_dp*multipoles(i)%charge_typ(isize)%charge(k)*r2
4819 END SUBROUTINE create_multi_type
4827 SUBROUTINE release_multi_type(multipoles)
4828 TYPE(multi_charge_type),
DIMENSION(:),
POINTER :: multipoles
4832 IF (
ASSOCIATED(multipoles))
THEN
4833 DO i = 1,
SIZE(multipoles)
4834 DO j = 1,
SIZE(multipoles(i)%charge_typ)
4835 DEALLOCATE (multipoles(i)%charge_typ(j)%charge)
4836 DEALLOCATE (multipoles(i)%charge_typ(j)%pos)
4838 DEALLOCATE (multipoles(i)%charge_typ)
4841 END SUBROUTINE release_multi_type
4851 SUBROUTINE reallocate_charge_type(charge_typ, istart, iend)
4852 TYPE(charge_mono_type),
DIMENSION(:),
POINTER :: charge_typ
4853 INTEGER,
INTENT(IN) :: istart, iend
4855 INTEGER :: i, isize, j, jsize, jsize1, jsize2
4856 TYPE(charge_mono_type),
DIMENSION(:),
POINTER :: charge_typ_bk
4858 IF (
ASSOCIATED(charge_typ))
THEN
4859 isize =
SIZE(charge_typ)
4860 ALLOCATE (charge_typ_bk(1:isize))
4862 jsize =
SIZE(charge_typ(j)%charge)
4863 ALLOCATE (charge_typ_bk(j)%charge(jsize))
4864 jsize1 =
SIZE(charge_typ(j)%pos, 1)
4865 jsize2 =
SIZE(charge_typ(j)%pos, 2)
4866 ALLOCATE (charge_typ_bk(j)%pos(jsize1, jsize2))
4867 charge_typ_bk(j)%pos = charge_typ(j)%pos
4868 charge_typ_bk(j)%charge = charge_typ(j)%charge
4870 DO j = 1,
SIZE(charge_typ)
4871 DEALLOCATE (charge_typ(j)%charge)
4872 DEALLOCATE (charge_typ(j)%pos)
4874 DEALLOCATE (charge_typ)
4876 ALLOCATE (charge_typ_bk(istart:iend))
4877 DO i = istart, isize
4878 jsize =
SIZE(charge_typ_bk(j)%charge)
4879 ALLOCATE (charge_typ(j)%charge(jsize))
4880 jsize1 =
SIZE(charge_typ_bk(j)%pos, 1)
4881 jsize2 =
SIZE(charge_typ_bk(j)%pos, 2)
4882 ALLOCATE (charge_typ(j)%pos(jsize1, jsize2))
4883 charge_typ(j)%pos = charge_typ_bk(j)%pos
4884 charge_typ(j)%charge = charge_typ_bk(j)%charge
4886 DO j = 1,
SIZE(charge_typ_bk)
4887 DEALLOCATE (charge_typ_bk(j)%charge)
4888 DEALLOCATE (charge_typ_bk(j)%pos)
4890 DEALLOCATE (charge_typ_bk)
4892 ALLOCATE (charge_typ(istart:iend))
4895 END SUBROUTINE reallocate_charge_type
4897 END SUBROUTINE debug_ewald_multipoles
4918 SUBROUTINE debug_ewald_multipoles_fields(ewald_env, ewald_pw, nonbond_env, cell, &
4919 particle_set, local_particles, radii, charges, dipoles, quadrupoles, task, iw, &
4920 atomic_kind_set, mm_section)
4927 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: radii, charges
4928 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
4929 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
4930 POINTER :: quadrupoles
4931 LOGICAL,
DIMENSION(3),
INTENT(IN) :: task
4932 INTEGER,
INTENT(IN) :: iw
4936 INTEGER :: i, iparticle_kind, j, k, &
4937 nparticle_local, nparticles
4938 REAL(kind=
dp) :: coord(3), dq, e_neut, e_self, efield1n(3), efield2n(3, 3), ene(2), &
4939 energy_glob, energy_local, enev(3, 2), o_tot_ene, pot, pv_glob(3, 3), pv_local(3, 3), &
4941 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: efield1, efield2, forces_glob, &
4943 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0, lcharges
4945 TYPE(
particle_type),
DIMENSION(:),
POINTER :: core_particle_set, shell_particle_set
4947 NULLIFY (lcharges, shell_particle_set, core_particle_set)
4951 nparticles =
SIZE(particle_set)
4953 DO iparticle_kind = 1,
SIZE(local_particles%n_el)
4954 nparticle_local = nparticle_local + local_particles%n_el(iparticle_kind)
4956 ALLOCATE (lcharges(nparticles))
4957 ALLOCATE (forces_glob(3, nparticles))
4958 ALLOCATE (forces_local(3, nparticle_local))
4959 ALLOCATE (efield0(nparticles))
4960 ALLOCATE (efield1(3, nparticles))
4961 ALLOCATE (efield2(9, nparticles))
4962 forces_glob = 0.0_dp
4963 forces_local = 0.0_dp
4969 energy_glob = 0.0_dp
4970 energy_local = 0.0_dp
4974 local_particles, energy_local, energy_glob, e_neut, e_self, task, .false., .true., .true., &
4975 .true., radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob, &
4976 efield0, efield1, efield2, iw, do_debug=.false.)
4977 o_tot_ene = energy_local + energy_glob + e_neut + e_self
4978 WRITE (iw, *)
"TOTAL ENERGY :: ========>", o_tot_ene
4982 DO i = 1, nparticles
4985 lcharges(i) = charges(i) + (-1.0_dp)**k*dq
4986 forces_glob = 0.0_dp
4987 forces_local = 0.0_dp
4990 energy_glob = 0.0_dp
4991 energy_local = 0.0_dp
4995 local_particles, energy_local, energy_glob, e_neut, e_self, &
4996 task, .false., .false., .false., .false., radii, &
4997 lcharges, dipoles, quadrupoles, iw=iw, do_debug=.false.)
4998 ene(k) = energy_local + energy_glob + e_neut + e_self
5000 pot = (ene(2) - ene(1))/(2.0_dp*dq)
5001 WRITE (iw,
'(A,I8,3(A,F15.9))')
"POTENTIAL FOR ATOM: ", i,
" NUMERICAL: ", pot,
" ANALYTICAL: ", efield0(i), &
5002 " ERROR: ", pot - efield0(i)
5003 tot_ene = tot_ene + 0.5_dp*efield0(i)*charges(i)
5005 WRITE (iw, *)
"ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
5006 WRITE (iw,
'(/,/,/)')
5009 DO i = 1, nparticles
5010 coord = particle_set(i)%r
5013 particle_set(i)%r(j) = coord(j) + (-1.0_dp)**k*dq
5016 CALL list_control(atomic_kind_set, particle_set, local_particles, &
5017 cell, nonbond_env, logger%para_env, mm_section, &
5018 shell_particle_set, core_particle_set)
5020 forces_glob = 0.0_dp
5021 forces_local = 0.0_dp
5024 energy_glob = 0.0_dp
5025 energy_local = 0.0_dp
5030 local_particles, energy_local, energy_glob, e_neut, e_self, &
5031 task, .false., .true., .true., .true., radii, &
5032 charges, dipoles, quadrupoles, forces_local, forces_glob, &
5033 pv_local, pv_glob, efield0, iw=iw, do_debug=.false.)
5035 particle_set(i)%r(j) = coord(j)
5037 efield1n(j) = -(ene(2) - ene(1))/(2.0_dp*dq)
5039 WRITE (iw,
'(/,A,I8)')
"FIELD FOR ATOM: ", i
5040 WRITE (iw,
'(A,3F15.9)')
" NUMERICAL: ", efield1n,
" ANALYTICAL: ", efield1(:, i), &
5041 " ERROR: ", efield1n - efield1(:, i)
5043 tot_ene = tot_ene - 0.5_dp*dot_product(efield1(:, i), dipoles(:, i))
5046 WRITE (iw, *)
"ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
5050 DO i = 1, nparticles
5051 coord = particle_set(i)%r
5054 particle_set(i)%r(j) = coord(j) + (-1.0_dp)**k*dq
5057 CALL list_control(atomic_kind_set, particle_set, local_particles, &
5058 cell, nonbond_env, logger%para_env, mm_section, &
5059 shell_particle_set, core_particle_set)
5061 forces_glob = 0.0_dp
5062 forces_local = 0.0_dp
5065 energy_glob = 0.0_dp
5066 energy_local = 0.0_dp
5071 local_particles, energy_local, energy_glob, e_neut, e_self, &
5072 task, .false., .true., .true., .true., radii, &
5073 charges, dipoles, quadrupoles, forces_local, forces_glob, &
5074 pv_local, pv_glob, efield1=efield1, iw=iw, do_debug=.false.)
5075 enev(:, k) = efield1(:, i)
5076 particle_set(i)%r(j) = coord(j)
5078 efield2n(:, j) = (enev(:, 2) - enev(:, 1))/(2.0_dp*dq)
5080 WRITE (iw,
'(/,A,I8)')
"FIELD GRADIENT FOR ATOM: ", i
5081 WRITE (iw,
'(A,9F15.9)')
" NUMERICAL: ", efield2n, &
5082 " ANALYTICAL: ", efield2(:, i), &
5083 " ERROR: ", reshape(efield2n, [9]) - efield2(:, i)
5085 END SUBROUTINE debug_ewald_multipoles_fields
5104 SUBROUTINE debug_ewald_multipoles_fields2(ewald_env, ewald_pw, nonbond_env, cell, &
5105 particle_set, local_particles, radii, charges, dipoles, quadrupoles, task, iw)
5112 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: radii, charges
5113 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
5114 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
5115 POINTER :: quadrupoles
5116 LOGICAL,
DIMENSION(3),
INTENT(IN) :: task
5117 INTEGER,
INTENT(IN) :: iw
5119 INTEGER :: i, ind, iparticle_kind, j, k, &
5120 nparticle_local, nparticles
5121 REAL(kind=
dp) :: e_neut, e_self, energy_glob, &
5122 energy_local, o_tot_ene, prod, &
5123 pv_glob(3, 3), pv_local(3, 3), tot_ene
5124 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: efield1, efield2, forces_glob, &
5126 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0
5132 nparticles =
SIZE(particle_set)
5134 DO iparticle_kind = 1,
SIZE(local_particles%n_el)
5135 nparticle_local = nparticle_local + local_particles%n_el(iparticle_kind)
5137 ALLOCATE (forces_glob(3, nparticles))
5138 ALLOCATE (forces_local(3, nparticle_local))
5139 ALLOCATE (efield0(nparticles))
5140 ALLOCATE (efield1(3, nparticles))
5141 ALLOCATE (efield2(9, nparticles))
5142 forces_glob = 0.0_dp
5143 forces_local = 0.0_dp
5149 energy_glob = 0.0_dp
5150 energy_local = 0.0_dp
5154 local_particles, energy_local, energy_glob, e_neut, e_self, task, .false., .true., .true., &
5155 .true., radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob, &
5156 efield0, efield1, efield2, iw, do_debug=.false.)
5157 o_tot_ene = energy_local + energy_glob + e_neut + e_self
5158 WRITE (iw, *)
"TOTAL ENERGY :: ========>", o_tot_ene
5163 DO i = 1, nparticles
5164 tot_ene = tot_ene + 0.5_dp*efield0(i)*charges(i)
5166 WRITE (iw, *)
"ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
5167 WRITE (iw,
'(/,/,/)')
5172 DO i = 1, nparticles
5173 tot_ene = tot_ene - 0.5_dp*dot_product(efield1(:, i), dipoles(:, i))
5175 WRITE (iw, *)
"ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
5176 WRITE (iw,
'(/,/,/)')
5181 DO i = 1, nparticles
5187 prod = prod + efield2(ind, i)*quadrupoles(j, k, i)
5190 tot_ene = tot_ene - 0.5_dp*(1.0_dp/3.0_dp)*prod
5192 WRITE (iw, *)
"ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
5193 WRITE (iw,
'(/,/,/)')
5196 END SUBROUTINE debug_ewald_multipoles_fields2
Define the atomic kind types and their sub types.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public aguado2003
integer, save, public laino2008
Handles all functions related to the CELL.
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
integer, parameter, public tang_toennies
integer, parameter, public no_damping
subroutine, public dg_get(dg, dg_rho0)
Get the dg_type.
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
subroutine, public ewald_pw_get(ewald_pw, pw_big_pool, pw_small_pool, rs_desc, poisson_env, dg)
get the ewald_pw environment to the correct program.
Treats the electrostatic for multipoles (up to quadrupoles)
recursive subroutine, public ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, task, do_correction_bonded, do_forces, do_stress, do_efield, radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob, efield0, efield1, efield2, iw, do_debug, atomic_kind_set, mm_section)
Computes the potential and the force for a lattice sum of multipoles (up to quadrupole)
subroutine, public list_control(atomic_kind_set, particle_set, local_particles, cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, core_particle_set, force_update, exclusions)
...
Define the neighbor list data types and the corresponding functionality.
subroutine, public fist_nonbond_env_get(fist_nonbond_env, potparm14, potparm, nonbonded, rlist_cut, rlist_lowsq, aup, lup, ei_scale14, vdw_scale14, shift_cutoff, do_electrostatics, r_last_update, r_last_update_pbc, rshell_last_update_pbc, rcore_last_update_pbc, cell_last_update, num_update, last_update, counter, natom_types, long_range_correction, ij_kind_full_fac, eam_data, nequip_data, deepmd_data, ace_data, charges)
sets a fist_nonbond_env
Defines the basic variable types.
integer, parameter, public dp
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition of mathematical constants and functions.
real(kind=dp), parameter, public oorootpi
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public sqrthalf
real(kind=dp), parameter, public fourpi
real(kind=dp), dimension(0:maxfac), parameter, public fac
complex(kind=dp), parameter, public z_zero
Interface to the message passing library MPI.
Parallel (pseudo)random number generator (RNG) for multiple streams and substreams of random numbers.
integer, parameter, public uniform
Define the data structure for the particle information.
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
subroutine, public structure_factor_deallocate(exp_igr)
...
subroutine, public structure_factor_allocate(bds, nparts, exp_igr, allocate_centre, allocate_shell_e, allocate_shell_centre, nshell)
...
subroutine, public structure_factor_evaluate(delta, lb, ex, ey, ez)
...
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Type for Gaussian Densities type = type of gaussian (PME) grid = grid number gcc = Gaussian contracti...
structure to store local (to a processor) ordered lists of integers.
to build arrays of pointers
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...