57#include "./base/base_uses.f90"
63 LOGICAL,
PRIVATE,
PARAMETER :: debug_this_module = .false.
64 LOGICAL,
PRIVATE,
PARAMETER :: debug_r_space = .false.
65 LOGICAL,
PRIVATE,
PARAMETER :: debug_g_space = .false.
66 LOGICAL,
PRIVATE,
PARAMETER :: debug_e_field = .false.
67 LOGICAL,
PRIVATE,
PARAMETER :: debug_e_field_en = .false.
68 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'ewalds_multipole'
118 cell, particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, &
119 task, do_correction_bonded, do_forces, do_stress, &
120 do_efield, radii, charges, dipoles, &
121 quadrupoles, forces_local, forces_glob, pv_local, pv_glob, efield0, efield1, &
122 efield2, iw, do_debug, atomic_kind_set, mm_section)
129 REAL(kind=
dp),
INTENT(INOUT) :: energy_local, energy_glob
130 REAL(kind=
dp),
INTENT(OUT) :: e_neut, e_self
131 LOGICAL,
DIMENSION(3),
INTENT(IN) :: task
132 LOGICAL,
INTENT(IN) :: do_correction_bonded, do_forces, &
134 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: radii, charges
135 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
136 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
137 POINTER :: quadrupoles
138 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT), &
139 OPTIONAL :: forces_local, forces_glob, pv_local, &
141 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT),
OPTIONAL :: efield0
142 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT), &
143 OPTIONAL :: efield1, efield2
144 INTEGER,
INTENT(IN) :: iw
145 LOGICAL,
INTENT(IN) :: do_debug
147 POINTER :: atomic_kind_set
150 CHARACTER(len=*),
PARAMETER :: routinen =
'ewald_multipole_evaluate'
152 INTEGER :: handle, i, j, size1, size2
153 LOGICAL :: check_debug, check_efield, check_forces, &
155 LOGICAL,
DIMENSION(3, 3) :: my_task
156 REAL(kind=
dp) :: e_bonded, e_bonded_t, e_rspace, &
157 e_rspace_t, energy_glob_t
158 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0_lr, efield0_sr
159 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: efield1_lr, efield1_sr, efield2_lr, &
165 CALL timeset(routinen, handle)
166 cpassert(
ASSOCIATED(nonbond_env))
167 check_debug = (debug_this_module .OR. debug_r_space .OR. debug_g_space .OR. debug_e_field .OR. debug_e_field_en) &
168 .EQV. debug_this_module
169 cpassert(check_debug)
170 check_forces = do_forces .EQV. (
PRESENT(forces_local) .AND.
PRESENT(forces_glob))
171 cpassert(check_forces)
172 check_efield = do_efield .EQV. (
PRESENT(efield0) .OR.
PRESENT(efield1) .OR.
PRESENT(efield2))
173 cpassert(check_efield)
175 IF (debug_this_module .AND. do_debug)
THEN
177 IF (debug_r_space)
THEN
178 CALL debug_ewald_multipoles(ewald_env, ewald_pw, nonbond_env, cell, &
179 particle_set, local_particles, iw, debug_r_space)
180 cpabort(
"Debug Multipole Requested: Real Part!")
183 IF (debug_e_field)
THEN
184 cpassert(
PRESENT(atomic_kind_set))
185 cpassert(
PRESENT(mm_section))
186 CALL debug_ewald_multipoles_fields(ewald_env, ewald_pw, nonbond_env, &
187 cell, particle_set, local_particles, radii, charges, dipoles, &
188 quadrupoles, task, iw, atomic_kind_set, mm_section)
189 cpabort(
"Debug Multipole Requested: POT+EFIELDS+GRAD!")
193 IF (debug_e_field_en)
THEN
194 CALL debug_ewald_multipoles_fields2(ewald_env, ewald_pw, nonbond_env, &
195 cell, particle_set, local_particles, radii, charges, dipoles, &
196 quadrupoles, task, iw)
197 cpabort(
"Debug Multipole Requested: POT+EFIELDS+GRAD to give the correct energy!")
207 do_task(1) = any(charges /= 0.0_dp)
209 do_task(2) = any(dipoles /= 0.0_dp)
211 do_task(3) = any(quadrupoles /= 0.0_dp)
217 my_task(j, i) = do_task(i) .AND. do_task(j)
218 my_task(i, j) = my_task(j, i)
223 NULLIFY (efield0_sr, efield0_lr, efield1_sr, efield1_lr, efield2_sr, efield2_lr)
225 IF (
PRESENT(efield0))
THEN
226 size1 =
SIZE(efield0)
227 ALLOCATE (efield0_sr(size1))
228 ALLOCATE (efield0_lr(size1))
232 IF (
PRESENT(efield1))
THEN
233 size1 =
SIZE(efield1, 1)
234 size2 =
SIZE(efield1, 2)
235 ALLOCATE (efield1_sr(size1, size2))
236 ALLOCATE (efield1_lr(size1, size2))
240 IF (
PRESENT(efield2))
THEN
241 size1 =
SIZE(efield2, 1)
242 size2 =
SIZE(efield2, 2)
243 ALLOCATE (efield2_sr(size1, size2))
244 ALLOCATE (efield2_lr(size1, size2))
252 IF ((.NOT. debug_g_space) .AND. (nonbond_env%do_nonbonded))
THEN
256 CALL ewald_multipole_sr(nonbond_env, ewald_env, atomic_kind_set, &
257 particle_set, cell, e_rspace, my_task, &
258 do_forces, do_efield, do_stress, radii, charges, dipoles, quadrupoles, &
259 forces_glob, pv_glob, efield0_sr, efield1_sr, efield2_sr)
260 energy_glob = energy_glob + e_rspace
262 IF (do_correction_bonded)
THEN
265 CALL ewald_multipole_bonded(nonbond_env, particle_set, ewald_env, &
266 cell, e_bonded, my_task, do_forces, do_efield, do_stress, &
267 charges, dipoles, quadrupoles, forces_glob, pv_glob, &
268 efield0_sr, efield1_sr, efield2_sr)
269 energy_glob = energy_glob + e_bonded
275 energy_local = 0.0_dp
276 IF (.NOT. debug_r_space)
THEN
278 CALL ewald_multipole_lr(ewald_env, ewald_pw, cell, particle_set, &
279 local_particles, energy_local, my_task, do_forces, do_efield, do_stress, &
280 charges, dipoles, quadrupoles, forces_local, pv_local, efield0_lr, efield1_lr, &
284 CALL ewald_multipole_self(ewald_env, cell, local_particles, e_self, &
285 e_neut, my_task, do_efield, radii, charges, dipoles, quadrupoles, &
286 efield0_lr, efield1_lr, efield2_lr)
291 energy_glob_t = energy_glob
292 e_rspace_t = e_rspace
293 e_bonded_t = e_bonded
294 CALL group%sum(energy_glob_t)
295 CALL group%sum(e_rspace_t)
296 CALL group%sum(e_bonded_t)
298 CALL ewald_multipole_print(iw, energy_local, e_rspace_t, e_bonded_t, e_self, e_neut)
302 IF (
PRESENT(efield0))
THEN
303 efield0 = efield0_sr + efield0_lr
304 CALL group%sum(efield0)
305 DEALLOCATE (efield0_sr)
306 DEALLOCATE (efield0_lr)
308 IF (
PRESENT(efield1))
THEN
309 efield1 = efield1_sr + efield1_lr
310 CALL group%sum(efield1)
311 DEALLOCATE (efield1_sr)
312 DEALLOCATE (efield1_lr)
314 IF (
PRESENT(efield2))
THEN
315 efield2 = efield2_sr + efield2_lr
316 CALL group%sum(efield2)
317 DEALLOCATE (efield2_sr)
318 DEALLOCATE (efield2_lr)
321 CALL timestop(handle)
348 SUBROUTINE ewald_multipole_sr(nonbond_env, ewald_env, atomic_kind_set, &
349 particle_set, cell, energy, task, &
350 do_forces, do_efield, do_stress, radii, charges, dipoles, quadrupoles, &
351 forces, pv, efield0, efield1, efield2)
355 POINTER :: atomic_kind_set
358 REAL(kind=
dp),
INTENT(INOUT) :: energy
359 LOGICAL,
DIMENSION(3, 3),
INTENT(IN) :: task
360 LOGICAL,
INTENT(IN) :: do_forces, do_efield, do_stress
361 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: radii, charges
362 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
363 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
364 POINTER :: quadrupoles
365 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT), &
366 OPTIONAL :: forces, pv
367 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0
368 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: efield1, efield2
370 CHARACTER(len=*),
PARAMETER :: routinen =
'ewald_multipole_SR'
372 INTEGER :: a, atom_a, atom_b, b, c, d, e, handle, i, iend, igrp, ikind, ilist, ipair, &
373 istart, itype_ij, itype_ji, jkind, k, kind_a, kind_b, kk, nkdamp_ij, nkdamp_ji, nkinds, &
375 INTEGER,
DIMENSION(:, :),
POINTER ::
list
376 LOGICAL :: do_efield0, do_efield1, do_efield2, &
378 REAL(kind=
dp) :: alpha, beta, ch_i, ch_j, dampa_ij, dampa_ji, dampaexpi, dampaexpj, &
379 dampfac_ij, dampfac_ji, dampfuncdiffi, dampfuncdiffj, dampfunci, dampfuncj, dampsumfi, &
380 dampsumfj, ef0_i, ef0_j, eloc,
fac, fac_ij, factorial, ir, irab2, ptens11, ptens12, &
381 ptens13, ptens21, ptens22, ptens23, ptens31, ptens32, ptens33, r, rab2, rab2_max, radius, &
382 rcut, tij, tmp, tmp1, tmp11, tmp12, tmp13, tmp2, tmp21, tmp22, tmp23, tmp31, tmp32, &
383 tmp33, tmp_ij, tmp_ji, xf
384 REAL(kind=
dp),
DIMENSION(0:5) :: f
385 REAL(kind=
dp),
DIMENSION(3) :: cell_v, cvi, damptij_a, damptji_a, dp_i, &
386 dp_j, ef1_i, ef1_j, fr, rab, tij_a
387 REAL(kind=
dp),
DIMENSION(3, 3) :: damptij_ab, damptji_ab, ef2_i, ef2_j, &
389 REAL(kind=
dp),
DIMENSION(3, 3, 3) :: tij_abc
390 REAL(kind=
dp),
DIMENSION(3, 3, 3, 3) :: tij_abcd
391 REAL(kind=
dp),
DIMENSION(3, 3, 3, 3, 3) :: tij_abcde
395 TYPE(
pos_type),
DIMENSION(:),
POINTER :: r_last_update, r_last_update_pbc
397 CALL timeset(routinen, handle)
398 NULLIFY (nonbonded, r_last_update, r_last_update_pbc)
399 do_efield0 = do_efield .AND.
ASSOCIATED(efield0)
400 do_efield1 = do_efield .AND.
ASSOCIATED(efield1)
401 do_efield2 = do_efield .AND.
ASSOCIATED(efield2)
403 ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp
404 ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp
405 ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp
409 r_last_update=r_last_update, r_last_update_pbc=r_last_update_pbc)
412 IF (debug_r_space)
THEN
413 rab2_max = huge(0.0_dp)
416 lists:
DO ilist = 1, nonbonded%nlists
417 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
418 npairs = neighbor_kind_pair%npairs
419 IF (npairs == 0) cycle
420 list => neighbor_kind_pair%list
421 cvi = neighbor_kind_pair%cell_vector
422 cell_v = matmul(cell%hmat, cvi)
423 kind_group_loop:
DO igrp = 1, neighbor_kind_pair%ngrp_kind
424 istart = neighbor_kind_pair%grp_kind_start(igrp)
425 iend = neighbor_kind_pair%grp_kind_end(igrp)
426 ikind = neighbor_kind_pair%ij_kind(1, igrp)
427 jkind = neighbor_kind_pair%ij_kind(2, igrp)
438 IF (
PRESENT(atomic_kind_set))
THEN
439 IF (
ASSOCIATED(atomic_kind_set(jkind)%damping))
THEN
440 damping_ij = atomic_kind_set(jkind)%damping%damp(ikind)
441 itype_ij = damping_ij%itype
442 nkdamp_ij = damping_ij%order
443 dampa_ij = damping_ij%bij
444 dampfac_ij = damping_ij%cij
447 IF (
ASSOCIATED(atomic_kind_set(ikind)%damping))
THEN
448 damping_ji = atomic_kind_set(ikind)%damping%damp(jkind)
449 itype_ji = damping_ji%itype
450 nkdamp_ji = damping_ji%order
451 dampa_ji = damping_ji%bij
452 dampfac_ji = damping_ji%cij
456 pairs:
DO ipair = istart, iend
457 IF (ipair <= neighbor_kind_pair%nscale)
THEN
460 fac_ij = neighbor_kind_pair%ei_scale(ipair)
461 IF (fac_ij <= 0) cycle
465 atom_a =
list(1, ipair)
466 atom_b =
list(2, ipair)
467 kind_a = particle_set(atom_a)%atomic_kind%kind_number
468 kind_b = particle_set(atom_b)%atomic_kind%kind_number
469 IF (atom_a == atom_b) fac_ij = 0.5_dp
470 rab = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
472 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
473 IF (rab2 <= rab2_max)
THEN
474 IF (
PRESENT(radii))
THEN
475 radius = sqrt(radii(atom_a)*radii(atom_a) + radii(atom_b)*radii(atom_b))
479 IF (radius > 0.0_dp)
THEN
485 tij_ab = huge(0.0_dp)
486 tij_abc = huge(0.0_dp)
487 tij_abcd = huge(0.0_dp)
488 tij_abcde = huge(0.0_dp)
495 IF (debug_this_module .AND. debug_r_space .AND. (.NOT. debug_g_space))
THEN
500 f(0) = erf(beta*r)*ir - erf(alpha*r)*ir
507 f(i) = irab2*(f(i - 1) + tmp1*((2.0_dp*alpha**2)**i)/(
fac*alpha) - tmp2*((2.0_dp*beta**2)**i)/(
fac*beta))
512 force_eval = do_stress
515 force_eval = do_forces .OR. do_efield1
517 IF (task(2, 2)) force_eval = force_eval .OR. do_efield0
518 IF (task(1, 2) .OR. force_eval)
THEN
519 force_eval = do_stress
520 tij_a = -rab*f(1)*fac_ij
521 IF (task(1, 2)) force_eval = force_eval .OR. do_forces
523 IF (task(1, 1)) force_eval = force_eval .OR. do_efield2
524 IF (task(3, 3)) force_eval = force_eval .OR. do_efield0
525 IF (task(2, 2) .OR. task(3, 1) .OR. force_eval)
THEN
526 force_eval = do_stress
529 tmp = rab(a)*rab(b)*fac_ij
530 tij_ab(a, b) = 3.0_dp*tmp*f(2)
531 IF (a == b) tij_ab(a, b) = tij_ab(a, b) - f(1)*fac_ij
534 IF (task(2, 2) .OR. task(3, 1)) force_eval = force_eval .OR. do_forces
536 IF (task(2, 2)) force_eval = force_eval .OR. do_efield2
537 IF (task(3, 3)) force_eval = force_eval .OR. do_efield1
538 IF (task(3, 2) .OR. force_eval)
THEN
539 force_eval = do_stress
543 tmp = rab(a)*rab(b)*rab(c)*fac_ij
544 tij_abc(a, b, c) = -15.0_dp*tmp*f(3)
545 tmp = 3.0_dp*f(2)*fac_ij
546 IF (a == b) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(c)
547 IF (a == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(b)
548 IF (b == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(a)
552 IF (task(3, 2)) force_eval = force_eval .OR. do_forces
554 IF (task(3, 3) .OR. force_eval)
THEN
555 force_eval = do_stress
560 tmp = rab(a)*rab(b)*rab(c)*rab(d)*fac_ij
561 tij_abcd(a, b, c, d) = 105.0_dp*tmp*f(4)
562 tmp1 = 15.0_dp*f(3)*fac_ij
563 tmp2 = 3.0_dp*f(2)*fac_ij
565 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(c)*rab(d)
566 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
569 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(d)
570 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
572 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(c)
574 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(d)
575 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
577 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(c)
578 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(b)
583 IF (task(3, 3)) force_eval = force_eval .OR. do_forces
586 force_eval = do_stress
592 tmp = rab(a)*rab(b)*rab(c)*rab(d)*rab(e)*fac_ij
593 tij_abcde(a, b, c, d, e) = -945.0_dp*tmp*f(5)
594 tmp1 = 105.0_dp*f(4)*fac_ij
595 tmp2 = 15.0_dp*f(3)*fac_ij
597 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(c)*rab(d)*rab(e)
598 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
599 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
600 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
603 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(d)*rab(e)
604 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
605 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
606 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
609 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(e)
610 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
611 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
612 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
615 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(d)
616 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
617 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
618 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
621 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(d)*rab(e)
622 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
625 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(e)
626 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
629 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(d)
630 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
632 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(e)
633 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(d)
634 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(c)
652 IF (debug_this_module)
THEN
660 IF (any(task(1, :)))
THEN
661 ch_j = charges(atom_a)
662 ch_i = charges(atom_b)
664 IF (any(task(2, :)))
THEN
665 dp_j = dipoles(:, atom_a)
666 dp_i = dipoles(:, atom_b)
668 IF (any(task(3, :)))
THEN
669 qp_j = quadrupoles(:, :, atom_a)
670 qp_i = quadrupoles(:, :, atom_b)
674 eloc = eloc + ch_i*tij*ch_j
676 IF (do_forces .OR. do_stress)
THEN
677 fr(1) = fr(1) - ch_j*tij_a(1)*ch_i
678 fr(2) = fr(2) - ch_j*tij_a(2)*ch_i
679 fr(3) = fr(3) - ch_j*tij_a(3)*ch_i
685 ef0_i = ef0_i + tij*ch_j
687 ef0_j = ef0_j + tij*ch_i
691 ef1_i(1) = ef1_i(1) - tij_a(1)*ch_j
692 ef1_i(2) = ef1_i(2) - tij_a(2)*ch_j
693 ef1_i(3) = ef1_i(3) - tij_a(3)*ch_j
695 ef1_j(1) = ef1_j(1) + tij_a(1)*ch_i
696 ef1_j(2) = ef1_j(2) + tij_a(2)*ch_i
697 ef1_j(3) = ef1_j(3) + tij_a(3)*ch_i
703 ef2_i(1, 1) = ef2_i(1, 1) - tij_ab(1, 1)*ch_j
704 ef2_i(2, 1) = ef2_i(2, 1) - tij_ab(2, 1)*ch_j
705 ef2_i(3, 1) = ef2_i(3, 1) - tij_ab(3, 1)*ch_j
706 ef2_i(1, 2) = ef2_i(1, 2) - tij_ab(1, 2)*ch_j
707 ef2_i(2, 2) = ef2_i(2, 2) - tij_ab(2, 2)*ch_j
708 ef2_i(3, 2) = ef2_i(3, 2) - tij_ab(3, 2)*ch_j
709 ef2_i(1, 3) = ef2_i(1, 3) - tij_ab(1, 3)*ch_j
710 ef2_i(2, 3) = ef2_i(2, 3) - tij_ab(2, 3)*ch_j
711 ef2_i(3, 3) = ef2_i(3, 3) - tij_ab(3, 3)*ch_j
713 ef2_j(1, 1) = ef2_j(1, 1) - tij_ab(1, 1)*ch_i
714 ef2_j(2, 1) = ef2_j(2, 1) - tij_ab(2, 1)*ch_i
715 ef2_j(3, 1) = ef2_j(3, 1) - tij_ab(3, 1)*ch_i
716 ef2_j(1, 2) = ef2_j(1, 2) - tij_ab(1, 2)*ch_i
717 ef2_j(2, 2) = ef2_j(2, 2) - tij_ab(2, 2)*ch_i
718 ef2_j(3, 2) = ef2_j(3, 2) - tij_ab(3, 2)*ch_i
719 ef2_j(1, 3) = ef2_j(1, 3) - tij_ab(1, 3)*ch_i
720 ef2_j(2, 3) = ef2_j(2, 3) - tij_ab(2, 3)*ch_i
721 ef2_j(3, 3) = ef2_j(3, 3) - tij_ab(3, 3)*ch_i
727 tmp = -(dp_i(1)*(tij_ab(1, 1)*dp_j(1) + &
728 tij_ab(2, 1)*dp_j(2) + &
729 tij_ab(3, 1)*dp_j(3)) + &
730 dp_i(2)*(tij_ab(1, 2)*dp_j(1) + &
731 tij_ab(2, 2)*dp_j(2) + &
732 tij_ab(3, 2)*dp_j(3)) + &
733 dp_i(3)*(tij_ab(1, 3)*dp_j(1) + &
734 tij_ab(2, 3)*dp_j(2) + &
735 tij_ab(3, 3)*dp_j(3)))
738 IF (do_forces .OR. do_stress)
THEN
740 fr(k) = fr(k) + dp_i(1)*(tij_abc(1, 1, k)*dp_j(1) + &
741 tij_abc(2, 1, k)*dp_j(2) + &
742 tij_abc(3, 1, k)*dp_j(3)) &
743 + dp_i(2)*(tij_abc(1, 2, k)*dp_j(1) + &
744 tij_abc(2, 2, k)*dp_j(2) + &
745 tij_abc(3, 2, k)*dp_j(3)) &
746 + dp_i(3)*(tij_abc(1, 3, k)*dp_j(1) + &
747 tij_abc(2, 3, k)*dp_j(2) + &
748 tij_abc(3, 3, k)*dp_j(3))
755 ef0_i = ef0_i - (tij_a(1)*dp_j(1) + &
759 ef0_j = ef0_j + (tij_a(1)*dp_i(1) + &
765 ef1_i(1) = ef1_i(1) + (tij_ab(1, 1)*dp_j(1) + &
766 tij_ab(2, 1)*dp_j(2) + &
767 tij_ab(3, 1)*dp_j(3))
768 ef1_i(2) = ef1_i(2) + (tij_ab(1, 2)*dp_j(1) + &
769 tij_ab(2, 2)*dp_j(2) + &
770 tij_ab(3, 2)*dp_j(3))
771 ef1_i(3) = ef1_i(3) + (tij_ab(1, 3)*dp_j(1) + &
772 tij_ab(2, 3)*dp_j(2) + &
773 tij_ab(3, 3)*dp_j(3))
775 ef1_j(1) = ef1_j(1) + (tij_ab(1, 1)*dp_i(1) + &
776 tij_ab(2, 1)*dp_i(2) + &
777 tij_ab(3, 1)*dp_i(3))
778 ef1_j(2) = ef1_j(2) + (tij_ab(1, 2)*dp_i(1) + &
779 tij_ab(2, 2)*dp_i(2) + &
780 tij_ab(3, 2)*dp_i(3))
781 ef1_j(3) = ef1_j(3) + (tij_ab(1, 3)*dp_i(1) + &
782 tij_ab(2, 3)*dp_i(2) + &
783 tij_ab(3, 3)*dp_i(3))
787 ef2_i(1, 1) = ef2_i(1, 1) + (tij_abc(1, 1, 1)*dp_j(1) + &
788 tij_abc(2, 1, 1)*dp_j(2) + &
789 tij_abc(3, 1, 1)*dp_j(3))
790 ef2_i(1, 2) = ef2_i(1, 2) + (tij_abc(1, 1, 2)*dp_j(1) + &
791 tij_abc(2, 1, 2)*dp_j(2) + &
792 tij_abc(3, 1, 2)*dp_j(3))
793 ef2_i(1, 3) = ef2_i(1, 3) + (tij_abc(1, 1, 3)*dp_j(1) + &
794 tij_abc(2, 1, 3)*dp_j(2) + &
795 tij_abc(3, 1, 3)*dp_j(3))
796 ef2_i(2, 1) = ef2_i(2, 1) + (tij_abc(1, 2, 1)*dp_j(1) + &
797 tij_abc(2, 2, 1)*dp_j(2) + &
798 tij_abc(3, 2, 1)*dp_j(3))
799 ef2_i(2, 2) = ef2_i(2, 2) + (tij_abc(1, 2, 2)*dp_j(1) + &
800 tij_abc(2, 2, 2)*dp_j(2) + &
801 tij_abc(3, 2, 2)*dp_j(3))
802 ef2_i(2, 3) = ef2_i(2, 3) + (tij_abc(1, 2, 3)*dp_j(1) + &
803 tij_abc(2, 2, 3)*dp_j(2) + &
804 tij_abc(3, 2, 3)*dp_j(3))
805 ef2_i(3, 1) = ef2_i(3, 1) + (tij_abc(1, 3, 1)*dp_j(1) + &
806 tij_abc(2, 3, 1)*dp_j(2) + &
807 tij_abc(3, 3, 1)*dp_j(3))
808 ef2_i(3, 2) = ef2_i(3, 2) + (tij_abc(1, 3, 2)*dp_j(1) + &
809 tij_abc(2, 3, 2)*dp_j(2) + &
810 tij_abc(3, 3, 2)*dp_j(3))
811 ef2_i(3, 3) = ef2_i(3, 3) + (tij_abc(1, 3, 3)*dp_j(1) + &
812 tij_abc(2, 3, 3)*dp_j(2) + &
813 tij_abc(3, 3, 3)*dp_j(3))
815 ef2_j(1, 1) = ef2_j(1, 1) - (tij_abc(1, 1, 1)*dp_i(1) + &
816 tij_abc(2, 1, 1)*dp_i(2) + &
817 tij_abc(3, 1, 1)*dp_i(3))
818 ef2_j(1, 2) = ef2_j(1, 2) - (tij_abc(1, 1, 2)*dp_i(1) + &
819 tij_abc(2, 1, 2)*dp_i(2) + &
820 tij_abc(3, 1, 2)*dp_i(3))
821 ef2_j(1, 3) = ef2_j(1, 3) - (tij_abc(1, 1, 3)*dp_i(1) + &
822 tij_abc(2, 1, 3)*dp_i(2) + &
823 tij_abc(3, 1, 3)*dp_i(3))
824 ef2_j(2, 1) = ef2_j(2, 1) - (tij_abc(1, 2, 1)*dp_i(1) + &
825 tij_abc(2, 2, 1)*dp_i(2) + &
826 tij_abc(3, 2, 1)*dp_i(3))
827 ef2_j(2, 2) = ef2_j(2, 2) - (tij_abc(1, 2, 2)*dp_i(1) + &
828 tij_abc(2, 2, 2)*dp_i(2) + &
829 tij_abc(3, 2, 2)*dp_i(3))
830 ef2_j(2, 3) = ef2_j(2, 3) - (tij_abc(1, 2, 3)*dp_i(1) + &
831 tij_abc(2, 2, 3)*dp_i(2) + &
832 tij_abc(3, 2, 3)*dp_i(3))
833 ef2_j(3, 1) = ef2_j(3, 1) - (tij_abc(1, 3, 1)*dp_i(1) + &
834 tij_abc(2, 3, 1)*dp_i(2) + &
835 tij_abc(3, 3, 1)*dp_i(3))
836 ef2_j(3, 2) = ef2_j(3, 2) - (tij_abc(1, 3, 2)*dp_i(1) + &
837 tij_abc(2, 3, 2)*dp_i(2) + &
838 tij_abc(3, 3, 2)*dp_i(3))
839 ef2_j(3, 3) = ef2_j(3, 3) - (tij_abc(1, 3, 3)*dp_i(1) + &
840 tij_abc(2, 3, 3)*dp_i(2) + &
841 tij_abc(3, 3, 3)*dp_i(3))
847 tmp = ch_j*(tij_a(1)*dp_i(1) + &
850 - ch_i*(tij_a(1)*dp_j(1) + &
855 IF (do_forces .OR. do_stress)
THEN
857 fr(k) = fr(k) - ch_j*(tij_ab(1, k)*dp_i(1) + &
858 tij_ab(2, k)*dp_i(2) + &
859 tij_ab(3, k)*dp_i(3)) &
860 + ch_i*(tij_ab(1, k)*dp_j(1) + &
861 tij_ab(2, k)*dp_j(2) + &
862 tij_ab(3, k)*dp_j(3))
869 tmp11 = qp_i(1, 1)*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
870 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
871 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
872 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
873 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
874 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
875 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
876 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
877 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
878 tmp21 = qp_i(2, 1)*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
879 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
880 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
881 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
882 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
883 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
884 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
885 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
886 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
887 tmp31 = qp_i(3, 1)*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
888 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
889 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
890 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
891 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
892 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
893 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
894 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
895 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
896 tmp22 = qp_i(2, 2)*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
897 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
898 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
899 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
900 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
901 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
902 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
903 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
904 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
905 tmp32 = qp_i(3, 2)*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
906 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
907 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
908 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
909 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
910 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
911 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
912 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
913 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
914 tmp33 = qp_i(3, 3)*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
915 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
916 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
917 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
918 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
919 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
920 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
921 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
922 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
926 tmp = tmp11 + tmp12 + tmp13 + &
927 tmp21 + tmp22 + tmp23 + &
928 tmp31 + tmp32 + tmp33
930 eloc = eloc +
fac*tmp
932 IF (do_forces .OR. do_stress)
THEN
934 tmp11 = qp_i(1, 1)*(tij_abcde(1, 1, 1, 1, k)*qp_j(1, 1) + &
935 tij_abcde(2, 1, 1, 1, k)*qp_j(2, 1) + &
936 tij_abcde(3, 1, 1, 1, k)*qp_j(3, 1) + &
937 tij_abcde(1, 2, 1, 1, k)*qp_j(1, 2) + &
938 tij_abcde(2, 2, 1, 1, k)*qp_j(2, 2) + &
939 tij_abcde(3, 2, 1, 1, k)*qp_j(3, 2) + &
940 tij_abcde(1, 3, 1, 1, k)*qp_j(1, 3) + &
941 tij_abcde(2, 3, 1, 1, k)*qp_j(2, 3) + &
942 tij_abcde(3, 3, 1, 1, k)*qp_j(3, 3))
943 tmp21 = qp_i(2, 1)*(tij_abcde(1, 1, 2, 1, k)*qp_j(1, 1) + &
944 tij_abcde(2, 1, 2, 1, k)*qp_j(2, 1) + &
945 tij_abcde(3, 1, 2, 1, k)*qp_j(3, 1) + &
946 tij_abcde(1, 2, 2, 1, k)*qp_j(1, 2) + &
947 tij_abcde(2, 2, 2, 1, k)*qp_j(2, 2) + &
948 tij_abcde(3, 2, 2, 1, k)*qp_j(3, 2) + &
949 tij_abcde(1, 3, 2, 1, k)*qp_j(1, 3) + &
950 tij_abcde(2, 3, 2, 1, k)*qp_j(2, 3) + &
951 tij_abcde(3, 3, 2, 1, k)*qp_j(3, 3))
952 tmp31 = qp_i(3, 1)*(tij_abcde(1, 1, 3, 1, k)*qp_j(1, 1) + &
953 tij_abcde(2, 1, 3, 1, k)*qp_j(2, 1) + &
954 tij_abcde(3, 1, 3, 1, k)*qp_j(3, 1) + &
955 tij_abcde(1, 2, 3, 1, k)*qp_j(1, 2) + &
956 tij_abcde(2, 2, 3, 1, k)*qp_j(2, 2) + &
957 tij_abcde(3, 2, 3, 1, k)*qp_j(3, 2) + &
958 tij_abcde(1, 3, 3, 1, k)*qp_j(1, 3) + &
959 tij_abcde(2, 3, 3, 1, k)*qp_j(2, 3) + &
960 tij_abcde(3, 3, 3, 1, k)*qp_j(3, 3))
961 tmp22 = qp_i(2, 2)*(tij_abcde(1, 1, 2, 2, k)*qp_j(1, 1) + &
962 tij_abcde(2, 1, 2, 2, k)*qp_j(2, 1) + &
963 tij_abcde(3, 1, 2, 2, k)*qp_j(3, 1) + &
964 tij_abcde(1, 2, 2, 2, k)*qp_j(1, 2) + &
965 tij_abcde(2, 2, 2, 2, k)*qp_j(2, 2) + &
966 tij_abcde(3, 2, 2, 2, k)*qp_j(3, 2) + &
967 tij_abcde(1, 3, 2, 2, k)*qp_j(1, 3) + &
968 tij_abcde(2, 3, 2, 2, k)*qp_j(2, 3) + &
969 tij_abcde(3, 3, 2, 2, k)*qp_j(3, 3))
970 tmp32 = qp_i(3, 2)*(tij_abcde(1, 1, 3, 2, k)*qp_j(1, 1) + &
971 tij_abcde(2, 1, 3, 2, k)*qp_j(2, 1) + &
972 tij_abcde(3, 1, 3, 2, k)*qp_j(3, 1) + &
973 tij_abcde(1, 2, 3, 2, k)*qp_j(1, 2) + &
974 tij_abcde(2, 2, 3, 2, k)*qp_j(2, 2) + &
975 tij_abcde(3, 2, 3, 2, k)*qp_j(3, 2) + &
976 tij_abcde(1, 3, 3, 2, k)*qp_j(1, 3) + &
977 tij_abcde(2, 3, 3, 2, k)*qp_j(2, 3) + &
978 tij_abcde(3, 3, 3, 2, k)*qp_j(3, 3))
979 tmp33 = qp_i(3, 3)*(tij_abcde(1, 1, 3, 3, k)*qp_j(1, 1) + &
980 tij_abcde(2, 1, 3, 3, k)*qp_j(2, 1) + &
981 tij_abcde(3, 1, 3, 3, k)*qp_j(3, 1) + &
982 tij_abcde(1, 2, 3, 3, k)*qp_j(1, 2) + &
983 tij_abcde(2, 2, 3, 3, k)*qp_j(2, 2) + &
984 tij_abcde(3, 2, 3, 3, k)*qp_j(3, 2) + &
985 tij_abcde(1, 3, 3, 3, k)*qp_j(1, 3) + &
986 tij_abcde(2, 3, 3, 3, k)*qp_j(2, 3) + &
987 tij_abcde(3, 3, 3, 3, k)*qp_j(3, 3))
991 fr(k) = fr(k) -
fac*(tmp11 + tmp12 + tmp13 + &
992 tmp21 + tmp22 + tmp23 + &
993 tmp31 + tmp32 + tmp33)
1000 IF (do_efield0)
THEN
1001 ef0_i = ef0_i +
fac*(tij_ab(1, 1)*qp_j(1, 1) + &
1002 tij_ab(2, 1)*qp_j(2, 1) + &
1003 tij_ab(3, 1)*qp_j(3, 1) + &
1004 tij_ab(1, 2)*qp_j(1, 2) + &
1005 tij_ab(2, 2)*qp_j(2, 2) + &
1006 tij_ab(3, 2)*qp_j(3, 2) + &
1007 tij_ab(1, 3)*qp_j(1, 3) + &
1008 tij_ab(2, 3)*qp_j(2, 3) + &
1009 tij_ab(3, 3)*qp_j(3, 3))
1011 ef0_j = ef0_j +
fac*(tij_ab(1, 1)*qp_i(1, 1) + &
1012 tij_ab(2, 1)*qp_i(2, 1) + &
1013 tij_ab(3, 1)*qp_i(3, 1) + &
1014 tij_ab(1, 2)*qp_i(1, 2) + &
1015 tij_ab(2, 2)*qp_i(2, 2) + &
1016 tij_ab(3, 2)*qp_i(3, 2) + &
1017 tij_ab(1, 3)*qp_i(1, 3) + &
1018 tij_ab(2, 3)*qp_i(2, 3) + &
1019 tij_ab(3, 3)*qp_i(3, 3))
1022 IF (do_efield1)
THEN
1023 ef1_i(1) = ef1_i(1) -
fac*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
1024 tij_abc(2, 1, 1)*qp_j(2, 1) + &
1025 tij_abc(3, 1, 1)*qp_j(3, 1) + &
1026 tij_abc(1, 2, 1)*qp_j(1, 2) + &
1027 tij_abc(2, 2, 1)*qp_j(2, 2) + &
1028 tij_abc(3, 2, 1)*qp_j(3, 2) + &
1029 tij_abc(1, 3, 1)*qp_j(1, 3) + &
1030 tij_abc(2, 3, 1)*qp_j(2, 3) + &
1031 tij_abc(3, 3, 1)*qp_j(3, 3))
1032 ef1_i(2) = ef1_i(2) -
fac*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
1033 tij_abc(2, 1, 2)*qp_j(2, 1) + &
1034 tij_abc(3, 1, 2)*qp_j(3, 1) + &
1035 tij_abc(1, 2, 2)*qp_j(1, 2) + &
1036 tij_abc(2, 2, 2)*qp_j(2, 2) + &
1037 tij_abc(3, 2, 2)*qp_j(3, 2) + &
1038 tij_abc(1, 3, 2)*qp_j(1, 3) + &
1039 tij_abc(2, 3, 2)*qp_j(2, 3) + &
1040 tij_abc(3, 3, 2)*qp_j(3, 3))
1041 ef1_i(3) = ef1_i(3) -
fac*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
1042 tij_abc(2, 1, 3)*qp_j(2, 1) + &
1043 tij_abc(3, 1, 3)*qp_j(3, 1) + &
1044 tij_abc(1, 2, 3)*qp_j(1, 2) + &
1045 tij_abc(2, 2, 3)*qp_j(2, 2) + &
1046 tij_abc(3, 2, 3)*qp_j(3, 2) + &
1047 tij_abc(1, 3, 3)*qp_j(1, 3) + &
1048 tij_abc(2, 3, 3)*qp_j(2, 3) + &
1049 tij_abc(3, 3, 3)*qp_j(3, 3))
1051 ef1_j(1) = ef1_j(1) +
fac*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
1052 tij_abc(2, 1, 1)*qp_i(2, 1) + &
1053 tij_abc(3, 1, 1)*qp_i(3, 1) + &
1054 tij_abc(1, 2, 1)*qp_i(1, 2) + &
1055 tij_abc(2, 2, 1)*qp_i(2, 2) + &
1056 tij_abc(3, 2, 1)*qp_i(3, 2) + &
1057 tij_abc(1, 3, 1)*qp_i(1, 3) + &
1058 tij_abc(2, 3, 1)*qp_i(2, 3) + &
1059 tij_abc(3, 3, 1)*qp_i(3, 3))
1060 ef1_j(2) = ef1_j(2) +
fac*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
1061 tij_abc(2, 1, 2)*qp_i(2, 1) + &
1062 tij_abc(3, 1, 2)*qp_i(3, 1) + &
1063 tij_abc(1, 2, 2)*qp_i(1, 2) + &
1064 tij_abc(2, 2, 2)*qp_i(2, 2) + &
1065 tij_abc(3, 2, 2)*qp_i(3, 2) + &
1066 tij_abc(1, 3, 2)*qp_i(1, 3) + &
1067 tij_abc(2, 3, 2)*qp_i(2, 3) + &
1068 tij_abc(3, 3, 2)*qp_i(3, 3))
1069 ef1_j(3) = ef1_j(3) +
fac*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
1070 tij_abc(2, 1, 3)*qp_i(2, 1) + &
1071 tij_abc(3, 1, 3)*qp_i(3, 1) + &
1072 tij_abc(1, 2, 3)*qp_i(1, 2) + &
1073 tij_abc(2, 2, 3)*qp_i(2, 2) + &
1074 tij_abc(3, 2, 3)*qp_i(3, 2) + &
1075 tij_abc(1, 3, 3)*qp_i(1, 3) + &
1076 tij_abc(2, 3, 3)*qp_i(2, 3) + &
1077 tij_abc(3, 3, 3)*qp_i(3, 3))
1080 IF (do_efield2)
THEN
1081 tmp11 =
fac*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
1082 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
1083 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
1084 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
1085 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
1086 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
1087 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
1088 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
1089 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
1090 tmp12 =
fac*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
1091 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
1092 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
1093 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
1094 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
1095 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
1096 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
1097 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
1098 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
1099 tmp13 =
fac*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
1100 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
1101 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
1102 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
1103 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
1104 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
1105 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
1106 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
1107 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
1108 tmp22 =
fac*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
1109 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
1110 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
1111 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
1112 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
1113 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
1114 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
1115 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
1116 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
1117 tmp23 =
fac*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
1118 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
1119 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
1120 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
1121 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
1122 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
1123 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
1124 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
1125 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
1126 tmp33 =
fac*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
1127 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
1128 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
1129 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
1130 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
1131 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
1132 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
1133 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
1134 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
1136 ef2_i(1, 1) = ef2_i(1, 1) - tmp11
1137 ef2_i(1, 2) = ef2_i(1, 2) - tmp12
1138 ef2_i(1, 3) = ef2_i(1, 3) - tmp13
1139 ef2_i(2, 1) = ef2_i(2, 1) - tmp12
1140 ef2_i(2, 2) = ef2_i(2, 2) - tmp22
1141 ef2_i(2, 3) = ef2_i(2, 3) - tmp23
1142 ef2_i(3, 1) = ef2_i(3, 1) - tmp13
1143 ef2_i(3, 2) = ef2_i(3, 2) - tmp23
1144 ef2_i(3, 3) = ef2_i(3, 3) - tmp33
1146 tmp11 =
fac*(tij_abcd(1, 1, 1, 1)*qp_i(1, 1) + &
1147 tij_abcd(2, 1, 1, 1)*qp_i(2, 1) + &
1148 tij_abcd(3, 1, 1, 1)*qp_i(3, 1) + &
1149 tij_abcd(1, 2, 1, 1)*qp_i(1, 2) + &
1150 tij_abcd(2, 2, 1, 1)*qp_i(2, 2) + &
1151 tij_abcd(3, 2, 1, 1)*qp_i(3, 2) + &
1152 tij_abcd(1, 3, 1, 1)*qp_i(1, 3) + &
1153 tij_abcd(2, 3, 1, 1)*qp_i(2, 3) + &
1154 tij_abcd(3, 3, 1, 1)*qp_i(3, 3))
1155 tmp12 =
fac*(tij_abcd(1, 1, 1, 2)*qp_i(1, 1) + &
1156 tij_abcd(2, 1, 1, 2)*qp_i(2, 1) + &
1157 tij_abcd(3, 1, 1, 2)*qp_i(3, 1) + &
1158 tij_abcd(1, 2, 1, 2)*qp_i(1, 2) + &
1159 tij_abcd(2, 2, 1, 2)*qp_i(2, 2) + &
1160 tij_abcd(3, 2, 1, 2)*qp_i(3, 2) + &
1161 tij_abcd(1, 3, 1, 2)*qp_i(1, 3) + &
1162 tij_abcd(2, 3, 1, 2)*qp_i(2, 3) + &
1163 tij_abcd(3, 3, 1, 2)*qp_i(3, 3))
1164 tmp13 =
fac*(tij_abcd(1, 1, 1, 3)*qp_i(1, 1) + &
1165 tij_abcd(2, 1, 1, 3)*qp_i(2, 1) + &
1166 tij_abcd(3, 1, 1, 3)*qp_i(3, 1) + &
1167 tij_abcd(1, 2, 1, 3)*qp_i(1, 2) + &
1168 tij_abcd(2, 2, 1, 3)*qp_i(2, 2) + &
1169 tij_abcd(3, 2, 1, 3)*qp_i(3, 2) + &
1170 tij_abcd(1, 3, 1, 3)*qp_i(1, 3) + &
1171 tij_abcd(2, 3, 1, 3)*qp_i(2, 3) + &
1172 tij_abcd(3, 3, 1, 3)*qp_i(3, 3))
1173 tmp22 =
fac*(tij_abcd(1, 1, 2, 2)*qp_i(1, 1) + &
1174 tij_abcd(2, 1, 2, 2)*qp_i(2, 1) + &
1175 tij_abcd(3, 1, 2, 2)*qp_i(3, 1) + &
1176 tij_abcd(1, 2, 2, 2)*qp_i(1, 2) + &
1177 tij_abcd(2, 2, 2, 2)*qp_i(2, 2) + &
1178 tij_abcd(3, 2, 2, 2)*qp_i(3, 2) + &
1179 tij_abcd(1, 3, 2, 2)*qp_i(1, 3) + &
1180 tij_abcd(2, 3, 2, 2)*qp_i(2, 3) + &
1181 tij_abcd(3, 3, 2, 2)*qp_i(3, 3))
1182 tmp23 =
fac*(tij_abcd(1, 1, 2, 3)*qp_i(1, 1) + &
1183 tij_abcd(2, 1, 2, 3)*qp_i(2, 1) + &
1184 tij_abcd(3, 1, 2, 3)*qp_i(3, 1) + &
1185 tij_abcd(1, 2, 2, 3)*qp_i(1, 2) + &
1186 tij_abcd(2, 2, 2, 3)*qp_i(2, 2) + &
1187 tij_abcd(3, 2, 2, 3)*qp_i(3, 2) + &
1188 tij_abcd(1, 3, 2, 3)*qp_i(1, 3) + &
1189 tij_abcd(2, 3, 2, 3)*qp_i(2, 3) + &
1190 tij_abcd(3, 3, 2, 3)*qp_i(3, 3))
1191 tmp33 =
fac*(tij_abcd(1, 1, 3, 3)*qp_i(1, 1) + &
1192 tij_abcd(2, 1, 3, 3)*qp_i(2, 1) + &
1193 tij_abcd(3, 1, 3, 3)*qp_i(3, 1) + &
1194 tij_abcd(1, 2, 3, 3)*qp_i(1, 2) + &
1195 tij_abcd(2, 2, 3, 3)*qp_i(2, 2) + &
1196 tij_abcd(3, 2, 3, 3)*qp_i(3, 2) + &
1197 tij_abcd(1, 3, 3, 3)*qp_i(1, 3) + &
1198 tij_abcd(2, 3, 3, 3)*qp_i(2, 3) + &
1199 tij_abcd(3, 3, 3, 3)*qp_i(3, 3))
1201 ef2_j(1, 1) = ef2_j(1, 1) - tmp11
1202 ef2_j(1, 2) = ef2_j(1, 2) - tmp12
1203 ef2_j(1, 3) = ef2_j(1, 3) - tmp13
1204 ef2_j(2, 1) = ef2_j(2, 1) - tmp12
1205 ef2_j(2, 2) = ef2_j(2, 2) - tmp22
1206 ef2_j(2, 3) = ef2_j(2, 3) - tmp23
1207 ef2_j(3, 1) = ef2_j(3, 1) - tmp13
1208 ef2_j(3, 2) = ef2_j(3, 2) - tmp23
1209 ef2_j(3, 3) = ef2_j(3, 3) - tmp33
1213 IF (task(3, 2))
THEN
1217 tmp_ij = dp_i(1)*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
1218 tij_abc(2, 1, 1)*qp_j(2, 1) + &
1219 tij_abc(3, 1, 1)*qp_j(3, 1) + &
1220 tij_abc(1, 2, 1)*qp_j(1, 2) + &
1221 tij_abc(2, 2, 1)*qp_j(2, 2) + &
1222 tij_abc(3, 2, 1)*qp_j(3, 2) + &
1223 tij_abc(1, 3, 1)*qp_j(1, 3) + &
1224 tij_abc(2, 3, 1)*qp_j(2, 3) + &
1225 tij_abc(3, 3, 1)*qp_j(3, 3)) + &
1226 dp_i(2)*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
1227 tij_abc(2, 1, 2)*qp_j(2, 1) + &
1228 tij_abc(3, 1, 2)*qp_j(3, 1) + &
1229 tij_abc(1, 2, 2)*qp_j(1, 2) + &
1230 tij_abc(2, 2, 2)*qp_j(2, 2) + &
1231 tij_abc(3, 2, 2)*qp_j(3, 2) + &
1232 tij_abc(1, 3, 2)*qp_j(1, 3) + &
1233 tij_abc(2, 3, 2)*qp_j(2, 3) + &
1234 tij_abc(3, 3, 2)*qp_j(3, 3)) + &
1235 dp_i(3)*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
1236 tij_abc(2, 1, 3)*qp_j(2, 1) + &
1237 tij_abc(3, 1, 3)*qp_j(3, 1) + &
1238 tij_abc(1, 2, 3)*qp_j(1, 2) + &
1239 tij_abc(2, 2, 3)*qp_j(2, 2) + &
1240 tij_abc(3, 2, 3)*qp_j(3, 2) + &
1241 tij_abc(1, 3, 3)*qp_j(1, 3) + &
1242 tij_abc(2, 3, 3)*qp_j(2, 3) + &
1243 tij_abc(3, 3, 3)*qp_j(3, 3))
1246 tmp_ji = dp_j(1)*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
1247 tij_abc(2, 1, 1)*qp_i(2, 1) + &
1248 tij_abc(3, 1, 1)*qp_i(3, 1) + &
1249 tij_abc(1, 2, 1)*qp_i(1, 2) + &
1250 tij_abc(2, 2, 1)*qp_i(2, 2) + &
1251 tij_abc(3, 2, 1)*qp_i(3, 2) + &
1252 tij_abc(1, 3, 1)*qp_i(1, 3) + &
1253 tij_abc(2, 3, 1)*qp_i(2, 3) + &
1254 tij_abc(3, 3, 1)*qp_i(3, 3)) + &
1255 dp_j(2)*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
1256 tij_abc(2, 1, 2)*qp_i(2, 1) + &
1257 tij_abc(3, 1, 2)*qp_i(3, 1) + &
1258 tij_abc(1, 2, 2)*qp_i(1, 2) + &
1259 tij_abc(2, 2, 2)*qp_i(2, 2) + &
1260 tij_abc(3, 2, 2)*qp_i(3, 2) + &
1261 tij_abc(1, 3, 2)*qp_i(1, 3) + &
1262 tij_abc(2, 3, 2)*qp_i(2, 3) + &
1263 tij_abc(3, 3, 2)*qp_i(3, 3)) + &
1264 dp_j(3)*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
1265 tij_abc(2, 1, 3)*qp_i(2, 1) + &
1266 tij_abc(3, 1, 3)*qp_i(3, 1) + &
1267 tij_abc(1, 2, 3)*qp_i(1, 2) + &
1268 tij_abc(2, 2, 3)*qp_i(2, 2) + &
1269 tij_abc(3, 2, 3)*qp_i(3, 2) + &
1270 tij_abc(1, 3, 3)*qp_i(1, 3) + &
1271 tij_abc(2, 3, 3)*qp_i(2, 3) + &
1272 tij_abc(3, 3, 3)*qp_i(3, 3))
1274 tmp =
fac*(tmp_ij - tmp_ji)
1276 IF (do_forces .OR. do_stress)
THEN
1279 tmp_ij = dp_i(1)*(tij_abcd(1, 1, 1, k)*qp_j(1, 1) + &
1280 tij_abcd(2, 1, 1, k)*qp_j(2, 1) + &
1281 tij_abcd(3, 1, 1, k)*qp_j(3, 1) + &
1282 tij_abcd(1, 2, 1, k)*qp_j(1, 2) + &
1283 tij_abcd(2, 2, 1, k)*qp_j(2, 2) + &
1284 tij_abcd(3, 2, 1, k)*qp_j(3, 2) + &
1285 tij_abcd(1, 3, 1, k)*qp_j(1, 3) + &
1286 tij_abcd(2, 3, 1, k)*qp_j(2, 3) + &
1287 tij_abcd(3, 3, 1, k)*qp_j(3, 3)) + &
1288 dp_i(2)*(tij_abcd(1, 1, 2, k)*qp_j(1, 1) + &
1289 tij_abcd(2, 1, 2, k)*qp_j(2, 1) + &
1290 tij_abcd(3, 1, 2, k)*qp_j(3, 1) + &
1291 tij_abcd(1, 2, 2, k)*qp_j(1, 2) + &
1292 tij_abcd(2, 2, 2, k)*qp_j(2, 2) + &
1293 tij_abcd(3, 2, 2, k)*qp_j(3, 2) + &
1294 tij_abcd(1, 3, 2, k)*qp_j(1, 3) + &
1295 tij_abcd(2, 3, 2, k)*qp_j(2, 3) + &
1296 tij_abcd(3, 3, 2, k)*qp_j(3, 3)) + &
1297 dp_i(3)*(tij_abcd(1, 1, 3, k)*qp_j(1, 1) + &
1298 tij_abcd(2, 1, 3, k)*qp_j(2, 1) + &
1299 tij_abcd(3, 1, 3, k)*qp_j(3, 1) + &
1300 tij_abcd(1, 2, 3, k)*qp_j(1, 2) + &
1301 tij_abcd(2, 2, 3, k)*qp_j(2, 2) + &
1302 tij_abcd(3, 2, 3, k)*qp_j(3, 2) + &
1303 tij_abcd(1, 3, 3, k)*qp_j(1, 3) + &
1304 tij_abcd(2, 3, 3, k)*qp_j(2, 3) + &
1305 tij_abcd(3, 3, 3, k)*qp_j(3, 3))
1308 tmp_ji = dp_j(1)*(tij_abcd(1, 1, 1, k)*qp_i(1, 1) + &
1309 tij_abcd(2, 1, 1, k)*qp_i(2, 1) + &
1310 tij_abcd(3, 1, 1, k)*qp_i(3, 1) + &
1311 tij_abcd(1, 2, 1, k)*qp_i(1, 2) + &
1312 tij_abcd(2, 2, 1, k)*qp_i(2, 2) + &
1313 tij_abcd(3, 2, 1, k)*qp_i(3, 2) + &
1314 tij_abcd(1, 3, 1, k)*qp_i(1, 3) + &
1315 tij_abcd(2, 3, 1, k)*qp_i(2, 3) + &
1316 tij_abcd(3, 3, 1, k)*qp_i(3, 3)) + &
1317 dp_j(2)*(tij_abcd(1, 1, 2, k)*qp_i(1, 1) + &
1318 tij_abcd(2, 1, 2, k)*qp_i(2, 1) + &
1319 tij_abcd(3, 1, 2, k)*qp_i(3, 1) + &
1320 tij_abcd(1, 2, 2, k)*qp_i(1, 2) + &
1321 tij_abcd(2, 2, 2, k)*qp_i(2, 2) + &
1322 tij_abcd(3, 2, 2, k)*qp_i(3, 2) + &
1323 tij_abcd(1, 3, 2, k)*qp_i(1, 3) + &
1324 tij_abcd(2, 3, 2, k)*qp_i(2, 3) + &
1325 tij_abcd(3, 3, 2, k)*qp_i(3, 3)) + &
1326 dp_j(3)*(tij_abcd(1, 1, 3, k)*qp_i(1, 1) + &
1327 tij_abcd(2, 1, 3, k)*qp_i(2, 1) + &
1328 tij_abcd(3, 1, 3, k)*qp_i(3, 1) + &
1329 tij_abcd(1, 2, 3, k)*qp_i(1, 2) + &
1330 tij_abcd(2, 2, 3, k)*qp_i(2, 2) + &
1331 tij_abcd(3, 2, 3, k)*qp_i(3, 2) + &
1332 tij_abcd(1, 3, 3, k)*qp_i(1, 3) + &
1333 tij_abcd(2, 3, 3, k)*qp_i(2, 3) + &
1334 tij_abcd(3, 3, 3, k)*qp_i(3, 3))
1336 fr(k) = fr(k) -
fac*(tmp_ij - tmp_ji)
1340 IF (task(3, 1))
THEN
1345 tmp_ij = ch_i*(tij_ab(1, 1)*qp_j(1, 1) + &
1346 tij_ab(2, 1)*qp_j(2, 1) + &
1347 tij_ab(3, 1)*qp_j(3, 1) + &
1348 tij_ab(1, 2)*qp_j(1, 2) + &
1349 tij_ab(2, 2)*qp_j(2, 2) + &
1350 tij_ab(3, 2)*qp_j(3, 2) + &
1351 tij_ab(1, 3)*qp_j(1, 3) + &
1352 tij_ab(2, 3)*qp_j(2, 3) + &
1353 tij_ab(3, 3)*qp_j(3, 3))
1356 tmp_ji = ch_j*(tij_ab(1, 1)*qp_i(1, 1) + &
1357 tij_ab(2, 1)*qp_i(2, 1) + &
1358 tij_ab(3, 1)*qp_i(3, 1) + &
1359 tij_ab(1, 2)*qp_i(1, 2) + &
1360 tij_ab(2, 2)*qp_i(2, 2) + &
1361 tij_ab(3, 2)*qp_i(3, 2) + &
1362 tij_ab(1, 3)*qp_i(1, 3) + &
1363 tij_ab(2, 3)*qp_i(2, 3) + &
1364 tij_ab(3, 3)*qp_i(3, 3))
1366 eloc = eloc +
fac*(tmp_ij + tmp_ji)
1367 IF (do_forces .OR. do_stress)
THEN
1370 tmp_ij = ch_i*(tij_abc(1, 1, k)*qp_j(1, 1) + &
1371 tij_abc(2, 1, k)*qp_j(2, 1) + &
1372 tij_abc(3, 1, k)*qp_j(3, 1) + &
1373 tij_abc(1, 2, k)*qp_j(1, 2) + &
1374 tij_abc(2, 2, k)*qp_j(2, 2) + &
1375 tij_abc(3, 2, k)*qp_j(3, 2) + &
1376 tij_abc(1, 3, k)*qp_j(1, 3) + &
1377 tij_abc(2, 3, k)*qp_j(2, 3) + &
1378 tij_abc(3, 3, k)*qp_j(3, 3))
1381 tmp_ji = ch_j*(tij_abc(1, 1, k)*qp_i(1, 1) + &
1382 tij_abc(2, 1, k)*qp_i(2, 1) + &
1383 tij_abc(3, 1, k)*qp_i(3, 1) + &
1384 tij_abc(1, 2, k)*qp_i(1, 2) + &
1385 tij_abc(2, 2, k)*qp_i(2, 2) + &
1386 tij_abc(3, 2, k)*qp_i(3, 2) + &
1387 tij_abc(1, 3, k)*qp_i(1, 3) + &
1388 tij_abc(2, 3, k)*qp_i(2, 3) + &
1389 tij_abc(3, 3, k)*qp_i(3, 3))
1391 fr(k) = fr(k) -
fac*(tmp_ij + tmp_ji)
1395 energy = energy + eloc
1397 forces(1, atom_a) = forces(1, atom_a) - fr(1)
1398 forces(2, atom_a) = forces(2, atom_a) - fr(2)
1399 forces(3, atom_a) = forces(3, atom_a) - fr(3)
1400 forces(1, atom_b) = forces(1, atom_b) + fr(1)
1401 forces(2, atom_b) = forces(2, atom_b) + fr(2)
1402 forces(3, atom_b) = forces(3, atom_b) + fr(3)
1407 IF (do_efield0)
THEN
1408 efield0(atom_a) = efield0(atom_a) + ef0_j
1410 efield0(atom_b) = efield0(atom_b) + ef0_i
1413 IF (do_efield1)
THEN
1414 efield1(1, atom_a) = efield1(1, atom_a) + ef1_j(1)
1415 efield1(2, atom_a) = efield1(2, atom_a) + ef1_j(2)
1416 efield1(3, atom_a) = efield1(3, atom_a) + ef1_j(3)
1418 efield1(1, atom_b) = efield1(1, atom_b) + ef1_i(1)
1419 efield1(2, atom_b) = efield1(2, atom_b) + ef1_i(2)
1420 efield1(3, atom_b) = efield1(3, atom_b) + ef1_i(3)
1423 IF (do_efield2)
THEN
1424 efield2(1, atom_a) = efield2(1, atom_a) + ef2_j(1, 1)
1425 efield2(2, atom_a) = efield2(2, atom_a) + ef2_j(1, 2)
1426 efield2(3, atom_a) = efield2(3, atom_a) + ef2_j(1, 3)
1427 efield2(4, atom_a) = efield2(4, atom_a) + ef2_j(2, 1)
1428 efield2(5, atom_a) = efield2(5, atom_a) + ef2_j(2, 2)
1429 efield2(6, atom_a) = efield2(6, atom_a) + ef2_j(2, 3)
1430 efield2(7, atom_a) = efield2(7, atom_a) + ef2_j(3, 1)
1431 efield2(8, atom_a) = efield2(8, atom_a) + ef2_j(3, 2)
1432 efield2(9, atom_a) = efield2(9, atom_a) + ef2_j(3, 3)
1434 efield2(1, atom_b) = efield2(1, atom_b) + ef2_i(1, 1)
1435 efield2(2, atom_b) = efield2(2, atom_b) + ef2_i(1, 2)
1436 efield2(3, atom_b) = efield2(3, atom_b) + ef2_i(1, 3)
1437 efield2(4, atom_b) = efield2(4, atom_b) + ef2_i(2, 1)
1438 efield2(5, atom_b) = efield2(5, atom_b) + ef2_i(2, 2)
1439 efield2(6, atom_b) = efield2(6, atom_b) + ef2_i(2, 3)
1440 efield2(7, atom_b) = efield2(7, atom_b) + ef2_i(3, 1)
1441 efield2(8, atom_b) = efield2(8, atom_b) + ef2_i(3, 2)
1442 efield2(9, atom_b) = efield2(9, atom_b) + ef2_i(3, 3)
1446 ptens11 = ptens11 + rab(1)*fr(1)
1447 ptens21 = ptens21 + rab(2)*fr(1)
1448 ptens31 = ptens31 + rab(3)*fr(1)
1449 ptens12 = ptens12 + rab(1)*fr(2)
1450 ptens22 = ptens22 + rab(2)*fr(2)
1451 ptens32 = ptens32 + rab(3)*fr(2)
1452 ptens13 = ptens13 + rab(1)*fr(3)
1453 ptens23 = ptens23 + rab(2)*fr(3)
1454 ptens33 = ptens33 + rab(3)*fr(3)
1461 tij_a = huge(0.0_dp)
1462 tij_ab = huge(0.0_dp)
1463 tij_abc = huge(0.0_dp)
1464 tij_abcd = huge(0.0_dp)
1465 tij_abcde = huge(0.0_dp)
1472 IF (debug_this_module .AND. debug_r_space .AND. (.NOT. debug_g_space))
THEN
1476 f(0) = erfc(alpha*r)*ir
1482 f(i) = irab2*(f(i - 1) + tmp*((2.0_dp*alpha**2)**i)/(
fac*alpha))
1487 force_eval = do_stress
1488 IF (task(1, 1))
THEN
1490 force_eval = do_forces .OR. do_efield1
1492 IF (task(2, 2)) force_eval = force_eval .OR. do_efield0
1493 IF (task(1, 2) .OR. force_eval)
THEN
1494 force_eval = do_stress
1495 tij_a = -rab*f(1)*fac_ij
1496 IF (task(1, 2)) force_eval = force_eval .OR. do_forces
1498 IF (task(1, 1)) force_eval = force_eval .OR. do_efield2
1499 IF (task(3, 3)) force_eval = force_eval .OR. do_efield0
1500 IF (task(2, 2) .OR. task(3, 1) .OR. force_eval)
THEN
1501 force_eval = do_stress
1504 tmp = rab(a)*rab(b)*fac_ij
1505 tij_ab(a, b) = 3.0_dp*tmp*f(2)
1506 IF (a == b) tij_ab(a, b) = tij_ab(a, b) - f(1)*fac_ij
1509 IF (task(2, 2) .OR. task(3, 1)) force_eval = force_eval .OR. do_forces
1511 IF (task(2, 2)) force_eval = force_eval .OR. do_efield2
1512 IF (task(3, 3)) force_eval = force_eval .OR. do_efield1
1513 IF (task(3, 2) .OR. force_eval)
THEN
1514 force_eval = do_stress
1518 tmp = rab(a)*rab(b)*rab(c)*fac_ij
1519 tij_abc(a, b, c) = -15.0_dp*tmp*f(3)
1520 tmp = 3.0_dp*f(2)*fac_ij
1521 IF (a == b) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(c)
1522 IF (a == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(b)
1523 IF (b == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(a)
1527 IF (task(3, 2)) force_eval = force_eval .OR. do_forces
1529 IF (task(3, 3) .OR. force_eval)
THEN
1530 force_eval = do_stress
1535 tmp = rab(a)*rab(b)*rab(c)*rab(d)*fac_ij
1536 tij_abcd(a, b, c, d) = 105.0_dp*tmp*f(4)
1537 tmp1 = 15.0_dp*f(3)*fac_ij
1538 tmp2 = 3.0_dp*f(2)*fac_ij
1540 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(c)*rab(d)
1541 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
1544 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(d)
1545 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
1547 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(c)
1549 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(d)
1550 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
1552 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(c)
1553 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(b)
1558 IF (task(3, 3)) force_eval = force_eval .OR. do_forces
1560 IF (force_eval)
THEN
1561 force_eval = do_stress
1567 tmp = rab(a)*rab(b)*rab(c)*rab(d)*rab(e)*fac_ij
1568 tij_abcde(a, b, c, d, e) = -945.0_dp*tmp*f(5)
1569 tmp1 = 105.0_dp*f(4)*fac_ij
1570 tmp2 = 15.0_dp*f(3)*fac_ij
1572 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(c)*rab(d)*rab(e)
1573 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
1574 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
1575 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
1578 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(d)*rab(e)
1579 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
1580 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
1581 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
1584 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(e)
1585 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
1586 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
1587 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
1590 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(d)
1591 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
1592 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
1593 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
1596 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(d)*rab(e)
1597 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
1600 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(e)
1601 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
1604 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(d)
1605 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
1607 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(e)
1608 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(d)
1609 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(c)
1627 IF (kind_a == ikind)
THEN
1629 SELECT CASE (itype_ij)
1634 DO kk = 1, nkdamp_ij
1636 factorial = factorial*real(kk, kind=
dp)
1637 dampsumfi = dampsumfi + (xf/factorial)
1639 dampaexpi = dexp(-dampa_ij*r)
1640 dampfunci = dampsumfi*dampaexpi*dampfac_ij
1641 dampfuncdiffi = -dampa_ij*dampaexpi* &
1642 dampfac_ij*(((dampa_ij*r)**nkdamp_ij)/ &
1646 dampfuncdiffi = 0.0_dp
1650 SELECT CASE (itype_ji)
1655 DO kk = 1, nkdamp_ji
1657 factorial = factorial*real(kk, kind=
dp)
1658 dampsumfj = dampsumfj + (xf/factorial)
1660 dampaexpj = dexp(-dampa_ji*r)
1661 dampfuncj = dampsumfj*dampaexpj*dampfac_ji
1662 dampfuncdiffj = -dampa_ji*dampaexpj* &
1663 dampfac_ji*(((dampa_ji*r)**nkdamp_ji)/ &
1667 dampfuncdiffj = 0.0_dp
1670 SELECT CASE (itype_ij)
1675 DO kk = 1, nkdamp_ij
1677 factorial = factorial*real(kk, kind=
dp)
1678 dampsumfj = dampsumfj + (xf/factorial)
1680 dampaexpj = dexp(-dampa_ij*r)
1681 dampfuncj = dampsumfj*dampaexpj*dampfac_ij
1682 dampfuncdiffj = -dampa_ij*dampaexpj* &
1683 dampfac_ij*(((dampa_ij*r)**nkdamp_ij)/ &
1687 dampfuncdiffj = 0.0_dp
1691 SELECT CASE (itype_ji)
1696 DO kk = 1, nkdamp_ji
1698 factorial = factorial*real(kk, kind=
dp)
1699 dampsumfi = dampsumfi + (xf/factorial)
1701 dampaexpi = dexp(-dampa_ji*r)
1702 dampfunci = dampsumfi*dampaexpi*dampfac_ji
1703 dampfuncdiffi = -dampa_ji*dampaexpi* &
1704 dampfac_ji*(((dampa_ji*r)**nkdamp_ji)/ &
1708 dampfuncdiffi = 0.0_dp
1712 damptij_a = -rab*dampfunci*fac_ij*irab2*ir
1713 damptji_a = -rab*dampfuncj*fac_ij*irab2*ir
1716 tmp = rab(a)*rab(b)*fac_ij
1717 damptij_ab(a, b) = tmp*(-dampfuncdiffi*irab2*irab2 + 3.0_dp*dampfunci*irab2*irab2*ir)
1718 damptji_ab(a, b) = tmp*(-dampfuncdiffj*irab2*irab2 + 3.0_dp*dampfuncj*irab2*irab2*ir)
1719 IF (a == b) damptij_ab(a, b) = damptij_ab(a, b) - dampfunci*fac_ij*irab2*ir
1720 IF (a == b) damptji_ab(a, b) = damptji_ab(a, b) - dampfuncj*fac_ij*irab2*ir
1726 IF (debug_this_module)
THEN
1734 IF (any(task(1, :)))
THEN
1735 ch_j = charges(atom_a)
1736 ch_i = charges(atom_b)
1738 IF (any(task(2, :)))
THEN
1739 dp_j = dipoles(:, atom_a)
1740 dp_i = dipoles(:, atom_b)
1742 IF (any(task(3, :)))
THEN
1743 qp_j = quadrupoles(:, :, atom_a)
1744 qp_i = quadrupoles(:, :, atom_b)
1746 IF (task(1, 1))
THEN
1748 eloc = eloc + ch_i*tij*ch_j
1750 IF (do_forces .OR. do_stress)
THEN
1751 fr(1) = fr(1) - ch_j*tij_a(1)*ch_i
1752 fr(2) = fr(2) - ch_j*tij_a(2)*ch_i
1753 fr(3) = fr(3) - ch_j*tij_a(3)*ch_i
1758 IF (do_efield0)
THEN
1759 ef0_i = ef0_i + tij*ch_j
1761 ef0_j = ef0_j + tij*ch_i
1764 IF (do_efield1)
THEN
1765 ef1_i(1) = ef1_i(1) - tij_a(1)*ch_j
1766 ef1_i(2) = ef1_i(2) - tij_a(2)*ch_j
1767 ef1_i(3) = ef1_i(3) - tij_a(3)*ch_j
1769 ef1_j(1) = ef1_j(1) + tij_a(1)*ch_i
1770 ef1_j(2) = ef1_j(2) + tij_a(2)*ch_i
1771 ef1_j(3) = ef1_j(3) + tij_a(3)*ch_i
1773 ef1_i(1) = ef1_i(1) + damptij_a(1)*ch_j
1774 ef1_i(2) = ef1_i(2) + damptij_a(2)*ch_j
1775 ef1_i(3) = ef1_i(3) + damptij_a(3)*ch_j
1777 ef1_j(1) = ef1_j(1) - damptji_a(1)*ch_i
1778 ef1_j(2) = ef1_j(2) - damptji_a(2)*ch_i
1779 ef1_j(3) = ef1_j(3) - damptji_a(3)*ch_i
1783 IF (do_efield2)
THEN
1784 ef2_i(1, 1) = ef2_i(1, 1) - tij_ab(1, 1)*ch_j
1785 ef2_i(2, 1) = ef2_i(2, 1) - tij_ab(2, 1)*ch_j
1786 ef2_i(3, 1) = ef2_i(3, 1) - tij_ab(3, 1)*ch_j
1787 ef2_i(1, 2) = ef2_i(1, 2) - tij_ab(1, 2)*ch_j
1788 ef2_i(2, 2) = ef2_i(2, 2) - tij_ab(2, 2)*ch_j
1789 ef2_i(3, 2) = ef2_i(3, 2) - tij_ab(3, 2)*ch_j
1790 ef2_i(1, 3) = ef2_i(1, 3) - tij_ab(1, 3)*ch_j
1791 ef2_i(2, 3) = ef2_i(2, 3) - tij_ab(2, 3)*ch_j
1792 ef2_i(3, 3) = ef2_i(3, 3) - tij_ab(3, 3)*ch_j
1794 ef2_j(1, 1) = ef2_j(1, 1) - tij_ab(1, 1)*ch_i
1795 ef2_j(2, 1) = ef2_j(2, 1) - tij_ab(2, 1)*ch_i
1796 ef2_j(3, 1) = ef2_j(3, 1) - tij_ab(3, 1)*ch_i
1797 ef2_j(1, 2) = ef2_j(1, 2) - tij_ab(1, 2)*ch_i
1798 ef2_j(2, 2) = ef2_j(2, 2) - tij_ab(2, 2)*ch_i
1799 ef2_j(3, 2) = ef2_j(3, 2) - tij_ab(3, 2)*ch_i
1800 ef2_j(1, 3) = ef2_j(1, 3) - tij_ab(1, 3)*ch_i
1801 ef2_j(2, 3) = ef2_j(2, 3) - tij_ab(2, 3)*ch_i
1802 ef2_j(3, 3) = ef2_j(3, 3) - tij_ab(3, 3)*ch_i
1806 IF (task(2, 2))
THEN
1808 tmp = -(dp_i(1)*(tij_ab(1, 1)*dp_j(1) + &
1809 tij_ab(2, 1)*dp_j(2) + &
1810 tij_ab(3, 1)*dp_j(3)) + &
1811 dp_i(2)*(tij_ab(1, 2)*dp_j(1) + &
1812 tij_ab(2, 2)*dp_j(2) + &
1813 tij_ab(3, 2)*dp_j(3)) + &
1814 dp_i(3)*(tij_ab(1, 3)*dp_j(1) + &
1815 tij_ab(2, 3)*dp_j(2) + &
1816 tij_ab(3, 3)*dp_j(3)))
1819 IF (do_forces .OR. do_stress)
THEN
1821 fr(k) = fr(k) + dp_i(1)*(tij_abc(1, 1, k)*dp_j(1) + &
1822 tij_abc(2, 1, k)*dp_j(2) + &
1823 tij_abc(3, 1, k)*dp_j(3)) &
1824 + dp_i(2)*(tij_abc(1, 2, k)*dp_j(1) + &
1825 tij_abc(2, 2, k)*dp_j(2) + &
1826 tij_abc(3, 2, k)*dp_j(3)) &
1827 + dp_i(3)*(tij_abc(1, 3, k)*dp_j(1) + &
1828 tij_abc(2, 3, k)*dp_j(2) + &
1829 tij_abc(3, 3, k)*dp_j(3))
1835 IF (do_efield0)
THEN
1836 ef0_i = ef0_i - (tij_a(1)*dp_j(1) + &
1837 tij_a(2)*dp_j(2) + &
1840 ef0_j = ef0_j + (tij_a(1)*dp_i(1) + &
1841 tij_a(2)*dp_i(2) + &
1845 IF (do_efield1)
THEN
1846 ef1_i(1) = ef1_i(1) + (tij_ab(1, 1)*dp_j(1) + &
1847 tij_ab(2, 1)*dp_j(2) + &
1848 tij_ab(3, 1)*dp_j(3))
1849 ef1_i(2) = ef1_i(2) + (tij_ab(1, 2)*dp_j(1) + &
1850 tij_ab(2, 2)*dp_j(2) + &
1851 tij_ab(3, 2)*dp_j(3))
1852 ef1_i(3) = ef1_i(3) + (tij_ab(1, 3)*dp_j(1) + &
1853 tij_ab(2, 3)*dp_j(2) + &
1854 tij_ab(3, 3)*dp_j(3))
1856 ef1_j(1) = ef1_j(1) + (tij_ab(1, 1)*dp_i(1) + &
1857 tij_ab(2, 1)*dp_i(2) + &
1858 tij_ab(3, 1)*dp_i(3))
1859 ef1_j(2) = ef1_j(2) + (tij_ab(1, 2)*dp_i(1) + &
1860 tij_ab(2, 2)*dp_i(2) + &
1861 tij_ab(3, 2)*dp_i(3))
1862 ef1_j(3) = ef1_j(3) + (tij_ab(1, 3)*dp_i(1) + &
1863 tij_ab(2, 3)*dp_i(2) + &
1864 tij_ab(3, 3)*dp_i(3))
1867 IF (do_efield2)
THEN
1868 ef2_i(1, 1) = ef2_i(1, 1) + (tij_abc(1, 1, 1)*dp_j(1) + &
1869 tij_abc(2, 1, 1)*dp_j(2) + &
1870 tij_abc(3, 1, 1)*dp_j(3))
1871 ef2_i(1, 2) = ef2_i(1, 2) + (tij_abc(1, 1, 2)*dp_j(1) + &
1872 tij_abc(2, 1, 2)*dp_j(2) + &
1873 tij_abc(3, 1, 2)*dp_j(3))
1874 ef2_i(1, 3) = ef2_i(1, 3) + (tij_abc(1, 1, 3)*dp_j(1) + &
1875 tij_abc(2, 1, 3)*dp_j(2) + &
1876 tij_abc(3, 1, 3)*dp_j(3))
1877 ef2_i(2, 1) = ef2_i(2, 1) + (tij_abc(1, 2, 1)*dp_j(1) + &
1878 tij_abc(2, 2, 1)*dp_j(2) + &
1879 tij_abc(3, 2, 1)*dp_j(3))
1880 ef2_i(2, 2) = ef2_i(2, 2) + (tij_abc(1, 2, 2)*dp_j(1) + &
1881 tij_abc(2, 2, 2)*dp_j(2) + &
1882 tij_abc(3, 2, 2)*dp_j(3))
1883 ef2_i(2, 3) = ef2_i(2, 3) + (tij_abc(1, 2, 3)*dp_j(1) + &
1884 tij_abc(2, 2, 3)*dp_j(2) + &
1885 tij_abc(3, 2, 3)*dp_j(3))
1886 ef2_i(3, 1) = ef2_i(3, 1) + (tij_abc(1, 3, 1)*dp_j(1) + &
1887 tij_abc(2, 3, 1)*dp_j(2) + &
1888 tij_abc(3, 3, 1)*dp_j(3))
1889 ef2_i(3, 2) = ef2_i(3, 2) + (tij_abc(1, 3, 2)*dp_j(1) + &
1890 tij_abc(2, 3, 2)*dp_j(2) + &
1891 tij_abc(3, 3, 2)*dp_j(3))
1892 ef2_i(3, 3) = ef2_i(3, 3) + (tij_abc(1, 3, 3)*dp_j(1) + &
1893 tij_abc(2, 3, 3)*dp_j(2) + &
1894 tij_abc(3, 3, 3)*dp_j(3))
1896 ef2_j(1, 1) = ef2_j(1, 1) - (tij_abc(1, 1, 1)*dp_i(1) + &
1897 tij_abc(2, 1, 1)*dp_i(2) + &
1898 tij_abc(3, 1, 1)*dp_i(3))
1899 ef2_j(1, 2) = ef2_j(1, 2) - (tij_abc(1, 1, 2)*dp_i(1) + &
1900 tij_abc(2, 1, 2)*dp_i(2) + &
1901 tij_abc(3, 1, 2)*dp_i(3))
1902 ef2_j(1, 3) = ef2_j(1, 3) - (tij_abc(1, 1, 3)*dp_i(1) + &
1903 tij_abc(2, 1, 3)*dp_i(2) + &
1904 tij_abc(3, 1, 3)*dp_i(3))
1905 ef2_j(2, 1) = ef2_j(2, 1) - (tij_abc(1, 2, 1)*dp_i(1) + &
1906 tij_abc(2, 2, 1)*dp_i(2) + &
1907 tij_abc(3, 2, 1)*dp_i(3))
1908 ef2_j(2, 2) = ef2_j(2, 2) - (tij_abc(1, 2, 2)*dp_i(1) + &
1909 tij_abc(2, 2, 2)*dp_i(2) + &
1910 tij_abc(3, 2, 2)*dp_i(3))
1911 ef2_j(2, 3) = ef2_j(2, 3) - (tij_abc(1, 2, 3)*dp_i(1) + &
1912 tij_abc(2, 2, 3)*dp_i(2) + &
1913 tij_abc(3, 2, 3)*dp_i(3))
1914 ef2_j(3, 1) = ef2_j(3, 1) - (tij_abc(1, 3, 1)*dp_i(1) + &
1915 tij_abc(2, 3, 1)*dp_i(2) + &
1916 tij_abc(3, 3, 1)*dp_i(3))
1917 ef2_j(3, 2) = ef2_j(3, 2) - (tij_abc(1, 3, 2)*dp_i(1) + &
1918 tij_abc(2, 3, 2)*dp_i(2) + &
1919 tij_abc(3, 3, 2)*dp_i(3))
1920 ef2_j(3, 3) = ef2_j(3, 3) - (tij_abc(1, 3, 3)*dp_i(1) + &
1921 tij_abc(2, 3, 3)*dp_i(2) + &
1922 tij_abc(3, 3, 3)*dp_i(3))
1926 IF (task(2, 1))
THEN
1928 tmp = ch_j*(tij_a(1)*dp_i(1) + &
1929 tij_a(2)*dp_i(2) + &
1931 - ch_i*(tij_a(1)*dp_j(1) + &
1932 tij_a(2)*dp_j(2) + &
1934 tmp = tmp - ch_j*(damptij_a(1)*dp_i(1) + &
1935 damptij_a(2)*dp_i(2) + &
1936 damptij_a(3)*dp_i(3)) &
1937 + ch_i*(damptji_a(1)*dp_j(1) + &
1938 damptji_a(2)*dp_j(2) + &
1939 damptji_a(3)*dp_j(3))
1942 IF (do_forces .OR. do_stress)
THEN
1944 fr(k) = fr(k) - ch_j*(tij_ab(1, k)*dp_i(1) + &
1945 tij_ab(2, k)*dp_i(2) + &
1946 tij_ab(3, k)*dp_i(3)) &
1947 + ch_i*(tij_ab(1, k)*dp_j(1) + &
1948 tij_ab(2, k)*dp_j(2) + &
1949 tij_ab(3, k)*dp_j(3))
1950 fr(k) = fr(k) + ch_j*(damptij_ab(1, k)*dp_i(1) + &
1951 damptij_ab(2, k)*dp_i(2) + &
1952 damptij_ab(3, k)*dp_i(3)) &
1953 - ch_i*(damptji_ab(1, k)*dp_j(1) + &
1954 damptji_ab(2, k)*dp_j(2) + &
1955 damptji_ab(3, k)*dp_j(3))
1959 IF (task(3, 3))
THEN
1962 tmp11 = qp_i(1, 1)*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
1963 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
1964 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
1965 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
1966 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
1967 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
1968 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
1969 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
1970 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
1971 tmp21 = qp_i(2, 1)*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
1972 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
1973 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
1974 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
1975 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
1976 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
1977 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
1978 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
1979 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
1980 tmp31 = qp_i(3, 1)*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
1981 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
1982 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
1983 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
1984 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
1985 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
1986 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
1987 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
1988 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
1989 tmp22 = qp_i(2, 2)*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
1990 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
1991 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
1992 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
1993 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
1994 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
1995 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
1996 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
1997 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
1998 tmp32 = qp_i(3, 2)*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
1999 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
2000 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
2001 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
2002 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
2003 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
2004 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
2005 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
2006 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
2007 tmp33 = qp_i(3, 3)*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
2008 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
2009 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
2010 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
2011 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
2012 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
2013 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
2014 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
2015 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
2019 tmp = tmp11 + tmp12 + tmp13 + &
2020 tmp21 + tmp22 + tmp23 + &
2021 tmp31 + tmp32 + tmp33
2023 eloc = eloc +
fac*tmp
2025 IF (do_forces .OR. do_stress)
THEN
2027 tmp11 = qp_i(1, 1)*(tij_abcde(1, 1, 1, 1, k)*qp_j(1, 1) + &
2028 tij_abcde(2, 1, 1, 1, k)*qp_j(2, 1) + &
2029 tij_abcde(3, 1, 1, 1, k)*qp_j(3, 1) + &
2030 tij_abcde(1, 2, 1, 1, k)*qp_j(1, 2) + &
2031 tij_abcde(2, 2, 1, 1, k)*qp_j(2, 2) + &
2032 tij_abcde(3, 2, 1, 1, k)*qp_j(3, 2) + &
2033 tij_abcde(1, 3, 1, 1, k)*qp_j(1, 3) + &
2034 tij_abcde(2, 3, 1, 1, k)*qp_j(2, 3) + &
2035 tij_abcde(3, 3, 1, 1, k)*qp_j(3, 3))
2036 tmp21 = qp_i(2, 1)*(tij_abcde(1, 1, 2, 1, k)*qp_j(1, 1) + &
2037 tij_abcde(2, 1, 2, 1, k)*qp_j(2, 1) + &
2038 tij_abcde(3, 1, 2, 1, k)*qp_j(3, 1) + &
2039 tij_abcde(1, 2, 2, 1, k)*qp_j(1, 2) + &
2040 tij_abcde(2, 2, 2, 1, k)*qp_j(2, 2) + &
2041 tij_abcde(3, 2, 2, 1, k)*qp_j(3, 2) + &
2042 tij_abcde(1, 3, 2, 1, k)*qp_j(1, 3) + &
2043 tij_abcde(2, 3, 2, 1, k)*qp_j(2, 3) + &
2044 tij_abcde(3, 3, 2, 1, k)*qp_j(3, 3))
2045 tmp31 = qp_i(3, 1)*(tij_abcde(1, 1, 3, 1, k)*qp_j(1, 1) + &
2046 tij_abcde(2, 1, 3, 1, k)*qp_j(2, 1) + &
2047 tij_abcde(3, 1, 3, 1, k)*qp_j(3, 1) + &
2048 tij_abcde(1, 2, 3, 1, k)*qp_j(1, 2) + &
2049 tij_abcde(2, 2, 3, 1, k)*qp_j(2, 2) + &
2050 tij_abcde(3, 2, 3, 1, k)*qp_j(3, 2) + &
2051 tij_abcde(1, 3, 3, 1, k)*qp_j(1, 3) + &
2052 tij_abcde(2, 3, 3, 1, k)*qp_j(2, 3) + &
2053 tij_abcde(3, 3, 3, 1, k)*qp_j(3, 3))
2054 tmp22 = qp_i(2, 2)*(tij_abcde(1, 1, 2, 2, k)*qp_j(1, 1) + &
2055 tij_abcde(2, 1, 2, 2, k)*qp_j(2, 1) + &
2056 tij_abcde(3, 1, 2, 2, k)*qp_j(3, 1) + &
2057 tij_abcde(1, 2, 2, 2, k)*qp_j(1, 2) + &
2058 tij_abcde(2, 2, 2, 2, k)*qp_j(2, 2) + &
2059 tij_abcde(3, 2, 2, 2, k)*qp_j(3, 2) + &
2060 tij_abcde(1, 3, 2, 2, k)*qp_j(1, 3) + &
2061 tij_abcde(2, 3, 2, 2, k)*qp_j(2, 3) + &
2062 tij_abcde(3, 3, 2, 2, k)*qp_j(3, 3))
2063 tmp32 = qp_i(3, 2)*(tij_abcde(1, 1, 3, 2, k)*qp_j(1, 1) + &
2064 tij_abcde(2, 1, 3, 2, k)*qp_j(2, 1) + &
2065 tij_abcde(3, 1, 3, 2, k)*qp_j(3, 1) + &
2066 tij_abcde(1, 2, 3, 2, k)*qp_j(1, 2) + &
2067 tij_abcde(2, 2, 3, 2, k)*qp_j(2, 2) + &
2068 tij_abcde(3, 2, 3, 2, k)*qp_j(3, 2) + &
2069 tij_abcde(1, 3, 3, 2, k)*qp_j(1, 3) + &
2070 tij_abcde(2, 3, 3, 2, k)*qp_j(2, 3) + &
2071 tij_abcde(3, 3, 3, 2, k)*qp_j(3, 3))
2072 tmp33 = qp_i(3, 3)*(tij_abcde(1, 1, 3, 3, k)*qp_j(1, 1) + &
2073 tij_abcde(2, 1, 3, 3, k)*qp_j(2, 1) + &
2074 tij_abcde(3, 1, 3, 3, k)*qp_j(3, 1) + &
2075 tij_abcde(1, 2, 3, 3, k)*qp_j(1, 2) + &
2076 tij_abcde(2, 2, 3, 3, k)*qp_j(2, 2) + &
2077 tij_abcde(3, 2, 3, 3, k)*qp_j(3, 2) + &
2078 tij_abcde(1, 3, 3, 3, k)*qp_j(1, 3) + &
2079 tij_abcde(2, 3, 3, 3, k)*qp_j(2, 3) + &
2080 tij_abcde(3, 3, 3, 3, k)*qp_j(3, 3))
2084 fr(k) = fr(k) -
fac*(tmp11 + tmp12 + tmp13 + &
2085 tmp21 + tmp22 + tmp23 + &
2086 tmp31 + tmp32 + tmp33)
2093 IF (do_efield0)
THEN
2094 ef0_i = ef0_i +
fac*(tij_ab(1, 1)*qp_j(1, 1) + &
2095 tij_ab(2, 1)*qp_j(2, 1) + &
2096 tij_ab(3, 1)*qp_j(3, 1) + &
2097 tij_ab(1, 2)*qp_j(1, 2) + &
2098 tij_ab(2, 2)*qp_j(2, 2) + &
2099 tij_ab(3, 2)*qp_j(3, 2) + &
2100 tij_ab(1, 3)*qp_j(1, 3) + &
2101 tij_ab(2, 3)*qp_j(2, 3) + &
2102 tij_ab(3, 3)*qp_j(3, 3))
2104 ef0_j = ef0_j +
fac*(tij_ab(1, 1)*qp_i(1, 1) + &
2105 tij_ab(2, 1)*qp_i(2, 1) + &
2106 tij_ab(3, 1)*qp_i(3, 1) + &
2107 tij_ab(1, 2)*qp_i(1, 2) + &
2108 tij_ab(2, 2)*qp_i(2, 2) + &
2109 tij_ab(3, 2)*qp_i(3, 2) + &
2110 tij_ab(1, 3)*qp_i(1, 3) + &
2111 tij_ab(2, 3)*qp_i(2, 3) + &
2112 tij_ab(3, 3)*qp_i(3, 3))
2115 IF (do_efield1)
THEN
2116 ef1_i(1) = ef1_i(1) -
fac*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
2117 tij_abc(2, 1, 1)*qp_j(2, 1) + &
2118 tij_abc(3, 1, 1)*qp_j(3, 1) + &
2119 tij_abc(1, 2, 1)*qp_j(1, 2) + &
2120 tij_abc(2, 2, 1)*qp_j(2, 2) + &
2121 tij_abc(3, 2, 1)*qp_j(3, 2) + &
2122 tij_abc(1, 3, 1)*qp_j(1, 3) + &
2123 tij_abc(2, 3, 1)*qp_j(2, 3) + &
2124 tij_abc(3, 3, 1)*qp_j(3, 3))
2125 ef1_i(2) = ef1_i(2) -
fac*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
2126 tij_abc(2, 1, 2)*qp_j(2, 1) + &
2127 tij_abc(3, 1, 2)*qp_j(3, 1) + &
2128 tij_abc(1, 2, 2)*qp_j(1, 2) + &
2129 tij_abc(2, 2, 2)*qp_j(2, 2) + &
2130 tij_abc(3, 2, 2)*qp_j(3, 2) + &
2131 tij_abc(1, 3, 2)*qp_j(1, 3) + &
2132 tij_abc(2, 3, 2)*qp_j(2, 3) + &
2133 tij_abc(3, 3, 2)*qp_j(3, 3))
2134 ef1_i(3) = ef1_i(3) -
fac*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
2135 tij_abc(2, 1, 3)*qp_j(2, 1) + &
2136 tij_abc(3, 1, 3)*qp_j(3, 1) + &
2137 tij_abc(1, 2, 3)*qp_j(1, 2) + &
2138 tij_abc(2, 2, 3)*qp_j(2, 2) + &
2139 tij_abc(3, 2, 3)*qp_j(3, 2) + &
2140 tij_abc(1, 3, 3)*qp_j(1, 3) + &
2141 tij_abc(2, 3, 3)*qp_j(2, 3) + &
2142 tij_abc(3, 3, 3)*qp_j(3, 3))
2144 ef1_j(1) = ef1_j(1) +
fac*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
2145 tij_abc(2, 1, 1)*qp_i(2, 1) + &
2146 tij_abc(3, 1, 1)*qp_i(3, 1) + &
2147 tij_abc(1, 2, 1)*qp_i(1, 2) + &
2148 tij_abc(2, 2, 1)*qp_i(2, 2) + &
2149 tij_abc(3, 2, 1)*qp_i(3, 2) + &
2150 tij_abc(1, 3, 1)*qp_i(1, 3) + &
2151 tij_abc(2, 3, 1)*qp_i(2, 3) + &
2152 tij_abc(3, 3, 1)*qp_i(3, 3))
2153 ef1_j(2) = ef1_j(2) +
fac*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
2154 tij_abc(2, 1, 2)*qp_i(2, 1) + &
2155 tij_abc(3, 1, 2)*qp_i(3, 1) + &
2156 tij_abc(1, 2, 2)*qp_i(1, 2) + &
2157 tij_abc(2, 2, 2)*qp_i(2, 2) + &
2158 tij_abc(3, 2, 2)*qp_i(3, 2) + &
2159 tij_abc(1, 3, 2)*qp_i(1, 3) + &
2160 tij_abc(2, 3, 2)*qp_i(2, 3) + &
2161 tij_abc(3, 3, 2)*qp_i(3, 3))
2162 ef1_j(3) = ef1_j(3) +
fac*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
2163 tij_abc(2, 1, 3)*qp_i(2, 1) + &
2164 tij_abc(3, 1, 3)*qp_i(3, 1) + &
2165 tij_abc(1, 2, 3)*qp_i(1, 2) + &
2166 tij_abc(2, 2, 3)*qp_i(2, 2) + &
2167 tij_abc(3, 2, 3)*qp_i(3, 2) + &
2168 tij_abc(1, 3, 3)*qp_i(1, 3) + &
2169 tij_abc(2, 3, 3)*qp_i(2, 3) + &
2170 tij_abc(3, 3, 3)*qp_i(3, 3))
2173 IF (do_efield2)
THEN
2174 tmp11 =
fac*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
2175 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
2176 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
2177 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
2178 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
2179 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
2180 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
2181 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
2182 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
2183 tmp12 =
fac*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
2184 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
2185 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
2186 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
2187 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
2188 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
2189 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
2190 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
2191 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
2192 tmp13 =
fac*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
2193 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
2194 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
2195 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
2196 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
2197 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
2198 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
2199 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
2200 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
2201 tmp22 =
fac*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
2202 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
2203 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
2204 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
2205 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
2206 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
2207 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
2208 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
2209 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
2210 tmp23 =
fac*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
2211 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
2212 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
2213 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
2214 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
2215 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
2216 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
2217 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
2218 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
2219 tmp33 =
fac*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
2220 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
2221 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
2222 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
2223 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
2224 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
2225 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
2226 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
2227 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
2229 ef2_i(1, 1) = ef2_i(1, 1) - tmp11
2230 ef2_i(1, 2) = ef2_i(1, 2) - tmp12
2231 ef2_i(1, 3) = ef2_i(1, 3) - tmp13
2232 ef2_i(2, 1) = ef2_i(2, 1) - tmp12
2233 ef2_i(2, 2) = ef2_i(2, 2) - tmp22
2234 ef2_i(2, 3) = ef2_i(2, 3) - tmp23
2235 ef2_i(3, 1) = ef2_i(3, 1) - tmp13
2236 ef2_i(3, 2) = ef2_i(3, 2) - tmp23
2237 ef2_i(3, 3) = ef2_i(3, 3) - tmp33
2239 tmp11 =
fac*(tij_abcd(1, 1, 1, 1)*qp_i(1, 1) + &
2240 tij_abcd(2, 1, 1, 1)*qp_i(2, 1) + &
2241 tij_abcd(3, 1, 1, 1)*qp_i(3, 1) + &
2242 tij_abcd(1, 2, 1, 1)*qp_i(1, 2) + &
2243 tij_abcd(2, 2, 1, 1)*qp_i(2, 2) + &
2244 tij_abcd(3, 2, 1, 1)*qp_i(3, 2) + &
2245 tij_abcd(1, 3, 1, 1)*qp_i(1, 3) + &
2246 tij_abcd(2, 3, 1, 1)*qp_i(2, 3) + &
2247 tij_abcd(3, 3, 1, 1)*qp_i(3, 3))
2248 tmp12 =
fac*(tij_abcd(1, 1, 1, 2)*qp_i(1, 1) + &
2249 tij_abcd(2, 1, 1, 2)*qp_i(2, 1) + &
2250 tij_abcd(3, 1, 1, 2)*qp_i(3, 1) + &
2251 tij_abcd(1, 2, 1, 2)*qp_i(1, 2) + &
2252 tij_abcd(2, 2, 1, 2)*qp_i(2, 2) + &
2253 tij_abcd(3, 2, 1, 2)*qp_i(3, 2) + &
2254 tij_abcd(1, 3, 1, 2)*qp_i(1, 3) + &
2255 tij_abcd(2, 3, 1, 2)*qp_i(2, 3) + &
2256 tij_abcd(3, 3, 1, 2)*qp_i(3, 3))
2257 tmp13 =
fac*(tij_abcd(1, 1, 1, 3)*qp_i(1, 1) + &
2258 tij_abcd(2, 1, 1, 3)*qp_i(2, 1) + &
2259 tij_abcd(3, 1, 1, 3)*qp_i(3, 1) + &
2260 tij_abcd(1, 2, 1, 3)*qp_i(1, 2) + &
2261 tij_abcd(2, 2, 1, 3)*qp_i(2, 2) + &
2262 tij_abcd(3, 2, 1, 3)*qp_i(3, 2) + &
2263 tij_abcd(1, 3, 1, 3)*qp_i(1, 3) + &
2264 tij_abcd(2, 3, 1, 3)*qp_i(2, 3) + &
2265 tij_abcd(3, 3, 1, 3)*qp_i(3, 3))
2266 tmp22 =
fac*(tij_abcd(1, 1, 2, 2)*qp_i(1, 1) + &
2267 tij_abcd(2, 1, 2, 2)*qp_i(2, 1) + &
2268 tij_abcd(3, 1, 2, 2)*qp_i(3, 1) + &
2269 tij_abcd(1, 2, 2, 2)*qp_i(1, 2) + &
2270 tij_abcd(2, 2, 2, 2)*qp_i(2, 2) + &
2271 tij_abcd(3, 2, 2, 2)*qp_i(3, 2) + &
2272 tij_abcd(1, 3, 2, 2)*qp_i(1, 3) + &
2273 tij_abcd(2, 3, 2, 2)*qp_i(2, 3) + &
2274 tij_abcd(3, 3, 2, 2)*qp_i(3, 3))
2275 tmp23 =
fac*(tij_abcd(1, 1, 2, 3)*qp_i(1, 1) + &
2276 tij_abcd(2, 1, 2, 3)*qp_i(2, 1) + &
2277 tij_abcd(3, 1, 2, 3)*qp_i(3, 1) + &
2278 tij_abcd(1, 2, 2, 3)*qp_i(1, 2) + &
2279 tij_abcd(2, 2, 2, 3)*qp_i(2, 2) + &
2280 tij_abcd(3, 2, 2, 3)*qp_i(3, 2) + &
2281 tij_abcd(1, 3, 2, 3)*qp_i(1, 3) + &
2282 tij_abcd(2, 3, 2, 3)*qp_i(2, 3) + &
2283 tij_abcd(3, 3, 2, 3)*qp_i(3, 3))
2284 tmp33 =
fac*(tij_abcd(1, 1, 3, 3)*qp_i(1, 1) + &
2285 tij_abcd(2, 1, 3, 3)*qp_i(2, 1) + &
2286 tij_abcd(3, 1, 3, 3)*qp_i(3, 1) + &
2287 tij_abcd(1, 2, 3, 3)*qp_i(1, 2) + &
2288 tij_abcd(2, 2, 3, 3)*qp_i(2, 2) + &
2289 tij_abcd(3, 2, 3, 3)*qp_i(3, 2) + &
2290 tij_abcd(1, 3, 3, 3)*qp_i(1, 3) + &
2291 tij_abcd(2, 3, 3, 3)*qp_i(2, 3) + &
2292 tij_abcd(3, 3, 3, 3)*qp_i(3, 3))
2294 ef2_j(1, 1) = ef2_j(1, 1) - tmp11
2295 ef2_j(1, 2) = ef2_j(1, 2) - tmp12
2296 ef2_j(1, 3) = ef2_j(1, 3) - tmp13
2297 ef2_j(2, 1) = ef2_j(2, 1) - tmp12
2298 ef2_j(2, 2) = ef2_j(2, 2) - tmp22
2299 ef2_j(2, 3) = ef2_j(2, 3) - tmp23
2300 ef2_j(3, 1) = ef2_j(3, 1) - tmp13
2301 ef2_j(3, 2) = ef2_j(3, 2) - tmp23
2302 ef2_j(3, 3) = ef2_j(3, 3) - tmp33
2306 IF (task(3, 2))
THEN
2310 tmp_ij = dp_i(1)*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
2311 tij_abc(2, 1, 1)*qp_j(2, 1) + &
2312 tij_abc(3, 1, 1)*qp_j(3, 1) + &
2313 tij_abc(1, 2, 1)*qp_j(1, 2) + &
2314 tij_abc(2, 2, 1)*qp_j(2, 2) + &
2315 tij_abc(3, 2, 1)*qp_j(3, 2) + &
2316 tij_abc(1, 3, 1)*qp_j(1, 3) + &
2317 tij_abc(2, 3, 1)*qp_j(2, 3) + &
2318 tij_abc(3, 3, 1)*qp_j(3, 3)) + &
2319 dp_i(2)*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
2320 tij_abc(2, 1, 2)*qp_j(2, 1) + &
2321 tij_abc(3, 1, 2)*qp_j(3, 1) + &
2322 tij_abc(1, 2, 2)*qp_j(1, 2) + &
2323 tij_abc(2, 2, 2)*qp_j(2, 2) + &
2324 tij_abc(3, 2, 2)*qp_j(3, 2) + &
2325 tij_abc(1, 3, 2)*qp_j(1, 3) + &
2326 tij_abc(2, 3, 2)*qp_j(2, 3) + &
2327 tij_abc(3, 3, 2)*qp_j(3, 3)) + &
2328 dp_i(3)*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
2329 tij_abc(2, 1, 3)*qp_j(2, 1) + &
2330 tij_abc(3, 1, 3)*qp_j(3, 1) + &
2331 tij_abc(1, 2, 3)*qp_j(1, 2) + &
2332 tij_abc(2, 2, 3)*qp_j(2, 2) + &
2333 tij_abc(3, 2, 3)*qp_j(3, 2) + &
2334 tij_abc(1, 3, 3)*qp_j(1, 3) + &
2335 tij_abc(2, 3, 3)*qp_j(2, 3) + &
2336 tij_abc(3, 3, 3)*qp_j(3, 3))
2339 tmp_ji = dp_j(1)*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
2340 tij_abc(2, 1, 1)*qp_i(2, 1) + &
2341 tij_abc(3, 1, 1)*qp_i(3, 1) + &
2342 tij_abc(1, 2, 1)*qp_i(1, 2) + &
2343 tij_abc(2, 2, 1)*qp_i(2, 2) + &
2344 tij_abc(3, 2, 1)*qp_i(3, 2) + &
2345 tij_abc(1, 3, 1)*qp_i(1, 3) + &
2346 tij_abc(2, 3, 1)*qp_i(2, 3) + &
2347 tij_abc(3, 3, 1)*qp_i(3, 3)) + &
2348 dp_j(2)*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
2349 tij_abc(2, 1, 2)*qp_i(2, 1) + &
2350 tij_abc(3, 1, 2)*qp_i(3, 1) + &
2351 tij_abc(1, 2, 2)*qp_i(1, 2) + &
2352 tij_abc(2, 2, 2)*qp_i(2, 2) + &
2353 tij_abc(3, 2, 2)*qp_i(3, 2) + &
2354 tij_abc(1, 3, 2)*qp_i(1, 3) + &
2355 tij_abc(2, 3, 2)*qp_i(2, 3) + &
2356 tij_abc(3, 3, 2)*qp_i(3, 3)) + &
2357 dp_j(3)*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
2358 tij_abc(2, 1, 3)*qp_i(2, 1) + &
2359 tij_abc(3, 1, 3)*qp_i(3, 1) + &
2360 tij_abc(1, 2, 3)*qp_i(1, 2) + &
2361 tij_abc(2, 2, 3)*qp_i(2, 2) + &
2362 tij_abc(3, 2, 3)*qp_i(3, 2) + &
2363 tij_abc(1, 3, 3)*qp_i(1, 3) + &
2364 tij_abc(2, 3, 3)*qp_i(2, 3) + &
2365 tij_abc(3, 3, 3)*qp_i(3, 3))
2367 tmp =
fac*(tmp_ij - tmp_ji)
2369 IF (do_forces .OR. do_stress)
THEN
2372 tmp_ij = dp_i(1)*(tij_abcd(1, 1, 1, k)*qp_j(1, 1) + &
2373 tij_abcd(2, 1, 1, k)*qp_j(2, 1) + &
2374 tij_abcd(3, 1, 1, k)*qp_j(3, 1) + &
2375 tij_abcd(1, 2, 1, k)*qp_j(1, 2) + &
2376 tij_abcd(2, 2, 1, k)*qp_j(2, 2) + &
2377 tij_abcd(3, 2, 1, k)*qp_j(3, 2) + &
2378 tij_abcd(1, 3, 1, k)*qp_j(1, 3) + &
2379 tij_abcd(2, 3, 1, k)*qp_j(2, 3) + &
2380 tij_abcd(3, 3, 1, k)*qp_j(3, 3)) + &
2381 dp_i(2)*(tij_abcd(1, 1, 2, k)*qp_j(1, 1) + &
2382 tij_abcd(2, 1, 2, k)*qp_j(2, 1) + &
2383 tij_abcd(3, 1, 2, k)*qp_j(3, 1) + &
2384 tij_abcd(1, 2, 2, k)*qp_j(1, 2) + &
2385 tij_abcd(2, 2, 2, k)*qp_j(2, 2) + &
2386 tij_abcd(3, 2, 2, k)*qp_j(3, 2) + &
2387 tij_abcd(1, 3, 2, k)*qp_j(1, 3) + &
2388 tij_abcd(2, 3, 2, k)*qp_j(2, 3) + &
2389 tij_abcd(3, 3, 2, k)*qp_j(3, 3)) + &
2390 dp_i(3)*(tij_abcd(1, 1, 3, k)*qp_j(1, 1) + &
2391 tij_abcd(2, 1, 3, k)*qp_j(2, 1) + &
2392 tij_abcd(3, 1, 3, k)*qp_j(3, 1) + &
2393 tij_abcd(1, 2, 3, k)*qp_j(1, 2) + &
2394 tij_abcd(2, 2, 3, k)*qp_j(2, 2) + &
2395 tij_abcd(3, 2, 3, k)*qp_j(3, 2) + &
2396 tij_abcd(1, 3, 3, k)*qp_j(1, 3) + &
2397 tij_abcd(2, 3, 3, k)*qp_j(2, 3) + &
2398 tij_abcd(3, 3, 3, k)*qp_j(3, 3))
2401 tmp_ji = dp_j(1)*(tij_abcd(1, 1, 1, k)*qp_i(1, 1) + &
2402 tij_abcd(2, 1, 1, k)*qp_i(2, 1) + &
2403 tij_abcd(3, 1, 1, k)*qp_i(3, 1) + &
2404 tij_abcd(1, 2, 1, k)*qp_i(1, 2) + &
2405 tij_abcd(2, 2, 1, k)*qp_i(2, 2) + &
2406 tij_abcd(3, 2, 1, k)*qp_i(3, 2) + &
2407 tij_abcd(1, 3, 1, k)*qp_i(1, 3) + &
2408 tij_abcd(2, 3, 1, k)*qp_i(2, 3) + &
2409 tij_abcd(3, 3, 1, k)*qp_i(3, 3)) + &
2410 dp_j(2)*(tij_abcd(1, 1, 2, k)*qp_i(1, 1) + &
2411 tij_abcd(2, 1, 2, k)*qp_i(2, 1) + &
2412 tij_abcd(3, 1, 2, k)*qp_i(3, 1) + &
2413 tij_abcd(1, 2, 2, k)*qp_i(1, 2) + &
2414 tij_abcd(2, 2, 2, k)*qp_i(2, 2) + &
2415 tij_abcd(3, 2, 2, k)*qp_i(3, 2) + &
2416 tij_abcd(1, 3, 2, k)*qp_i(1, 3) + &
2417 tij_abcd(2, 3, 2, k)*qp_i(2, 3) + &
2418 tij_abcd(3, 3, 2, k)*qp_i(3, 3)) + &
2419 dp_j(3)*(tij_abcd(1, 1, 3, k)*qp_i(1, 1) + &
2420 tij_abcd(2, 1, 3, k)*qp_i(2, 1) + &
2421 tij_abcd(3, 1, 3, k)*qp_i(3, 1) + &
2422 tij_abcd(1, 2, 3, k)*qp_i(1, 2) + &
2423 tij_abcd(2, 2, 3, k)*qp_i(2, 2) + &
2424 tij_abcd(3, 2, 3, k)*qp_i(3, 2) + &
2425 tij_abcd(1, 3, 3, k)*qp_i(1, 3) + &
2426 tij_abcd(2, 3, 3, k)*qp_i(2, 3) + &
2427 tij_abcd(3, 3, 3, k)*qp_i(3, 3))
2429 fr(k) = fr(k) -
fac*(tmp_ij - tmp_ji)
2433 IF (task(3, 1))
THEN
2438 tmp_ij = ch_i*(tij_ab(1, 1)*qp_j(1, 1) + &
2439 tij_ab(2, 1)*qp_j(2, 1) + &
2440 tij_ab(3, 1)*qp_j(3, 1) + &
2441 tij_ab(1, 2)*qp_j(1, 2) + &
2442 tij_ab(2, 2)*qp_j(2, 2) + &
2443 tij_ab(3, 2)*qp_j(3, 2) + &
2444 tij_ab(1, 3)*qp_j(1, 3) + &
2445 tij_ab(2, 3)*qp_j(2, 3) + &
2446 tij_ab(3, 3)*qp_j(3, 3))
2449 tmp_ji = ch_j*(tij_ab(1, 1)*qp_i(1, 1) + &
2450 tij_ab(2, 1)*qp_i(2, 1) + &
2451 tij_ab(3, 1)*qp_i(3, 1) + &
2452 tij_ab(1, 2)*qp_i(1, 2) + &
2453 tij_ab(2, 2)*qp_i(2, 2) + &
2454 tij_ab(3, 2)*qp_i(3, 2) + &
2455 tij_ab(1, 3)*qp_i(1, 3) + &
2456 tij_ab(2, 3)*qp_i(2, 3) + &
2457 tij_ab(3, 3)*qp_i(3, 3))
2459 eloc = eloc +
fac*(tmp_ij + tmp_ji)
2460 IF (do_forces .OR. do_stress)
THEN
2463 tmp_ij = ch_i*(tij_abc(1, 1, k)*qp_j(1, 1) + &
2464 tij_abc(2, 1, k)*qp_j(2, 1) + &
2465 tij_abc(3, 1, k)*qp_j(3, 1) + &
2466 tij_abc(1, 2, k)*qp_j(1, 2) + &
2467 tij_abc(2, 2, k)*qp_j(2, 2) + &
2468 tij_abc(3, 2, k)*qp_j(3, 2) + &
2469 tij_abc(1, 3, k)*qp_j(1, 3) + &
2470 tij_abc(2, 3, k)*qp_j(2, 3) + &
2471 tij_abc(3, 3, k)*qp_j(3, 3))
2474 tmp_ji = ch_j*(tij_abc(1, 1, k)*qp_i(1, 1) + &
2475 tij_abc(2, 1, k)*qp_i(2, 1) + &
2476 tij_abc(3, 1, k)*qp_i(3, 1) + &
2477 tij_abc(1, 2, k)*qp_i(1, 2) + &
2478 tij_abc(2, 2, k)*qp_i(2, 2) + &
2479 tij_abc(3, 2, k)*qp_i(3, 2) + &
2480 tij_abc(1, 3, k)*qp_i(1, 3) + &
2481 tij_abc(2, 3, k)*qp_i(2, 3) + &
2482 tij_abc(3, 3, k)*qp_i(3, 3))
2484 fr(k) = fr(k) -
fac*(tmp_ij + tmp_ji)
2488 energy = energy + eloc
2490 forces(1, atom_a) = forces(1, atom_a) - fr(1)
2491 forces(2, atom_a) = forces(2, atom_a) - fr(2)
2492 forces(3, atom_a) = forces(3, atom_a) - fr(3)
2493 forces(1, atom_b) = forces(1, atom_b) + fr(1)
2494 forces(2, atom_b) = forces(2, atom_b) + fr(2)
2495 forces(3, atom_b) = forces(3, atom_b) + fr(3)
2500 IF (do_efield0)
THEN
2501 efield0(atom_a) = efield0(atom_a) + ef0_j
2503 efield0(atom_b) = efield0(atom_b) + ef0_i
2506 IF (do_efield1)
THEN
2507 efield1(1, atom_a) = efield1(1, atom_a) + ef1_j(1)
2508 efield1(2, atom_a) = efield1(2, atom_a) + ef1_j(2)
2509 efield1(3, atom_a) = efield1(3, atom_a) + ef1_j(3)
2511 efield1(1, atom_b) = efield1(1, atom_b) + ef1_i(1)
2512 efield1(2, atom_b) = efield1(2, atom_b) + ef1_i(2)
2513 efield1(3, atom_b) = efield1(3, atom_b) + ef1_i(3)
2516 IF (do_efield2)
THEN
2517 efield2(1, atom_a) = efield2(1, atom_a) + ef2_j(1, 1)
2518 efield2(2, atom_a) = efield2(2, atom_a) + ef2_j(1, 2)
2519 efield2(3, atom_a) = efield2(3, atom_a) + ef2_j(1, 3)
2520 efield2(4, atom_a) = efield2(4, atom_a) + ef2_j(2, 1)
2521 efield2(5, atom_a) = efield2(5, atom_a) + ef2_j(2, 2)
2522 efield2(6, atom_a) = efield2(6, atom_a) + ef2_j(2, 3)
2523 efield2(7, atom_a) = efield2(7, atom_a) + ef2_j(3, 1)
2524 efield2(8, atom_a) = efield2(8, atom_a) + ef2_j(3, 2)
2525 efield2(9, atom_a) = efield2(9, atom_a) + ef2_j(3, 3)
2527 efield2(1, atom_b) = efield2(1, atom_b) + ef2_i(1, 1)
2528 efield2(2, atom_b) = efield2(2, atom_b) + ef2_i(1, 2)
2529 efield2(3, atom_b) = efield2(3, atom_b) + ef2_i(1, 3)
2530 efield2(4, atom_b) = efield2(4, atom_b) + ef2_i(2, 1)
2531 efield2(5, atom_b) = efield2(5, atom_b) + ef2_i(2, 2)
2532 efield2(6, atom_b) = efield2(6, atom_b) + ef2_i(2, 3)
2533 efield2(7, atom_b) = efield2(7, atom_b) + ef2_i(3, 1)
2534 efield2(8, atom_b) = efield2(8, atom_b) + ef2_i(3, 2)
2535 efield2(9, atom_b) = efield2(9, atom_b) + ef2_i(3, 3)
2539 ptens11 = ptens11 + rab(1)*fr(1)
2540 ptens21 = ptens21 + rab(2)*fr(1)
2541 ptens31 = ptens31 + rab(3)*fr(1)
2542 ptens12 = ptens12 + rab(1)*fr(2)
2543 ptens22 = ptens22 + rab(2)*fr(2)
2544 ptens32 = ptens32 + rab(3)*fr(2)
2545 ptens13 = ptens13 + rab(1)*fr(3)
2546 ptens23 = ptens23 + rab(2)*fr(3)
2547 ptens33 = ptens33 + rab(3)*fr(3)
2553 END DO kind_group_loop
2556 pv(1, 1) = pv(1, 1) + ptens11
2557 pv(1, 2) = pv(1, 2) + (ptens12 + ptens21)*0.5_dp
2558 pv(1, 3) = pv(1, 3) + (ptens13 + ptens31)*0.5_dp
2560 pv(2, 2) = pv(2, 2) + ptens22
2561 pv(2, 3) = pv(2, 3) + (ptens23 + ptens32)*0.5_dp
2564 pv(3, 3) = pv(3, 3) + ptens33
2567 CALL timestop(handle)
2568 END SUBROUTINE ewald_multipole_sr
2592 SUBROUTINE ewald_multipole_bonded(nonbond_env, particle_set, ewald_env, &
2593 cell, energy, task, do_forces, do_efield, do_stress, charges, &
2594 dipoles, quadrupoles, forces, pv, efield0, efield1, efield2)
2600 REAL(kind=
dp),
INTENT(INOUT) :: energy
2601 LOGICAL,
DIMENSION(3, 3),
INTENT(IN) :: task
2602 LOGICAL,
INTENT(IN) :: do_forces, do_efield, do_stress
2603 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: charges
2604 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
2605 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
2606 POINTER :: quadrupoles
2607 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT), &
2608 OPTIONAL :: forces, pv
2609 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0
2610 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: efield1, efield2
2612 CHARACTER(len=*),
PARAMETER :: routinen =
'ewald_multipole_bonded'
2614 INTEGER :: a, atom_a, atom_b, b, c, d, e, handle, &
2615 i, iend, igrp, ilist, ipair, istart, &
2617 INTEGER,
DIMENSION(:, :),
POINTER ::
list
2618 LOGICAL :: do_efield0, do_efield1, do_efield2, &
2620 REAL(kind=
dp) :: alpha, ch_i, ch_j, ef0_i, ef0_j, eloc,
fac, fac_ij, ir, irab2, ptens11, &
2621 ptens12, ptens13, ptens21, ptens22, ptens23, ptens31, ptens32, ptens33, r, rab2, tij, &
2622 tmp, tmp1, tmp11, tmp12, tmp13, tmp2, tmp21, tmp22, tmp23, tmp31, tmp32, tmp33, tmp_ij, &
2624 REAL(kind=
dp),
DIMENSION(0:5) :: f
2625 REAL(kind=
dp),
DIMENSION(3) :: dp_i, dp_j, ef1_i, ef1_j, fr, rab, tij_a
2626 REAL(kind=
dp),
DIMENSION(3, 3) :: ef2_i, ef2_j, qp_i, qp_j, tij_ab
2627 REAL(kind=
dp),
DIMENSION(3, 3, 3) :: tij_abc
2628 REAL(kind=
dp),
DIMENSION(3, 3, 3, 3) :: tij_abcd
2629 REAL(kind=
dp),
DIMENSION(3, 3, 3, 3, 3) :: tij_abcde
2633 CALL timeset(routinen, handle)
2634 do_efield0 = do_efield .AND.
ASSOCIATED(efield0)
2635 do_efield1 = do_efield .AND.
ASSOCIATED(efield1)
2636 do_efield2 = do_efield .AND.
ASSOCIATED(efield2)
2638 ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp
2639 ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp
2640 ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp
2646 lists:
DO ilist = 1, nonbonded%nlists
2647 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
2648 nscale = neighbor_kind_pair%nscale
2649 IF (nscale == 0) cycle
2650 list => neighbor_kind_pair%list
2651 kind_group_loop:
DO igrp = 1, neighbor_kind_pair%ngrp_kind
2652 istart = neighbor_kind_pair%grp_kind_start(igrp)
2653 IF (istart > nscale) cycle
2654 iend = min(neighbor_kind_pair%grp_kind_end(igrp), nscale)
2655 pairs:
DO ipair = istart, iend
2657 fac_ij = -1.0_dp + neighbor_kind_pair%ei_scale(ipair)
2658 IF (fac_ij >= 0) cycle
2660 atom_a =
list(1, ipair)
2661 atom_b =
list(2, ipair)
2663 rab = particle_set(atom_b)%r - particle_set(atom_a)%r
2664 rab =
pbc(rab, cell)
2665 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
2669 tij_a = huge(0.0_dp)
2670 tij_ab = huge(0.0_dp)
2671 tij_abc = huge(0.0_dp)
2672 tij_abcd = huge(0.0_dp)
2673 tij_abcde = huge(0.0_dp)
2679 IF (debug_this_module .AND. debug_r_space .AND. (.NOT. debug_g_space))
THEN
2683 f(0) = erf(alpha*r)*ir
2689 f(i) = irab2*(f(i - 1) - tmp*((2.0_dp*alpha**2)**i)/(
fac*alpha))
2694 force_eval = do_stress
2695 IF (task(1, 1))
THEN
2697 force_eval = do_forces .OR. do_efield1
2699 IF (task(2, 2)) force_eval = force_eval .OR. do_efield0
2700 IF (task(1, 2) .OR. force_eval)
THEN
2701 force_eval = do_stress
2702 tij_a = -rab*f(1)*fac_ij
2703 IF (task(1, 2)) force_eval = force_eval .OR. do_forces
2705 IF (task(1, 1)) force_eval = force_eval .OR. do_efield2
2706 IF (task(3, 3)) force_eval = force_eval .OR. do_efield0
2707 IF (task(2, 2) .OR. task(3, 1) .OR. force_eval)
THEN
2708 force_eval = do_stress
2711 tmp = rab(a)*rab(b)*fac_ij
2712 tij_ab(a, b) = 3.0_dp*tmp*f(2)
2713 IF (a == b) tij_ab(a, b) = tij_ab(a, b) - f(1)*fac_ij
2716 IF (task(2, 2) .OR. task(3, 1)) force_eval = force_eval .OR. do_forces
2718 IF (task(2, 2)) force_eval = force_eval .OR. do_efield2
2719 IF (task(3, 3)) force_eval = force_eval .OR. do_efield1
2720 IF (task(3, 2) .OR. force_eval)
THEN
2721 force_eval = do_stress
2725 tmp = rab(a)*rab(b)*rab(c)*fac_ij
2726 tij_abc(a, b, c) = -15.0_dp*tmp*f(3)
2727 tmp = 3.0_dp*f(2)*fac_ij
2728 IF (a == b) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(c)
2729 IF (a == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(b)
2730 IF (b == c) tij_abc(a, b, c) = tij_abc(a, b, c) + tmp*rab(a)
2734 IF (task(3, 2)) force_eval = force_eval .OR. do_forces
2736 IF (task(3, 3) .OR. force_eval)
THEN
2737 force_eval = do_stress
2742 tmp = rab(a)*rab(b)*rab(c)*rab(d)*fac_ij
2743 tij_abcd(a, b, c, d) = 105.0_dp*tmp*f(4)
2744 tmp1 = 15.0_dp*f(3)*fac_ij
2745 tmp2 = 3.0_dp*f(2)*fac_ij
2747 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(c)*rab(d)
2748 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
2751 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(d)
2752 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
2754 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(b)*rab(c)
2756 tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(d)
2757 IF (a == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) + tmp2
2759 IF (b == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(c)
2760 IF (c == d) tij_abcd(a, b, c, d) = tij_abcd(a, b, c, d) - tmp1*rab(a)*rab(b)
2765 IF (task(3, 3)) force_eval = force_eval .OR. do_forces
2767 IF (force_eval)
THEN
2768 force_eval = do_stress
2774 tmp = rab(a)*rab(b)*rab(c)*rab(d)*rab(e)*fac_ij
2775 tij_abcde(a, b, c, d, e) = -945.0_dp*tmp*f(5)
2776 tmp1 = 105.0_dp*f(4)*fac_ij
2777 tmp2 = 15.0_dp*f(3)*fac_ij
2779 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(c)*rab(d)*rab(e)
2780 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
2781 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
2782 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
2785 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(d)*rab(e)
2786 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
2787 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
2788 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
2791 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(e)
2792 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(e)
2793 IF (b == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
2794 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
2797 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(b)*rab(c)*rab(d)
2798 IF (b == c) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(d)
2799 IF (b == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(c)
2800 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(b)
2803 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(d)*rab(e)
2804 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
2807 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(e)
2808 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
2811 tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(c)*rab(d)
2812 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) - tmp2*rab(a)
2814 IF (c == d) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(e)
2815 IF (c == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(d)
2816 IF (d == e) tij_abcde(a, b, c, d, e) = tij_abcde(a, b, c, d, e) + tmp1*rab(a)*rab(b)*rab(c)
2834 IF (debug_this_module)
THEN
2842 IF (any(task(1, :)))
THEN
2843 ch_j = charges(atom_a)
2844 ch_i = charges(atom_b)
2846 IF (any(task(2, :)))
THEN
2847 dp_j = dipoles(:, atom_a)
2848 dp_i = dipoles(:, atom_b)
2850 IF (any(task(3, :)))
THEN
2851 qp_j = quadrupoles(:, :, atom_a)
2852 qp_i = quadrupoles(:, :, atom_b)
2854 IF (task(1, 1))
THEN
2856 eloc = eloc + ch_i*tij*ch_j
2858 IF (do_forces .OR. do_stress)
THEN
2859 fr(1) = fr(1) - ch_j*tij_a(1)*ch_i
2860 fr(2) = fr(2) - ch_j*tij_a(2)*ch_i
2861 fr(3) = fr(3) - ch_j*tij_a(3)*ch_i
2866 IF (do_efield0)
THEN
2867 ef0_i = ef0_i + tij*ch_j
2869 ef0_j = ef0_j + tij*ch_i
2872 IF (do_efield1)
THEN
2873 ef1_i(1) = ef1_i(1) - tij_a(1)*ch_j
2874 ef1_i(2) = ef1_i(2) - tij_a(2)*ch_j
2875 ef1_i(3) = ef1_i(3) - tij_a(3)*ch_j
2877 ef1_j(1) = ef1_j(1) + tij_a(1)*ch_i
2878 ef1_j(2) = ef1_j(2) + tij_a(2)*ch_i
2879 ef1_j(3) = ef1_j(3) + tij_a(3)*ch_i
2884 IF (do_efield2)
THEN
2885 ef2_i(1, 1) = ef2_i(1, 1) - tij_ab(1, 1)*ch_j
2886 ef2_i(2, 1) = ef2_i(2, 1) - tij_ab(2, 1)*ch_j
2887 ef2_i(3, 1) = ef2_i(3, 1) - tij_ab(3, 1)*ch_j
2888 ef2_i(1, 2) = ef2_i(1, 2) - tij_ab(1, 2)*ch_j
2889 ef2_i(2, 2) = ef2_i(2, 2) - tij_ab(2, 2)*ch_j
2890 ef2_i(3, 2) = ef2_i(3, 2) - tij_ab(3, 2)*ch_j
2891 ef2_i(1, 3) = ef2_i(1, 3) - tij_ab(1, 3)*ch_j
2892 ef2_i(2, 3) = ef2_i(2, 3) - tij_ab(2, 3)*ch_j
2893 ef2_i(3, 3) = ef2_i(3, 3) - tij_ab(3, 3)*ch_j
2895 ef2_j(1, 1) = ef2_j(1, 1) - tij_ab(1, 1)*ch_i
2896 ef2_j(2, 1) = ef2_j(2, 1) - tij_ab(2, 1)*ch_i
2897 ef2_j(3, 1) = ef2_j(3, 1) - tij_ab(3, 1)*ch_i
2898 ef2_j(1, 2) = ef2_j(1, 2) - tij_ab(1, 2)*ch_i
2899 ef2_j(2, 2) = ef2_j(2, 2) - tij_ab(2, 2)*ch_i
2900 ef2_j(3, 2) = ef2_j(3, 2) - tij_ab(3, 2)*ch_i
2901 ef2_j(1, 3) = ef2_j(1, 3) - tij_ab(1, 3)*ch_i
2902 ef2_j(2, 3) = ef2_j(2, 3) - tij_ab(2, 3)*ch_i
2903 ef2_j(3, 3) = ef2_j(3, 3) - tij_ab(3, 3)*ch_i
2907 IF (task(2, 2))
THEN
2909 tmp = -(dp_i(1)*(tij_ab(1, 1)*dp_j(1) + &
2910 tij_ab(2, 1)*dp_j(2) + &
2911 tij_ab(3, 1)*dp_j(3)) + &
2912 dp_i(2)*(tij_ab(1, 2)*dp_j(1) + &
2913 tij_ab(2, 2)*dp_j(2) + &
2914 tij_ab(3, 2)*dp_j(3)) + &
2915 dp_i(3)*(tij_ab(1, 3)*dp_j(1) + &
2916 tij_ab(2, 3)*dp_j(2) + &
2917 tij_ab(3, 3)*dp_j(3)))
2920 IF (do_forces .OR. do_stress)
THEN
2922 fr(k) = fr(k) + dp_i(1)*(tij_abc(1, 1, k)*dp_j(1) + &
2923 tij_abc(2, 1, k)*dp_j(2) + &
2924 tij_abc(3, 1, k)*dp_j(3)) &
2925 + dp_i(2)*(tij_abc(1, 2, k)*dp_j(1) + &
2926 tij_abc(2, 2, k)*dp_j(2) + &
2927 tij_abc(3, 2, k)*dp_j(3)) &
2928 + dp_i(3)*(tij_abc(1, 3, k)*dp_j(1) + &
2929 tij_abc(2, 3, k)*dp_j(2) + &
2930 tij_abc(3, 3, k)*dp_j(3))
2936 IF (do_efield0)
THEN
2937 ef0_i = ef0_i - (tij_a(1)*dp_j(1) + &
2938 tij_a(2)*dp_j(2) + &
2941 ef0_j = ef0_j + (tij_a(1)*dp_i(1) + &
2942 tij_a(2)*dp_i(2) + &
2946 IF (do_efield1)
THEN
2947 ef1_i(1) = ef1_i(1) + (tij_ab(1, 1)*dp_j(1) + &
2948 tij_ab(2, 1)*dp_j(2) + &
2949 tij_ab(3, 1)*dp_j(3))
2950 ef1_i(2) = ef1_i(2) + (tij_ab(1, 2)*dp_j(1) + &
2951 tij_ab(2, 2)*dp_j(2) + &
2952 tij_ab(3, 2)*dp_j(3))
2953 ef1_i(3) = ef1_i(3) + (tij_ab(1, 3)*dp_j(1) + &
2954 tij_ab(2, 3)*dp_j(2) + &
2955 tij_ab(3, 3)*dp_j(3))
2957 ef1_j(1) = ef1_j(1) + (tij_ab(1, 1)*dp_i(1) + &
2958 tij_ab(2, 1)*dp_i(2) + &
2959 tij_ab(3, 1)*dp_i(3))
2960 ef1_j(2) = ef1_j(2) + (tij_ab(1, 2)*dp_i(1) + &
2961 tij_ab(2, 2)*dp_i(2) + &
2962 tij_ab(3, 2)*dp_i(3))
2963 ef1_j(3) = ef1_j(3) + (tij_ab(1, 3)*dp_i(1) + &
2964 tij_ab(2, 3)*dp_i(2) + &
2965 tij_ab(3, 3)*dp_i(3))
2968 IF (do_efield2)
THEN
2969 ef2_i(1, 1) = ef2_i(1, 1) + (tij_abc(1, 1, 1)*dp_j(1) + &
2970 tij_abc(2, 1, 1)*dp_j(2) + &
2971 tij_abc(3, 1, 1)*dp_j(3))
2972 ef2_i(1, 2) = ef2_i(1, 2) + (tij_abc(1, 1, 2)*dp_j(1) + &
2973 tij_abc(2, 1, 2)*dp_j(2) + &
2974 tij_abc(3, 1, 2)*dp_j(3))
2975 ef2_i(1, 3) = ef2_i(1, 3) + (tij_abc(1, 1, 3)*dp_j(1) + &
2976 tij_abc(2, 1, 3)*dp_j(2) + &
2977 tij_abc(3, 1, 3)*dp_j(3))
2978 ef2_i(2, 1) = ef2_i(2, 1) + (tij_abc(1, 2, 1)*dp_j(1) + &
2979 tij_abc(2, 2, 1)*dp_j(2) + &
2980 tij_abc(3, 2, 1)*dp_j(3))
2981 ef2_i(2, 2) = ef2_i(2, 2) + (tij_abc(1, 2, 2)*dp_j(1) + &
2982 tij_abc(2, 2, 2)*dp_j(2) + &
2983 tij_abc(3, 2, 2)*dp_j(3))
2984 ef2_i(2, 3) = ef2_i(2, 3) + (tij_abc(1, 2, 3)*dp_j(1) + &
2985 tij_abc(2, 2, 3)*dp_j(2) + &
2986 tij_abc(3, 2, 3)*dp_j(3))
2987 ef2_i(3, 1) = ef2_i(3, 1) + (tij_abc(1, 3, 1)*dp_j(1) + &
2988 tij_abc(2, 3, 1)*dp_j(2) + &
2989 tij_abc(3, 3, 1)*dp_j(3))
2990 ef2_i(3, 2) = ef2_i(3, 2) + (tij_abc(1, 3, 2)*dp_j(1) + &
2991 tij_abc(2, 3, 2)*dp_j(2) + &
2992 tij_abc(3, 3, 2)*dp_j(3))
2993 ef2_i(3, 3) = ef2_i(3, 3) + (tij_abc(1, 3, 3)*dp_j(1) + &
2994 tij_abc(2, 3, 3)*dp_j(2) + &
2995 tij_abc(3, 3, 3)*dp_j(3))
2997 ef2_j(1, 1) = ef2_j(1, 1) - (tij_abc(1, 1, 1)*dp_i(1) + &
2998 tij_abc(2, 1, 1)*dp_i(2) + &
2999 tij_abc(3, 1, 1)*dp_i(3))
3000 ef2_j(1, 2) = ef2_j(1, 2) - (tij_abc(1, 1, 2)*dp_i(1) + &
3001 tij_abc(2, 1, 2)*dp_i(2) + &
3002 tij_abc(3, 1, 2)*dp_i(3))
3003 ef2_j(1, 3) = ef2_j(1, 3) - (tij_abc(1, 1, 3)*dp_i(1) + &
3004 tij_abc(2, 1, 3)*dp_i(2) + &
3005 tij_abc(3, 1, 3)*dp_i(3))
3006 ef2_j(2, 1) = ef2_j(2, 1) - (tij_abc(1, 2, 1)*dp_i(1) + &
3007 tij_abc(2, 2, 1)*dp_i(2) + &
3008 tij_abc(3, 2, 1)*dp_i(3))
3009 ef2_j(2, 2) = ef2_j(2, 2) - (tij_abc(1, 2, 2)*dp_i(1) + &
3010 tij_abc(2, 2, 2)*dp_i(2) + &
3011 tij_abc(3, 2, 2)*dp_i(3))
3012 ef2_j(2, 3) = ef2_j(2, 3) - (tij_abc(1, 2, 3)*dp_i(1) + &
3013 tij_abc(2, 2, 3)*dp_i(2) + &
3014 tij_abc(3, 2, 3)*dp_i(3))
3015 ef2_j(3, 1) = ef2_j(3, 1) - (tij_abc(1, 3, 1)*dp_i(1) + &
3016 tij_abc(2, 3, 1)*dp_i(2) + &
3017 tij_abc(3, 3, 1)*dp_i(3))
3018 ef2_j(3, 2) = ef2_j(3, 2) - (tij_abc(1, 3, 2)*dp_i(1) + &
3019 tij_abc(2, 3, 2)*dp_i(2) + &
3020 tij_abc(3, 3, 2)*dp_i(3))
3021 ef2_j(3, 3) = ef2_j(3, 3) - (tij_abc(1, 3, 3)*dp_i(1) + &
3022 tij_abc(2, 3, 3)*dp_i(2) + &
3023 tij_abc(3, 3, 3)*dp_i(3))
3027 IF (task(2, 1))
THEN
3029 tmp = ch_j*(tij_a(1)*dp_i(1) + &
3030 tij_a(2)*dp_i(2) + &
3032 - ch_i*(tij_a(1)*dp_j(1) + &
3033 tij_a(2)*dp_j(2) + &
3037 IF (do_forces .OR. do_stress)
THEN
3039 fr(k) = fr(k) - ch_j*(tij_ab(1, k)*dp_i(1) + &
3040 tij_ab(2, k)*dp_i(2) + &
3041 tij_ab(3, k)*dp_i(3)) &
3042 + ch_i*(tij_ab(1, k)*dp_j(1) + &
3043 tij_ab(2, k)*dp_j(2) + &
3044 tij_ab(3, k)*dp_j(3))
3048 IF (task(3, 3))
THEN
3051 tmp11 = qp_i(1, 1)*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
3052 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
3053 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
3054 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
3055 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
3056 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
3057 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
3058 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
3059 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
3060 tmp21 = qp_i(2, 1)*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
3061 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
3062 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
3063 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
3064 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
3065 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
3066 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
3067 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
3068 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
3069 tmp31 = qp_i(3, 1)*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
3070 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
3071 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
3072 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
3073 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
3074 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
3075 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
3076 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
3077 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
3078 tmp22 = qp_i(2, 2)*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
3079 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
3080 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
3081 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
3082 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
3083 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
3084 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
3085 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
3086 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
3087 tmp32 = qp_i(3, 2)*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
3088 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
3089 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
3090 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
3091 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
3092 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
3093 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
3094 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
3095 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
3096 tmp33 = qp_i(3, 3)*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
3097 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
3098 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
3099 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
3100 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
3101 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
3102 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
3103 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
3104 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
3108 tmp = tmp11 + tmp12 + tmp13 + &
3109 tmp21 + tmp22 + tmp23 + &
3110 tmp31 + tmp32 + tmp33
3112 eloc = eloc +
fac*tmp
3114 IF (do_forces .OR. do_stress)
THEN
3116 tmp11 = qp_i(1, 1)*(tij_abcde(1, 1, 1, 1, k)*qp_j(1, 1) + &
3117 tij_abcde(2, 1, 1, 1, k)*qp_j(2, 1) + &
3118 tij_abcde(3, 1, 1, 1, k)*qp_j(3, 1) + &
3119 tij_abcde(1, 2, 1, 1, k)*qp_j(1, 2) + &
3120 tij_abcde(2, 2, 1, 1, k)*qp_j(2, 2) + &
3121 tij_abcde(3, 2, 1, 1, k)*qp_j(3, 2) + &
3122 tij_abcde(1, 3, 1, 1, k)*qp_j(1, 3) + &
3123 tij_abcde(2, 3, 1, 1, k)*qp_j(2, 3) + &
3124 tij_abcde(3, 3, 1, 1, k)*qp_j(3, 3))
3125 tmp21 = qp_i(2, 1)*(tij_abcde(1, 1, 2, 1, k)*qp_j(1, 1) + &
3126 tij_abcde(2, 1, 2, 1, k)*qp_j(2, 1) + &
3127 tij_abcde(3, 1, 2, 1, k)*qp_j(3, 1) + &
3128 tij_abcde(1, 2, 2, 1, k)*qp_j(1, 2) + &
3129 tij_abcde(2, 2, 2, 1, k)*qp_j(2, 2) + &
3130 tij_abcde(3, 2, 2, 1, k)*qp_j(3, 2) + &
3131 tij_abcde(1, 3, 2, 1, k)*qp_j(1, 3) + &
3132 tij_abcde(2, 3, 2, 1, k)*qp_j(2, 3) + &
3133 tij_abcde(3, 3, 2, 1, k)*qp_j(3, 3))
3134 tmp31 = qp_i(3, 1)*(tij_abcde(1, 1, 3, 1, k)*qp_j(1, 1) + &
3135 tij_abcde(2, 1, 3, 1, k)*qp_j(2, 1) + &
3136 tij_abcde(3, 1, 3, 1, k)*qp_j(3, 1) + &
3137 tij_abcde(1, 2, 3, 1, k)*qp_j(1, 2) + &
3138 tij_abcde(2, 2, 3, 1, k)*qp_j(2, 2) + &
3139 tij_abcde(3, 2, 3, 1, k)*qp_j(3, 2) + &
3140 tij_abcde(1, 3, 3, 1, k)*qp_j(1, 3) + &
3141 tij_abcde(2, 3, 3, 1, k)*qp_j(2, 3) + &
3142 tij_abcde(3, 3, 3, 1, k)*qp_j(3, 3))
3143 tmp22 = qp_i(2, 2)*(tij_abcde(1, 1, 2, 2, k)*qp_j(1, 1) + &
3144 tij_abcde(2, 1, 2, 2, k)*qp_j(2, 1) + &
3145 tij_abcde(3, 1, 2, 2, k)*qp_j(3, 1) + &
3146 tij_abcde(1, 2, 2, 2, k)*qp_j(1, 2) + &
3147 tij_abcde(2, 2, 2, 2, k)*qp_j(2, 2) + &
3148 tij_abcde(3, 2, 2, 2, k)*qp_j(3, 2) + &
3149 tij_abcde(1, 3, 2, 2, k)*qp_j(1, 3) + &
3150 tij_abcde(2, 3, 2, 2, k)*qp_j(2, 3) + &
3151 tij_abcde(3, 3, 2, 2, k)*qp_j(3, 3))
3152 tmp32 = qp_i(3, 2)*(tij_abcde(1, 1, 3, 2, k)*qp_j(1, 1) + &
3153 tij_abcde(2, 1, 3, 2, k)*qp_j(2, 1) + &
3154 tij_abcde(3, 1, 3, 2, k)*qp_j(3, 1) + &
3155 tij_abcde(1, 2, 3, 2, k)*qp_j(1, 2) + &
3156 tij_abcde(2, 2, 3, 2, k)*qp_j(2, 2) + &
3157 tij_abcde(3, 2, 3, 2, k)*qp_j(3, 2) + &
3158 tij_abcde(1, 3, 3, 2, k)*qp_j(1, 3) + &
3159 tij_abcde(2, 3, 3, 2, k)*qp_j(2, 3) + &
3160 tij_abcde(3, 3, 3, 2, k)*qp_j(3, 3))
3161 tmp33 = qp_i(3, 3)*(tij_abcde(1, 1, 3, 3, k)*qp_j(1, 1) + &
3162 tij_abcde(2, 1, 3, 3, k)*qp_j(2, 1) + &
3163 tij_abcde(3, 1, 3, 3, k)*qp_j(3, 1) + &
3164 tij_abcde(1, 2, 3, 3, k)*qp_j(1, 2) + &
3165 tij_abcde(2, 2, 3, 3, k)*qp_j(2, 2) + &
3166 tij_abcde(3, 2, 3, 3, k)*qp_j(3, 2) + &
3167 tij_abcde(1, 3, 3, 3, k)*qp_j(1, 3) + &
3168 tij_abcde(2, 3, 3, 3, k)*qp_j(2, 3) + &
3169 tij_abcde(3, 3, 3, 3, k)*qp_j(3, 3))
3173 fr(k) = fr(k) -
fac*(tmp11 + tmp12 + tmp13 + &
3174 tmp21 + tmp22 + tmp23 + &
3175 tmp31 + tmp32 + tmp33)
3182 IF (do_efield0)
THEN
3183 ef0_i = ef0_i +
fac*(tij_ab(1, 1)*qp_j(1, 1) + &
3184 tij_ab(2, 1)*qp_j(2, 1) + &
3185 tij_ab(3, 1)*qp_j(3, 1) + &
3186 tij_ab(1, 2)*qp_j(1, 2) + &
3187 tij_ab(2, 2)*qp_j(2, 2) + &
3188 tij_ab(3, 2)*qp_j(3, 2) + &
3189 tij_ab(1, 3)*qp_j(1, 3) + &
3190 tij_ab(2, 3)*qp_j(2, 3) + &
3191 tij_ab(3, 3)*qp_j(3, 3))
3193 ef0_j = ef0_j +
fac*(tij_ab(1, 1)*qp_i(1, 1) + &
3194 tij_ab(2, 1)*qp_i(2, 1) + &
3195 tij_ab(3, 1)*qp_i(3, 1) + &
3196 tij_ab(1, 2)*qp_i(1, 2) + &
3197 tij_ab(2, 2)*qp_i(2, 2) + &
3198 tij_ab(3, 2)*qp_i(3, 2) + &
3199 tij_ab(1, 3)*qp_i(1, 3) + &
3200 tij_ab(2, 3)*qp_i(2, 3) + &
3201 tij_ab(3, 3)*qp_i(3, 3))
3204 IF (do_efield1)
THEN
3205 ef1_i(1) = ef1_i(1) -
fac*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
3206 tij_abc(2, 1, 1)*qp_j(2, 1) + &
3207 tij_abc(3, 1, 1)*qp_j(3, 1) + &
3208 tij_abc(1, 2, 1)*qp_j(1, 2) + &
3209 tij_abc(2, 2, 1)*qp_j(2, 2) + &
3210 tij_abc(3, 2, 1)*qp_j(3, 2) + &
3211 tij_abc(1, 3, 1)*qp_j(1, 3) + &
3212 tij_abc(2, 3, 1)*qp_j(2, 3) + &
3213 tij_abc(3, 3, 1)*qp_j(3, 3))
3214 ef1_i(2) = ef1_i(2) -
fac*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
3215 tij_abc(2, 1, 2)*qp_j(2, 1) + &
3216 tij_abc(3, 1, 2)*qp_j(3, 1) + &
3217 tij_abc(1, 2, 2)*qp_j(1, 2) + &
3218 tij_abc(2, 2, 2)*qp_j(2, 2) + &
3219 tij_abc(3, 2, 2)*qp_j(3, 2) + &
3220 tij_abc(1, 3, 2)*qp_j(1, 3) + &
3221 tij_abc(2, 3, 2)*qp_j(2, 3) + &
3222 tij_abc(3, 3, 2)*qp_j(3, 3))
3223 ef1_i(3) = ef1_i(3) -
fac*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
3224 tij_abc(2, 1, 3)*qp_j(2, 1) + &
3225 tij_abc(3, 1, 3)*qp_j(3, 1) + &
3226 tij_abc(1, 2, 3)*qp_j(1, 2) + &
3227 tij_abc(2, 2, 3)*qp_j(2, 2) + &
3228 tij_abc(3, 2, 3)*qp_j(3, 2) + &
3229 tij_abc(1, 3, 3)*qp_j(1, 3) + &
3230 tij_abc(2, 3, 3)*qp_j(2, 3) + &
3231 tij_abc(3, 3, 3)*qp_j(3, 3))
3233 ef1_j(1) = ef1_j(1) +
fac*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
3234 tij_abc(2, 1, 1)*qp_i(2, 1) + &
3235 tij_abc(3, 1, 1)*qp_i(3, 1) + &
3236 tij_abc(1, 2, 1)*qp_i(1, 2) + &
3237 tij_abc(2, 2, 1)*qp_i(2, 2) + &
3238 tij_abc(3, 2, 1)*qp_i(3, 2) + &
3239 tij_abc(1, 3, 1)*qp_i(1, 3) + &
3240 tij_abc(2, 3, 1)*qp_i(2, 3) + &
3241 tij_abc(3, 3, 1)*qp_i(3, 3))
3242 ef1_j(2) = ef1_j(2) +
fac*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
3243 tij_abc(2, 1, 2)*qp_i(2, 1) + &
3244 tij_abc(3, 1, 2)*qp_i(3, 1) + &
3245 tij_abc(1, 2, 2)*qp_i(1, 2) + &
3246 tij_abc(2, 2, 2)*qp_i(2, 2) + &
3247 tij_abc(3, 2, 2)*qp_i(3, 2) + &
3248 tij_abc(1, 3, 2)*qp_i(1, 3) + &
3249 tij_abc(2, 3, 2)*qp_i(2, 3) + &
3250 tij_abc(3, 3, 2)*qp_i(3, 3))
3251 ef1_j(3) = ef1_j(3) +
fac*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
3252 tij_abc(2, 1, 3)*qp_i(2, 1) + &
3253 tij_abc(3, 1, 3)*qp_i(3, 1) + &
3254 tij_abc(1, 2, 3)*qp_i(1, 2) + &
3255 tij_abc(2, 2, 3)*qp_i(2, 2) + &
3256 tij_abc(3, 2, 3)*qp_i(3, 2) + &
3257 tij_abc(1, 3, 3)*qp_i(1, 3) + &
3258 tij_abc(2, 3, 3)*qp_i(2, 3) + &
3259 tij_abc(3, 3, 3)*qp_i(3, 3))
3262 IF (do_efield2)
THEN
3263 tmp11 =
fac*(tij_abcd(1, 1, 1, 1)*qp_j(1, 1) + &
3264 tij_abcd(2, 1, 1, 1)*qp_j(2, 1) + &
3265 tij_abcd(3, 1, 1, 1)*qp_j(3, 1) + &
3266 tij_abcd(1, 2, 1, 1)*qp_j(1, 2) + &
3267 tij_abcd(2, 2, 1, 1)*qp_j(2, 2) + &
3268 tij_abcd(3, 2, 1, 1)*qp_j(3, 2) + &
3269 tij_abcd(1, 3, 1, 1)*qp_j(1, 3) + &
3270 tij_abcd(2, 3, 1, 1)*qp_j(2, 3) + &
3271 tij_abcd(3, 3, 1, 1)*qp_j(3, 3))
3272 tmp12 =
fac*(tij_abcd(1, 1, 1, 2)*qp_j(1, 1) + &
3273 tij_abcd(2, 1, 1, 2)*qp_j(2, 1) + &
3274 tij_abcd(3, 1, 1, 2)*qp_j(3, 1) + &
3275 tij_abcd(1, 2, 1, 2)*qp_j(1, 2) + &
3276 tij_abcd(2, 2, 1, 2)*qp_j(2, 2) + &
3277 tij_abcd(3, 2, 1, 2)*qp_j(3, 2) + &
3278 tij_abcd(1, 3, 1, 2)*qp_j(1, 3) + &
3279 tij_abcd(2, 3, 1, 2)*qp_j(2, 3) + &
3280 tij_abcd(3, 3, 1, 2)*qp_j(3, 3))
3281 tmp13 =
fac*(tij_abcd(1, 1, 1, 3)*qp_j(1, 1) + &
3282 tij_abcd(2, 1, 1, 3)*qp_j(2, 1) + &
3283 tij_abcd(3, 1, 1, 3)*qp_j(3, 1) + &
3284 tij_abcd(1, 2, 1, 3)*qp_j(1, 2) + &
3285 tij_abcd(2, 2, 1, 3)*qp_j(2, 2) + &
3286 tij_abcd(3, 2, 1, 3)*qp_j(3, 2) + &
3287 tij_abcd(1, 3, 1, 3)*qp_j(1, 3) + &
3288 tij_abcd(2, 3, 1, 3)*qp_j(2, 3) + &
3289 tij_abcd(3, 3, 1, 3)*qp_j(3, 3))
3290 tmp22 =
fac*(tij_abcd(1, 1, 2, 2)*qp_j(1, 1) + &
3291 tij_abcd(2, 1, 2, 2)*qp_j(2, 1) + &
3292 tij_abcd(3, 1, 2, 2)*qp_j(3, 1) + &
3293 tij_abcd(1, 2, 2, 2)*qp_j(1, 2) + &
3294 tij_abcd(2, 2, 2, 2)*qp_j(2, 2) + &
3295 tij_abcd(3, 2, 2, 2)*qp_j(3, 2) + &
3296 tij_abcd(1, 3, 2, 2)*qp_j(1, 3) + &
3297 tij_abcd(2, 3, 2, 2)*qp_j(2, 3) + &
3298 tij_abcd(3, 3, 2, 2)*qp_j(3, 3))
3299 tmp23 =
fac*(tij_abcd(1, 1, 2, 3)*qp_j(1, 1) + &
3300 tij_abcd(2, 1, 2, 3)*qp_j(2, 1) + &
3301 tij_abcd(3, 1, 2, 3)*qp_j(3, 1) + &
3302 tij_abcd(1, 2, 2, 3)*qp_j(1, 2) + &
3303 tij_abcd(2, 2, 2, 3)*qp_j(2, 2) + &
3304 tij_abcd(3, 2, 2, 3)*qp_j(3, 2) + &
3305 tij_abcd(1, 3, 2, 3)*qp_j(1, 3) + &
3306 tij_abcd(2, 3, 2, 3)*qp_j(2, 3) + &
3307 tij_abcd(3, 3, 2, 3)*qp_j(3, 3))
3308 tmp33 =
fac*(tij_abcd(1, 1, 3, 3)*qp_j(1, 1) + &
3309 tij_abcd(2, 1, 3, 3)*qp_j(2, 1) + &
3310 tij_abcd(3, 1, 3, 3)*qp_j(3, 1) + &
3311 tij_abcd(1, 2, 3, 3)*qp_j(1, 2) + &
3312 tij_abcd(2, 2, 3, 3)*qp_j(2, 2) + &
3313 tij_abcd(3, 2, 3, 3)*qp_j(3, 2) + &
3314 tij_abcd(1, 3, 3, 3)*qp_j(1, 3) + &
3315 tij_abcd(2, 3, 3, 3)*qp_j(2, 3) + &
3316 tij_abcd(3, 3, 3, 3)*qp_j(3, 3))
3318 ef2_i(1, 1) = ef2_i(1, 1) - tmp11
3319 ef2_i(1, 2) = ef2_i(1, 2) - tmp12
3320 ef2_i(1, 3) = ef2_i(1, 3) - tmp13
3321 ef2_i(2, 1) = ef2_i(2, 1) - tmp12
3322 ef2_i(2, 2) = ef2_i(2, 2) - tmp22
3323 ef2_i(2, 3) = ef2_i(2, 3) - tmp23
3324 ef2_i(3, 1) = ef2_i(3, 1) - tmp13
3325 ef2_i(3, 2) = ef2_i(3, 2) - tmp23
3326 ef2_i(3, 3) = ef2_i(3, 3) - tmp33
3328 tmp11 =
fac*(tij_abcd(1, 1, 1, 1)*qp_i(1, 1) + &
3329 tij_abcd(2, 1, 1, 1)*qp_i(2, 1) + &
3330 tij_abcd(3, 1, 1, 1)*qp_i(3, 1) + &
3331 tij_abcd(1, 2, 1, 1)*qp_i(1, 2) + &
3332 tij_abcd(2, 2, 1, 1)*qp_i(2, 2) + &
3333 tij_abcd(3, 2, 1, 1)*qp_i(3, 2) + &
3334 tij_abcd(1, 3, 1, 1)*qp_i(1, 3) + &
3335 tij_abcd(2, 3, 1, 1)*qp_i(2, 3) + &
3336 tij_abcd(3, 3, 1, 1)*qp_i(3, 3))
3337 tmp12 =
fac*(tij_abcd(1, 1, 1, 2)*qp_i(1, 1) + &
3338 tij_abcd(2, 1, 1, 2)*qp_i(2, 1) + &
3339 tij_abcd(3, 1, 1, 2)*qp_i(3, 1) + &
3340 tij_abcd(1, 2, 1, 2)*qp_i(1, 2) + &
3341 tij_abcd(2, 2, 1, 2)*qp_i(2, 2) + &
3342 tij_abcd(3, 2, 1, 2)*qp_i(3, 2) + &
3343 tij_abcd(1, 3, 1, 2)*qp_i(1, 3) + &
3344 tij_abcd(2, 3, 1, 2)*qp_i(2, 3) + &
3345 tij_abcd(3, 3, 1, 2)*qp_i(3, 3))
3346 tmp13 =
fac*(tij_abcd(1, 1, 1, 3)*qp_i(1, 1) + &
3347 tij_abcd(2, 1, 1, 3)*qp_i(2, 1) + &
3348 tij_abcd(3, 1, 1, 3)*qp_i(3, 1) + &
3349 tij_abcd(1, 2, 1, 3)*qp_i(1, 2) + &
3350 tij_abcd(2, 2, 1, 3)*qp_i(2, 2) + &
3351 tij_abcd(3, 2, 1, 3)*qp_i(3, 2) + &
3352 tij_abcd(1, 3, 1, 3)*qp_i(1, 3) + &
3353 tij_abcd(2, 3, 1, 3)*qp_i(2, 3) + &
3354 tij_abcd(3, 3, 1, 3)*qp_i(3, 3))
3355 tmp22 =
fac*(tij_abcd(1, 1, 2, 2)*qp_i(1, 1) + &
3356 tij_abcd(2, 1, 2, 2)*qp_i(2, 1) + &
3357 tij_abcd(3, 1, 2, 2)*qp_i(3, 1) + &
3358 tij_abcd(1, 2, 2, 2)*qp_i(1, 2) + &
3359 tij_abcd(2, 2, 2, 2)*qp_i(2, 2) + &
3360 tij_abcd(3, 2, 2, 2)*qp_i(3, 2) + &
3361 tij_abcd(1, 3, 2, 2)*qp_i(1, 3) + &
3362 tij_abcd(2, 3, 2, 2)*qp_i(2, 3) + &
3363 tij_abcd(3, 3, 2, 2)*qp_i(3, 3))
3364 tmp23 =
fac*(tij_abcd(1, 1, 2, 3)*qp_i(1, 1) + &
3365 tij_abcd(2, 1, 2, 3)*qp_i(2, 1) + &
3366 tij_abcd(3, 1, 2, 3)*qp_i(3, 1) + &
3367 tij_abcd(1, 2, 2, 3)*qp_i(1, 2) + &
3368 tij_abcd(2, 2, 2, 3)*qp_i(2, 2) + &
3369 tij_abcd(3, 2, 2, 3)*qp_i(3, 2) + &
3370 tij_abcd(1, 3, 2, 3)*qp_i(1, 3) + &
3371 tij_abcd(2, 3, 2, 3)*qp_i(2, 3) + &
3372 tij_abcd(3, 3, 2, 3)*qp_i(3, 3))
3373 tmp33 =
fac*(tij_abcd(1, 1, 3, 3)*qp_i(1, 1) + &
3374 tij_abcd(2, 1, 3, 3)*qp_i(2, 1) + &
3375 tij_abcd(3, 1, 3, 3)*qp_i(3, 1) + &
3376 tij_abcd(1, 2, 3, 3)*qp_i(1, 2) + &
3377 tij_abcd(2, 2, 3, 3)*qp_i(2, 2) + &
3378 tij_abcd(3, 2, 3, 3)*qp_i(3, 2) + &
3379 tij_abcd(1, 3, 3, 3)*qp_i(1, 3) + &
3380 tij_abcd(2, 3, 3, 3)*qp_i(2, 3) + &
3381 tij_abcd(3, 3, 3, 3)*qp_i(3, 3))
3383 ef2_j(1, 1) = ef2_j(1, 1) - tmp11
3384 ef2_j(1, 2) = ef2_j(1, 2) - tmp12
3385 ef2_j(1, 3) = ef2_j(1, 3) - tmp13
3386 ef2_j(2, 1) = ef2_j(2, 1) - tmp12
3387 ef2_j(2, 2) = ef2_j(2, 2) - tmp22
3388 ef2_j(2, 3) = ef2_j(2, 3) - tmp23
3389 ef2_j(3, 1) = ef2_j(3, 1) - tmp13
3390 ef2_j(3, 2) = ef2_j(3, 2) - tmp23
3391 ef2_j(3, 3) = ef2_j(3, 3) - tmp33
3395 IF (task(3, 2))
THEN
3399 tmp_ij = dp_i(1)*(tij_abc(1, 1, 1)*qp_j(1, 1) + &
3400 tij_abc(2, 1, 1)*qp_j(2, 1) + &
3401 tij_abc(3, 1, 1)*qp_j(3, 1) + &
3402 tij_abc(1, 2, 1)*qp_j(1, 2) + &
3403 tij_abc(2, 2, 1)*qp_j(2, 2) + &
3404 tij_abc(3, 2, 1)*qp_j(3, 2) + &
3405 tij_abc(1, 3, 1)*qp_j(1, 3) + &
3406 tij_abc(2, 3, 1)*qp_j(2, 3) + &
3407 tij_abc(3, 3, 1)*qp_j(3, 3)) + &
3408 dp_i(2)*(tij_abc(1, 1, 2)*qp_j(1, 1) + &
3409 tij_abc(2, 1, 2)*qp_j(2, 1) + &
3410 tij_abc(3, 1, 2)*qp_j(3, 1) + &
3411 tij_abc(1, 2, 2)*qp_j(1, 2) + &
3412 tij_abc(2, 2, 2)*qp_j(2, 2) + &
3413 tij_abc(3, 2, 2)*qp_j(3, 2) + &
3414 tij_abc(1, 3, 2)*qp_j(1, 3) + &
3415 tij_abc(2, 3, 2)*qp_j(2, 3) + &
3416 tij_abc(3, 3, 2)*qp_j(3, 3)) + &
3417 dp_i(3)*(tij_abc(1, 1, 3)*qp_j(1, 1) + &
3418 tij_abc(2, 1, 3)*qp_j(2, 1) + &
3419 tij_abc(3, 1, 3)*qp_j(3, 1) + &
3420 tij_abc(1, 2, 3)*qp_j(1, 2) + &
3421 tij_abc(2, 2, 3)*qp_j(2, 2) + &
3422 tij_abc(3, 2, 3)*qp_j(3, 2) + &
3423 tij_abc(1, 3, 3)*qp_j(1, 3) + &
3424 tij_abc(2, 3, 3)*qp_j(2, 3) + &
3425 tij_abc(3, 3, 3)*qp_j(3, 3))
3428 tmp_ji = dp_j(1)*(tij_abc(1, 1, 1)*qp_i(1, 1) + &
3429 tij_abc(2, 1, 1)*qp_i(2, 1) + &
3430 tij_abc(3, 1, 1)*qp_i(3, 1) + &
3431 tij_abc(1, 2, 1)*qp_i(1, 2) + &
3432 tij_abc(2, 2, 1)*qp_i(2, 2) + &
3433 tij_abc(3, 2, 1)*qp_i(3, 2) + &
3434 tij_abc(1, 3, 1)*qp_i(1, 3) + &
3435 tij_abc(2, 3, 1)*qp_i(2, 3) + &
3436 tij_abc(3, 3, 1)*qp_i(3, 3)) + &
3437 dp_j(2)*(tij_abc(1, 1, 2)*qp_i(1, 1) + &
3438 tij_abc(2, 1, 2)*qp_i(2, 1) + &
3439 tij_abc(3, 1, 2)*qp_i(3, 1) + &
3440 tij_abc(1, 2, 2)*qp_i(1, 2) + &
3441 tij_abc(2, 2, 2)*qp_i(2, 2) + &
3442 tij_abc(3, 2, 2)*qp_i(3, 2) + &
3443 tij_abc(1, 3, 2)*qp_i(1, 3) + &
3444 tij_abc(2, 3, 2)*qp_i(2, 3) + &
3445 tij_abc(3, 3, 2)*qp_i(3, 3)) + &
3446 dp_j(3)*(tij_abc(1, 1, 3)*qp_i(1, 1) + &
3447 tij_abc(2, 1, 3)*qp_i(2, 1) + &
3448 tij_abc(3, 1, 3)*qp_i(3, 1) + &
3449 tij_abc(1, 2, 3)*qp_i(1, 2) + &
3450 tij_abc(2, 2, 3)*qp_i(2, 2) + &
3451 tij_abc(3, 2, 3)*qp_i(3, 2) + &
3452 tij_abc(1, 3, 3)*qp_i(1, 3) + &
3453 tij_abc(2, 3, 3)*qp_i(2, 3) + &
3454 tij_abc(3, 3, 3)*qp_i(3, 3))
3456 tmp =
fac*(tmp_ij - tmp_ji)
3458 IF (do_forces .OR. do_stress)
THEN
3461 tmp_ij = dp_i(1)*(tij_abcd(1, 1, 1, k)*qp_j(1, 1) + &
3462 tij_abcd(2, 1, 1, k)*qp_j(2, 1) + &
3463 tij_abcd(3, 1, 1, k)*qp_j(3, 1) + &
3464 tij_abcd(1, 2, 1, k)*qp_j(1, 2) + &
3465 tij_abcd(2, 2, 1, k)*qp_j(2, 2) + &
3466 tij_abcd(3, 2, 1, k)*qp_j(3, 2) + &
3467 tij_abcd(1, 3, 1, k)*qp_j(1, 3) + &
3468 tij_abcd(2, 3, 1, k)*qp_j(2, 3) + &
3469 tij_abcd(3, 3, 1, k)*qp_j(3, 3)) + &
3470 dp_i(2)*(tij_abcd(1, 1, 2, k)*qp_j(1, 1) + &
3471 tij_abcd(2, 1, 2, k)*qp_j(2, 1) + &
3472 tij_abcd(3, 1, 2, k)*qp_j(3, 1) + &
3473 tij_abcd(1, 2, 2, k)*qp_j(1, 2) + &
3474 tij_abcd(2, 2, 2, k)*qp_j(2, 2) + &
3475 tij_abcd(3, 2, 2, k)*qp_j(3, 2) + &
3476 tij_abcd(1, 3, 2, k)*qp_j(1, 3) + &
3477 tij_abcd(2, 3, 2, k)*qp_j(2, 3) + &
3478 tij_abcd(3, 3, 2, k)*qp_j(3, 3)) + &
3479 dp_i(3)*(tij_abcd(1, 1, 3, k)*qp_j(1, 1) + &
3480 tij_abcd(2, 1, 3, k)*qp_j(2, 1) + &
3481 tij_abcd(3, 1, 3, k)*qp_j(3, 1) + &
3482 tij_abcd(1, 2, 3, k)*qp_j(1, 2) + &
3483 tij_abcd(2, 2, 3, k)*qp_j(2, 2) + &
3484 tij_abcd(3, 2, 3, k)*qp_j(3, 2) + &
3485 tij_abcd(1, 3, 3, k)*qp_j(1, 3) + &
3486 tij_abcd(2, 3, 3, k)*qp_j(2, 3) + &
3487 tij_abcd(3, 3, 3, k)*qp_j(3, 3))
3490 tmp_ji = dp_j(1)*(tij_abcd(1, 1, 1, k)*qp_i(1, 1) + &
3491 tij_abcd(2, 1, 1, k)*qp_i(2, 1) + &
3492 tij_abcd(3, 1, 1, k)*qp_i(3, 1) + &
3493 tij_abcd(1, 2, 1, k)*qp_i(1, 2) + &
3494 tij_abcd(2, 2, 1, k)*qp_i(2, 2) + &
3495 tij_abcd(3, 2, 1, k)*qp_i(3, 2) + &
3496 tij_abcd(1, 3, 1, k)*qp_i(1, 3) + &
3497 tij_abcd(2, 3, 1, k)*qp_i(2, 3) + &
3498 tij_abcd(3, 3, 1, k)*qp_i(3, 3)) + &
3499 dp_j(2)*(tij_abcd(1, 1, 2, k)*qp_i(1, 1) + &
3500 tij_abcd(2, 1, 2, k)*qp_i(2, 1) + &
3501 tij_abcd(3, 1, 2, k)*qp_i(3, 1) + &
3502 tij_abcd(1, 2, 2, k)*qp_i(1, 2) + &
3503 tij_abcd(2, 2, 2, k)*qp_i(2, 2) + &
3504 tij_abcd(3, 2, 2, k)*qp_i(3, 2) + &
3505 tij_abcd(1, 3, 2, k)*qp_i(1, 3) + &
3506 tij_abcd(2, 3, 2, k)*qp_i(2, 3) + &
3507 tij_abcd(3, 3, 2, k)*qp_i(3, 3)) + &
3508 dp_j(3)*(tij_abcd(1, 1, 3, k)*qp_i(1, 1) + &
3509 tij_abcd(2, 1, 3, k)*qp_i(2, 1) + &
3510 tij_abcd(3, 1, 3, k)*qp_i(3, 1) + &
3511 tij_abcd(1, 2, 3, k)*qp_i(1, 2) + &
3512 tij_abcd(2, 2, 3, k)*qp_i(2, 2) + &
3513 tij_abcd(3, 2, 3, k)*qp_i(3, 2) + &
3514 tij_abcd(1, 3, 3, k)*qp_i(1, 3) + &
3515 tij_abcd(2, 3, 3, k)*qp_i(2, 3) + &
3516 tij_abcd(3, 3, 3, k)*qp_i(3, 3))
3518 fr(k) = fr(k) -
fac*(tmp_ij - tmp_ji)
3522 IF (task(3, 1))
THEN
3527 tmp_ij = ch_i*(tij_ab(1, 1)*qp_j(1, 1) + &
3528 tij_ab(2, 1)*qp_j(2, 1) + &
3529 tij_ab(3, 1)*qp_j(3, 1) + &
3530 tij_ab(1, 2)*qp_j(1, 2) + &
3531 tij_ab(2, 2)*qp_j(2, 2) + &
3532 tij_ab(3, 2)*qp_j(3, 2) + &
3533 tij_ab(1, 3)*qp_j(1, 3) + &
3534 tij_ab(2, 3)*qp_j(2, 3) + &
3535 tij_ab(3, 3)*qp_j(3, 3))
3538 tmp_ji = ch_j*(tij_ab(1, 1)*qp_i(1, 1) + &
3539 tij_ab(2, 1)*qp_i(2, 1) + &
3540 tij_ab(3, 1)*qp_i(3, 1) + &
3541 tij_ab(1, 2)*qp_i(1, 2) + &
3542 tij_ab(2, 2)*qp_i(2, 2) + &
3543 tij_ab(3, 2)*qp_i(3, 2) + &
3544 tij_ab(1, 3)*qp_i(1, 3) + &
3545 tij_ab(2, 3)*qp_i(2, 3) + &
3546 tij_ab(3, 3)*qp_i(3, 3))
3548 eloc = eloc +
fac*(tmp_ij + tmp_ji)
3549 IF (do_forces .OR. do_stress)
THEN
3552 tmp_ij = ch_i*(tij_abc(1, 1, k)*qp_j(1, 1) + &
3553 tij_abc(2, 1, k)*qp_j(2, 1) + &
3554 tij_abc(3, 1, k)*qp_j(3, 1) + &
3555 tij_abc(1, 2, k)*qp_j(1, 2) + &
3556 tij_abc(2, 2, k)*qp_j(2, 2) + &
3557 tij_abc(3, 2, k)*qp_j(3, 2) + &
3558 tij_abc(1, 3, k)*qp_j(1, 3) + &
3559 tij_abc(2, 3, k)*qp_j(2, 3) + &
3560 tij_abc(3, 3, k)*qp_j(3, 3))
3563 tmp_ji = ch_j*(tij_abc(1, 1, k)*qp_i(1, 1) + &
3564 tij_abc(2, 1, k)*qp_i(2, 1) + &
3565 tij_abc(3, 1, k)*qp_i(3, 1) + &
3566 tij_abc(1, 2, k)*qp_i(1, 2) + &
3567 tij_abc(2, 2, k)*qp_i(2, 2) + &
3568 tij_abc(3, 2, k)*qp_i(3, 2) + &
3569 tij_abc(1, 3, k)*qp_i(1, 3) + &
3570 tij_abc(2, 3, k)*qp_i(2, 3) + &
3571 tij_abc(3, 3, k)*qp_i(3, 3))
3573 fr(k) = fr(k) -
fac*(tmp_ij + tmp_ji)
3577 energy = energy + eloc
3579 forces(1, atom_a) = forces(1, atom_a) - fr(1)
3580 forces(2, atom_a) = forces(2, atom_a) - fr(2)
3581 forces(3, atom_a) = forces(3, atom_a) - fr(3)
3582 forces(1, atom_b) = forces(1, atom_b) + fr(1)
3583 forces(2, atom_b) = forces(2, atom_b) + fr(2)
3584 forces(3, atom_b) = forces(3, atom_b) + fr(3)
3589 IF (do_efield0)
THEN
3590 efield0(atom_a) = efield0(atom_a) + ef0_j
3592 efield0(atom_b) = efield0(atom_b) + ef0_i
3595 IF (do_efield1)
THEN
3596 efield1(1, atom_a) = efield1(1, atom_a) + ef1_j(1)
3597 efield1(2, atom_a) = efield1(2, atom_a) + ef1_j(2)
3598 efield1(3, atom_a) = efield1(3, atom_a) + ef1_j(3)
3600 efield1(1, atom_b) = efield1(1, atom_b) + ef1_i(1)
3601 efield1(2, atom_b) = efield1(2, atom_b) + ef1_i(2)
3602 efield1(3, atom_b) = efield1(3, atom_b) + ef1_i(3)
3605 IF (do_efield2)
THEN
3606 efield2(1, atom_a) = efield2(1, atom_a) + ef2_j(1, 1)
3607 efield2(2, atom_a) = efield2(2, atom_a) + ef2_j(1, 2)
3608 efield2(3, atom_a) = efield2(3, atom_a) + ef2_j(1, 3)
3609 efield2(4, atom_a) = efield2(4, atom_a) + ef2_j(2, 1)
3610 efield2(5, atom_a) = efield2(5, atom_a) + ef2_j(2, 2)
3611 efield2(6, atom_a) = efield2(6, atom_a) + ef2_j(2, 3)
3612 efield2(7, atom_a) = efield2(7, atom_a) + ef2_j(3, 1)
3613 efield2(8, atom_a) = efield2(8, atom_a) + ef2_j(3, 2)
3614 efield2(9, atom_a) = efield2(9, atom_a) + ef2_j(3, 3)
3616 efield2(1, atom_b) = efield2(1, atom_b) + ef2_i(1, 1)
3617 efield2(2, atom_b) = efield2(2, atom_b) + ef2_i(1, 2)
3618 efield2(3, atom_b) = efield2(3, atom_b) + ef2_i(1, 3)
3619 efield2(4, atom_b) = efield2(4, atom_b) + ef2_i(2, 1)
3620 efield2(5, atom_b) = efield2(5, atom_b) + ef2_i(2, 2)
3621 efield2(6, atom_b) = efield2(6, atom_b) + ef2_i(2, 3)
3622 efield2(7, atom_b) = efield2(7, atom_b) + ef2_i(3, 1)
3623 efield2(8, atom_b) = efield2(8, atom_b) + ef2_i(3, 2)
3624 efield2(9, atom_b) = efield2(9, atom_b) + ef2_i(3, 3)
3628 ptens11 = ptens11 + rab(1)*fr(1)
3629 ptens21 = ptens21 + rab(2)*fr(1)
3630 ptens31 = ptens31 + rab(3)*fr(1)
3631 ptens12 = ptens12 + rab(1)*fr(2)
3632 ptens22 = ptens22 + rab(2)*fr(2)
3633 ptens32 = ptens32 + rab(3)*fr(2)
3634 ptens13 = ptens13 + rab(1)*fr(3)
3635 ptens23 = ptens23 + rab(2)*fr(3)
3636 ptens33 = ptens33 + rab(3)*fr(3)
3640 END DO kind_group_loop
3643 pv(1, 1) = pv(1, 1) + ptens11
3644 pv(1, 2) = pv(1, 2) + (ptens12 + ptens21)*0.5_dp
3645 pv(1, 3) = pv(1, 3) + (ptens13 + ptens31)*0.5_dp
3647 pv(2, 2) = pv(2, 2) + ptens22
3648 pv(2, 3) = pv(2, 3) + (ptens23 + ptens32)*0.5_dp
3651 pv(3, 3) = pv(3, 3) + ptens33
3654 CALL timestop(handle)
3655 END SUBROUTINE ewald_multipole_bonded
3680 SUBROUTINE ewald_multipole_lr(ewald_env, ewald_pw, cell, particle_set, &
3681 local_particles, energy, task, do_forces, do_efield, do_stress, &
3682 charges, dipoles, quadrupoles, forces, pv, efield0, efield1, efield2)
3688 REAL(kind=
dp),
INTENT(INOUT) :: energy
3689 LOGICAL,
DIMENSION(3, 3),
INTENT(IN) :: task
3690 LOGICAL,
INTENT(IN) :: do_forces, do_efield, do_stress
3691 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: charges
3692 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
3693 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
3694 POINTER :: quadrupoles
3695 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT), &
3696 OPTIONAL :: forces, pv
3697 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0
3698 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: efield1, efield2
3700 CHARACTER(len=*),
PARAMETER :: routinen =
'ewald_multipole_LR'
3702 COMPLEX(KIND=dp) :: atm_factor, atm_factor_st(3), cnjg_fac, &
3704 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: summe_ef
3705 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:, :) :: summe_st
3706 INTEGER :: gpt, handle, iparticle, iparticle_kind, iparticle_local, lp, mp, nnodes, &
3707 node, np, nparticle_kind, nparticle_local
3708 INTEGER,
DIMENSION(:, :),
POINTER :: bds
3709 LOGICAL :: do_efield0, do_efield1, do_efield2
3710 REAL(kind=
dp) :: alpha, denom, dipole_t(3), f0, factor, &
3711 four_alpha_sq, gauss, pref, q_t, tmp, &
3713 REAL(kind=
dp),
DIMENSION(3) :: tmp_v, vec
3714 REAL(kind=
dp),
DIMENSION(3, 3) :: pv_tmp
3715 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: rho0
3723 CALL timeset(routinen, handle)
3724 do_efield0 = do_efield .AND.
ASSOCIATED(efield0)
3725 do_efield1 = do_efield .AND.
ASSOCIATED(efield1)
3726 do_efield2 = do_efield .AND.
ASSOCIATED(efield2)
3730 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, dg=dg)
3731 CALL dg_get(dg, dg_rho0=dg_rho0)
3732 rho0 => dg_rho0%density%array
3733 pw_grid => pw_pool%pw_grid
3734 bds => pw_grid%bounds
3737 nparticle_kind =
SIZE(local_particles%n_el)
3739 DO iparticle_kind = 1, nparticle_kind
3740 nnodes = nnodes + local_particles%n_el(iparticle_kind)
3744 ALLOCATE (summe_ef(1:pw_grid%ngpts_cut))
3749 ALLOCATE (summe_st(3, 1:pw_grid%ngpts_cut))
3754 four_alpha_sq = 4.0_dp*alpha**2
3760 DO iparticle_kind = 1, nparticle_kind
3761 nparticle_local = local_particles%n_el(iparticle_kind)
3762 DO iparticle_local = 1, nparticle_local
3764 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
3765 vec = matmul(cell%h_inv, particle_set(iparticle)%r)
3767 exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
3770 IF (any(task(1, :)))
THEN
3771 q_t = q_t + charges(iparticle)
3773 IF (any(task(2, :)))
THEN
3774 dipole_t = dipole_t + dipoles(:, iparticle)
3776 IF (any(task(3, :)))
THEN
3777 trq_t = trq_t + quadrupoles(1, 1, iparticle) + &
3778 quadrupoles(2, 2, iparticle) + &
3779 quadrupoles(3, 3, iparticle)
3785 DO gpt = 1, pw_grid%ngpts_cut_local
3786 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
3787 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
3788 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
3796 DO iparticle_kind = 1, nparticle_kind
3797 nparticle_local = local_particles%n_el(iparticle_kind)
3798 DO iparticle_local = 1, nparticle_local
3800 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
3802 CALL get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges, &
3803 dipoles, quadrupoles)
3804 summe_tmp = exp_igr%ex(lp, node)*exp_igr%ey(mp, node)*exp_igr%ez(np, node)
3805 summe_ef(gpt) = summe_ef(gpt) + atm_factor*summe_tmp
3809 CALL get_atom_factor_stress(atm_factor_st, pw_grid, gpt, iparticle, task, &
3810 dipoles, quadrupoles)
3811 summe_st(1:3, gpt) = summe_st(1:3, gpt) + atm_factor_st(1:3)*summe_tmp
3817 CALL group%sum(dipole_t)
3818 CALL group%sum(trq_t)
3819 CALL group%sum(summe_ef)
3820 IF (do_stress)
CALL group%sum(summe_st)
3823 DO gpt = 1, pw_grid%ngpts_cut_local
3825 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
3826 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
3827 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
3833 IF (pw_grid%gsq(gpt) == 0.0_dp)
THEN
3835 energy = energy + (1.0_dp/6.0_dp)*dot_product(dipole_t, dipole_t) &
3836 - (1.0_dp/9.0_dp)*q_t*trq_t
3839 pv_tmp(1, 1) = pv_tmp(1, 1) + (1.0_dp/6.0_dp)*dot_product(dipole_t, dipole_t)
3840 pv_tmp(2, 2) = pv_tmp(2, 2) + (1.0_dp/6.0_dp)*dot_product(dipole_t, dipole_t)
3841 pv_tmp(3, 3) = pv_tmp(3, 3) + (1.0_dp/6.0_dp)*dot_product(dipole_t, dipole_t)
3844 IF (do_efield .AND. (debug_e_field_en .OR. (.NOT. debug_this_module)))
THEN
3849 DO iparticle_kind = 1, nparticle_kind
3850 nparticle_local = local_particles%n_el(iparticle_kind)
3851 DO iparticle_local = 1, nparticle_local
3853 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
3856 IF (do_efield0)
THEN
3857 efield0(iparticle) = efield0(iparticle)
3860 IF (do_efield1)
THEN
3861 efield1(1:3, iparticle) = efield1(1:3, iparticle) - (1.0_dp/6.0_dp)*dipole_t(1:3)
3864 IF (do_efield2)
THEN
3865 efield2(1, iparticle) = efield2(1, iparticle) - (1.0_dp/(18.0_dp))*q_t
3866 efield2(5, iparticle) = efield2(5, iparticle) - (1.0_dp/(18.0_dp))*q_t
3867 efield2(9, iparticle) = efield2(9, iparticle) - (1.0_dp/(18.0_dp))*q_t
3874 gauss = (rho0(lp, mp, np)*pw_grid%vol)**2/pw_grid%gsq(gpt)
3875 factor = gauss*real(summe_ef(gpt)*conjg(summe_ef(gpt)), kind=
dp)
3876 energy = energy + factor
3878 IF (do_forces .OR. do_efield)
THEN
3880 DO iparticle_kind = 1, nparticle_kind
3881 nparticle_local = local_particles%n_el(iparticle_kind)
3882 DO iparticle_local = 1, nparticle_local
3884 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
3885 fac = exp_igr%ex(lp, node)*exp_igr%ey(mp, node)*exp_igr%ez(np, node)
3886 cnjg_fac = conjg(
fac)
3890 CALL get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges, &
3891 dipoles, quadrupoles)
3893 tmp = gauss*aimag(summe_ef(gpt)*(cnjg_fac*conjg(atm_factor)))
3894 forces(1, node) = forces(1, node) + tmp*pw_grid%g(1, gpt)
3895 forces(2, node) = forces(2, node) + tmp*pw_grid%g(2, gpt)
3896 forces(3, node) = forces(3, node) + tmp*pw_grid%g(3, gpt)
3902 IF (do_efield0)
THEN
3903 efield0(iparticle) = efield0(iparticle) + gauss*real(
fac*conjg(summe_ef(gpt)), kind=
dp)
3906 IF (do_efield1)
THEN
3907 tmp = aimag(
fac*conjg(summe_ef(gpt)))*gauss
3908 efield1(1, iparticle) = efield1(1, iparticle) - tmp*pw_grid%g(1, gpt)
3909 efield1(2, iparticle) = efield1(2, iparticle) - tmp*pw_grid%g(2, gpt)
3910 efield1(3, iparticle) = efield1(3, iparticle) - tmp*pw_grid%g(3, gpt)
3913 IF (do_efield2)
THEN
3914 tmp_v(1) = real(
fac*conjg(summe_ef(gpt)), kind=
dp)*pw_grid%g(1, gpt)*gauss
3915 tmp_v(2) = real(
fac*conjg(summe_ef(gpt)), kind=
dp)*pw_grid%g(2, gpt)*gauss
3916 tmp_v(3) = real(
fac*conjg(summe_ef(gpt)), kind=
dp)*pw_grid%g(3, gpt)*gauss
3918 efield2(1, iparticle) = efield2(1, iparticle) + tmp_v(1)*pw_grid%g(1, gpt)
3919 efield2(2, iparticle) = efield2(2, iparticle) + tmp_v(1)*pw_grid%g(2, gpt)
3920 efield2(3, iparticle) = efield2(3, iparticle) + tmp_v(1)*pw_grid%g(3, gpt)
3921 efield2(4, iparticle) = efield2(4, iparticle) + tmp_v(2)*pw_grid%g(1, gpt)
3922 efield2(5, iparticle) = efield2(5, iparticle) + tmp_v(2)*pw_grid%g(2, gpt)
3923 efield2(6, iparticle) = efield2(6, iparticle) + tmp_v(2)*pw_grid%g(3, gpt)
3924 efield2(7, iparticle) = efield2(7, iparticle) + tmp_v(3)*pw_grid%g(1, gpt)
3925 efield2(8, iparticle) = efield2(8, iparticle) + tmp_v(3)*pw_grid%g(2, gpt)
3926 efield2(9, iparticle) = efield2(9, iparticle) + tmp_v(3)*pw_grid%g(3, gpt)
3937 denom = 1.0_dp/four_alpha_sq + 1.0_dp/pw_grid%gsq(gpt)
3938 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)
3939 pv_tmp(1, 2) = pv_tmp(1, 2) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(2, gpt)*denom)
3940 pv_tmp(1, 3) = pv_tmp(1, 3) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(3, gpt)*denom)
3941 pv_tmp(2, 1) = pv_tmp(2, 1) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(1, gpt)*denom)
3942 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)
3943 pv_tmp(2, 3) = pv_tmp(2, 3) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(3, gpt)*denom)
3944 pv_tmp(3, 1) = pv_tmp(3, 1) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(1, gpt)*denom)
3945 pv_tmp(3, 2) = pv_tmp(3, 2) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(2, gpt)*denom)
3946 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)
3949 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)
3950 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)
3951 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)
3952 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)
3953 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)
3954 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)
3955 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)
3956 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)
3957 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)
3960 pref =
fourpi/pw_grid%vol
3961 energy = energy*pref
3964 DEALLOCATE (summe_ef)
3966 pv_tmp = pv_tmp*pref
3968 pv(1, 1) = pv(1, 1) + pv_tmp(1, 1)
3969 pv(1, 2) = pv(1, 2) + (pv_tmp(1, 2) + pv_tmp(2, 1))*0.5_dp
3970 pv(1, 3) = pv(1, 3) + (pv_tmp(1, 3) + pv_tmp(3, 1))*0.5_dp
3972 pv(2, 2) = pv(2, 2) + pv_tmp(2, 2)
3973 pv(2, 3) = pv(2, 3) + (pv_tmp(2, 3) + pv_tmp(3, 2))*0.5_dp
3976 pv(3, 3) = pv(3, 3) + pv_tmp(3, 3)
3977 DEALLOCATE (summe_st)
3980 forces = 2.0_dp*forces*pref
3982 IF (do_efield0)
THEN
3983 efield0 = 2.0_dp*efield0*pref
3985 IF (do_efield1)
THEN
3986 efield1 = 2.0_dp*efield1*pref
3988 IF (do_efield2)
THEN
3989 efield2 = 2.0_dp*efield2*pref
3991 CALL timestop(handle)
3993 END SUBROUTINE ewald_multipole_lr
4009 SUBROUTINE get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges, &
4010 dipoles, quadrupoles)
4011 COMPLEX(KIND=dp),
INTENT(OUT) :: atm_factor
4013 INTEGER,
INTENT(IN) :: gpt
4014 INTEGER :: iparticle
4015 LOGICAL,
DIMENSION(3, 3),
INTENT(IN) :: task
4016 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: charges
4017 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
4018 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
4019 POINTER :: quadrupoles
4021 COMPLEX(KIND=dp) :: tmp
4025 IF (task(1, 1))
THEN
4027 atm_factor = atm_factor + charges(iparticle)
4029 IF (task(2, 2))
THEN
4033 tmp = tmp + dipoles(i, iparticle)*pw_grid%g(i, gpt)
4035 atm_factor = atm_factor + tmp*cmplx(0.0_dp, -1.0_dp, kind=
dp)
4037 IF (task(3, 3))
THEN
4042 tmp = tmp + quadrupoles(j, i, iparticle)*pw_grid%g(j, gpt)*pw_grid%g(i, gpt)
4045 atm_factor = atm_factor - 1.0_dp/3.0_dp*tmp
4048 END SUBROUTINE get_atom_factor
4063 SUBROUTINE get_atom_factor_stress(atm_factor, pw_grid, gpt, iparticle, task, &
4064 dipoles, quadrupoles)
4065 COMPLEX(KIND=dp),
INTENT(OUT) :: atm_factor(3)
4067 INTEGER,
INTENT(IN) :: gpt
4068 INTEGER :: iparticle
4069 LOGICAL,
DIMENSION(3, 3),
INTENT(IN) :: task
4070 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
4071 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
4072 POINTER :: quadrupoles
4077 IF (any(task(2, :)))
THEN
4079 atm_factor = dipoles(:, iparticle)*cmplx(0.0_dp, -1.0_dp, kind=
dp)
4081 IF (any(task(3, :)))
THEN
4084 atm_factor(1) = atm_factor(1) - 1.0_dp/3.0_dp* &
4085 (quadrupoles(1, i, iparticle)*pw_grid%g(i, gpt) + &
4086 quadrupoles(i, 1, iparticle)*pw_grid%g(i, gpt))
4087 atm_factor(2) = atm_factor(2) - 1.0_dp/3.0_dp* &
4088 (quadrupoles(2, i, iparticle)*pw_grid%g(i, gpt) + &
4089 quadrupoles(i, 2, iparticle)*pw_grid%g(i, gpt))
4090 atm_factor(3) = atm_factor(3) - 1.0_dp/3.0_dp* &
4091 (quadrupoles(3, i, iparticle)*pw_grid%g(i, gpt) + &
4092 quadrupoles(i, 3, iparticle)*pw_grid%g(i, gpt))
4096 END SUBROUTINE get_atom_factor_stress
4117 SUBROUTINE ewald_multipole_self(ewald_env, cell, local_particles, e_self, &
4118 e_neut, task, do_efield, radii, charges, dipoles, quadrupoles, efield0, &
4123 REAL(kind=
dp),
INTENT(OUT) :: e_self, e_neut
4124 LOGICAL,
DIMENSION(3, 3),
INTENT(IN) :: task
4125 LOGICAL,
INTENT(IN) :: do_efield
4126 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: radii, charges
4127 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
4128 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
4129 POINTER :: quadrupoles
4130 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0
4131 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: efield1, efield2
4133 REAL(kind=
dp),
PARAMETER :: f23 = 2.0_dp/3.0_dp, &
4134 f415 = 4.0_dp/15.0_dp
4136 INTEGER :: ewald_type, i, iparticle, &
4137 iparticle_kind, iparticle_local, j, &
4139 LOGICAL :: do_efield0, do_efield1, do_efield2, &
4141 REAL(kind=
dp) :: alpha, ch_qu_self, ch_qu_self_tmp, &
4142 dipole_self, fac1, fac2, fac3, fac4, &
4143 q, q_neutg, q_self, q_sum, qu_qu_self, &
4147 CALL ewald_env_get(ewald_env, ewald_type=ewald_type, alpha=alpha, &
4150 do_efield0 = do_efield .AND.
ASSOCIATED(efield0)
4151 do_efield1 = do_efield .AND.
ASSOCIATED(efield1)
4152 do_efield2 = do_efield .AND.
ASSOCIATED(efield2)
4155 dipole_self = 0.0_dp
4159 fac2 = 6.0_dp*(f23**2)*(alpha**3)*
oorootpi
4160 fac3 = (2.0_dp*
oorootpi)*f23*alpha**3
4161 fac4 = (4.0_dp*
oorootpi)*f415*alpha**5
4162 lradii =
PRESENT(radii)
4165 DO iparticle_kind = 1,
SIZE(local_particles%n_el)
4166 nparticle_local = local_particles%n_el(iparticle_kind)
4167 DO iparticle_local = 1, nparticle_local
4168 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
4169 IF (any(task(1, :)))
THEN
4171 q = charges(iparticle)
4172 IF (lradii) radius = radii(iparticle)
4173 IF (radius > 0)
THEN
4174 q_neutg = q_neutg + 2.0_dp*q*radius**2
4176 q_self = q_self + q*q
4179 IF (do_efield0)
THEN
4180 efield0(iparticle) = efield0(iparticle) - q*fac1
4183 IF (task(1, 3))
THEN
4185 ch_qu_self_tmp = 0.0_dp
4187 ch_qu_self_tmp = ch_qu_self_tmp + quadrupoles(i, i, iparticle)*q
4189 ch_qu_self = ch_qu_self + ch_qu_self_tmp
4191 IF (do_efield2)
THEN
4192 efield2(1, iparticle) = efield2(1, iparticle) + fac2*q
4193 efield2(5, iparticle) = efield2(5, iparticle) + fac2*q
4194 efield2(9, iparticle) = efield2(9, iparticle) + fac2*q
4198 IF (any(task(2, :)))
THEN
4201 dipole_self = dipole_self + dipoles(i, iparticle)**2
4204 IF (do_efield1)
THEN
4205 efield1(1, iparticle) = efield1(1, iparticle) + fac3*dipoles(1, iparticle)
4206 efield1(2, iparticle) = efield1(2, iparticle) + fac3*dipoles(2, iparticle)
4207 efield1(3, iparticle) = efield1(3, iparticle) + fac3*dipoles(3, iparticle)
4210 IF (any(task(3, :)))
THEN
4214 qu_qu_self = qu_qu_self + quadrupoles(j, i, iparticle)**2
4218 IF (do_efield2)
THEN
4219 efield2(1, iparticle) = efield2(1, iparticle) + fac4*quadrupoles(1, 1, iparticle)
4220 efield2(2, iparticle) = efield2(2, iparticle) + fac4*quadrupoles(2, 1, iparticle)
4221 efield2(3, iparticle) = efield2(3, iparticle) + fac4*quadrupoles(3, 1, iparticle)
4222 efield2(4, iparticle) = efield2(4, iparticle) + fac4*quadrupoles(1, 2, iparticle)
4223 efield2(5, iparticle) = efield2(5, iparticle) + fac4*quadrupoles(2, 2, iparticle)
4224 efield2(6, iparticle) = efield2(6, iparticle) + fac4*quadrupoles(3, 2, iparticle)
4225 efield2(7, iparticle) = efield2(7, iparticle) + fac4*quadrupoles(1, 3, iparticle)
4226 efield2(8, iparticle) = efield2(8, iparticle) + fac4*quadrupoles(2, 3, iparticle)
4227 efield2(9, iparticle) = efield2(9, iparticle) + fac4*quadrupoles(3, 3, iparticle)
4233 CALL group%sum(q_neutg)
4234 CALL group%sum(q_self)
4235 CALL group%sum(q_sum)
4236 CALL group%sum(dipole_self)
4237 CALL group%sum(ch_qu_self)
4238 CALL group%sum(qu_qu_self)
4240 e_self = -(q_self + f23*(dipole_self - f23*ch_qu_self + f415*qu_qu_self*alpha**2)*alpha**2)*alpha*
oorootpi
4241 fac1 =
pi/(2.0_dp*cell%deth)
4242 e_neut = -q_sum*fac1*(q_sum/alpha**2 - q_neutg)
4245 DO iparticle_kind = 1,
SIZE(local_particles%n_el)
4246 nparticle_local = local_particles%n_el(iparticle_kind)
4247 DO iparticle_local = 1, nparticle_local
4248 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
4249 IF (any(task(1, :)))
THEN
4251 IF (do_efield0)
THEN
4252 efield0(iparticle) = efield0(iparticle) - q_sum*2.0_dp*fac1/alpha**2
4253 IF (lradii) radius = radii(iparticle)
4254 IF (radius > 0)
THEN
4255 q = charges(iparticle)
4256 efield0(iparticle) = efield0(iparticle) + fac1*radius**2*(q_sum + q)
4263 END SUBROUTINE ewald_multipole_self
4275 SUBROUTINE ewald_multipole_print(iw, e_gspace, e_rspace, e_bonded, e_self, e_neut)
4277 INTEGER,
INTENT(IN) :: iw
4278 REAL(kind=
dp),
INTENT(IN) :: e_gspace, e_rspace, e_bonded, e_self, &
4282 WRITE (iw,
'( A, A )')
' *********************************', &
4283 '**********************************************'
4284 WRITE (iw,
'( A, A, T35, A, T56, E25.15 )')
' INITIAL GSPACE ENERGY', &
4285 '[hartree]',
'= ', e_gspace
4286 WRITE (iw,
'( A, A, T35, A, T56, E25.15 )')
' INITIAL RSPACE ENERGY', &
4287 '[hartree]',
'= ', e_rspace
4288 WRITE (iw,
'( A, A, T35, A, T56, E25.15 )')
' BONDED CORRECTION', &
4289 '[hartree]',
'= ', e_bonded
4290 WRITE (iw,
'( A, A, T35, A, T56, E25.15 )')
' SELF ENERGY CORRECTION', &
4291 '[hartree]',
'= ', e_self
4292 WRITE (iw,
'( A, A, T35, A, T56, E25.15 )')
' NEUTRALIZ. BCKGR. ENERGY', &
4293 '[hartree]',
'= ', e_neut
4294 WRITE (iw,
'( A, A, T35, A, T56, E25.15 )')
' TOTAL ELECTROSTATIC EN.', &
4295 '[hartree]',
'= ', e_rspace + e_bonded + e_gspace + e_self + e_neut
4296 WRITE (iw,
'( A, A )')
' *********************************', &
4297 '**********************************************'
4299 END SUBROUTINE ewald_multipole_print
4314 SUBROUTINE debug_ewald_multipoles(ewald_env, ewald_pw, nonbond_env, cell, &
4315 particle_set, local_particles, iw, debug_r_space)
4316 TYPE charge_mono_type
4317 REAL(kind=
dp),
DIMENSION(:), &
4319 REAL(kind=
dp),
DIMENSION(:, :), &
4321 END TYPE charge_mono_type
4322 TYPE multi_charge_type
4323 TYPE(charge_mono_type),
DIMENSION(:), &
4324 POINTER :: charge_typ
4325 END TYPE multi_charge_type
4331 POINTER :: particle_set
4333 INTEGER,
INTENT(IN) :: iw
4334 LOGICAL,
INTENT(IN) :: debug_r_space
4336 INTEGER :: nparticles
4337 LOGICAL,
DIMENSION(3) :: task
4338 REAL(kind=
dp) :: e_neut, e_self, g_energy, &
4339 r_energy, debug_energy
4340 REAL(kind=
dp),
POINTER,
DIMENSION(:) :: charges
4341 REAL(kind=
dp),
POINTER, &
4342 DIMENSION(:, :) :: dipoles, g_forces, g_pv, &
4343 r_forces, r_pv, e_field1, &
4345 REAL(kind=
dp),
POINTER, &
4346 DIMENSION(:, :, :) :: quadrupoles
4348 TYPE(multi_charge_type),
DIMENSION(:), &
4349 POINTER :: multipoles
4351 NULLIFY (multipoles, charges, dipoles, g_forces, g_pv, &
4352 r_forces, r_pv, e_field1, e_field2)
4357 nparticles =
SIZE(particle_set)
4360 ALLOCATE (charges(nparticles))
4361 ALLOCATE (dipoles(3, nparticles))
4362 ALLOCATE (quadrupoles(3, 3, nparticles))
4365 ALLOCATE (r_forces(3, nparticles))
4366 ALLOCATE (g_forces(3, nparticles))
4367 ALLOCATE (e_field1(3, nparticles))
4368 ALLOCATE (e_field2(3, nparticles))
4369 ALLOCATE (g_pv(3, 3))
4370 ALLOCATE (r_pv(3, 3))
4376 quadrupoles = 0.0_dp
4388 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2,
"CHARGE", echarge=-1.0_dp, &
4389 random_stream=random_stream, charges=charges)
4390 CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles,
"CHARGE", echarge=1.0_dp, &
4391 random_stream=random_stream, charges=charges)
4392 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
4395 WRITE (iw, *)
"DEBUG ENERGY (CHARGE-CHARGE): ", debug_energy
4397 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
4398 task, do_correction_bonded=.false., do_forces=.true., do_stress=.true., do_efield=.false., &
4399 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
4400 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.false.)
4401 CALL release_multi_type(multipoles)
4408 quadrupoles = 0.0_dp
4420 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2,
"CHARGE", echarge=-1.0_dp, &
4421 random_stream=random_stream, charges=charges)
4422 CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles,
"DIPOLE", echarge=0.5_dp, &
4423 random_stream=random_stream, dipoles=dipoles)
4424 WRITE (iw,
'("CHARGES",F15.9)') charges
4425 WRITE (iw,
'("DIPOLES",3F15.9)') dipoles
4426 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
4429 WRITE (iw, *)
"DEBUG ENERGY (CHARGE-DIPOLE): ", debug_energy
4431 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
4432 task, do_correction_bonded=.false., do_forces=.true., do_stress=.true., do_efield=.false., &
4433 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
4434 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.false.)
4435 CALL release_multi_type(multipoles)
4441 quadrupoles = 0.0_dp
4453 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2,
"DIPOLE", echarge=10000.0_dp, &
4454 random_stream=random_stream, dipoles=dipoles)
4455 CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles,
"DIPOLE", echarge=20000._dp, &
4456 random_stream=random_stream, dipoles=dipoles)
4457 WRITE (iw,
'("DIPOLES",3F15.9)') dipoles
4458 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
4461 WRITE (iw, *)
"DEBUG ENERGY (DIPOLE-DIPOLE): ", debug_energy
4463 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
4464 task, do_correction_bonded=.false., do_forces=.true., do_stress=.true., do_efield=.false., &
4465 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
4466 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.false.)
4467 CALL release_multi_type(multipoles)
4474 quadrupoles = 0.0_dp
4486 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2,
"CHARGE", echarge=-1.0_dp, &
4487 random_stream=random_stream, charges=charges)
4488 CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles,
"QUADRUPOLE", echarge=10.0_dp, &
4489 random_stream=random_stream, quadrupoles=quadrupoles)
4490 WRITE (iw,
'("CHARGES",F15.9)') charges
4491 WRITE (iw,
'("QUADRUPOLES",9F15.9)') quadrupoles
4492 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
4495 WRITE (iw, *)
"DEBUG ENERGY (CHARGE-QUADRUPOLE): ", debug_energy
4497 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
4498 task, do_correction_bonded=.false., do_forces=.true., do_stress=.true., do_efield=.false., &
4499 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
4500 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.false.)
4501 CALL release_multi_type(multipoles)
4508 quadrupoles = 0.0_dp
4520 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2,
"DIPOLE", echarge=10000.0_dp, &
4521 random_stream=random_stream, dipoles=dipoles)
4522 CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles,
"QUADRUPOLE", echarge=10000.0_dp, &
4523 random_stream=random_stream, quadrupoles=quadrupoles)
4524 WRITE (iw,
'("DIPOLES",3F15.9)') dipoles
4525 WRITE (iw,
'("QUADRUPOLES",9F15.9)') quadrupoles
4526 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
4529 WRITE (iw, *)
"DEBUG ENERGY (DIPOLE-QUADRUPOLE): ", debug_energy
4531 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
4532 task, do_correction_bonded=.false., do_forces=.true., do_stress=.true., do_efield=.false., &
4533 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
4534 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.false.)
4535 CALL release_multi_type(multipoles)
4541 quadrupoles = 0.0_dp
4553 CALL create_multi_type(multipoles, nparticles, 1, nparticles/2,
"QUADRUPOLE", echarge=-20000.0_dp, &
4554 random_stream=random_stream, quadrupoles=quadrupoles)
4555 CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles,
"QUADRUPOLE", echarge=10000.0_dp, &
4556 random_stream=random_stream, quadrupoles=quadrupoles)
4557 WRITE (iw,
'("QUADRUPOLES",9F15.9)') quadrupoles
4558 CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
4561 WRITE (iw, *)
"DEBUG ENERGY (QUADRUPOLE-QUADRUPOLE): ", debug_energy
4563 particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
4564 task, do_correction_bonded=.false., do_forces=.true., do_stress=.true., do_efield=.false., &
4565 charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
4566 forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.false.)
4567 CALL release_multi_type(multipoles)
4569 DEALLOCATE (charges)
4570 DEALLOCATE (dipoles)
4571 DEALLOCATE (quadrupoles)
4572 DEALLOCATE (r_forces)
4573 DEALLOCATE (g_forces)
4574 DEALLOCATE (e_field1)
4575 DEALLOCATE (e_field2)
4591 SUBROUTINE debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, &
4592 energy, debug_r_space)
4596 TYPE(multi_charge_type),
DIMENSION(:),
POINTER :: multipoles
4597 REAL(kind=
dp),
INTENT(OUT) :: energy
4598 LOGICAL,
INTENT(IN) :: debug_r_space
4600 INTEGER :: atom_a, atom_b, icell, iend, igrp, &
4601 ikind, ilist, ipair, istart, jcell, &
4602 jkind, k, k1, kcell, l, l1, ncells, &
4604 INTEGER,
DIMENSION(:, :),
POINTER ::
list
4605 REAL(kind=
dp) :: fac_ij, q, r, rab2, rab2_max
4606 REAL(kind=
dp),
DIMENSION(3) :: cell_v, cvi, rab, rab0, rm
4609 TYPE(
pos_type),
DIMENSION(:),
POINTER :: r_last_update, r_last_update_pbc
4613 r_last_update=r_last_update, r_last_update_pbc=r_last_update_pbc)
4614 rab2_max = huge(0.0_dp)
4615 IF (debug_r_space)
THEN
4618 lists:
DO ilist = 1, nonbonded%nlists
4619 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
4620 npairs = neighbor_kind_pair%npairs
4621 IF (npairs == 0) cycle
4622 list => neighbor_kind_pair%list
4623 cvi = neighbor_kind_pair%cell_vector
4624 cell_v = matmul(cell%hmat, cvi)
4625 kind_group_loop:
DO igrp = 1, neighbor_kind_pair%ngrp_kind
4626 istart = neighbor_kind_pair%grp_kind_start(igrp)
4627 iend = neighbor_kind_pair%grp_kind_end(igrp)
4628 ikind = neighbor_kind_pair%ij_kind(1, igrp)
4629 jkind = neighbor_kind_pair%ij_kind(2, igrp)
4630 pairs:
DO ipair = istart, iend
4632 atom_a =
list(1, ipair)
4633 atom_b =
list(2, ipair)
4634 IF (atom_a == atom_b) fac_ij = 0.5_dp
4635 rab = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
4637 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
4638 IF (rab2 <= rab2_max)
THEN
4640 DO k = 1,
SIZE(multipoles(atom_a)%charge_typ)
4641 DO k1 = 1,
SIZE(multipoles(atom_a)%charge_typ(k)%charge)
4643 DO l = 1,
SIZE(multipoles(atom_b)%charge_typ)
4644 DO l1 = 1,
SIZE(multipoles(atom_b)%charge_typ(l)%charge)
4646 rm = rab + multipoles(atom_b)%charge_typ(l)%pos(:, l1) - multipoles(atom_a)%charge_typ(k)%pos(:, k1)
4647 r = sqrt(dot_product(rm, rm))
4648 q = multipoles(atom_b)%charge_typ(l)%charge(l1)*multipoles(atom_a)%charge_typ(k)%charge(k1)
4649 energy = energy + q/r*fac_ij
4658 END DO kind_group_loop
4664 DO atom_a = 1,
SIZE(particle_set)
4665 DO atom_b = atom_a,
SIZE(particle_set)
4667 IF (atom_a == atom_b) fac_ij = 0.5_dp
4668 rab0 = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
4670 DO icell = -ncells, ncells
4671 DO jcell = -ncells, ncells
4672 DO kcell = -ncells, ncells
4673 cell_v = matmul(cell%hmat, real((/icell, jcell, kcell/), kind=
dp))
4674 IF (all(cell_v == 0.0_dp) .AND. (atom_a == atom_b)) cycle
4676 rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
4677 IF (rab2 <= rab2_max)
THEN
4679 DO k = 1,
SIZE(multipoles(atom_a)%charge_typ)
4680 DO k1 = 1,
SIZE(multipoles(atom_a)%charge_typ(k)%charge)
4682 DO l = 1,
SIZE(multipoles(atom_b)%charge_typ)
4683 DO l1 = 1,
SIZE(multipoles(atom_b)%charge_typ(l)%charge)
4685 rm = rab + multipoles(atom_b)%charge_typ(l)%pos(:, l1) - multipoles(atom_a)%charge_typ(k)%pos(:, k1)
4686 r = sqrt(dot_product(rm, rm))
4687 q = multipoles(atom_b)%charge_typ(l)%charge(l1)*multipoles(atom_a)%charge_typ(k)%charge(k1)
4688 energy = energy + q/r*fac_ij
4702 END SUBROUTINE debug_ewald_multipole_low
4719 SUBROUTINE create_multi_type(multipoles, idim, istart, iend, label, echarge, &
4720 random_stream, charges, dipoles, quadrupoles)
4721 TYPE(multi_charge_type),
DIMENSION(:),
POINTER :: multipoles
4722 INTEGER,
INTENT(IN) :: idim, istart, iend
4723 CHARACTER(LEN=*),
INTENT(IN) :: label
4724 REAL(kind=
dp),
INTENT(IN) :: echarge
4726 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: charges
4727 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
4728 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
4729 POINTER :: quadrupoles
4731 INTEGER :: i, isize, k, l, m
4732 REAL(kind=
dp) :: dx, r2, rvec(3), rvec1(3), rvec2(3)
4734 IF (
ASSOCIATED(multipoles))
THEN
4735 cpassert(
SIZE(multipoles) == idim)
4737 ALLOCATE (multipoles(idim))
4739 NULLIFY (multipoles(i)%charge_typ)
4743 IF (
ASSOCIATED(multipoles(i)%charge_typ))
THEN
4745 isize =
SIZE(multipoles(i)%charge_typ) + 1
4749 CALL reallocate_charge_type(multipoles(i)%charge_typ, 1, isize)
4752 cpassert(
PRESENT(charges))
4753 cpassert(
ASSOCIATED(charges))
4754 ALLOCATE (multipoles(i)%charge_typ(isize)%charge(1))
4755 ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 1))
4757 multipoles(i)%charge_typ(isize)%charge(1) = echarge
4758 multipoles(i)%charge_typ(isize)%pos(1:3, 1) = 0.0_dp
4759 charges(i) = charges(i) + echarge
4762 cpassert(
PRESENT(dipoles))
4763 cpassert(
ASSOCIATED(dipoles))
4764 ALLOCATE (multipoles(i)%charge_typ(isize)%charge(2))
4765 ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 2))
4766 CALL random_stream%fill(rvec)
4767 rvec = rvec/(2.0_dp*sqrt(dot_product(rvec, rvec)))*dx
4768 multipoles(i)%charge_typ(isize)%charge(1) = echarge
4769 multipoles(i)%charge_typ(isize)%pos(1:3, 1) = rvec
4770 multipoles(i)%charge_typ(isize)%charge(2) = -echarge
4771 multipoles(i)%charge_typ(isize)%pos(1:3, 2) = -rvec
4773 dipoles(:, i) = dipoles(:, i) + 2.0_dp*echarge*rvec
4776 cpassert(
PRESENT(quadrupoles))
4777 cpassert(
ASSOCIATED(quadrupoles))
4778 ALLOCATE (multipoles(i)%charge_typ(isize)%charge(4))
4779 ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 4))
4780 CALL random_stream%fill(rvec1)
4781 CALL random_stream%fill(rvec2)
4782 rvec1 = rvec1/sqrt(dot_product(rvec1, rvec1))
4783 rvec2 = rvec2 - dot_product(rvec2, rvec1)*rvec1
4784 rvec2 = rvec2/sqrt(dot_product(rvec2, rvec2))
4786 rvec1 = rvec1/2.0_dp*dx
4787 rvec2 = rvec2/2.0_dp*dx
4795 multipoles(i)%charge_typ(isize)%charge(1) = -echarge
4796 multipoles(i)%charge_typ(isize)%pos(1:3, 1) = rvec1 + rvec2
4797 multipoles(i)%charge_typ(isize)%charge(2) = echarge
4798 multipoles(i)%charge_typ(isize)%pos(1:3, 2) = rvec1 - rvec2
4799 multipoles(i)%charge_typ(isize)%charge(3) = -echarge
4800 multipoles(i)%charge_typ(isize)%pos(1:3, 3) = -rvec1 - rvec2
4801 multipoles(i)%charge_typ(isize)%charge(4) = echarge
4802 multipoles(i)%charge_typ(isize)%pos(1:3, 4) = -rvec1 + rvec2
4805 r2 = dot_product(multipoles(i)%charge_typ(isize)%pos(:, k), multipoles(i)%charge_typ(isize)%pos(:, k))
4808 quadrupoles(m, l, i) = quadrupoles(m, l, i) + 3.0_dp*0.5_dp*multipoles(i)%charge_typ(isize)%charge(k)* &
4809 multipoles(i)%charge_typ(isize)%pos(l, k)* &
4810 multipoles(i)%charge_typ(isize)%pos(m, k)
4811 IF (m == l) quadrupoles(m, l, i) = quadrupoles(m, l, i) - 0.5_dp*multipoles(i)%charge_typ(isize)%charge(k)*r2
4818 END SUBROUTINE create_multi_type
4826 SUBROUTINE release_multi_type(multipoles)
4827 TYPE(multi_charge_type),
DIMENSION(:),
POINTER :: multipoles
4831 IF (
ASSOCIATED(multipoles))
THEN
4832 DO i = 1,
SIZE(multipoles)
4833 DO j = 1,
SIZE(multipoles(i)%charge_typ)
4834 DEALLOCATE (multipoles(i)%charge_typ(j)%charge)
4835 DEALLOCATE (multipoles(i)%charge_typ(j)%pos)
4837 DEALLOCATE (multipoles(i)%charge_typ)
4840 END SUBROUTINE release_multi_type
4850 SUBROUTINE reallocate_charge_type(charge_typ, istart, iend)
4851 TYPE(charge_mono_type),
DIMENSION(:),
POINTER :: charge_typ
4852 INTEGER,
INTENT(IN) :: istart, iend
4854 INTEGER :: i, isize, j, jsize, jsize1, jsize2
4855 TYPE(charge_mono_type),
DIMENSION(:),
POINTER :: charge_typ_bk
4857 IF (
ASSOCIATED(charge_typ))
THEN
4858 isize =
SIZE(charge_typ)
4859 ALLOCATE (charge_typ_bk(1:isize))
4861 jsize =
SIZE(charge_typ(j)%charge)
4862 ALLOCATE (charge_typ_bk(j)%charge(jsize))
4863 jsize1 =
SIZE(charge_typ(j)%pos, 1)
4864 jsize2 =
SIZE(charge_typ(j)%pos, 2)
4865 ALLOCATE (charge_typ_bk(j)%pos(jsize1, jsize2))
4866 charge_typ_bk(j)%pos = charge_typ(j)%pos
4867 charge_typ_bk(j)%charge = charge_typ(j)%charge
4869 DO j = 1,
SIZE(charge_typ)
4870 DEALLOCATE (charge_typ(j)%charge)
4871 DEALLOCATE (charge_typ(j)%pos)
4873 DEALLOCATE (charge_typ)
4875 ALLOCATE (charge_typ_bk(istart:iend))
4876 DO i = istart, isize
4877 jsize =
SIZE(charge_typ_bk(j)%charge)
4878 ALLOCATE (charge_typ(j)%charge(jsize))
4879 jsize1 =
SIZE(charge_typ_bk(j)%pos, 1)
4880 jsize2 =
SIZE(charge_typ_bk(j)%pos, 2)
4881 ALLOCATE (charge_typ(j)%pos(jsize1, jsize2))
4882 charge_typ(j)%pos = charge_typ_bk(j)%pos
4883 charge_typ(j)%charge = charge_typ_bk(j)%charge
4885 DO j = 1,
SIZE(charge_typ_bk)
4886 DEALLOCATE (charge_typ_bk(j)%charge)
4887 DEALLOCATE (charge_typ_bk(j)%pos)
4889 DEALLOCATE (charge_typ_bk)
4891 ALLOCATE (charge_typ(istart:iend))
4894 END SUBROUTINE reallocate_charge_type
4896 END SUBROUTINE debug_ewald_multipoles
4917 SUBROUTINE debug_ewald_multipoles_fields(ewald_env, ewald_pw, nonbond_env, cell, &
4918 particle_set, local_particles, radii, charges, dipoles, quadrupoles, task, iw, &
4919 atomic_kind_set, mm_section)
4926 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: radii, charges
4927 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
4928 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
4929 POINTER :: quadrupoles
4930 LOGICAL,
DIMENSION(3),
INTENT(IN) :: task
4931 INTEGER,
INTENT(IN) :: iw
4935 INTEGER :: i, iparticle_kind, j, k, &
4936 nparticle_local, nparticles
4937 REAL(kind=
dp) :: coord(3), dq, e_neut, e_self, efield1n(3), efield2n(3, 3), ene(2), &
4938 energy_glob, energy_local, enev(3, 2), o_tot_ene, pot, pv_glob(3, 3), pv_local(3, 3), &
4940 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: efield1, efield2, forces_glob, &
4942 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0, lcharges
4944 TYPE(
particle_type),
DIMENSION(:),
POINTER :: core_particle_set, shell_particle_set
4946 NULLIFY (lcharges, shell_particle_set, core_particle_set)
4950 nparticles =
SIZE(particle_set)
4952 DO iparticle_kind = 1,
SIZE(local_particles%n_el)
4953 nparticle_local = nparticle_local + local_particles%n_el(iparticle_kind)
4955 ALLOCATE (lcharges(nparticles))
4956 ALLOCATE (forces_glob(3, nparticles))
4957 ALLOCATE (forces_local(3, nparticle_local))
4958 ALLOCATE (efield0(nparticles))
4959 ALLOCATE (efield1(3, nparticles))
4960 ALLOCATE (efield2(9, nparticles))
4961 forces_glob = 0.0_dp
4962 forces_local = 0.0_dp
4968 energy_glob = 0.0_dp
4969 energy_local = 0.0_dp
4973 local_particles, energy_local, energy_glob, e_neut, e_self, task, .false., .true., .true., &
4974 .true., radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob, &
4975 efield0, efield1, efield2, iw, do_debug=.false.)
4976 o_tot_ene = energy_local + energy_glob + e_neut + e_self
4977 WRITE (iw, *)
"TOTAL ENERGY :: ========>", o_tot_ene
4981 DO i = 1, nparticles
4984 lcharges(i) = charges(i) + (-1.0_dp)**k*dq
4985 forces_glob = 0.0_dp
4986 forces_local = 0.0_dp
4989 energy_glob = 0.0_dp
4990 energy_local = 0.0_dp
4994 local_particles, energy_local, energy_glob, e_neut, e_self, &
4995 task, .false., .false., .false., .false., radii, &
4996 lcharges, dipoles, quadrupoles, iw=iw, do_debug=.false.)
4997 ene(k) = energy_local + energy_glob + e_neut + e_self
4999 pot = (ene(2) - ene(1))/(2.0_dp*dq)
5000 WRITE (iw,
'(A,I8,3(A,F15.9))')
"POTENTIAL FOR ATOM: ", i,
" NUMERICAL: ", pot,
" ANALYTICAL: ", efield0(i), &
5001 " ERROR: ", pot - efield0(i)
5002 tot_ene = tot_ene + 0.5_dp*efield0(i)*charges(i)
5004 WRITE (iw, *)
"ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
5005 WRITE (iw,
'(/,/,/)')
5008 DO i = 1, nparticles
5009 coord = particle_set(i)%r
5012 particle_set(i)%r(j) = coord(j) + (-1.0_dp)**k*dq
5015 CALL list_control(atomic_kind_set, particle_set, local_particles, &
5016 cell, nonbond_env, logger%para_env, mm_section, &
5017 shell_particle_set, core_particle_set)
5019 forces_glob = 0.0_dp
5020 forces_local = 0.0_dp
5023 energy_glob = 0.0_dp
5024 energy_local = 0.0_dp
5029 local_particles, energy_local, energy_glob, e_neut, e_self, &
5030 task, .false., .true., .true., .true., radii, &
5031 charges, dipoles, quadrupoles, forces_local, forces_glob, &
5032 pv_local, pv_glob, efield0, iw=iw, do_debug=.false.)
5034 particle_set(i)%r(j) = coord(j)
5036 efield1n(j) = -(ene(2) - ene(1))/(2.0_dp*dq)
5038 WRITE (iw,
'(/,A,I8)')
"FIELD FOR ATOM: ", i
5039 WRITE (iw,
'(A,3F15.9)')
" NUMERICAL: ", efield1n,
" ANALYTICAL: ", efield1(:, i), &
5040 " ERROR: ", efield1n - efield1(:, i)
5042 tot_ene = tot_ene - 0.5_dp*dot_product(efield1(:, i), dipoles(:, i))
5045 WRITE (iw, *)
"ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
5049 DO i = 1, nparticles
5050 coord = particle_set(i)%r
5053 particle_set(i)%r(j) = coord(j) + (-1.0_dp)**k*dq
5056 CALL list_control(atomic_kind_set, particle_set, local_particles, &
5057 cell, nonbond_env, logger%para_env, mm_section, &
5058 shell_particle_set, core_particle_set)
5060 forces_glob = 0.0_dp
5061 forces_local = 0.0_dp
5064 energy_glob = 0.0_dp
5065 energy_local = 0.0_dp
5070 local_particles, energy_local, energy_glob, e_neut, e_self, &
5071 task, .false., .true., .true., .true., radii, &
5072 charges, dipoles, quadrupoles, forces_local, forces_glob, &
5073 pv_local, pv_glob, efield1=efield1, iw=iw, do_debug=.false.)
5074 enev(:, k) = efield1(:, i)
5075 particle_set(i)%r(j) = coord(j)
5077 efield2n(:, j) = (enev(:, 2) - enev(:, 1))/(2.0_dp*dq)
5079 WRITE (iw,
'(/,A,I8)')
"FIELD GRADIENT FOR ATOM: ", i
5080 WRITE (iw,
'(A,9F15.9)')
" NUMERICAL: ", efield2n, &
5081 " ANALYTICAL: ", efield2(:, i), &
5082 " ERROR: ", reshape(efield2n, (/9/)) - efield2(:, i)
5084 END SUBROUTINE debug_ewald_multipoles_fields
5103 SUBROUTINE debug_ewald_multipoles_fields2(ewald_env, ewald_pw, nonbond_env, cell, &
5104 particle_set, local_particles, radii, charges, dipoles, quadrupoles, task, iw)
5111 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: radii, charges
5112 REAL(kind=
dp),
DIMENSION(:, :),
OPTIONAL,
POINTER :: dipoles
5113 REAL(kind=
dp),
DIMENSION(:, :, :),
OPTIONAL, &
5114 POINTER :: quadrupoles
5115 LOGICAL,
DIMENSION(3),
INTENT(IN) :: task
5116 INTEGER,
INTENT(IN) :: iw
5118 INTEGER :: i, ind, iparticle_kind, j, k, &
5119 nparticle_local, nparticles
5120 REAL(kind=
dp) :: e_neut, e_self, energy_glob, &
5121 energy_local, o_tot_ene, prod, &
5122 pv_glob(3, 3), pv_local(3, 3), tot_ene
5123 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: efield1, efield2, forces_glob, &
5125 REAL(kind=
dp),
DIMENSION(:),
POINTER :: efield0
5131 nparticles =
SIZE(particle_set)
5133 DO iparticle_kind = 1,
SIZE(local_particles%n_el)
5134 nparticle_local = nparticle_local + local_particles%n_el(iparticle_kind)
5136 ALLOCATE (forces_glob(3, nparticles))
5137 ALLOCATE (forces_local(3, nparticle_local))
5138 ALLOCATE (efield0(nparticles))
5139 ALLOCATE (efield1(3, nparticles))
5140 ALLOCATE (efield2(9, nparticles))
5141 forces_glob = 0.0_dp
5142 forces_local = 0.0_dp
5148 energy_glob = 0.0_dp
5149 energy_local = 0.0_dp
5153 local_particles, energy_local, energy_glob, e_neut, e_self, task, .false., .true., .true., &
5154 .true., radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob, &
5155 efield0, efield1, efield2, iw, do_debug=.false.)
5156 o_tot_ene = energy_local + energy_glob + e_neut + e_self
5157 WRITE (iw, *)
"TOTAL ENERGY :: ========>", o_tot_ene
5162 DO i = 1, nparticles
5163 tot_ene = tot_ene + 0.5_dp*efield0(i)*charges(i)
5165 WRITE (iw, *)
"ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
5166 WRITE (iw,
'(/,/,/)')
5171 DO i = 1, nparticles
5172 tot_ene = tot_ene - 0.5_dp*dot_product(efield1(:, i), dipoles(:, i))
5174 WRITE (iw, *)
"ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
5175 WRITE (iw,
'(/,/,/)')
5180 DO i = 1, nparticles
5186 prod = prod + efield2(ind, i)*quadrupoles(j, k, i)
5189 tot_ene = tot_ene - 0.5_dp*(1.0_dp/3.0_dp)*prod
5191 WRITE (iw, *)
"ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
5192 WRITE (iw,
'(/,/,/)')
5195 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, quip_data, nequip_data, allegro_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 ...