LCOV - code coverage report
Current view: top level - src/motion - helium_nnp.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 63 128 49.2 %
Date: 2024-04-20 06:29:22 Functions: 2 5 40.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 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 interaction potential
      10             : !> \author Laura Duran
      11             : !> \date   2023-02-17
      12             : ! **************************************************************************************************
      13             : MODULE helium_nnp
      14             : 
      15             :    USE bibliography,                    ONLY: Behler2007,&
      16             :                                               Behler2011,&
      17             :                                               Schran2020a,&
      18             :                                               Schran2020b,&
      19             :                                               cite_reference
      20             :    USE cell_methods,                    ONLY: cell_create
      21             :    USE cell_types,                      ONLY: cell_release,&
      22             :                                               cell_type,&
      23             :                                               pbc
      24             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      25             :                                               cp_logger_type
      26             :    USE cp_output_handling,              ONLY: cp_p_file,&
      27             :                                               cp_print_key_finished_output,&
      28             :                                               cp_print_key_should_output,&
      29             :                                               cp_print_key_unit_nr
      30             :    USE cp_units,                        ONLY: cp_unit_from_cp2k
      31             :    USE helium_types,                    ONLY: helium_solvent_type
      32             :    USE input_section_types,             ONLY: section_vals_get,&
      33             :                                               section_vals_get_subs_vals,&
      34             :                                               section_vals_type,&
      35             :                                               section_vals_val_get
      36             :    USE kinds,                           ONLY: default_path_length,&
      37             :                                               default_string_length,&
      38             :                                               dp
      39             :    USE nnp_environment,                 ONLY: nnp_init_model
      40             :    USE nnp_environment_types,           ONLY: nnp_env_get,&
      41             :                                               nnp_env_set,&
      42             :                                               nnp_type
      43             :    USE periodic_table,                  ONLY: get_ptable_info
      44             :    USE physcon,                         ONLY: angstrom
      45             : #include "../base/base_uses.f90"
      46             : 
      47             :    IMPLICIT NONE
      48             : 
      49             :    PRIVATE
      50             : 
      51             :    LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
      52             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_nnp'
      53             : 
      54             :    PUBLIC :: helium_init_nnp, &
      55             :              helium_nnp_print
      56             : 
      57             : CONTAINS
      58             : 
      59             : ! ***************************************************************************
      60             : !> \brief  Read and initialize all the information for neural network potentials
      61             : !> \param helium ...
      62             : !> \param nnp ...
      63             : !> \param input ...
      64             : !> \date   2023-02-21
      65             : !> \author lduran
      66             : ! **************************************************************************************************
      67           1 :    SUBROUTINE helium_init_nnp(helium, nnp, input)
      68             :       TYPE(helium_solvent_type), INTENT(INOUT)           :: helium
      69             :       TYPE(nnp_type), POINTER                            :: nnp
      70             :       TYPE(section_vals_type), POINTER                   :: input
      71             : 
      72             :       CHARACTER(len=default_path_length)                 :: msg_str
      73             :       CHARACTER(len=default_string_length)               :: elem
      74             :       INTEGER                                            :: i, ig, is, j
      75             :       INTEGER, DIMENSION(3)                              :: periodicity
      76             :       LOGICAL                                            :: found
      77             :       TYPE(cell_type), POINTER                           :: he_cell
      78             :       TYPE(cp_logger_type), POINTER                      :: logger
      79             :       TYPE(section_vals_type), POINTER                   :: sr_cutoff_section
      80             : 
      81           1 :       CALL cite_reference(Behler2007)
      82           1 :       CALL cite_reference(Behler2011)
      83           1 :       CALL cite_reference(Schran2020a)
      84           1 :       CALL cite_reference(Schran2020b)
      85             : 
      86           1 :       NULLIFY (logger)
      87           1 :       logger => cp_get_default_logger()
      88             : 
      89           1 :       CALL nnp_env_set(nnp_env=nnp, nnp_input=input)
      90             : 
      91           1 :       nnp%num_atoms = helium%solute_atoms + 1
      92             : 
      93           1 :       CALL nnp_init_model(nnp, "HELIUM NNP")
      94             : 
      95           1 :       periodicity = 0
      96           1 :       IF (helium%periodic) periodicity = 1
      97           1 :       NULLIFY (he_cell)
      98             :       CALL cell_create(he_cell, hmat=helium%cell_m, &
      99           1 :                        periodic=periodicity, tag="HELIUM NNP")
     100           1 :       CALL nnp_env_set(nnp, cell=he_cell)
     101           1 :       CALL cell_release(he_cell)
     102             : 
     103             :       ! Set up arrays for calculation:
     104           3 :       ALLOCATE (nnp%ele_ind(nnp%num_atoms))
     105           3 :       ALLOCATE (nnp%nuc_atoms(nnp%num_atoms))
     106           3 :       ALLOCATE (nnp%coord(3, nnp%num_atoms))
     107           3 :       ALLOCATE (nnp%nnp_forces(3, nnp%num_atoms))
     108           3 :       ALLOCATE (nnp%atoms(nnp%num_atoms))
     109           3 :       ALLOCATE (nnp%sort(nnp%num_atoms - 1))
     110             : 
     111             :       !fill arrays, assume that order will not change during simulation
     112           1 :       ig = 1
     113           1 :       is = 1
     114           4 :       DO i = 1, nnp%n_ele
     115           3 :          IF (nnp%ele(i) == 'He') THEN
     116           1 :             nnp%atoms(ig) = 'He'
     117           1 :             CALL get_ptable_info(nnp%atoms(ig), number=nnp%nuc_atoms(ig))
     118           1 :             nnp%ele_ind(ig) = i
     119           1 :             ig = ig + 1
     120             :          END IF
     121          13 :          DO j = 1, helium%solute_atoms
     122          12 :             IF (nnp%ele(i) == helium%solute_element(j)) THEN
     123           3 :                nnp%atoms(ig) = nnp%ele(i)
     124           3 :                CALL get_ptable_info(nnp%atoms(ig), number=nnp%nuc_atoms(ig))
     125           3 :                nnp%ele_ind(ig) = i
     126           3 :                nnp%sort(is) = j
     127           3 :                ig = ig + 1
     128           3 :                is = is + 1
     129             :             END IF
     130             :          END DO
     131             :       END DO
     132             : 
     133           3 :       ALLOCATE (helium%nnp_sr_cut(nnp%n_ele))
     134           4 :       helium%nnp_sr_cut = 0.0_dp
     135             : 
     136           1 :       sr_cutoff_section => section_vals_get_subs_vals(nnp%nnp_input, "SR_CUTOFF")
     137           1 :       CALL section_vals_get(sr_cutoff_section, n_repetition=is)
     138           3 :       DO i = 1, is
     139           2 :          CALL section_vals_val_get(sr_cutoff_section, "ELEMENT", c_val=elem, i_rep_section=i)
     140           2 :          found = .FALSE.
     141           8 :          DO ig = 1, nnp%n_ele
     142           8 :             IF (TRIM(nnp%ele(ig)) == TRIM(elem)) THEN
     143           2 :                found = .TRUE.
     144             :                CALL section_vals_val_get(sr_cutoff_section, "RADIUS", r_val=helium%nnp_sr_cut(ig), &
     145           2 :                                          i_rep_section=i)
     146             :             END IF
     147             :          END DO
     148           3 :          IF (.NOT. found) THEN
     149           0 :             msg_str = "SR_CUTOFF for element "//TRIM(elem)//" defined but not found in NNP"
     150           0 :             CPWARN(msg_str)
     151             :          END IF
     152             :       END DO
     153           4 :       helium%nnp_sr_cut(:) = helium%nnp_sr_cut(:)**2
     154             : 
     155           1 :       RETURN
     156             : 
     157           1 :    END SUBROUTINE helium_init_nnp
     158             : 
     159             : ! **************************************************************************************************
     160             : !> \brief Print properties according to the requests in input file
     161             : !> \param nnp ...
     162             : !> \param print_section ...
     163             : !> \param ind_he ...
     164             : !> \date   2023-07-31
     165             : !> \author Laura Duran
     166             : ! **************************************************************************************************
     167       55172 :    SUBROUTINE helium_nnp_print(nnp, print_section, ind_he)
     168             :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     169             :       TYPE(section_vals_type), INTENT(IN), POINTER       :: print_section
     170             :       INTEGER, INTENT(IN)                                :: ind_he
     171             : 
     172             :       INTEGER                                            :: unit_nr
     173             :       LOGICAL                                            :: file_is_new
     174             :       TYPE(cp_logger_type), POINTER                      :: logger
     175             :       TYPE(section_vals_type), POINTER                   :: print_key
     176             : 
     177       55172 :       NULLIFY (logger, print_key)
     178       55172 :       logger => cp_get_default_logger()
     179             : 
     180       55172 :       print_key => section_vals_get_subs_vals(print_section, "ENERGIES")
     181       55172 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     182             :          unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
     183           0 :                                         middle_name="helium-nnp-energies", is_new_file=file_is_new)
     184           0 :          IF (unit_nr > 0) CALL helium_nnp_print_energies(nnp, unit_nr, file_is_new)
     185           0 :          CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     186             :       END IF
     187             : 
     188       55172 :       print_key => section_vals_get_subs_vals(print_section, "FORCES_SIGMA")
     189       55172 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     190             :          unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
     191           0 :                                         middle_name="helium-nnp-forces-std", is_new_file=file_is_new)
     192           0 :          IF (unit_nr > 0) CALL helium_nnp_print_force_sigma(nnp, unit_nr, file_is_new)
     193           0 :          CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     194             :       END IF
     195             : 
     196       55172 :       CALL logger%para_env%sum(nnp%output_expol)
     197       55172 :       IF (nnp%output_expol) THEN
     198           0 :          print_key => section_vals_get_subs_vals(print_section, "EXTRAPOLATION")
     199           0 :          IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     200             :             unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
     201           0 :                                            middle_name="-NNP-He-extrapolation")
     202           0 :             IF (unit_nr > 0) CALL helium_nnp_print_expol(nnp, unit_nr, ind_he)
     203           0 :             CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     204             :          END IF
     205             :       END IF
     206             : 
     207       55172 :    END SUBROUTINE helium_nnp_print
     208             : 
     209             : ! **************************************************************************************************
     210             : !> \brief Print NNP energies and standard deviation sigma
     211             : !> \param nnp ...
     212             : !> \param unit_nr ...
     213             : !> \param file_is_new ...
     214             : !> \date   2023-07-31
     215             : !> \author Laura Duran
     216             : ! **************************************************************************************************
     217           0 :    SUBROUTINE helium_nnp_print_energies(nnp, unit_nr, file_is_new)
     218             :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     219             :       INTEGER, INTENT(IN)                                :: unit_nr
     220             :       LOGICAL, INTENT(IN)                                :: file_is_new
     221             : 
     222             :       CHARACTER(len=default_path_length)                 :: fmt_string
     223             :       INTEGER                                            :: i
     224             :       REAL(KIND=dp)                                      :: std
     225             : 
     226           0 :       IF (file_is_new) THEN
     227           0 :          WRITE (unit_nr, "(A1,1X,A20)", ADVANCE='no') "#", "NNP Average [a.u.],"
     228           0 :          WRITE (unit_nr, "(A20)", ADVANCE='no') "NNP sigma [a.u.]"
     229           0 :          DO i = 1, nnp%n_committee
     230           0 :             WRITE (unit_nr, "(A17,I3)", ADVANCE='no') "NNP", i
     231             :          END DO
     232           0 :          WRITE (unit_nr, "(A)") ""
     233             :       END IF
     234             : 
     235           0 :       fmt_string = "(2X,2(F20.9))"
     236           0 :       WRITE (fmt_string, "(A,I3,A)") "(2X", nnp%n_committee + 2, "(F20.9))"
     237           0 :       std = SUM((SUM(nnp%atomic_energy, 1) - nnp%nnp_potential_energy)**2)
     238           0 :       std = std/REAL(nnp%n_committee, dp)
     239           0 :       std = SQRT(std)
     240           0 :       WRITE (unit_nr, fmt_string) nnp%nnp_potential_energy, std, SUM(nnp%atomic_energy, 1)
     241             : 
     242           0 :    END SUBROUTINE helium_nnp_print_energies
     243             : 
     244             : ! **************************************************************************************************
     245             : !> \brief Print standard deviation sigma of NNP forces
     246             : !> \param nnp ...
     247             : !> \param unit_nr ...
     248             : !> \param file_is_new ...
     249             : !> \date   2023-07-31
     250             : !> \author Laura Duran
     251             : ! **************************************************************************************************
     252           0 :    SUBROUTINE helium_nnp_print_force_sigma(nnp, unit_nr, file_is_new)
     253             :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     254             :       INTEGER, INTENT(IN)                                :: unit_nr
     255             :       LOGICAL, INTENT(IN)                                :: file_is_new
     256             : 
     257             :       INTEGER                                            :: i, ig, j
     258             :       REAL(KIND=dp), DIMENSION(3)                        :: var
     259             : 
     260           0 :       IF (unit_nr > 0) THEN
     261           0 :          IF (file_is_new) THEN
     262           0 :             WRITE (unit_nr, "(A,1X,A)") "#   NNP sigma of forces [a.u.]    x, y, z coordinates"
     263             :          END IF
     264             : 
     265           0 :          ig = 1
     266           0 :          DO i = 1, nnp%num_atoms
     267           0 :             IF (nnp%ele(i) == 'He') THEN
     268           0 :                var = 0.0_dp
     269           0 :                DO j = 1, nnp%n_committee
     270           0 :                   var = var + (nnp%committee_forces(:, i, j) - nnp%nnp_forces(:, i))**2
     271             :                END DO
     272           0 :                WRITE (unit_nr, "(A4,1X,3E20.10)") nnp%ele(i), var
     273             :             END IF
     274           0 :             ig = ig + 1
     275             :          END DO
     276             :       END IF
     277             : 
     278           0 :    END SUBROUTINE helium_nnp_print_force_sigma
     279             : 
     280             : ! **************************************************************************************************
     281             : !> \brief Print structures with extrapolation warning
     282             : !> \param nnp ...
     283             : !> \param unit_nr ...
     284             : !> \param ind_he ...
     285             : !> \date   2023-10-11
     286             : !> \author Harald Forbert (harald.forbert@rub.de)
     287             : ! **************************************************************************************************
     288           0 :    SUBROUTINE helium_nnp_print_expol(nnp, unit_nr, ind_he)
     289             :       TYPE(nnp_type), INTENT(INOUT)                      :: nnp
     290             :       INTEGER, INTENT(IN)                                :: unit_nr, ind_he
     291             : 
     292             :       CHARACTER(len=default_path_length)                 :: fmt_string
     293             :       INTEGER                                            :: i
     294             :       REAL(KIND=dp)                                      :: mass, unit_conv
     295             :       REAL(KIND=dp), DIMENSION(3)                        :: com
     296             :       TYPE(cell_type), POINTER                           :: cell
     297             : 
     298           0 :       NULLIFY (cell)
     299           0 :       CALL nnp_env_get(nnp_env=nnp, cell=cell)
     300             : 
     301           0 :       nnp%expol = nnp%expol + 1
     302           0 :       WRITE (unit_nr, *) nnp%num_atoms
     303           0 :       WRITE (unit_nr, "(A,1X,I6)") "HELIUM-NNP extrapolation point N =", nnp%expol
     304             : 
     305             :       ! move to COM of solute and wrap the box
     306             :       ! coord not needed afterwards, therefore manipulation ok
     307           0 :       com = 0.0_dp
     308           0 :       mass = 0.0_dp
     309           0 :       DO i = 1, nnp%num_atoms
     310           0 :          IF (i == ind_he) CYCLE
     311           0 :          CALL get_ptable_info(nnp%atoms(i), amass=unit_conv)
     312           0 :          com(:) = com(:) + nnp%coord(:, i)*unit_conv
     313           0 :          mass = mass + unit_conv
     314             :       END DO
     315           0 :       com(:) = com(:)/mass
     316             : 
     317           0 :       DO i = 1, nnp%num_atoms
     318           0 :          nnp%coord(:, i) = nnp%coord(:, i) - com(:)
     319           0 :          nnp%coord(:, i) = pbc(nnp%coord(:, i), cell)
     320             :       END DO
     321             : 
     322           0 :       unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM("angstrom"))
     323           0 :       fmt_string = "(A4,1X,3F20.10)"
     324           0 :       DO i = 1, nnp%num_atoms
     325             :          WRITE (unit_nr, fmt_string) &
     326           0 :             nnp%atoms(i), &
     327           0 :             nnp%coord(1, i)*unit_conv, &
     328           0 :             nnp%coord(2, i)*unit_conv, &
     329           0 :             nnp%coord(3, i)*unit_conv
     330             :       END DO
     331             : 
     332           0 :    END SUBROUTINE helium_nnp_print_expol
     333             : 
     334             : END MODULE helium_nnp

Generated by: LCOV version 1.15