LCOV - code coverage report
Current view: top level - src - nnp_force.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 81.7 % 263 215
Test Date: 2025-07-25 12:55:17 Functions: 72.7 % 11 8

            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              : !> \brief  Methods dealing with Neural Network potentials
      10              : !> \author Christoph Schran (christoph.schran@rub.de)
      11              : !> \date   2020-10-10
      12              : ! **************************************************************************************************
      13              : MODULE nnp_force
      14              : 
      15              :    USE atomic_kind_types,               ONLY: atomic_kind_type
      16              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      17              :                                               cp_logger_type
      18              :    USE cp_output_handling,              ONLY: cp_p_file,&
      19              :                                               cp_print_key_finished_output,&
      20              :                                               cp_print_key_should_output,&
      21              :                                               cp_print_key_unit_nr
      22              :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      23              :                                               cp_subsys_type
      24              :    USE cp_units,                        ONLY: cp_unit_from_cp2k
      25              :    USE distribution_1d_types,           ONLY: distribution_1d_type
      26              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      27              :                                               section_vals_type,&
      28              :                                               section_vals_val_get
      29              :    USE kinds,                           ONLY: default_path_length,&
      30              :                                               default_string_length,&
      31              :                                               dp
      32              :    USE nnp_acsf,                        ONLY: nnp_calc_acsf
      33              :    USE nnp_environment_types,           ONLY: nnp_env_get,&
      34              :                                               nnp_type
      35              :    USE nnp_model,                       ONLY: nnp_gradients,&
      36              :                                               nnp_predict
      37              :    USE particle_types,                  ONLY: particle_type
      38              :    USE periodic_table,                  ONLY: get_ptable_info
      39              :    USE physcon,                         ONLY: angstrom
      40              :    USE virial_types,                    ONLY: virial_type
      41              : #include "./base/base_uses.f90"
      42              : 
      43              :    IMPLICIT NONE
      44              : 
      45              :    PRIVATE
      46              : 
      47              :    LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
      48              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'nnp_force'
      49              : 
      50              :    PUBLIC :: nnp_calc_energy_force
      51              : 
      52              : CONTAINS
      53              : 
      54              : ! **************************************************************************************************
      55              : !> \brief Calculate the energy and force for a given configuration with the NNP
      56              : !> \param nnp ...
      57              : !> \param calc_forces ...
      58              : !> \date   2020-10-10
      59              : !> \author Christoph Schran (christoph.schran@rub.de)
      60              : ! **************************************************************************************************
      61          308 :    SUBROUTINE nnp_calc_energy_force(nnp, calc_forces)
      62              :       TYPE(nnp_type), INTENT(INOUT), POINTER             :: nnp
      63              :       LOGICAL, INTENT(IN)                                :: calc_forces
      64              : 
      65              :       CHARACTER(len=*), PARAMETER :: routineN = 'nnp_calc_energy_force'
      66              : 
      67              :       INTEGER                                            :: handle, i, i_com, ig, ind, istart, j, k, &
      68              :                                                             m, mecalc
      69          308 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: allcalc
      70              :       LOGICAL                                            :: calc_stress
      71          308 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: denergydsym
      72          308 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: dsymdxyz, stress
      73          308 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      74              :       TYPE(cp_logger_type), POINTER                      :: logger
      75              :       TYPE(cp_subsys_type), POINTER                      :: subsys
      76              :       TYPE(distribution_1d_type), POINTER                :: local_particles
      77          308 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      78              :       TYPE(section_vals_type), POINTER                   :: print_section
      79              :       TYPE(virial_type), POINTER                         :: virial
      80              : 
      81          308 :       CALL timeset(routineN, handle)
      82              : 
      83          308 :       NULLIFY (particle_set, logger, local_particles, subsys, &
      84          308 :                atomic_kind_set)
      85          308 :       logger => cp_get_default_logger()
      86              : 
      87          308 :       CPASSERT(ASSOCIATED(nnp))
      88              :       CALL nnp_env_get(nnp_env=nnp, particle_set=particle_set, &
      89              :                        subsys=subsys, local_particles=local_particles, &
      90          308 :                        atomic_kind_set=atomic_kind_set)
      91              : 
      92              :       CALL cp_subsys_get(subsys, &
      93          308 :                          virial=virial)
      94              : 
      95          308 :       calc_stress = virial%pv_availability .AND. (.NOT. virial%pv_numer)
      96           22 :       IF (calc_stress .AND. .NOT. calc_forces) CPABORT('Stress cannot be calculated without forces')
      97              : 
      98              :       ! Initialize energy and gradient:
      99       412874 :       nnp%atomic_energy(:, :) = 0.0_dp
     100       696234 :       IF (calc_forces) nnp%myforce(:, :, :) = 0.0_dp
     101       696234 :       IF (calc_forces) nnp%committee_forces(:, :, :) = 0.0_dp
     102         2596 :       IF (calc_stress) nnp%committee_stress(:, :, :) = 0.0_dp
     103              : 
     104              :       !fill coord array
     105          308 :       ig = 1
     106          924 :       DO i = 1, nnp%n_ele
     107       110880 :          DO j = 1, nnp%num_atoms
     108       110572 :             IF (nnp%ele(i) == particle_set(j)%atomic_kind%element_symbol) THEN
     109       219912 :                DO m = 1, 3
     110       219912 :                   nnp%coord(m, ig) = particle_set(j)%r(m)
     111              :                END DO
     112        54978 :                nnp%atoms(ig) = nnp%ele(i)
     113        54978 :                CALL get_ptable_info(nnp%atoms(ig), number=nnp%nuc_atoms(ig))
     114        54978 :                nnp%ele_ind(ig) = i
     115        54978 :                nnp%sort(ig) = j
     116        54978 :                nnp%sort_inv(j) = ig
     117        54978 :                ig = ig + 1
     118              :             END IF
     119              :          END DO
     120              :       END DO
     121              : 
     122              :       ! parallization:
     123              :       mecalc = nnp%num_atoms/logger%para_env%num_pe + &
     124              :                MIN(MOD(nnp%num_atoms, logger%para_env%num_pe)/ &
     125          308 :                    (logger%para_env%mepos + 1), 1)
     126          924 :       ALLOCATE (allcalc(logger%para_env%num_pe))
     127          924 :       allcalc(:) = 0
     128          308 :       CALL logger%para_env%allgather(mecalc, allcalc)
     129          308 :       istart = 1
     130          462 :       DO i = 2, logger%para_env%mepos + 1
     131          462 :          istart = istart + allcalc(i - 1)
     132              :       END DO
     133              : 
     134              :       ! reset extrapolation status
     135          308 :       nnp%output_expol = .FALSE.
     136              : 
     137              :       ! calc atomic contribution to energy and force
     138        27797 :       DO i = istart, istart + mecalc - 1
     139              : 
     140              :          ! determine index of atom type and offset
     141        27489 :          ind = nnp%ele_ind(i)
     142              : 
     143              :          ! reset input nodes of ele(ind):
     144       797181 :          nnp%arc(ind)%layer(1)%node(:) = 0.0_dp
     145              : 
     146              :          ! compute sym fnct values
     147        27489 :          IF (calc_forces) THEN
     148              :             !reset input grads of ele(ind):
     149       368445 :             nnp%arc(ind)%layer(1)%node_grad(:) = 0.0_dp
     150        50820 :             ALLOCATE (dsymdxyz(3, nnp%arc(ind)%n_nodes(1), nnp%num_atoms))
     151    274955604 :             dsymdxyz(:, :, :) = 0.0_dp
     152        12705 :             IF (calc_stress) THEN
     153         6336 :                ALLOCATE (stress(3, 3, nnp%arc(ind)%n_nodes(1)))
     154       770880 :                stress(:, :, :) = 0.0_dp
     155         2112 :                CALL nnp_calc_acsf(nnp, i, dsymdxyz, stress)
     156              :             ELSE
     157        10593 :                CALL nnp_calc_acsf(nnp, i, dsymdxyz)
     158              :             END IF
     159              :          ELSE
     160        14784 :             CALL nnp_calc_acsf(nnp, i)
     161              :          END IF
     162              : 
     163       232617 :          DO i_com = 1, nnp%n_committee
     164              :             ! predict energy
     165       205128 :             CALL nnp_predict(nnp%arc(ind), nnp, i_com)
     166              :             nnp%atomic_energy(i, i_com) = nnp%arc(ind)%layer(nnp%n_layer)%node(1) + &
     167       205128 :                                           nnp%atom_energies(ind)
     168              : 
     169              :             ! predict forces
     170       232617 :             IF (calc_forces) THEN
     171       260568 :                ALLOCATE (denergydsym(nnp%arc(ind)%n_nodes(1)))
     172      2518824 :                denergydsym(:) = 0.0_dp
     173        86856 :                CALL nnp_gradients(nnp%arc(ind), nnp, i_com, denergydsym)
     174      2518824 :                DO j = 1, nnp%arc(ind)%n_nodes(1)
     175    467972736 :                   DO k = 1, nnp%num_atoms
     176   1864595040 :                      DO m = 1, 3
     177   1862163072 :                         nnp%myforce(m, k, i_com) = nnp%myforce(m, k, i_com) - denergydsym(j)*dsymdxyz(m, j, k)
     178              :                      END DO
     179              :                   END DO
     180      2518824 :                   IF (calc_stress) THEN
     181              :                      nnp%committee_stress(:, :, i_com) = nnp%committee_stress(:, :, i_com) - &
     182      6150144 :                                                          denergydsym(j)*stress(:, :, j)
     183              :                   END IF
     184              :                END DO
     185        86856 :                DEALLOCATE (denergydsym)
     186              :             END IF
     187              :          END DO
     188              : 
     189              :          !deallocate memory
     190        27797 :          IF (calc_forces) THEN
     191        12705 :             DEALLOCATE (dsymdxyz)
     192        12705 :             IF (calc_stress) THEN
     193         2112 :                DEALLOCATE (stress)
     194              :             END IF
     195              :          END IF
     196              : 
     197              :       END DO ! loop over num_atoms
     198              : 
     199              :       ! calculate energy:
     200          308 :       CALL logger%para_env%sum(nnp%atomic_energy(:, :))
     201       412874 :       nnp%committee_energy(:) = SUM(nnp%atomic_energy, 1)
     202         2618 :       nnp%nnp_potential_energy = SUM(nnp%committee_energy)/REAL(nnp%n_committee, dp)
     203              : 
     204          308 :       IF (calc_forces) THEN
     205              :          ! bring myforce to force array
     206        25564 :          DO j = 1, nnp%num_atoms
     207       101794 :             DO k = 1, 3
     208       622776 :                nnp%committee_forces(k, (nnp%sort(j)), :) = nnp%myforce(k, j, :)
     209              :             END DO
     210              :          END DO
     211          154 :          CALL logger%para_env%sum(nnp%committee_forces)
     212       622930 :          nnp%nnp_forces(:, :) = SUM(nnp%committee_forces, DIM=3)/REAL(nnp%n_committee, dp)
     213        25564 :          DO j = 1, nnp%num_atoms
     214       178024 :             particle_set(j)%f(:) = nnp%nnp_forces(:, j)
     215              :          END DO
     216              :       END IF
     217              : 
     218          308 :       IF (calc_stress) THEN
     219           22 :          CALL logger%para_env%sum(nnp%committee_stress)
     220         2134 :          virial%pv_virial = SUM(nnp%committee_stress, DIM=3)/REAL(nnp%n_committee, dp)
     221              :       END IF
     222              : 
     223              :       ! Bias the standard deviation of committee disagreement
     224          308 :       IF (nnp%bias) THEN
     225           44 :          CALL nnp_bias_sigma(nnp, calc_forces)
     226           44 :          nnp%nnp_potential_energy = nnp%nnp_potential_energy + nnp%bias_energy
     227           44 :          IF (calc_forces) THEN
     228         8492 :             DO j = 1, nnp%num_atoms
     229        59180 :                particle_set(j)%f(:) = particle_set(j)%f(:) + nnp%bias_forces(:, j)
     230              :             END DO
     231              :          END IF
     232              :          ! print properties if requested
     233           44 :          print_section => section_vals_get_subs_vals(nnp%nnp_input, "BIAS%PRINT")
     234           44 :          CALL nnp_print_bias(nnp, print_section)
     235              :       END IF
     236              : 
     237              :       ! print properties if requested
     238          308 :       print_section => section_vals_get_subs_vals(nnp%nnp_input, "PRINT")
     239          308 :       CALL nnp_print(nnp, print_section)
     240              : 
     241          308 :       DEALLOCATE (allcalc)
     242              : 
     243          308 :       CALL timestop(handle)
     244              : 
     245          616 :    END SUBROUTINE nnp_calc_energy_force
     246              : 
     247              : ! **************************************************************************************************
     248              : !> \brief Calculate bias potential and force based on standard deviation of committee disagreement
     249              : !> \param nnp ...
     250              : !> \param calc_forces ...
     251              : !> \date   2020-10-10
     252              : !> \author Christoph Schran (christoph.schran@rub.de)
     253              : ! **************************************************************************************************
     254           44 :    SUBROUTINE nnp_bias_sigma(nnp, calc_forces)
     255              :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     256              :       LOGICAL, INTENT(IN)                                :: calc_forces
     257              : 
     258              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'nnp_bias_sigma'
     259              : 
     260              :       INTEGER                                            :: handle, i
     261              :       REAL(KIND=dp)                                      :: avrg, pref, sigma
     262              : 
     263           44 :       CALL timeset(routineN, handle)
     264              : 
     265              :       ! init
     266           44 :       sigma = 0.0_dp
     267           44 :       nnp%bias_energy = 0.0_dp
     268        33836 :       IF (calc_forces) nnp%bias_forces = 0.0_dp
     269              : 
     270              :       ! Subtract reference energy of each committee member, if requested
     271           44 :       IF (nnp%bias_align) THEN
     272              :          ! committee energy not used afterward, therefore overwritten
     273          396 :          nnp%committee_energy(:) = nnp%committee_energy(:) - nnp%bias_e_avrg(:)
     274              :       END IF
     275              : 
     276              :       ! <E> = 1/n sum(E_i)
     277              :       ! sigma = sqrt(1/n sum((E_i - <E>)**2))
     278              :       !       = sqrt(1/n sum(dE_i**2))
     279          396 :       avrg = SUM(nnp%committee_energy)/REAL(nnp%n_committee, dp)
     280          396 :       DO i = 1, nnp%n_committee
     281          396 :          sigma = sigma + (nnp%committee_energy(i) - avrg)**2
     282              :       END DO
     283           44 :       sigma = SQRT(sigma/REAL(nnp%n_committee, dp))
     284           44 :       nnp%bias_sigma = sigma
     285              : 
     286           44 :       IF (sigma > nnp%bias_sigma0) THEN
     287              :          ! E_b = 0.5 * kb * (sigma - sigma_0)**2
     288           44 :          nnp%bias_energy = 0.5_dp*nnp%bias_kb*(sigma - nnp%bias_sigma0)**2
     289              : 
     290           44 :          IF (calc_forces) THEN
     291              :             ! nabla(E_b) = kb*(sigma - sigma_0)*nabla(sigma)
     292              :             ! nabla(sigma) = 1/sigma * 1/n sum(dE_i* nabla(dE_i))
     293              :             ! nabla(dE_i) = nabla(E_i) - nabla(<E>)
     294           44 :             pref = nnp%bias_kb*(1.0_dp - nnp%bias_sigma0/sigma)
     295          396 :             DO i = 1, nnp%n_committee
     296              :                nnp%bias_forces(:, :) = nnp%bias_forces(:, :) + &
     297              :                                        (nnp%committee_energy(i) - avrg)* &
     298       270732 :                                        (nnp%committee_forces(:, :, i) - nnp%nnp_forces(:, :))
     299              :             END DO
     300           44 :             pref = pref/REAL(nnp%n_committee, dp)
     301        33836 :             nnp%bias_forces(:, :) = nnp%bias_forces(:, :)*pref
     302              :          END IF
     303              :       END IF
     304              : 
     305           44 :       CALL timestop(handle)
     306              : 
     307           44 :    END SUBROUTINE nnp_bias_sigma
     308              : 
     309              : ! **************************************************************************************************
     310              : !> \brief Print properties according to the requests in input file
     311              : !> \param nnp ...
     312              : !> \param print_section ...
     313              : !> \date   2020-10-10
     314              : !> \author Christoph Schran (christoph.schran@rub.de)
     315              : ! **************************************************************************************************
     316          616 :    SUBROUTINE nnp_print(nnp, print_section)
     317              :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     318              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: print_section
     319              : 
     320              :       INTEGER                                            :: unit_nr
     321              :       LOGICAL                                            :: explicit, file_is_new
     322              :       TYPE(cp_logger_type), POINTER                      :: logger
     323              :       TYPE(section_vals_type), POINTER                   :: print_key
     324              : 
     325          308 :       NULLIFY (logger, print_key)
     326          308 :       logger => cp_get_default_logger()
     327              : 
     328          308 :       print_key => section_vals_get_subs_vals(print_section, "ENERGIES")
     329          308 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     330              :          unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
     331          308 :                                         middle_name="nnp-energies", is_new_file=file_is_new)
     332          308 :          IF (unit_nr > 0) CALL nnp_print_energies(nnp, unit_nr, file_is_new)
     333          308 :          CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     334              :       END IF
     335              : 
     336          308 :       print_key => section_vals_get_subs_vals(print_section, "FORCES")
     337          308 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     338          286 :          IF (unit_nr > 0) CALL nnp_print_forces(nnp, print_key)
     339              :       END IF
     340              : 
     341          308 :       print_key => section_vals_get_subs_vals(print_section, "FORCES_SIGMA")
     342          308 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     343              :          unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
     344          264 :                                         middle_name="nnp-forces-std")
     345          264 :          IF (unit_nr > 0) CALL nnp_print_force_sigma(nnp, unit_nr)
     346          264 :          CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     347              :       END IF
     348              : 
     349              :       ! Output structures with extrapolation warning on any processor
     350          308 :       CALL logger%para_env%sum(nnp%output_expol)
     351          308 :       IF (nnp%output_expol) THEN
     352           22 :          print_key => section_vals_get_subs_vals(print_section, "EXTRAPOLATION")
     353           22 :          IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     354              :             unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
     355           22 :                                            middle_name="nnp-extrapolation")
     356           22 :             IF (unit_nr > 0) CALL nnp_print_expol(nnp, unit_nr)
     357           22 :             CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     358              :          END IF
     359              :       END IF
     360              : 
     361          308 :       print_key => section_vals_get_subs_vals(print_section, "SUM_FORCE")
     362              : 
     363              :       CALL section_vals_val_get(print_section, "SUM_FORCE%ATOM_LIST", &
     364          308 :                                 explicit=explicit)
     365          308 :       IF (explicit) THEN
     366            0 :          IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     367              :             unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".dat", &
     368            0 :                                            middle_name="nnp-sumforce", is_new_file=file_is_new)
     369            0 :             IF (unit_nr > 0) CALL nnp_print_sumforces(nnp, print_section, unit_nr, file_is_new)
     370            0 :             CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     371              :          END IF
     372              :       END IF
     373              : 
     374          308 :    END SUBROUTINE nnp_print
     375              : 
     376              : ! **************************************************************************************************
     377              : !> \brief Print NNP energies and standard deviation sigma
     378              : !> \param nnp ...
     379              : !> \param unit_nr ...
     380              : !> \param file_is_new ...
     381              : !> \date   2020-10-10
     382              : !> \author Christoph Schran (christoph.schran@rub.de)
     383              : ! **************************************************************************************************
     384          154 :    SUBROUTINE nnp_print_energies(nnp, unit_nr, file_is_new)
     385              :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     386              :       INTEGER, INTENT(IN)                                :: unit_nr
     387              :       LOGICAL, INTENT(IN)                                :: file_is_new
     388              : 
     389              :       CHARACTER(LEN=12)                                  :: fmt_string
     390              :       INTEGER                                            :: i
     391              :       REAL(KIND=dp)                                      :: std
     392              : 
     393          154 :       IF (file_is_new) THEN
     394            6 :          WRITE (unit_nr, "(A1,1X,A20)", ADVANCE='no') "#", "NNP Average [a.u.],"
     395            6 :          WRITE (unit_nr, "(A20)", ADVANCE='no') "NNP sigma [a.u.]"
     396           47 :          DO i = 1, nnp%n_committee
     397           47 :             WRITE (unit_nr, "(A17,I3)", ADVANCE='no') "NNP", i
     398              :          END DO
     399            6 :          WRITE (unit_nr, "(A)") ""
     400              :       END IF
     401              : 
     402          154 :       fmt_string = "(2X,  F20.9)"
     403          154 :       WRITE (UNIT=fmt_string(5:6), FMT="(I2)") nnp%n_committee + 2
     404       206437 :       std = SUM((SUM(nnp%atomic_energy, 1) - nnp%nnp_potential_energy)**2)
     405          154 :       std = std/REAL(nnp%n_committee, dp)
     406          154 :       std = SQRT(std)
     407       206437 :       WRITE (unit_nr, fmt_string) nnp%nnp_potential_energy, std, SUM(nnp%atomic_energy, 1)
     408              : 
     409          154 :    END SUBROUTINE nnp_print_energies
     410              : 
     411              : ! **************************************************************************************************
     412              : !> \brief Print nnp forces
     413              : !> \param nnp ...
     414              : !> \param print_key ...
     415              : !> \date   2020-10-10
     416              : !> \author Christoph Schran (christoph.schran@rub.de)
     417              : ! **************************************************************************************************
     418            0 :    SUBROUTINE nnp_print_forces(nnp, print_key)
     419              :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     420              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: print_key
     421              : 
     422              :       CHARACTER(len=default_path_length)                 :: fmt_string, middle_name
     423              :       INTEGER                                            :: i, j, unit_nr
     424              :       TYPE(cp_logger_type), POINTER                      :: logger
     425              : 
     426            0 :       NULLIFY (logger)
     427            0 :       logger => cp_get_default_logger()
     428              : 
     429              :       ! TODO use write_particle_coordinates from particle_methods.F
     430            0 :       DO i = 1, nnp%n_committee
     431            0 :          WRITE (fmt_string, *) i
     432            0 :          WRITE (middle_name, "(A,A)") "nnp-forces-", ADJUSTL(TRIM(fmt_string))
     433              :          unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
     434            0 :                                         middle_name=TRIM(middle_name))
     435            0 :          IF (unit_nr > 0) THEN
     436            0 :             WRITE (unit_nr, *) nnp%num_atoms
     437            0 :             WRITE (unit_nr, "(A,1X,A,A,F20.9)") "NNP forces [a.u.] of committee member", &
     438            0 :                ADJUSTL(TRIM(fmt_string)), "energy [a.u.]=", nnp%committee_energy(i)
     439              : 
     440            0 :             fmt_string = "(A4,1X,3F20.10)"
     441            0 :             DO j = 1, nnp%num_atoms
     442            0 :                WRITE (unit_nr, fmt_string) nnp%atoms(nnp%sort_inv(j)), nnp%committee_forces(:, j, i)
     443              :             END DO
     444              :          END IF
     445            0 :          CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     446              :       END DO
     447              : 
     448            0 :    END SUBROUTINE nnp_print_forces
     449              : 
     450              : ! **************************************************************************************************
     451              : !> \brief Print standard deviation sigma of NNP forces
     452              : !> \param nnp ...
     453              : !> \param unit_nr ...
     454              : !> \date   2020-10-10
     455              : !> \author Christoph Schran (christoph.schran@rub.de)
     456              : ! **************************************************************************************************
     457          132 :    SUBROUTINE nnp_print_force_sigma(nnp, unit_nr)
     458              :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     459              :       INTEGER, INTENT(IN)                                :: unit_nr
     460              : 
     461              :       INTEGER                                            :: i, j
     462              :       REAL(KIND=dp), DIMENSION(3)                        :: var
     463              : 
     464          132 :       IF (unit_nr > 0) THEN
     465          132 :          WRITE (unit_nr, *) nnp%num_atoms
     466          132 :          WRITE (unit_nr, "(A,1X,A)") "NNP sigma of forces [a.u.]"
     467              : 
     468        25476 :          DO i = 1, nnp%num_atoms
     469        25344 :             var = 0.0_dp
     470       228096 :             DO j = 1, nnp%n_committee
     471       836352 :                var = var + (nnp%committee_forces(:, i, j) - nnp%nnp_forces(:, i))**2
     472              :             END DO
     473       101376 :             var = var/REAL(nnp%n_committee, dp)
     474       101376 :             var = SQRT(var)
     475        25476 :             WRITE (unit_nr, "(A4,1X,3F20.10)") nnp%atoms(nnp%sort_inv(i)), var
     476              :          END DO
     477              :       END IF
     478              : 
     479          132 :    END SUBROUTINE nnp_print_force_sigma
     480              : 
     481              : ! **************************************************************************************************
     482              : !> \brief Print structures with extrapolation warning
     483              : !> \param nnp ...
     484              : !> \param unit_nr ...
     485              : !> \date   2020-10-10
     486              : !> \author Christoph Schran (christoph.schran@rub.de)
     487              : ! **************************************************************************************************
     488           11 :    SUBROUTINE nnp_print_expol(nnp, unit_nr)
     489              :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     490              :       INTEGER, INTENT(IN)                                :: unit_nr
     491              : 
     492              :       CHARACTER(len=default_path_length)                 :: fmt_string
     493              :       INTEGER                                            :: i
     494              :       REAL(KIND=dp)                                      :: mass, unit_conv
     495              :       REAL(KIND=dp), DIMENSION(3)                        :: com
     496              : 
     497           11 :       nnp%expol = nnp%expol + 1
     498           11 :       WRITE (unit_nr, *) nnp%num_atoms
     499           11 :       WRITE (unit_nr, "(A,1X,I6)") "NNP extrapolation point N =", nnp%expol
     500              : 
     501              :       ! move to COM of solute and wrap the box
     502              :       ! coord not needed afterwards, therefore manipulation ok
     503           11 :       com = 0.0_dp
     504           11 :       mass = 0.0_dp
     505           44 :       DO i = 1, nnp%num_atoms
     506           33 :          CALL get_ptable_info(nnp%atoms(i), amass=unit_conv)
     507          132 :          com(:) = com(:) + nnp%coord(:, i)*unit_conv
     508           77 :          mass = mass + unit_conv
     509              :       END DO
     510           44 :       com(:) = com(:)/mass
     511              : 
     512           44 :       DO i = 1, nnp%num_atoms
     513          143 :          nnp%coord(:, i) = nnp%coord(:, i) - com(:)
     514              :       END DO
     515              : 
     516              :       ! write out coordinates
     517           11 :       unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM("angstrom"))
     518           11 :       fmt_string = "(A4,1X,3F20.10)"
     519           44 :       DO i = 1, nnp%num_atoms
     520              :          WRITE (unit_nr, fmt_string) &
     521           33 :             nnp%atoms(nnp%sort_inv(i)), &
     522           33 :             nnp%coord(1, nnp%sort_inv(i))*unit_conv, &
     523           33 :             nnp%coord(2, nnp%sort_inv(i))*unit_conv, &
     524           77 :             nnp%coord(3, nnp%sort_inv(i))*unit_conv
     525              :       END DO
     526              : 
     527           11 :    END SUBROUTINE nnp_print_expol
     528              : 
     529              : ! **************************************************************************************************
     530              : !> \brief Print properties number according the requests in input file
     531              : !> \param nnp ...
     532              : !> \param print_section ...
     533              : !> \date   2020-10-10
     534              : !> \author Christoph Schran (christoph.schran@rub.de)
     535              : ! **************************************************************************************************
     536           44 :    SUBROUTINE nnp_print_bias(nnp, print_section)
     537              :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     538              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: print_section
     539              : 
     540              :       INTEGER                                            :: unit_nr
     541              :       LOGICAL                                            :: file_is_new
     542              :       TYPE(cp_logger_type), POINTER                      :: logger
     543              :       TYPE(section_vals_type), POINTER                   :: print_key
     544              : 
     545           44 :       NULLIFY (logger, print_key)
     546           44 :       logger => cp_get_default_logger()
     547              : 
     548           44 :       print_key => section_vals_get_subs_vals(print_section, "BIAS_ENERGY")
     549           44 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     550              :          unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
     551           44 :                                         middle_name="nnp-bias-energy", is_new_file=file_is_new)
     552           44 :          IF (unit_nr > 0) CALL nnp_print_bias_energy(nnp, unit_nr, file_is_new)
     553           44 :          CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     554              :       END IF
     555              : 
     556           44 :       print_key => section_vals_get_subs_vals(print_section, "BIAS_FORCES")
     557           44 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     558              :          unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
     559            0 :                                         middle_name="nnp-bias-forces")
     560            0 :          IF (unit_nr > 0) CALL nnp_print_bias_forces(nnp, unit_nr)
     561            0 :          CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     562              :       END IF
     563              : 
     564           44 :    END SUBROUTINE nnp_print_bias
     565              : 
     566              : ! **************************************************************************************************
     567              : !> \brief Print NNP bias energies
     568              : !> \param nnp ...
     569              : !> \param unit_nr ...
     570              : !> \param file_is_new ...
     571              : !> \date   2020-10-10
     572              : !> \author Christoph Schran (christoph.schran@rub.de)
     573              : ! **************************************************************************************************
     574           22 :    SUBROUTINE nnp_print_bias_energy(nnp, unit_nr, file_is_new)
     575              :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     576              :       INTEGER, INTENT(IN)                                :: unit_nr
     577              :       LOGICAL, INTENT(IN)                                :: file_is_new
     578              : 
     579              :       CHARACTER(len=default_path_length)                 :: fmt_string
     580              :       INTEGER                                            :: i
     581              : 
     582           22 :       IF (file_is_new) THEN
     583            1 :          WRITE (unit_nr, "(A1)", ADVANCE='no') "#"
     584            1 :          WRITE (unit_nr, "(2(2X,A19))", ADVANCE='no') "Sigma [a.u.]", "Bias energy [a.u.]"
     585            9 :          DO i = 1, nnp%n_committee
     586            9 :             IF (nnp%bias_align) THEN
     587            8 :                WRITE (unit_nr, "(2X,A16,I3)", ADVANCE='no') "shifted E_NNP", i
     588              :             ELSE
     589            0 :                WRITE (unit_nr, "(2X,A16,I3)", ADVANCE='no') "E_NNP", i
     590              :             END IF
     591              :          END DO
     592            1 :          WRITE (unit_nr, "(A)") ""
     593              : 
     594              :       END IF
     595              : 
     596           22 :       WRITE (fmt_string, "(A,I3,A)") "(2X,", nnp%n_committee + 2, "(F20.9,1X))"
     597           22 :       WRITE (unit_nr, fmt_string) nnp%bias_sigma, nnp%bias_energy, nnp%committee_energy
     598              : 
     599           22 :    END SUBROUTINE nnp_print_bias_energy
     600              : 
     601              : ! **************************************************************************************************
     602              : !> \brief Print NNP bias forces
     603              : !> \param nnp ...
     604              : !> \param unit_nr ...
     605              : !> \date   2020-10-10
     606              : !> \author Christoph Schran (christoph.schran@rub.de)
     607              : ! **************************************************************************************************
     608            0 :    SUBROUTINE nnp_print_bias_forces(nnp, unit_nr)
     609              :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     610              :       INTEGER, INTENT(IN)                                :: unit_nr
     611              : 
     612              :       CHARACTER(len=default_path_length)                 :: fmt_string
     613              :       INTEGER                                            :: i
     614              : 
     615              :       ! TODO use write_particle_coordinates from particle_methods.F
     616            0 :       WRITE (unit_nr, *) nnp%num_atoms
     617            0 :       WRITE (unit_nr, "(A,F20.9)") "NNP bias forces [a.u.] for bias energy [a.u]=", nnp%bias_energy
     618              : 
     619            0 :       fmt_string = "(A4,1X,3F20.10)"
     620            0 :       DO i = 1, nnp%num_atoms
     621            0 :          WRITE (unit_nr, fmt_string) nnp%atoms(nnp%sort_inv(i)), nnp%bias_forces(:, i)
     622              :       END DO
     623              : 
     624            0 :    END SUBROUTINE nnp_print_bias_forces
     625              : 
     626              : ! **************************************************************************************************
     627              : !> \brief Print NNP summed forces
     628              : !> \param nnp ...
     629              : !> \param print_section ...
     630              : !> \param unit_nr ...
     631              : !> \param file_is_new ...
     632              : !> \date   2020-10-10
     633              : !> \author Christoph Schran (christoph.schran@rub.de)
     634              : ! **************************************************************************************************
     635            0 :    SUBROUTINE nnp_print_sumforces(nnp, print_section, unit_nr, file_is_new)
     636              :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     637              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: print_section
     638              :       INTEGER, INTENT(IN)                                :: unit_nr
     639              :       LOGICAL, INTENT(IN)                                :: file_is_new
     640              : 
     641              :       CHARACTER(len=default_path_length)                 :: fmt_string
     642              :       CHARACTER(LEN=default_string_length), &
     643            0 :          DIMENSION(:), POINTER                           :: atomlist
     644              :       INTEGER                                            :: i, ig, j, n
     645              :       REAL(KIND=dp), DIMENSION(3)                        :: rvec
     646              : 
     647            0 :       NULLIFY (atomlist)
     648            0 :       IF (file_is_new) THEN
     649            0 :          WRITE (unit_nr, "(A)") "# Summed forces [a.u.]"
     650              :       END IF
     651              : 
     652            0 :       rvec = 0.0_dp
     653              : 
     654              :       ! get atoms to sum over:
     655              :       CALL section_vals_val_get(print_section, "SUM_FORCE%ATOM_LIST", &
     656            0 :                                 c_vals=atomlist)
     657            0 :       IF (ASSOCIATED(atomlist)) THEN
     658            0 :          n = SIZE(atomlist)
     659            0 :          DO i = 1, nnp%num_atoms
     660            0 :             DO j = 1, n
     661            0 :                ig = nnp%sort_inv(i)
     662            0 :                IF (TRIM(ADJUSTL(atomlist(j))) == TRIM(ADJUSTL(nnp%atoms(ig)))) THEN
     663            0 :                   rvec(:) = rvec(:) + nnp%nnp_forces(:, i)
     664              :                END IF
     665              :             END DO
     666              :          END DO
     667              :       END IF
     668              : 
     669            0 :       fmt_string = "(3(F20.10,1X))"
     670            0 :       WRITE (unit_nr, fmt_string) rvec
     671              : 
     672            0 :    END SUBROUTINE nnp_print_sumforces
     673              : 
     674              : END MODULE nnp_force
        

Generated by: LCOV version 2.0-1