LCOV - code coverage report
Current view: top level - src - fist_force.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 84.6 % 428 362
Test Date: 2025-07-25 12:55:17 Functions: 66.7 % 3 2

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \par History
      10              : !>      Torsions added(DG)05-Dec-2000
      11              : !>      Variable names changed(DG)05-Dec-2000
      12              : !>      CJM SEPT-12-2002: int_env is now passed
      13              : !>      CJM NOV-30-2003: only uses fist_env
      14              : !> \author CJM & JGH
      15              : ! **************************************************************************************************
      16              : MODULE fist_force
      17              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      18              :                                               get_atomic_kind
      19              :    USE atprop_types,                    ONLY: atprop_type
      20              :    USE cell_methods,                    ONLY: init_cell
      21              :    USE cell_types,                      ONLY: cell_type
      22              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      23              :                                               cp_logger_type
      24              :    USE cp_output_handling,              ONLY: cp_p_file,&
      25              :                                               cp_print_key_finished_output,&
      26              :                                               cp_print_key_should_output,&
      27              :                                               cp_print_key_unit_nr
      28              :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      29              :                                               cp_subsys_type
      30              :    USE cp_units,                        ONLY: cp_unit_from_cp2k
      31              :    USE distribution_1d_types,           ONLY: distribution_1d_type
      32              :    USE ewald_environment_types,         ONLY: ewald_env_get,&
      33              :                                               ewald_environment_type
      34              :    USE ewald_pw_methods,                ONLY: ewald_pw_grid_update
      35              :    USE ewald_pw_types,                  ONLY: ewald_pw_type
      36              :    USE ewalds,                          ONLY: ewald_evaluate,&
      37              :                                               ewald_print,&
      38              :                                               ewald_self,&
      39              :                                               ewald_self_atom
      40              :    USE ewalds_multipole,                ONLY: ewald_multipole_evaluate
      41              :    USE exclusion_types,                 ONLY: exclusion_type
      42              :    USE fist_efield_methods,             ONLY: fist_dipole,&
      43              :                                               fist_efield_energy_force
      44              :    USE fist_efield_types,               ONLY: fist_efield_type
      45              :    USE fist_energy_types,               ONLY: fist_energy_type
      46              :    USE fist_environment_types,          ONLY: fist_env_get,&
      47              :                                               fist_environment_type
      48              :    USE fist_intra_force,                ONLY: force_intra_control
      49              :    USE fist_neighbor_list_control,      ONLY: list_control
      50              :    USE fist_nonbond_env_types,          ONLY: fist_nonbond_env_type
      51              :    USE fist_nonbond_force,              ONLY: bonded_correct_gaussian,&
      52              :                                               force_nonbond
      53              :    USE fist_pol_scf,                    ONLY: fist_pol_evaluate
      54              :    USE input_constants,                 ONLY: do_fist_pol_none
      55              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      56              :                                               section_vals_type
      57              :    USE kinds,                           ONLY: dp
      58              :    USE manybody_eam,                    ONLY: density_nonbond
      59              :    USE manybody_potential,              ONLY: energy_manybody,&
      60              :                                               force_nonbond_manybody
      61              :    USE message_passing,                 ONLY: mp_para_env_type
      62              :    USE molecule_kind_types,             ONLY: molecule_kind_type
      63              :    USE molecule_types,                  ONLY: molecule_type
      64              :    USE multipole_types,                 ONLY: multipole_type
      65              :    USE particle_types,                  ONLY: particle_type
      66              :    USE pme,                             ONLY: pme_evaluate
      67              :    USE pw_poisson_types,                ONLY: do_ewald_ewald,&
      68              :                                               do_ewald_none,&
      69              :                                               do_ewald_pme,&
      70              :                                               do_ewald_spme
      71              :    USE shell_potential_types,           ONLY: shell_kind_type
      72              :    USE spme,                            ONLY: spme_evaluate
      73              :    USE virial_types,                    ONLY: virial_type
      74              : #include "./base/base_uses.f90"
      75              : 
      76              :    IMPLICIT NONE
      77              :    PRIVATE
      78              :    PUBLIC :: fist_calc_energy_force
      79              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_force'
      80              : 
      81              :    TYPE debug_variables_type
      82              :       REAL(KIND=dp)                           :: pot_manybody = 0.0_dp, pot_bend = 0.0_dp, pot_torsion = 0.0_dp
      83              :       REAL(KIND=dp)                           :: pot_nonbond = 0.0_dp, pot_g = 0.0_dp, pot_bond = 0.0_dp
      84              :       REAL(KIND=dp)                           :: pot_imptors = 0.0_dp, pot_urey_bradley = 0.0_dp
      85              :       REAL(KIND=dp)                           :: pot_opbend = 0.0_dp
      86              :       REAL(KIND=dp), DIMENSION(:, :), POINTER :: f_nonbond => NULL(), f_g => NULL(), f_bond => NULL(), f_bend => NULL(), &
      87              :                                                  f_torsion => NULL(), f_imptors => NULL(), f_ub => NULL(), &
      88              :                                                  f_opbend => NULL()
      89              :       REAL(KIND=dp), DIMENSION(3, 3)          :: pv_nonbond = 0.0_dp, pv_g = 0.0_dp, pv_bond = 0.0_dp, pv_ub = 0.0_dp, &
      90              :                                                  pv_bend = 0.0_dp, pv_torsion = 0.0_dp, pv_imptors = 0.0_dp, &
      91              :                                                  pv_opbend = 0.0_dp
      92              :    END TYPE debug_variables_type
      93              : 
      94              : CONTAINS
      95              : 
      96              : ! **************************************************************************************************
      97              : !> \brief Calculates the total potential energy, total force, and the
      98              : !>      total pressure tensor from the potentials
      99              : !> \param fist_env ...
     100              : !> \param debug ...
     101              : !> \par History
     102              : !>      Harald Forbert(Dec-2000): Changes for multiple linked lists
     103              : !>      cjm, 20-Feb-2001: box_ref used to initialize ewald.  Now
     104              : !>                        have consistent restarts with npt and ewald
     105              : !>      JGH(15-Mar-2001): box_change replaces ensemble parameter
     106              : !>                          Call ewald_setup if box has changed
     107              : !>                          Consistent setup for PME and SPME
     108              : !>      cjm, 28-Feb-2006: box_change is gone
     109              : !> \author CJM & JGH
     110              : ! **************************************************************************************************
     111        79152 :    SUBROUTINE fist_calc_energy_force(fist_env, debug)
     112              :       TYPE(fist_environment_type), POINTER               :: fist_env
     113              :       TYPE(debug_variables_type), INTENT(INOUT), &
     114              :          OPTIONAL                                        :: debug
     115              : 
     116              :       CHARACTER(len=*), PARAMETER :: routineN = 'fist_calc_energy_force'
     117              : 
     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, &
     125              :          xdum2, xdum3
     126        79152 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ef_f, f_nonbond, f_total, fcore_nonbond, &
     127        79152 :          fcore_shell_bonded, fcore_total, fg_coulomb, fgcore_coulomb, fgshell_coulomb, &
     128        79152 :          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        79152 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: e_coulomb
     133        79152 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     134              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     135              :       TYPE(atprop_type), POINTER                         :: atprop_env
     136              :       TYPE(cell_type), POINTER                           :: cell
     137              :       TYPE(cp_logger_type), POINTER                      :: logger
     138              :       TYPE(cp_subsys_type), POINTER                      :: subsys
     139              :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
     140              :       TYPE(ewald_environment_type), POINTER              :: ewald_env
     141              :       TYPE(ewald_pw_type), POINTER                       :: ewald_pw
     142        79152 :       TYPE(exclusion_type), DIMENSION(:), POINTER        :: exclusions
     143              :       TYPE(fist_efield_type), POINTER                    :: efield
     144              :       TYPE(fist_energy_type), POINTER                    :: thermo
     145              :       TYPE(fist_nonbond_env_type), POINTER               :: fist_nonbond_env
     146        79152 :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
     147        79152 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     148              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     149              :       TYPE(multipole_type), POINTER                      :: multipoles
     150        79152 :       TYPE(particle_type), DIMENSION(:), POINTER         :: core_particle_set, particle_set, &
     151        79152 :                                                             shell_particle_set
     152              :       TYPE(section_vals_type), POINTER                   :: force_env_section, mm_section, &
     153              :                                                             print_section
     154              :       TYPE(shell_kind_type), POINTER                     :: shell
     155              :       TYPE(virial_type), POINTER                         :: virial
     156              : 
     157        79152 :       CALL timeset(routineN, handle)
     158        79152 :       NULLIFY (logger)
     159        79152 :       NULLIFY (subsys, virial, atprop_env, para_env, force_env_section)
     160        79152 :       logger => cp_get_default_logger()
     161              : 
     162              :       CALL fist_env_get(fist_env, &
     163              :                         subsys=subsys, &
     164              :                         para_env=para_env, &
     165        79152 :                         input=force_env_section)
     166              : 
     167              :       CALL cp_subsys_get(subsys, &
     168              :                          virial=virial, &
     169        79152 :                          atprop=atprop_env)
     170              : 
     171        79152 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     172        79152 :       NULLIFY (atomic_kind, atomic_kind_set, cell, ewald_pw, ewald_env, &
     173        79152 :                fist_nonbond_env, mm_section, local_molecules, local_particles, &
     174        79152 :                molecule_kind_set, molecule_set, particle_set, print_section, &
     175        79152 :                shell, shell_particle_set, core_particle_set, thermo, multipoles, &
     176        79152 :                e_coulomb)
     177              : 
     178        79152 :       mm_section => section_vals_get_subs_vals(force_env_section, "MM")
     179              :       iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%DERIVATIVES", &
     180        79152 :                                 extension=".mmLog")
     181              :       iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%EWALD_INFO", &
     182        79152 :                                  extension=".mmLog")
     183              : 
     184              :       CALL fist_env_get(fist_env, ewald_pw=ewald_pw, ewald_env=ewald_env, &
     185              :                         local_particles=local_particles, particle_set=particle_set, &
     186              :                         atomic_kind_set=atomic_kind_set, molecule_set=molecule_set, &
     187              :                         local_molecules=local_molecules, thermo=thermo, &
     188              :                         molecule_kind_set=molecule_kind_set, fist_nonbond_env=fist_nonbond_env, &
     189              :                         cell=cell, shell_model=shell_present, shell_model_ad=shell_model_ad, &
     190        79152 :                         multipoles=multipoles, exclusions=exclusions, efield=efield)
     191              : 
     192              :       CALL ewald_env_get(ewald_env, ewald_type=ewald_type, do_multipoles=do_multipoles, &
     193        79152 :                          do_ipol=do_ipol)
     194              : 
     195              :       ! Initialize ewald grids
     196        79152 :       CALL init_cell(cell)
     197        79152 :       CALL ewald_pw_grid_update(ewald_pw, ewald_env, cell%hmat)
     198              : 
     199        79152 :       natoms = SIZE(particle_set)
     200        79152 :       nlocal_particles = 0
     201        79152 :       nparticle_kind = SIZE(atomic_kind_set)
     202       309644 :       DO ikind = 1, nparticle_kind
     203       309644 :          nlocal_particles = nlocal_particles + local_particles%n_el(ikind)
     204              :       END DO
     205              : 
     206       237456 :       ALLOCATE (f_nonbond(3, natoms))
     207     33130492 :       f_nonbond = 0.0_dp
     208              : 
     209        79152 :       nshell = 0
     210        79152 :       IF (shell_present) THEN
     211              :          CALL fist_env_get(fist_env, shell_particle_set=shell_particle_set, &
     212         9468 :                            core_particle_set=core_particle_set)
     213         9468 :          CPASSERT(ASSOCIATED(shell_particle_set))
     214         9468 :          CPASSERT(ASSOCIATED(core_particle_set))
     215         9468 :          nshell = SIZE(shell_particle_set)
     216        28404 :          ALLOCATE (fshell_nonbond(3, nshell))
     217      3062484 :          fshell_nonbond = 0.0_dp
     218        28404 :          ALLOCATE (fcore_nonbond(3, nshell))
     219      3062484 :          fcore_nonbond = 0.0_dp
     220              :       ELSE
     221        69684 :          NULLIFY (shell_particle_set, core_particle_set)
     222              :       END IF
     223              : 
     224        79152 :       IF (fist_nonbond_env%do_nonbonded) THEN
     225              :          ! First check with list_control to update neighbor lists
     226        79000 :          IF (ASSOCIATED(exclusions)) THEN
     227              :             CALL list_control(atomic_kind_set, particle_set, local_particles, &
     228              :                               cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, &
     229        74864 :                               core_particle_set, exclusions=exclusions)
     230              :          ELSE
     231              :             CALL list_control(atomic_kind_set, particle_set, local_particles, &
     232              :                               cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, &
     233         4136 :                               core_particle_set)
     234              :          END IF
     235              :       END IF
     236              : 
     237              :       ! Initialize force, energy and pressure tensor arrays
     238      8341987 :       DO i = 1, natoms
     239      8262835 :          particle_set(i)%f(1) = 0.0_dp
     240      8262835 :          particle_set(i)%f(2) = 0.0_dp
     241      8341987 :          particle_set(i)%f(3) = 0.0_dp
     242              :       END DO
     243        79152 :       IF (nshell > 0) THEN
     244       772722 :          DO i = 1, nshell
     245       763254 :             shell_particle_set(i)%f(1) = 0.0_dp
     246       763254 :             shell_particle_set(i)%f(2) = 0.0_dp
     247       763254 :             shell_particle_set(i)%f(3) = 0.0_dp
     248       763254 :             core_particle_set(i)%f(1) = 0.0_dp
     249       763254 :             core_particle_set(i)%f(2) = 0.0_dp
     250       772722 :             core_particle_set(i)%f(3) = 0.0_dp
     251              :          END DO
     252              :       END IF
     253              : 
     254        79152 :       IF (use_virial) THEN
     255        17288 :          pv_bc = 0.0_dp
     256        17288 :          pv_bond = 0.0_dp
     257        17288 :          pv_bend = 0.0_dp
     258        17288 :          pv_torsion = 0.0_dp
     259        17288 :          pv_imptors = 0.0_dp
     260        17288 :          pv_opbend = 0.0_dp
     261        17288 :          pv_urey_bradley = 0.0_dp
     262        17288 :          pv_nonbond = 0.0_dp
     263        17288 :          pv_g = 0.0_dp
     264       224744 :          virial%pv_virial = 0.0_dp
     265              :       END IF
     266              : 
     267        79152 :       pot_nonbond = 0.0_dp
     268        79152 :       pot_manybody = 0.0_dp
     269        79152 :       pot_bond = 0.0_dp
     270        79152 :       pot_bend = 0.0_dp
     271        79152 :       pot_torsion = 0.0_dp
     272        79152 :       pot_opbend = 0.0_dp
     273        79152 :       pot_imptors = 0.0_dp
     274        79152 :       pot_urey_bradley = 0.0_dp
     275        79152 :       pot_shell = 0.0_dp
     276        79152 :       vg_coulomb = 0.0_dp
     277        79152 :       thermo%pot = 0.0_dp
     278        79152 :       thermo%harm_shell = 0.0_dp
     279              : 
     280              :       ! Get real-space non-bonded forces:
     281        79152 :       IF (iw > 0) THEN
     282          285 :          WRITE (iw, '(A)') " FIST:: FORCES IN INPUT..."
     283        48061 :          WRITE (iw, '(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, SIZE(particle_set))
     284              :       END IF
     285              : 
     286        79152 :       IF (fist_nonbond_env%do_nonbonded) THEN
     287              :          ! Compute density for EAM
     288        79000 :          CALL density_nonbond(fist_nonbond_env, particle_set, cell, para_env)
     289              : 
     290              :          ! Compute embedding function and manybody energy
     291              :          CALL energy_manybody(fist_nonbond_env, atomic_kind_set, local_particles, particle_set, &
     292        79000 :                               cell, pot_manybody, para_env, mm_section, use_virial)
     293              : 
     294              :          ! Nonbond contribution + Manybody Forces
     295        79000 :          IF (shell_present) THEN
     296              :             CALL force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
     297              :                                pot_nonbond, f_nonbond, pv_nonbond, &
     298              :                                fshell_nonbond=fshell_nonbond, fcore_nonbond=fcore_nonbond, &
     299              :                                atprop_env=atprop_env, &
     300         9380 :                                atomic_kind_set=atomic_kind_set, use_virial=use_virial)
     301              :          ELSE
     302              :             CALL force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
     303              :                                pot_nonbond, f_nonbond, pv_nonbond, atprop_env=atprop_env, &
     304        69620 :                                atomic_kind_set=atomic_kind_set, use_virial=use_virial)
     305              :             CALL force_nonbond_manybody(fist_nonbond_env, particle_set, cell, f_nonbond, pv_nonbond, &
     306        69620 :                                         use_virial=use_virial)
     307              :          END IF
     308              :       END IF
     309              : 
     310        79152 :       IF (iw > 0) THEN
     311          285 :          WRITE (iw, '(A)') " FIST:: NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
     312          285 :          WRITE (iw, '(3f15.9)') f_nonbond
     313          285 :          IF (shell_present .AND. shell_model_ad) THEN
     314           44 :             WRITE (iw, '(A)') " FIST:: SHELL NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
     315           44 :             WRITE (iw, '(3f15.9)') fshell_nonbond
     316           44 :             WRITE (iw, '(A)') " FIST:: CORE NONBOND + R-SPACE ELECTROSTATIC FORCES ..."
     317           44 :             WRITE (iw, '(3f15.9)') fcore_nonbond
     318              :          END IF
     319              :       END IF
     320              : 
     321              :       ! Get g-space non-bonded forces:
     322        79152 :       IF (ewald_type /= do_ewald_none) THEN
     323              :          ! Determine size of the forces array
     324              :          SELECT CASE (ewald_type)
     325              :          CASE (do_ewald_ewald)
     326        31474 :             fg_coulomb_size = nlocal_particles
     327              :          CASE DEFAULT
     328        61433 :             fg_coulomb_size = natoms
     329              :          END SELECT
     330              :          ! Allocate and zeroing arrays
     331       183753 :          ALLOCATE (fg_coulomb(3, fg_coulomb_size))
     332     27719913 :          fg_coulomb = 0.0_dp
     333        61433 :          IF (shell_present) THEN
     334        28362 :             ALLOCATE (fgshell_coulomb(3, nshell))
     335        28362 :             ALLOCATE (fgcore_coulomb(3, nshell))
     336      3062390 :             fgshell_coulomb = 0.0_dp
     337      3062390 :             fgcore_coulomb = 0.0_dp
     338              :          END IF
     339        61433 :          IF (shell_present .AND. do_multipoles) THEN
     340            0 :             CPABORT("Multipoles and Core-Shell model not implemented.")
     341              :          END IF
     342              :          ! If not multipole: Compute self-interaction and neutralizing background
     343              :          ! for multipoles is handled separately..
     344        61433 :          IF (.NOT. do_multipoles) THEN
     345              :             CALL ewald_self(ewald_env, cell, atomic_kind_set, local_particles, &
     346        60275 :                             thermo%e_self, thermo%e_neut, fist_nonbond_env%charges)
     347        60275 :             IF (atprop_env%energy) THEN
     348              :                CALL ewald_self_atom(ewald_env, atomic_kind_set, local_particles, &
     349          636 :                                     atprop_env%atener, fist_nonbond_env%charges)
     350         6054 :                atprop_env%atener = atprop_env%atener + thermo%e_neut/SIZE(atprop_env%atener)
     351              :             END IF
     352              :          END IF
     353              : 
     354              :          ! Polarizable force-field
     355        61433 :          IF (do_ipol /= do_fist_pol_none) THEN
     356              :             ! Check if an array of chagres was provided and in case abort due to lack of implementation
     357          104 :             IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
     358            0 :                CPABORT("Polarizable force field and array charges not implemented!")
     359              :             END IF
     360              :             ! Converge the dipoles self-consistently
     361              :             CALL fist_pol_evaluate(atomic_kind_set, multipoles, ewald_env, &
     362              :                                    ewald_pw, fist_nonbond_env, cell, particle_set, &
     363              :                                    local_particles, thermo, vg_coulomb, pot_nonbond, f_nonbond, &
     364          104 :                                    fg_coulomb, use_virial, pv_g, pv_nonbond, mm_section, do_ipol)
     365              :          ELSE
     366              :             ! Non-Polarizable force-field
     367        29855 :             SELECT CASE (ewald_type)
     368              :             CASE (do_ewald_ewald)
     369              :                ! Parallelisation over atoms --> allocate local atoms
     370        29855 :                IF (shell_present) THEN
     371              :                   ! Check if an array of chagres was provided and in case abort due to lack of implementation
     372            0 :                   IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
     373            0 :                      CPABORT("Core-Shell and array charges not implemented!")
     374              :                   END IF
     375            0 :                   IF (do_multipoles) THEN
     376            0 :                      CPABORT("Multipole Ewald and CORE-SHELL not yet implemented within Ewald sum!")
     377              :                   ELSE
     378            0 :                      CPABORT("Core-Shell model not yet implemented within Ewald sums.")
     379              :                   END IF
     380              :                ELSE
     381        29855 :                   IF (do_multipoles) THEN
     382              :                      ! Check if an array of chagres was provided and in case abort due to lack of implementation
     383         1054 :                      IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
     384            0 :                         CPABORT("Multipole Ewald and array charges not implemented!")
     385              :                      END IF
     386              :                      CALL ewald_multipole_evaluate(ewald_env, ewald_pw, fist_nonbond_env, cell, &
     387              :                                                    particle_set, local_particles, vg_coulomb, pot_nonbond, thermo%e_neut, &
     388              :                                                    thermo%e_self, multipoles%task, do_correction_bonded=.TRUE., do_forces=.TRUE., &
     389              :                                                    do_stress=use_virial, do_efield=.FALSE., radii=multipoles%radii, &
     390              :                                                    charges=multipoles%charges, dipoles=multipoles%dipoles, &
     391              :                                                    quadrupoles=multipoles%quadrupoles, forces_local=fg_coulomb, &
     392              :                                                    forces_glob=f_nonbond, pv_local=pv_g, pv_glob=pv_nonbond, iw=iw2, &
     393         1054 :                                                    do_debug=.TRUE., atomic_kind_set=atomic_kind_set, mm_section=mm_section)
     394              :                   ELSE
     395        28801 :                      IF (atprop_env%energy) THEN
     396          505 :                         ALLOCATE (e_coulomb(fg_coulomb_size))
     397              :                      END IF
     398              :                      CALL ewald_evaluate(ewald_env, ewald_pw, cell, atomic_kind_set, particle_set, &
     399              :                                          local_particles, fg_coulomb, vg_coulomb, pv_g, use_virial=use_virial, &
     400        28801 :                                          charges=fist_nonbond_env%charges, e_coulomb=e_coulomb)
     401              :                   END IF
     402              :                END IF
     403              :             CASE (do_ewald_pme)
     404              :                ! Parallelisation over grids --> allocate all atoms
     405         1818 :                IF (shell_present) THEN
     406              :                   ! Check if an array of chagres was provided and in case abort due to lack of implementation
     407            0 :                   IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
     408            0 :                      CPABORT("Core-Shell and array charges not implemented!")
     409              :                   END IF
     410            0 :                   IF (do_multipoles) THEN
     411            0 :                      CPABORT("Multipole Ewald and CORE-SHELL not yet implemented within a PME scheme!")
     412              :                   ELSE
     413              :                      CALL pme_evaluate(ewald_env, ewald_pw, cell, particle_set, vg_coulomb, fg_coulomb, &
     414              :                                        pv_g, shell_particle_set=shell_particle_set, core_particle_set=core_particle_set, &
     415              :                                        fgshell_coulomb=fgshell_coulomb, fgcore_coulomb=fgcore_coulomb, use_virial=use_virial, &
     416            0 :                                        atprop=atprop_env)
     417            0 :                      CALL para_env%sum(fgshell_coulomb)
     418            0 :                      CALL para_env%sum(fgcore_coulomb)
     419              :                   END IF
     420              :                ELSE
     421         1818 :                   IF (do_multipoles) THEN
     422            0 :                      CPABORT("Multipole Ewald not yet implemented within a PME scheme!")
     423              :                   ELSE
     424              :                      CALL pme_evaluate(ewald_env, ewald_pw, cell, particle_set, vg_coulomb, fg_coulomb, &
     425              :                                        pv_g, use_virial=use_virial, charges=fist_nonbond_env%charges, &
     426         1818 :                                        atprop=atprop_env)
     427              :                   END IF
     428              :                END IF
     429         1818 :                CALL para_env%sum(fg_coulomb)
     430              :             CASE (do_ewald_spme)
     431              :                ! Parallelisation over grids --> allocate all atoms
     432        29656 :                IF (shell_present) THEN
     433              :                   ! Check if an array of charges was provided and in case abort due to lack of implementation
     434         9454 :                   IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
     435            0 :                      CPABORT("Core-Shell and array charges not implemented!")
     436              :                   END IF
     437         9454 :                   IF (do_multipoles) THEN
     438            0 :                      CPABORT("Multipole Ewald and CORE-SHELL not yet implemented within a SPME scheme!")
     439              :                   ELSE
     440              :                      CALL spme_evaluate(ewald_env, ewald_pw, cell, particle_set, fg_coulomb, vg_coulomb, &
     441              :                                         pv_g, shell_particle_set=shell_particle_set, core_particle_set=core_particle_set, &
     442              :                                         fgshell_coulomb=fgshell_coulomb, fgcore_coulomb=fgcore_coulomb, use_virial=use_virial, &
     443         9454 :                                         atprop=atprop_env)
     444         9454 :                      CALL para_env%sum(fgshell_coulomb)
     445         9454 :                      CALL para_env%sum(fgcore_coulomb)
     446              :                   END IF
     447              :                ELSE
     448        20202 :                   IF (do_multipoles) THEN
     449            0 :                      CPABORT("Multipole Ewald not yet implemented within a SPME scheme!")
     450              :                   ELSE
     451              :                      CALL spme_evaluate(ewald_env, ewald_pw, cell, particle_set, fg_coulomb, vg_coulomb, &
     452              :                                         pv_g, use_virial=use_virial, charges=fist_nonbond_env%charges, &
     453        20202 :                                         atprop=atprop_env)
     454              :                   END IF
     455              :                END IF
     456        90985 :                CALL para_env%sum(fg_coulomb)
     457              :             END SELECT
     458              :          END IF
     459              : 
     460              :          ! Subtract the interaction between screening charges. This is a
     461              :          ! correction in real-space and depends on the neighborlists. Therefore
     462              :          ! it is only executed if fist_nonbond_env%do_nonbonded is set.
     463        61433 :          IF (fist_nonbond_env%do_nonbonded) THEN
     464              :             ! Correction for core-shell model
     465        61345 :             IF (shell_present) THEN
     466              :                CALL bonded_correct_gaussian(fist_nonbond_env, atomic_kind_set, &
     467              :                                             local_particles, particle_set, ewald_env, thermo%e_bonded, &
     468              :                                             pv_bc, shell_particle_set=shell_particle_set, &
     469              :                                             core_particle_set=core_particle_set, atprop_env=atprop_env, cell=cell, &
     470         9366 :                                             use_virial=use_virial)
     471              :             ELSE
     472        51979 :                IF (.NOT. do_multipoles) THEN
     473              :                   CALL bonded_correct_gaussian(fist_nonbond_env, &
     474              :                                                atomic_kind_set, local_particles, particle_set, &
     475              :                                                ewald_env, thermo%e_bonded, pv_bc=pv_bc, atprop_env=atprop_env, cell=cell, &
     476        50821 :                                                use_virial=use_virial)
     477              :                END IF
     478              :             END IF
     479              :          END IF
     480              : 
     481        61433 :          IF (.NOT. do_multipoles) THEN
     482              :             ! Multipole code has its own printing routine.
     483              :             CALL ewald_print(iw2, pot_nonbond, vg_coulomb, thermo%e_self, thermo%e_neut, &
     484        60275 :                              thermo%e_bonded)
     485              :          END IF
     486              :       ELSE
     487        17719 :          IF (use_virial) THEN
     488         3182 :             pv_g = 0.0_dp
     489         3182 :             pv_bc = 0.0_dp
     490              :          END IF
     491        17719 :          thermo%e_neut = 0.0_dp
     492              :       END IF
     493              : 
     494        79152 :       IF (iw > 0) THEN
     495          285 :          IF (ALLOCATED(fg_coulomb)) THEN
     496          206 :             WRITE (iw, '(A)') " FIST:: NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
     497          206 :             WRITE (iw, '(3f15.9)') ((fg_coulomb(j, i), j=1, 3), i=1, SIZE(fg_coulomb, 2))
     498              :          END IF
     499          285 :          IF (shell_present .AND. shell_model_ad .AND. ALLOCATED(fgshell_coulomb)) THEN
     500           44 :             WRITE (iw, '(A)') " FIST:: SHELL NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
     501           44 :             WRITE (iw, '(3f15.9)') ((fgshell_coulomb(j, i), j=1, 3), i=1, SIZE(fgshell_coulomb, 2))
     502           44 :             WRITE (iw, '(A)') " FIST:: CORE NONBONDED ELECTROSTATIC FORCES IN G-SPACE..."
     503           44 :             WRITE (iw, '(3f15.9)') ((fgcore_coulomb(j, i), j=1, 3), i=1, SIZE(fgcore_coulomb, 2))
     504              :          END IF
     505              :       END IF
     506              : 
     507              :       ! calculate the action of an (external) electric field
     508        79152 :       IF (efield%apply_field) THEN
     509          348 :          ALLOCATE (ef_f(3, natoms))
     510              :          CALL fist_efield_energy_force(ef_ener, ef_f, ef_pv, atomic_kind_set, particle_set, cell, &
     511          116 :                                        efield, use_virial=use_virial, iunit=iw, charges=fist_nonbond_env%charges)
     512              :       END IF
     513              : 
     514              :       ! Get intramolecular forces
     515        79152 :       IF (PRESENT(debug)) THEN
     516              :          CALL force_intra_control(molecule_set, molecule_kind_set, &
     517              :                                   local_molecules, particle_set, shell_particle_set, &
     518              :                                   core_particle_set, pot_bond, pot_bend, pot_urey_bradley, &
     519              :                                   pot_torsion, pot_imptors, pot_opbend, pot_shell, pv_bond, pv_bend, &
     520              :                                   pv_urey_bradley, pv_torsion, pv_imptors, pv_opbend, &
     521              :                                   debug%f_bond, debug%f_bend, debug%f_torsion, debug%f_ub, &
     522            0 :                                   debug%f_imptors, debug%f_opbend, cell, use_virial, atprop_env)
     523              : 
     524              :       ELSE
     525              :          CALL force_intra_control(molecule_set, molecule_kind_set, &
     526              :                                   local_molecules, particle_set, shell_particle_set, &
     527              :                                   core_particle_set, pot_bond, pot_bend, pot_urey_bradley, &
     528              :                                   pot_torsion, pot_imptors, pot_opbend, pot_shell, pv_bond, pv_bend, &
     529              :                                   pv_urey_bradley, pv_torsion, pv_imptors, pv_opbend, &
     530        79152 :                                   cell=cell, use_virial=use_virial, atprop_env=atprop_env)
     531              :       END IF
     532              : 
     533              :       ! Perform global sum of the intra-molecular (bonded) forces for the core-shell atoms
     534        79152 :       IF (shell_present .AND. shell_model_ad) THEN
     535        28404 :          ALLOCATE (fcore_shell_bonded(3, nshell))
     536      3062484 :          fcore_shell_bonded = 0.0_dp
     537       808626 :          DO i = 1, natoms
     538       799158 :             shell_index = particle_set(i)%shell_index
     539       808626 :             IF (shell_index /= 0) THEN
     540      3053016 :                fcore_shell_bonded(1:3, shell_index) = particle_set(i)%f(1:3)
     541              :             END IF
     542              :          END DO
     543         9468 :          CALL para_env%sum(fcore_shell_bonded)
     544       808626 :          DO i = 1, natoms
     545       799158 :             shell_index = particle_set(i)%shell_index
     546       808626 :             IF (shell_index /= 0) THEN
     547      3053016 :                particle_set(i)%f(1:3) = fcore_shell_bonded(1:3, shell_index)
     548              :             END IF
     549              :          END DO
     550         9468 :          DEALLOCATE (fcore_shell_bonded)
     551              :       END IF
     552              : 
     553        79152 :       IF (iw > 0) THEN
     554          285 :          xdum1 = cp_unit_from_cp2k(pot_bond, "kcalmol")
     555          285 :          xdum2 = cp_unit_from_cp2k(pot_bend, "kcalmol")
     556          285 :          xdum3 = cp_unit_from_cp2k(pot_urey_bradley, "kcalmol")
     557          285 :          WRITE (iw, '(A)') " FIST energy contributions in kcal/mol:"
     558              :          WRITE (iw, '(1x,"BOND    = ",f13.4,'// &
     559              :                 '2x,"ANGLE   = ",f13.4,'// &
     560          285 :                 '2x,"UBRAD   = ",f13.4)') xdum1, xdum2, xdum3
     561          285 :          xdum1 = cp_unit_from_cp2k(pot_torsion, "kcalmol")
     562          285 :          xdum2 = cp_unit_from_cp2k(pot_imptors, "kcalmol")
     563          285 :          xdum3 = cp_unit_from_cp2k(pot_opbend, "kcalmol")
     564              :          WRITE (iw, '(1x,"TORSION = ",f13.4,'// &
     565              :                 '2x,"IMPTORS = ",f13.4,'// &
     566          285 :                 '2x,"OPBEND  = ",f13.4)') xdum1, xdum2, xdum3
     567              : 
     568          285 :          WRITE (iw, '(A)') " FIST:: CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
     569        48061 :          WRITE (iw, '(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, SIZE(particle_set))
     570          285 :          IF (shell_present .AND. shell_model_ad) THEN
     571           44 :             WRITE (iw, '(A)') " FIST:: SHELL CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
     572        16940 :             WRITE (iw, '(3f15.9)') ((shell_particle_set(i)%f(j), j=1, 3), i=1, SIZE(shell_particle_set))
     573           44 :             WRITE (iw, '(A)') " FIST:: CORE CORRECTED BONDED ELECTROSTATIC FORCES + INTERNAL FORCES..."
     574        16940 :             WRITE (iw, '(3f15.9)') ((core_particle_set(i)%f(j), j=1, 3), i=1, SIZE(core_particle_set))
     575              :          END IF
     576              :       END IF
     577              : 
     578              :       ! add up all the potential energies
     579              :       thermo%pot = pot_nonbond + pot_bond + pot_bend + pot_torsion + pot_opbend + &
     580        79152 :                    pot_imptors + pot_urey_bradley + pot_manybody + pot_shell
     581              : 
     582        79152 :       CALL para_env%sum(thermo%pot)
     583              : 
     584        79152 :       IF (shell_present) THEN
     585         9468 :          thermo%harm_shell = pot_shell
     586         9468 :          CALL para_env%sum(thermo%harm_shell)
     587              :       END IF
     588              :       ! add g-space contributions if needed
     589        79152 :       IF (ewald_type /= do_ewald_none) THEN
     590              :          ! e_self, e_neut, and ebonded are already summed over all processors
     591              :          ! vg_coulomb is not calculated in parallel
     592        61433 :          thermo%e_gspace = vg_coulomb
     593        61433 :          thermo%pot = thermo%pot + thermo%e_self + thermo%e_neut
     594        61433 :          thermo%pot = thermo%pot + vg_coulomb + thermo%e_bonded
     595              :          ! add the induction energy of the dipoles for polarizable models
     596        61433 :          IF (do_ipol /= do_fist_pol_none) thermo%pot = thermo%pot + thermo%e_induction
     597              :       END IF
     598              : 
     599              :       ! add up all the forces
     600              :       !
     601              :       ! nonbonded forces might be calculated for atoms not on this node
     602              :       ! ewald forces are strictly local -> sum only over pnode
     603              :       ! We first sum the forces in f_nonbond, this allows for a more efficient
     604              :       ! global sum in the parallel code and in the end copy them back to part
     605       237456 :       ALLOCATE (f_total(3, natoms))
     606     33130492 :       f_total = 0.0_dp
     607      8341987 :       DO i = 1, natoms
     608      8262835 :          f_total(1, i) = particle_set(i)%f(1) + f_nonbond(1, i)
     609      8262835 :          f_total(2, i) = particle_set(i)%f(2) + f_nonbond(2, i)
     610      8341987 :          f_total(3, i) = particle_set(i)%f(3) + f_nonbond(3, i)
     611              :       END DO
     612        79152 :       IF (shell_present) THEN
     613        28404 :          ALLOCATE (fshell_total(3, nshell))
     614      3062484 :          fshell_total = 0.0_dp
     615        28404 :          ALLOCATE (fcore_total(3, nshell))
     616      3062484 :          fcore_total = 0.0_dp
     617       772722 :          DO i = 1, nshell
     618       763254 :             fcore_total(1, i) = core_particle_set(i)%f(1) + fcore_nonbond(1, i)
     619       763254 :             fcore_total(2, i) = core_particle_set(i)%f(2) + fcore_nonbond(2, i)
     620       763254 :             fcore_total(3, i) = core_particle_set(i)%f(3) + fcore_nonbond(3, i)
     621       763254 :             fshell_total(1, i) = shell_particle_set(i)%f(1) + fshell_nonbond(1, i)
     622       763254 :             fshell_total(2, i) = shell_particle_set(i)%f(2) + fshell_nonbond(2, i)
     623       772722 :             fshell_total(3, i) = shell_particle_set(i)%f(3) + fshell_nonbond(3, i)
     624              :          END DO
     625              :       END IF
     626              : 
     627        79152 :       IF (iw > 0) THEN
     628          285 :          WRITE (iw, '(A)') " FIST::(1)INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
     629          285 :          WRITE (iw, '(3f15.9)') ((f_total(j, i), j=1, 3), i=1, natoms)
     630          285 :          IF (shell_present .AND. shell_model_ad) THEN
     631           44 :             WRITE (iw, '(A)') " FIST::(1)SHELL INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
     632           44 :             WRITE (iw, '(3f15.9)') ((fshell_total(j, i), j=1, 3), i=1, nshell)
     633           44 :             WRITE (iw, '(A)') " FIST::(1)CORE INTERNAL + ELECTROSTATIC BONDED + NONBONDED"
     634           44 :             WRITE (iw, '(3f15.9)') ((fcore_total(j, i), j=1, 3), i=1, nshell)
     635              :          END IF
     636              :       END IF
     637              : 
     638              :       ! Adding in the reciprocal forces: EWALD is a special case because of distrubted data
     639        79152 :       IF (ewald_type == do_ewald_ewald) THEN
     640              :          node = 0
     641       131819 :          DO iparticle_kind = 1, nparticle_kind
     642       101860 :             nparticle_local = local_particles%n_el(iparticle_kind)
     643       565241 :             DO iparticle_local = 1, nparticle_local
     644       433422 :                ii = local_particles%list(iparticle_kind)%array(iparticle_local)
     645       433422 :                node = node + 1
     646       433422 :                f_total(1, ii) = f_total(1, ii) + fg_coulomb(1, node)
     647       433422 :                f_total(2, ii) = f_total(2, ii) + fg_coulomb(2, node)
     648       433422 :                f_total(3, ii) = f_total(3, ii) + fg_coulomb(3, node)
     649       433422 :                IF (PRESENT(debug)) THEN
     650            0 :                   debug%f_g(1, ii) = debug%f_g(1, ii) + fg_coulomb(1, node)
     651            0 :                   debug%f_g(2, ii) = debug%f_g(2, ii) + fg_coulomb(2, node)
     652            0 :                   debug%f_g(3, ii) = debug%f_g(3, ii) + fg_coulomb(3, node)
     653              :                END IF
     654       535282 :                IF (atprop_env%energy) THEN
     655          303 :                   atprop_env%atener(ii) = atprop_env%atener(ii) + e_coulomb(node)
     656              :                END IF
     657              :             END DO
     658              :          END DO
     659        29959 :          IF (atprop_env%energy) THEN
     660          202 :             DEALLOCATE (e_coulomb)
     661              :          END IF
     662              :       END IF
     663              : 
     664        79152 :       IF (iw > 0) THEN
     665          285 :          WRITE (iw, '(A)') " FIST::(2)TOTAL FORCES(1)+ ELECTROSTATIC FORCES"
     666          285 :          WRITE (iw, '(3f15.9)') ((f_total(j, i), j=1, 3), i=1, natoms)
     667          285 :          IF (shell_present .AND. shell_model_ad) THEN
     668           44 :             WRITE (iw, '(A)') " FIST::(2)SHELL TOTAL FORCES(1)+ ELECTROSTATIC FORCES "
     669           44 :             WRITE (iw, '(3f15.9)') ((fshell_total(j, i), j=1, 3), i=1, nshell)
     670           44 :             WRITE (iw, '(A)') " FIST::(2)CORE TOTAL FORCES(1)+ ELECTROSTATIC FORCES"
     671           44 :             WRITE (iw, '(3f15.9)') ((fcore_total(j, i), j=1, 3), i=1, nshell)
     672              :          END IF
     673              :       END IF
     674              : 
     675        79152 :       IF (use_virial) THEN
     676              :          ! Add up all the pressure tensors
     677        17288 :          IF (ewald_type == do_ewald_none) THEN
     678        41366 :             virial%pv_virial = pv_nonbond + pv_bond + pv_bend + pv_torsion + pv_imptors + pv_urey_bradley
     679         3182 :             CALL para_env%sum(virial%pv_virial)
     680              :          ELSE
     681        14106 :             ident = 0.0_dp
     682        56424 :             DO i = 1, 3
     683        56424 :                ident(i, i) = 1.0_dp
     684              :             END DO
     685              : 
     686       183378 :             virial%pv_virial = pv_nonbond + pv_bond + pv_bend + pv_torsion + pv_imptors + pv_urey_bradley + pv_bc
     687        14106 :             CALL para_env%sum(virial%pv_virial)
     688              : 
     689       183378 :             virial%pv_virial = virial%pv_virial + ident*thermo%e_neut
     690       183378 :             virial%pv_virial = virial%pv_virial + pv_g
     691              :          END IF
     692              :       END IF
     693              : 
     694              :       ! Sum total forces
     695        79152 :       CALL para_env%sum(f_total)
     696        79152 :       IF (shell_present .AND. shell_model_ad) THEN
     697         9468 :          CALL para_env%sum(fcore_total)
     698         9468 :          CALL para_env%sum(fshell_total)
     699              :       END IF
     700              : 
     701              :       ! contributions from fields (currently all quantities are fully replicated!)
     702        79152 :       IF (efield%apply_field) THEN
     703          116 :          thermo%pot = thermo%pot + ef_ener
     704         1508 :          f_total(1:3, 1:natoms) = f_total(1:3, 1:natoms) + ef_f(1:3, 1:natoms)
     705          116 :          IF (use_virial) THEN
     706           26 :             virial%pv_virial(1:3, 1:3) = virial%pv_virial(1:3, 1:3) + ef_pv(1:3, 1:3)
     707              :          END IF
     708          116 :          DEALLOCATE (ef_f)
     709              :       END IF
     710              : 
     711              :       ! Assign to particle_set
     712        31474 :       SELECT CASE (ewald_type)
     713              :       CASE (do_ewald_spme, do_ewald_pme)
     714        31474 :          IF (shell_present .AND. shell_model_ad) THEN
     715       808584 :             DO i = 1, natoms
     716       799130 :                shell_index = particle_set(i)%shell_index
     717       808584 :                IF (shell_index == 0) THEN
     718       143584 :                   particle_set(i)%f(1:3) = f_total(1:3, i) + fg_coulomb(1:3, i)
     719              :                ELSE
     720       763234 :                   atomic_kind => particle_set(i)%atomic_kind
     721       763234 :                   CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell, mass=mass)
     722       763234 :                   fc = shell%mass_core/mass
     723      3052936 :                   fcore_total(1:3, shell_index) = fcore_total(1:3, shell_index) + particle_set(i)%f(1:3)*fc
     724       763234 :                   fs = shell%mass_shell/mass
     725      3052936 :                   fshell_total(1:3, shell_index) = fshell_total(1:3, shell_index) + particle_set(i)%f(1:3)*fs
     726              :                END IF
     727              :             END DO
     728              : 
     729       772688 :             DO i = 1, nshell
     730       763234 :                core_particle_set(i)%f(1) = fcore_total(1, i) + fgcore_coulomb(1, i)
     731       763234 :                core_particle_set(i)%f(2) = fcore_total(2, i) + fgcore_coulomb(2, i)
     732       763234 :                core_particle_set(i)%f(3) = fcore_total(3, i) + fgcore_coulomb(3, i)
     733       763234 :                shell_particle_set(i)%f(1) = fshell_total(1, i) + fgshell_coulomb(1, i)
     734       763234 :                shell_particle_set(i)%f(2) = fshell_total(2, i) + fgshell_coulomb(2, i)
     735       772688 :                shell_particle_set(i)%f(3) = fshell_total(3, i) + fgshell_coulomb(3, i)
     736              :             END DO
     737              : 
     738        22020 :          ELSEIF (shell_present .AND. .NOT. shell_model_ad) THEN
     739            0 :             CPABORT("Non adiabatic shell-model not implemented.")
     740              :          ELSE
     741      5704088 :             DO i = 1, natoms
     742      5682068 :                particle_set(i)%f(1) = f_total(1, i) + fg_coulomb(1, i)
     743      5682068 :                particle_set(i)%f(2) = f_total(2, i) + fg_coulomb(2, i)
     744      5704088 :                particle_set(i)%f(3) = f_total(3, i) + fg_coulomb(3, i)
     745              :             END DO
     746              :          END IF
     747              :       CASE (do_ewald_ewald, do_ewald_none)
     748        79152 :          IF (shell_present .AND. shell_model_ad) THEN
     749           42 :             DO i = 1, natoms
     750           28 :                shell_index = particle_set(i)%shell_index
     751           42 :                IF (shell_index == 0) THEN
     752           32 :                   particle_set(i)%f(1:3) = f_total(1:3, i)
     753              :                ELSE
     754           20 :                   atomic_kind => particle_set(i)%atomic_kind
     755           20 :                   CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell, mass=mass)
     756           20 :                   fc = shell%mass_core/mass
     757           80 :                   fcore_total(1:3, shell_index) = fcore_total(1:3, shell_index) + particle_set(i)%f(1:3)*fc
     758           20 :                   fs = shell%mass_shell/mass
     759           80 :                   fshell_total(1:3, shell_index) = fshell_total(1:3, shell_index) + particle_set(i)%f(1:3)*fs
     760              :                END IF
     761              :             END DO
     762           34 :             DO i = 1, nshell
     763           20 :                core_particle_set(i)%f(1) = fcore_total(1, i)
     764           20 :                core_particle_set(i)%f(2) = fcore_total(2, i)
     765           20 :                core_particle_set(i)%f(3) = fcore_total(3, i)
     766           20 :                shell_particle_set(i)%f(1) = fshell_total(1, i)
     767           20 :                shell_particle_set(i)%f(2) = fshell_total(2, i)
     768           34 :                shell_particle_set(i)%f(3) = fshell_total(3, i)
     769              :             END DO
     770        47664 :          ELSEIF (shell_present .AND. .NOT. shell_model_ad) THEN
     771            0 :             CPABORT("Non adiabatic shell-model not implemented.")
     772              :          ELSE
     773      1829273 :             DO i = 1, natoms
     774      1781609 :                particle_set(i)%f(1) = f_total(1, i)
     775      1781609 :                particle_set(i)%f(2) = f_total(2, i)
     776      1829273 :                particle_set(i)%f(3) = f_total(3, i)
     777              :             END DO
     778              :          END IF
     779              :       END SELECT
     780              : 
     781        79152 :       IF (iw > 0) THEN
     782          285 :          WRITE (iw, '(A)') " FIST::(3)TOTAL FORCES - THE END..."
     783        48061 :          WRITE (iw, '(3f15.9)') ((particle_set(i)%f(j), j=1, 3), i=1, natoms)
     784          285 :          IF (shell_present .AND. shell_model_ad) THEN
     785           44 :             WRITE (iw, '(A)') " FIST::(3)SHELL TOTAL FORCES - THE END..."
     786        16940 :             WRITE (iw, '(3f15.9)') ((shell_particle_set(i)%f(j), j=1, 3), i=1, nshell)
     787           44 :             WRITE (iw, '(A)') " FIST::(3)CORE TOTAL FORCES - THE END..."
     788        16940 :             WRITE (iw, '(3f15.9)') ((core_particle_set(i)%f(j), j=1, 3), i=1, nshell)
     789              :          END IF
     790          285 :          WRITE (iw, '(A,f15.9)') "Energy after FIST calculation.. exiting now ::", thermo%pot
     791              :       END IF
     792              :       !
     793              :       ! If we are doing debugging, check if variables are present and assign
     794              :       !
     795        79152 :       IF (PRESENT(debug)) THEN
     796            0 :          CALL para_env%sum(pot_manybody)
     797            0 :          debug%pot_manybody = pot_manybody
     798            0 :          CALL para_env%sum(pot_nonbond)
     799            0 :          debug%pot_nonbond = pot_nonbond
     800            0 :          CALL para_env%sum(pot_bond)
     801            0 :          debug%pot_bond = pot_bond
     802            0 :          CALL para_env%sum(pot_bend)
     803            0 :          debug%pot_bend = pot_bend
     804            0 :          CALL para_env%sum(pot_torsion)
     805            0 :          debug%pot_torsion = pot_torsion
     806            0 :          CALL para_env%sum(pot_imptors)
     807            0 :          debug%pot_imptors = pot_imptors
     808            0 :          CALL para_env%sum(pot_opbend)
     809            0 :          debug%pot_opbend = pot_opbend
     810            0 :          CALL para_env%sum(pot_urey_bradley)
     811            0 :          debug%pot_urey_bradley = pot_urey_bradley
     812            0 :          CALL para_env%sum(f_nonbond)
     813            0 :          debug%f_nonbond = f_nonbond
     814            0 :          IF (use_virial) THEN
     815            0 :             CALL para_env%sum(pv_nonbond)
     816            0 :             debug%pv_nonbond = pv_nonbond
     817            0 :             CALL para_env%sum(pv_bond)
     818            0 :             debug%pv_bond = pv_bond
     819            0 :             CALL para_env%sum(pv_bend)
     820            0 :             debug%pv_bend = pv_bend
     821            0 :             CALL para_env%sum(pv_torsion)
     822            0 :             debug%pv_torsion = pv_torsion
     823            0 :             CALL para_env%sum(pv_imptors)
     824            0 :             debug%pv_imptors = pv_imptors
     825            0 :             CALL para_env%sum(pv_urey_bradley)
     826            0 :             debug%pv_ub = pv_urey_bradley
     827              :          END IF
     828            0 :          SELECT CASE (ewald_type)
     829              :          CASE (do_ewald_spme, do_ewald_pme)
     830            0 :             debug%pot_g = vg_coulomb
     831            0 :             debug%pv_g = pv_g
     832            0 :             debug%f_g = fg_coulomb
     833              :          CASE (do_ewald_ewald)
     834            0 :             debug%pot_g = vg_coulomb
     835            0 :             IF (use_virial) debug%pv_g = pv_g
     836              :          CASE default
     837            0 :             debug%pot_g = 0.0_dp
     838            0 :             debug%f_g = 0.0_dp
     839            0 :             IF (use_virial) debug%pv_g = 0.0_dp
     840              :          END SELECT
     841              :       END IF
     842              : 
     843              :       ! print properties if requested
     844        79152 :       print_section => section_vals_get_subs_vals(mm_section, "PRINT")
     845        79152 :       CALL print_fist(fist_env, print_section, atomic_kind_set, particle_set, cell)
     846              : 
     847              :       ! deallocating all local variables
     848        79152 :       IF (ALLOCATED(fg_coulomb)) THEN
     849        61433 :          DEALLOCATE (fg_coulomb)
     850              :       END IF
     851        79152 :       IF (ALLOCATED(f_total)) THEN
     852        79152 :          DEALLOCATE (f_total)
     853              :       END IF
     854        79152 :       DEALLOCATE (f_nonbond)
     855        79152 :       IF (shell_present) THEN
     856         9468 :          DEALLOCATE (fshell_total)
     857         9468 :          IF (ewald_type /= do_ewald_none) THEN
     858         9454 :             DEALLOCATE (fgshell_coulomb)
     859              :          END IF
     860         9468 :          DEALLOCATE (fshell_nonbond)
     861              :       END IF
     862        79152 :       CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%DERIVATIVES")
     863        79152 :       CALL cp_print_key_finished_output(iw2, logger, mm_section, "PRINT%EWALD_INFO")
     864        79152 :       CALL timestop(handle)
     865              : 
     866       237456 :    END SUBROUTINE fist_calc_energy_force
     867              : 
     868              : ! **************************************************************************************************
     869              : !> \brief Print properties number according the requests in input file
     870              : !> \param fist_env ...
     871              : !> \param print_section ...
     872              : !> \param atomic_kind_set ...
     873              : !> \param particle_set ...
     874              : !> \param cell ...
     875              : !> \par History
     876              : !>      [01.2006] created
     877              : !> \author Teodoro Laino
     878              : ! **************************************************************************************************
     879        79152 :    SUBROUTINE print_fist(fist_env, print_section, atomic_kind_set, particle_set, cell)
     880              :       TYPE(fist_environment_type), POINTER               :: fist_env
     881              :       TYPE(section_vals_type), POINTER                   :: print_section
     882              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     883              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     884              :       TYPE(cell_type), POINTER                           :: cell
     885              : 
     886              :       INTEGER                                            :: unit_nr
     887              :       TYPE(cp_logger_type), POINTER                      :: logger
     888              :       TYPE(fist_nonbond_env_type), POINTER               :: fist_nonbond_env
     889              :       TYPE(section_vals_type), POINTER                   :: print_key
     890              : 
     891        79152 :       NULLIFY (logger, print_key, fist_nonbond_env)
     892        79152 :       logger => cp_get_default_logger()
     893        79152 :       print_key => section_vals_get_subs_vals(print_section, "dipole")
     894        79152 :       CALL fist_env_get(fist_env, fist_nonbond_env=fist_nonbond_env)
     895        79152 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), &
     896              :                 cp_p_file)) THEN
     897              :          unit_nr = cp_print_key_unit_nr(logger, print_section, "dipole", &
     898        21204 :                                         extension=".data", middle_name="MM_DIPOLE", log_filename=.FALSE.)
     899              :          CALL fist_dipole(fist_env, print_section, atomic_kind_set, particle_set, &
     900        21204 :                           cell, unit_nr, fist_nonbond_env%charges)
     901        21204 :          CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     902              :       END IF
     903              : 
     904        79152 :    END SUBROUTINE print_fist
     905              : 
     906            0 : END MODULE fist_force
        

Generated by: LCOV version 2.0-1