113 TYPE(debug_variables_type),
INTENT(INOUT), &
116 CHARACTER(len=*),
PARAMETER :: routinen =
'fist_calc_energy_force'
118 INTEGER :: do_ipol, ewald_type, fg_coulomb_size, handle, i, ii, ikind, iparticle_kind, &
119 iparticle_local, iw, iw2, j, natoms, nlocal_particles, node, nparticle_kind, &
120 nparticle_local, nshell, shell_index
121 LOGICAL :: do_multipoles, shell_model_ad, &
122 shell_present, use_virial
123 REAL(kind=
dp) :: ef_ener, fc, fs, mass, pot_bend, pot_bond, pot_imptors, pot_manybody, &
124 pot_nonbond, pot_opbend, pot_shell, pot_torsion, pot_urey_bradley, vg_coulomb, xdum1, &
126 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: ef_f, f_nonbond, f_total, fcore_nonbond, &
127 fcore_shell_bonded, fcore_total, fg_coulomb, fgcore_coulomb, fgshell_coulomb, &
128 fshell_nonbond, fshell_total
129 REAL(kind=
dp),
DIMENSION(3, 3) :: ef_pv, ident, pv_bc, pv_bend, pv_bond, &
130 pv_g, pv_imptors, pv_nonbond, &
131 pv_opbend, pv_torsion, pv_urey_bradley
132 REAL(kind=
dp),
DIMENSION(:),
POINTER :: e_coulomb
133 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: pv_coulomb
151 TYPE(
particle_type),
DIMENSION(:),
POINTER :: core_particle_set, particle_set, &
158 CALL timeset(routinen, handle)
160 NULLIFY (subsys, virial, atprop_env, para_env, force_env_section)
166 input=force_env_section)
172 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
173 NULLIFY (atomic_kind, atomic_kind_set, cell, ewald_pw, ewald_env, &
174 fist_nonbond_env, mm_section, local_molecules, local_particles, &
175 molecule_kind_set, molecule_set, particle_set, print_section, &
176 shell, shell_particle_set, core_particle_set, thermo, multipoles, &
177 e_coulomb, pv_coulomb)
185 CALL fist_env_get(fist_env, ewald_pw=ewald_pw, ewald_env=ewald_env, &
186 local_particles=local_particles, particle_set=particle_set, &
187 atomic_kind_set=atomic_kind_set, molecule_set=molecule_set, &
188 local_molecules=local_molecules, thermo=thermo, &
189 molecule_kind_set=molecule_kind_set, fist_nonbond_env=fist_nonbond_env, &
190 cell=cell, shell_model=shell_present, shell_model_ad=shell_model_ad, &
191 multipoles=multipoles, exclusions=exclusions, efield=efield)
193 CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, &
200 natoms =
SIZE(particle_set)
202 nparticle_kind =
SIZE(atomic_kind_set)
203 DO ikind = 1, nparticle_kind
204 nlocal_particles = nlocal_particles + local_particles%n_el(ikind)
207 ALLOCATE (f_nonbond(3, natoms))
211 IF (shell_present)
THEN
212 CALL fist_env_get(fist_env, shell_particle_set=shell_particle_set, &
213 core_particle_set=core_particle_set)
214 cpassert(
ASSOCIATED(shell_particle_set))
215 cpassert(
ASSOCIATED(core_particle_set))
216 nshell =
SIZE(shell_particle_set)
217 ALLOCATE (fshell_nonbond(3, nshell))
218 fshell_nonbond = 0.0_dp
219 ALLOCATE (fcore_nonbond(3, nshell))
220 fcore_nonbond = 0.0_dp
222 NULLIFY (shell_particle_set, core_particle_set)
225 IF (fist_nonbond_env%do_nonbonded)
THEN
227 IF (
ASSOCIATED(exclusions))
THEN
228 CALL list_control(atomic_kind_set, particle_set, local_particles, &
229 cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, &
230 core_particle_set, exclusions=exclusions)
232 CALL list_control(atomic_kind_set, particle_set, local_particles, &
233 cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, &
240 particle_set(i)%f(1) = 0.0_dp
241 particle_set(i)%f(2) = 0.0_dp
242 particle_set(i)%f(3) = 0.0_dp
246 shell_particle_set(i)%f(1) = 0.0_dp
247 shell_particle_set(i)%f(2) = 0.0_dp
248 shell_particle_set(i)%f(3) = 0.0_dp
249 core_particle_set(i)%f(1) = 0.0_dp
250 core_particle_set(i)%f(2) = 0.0_dp
251 core_particle_set(i)%f(3) = 0.0_dp
262 pv_urey_bradley = 0.0_dp
265 virial%pv_virial = 0.0_dp
269 pot_manybody = 0.0_dp
275 pot_urey_bradley = 0.0_dp
279 thermo%harm_shell = 0.0_dp
283 WRITE (iw,
'(A)')
" FIST:: FORCES IN INPUT..."
284 WRITE (iw,
'(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1,
SIZE(particle_set))
287 IF (fist_nonbond_env%do_nonbonded)
THEN
292 CALL energy_manybody(fist_nonbond_env, atomic_kind_set, local_particles, particle_set, &
293 cell, pot_manybody, para_env, mm_section)
296 IF (shell_present)
THEN
297 CALL force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
298 pot_nonbond, f_nonbond, pv_nonbond, &
299 fshell_nonbond=fshell_nonbond, fcore_nonbond=fcore_nonbond, &
300 atprop_env=atprop_env, &
301 atomic_kind_set=atomic_kind_set, use_virial=use_virial)
303 CALL force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
304 pot_nonbond, f_nonbond, pv_nonbond, atprop_env=atprop_env, &
305 atomic_kind_set=atomic_kind_set, use_virial=use_virial)
307 use_virial=use_virial)
312 WRITE (iw,
'(A)')
" FIST:: NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
313 WRITE (iw,
'(3f15.9)') f_nonbond
314 IF (shell_present .AND. shell_model_ad)
THEN
315 WRITE (iw,
'(A)')
" FIST:: SHELL NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
316 WRITE (iw,
'(3f15.9)') fshell_nonbond
317 WRITE (iw,
'(A)')
" FIST:: CORE NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
318 WRITE (iw,
'(3f15.9)') fcore_nonbond
325 SELECT CASE (ewald_type)
327 fg_coulomb_size = nlocal_particles
329 fg_coulomb_size = natoms
332 ALLOCATE (fg_coulomb(3, fg_coulomb_size))
334 IF (shell_present)
THEN
335 ALLOCATE (fgshell_coulomb(3, nshell))
336 ALLOCATE (fgcore_coulomb(3, nshell))
337 fgshell_coulomb = 0.0_dp
338 fgcore_coulomb = 0.0_dp
340 IF (shell_present .AND. do_multipoles)
THEN
341 cpabort(
"Multipoles and Core-Shell model not implemented.")
345 IF (.NOT. do_multipoles)
THEN
346 CALL ewald_self(ewald_env, cell, atomic_kind_set, local_particles, &
347 thermo%e_self, thermo%e_neut, fist_nonbond_env%charges)
348 IF (atprop_env%energy)
THEN
350 atprop_env%atener, fist_nonbond_env%charges)
351 atprop_env%atener = atprop_env%atener + thermo%e_neut/
SIZE(atprop_env%atener)
358 IF (
ASSOCIATED(fist_nonbond_env%charges))
THEN
359 cpabort(
"Polarizable force field and array charges not implemented!")
363 ewald_pw, fist_nonbond_env, cell, particle_set, &
364 local_particles, thermo, vg_coulomb, pot_nonbond, f_nonbond, &
365 fg_coulomb, use_virial, pv_g, pv_nonbond, mm_section, do_ipol)
368 SELECT CASE (ewald_type)
371 IF (shell_present)
THEN
373 IF (
ASSOCIATED(fist_nonbond_env%charges))
THEN
374 cpabort(
"Core-Shell and array charges not implemented!")
376 IF (do_multipoles)
THEN
377 cpabort(
"Multipole Ewald and CORE-SHELL not yet implemented within Ewald sum!")
379 cpabort(
"Core-Shell model not yet implemented within Ewald sums.")
382 IF (do_multipoles)
THEN
384 IF (
ASSOCIATED(fist_nonbond_env%charges))
THEN
385 cpabort(
"Multipole Ewald and array charges not implemented!")
388 particle_set, local_particles, vg_coulomb, pot_nonbond, thermo%e_neut, &
389 thermo%e_self, multipoles%task, do_correction_bonded=.true., do_forces=.true., &
390 do_stress=use_virial, do_efield=.false., radii=multipoles%radii, &
391 charges=multipoles%charges, dipoles=multipoles%dipoles, &
392 quadrupoles=multipoles%quadrupoles, forces_local=fg_coulomb, &
393 forces_glob=f_nonbond, pv_local=pv_g, pv_glob=pv_nonbond, iw=iw2, &
394 do_debug=.true., atomic_kind_set=atomic_kind_set, mm_section=mm_section)
396 IF (atprop_env%energy)
THEN
397 ALLOCATE (e_coulomb(fg_coulomb_size))
399 IF (atprop_env%stress)
THEN
400 ALLOCATE (pv_coulomb(3, 3, fg_coulomb_size))
402 CALL ewald_evaluate(ewald_env, ewald_pw, cell, atomic_kind_set, particle_set, &
403 local_particles, fg_coulomb, vg_coulomb, pv_g, use_virial=use_virial, &
404 charges=fist_nonbond_env%charges, e_coulomb=e_coulomb, &
405 pv_coulomb=pv_coulomb)
410 IF (shell_present)
THEN
412 IF (
ASSOCIATED(fist_nonbond_env%charges))
THEN
413 cpabort(
"Core-Shell and array charges not implemented!")
415 IF (do_multipoles)
THEN
416 cpabort(
"Multipole Ewald and CORE-SHELL not yet implemented within a PME scheme!")
418 CALL pme_evaluate(ewald_env, ewald_pw, cell, particle_set, vg_coulomb, fg_coulomb, &
419 pv_g, shell_particle_set=shell_particle_set, core_particle_set=core_particle_set, &
420 fgshell_coulomb=fgshell_coulomb, fgcore_coulomb=fgcore_coulomb, use_virial=use_virial, &
422 CALL para_env%sum(fgshell_coulomb)
423 CALL para_env%sum(fgcore_coulomb)
426 IF (do_multipoles)
THEN
427 cpabort(
"Multipole Ewald not yet implemented within a PME scheme!")
429 CALL pme_evaluate(ewald_env, ewald_pw, cell, particle_set, vg_coulomb, fg_coulomb, &
430 pv_g, use_virial=use_virial, charges=fist_nonbond_env%charges, &
434 CALL para_env%sum(fg_coulomb)
437 IF (shell_present)
THEN
439 IF (
ASSOCIATED(fist_nonbond_env%charges))
THEN
440 cpabort(
"Core-Shell and array charges not implemented!")
442 IF (do_multipoles)
THEN
443 cpabort(
"Multipole Ewald and CORE-SHELL not yet implemented within a SPME scheme!")
445 CALL spme_evaluate(ewald_env, ewald_pw, cell, particle_set, fg_coulomb, vg_coulomb, &
446 pv_g, shell_particle_set=shell_particle_set, core_particle_set=core_particle_set, &
447 fgshell_coulomb=fgshell_coulomb, fgcore_coulomb=fgcore_coulomb, use_virial=use_virial, &
449 CALL para_env%sum(fgshell_coulomb)
450 CALL para_env%sum(fgcore_coulomb)
453 IF (do_multipoles)
THEN
454 cpabort(
"Multipole Ewald not yet implemented within a SPME scheme!")
456 CALL spme_evaluate(ewald_env, ewald_pw, cell, particle_set, fg_coulomb, vg_coulomb, &
457 pv_g, use_virial=use_virial, charges=fist_nonbond_env%charges, &
461 CALL para_env%sum(fg_coulomb)
468 IF (fist_nonbond_env%do_nonbonded)
THEN
470 IF (shell_present)
THEN
472 local_particles, particle_set, ewald_env, thermo%e_bonded, &
473 pv_bc, shell_particle_set=shell_particle_set, &
474 core_particle_set=core_particle_set, atprop_env=atprop_env, cell=cell, &
475 use_virial=use_virial)
477 IF (.NOT. do_multipoles)
THEN
479 atomic_kind_set, local_particles, particle_set, &
480 ewald_env, thermo%e_bonded, pv_bc=pv_bc, atprop_env=atprop_env, cell=cell, &
481 use_virial=use_virial)
486 IF (.NOT. do_multipoles)
THEN
488 CALL ewald_print(iw2, pot_nonbond, vg_coulomb, thermo%e_self, thermo%e_neut, &
496 thermo%e_neut = 0.0_dp
500 IF (
ALLOCATED(fg_coulomb))
THEN
501 WRITE (iw,
'(A)')
" FIST:: NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
502 WRITE (iw,
'(3f15.9)') ((fg_coulomb(j, i), j=1, 3), i=1,
SIZE(fg_coulomb, 2))
504 IF (shell_present .AND. shell_model_ad .AND.
ALLOCATED(fgshell_coulomb))
THEN
505 WRITE (iw,
'(A)')
" FIST:: SHELL NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
506 WRITE (iw,
'(3f15.9)') ((fgshell_coulomb(j, i), j=1, 3), i=1,
SIZE(fgshell_coulomb, 2))
507 WRITE (iw,
'(A)')
" FIST:: CORE NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
508 WRITE (iw,
'(3f15.9)') ((fgcore_coulomb(j, i), j=1, 3), i=1,
SIZE(fgcore_coulomb, 2))
513 IF (efield%apply_field)
THEN
514 ALLOCATE (ef_f(3, natoms))
516 efield, use_virial=use_virial, iunit=iw, charges=fist_nonbond_env%charges)
520 IF (
PRESENT(debug))
THEN
522 local_molecules, particle_set, shell_particle_set, &
523 core_particle_set, pot_bond, pot_bend, pot_urey_bradley, &
524 pot_torsion, pot_imptors, pot_opbend, pot_shell, pv_bond, pv_bend, &
525 pv_urey_bradley, pv_torsion, pv_imptors, pv_opbend, &
526 debug%f_bond, debug%f_bend, debug%f_torsion, debug%f_ub, &
527 debug%f_imptors, debug%f_opbend, cell, use_virial, atprop_env)
531 local_molecules, particle_set, shell_particle_set, &
532 core_particle_set, pot_bond, pot_bend, pot_urey_bradley, &
533 pot_torsion, pot_imptors, pot_opbend, pot_shell, pv_bond, pv_bend, &
534 pv_urey_bradley, pv_torsion, pv_imptors, pv_opbend, &
535 cell=cell, use_virial=use_virial, atprop_env=atprop_env)
539 IF (shell_present .AND. shell_model_ad)
THEN
540 ALLOCATE (fcore_shell_bonded(3, nshell))
541 fcore_shell_bonded = 0.0_dp
543 shell_index = particle_set(i)%shell_index
544 IF (shell_index /= 0)
THEN
545 fcore_shell_bonded(1:3, shell_index) = particle_set(i)%f(1:3)
548 CALL para_env%sum(fcore_shell_bonded)
550 shell_index = particle_set(i)%shell_index
551 IF (shell_index /= 0)
THEN
552 particle_set(i)%f(1:3) = fcore_shell_bonded(1:3, shell_index)
555 DEALLOCATE (fcore_shell_bonded)
562 WRITE (iw,
'(A)')
" FIST energy contributions in kcal/mol:"
563 WRITE (iw,
'(1x,"BOND = ",f13.4,'// &
564 '2x,"ANGLE = ",f13.4,'// &
565 '2x,"UBRAD = ",f13.4)') xdum1, xdum2, xdum3
569 WRITE (iw,
'(1x,"TORSION = ",f13.4,'// &
570 '2x,"IMPTORS = ",f13.4,'// &
571 '2x,"OPBEND = ",f13.4)') xdum1, xdum2, xdum3
573 WRITE (iw,
'(A)')
" FIST:: CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
574 WRITE (iw,
'(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1,
SIZE(particle_set))
575 IF (shell_present .AND. shell_model_ad)
THEN
576 WRITE (iw,
'(A)')
" FIST:: SHELL CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
577 WRITE (iw,
'(3f15.9)') ((shell_particle_set(i)%f(j), j=1, 3), i=1,
SIZE(shell_particle_set))
578 WRITE (iw,
'(A)')
" FIST:: CORE CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
579 WRITE (iw,
'(3f15.9)') ((core_particle_set(i)%f(j), j=1, 3), i=1,
SIZE(core_particle_set))
584 thermo%pot = pot_nonbond + pot_bond + pot_bend + pot_torsion + pot_opbend + &
585 pot_imptors + pot_urey_bradley + pot_manybody + pot_shell
587 CALL para_env%sum(thermo%pot)
589 IF (shell_present)
THEN
590 thermo%harm_shell = pot_shell
591 CALL para_env%sum(thermo%harm_shell)
597 thermo%e_gspace = vg_coulomb
598 thermo%pot = thermo%pot + thermo%e_self + thermo%e_neut
599 thermo%pot = thermo%pot + vg_coulomb + thermo%e_bonded
601 IF (do_ipol /=
do_fist_pol_none) thermo%pot = thermo%pot + thermo%e_induction
610 ALLOCATE (f_total(3, natoms))
613 f_total(1, i) = particle_set(i)%f(1) + f_nonbond(1, i)
614 f_total(2, i) = particle_set(i)%f(2) + f_nonbond(2, i)
615 f_total(3, i) = particle_set(i)%f(3) + f_nonbond(3, i)
617 IF (shell_present)
THEN
618 ALLOCATE (fshell_total(3, nshell))
619 fshell_total = 0.0_dp
620 ALLOCATE (fcore_total(3, nshell))
623 fcore_total(1, i) = core_particle_set(i)%f(1) + fcore_nonbond(1, i)
624 fcore_total(2, i) = core_particle_set(i)%f(2) + fcore_nonbond(2, i)
625 fcore_total(3, i) = core_particle_set(i)%f(3) + fcore_nonbond(3, i)
626 fshell_total(1, i) = shell_particle_set(i)%f(1) + fshell_nonbond(1, i)
627 fshell_total(2, i) = shell_particle_set(i)%f(2) + fshell_nonbond(2, i)
628 fshell_total(3, i) = shell_particle_set(i)%f(3) + fshell_nonbond(3, i)
633 WRITE (iw,
'(A)')
" FIST::(1)INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
634 WRITE (iw,
'(3f15.9)') ((f_total(j, i), j=1, 3), i=1, natoms)
635 IF (shell_present .AND. shell_model_ad)
THEN
636 WRITE (iw,
'(A)')
" FIST::(1)SHELL INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
637 WRITE (iw,
'(3f15.9)') ((fshell_total(j, i), j=1, 3), i=1, nshell)
638 WRITE (iw,
'(A)')
" FIST::(1)CORE INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
639 WRITE (iw,
'(3f15.9)') ((fcore_total(j, i), j=1, 3), i=1, nshell)
646 DO iparticle_kind = 1, nparticle_kind
647 nparticle_local = local_particles%n_el(iparticle_kind)
648 DO iparticle_local = 1, nparticle_local
649 ii = local_particles%list(iparticle_kind)%array(iparticle_local)
651 f_total(1, ii) = f_total(1, ii) + fg_coulomb(1, node)
652 f_total(2, ii) = f_total(2, ii) + fg_coulomb(2, node)
653 f_total(3, ii) = f_total(3, ii) + fg_coulomb(3, node)
654 IF (
PRESENT(debug))
THEN
655 debug%f_g(1, ii) = debug%f_g(1, ii) + fg_coulomb(1, node)
656 debug%f_g(2, ii) = debug%f_g(2, ii) + fg_coulomb(2, node)
657 debug%f_g(3, ii) = debug%f_g(3, ii) + fg_coulomb(3, node)
659 IF (atprop_env%energy)
THEN
660 atprop_env%atener(ii) = atprop_env%atener(ii) + e_coulomb(node)
662 IF (atprop_env%stress)
THEN
663 atprop_env%atstress(:, :, ii) = atprop_env%atstress(:, :, ii) + pv_coulomb(:, :, node)
667 IF (atprop_env%energy)
THEN
668 DEALLOCATE (e_coulomb)
670 IF (atprop_env%stress)
THEN
671 DEALLOCATE (pv_coulomb)
676 WRITE (iw,
'(A)')
" FIST::(2)TOTAL FORCES(1)+ ELECTROSTATIC FORCES"
677 WRITE (iw,
'(3f15.9)') ((f_total(j, i), j=1, 3), i=1, natoms)
678 IF (shell_present .AND. shell_model_ad)
THEN
679 WRITE (iw,
'(A)')
" FIST::(2)SHELL TOTAL FORCES(1)+ ELECTROSTATIC FORCES "
680 WRITE (iw,
'(3f15.9)') ((fshell_total(j, i), j=1, 3), i=1, nshell)
681 WRITE (iw,
'(A)')
" FIST::(2)CORE TOTAL FORCES(1)+ ELECTROSTATIC FORCES"
682 WRITE (iw,
'(3f15.9)') ((fcore_total(j, i), j=1, 3), i=1, nshell)
689 virial%pv_virial = pv_nonbond + pv_bond + pv_bend + pv_torsion + pv_imptors + pv_urey_bradley
690 CALL para_env%sum(virial%pv_virial)
697 virial%pv_virial = pv_nonbond + pv_bond + pv_bend + pv_torsion + pv_imptors + pv_urey_bradley + pv_bc
698 CALL para_env%sum(virial%pv_virial)
700 virial%pv_virial = virial%pv_virial + ident*thermo%e_neut
701 virial%pv_virial = virial%pv_virial + pv_g
706 CALL para_env%sum(f_total)
707 IF (shell_present .AND. shell_model_ad)
THEN
708 CALL para_env%sum(fcore_total)
709 CALL para_env%sum(fshell_total)
713 IF (efield%apply_field)
THEN
714 thermo%pot = thermo%pot + ef_ener
715 f_total(1:3, 1:natoms) = f_total(1:3, 1:natoms) + ef_f(1:3, 1:natoms)
717 virial%pv_virial(1:3, 1:3) = virial%pv_virial(1:3, 1:3) + ef_pv(1:3, 1:3)
723 SELECT CASE (ewald_type)
725 IF (shell_present .AND. shell_model_ad)
THEN
727 shell_index = particle_set(i)%shell_index
728 IF (shell_index == 0)
THEN
729 particle_set(i)%f(1:3) = f_total(1:3, i) + fg_coulomb(1:3, i)
731 atomic_kind => particle_set(i)%atomic_kind
733 fc = shell%mass_core/mass
734 fcore_total(1:3, shell_index) = fcore_total(1:3, shell_index) + particle_set(i)%f(1:3)*fc
735 fs = shell%mass_shell/mass
736 fshell_total(1:3, shell_index) = fshell_total(1:3, shell_index) + particle_set(i)%f(1:3)*fs
741 core_particle_set(i)%f(1) = fcore_total(1, i) + fgcore_coulomb(1, i)
742 core_particle_set(i)%f(2) = fcore_total(2, i) + fgcore_coulomb(2, i)
743 core_particle_set(i)%f(3) = fcore_total(3, i) + fgcore_coulomb(3, i)
744 shell_particle_set(i)%f(1) = fshell_total(1, i) + fgshell_coulomb(1, i)
745 shell_particle_set(i)%f(2) = fshell_total(2, i) + fgshell_coulomb(2, i)
746 shell_particle_set(i)%f(3) = fshell_total(3, i) + fgshell_coulomb(3, i)
749 ELSEIF (shell_present .AND. .NOT. shell_model_ad)
THEN
750 cpabort(
"Non adiabatic shell-model not implemented.")
753 particle_set(i)%f(1) = f_total(1, i) + fg_coulomb(1, i)
754 particle_set(i)%f(2) = f_total(2, i) + fg_coulomb(2, i)
755 particle_set(i)%f(3) = f_total(3, i) + fg_coulomb(3, i)
759 IF (shell_present .AND. shell_model_ad)
THEN
761 shell_index = particle_set(i)%shell_index
762 IF (shell_index == 0)
THEN
763 particle_set(i)%f(1:3) = f_total(1:3, i)
765 atomic_kind => particle_set(i)%atomic_kind
767 fc = shell%mass_core/mass
768 fcore_total(1:3, shell_index) = fcore_total(1:3, shell_index) + particle_set(i)%f(1:3)*fc
769 fs = shell%mass_shell/mass
770 fshell_total(1:3, shell_index) = fshell_total(1:3, shell_index) + particle_set(i)%f(1:3)*fs
774 core_particle_set(i)%f(1) = fcore_total(1, i)
775 core_particle_set(i)%f(2) = fcore_total(2, i)
776 core_particle_set(i)%f(3) = fcore_total(3, i)
777 shell_particle_set(i)%f(1) = fshell_total(1, i)
778 shell_particle_set(i)%f(2) = fshell_total(2, i)
779 shell_particle_set(i)%f(3) = fshell_total(3, i)
781 ELSEIF (shell_present .AND. .NOT. shell_model_ad)
THEN
782 cpabort(
"Non adiabatic shell-model not implemented.")
785 particle_set(i)%f(1) = f_total(1, i)
786 particle_set(i)%f(2) = f_total(2, i)
787 particle_set(i)%f(3) = f_total(3, i)
793 WRITE (iw,
'(A)')
" FIST::(3)TOTAL FORCES - THE END..."
794 WRITE (iw,
'(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, natoms)
795 IF (shell_present .AND. shell_model_ad)
THEN
796 WRITE (iw,
'(A)')
" FIST::(3)SHELL TOTAL FORCES - THE END..."
797 WRITE (iw,
'(3f15.9)') ((shell_particle_set(i)%f(j), j=1, 3), i=1, nshell)
798 WRITE (iw,
'(A)')
" FIST::(3)CORE TOTAL FORCES - THE END..."
799 WRITE (iw,
'(3f15.9)') ((core_particle_set(i)%f(j), j=1, 3), i=1, nshell)
801 WRITE (iw,
'(A,f15.9)')
"Energy after FIST calculation.. exiting now ::", thermo%pot
806 IF (
PRESENT(debug))
THEN
807 CALL para_env%sum(pot_manybody)
808 debug%pot_manybody = pot_manybody
809 CALL para_env%sum(pot_nonbond)
810 debug%pot_nonbond = pot_nonbond
811 CALL para_env%sum(pot_bond)
812 debug%pot_bond = pot_bond
813 CALL para_env%sum(pot_bend)
814 debug%pot_bend = pot_bend
815 CALL para_env%sum(pot_torsion)
816 debug%pot_torsion = pot_torsion
817 CALL para_env%sum(pot_imptors)
818 debug%pot_imptors = pot_imptors
819 CALL para_env%sum(pot_opbend)
820 debug%pot_opbend = pot_opbend
821 CALL para_env%sum(pot_urey_bradley)
822 debug%pot_urey_bradley = pot_urey_bradley
823 CALL para_env%sum(f_nonbond)
824 debug%f_nonbond = f_nonbond
826 CALL para_env%sum(pv_nonbond)
827 debug%pv_nonbond = pv_nonbond
828 CALL para_env%sum(pv_bond)
829 debug%pv_bond = pv_bond
830 CALL para_env%sum(pv_bend)
831 debug%pv_bend = pv_bend
832 CALL para_env%sum(pv_torsion)
833 debug%pv_torsion = pv_torsion
834 CALL para_env%sum(pv_imptors)
835 debug%pv_imptors = pv_imptors
836 CALL para_env%sum(pv_urey_bradley)
837 debug%pv_ub = pv_urey_bradley
839 SELECT CASE (ewald_type)
841 debug%pot_g = vg_coulomb
843 debug%f_g = fg_coulomb
845 debug%pot_g = vg_coulomb
846 IF (use_virial) debug%pv_g = pv_g
850 IF (use_virial) debug%pv_g = 0.0_dp
856 CALL print_fist(fist_env, print_section, atomic_kind_set, particle_set, cell)
859 IF (
ALLOCATED(fg_coulomb))
THEN
860 DEALLOCATE (fg_coulomb)
862 IF (
ALLOCATED(f_total))
THEN
865 DEALLOCATE (f_nonbond)
866 IF (shell_present)
THEN
867 DEALLOCATE (fshell_total)
869 DEALLOCATE (fgshell_coulomb)
871 DEALLOCATE (fshell_nonbond)
875 CALL timestop(handle)