LCOV - code coverage report
Current view: top level - src - force_fields_input.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:2fce0f8) Lines: 1315 1405 93.6 %
Date: 2024-12-21 06:28:57 Functions: 38 41 92.7 %

          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             : !> \par History
      10             : !>      Subroutine input_torsions changed (DG) 05-Dec-2000
      11             : !>      Output formats changed (DG) 05-Dec-2000
      12             : !>      JGH (26-01-2002) : force field parameters stored in tables, not in
      13             : !>        matrices. Input changed to have parameters labeled by the position
      14             : !>        and not atom pairs (triples etc)
      15             : !>      Teo (11.2005) : Moved all information on force field  pair_potential to
      16             : !>                      a much lighter memory structure
      17             : !>      Teo 09.2006   : Split all routines force_field I/O in a separate file
      18             : !> \author CJM
      19             : ! **************************************************************************************************
      20             : MODULE force_fields_input
      21             :    USE bibliography,                    ONLY: Clabaut2020,&
      22             :                                               Clabaut2021,&
      23             :                                               Siepmann1995,&
      24             :                                               Tersoff1988,&
      25             :                                               Tosi1964a,&
      26             :                                               Tosi1964b,&
      27             :                                               Yamada2000,&
      28             :                                               cite_reference
      29             :    USE cp_files,                        ONLY: discover_file
      30             :    USE cp_linked_list_input,            ONLY: cp_sll_val_next,&
      31             :                                               cp_sll_val_type
      32             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      33             :                                               cp_logger_type,&
      34             :                                               cp_to_string
      35             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      36             :                                               cp_print_key_unit_nr
      37             :    USE cp_parser_methods,               ONLY: parser_get_next_line
      38             :    USE cp_parser_types,                 ONLY: cp_parser_type,&
      39             :                                               parser_create,&
      40             :                                               parser_release
      41             :    USE cp_units,                        ONLY: cp_unit_to_cp2k
      42             :    USE damping_dipole_types,            ONLY: damping_info_type
      43             :    USE force_field_kind_types,          ONLY: do_ff_amber,&
      44             :                                               do_ff_charmm,&
      45             :                                               do_ff_g87,&
      46             :                                               do_ff_g96,&
      47             :                                               do_ff_opls,&
      48             :                                               do_ff_undef,&
      49             :                                               legendre_data_type
      50             :    USE force_field_types,               ONLY: force_field_type,&
      51             :                                               input_info_type
      52             :    USE force_fields_util,               ONLY: get_generic_info
      53             :    USE input_section_types,             ONLY: section_vals_get,&
      54             :                                               section_vals_get_subs_vals,&
      55             :                                               section_vals_list_get,&
      56             :                                               section_vals_type,&
      57             :                                               section_vals_val_get
      58             :    USE input_val_types,                 ONLY: val_get,&
      59             :                                               val_type
      60             :    USE kinds,                           ONLY: default_path_length,&
      61             :                                               default_string_length,&
      62             :                                               dp
      63             :    USE mathconstants,                   ONLY: pi
      64             :    USE mathlib,                         ONLY: invert_matrix
      65             :    USE memory_utilities,                ONLY: reallocate
      66             :    USE message_passing,                 ONLY: mp_para_env_type
      67             :    USE pair_potential_types,            ONLY: &
      68             :         allegro_pot_type, allegro_type, b4_type, bm_type, deepmd_type, &
      69             :         do_potential_single_allocation, ea_type, eam_pot_type, ft_pot_type, ft_type, ftd_type, &
      70             :         gal21_type, gal_type, gp_type, gw_type, ip_type, ipbv_pot_type, lj_charmm_type, &
      71             :         nequip_pot_type, nequip_type, no_potential_single_allocation, pair_potential_p_type, &
      72             :         pair_potential_reallocate, potential_single_allocation, quip_type, siepmann_type, &
      73             :         tab_pot_type, tab_type, tersoff_type, wl_type
      74             :    USE shell_potential_types,           ONLY: shell_p_create,&
      75             :                                               shell_p_type
      76             :    USE string_utilities,                ONLY: uppercase
      77             :    USE torch_api,                       ONLY: torch_allow_tf32,&
      78             :                                               torch_model_read_metadata
      79             : #include "./base/base_uses.f90"
      80             : 
      81             :    IMPLICIT NONE
      82             : 
      83             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_fields_input'
      84             : 
      85             :    PRIVATE
      86             :    PUBLIC :: read_force_field_section, &
      87             :              read_lj_section, &
      88             :              read_wl_section, &
      89             :              read_gd_section, &
      90             :              read_gp_section, &
      91             :              read_chrg_section
      92             : 
      93             : CONTAINS
      94             : 
      95             : ! **************************************************************************************************
      96             : !> \brief Reads the force_field input section
      97             : !> \param ff_section ...
      98             : !> \param mm_section ...
      99             : !> \param ff_type ...
     100             : !> \param para_env ...
     101             : !> \author teo
     102             : ! **************************************************************************************************
     103       36946 :    SUBROUTINE read_force_field_section1(ff_section, mm_section, ff_type, para_env)
     104             :       TYPE(section_vals_type), POINTER                   :: ff_section, mm_section
     105             :       TYPE(force_field_type), INTENT(INOUT)              :: ff_type
     106             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     107             : 
     108             :       CHARACTER(LEN=default_string_length), &
     109        2639 :          DIMENSION(:), POINTER                           :: atm_names
     110             :       INTEGER :: nallegro, nb4, nbends, nbm, nbmhft, nbmhftd, nbonds, nchg, ndeepmd, neam, ngal, &
     111             :          ngal21, ngd, ngp, nimpr, nipbv, nlj, nnequip, nopbend, nquip, nshell, nsiepmann, ntab, &
     112             :          ntersoff, ntors, ntot, nubs, nwl
     113             :       LOGICAL                                            :: explicit, unique_spline
     114             :       REAL(KIND=dp)                                      :: min_eps_spline_allowed
     115             :       TYPE(input_info_type), POINTER                     :: inp_info
     116             :       TYPE(section_vals_type), POINTER                   :: tmp_section, tmp_section2
     117             : 
     118             :       INTEGER::i
     119             : 
     120        2639 :       NULLIFY (tmp_section, tmp_section2)
     121        2639 :       inp_info => ff_type%inp_info
     122        2639 :       CALL section_vals_val_get(ff_section, "PARMTYPE", i_val=ff_type%ff_type)
     123        2639 :       CALL section_vals_val_get(ff_section, "EI_SCALE14", r_val=ff_type%ei_scale14)
     124        2639 :       CALL section_vals_val_get(ff_section, "VDW_SCALE14", r_val=ff_type%vdw_scale14)
     125        2639 :       CALL section_vals_val_get(ff_section, "SPLINE%RCUT_NB", r_val=ff_type%rcut_nb)
     126        2639 :       CALL section_vals_val_get(ff_section, "SPLINE%R0_NB", r_val=ff_type%rlow_nb)
     127        2639 :       CALL section_vals_val_get(ff_section, "SPLINE%EPS_SPLINE", r_val=ff_type%eps_spline)
     128        2639 :       CALL section_vals_val_get(ff_section, "SPLINE%EMAX_SPLINE", r_val=ff_type%emax_spline)
     129        2639 :       CALL section_vals_val_get(ff_section, "SPLINE%EMAX_ACCURACY", r_val=ff_type%max_energy)
     130        2639 :       CALL section_vals_val_get(ff_section, "SPLINE%NPOINTS", i_val=ff_type%npoints)
     131        2639 :       CALL section_vals_val_get(ff_section, "IGNORE_MISSING_CRITICAL_PARAMS", l_val=ff_type%ignore_missing_critical)
     132        2639 :       CPASSERT(ff_type%max_energy <= ff_type%emax_spline)
     133             :       ! Read the parameter file name only if the force field type requires it..
     134        3543 :       SELECT CASE (ff_type%ff_type)
     135             :       CASE (do_ff_charmm, do_ff_amber, do_ff_g96, do_ff_g87)
     136         904 :          CALL section_vals_val_get(ff_section, "PARM_FILE_NAME", c_val=ff_type%ff_file_name)
     137         904 :          IF (TRIM(ff_type%ff_file_name) == "") &
     138           0 :             CPABORT("Force Field Parameter's filename is empty! Please check your input file.")
     139             :       CASE (do_ff_undef)
     140             :          ! Do Nothing
     141             :       CASE DEFAULT
     142        2639 :          CPABORT("Force field type not implemented")
     143             :       END SELECT
     144             :       ! Numerical Accuracy:
     145             :       ! the factors here should depend on the energy and on the shape of each potential mapped
     146             :       ! with splines. this would make everything un-necessarily complicated. Let's just be safe
     147             :       ! and assume that we cannot achieve an accuracy on the spline 2 orders of magnitude more
     148             :       ! than the smallest representable number (taking into account also the max_energy for the
     149             :       ! spline generation
     150        2639 :       min_eps_spline_allowed = 20.0_dp*MAX(ff_type%max_energy, 10.0_dp)*EPSILON(0.0_dp)
     151        2639 :       IF (ff_type%eps_spline < min_eps_spline_allowed) THEN
     152             :          CALL cp_warn(__LOCATION__, &
     153             :                       "Requested spline accuracy ("//TRIM(cp_to_string(ff_type%eps_spline))//" ) "// &
     154             :                       "is smaller than the minimum value allowed ("//TRIM(cp_to_string(min_eps_spline_allowed))// &
     155             :                       " ) with the present machine precision ("//TRIM(cp_to_string(EPSILON(0.0_dp)))//" ). "// &
     156           0 :                       "New EPS_SPLINE value ("//TRIM(cp_to_string(min_eps_spline_allowed))//" ). ")
     157           0 :          ff_type%eps_spline = min_eps_spline_allowed
     158             :       END IF
     159        2639 :       CALL section_vals_val_get(ff_section, "SHIFT_CUTOFF", l_val=ff_type%shift_cutoff)
     160        2639 :       CALL section_vals_val_get(ff_section, "SPLINE%UNIQUE_SPLINE", l_val=unique_spline)
     161             :       ! Single spline
     162        2639 :       potential_single_allocation = no_potential_single_allocation
     163        2639 :       IF (unique_spline) potential_single_allocation = do_potential_single_allocation
     164             : 
     165        2639 :       CALL section_vals_val_get(ff_section, "MULTIPLE_POTENTIAL", l_val=ff_type%multiple_potential)
     166        2639 :       CALL section_vals_val_get(ff_section, "DO_NONBONDED", l_val=ff_type%do_nonbonded)
     167        2639 :       CALL section_vals_val_get(ff_section, "DO_ELECTROSTATICS", l_val=ff_type%do_electrostatics)
     168        2639 :       tmp_section => section_vals_get_subs_vals(ff_section, "NONBONDED")
     169        2639 :       CALL section_vals_get(tmp_section, explicit=explicit)
     170        2639 :       IF (explicit .AND. ff_type%do_nonbonded) THEN
     171        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "LENNARD-JONES")
     172        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nlj)
     173        1735 :          ntot = 0
     174        1735 :          IF (explicit) THEN
     175         976 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nlj, lj_charmm=.TRUE.)
     176         976 :             CALL read_lj_section(inp_info%nonbonded, tmp_section2, ntot)
     177             :          END IF
     178             : 
     179        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "WILLIAMS")
     180        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nwl)
     181        1735 :          ntot = nlj
     182        1735 :          IF (explicit) THEN
     183         373 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nwl, williams=.TRUE.)
     184         373 :             CALL read_wl_section(inp_info%nonbonded, tmp_section2, ntot)
     185             :          END IF
     186             : 
     187        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "EAM")
     188        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=neam)
     189        1735 :          ntot = nlj + nwl
     190        1735 :          IF (explicit) THEN
     191          12 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + neam, eam=.TRUE.)
     192          12 :             CALL read_eam_section(inp_info%nonbonded, tmp_section2, ntot, para_env, mm_section)
     193             :          END IF
     194             : 
     195        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "GOODWIN")
     196        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ngd)
     197        1735 :          ntot = nlj + nwl + neam
     198        1735 :          IF (explicit) THEN
     199           0 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ngd, goodwin=.TRUE.)
     200           0 :             CALL read_gd_section(inp_info%nonbonded, tmp_section2, ntot)
     201             :          END IF
     202             : 
     203        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "IPBV")
     204        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nipbv)
     205        1735 :          ntot = nlj + nwl + neam + ngd
     206        1735 :          IF (explicit) THEN
     207          16 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nipbv, ipbv=.TRUE.)
     208          16 :             CALL read_ipbv_section(inp_info%nonbonded, tmp_section2, ntot)
     209             :          END IF
     210             : 
     211        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "BMHFT")
     212        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nbmhft)
     213        1735 :          ntot = nlj + nwl + neam + ngd + nipbv
     214        1735 :          IF (explicit) THEN
     215           4 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nbmhft, bmhft=.TRUE.)
     216           4 :             CALL read_bmhft_section(inp_info%nonbonded, tmp_section2, ntot)
     217             :          END IF
     218             : 
     219        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "BMHFTD")
     220        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nbmhftd)
     221        1735 :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft
     222        1735 :          IF (explicit) THEN
     223          18 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nbmhftd, bmhftd=.TRUE.)
     224          18 :             CALL read_bmhftd_section(inp_info%nonbonded, tmp_section2, ntot)
     225             :          END IF
     226             : 
     227        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "BUCK4RANGES")
     228        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nb4)
     229        1735 :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd
     230        1735 :          IF (explicit) THEN
     231         262 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nb4, buck4r=.TRUE.)
     232         262 :             CALL read_b4_section(inp_info%nonbonded, tmp_section2, ntot)
     233             :          END IF
     234             : 
     235        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "BUCKMORSE")
     236        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nbm)
     237        1735 :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4
     238        1735 :          IF (explicit) THEN
     239           8 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nbm, buckmo=.TRUE.)
     240           8 :             CALL read_bm_section(inp_info%nonbonded, tmp_section2, ntot)
     241             :          END IF
     242             : 
     243        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "GENPOT")
     244        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ngp)
     245        1735 :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm
     246        1735 :          IF (explicit) THEN
     247         308 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ngp, gp=.TRUE.)
     248         308 :             CALL read_gp_section(inp_info%nonbonded, tmp_section2, ntot)
     249             :          END IF
     250        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "TERSOFF")
     251        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ntersoff)
     252        1735 :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp
     253        1735 :          IF (explicit) THEN
     254          40 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ntersoff, tersoff=.TRUE.)
     255          40 :             CALL read_tersoff_section(inp_info%nonbonded, tmp_section2, ntot, tmp_section2)
     256             :          END IF
     257             : 
     258        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "GAL19")
     259        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ngal)
     260        1735 :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff
     261        1735 :          IF (explicit) THEN
     262           1 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ngal, gal=.TRUE.)
     263           1 :             CALL read_gal_section(inp_info%nonbonded, tmp_section2, ntot, tmp_section2)
     264             :          END IF
     265             : 
     266        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "GAL21")
     267        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ngal21)
     268        1735 :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + ngal
     269        1735 :          IF (explicit) THEN
     270           1 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ngal21, gal21=.TRUE.)
     271           1 :             CALL read_gal21_section(inp_info%nonbonded, tmp_section2, ntot, tmp_section2)
     272             :          END IF
     273             : 
     274        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "SIEPMANN")
     275        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nsiepmann)
     276        1735 :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + ngal + ngal21
     277        1735 :          IF (explicit) THEN
     278           5 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nsiepmann, siepmann=.TRUE.)
     279           5 :             CALL read_siepmann_section(inp_info%nonbonded, tmp_section2, ntot, tmp_section2)
     280             :          END IF
     281             : 
     282        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "quip")
     283        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nquip)
     284        1735 :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + ngal + ngal21 + nsiepmann
     285        1735 :          IF (explicit) THEN
     286           2 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nquip, quip=.TRUE.)
     287           2 :             CALL read_quip_section(inp_info%nonbonded, tmp_section2, ntot)
     288             :          END IF
     289             : 
     290        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "nequip")
     291        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nnequip)
     292             :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + ngal + ngal21 + nsiepmann + &
     293        1735 :                 nquip
     294        1735 :          IF (explicit) THEN
     295             :             ! avoid repeating the nequip section for each pair
     296           4 :             CALL section_vals_val_get(tmp_section2, "ATOMS", c_vals=atm_names)
     297           4 :             nnequip = nnequip - 1 + SIZE(atm_names) + (SIZE(atm_names)*SIZE(atm_names) - SIZE(atm_names))/2
     298           4 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nnequip, nequip=.TRUE.)
     299           4 :             CALL read_nequip_section(inp_info%nonbonded, tmp_section2, ntot)
     300             :          END IF
     301             : 
     302        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "allegro")
     303        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nallegro)
     304             :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + ngal + ngal21 + nsiepmann + &
     305        1735 :                 nquip + nnequip
     306        1735 :          IF (explicit) THEN
     307             :             ! avoid repeating the allegro section for each pair
     308           4 :             CALL section_vals_val_get(tmp_section2, "ATOMS", c_vals=atm_names)
     309           4 :             nallegro = nallegro - 1 + SIZE(atm_names) + (SIZE(atm_names)*SIZE(atm_names) - SIZE(atm_names))/2
     310           4 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nallegro, allegro=.TRUE.)
     311           4 :             CALL read_allegro_section(inp_info%nonbonded, tmp_section2, ntot)
     312             :          END IF
     313             : 
     314        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "TABPOT")
     315        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ntab)
     316             :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + ngal + ngal21 + nsiepmann + &
     317        1735 :                 nquip + nnequip + nallegro
     318        1735 :          IF (explicit) THEN
     319           8 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ntab, tab=.TRUE.)
     320           8 :             CALL read_tabpot_section(inp_info%nonbonded, tmp_section2, ntot, para_env, mm_section)
     321             :          END IF
     322             : 
     323        1735 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "DEEPMD")
     324        1735 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ndeepmd)
     325             :          ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + ngal + ngal21 + nsiepmann + &
     326        1735 :                 nquip + nnequip + nallegro + ntab
     327        1735 :          IF (explicit) THEN
     328             :             ! avoid repeating the deepmd section for each pair
     329           2 :             CALL section_vals_val_get(tmp_section2, "ATOMS", c_vals=atm_names)
     330           2 :             ndeepmd = ndeepmd - 1 + SIZE(atm_names) + (SIZE(atm_names)*SIZE(atm_names) - SIZE(atm_names))/2
     331           2 :             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ndeepmd, deepmd=.TRUE.)
     332           2 :             CALL read_deepmd_section(inp_info%nonbonded, tmp_section2, ntot)
     333             :          END IF
     334             : 
     335             :       END IF
     336             : 
     337        2639 :       tmp_section => section_vals_get_subs_vals(ff_section, "NONBONDED14")
     338        2639 :       CALL section_vals_get(tmp_section, explicit=explicit)
     339        2639 :       IF (explicit .AND. ff_type%do_nonbonded) THEN
     340         274 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "LENNARD-JONES")
     341         274 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nlj)
     342         274 :          ntot = 0
     343         274 :          IF (explicit) THEN
     344          12 :             CALL pair_potential_reallocate(inp_info%nonbonded14, 1, ntot + nlj, lj_charmm=.TRUE.)
     345          12 :             CALL read_lj_section(inp_info%nonbonded14, tmp_section2, ntot)
     346             :          END IF
     347         274 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "WILLIAMS")
     348         274 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=nwl)
     349         274 :          ntot = nlj
     350         274 :          IF (explicit) THEN
     351           0 :             CALL pair_potential_reallocate(inp_info%nonbonded14, 1, ntot + nwl, williams=.TRUE.)
     352           0 :             CALL read_wl_section(inp_info%nonbonded14, tmp_section2, ntot)
     353             :          END IF
     354         274 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "GOODWIN")
     355         274 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ngd)
     356         274 :          ntot = nlj + nwl
     357         274 :          IF (explicit) THEN
     358           0 :             CALL pair_potential_reallocate(inp_info%nonbonded14, 1, ntot + ngd, goodwin=.TRUE.)
     359           0 :             CALL read_gd_section(inp_info%nonbonded14, tmp_section2, ntot)
     360             :          END IF
     361         274 :          tmp_section2 => section_vals_get_subs_vals(tmp_section, "GENPOT")
     362         274 :          CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ngp)
     363         274 :          ntot = nlj + nwl + ngd
     364         274 :          IF (explicit) THEN
     365         262 :             CALL pair_potential_reallocate(inp_info%nonbonded14, 1, ntot + ngp, gp=.TRUE.)
     366         262 :             CALL read_gp_section(inp_info%nonbonded14, tmp_section2, ntot)
     367             :          END IF
     368             :       END IF
     369             : 
     370        2639 :       tmp_section => section_vals_get_subs_vals(ff_section, "CHARGE")
     371        2639 :       CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nchg)
     372        2639 :       IF (explicit) THEN
     373        2079 :          ntot = 0
     374        2079 :          CALL reallocate(inp_info%charge_atm, 1, nchg)
     375        2079 :          CALL reallocate(inp_info%charge, 1, nchg)
     376        2079 :          CALL read_chrg_section(inp_info%charge_atm, inp_info%charge, tmp_section, ntot)
     377             :       END IF
     378        2639 :       tmp_section => section_vals_get_subs_vals(ff_section, "DIPOLE")
     379        2639 :       CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nchg)
     380        2639 :       IF (explicit) THEN
     381          34 :          ntot = 0
     382          34 :          CALL reallocate(inp_info%apol_atm, 1, nchg)
     383          34 :          CALL reallocate(inp_info%apol, 1, nchg)
     384             :          CALL read_apol_section(inp_info%apol_atm, inp_info%apol, inp_info%damping_list, &
     385          34 :                                 tmp_section, ntot)
     386             :       END IF
     387        2639 :       tmp_section => section_vals_get_subs_vals(ff_section, "QUADRUPOLE")
     388        2639 :       CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nchg)
     389        2639 :       IF (explicit) THEN
     390           0 :          ntot = 0
     391           0 :          CALL reallocate(inp_info%cpol_atm, 1, nchg)
     392           0 :          CALL reallocate(inp_info%cpol, 1, nchg)
     393           0 :          CALL read_cpol_section(inp_info%cpol_atm, inp_info%cpol, tmp_section, ntot)
     394             :       END IF
     395        2639 :       tmp_section => section_vals_get_subs_vals(ff_section, "SHELL")
     396        2639 :       CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nshell)
     397        2639 :       IF (explicit) THEN
     398         258 :          ntot = 0
     399         258 :          CALL shell_p_create(inp_info%shell_list, nshell)
     400         258 :          CALL read_shell_section(inp_info%shell_list, tmp_section, ntot)
     401             :       END IF
     402             : 
     403        2639 :       tmp_section => section_vals_get_subs_vals(ff_section, "BOND")
     404        2639 :       CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nbonds)
     405        2639 :       IF (explicit) THEN
     406         975 :          ntot = 0
     407         975 :          CALL reallocate(inp_info%bond_kind, 1, nbonds)
     408         975 :          CALL reallocate(inp_info%bond_a, 1, nbonds)
     409         975 :          CALL reallocate(inp_info%bond_b, 1, nbonds)
     410         975 :          CALL reallocate(inp_info%bond_k, 1, 3, 1, nbonds)
     411         975 :          CALL reallocate(inp_info%bond_r0, 1, nbonds)
     412         975 :          CALL reallocate(inp_info%bond_cs, 1, nbonds)
     413             :          CALL read_bonds_section(inp_info%bond_kind, inp_info%bond_a, inp_info%bond_b, inp_info%bond_k, &
     414         975 :                                  inp_info%bond_r0, inp_info%bond_cs, tmp_section, ntot)
     415             :       END IF
     416        2639 :       tmp_section => section_vals_get_subs_vals(ff_section, "BEND")
     417        2639 :       CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nbends)
     418        2639 :       IF (explicit) THEN
     419         939 :          ntot = 0
     420         939 :          CALL reallocate(inp_info%bend_kind, 1, nbends)
     421         939 :          CALL reallocate(inp_info%bend_a, 1, nbends)
     422         939 :          CALL reallocate(inp_info%bend_b, 1, nbends)
     423         939 :          CALL reallocate(inp_info%bend_c, 1, nbends)
     424         939 :          CALL reallocate(inp_info%bend_k, 1, nbends)
     425         939 :          CALL reallocate(inp_info%bend_theta0, 1, nbends)
     426         939 :          CALL reallocate(inp_info%bend_cb, 1, nbends)
     427         939 :          CALL reallocate(inp_info%bend_r012, 1, nbends)
     428         939 :          CALL reallocate(inp_info%bend_r032, 1, nbends)
     429         939 :          CALL reallocate(inp_info%bend_kbs12, 1, nbends)
     430         939 :          CALL reallocate(inp_info%bend_kbs32, 1, nbends)
     431         939 :          CALL reallocate(inp_info%bend_kss, 1, nbends)
     432         939 :          IF (ASSOCIATED(inp_info%bend_legendre)) THEN
     433           0 :             DO i = 1, SIZE(inp_info%bend_legendre)
     434           0 :                IF (ASSOCIATED(inp_info%bend_legendre(i)%coeffs)) THEN
     435           0 :                   DEALLOCATE (inp_info%bend_legendre(i)%coeffs)
     436           0 :                   NULLIFY (inp_info%bend_legendre(i)%coeffs)
     437             :                END IF
     438             :             END DO
     439           0 :             DEALLOCATE (inp_info%bend_legendre)
     440             :             NULLIFY (inp_info%bend_legendre)
     441             :          END IF
     442        4936 :          ALLOCATE (inp_info%bend_legendre(1:nbends))
     443        3058 :          DO i = 1, SIZE(inp_info%bend_legendre(1:nbends))
     444        2119 :             NULLIFY (inp_info%bend_legendre(i)%coeffs)
     445        3058 :             inp_info%bend_legendre(i)%order = 0
     446             :          END DO
     447             :          CALL read_bends_section(inp_info%bend_kind, inp_info%bend_a, inp_info%bend_b, inp_info%bend_c, &
     448             :                                  inp_info%bend_k, inp_info%bend_theta0, inp_info%bend_cb, &
     449             :                                  inp_info%bend_r012, inp_info%bend_r032, inp_info%bend_kbs12, &
     450             :                                  inp_info%bend_kbs32, inp_info%bend_kss, &
     451         939 :                                  inp_info%bend_legendre, tmp_section, ntot)
     452             :       END IF
     453        2639 :       tmp_section => section_vals_get_subs_vals(ff_section, "BEND")
     454        2639 :       CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nubs)
     455        2639 :       IF (explicit) THEN
     456         939 :          ntot = 0
     457         939 :          CALL reallocate(inp_info%ub_kind, 1, nubs)
     458         939 :          CALL reallocate(inp_info%ub_a, 1, nubs)
     459         939 :          CALL reallocate(inp_info%ub_b, 1, nubs)
     460         939 :          CALL reallocate(inp_info%ub_c, 1, nubs)
     461         939 :          CALL reallocate(inp_info%ub_k, 1, 3, 1, nubs)
     462         939 :          CALL reallocate(inp_info%ub_r0, 1, nubs)
     463             :          CALL read_ubs_section(inp_info%ub_kind, inp_info%ub_a, inp_info%ub_b, inp_info%ub_c, &
     464         939 :                                inp_info%ub_k, inp_info%ub_r0, tmp_section, ntot)
     465             :       END IF
     466        2639 :       tmp_section => section_vals_get_subs_vals(ff_section, "TORSION")
     467        2639 :       CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=ntors)
     468        2639 :       IF (explicit) THEN
     469           6 :          ntot = 0
     470           6 :          CALL reallocate(inp_info%torsion_kind, 1, ntors)
     471           6 :          CALL reallocate(inp_info%torsion_a, 1, ntors)
     472           6 :          CALL reallocate(inp_info%torsion_b, 1, ntors)
     473           6 :          CALL reallocate(inp_info%torsion_c, 1, ntors)
     474           6 :          CALL reallocate(inp_info%torsion_d, 1, ntors)
     475           6 :          CALL reallocate(inp_info%torsion_k, 1, ntors)
     476           6 :          CALL reallocate(inp_info%torsion_m, 1, ntors)
     477           6 :          CALL reallocate(inp_info%torsion_phi0, 1, ntors)
     478             :          CALL read_torsions_section(inp_info%torsion_kind, inp_info%torsion_a, inp_info%torsion_b, &
     479             :                                     inp_info%torsion_c, inp_info%torsion_d, inp_info%torsion_k, inp_info%torsion_phi0, &
     480           6 :                                     inp_info%torsion_m, tmp_section, ntot)
     481             :       END IF
     482             : 
     483        2639 :       tmp_section => section_vals_get_subs_vals(ff_section, "IMPROPER")
     484        2639 :       CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nimpr)
     485        2639 :       IF (explicit) THEN
     486           8 :          ntot = 0
     487           8 :          CALL reallocate(inp_info%impr_kind, 1, nimpr)
     488           8 :          CALL reallocate(inp_info%impr_a, 1, nimpr)
     489           8 :          CALL reallocate(inp_info%impr_b, 1, nimpr)
     490           8 :          CALL reallocate(inp_info%impr_c, 1, nimpr)
     491           8 :          CALL reallocate(inp_info%impr_d, 1, nimpr)
     492           8 :          CALL reallocate(inp_info%impr_k, 1, nimpr)
     493           8 :          CALL reallocate(inp_info%impr_phi0, 1, nimpr)
     494             :          CALL read_improper_section(inp_info%impr_kind, inp_info%impr_a, inp_info%impr_b, &
     495             :                                     inp_info%impr_c, inp_info%impr_d, inp_info%impr_k, inp_info%impr_phi0, &
     496           8 :                                     tmp_section, ntot)
     497             :       END IF
     498             : 
     499        2639 :       tmp_section => section_vals_get_subs_vals(ff_section, "OPBEND")
     500        2639 :       CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nopbend)
     501        2639 :       IF (explicit) THEN
     502           2 :          ntot = 0
     503           2 :          CALL reallocate(inp_info%opbend_kind, 1, nopbend)
     504           2 :          CALL reallocate(inp_info%opbend_a, 1, nopbend)
     505           2 :          CALL reallocate(inp_info%opbend_b, 1, nopbend)
     506           2 :          CALL reallocate(inp_info%opbend_c, 1, nopbend)
     507           2 :          CALL reallocate(inp_info%opbend_d, 1, nopbend)
     508           2 :          CALL reallocate(inp_info%opbend_k, 1, nopbend)
     509           2 :          CALL reallocate(inp_info%opbend_phi0, 1, nopbend)
     510             :          CALL read_opbend_section(inp_info%opbend_kind, inp_info%opbend_a, inp_info%opbend_b, &
     511             :                                   inp_info%opbend_c, inp_info%opbend_d, inp_info%opbend_k, inp_info%opbend_phi0, &
     512           2 :                                   tmp_section, ntot)
     513             :       END IF
     514             : 
     515        2639 :    END SUBROUTINE read_force_field_section1
     516             : 
     517             : ! **************************************************************************************************
     518             : !> \brief Set up of the IPBV force fields
     519             : !> \param at1 ...
     520             : !> \param at2 ...
     521             : !> \param ipbv ...
     522             : !> \author teo
     523             : ! **************************************************************************************************
     524          48 :    SUBROUTINE set_IPBV_ff(at1, at2, ipbv)
     525             :       CHARACTER(LEN=*), INTENT(IN)                       :: at1, at2
     526             :       TYPE(ipbv_pot_type), POINTER                       :: ipbv
     527             : 
     528          48 :       IF ((at1(1:1) == 'O') .AND. (at2(1:1) == 'O')) THEN
     529          16 :          ipbv%rcore = 0.9_dp ! a.u.
     530          16 :          ipbv%m = -1.2226442563398141E+11_dp ! Kelvin/a.u.
     531          16 :          ipbv%b = 1.1791292385486696E+11_dp ! Hartree
     532             : 
     533             :          ! Hartree*a.u.^2
     534          16 :          ipbv%a(2) = 4.786380682394_dp
     535          16 :          ipbv%a(3) = -1543.407053545_dp
     536          16 :          ipbv%a(4) = 88783.31188529_dp
     537          16 :          ipbv%a(5) = -2361200.155376_dp
     538          16 :          ipbv%a(6) = 35940504.84679_dp
     539          16 :          ipbv%a(7) = -339762743.6358_dp
     540          16 :          ipbv%a(8) = 2043874926.466_dp
     541          16 :          ipbv%a(9) = -7654856796.383_dp
     542          16 :          ipbv%a(10) = 16195251405.65_dp
     543          16 :          ipbv%a(11) = -13140392992.18_dp
     544          16 :          ipbv%a(12) = -9285572894.245_dp
     545          16 :          ipbv%a(13) = 8756947519.029_dp
     546          16 :          ipbv%a(14) = 15793297761.67_dp
     547          16 :          ipbv%a(15) = 12917180227.21_dp
     548          32 :       ELSEIF (((at1(1:1) == 'O') .AND. (at2(1:1) == 'H')) .OR. &
     549             :               ((at1(1:1) == 'H') .AND. (at2(1:1) == 'O'))) THEN
     550          16 :          ipbv%rcore = 2.95_dp ! a.u.
     551             : 
     552          16 :          ipbv%m = -0.004025691139759147_dp ! Hartree/a.u.
     553          16 :          ipbv%b = -2.193731138097428_dp ! Hartree
     554             :          ! Hartree*a.u.^2
     555          16 :          ipbv%a(2) = -195.7716013277_dp
     556          16 :          ipbv%a(3) = 15343.78613395_dp
     557          16 :          ipbv%a(4) = -530864.4586516_dp
     558          16 :          ipbv%a(5) = 10707934.39058_dp
     559          16 :          ipbv%a(6) = -140099704.7890_dp
     560          16 :          ipbv%a(7) = 1250943273.785_dp
     561          16 :          ipbv%a(8) = -7795458330.676_dp
     562          16 :          ipbv%a(9) = 33955897217.31_dp
     563          16 :          ipbv%a(10) = -101135640744.0_dp
     564          16 :          ipbv%a(11) = 193107995718.7_dp
     565          16 :          ipbv%a(12) = -193440560940.0_dp
     566          16 :          ipbv%a(13) = -4224406093.918E0_dp
     567          16 :          ipbv%a(14) = 217192386506.5E0_dp
     568          16 :          ipbv%a(15) = -157581228915.5_dp
     569          16 :       ELSEIF ((at1(1:1) == 'H') .AND. (at2(1:1) == 'H')) THEN
     570          16 :          ipbv%rcore = 3.165_dp ! a.u.
     571          16 :          ipbv%m = 0.002639704108787555_dp ! Hartree/a.u.
     572          16 :          ipbv%b = -0.2735482611857583_dp ! Hartree
     573             :          ! Hartree*a.u.^2
     574          16 :          ipbv%a(2) = -26.29456010782_dp
     575          16 :          ipbv%a(3) = 2373.352548248_dp
     576          16 :          ipbv%a(4) = -93880.43551360_dp
     577          16 :          ipbv%a(5) = 2154624.884809_dp
     578          16 :          ipbv%a(6) = -31965151.34955_dp
     579          16 :          ipbv%a(7) = 322781785.3278_dp
     580          16 :          ipbv%a(8) = -2271097368.668_dp
     581          16 :          ipbv%a(9) = 11169163192.90_dp
     582          16 :          ipbv%a(10) = -37684457778.47_dp
     583          16 :          ipbv%a(11) = 82562104256.03_dp
     584          16 :          ipbv%a(12) = -100510435213.4_dp
     585          16 :          ipbv%a(13) = 24570342714.65E0_dp
     586          16 :          ipbv%a(14) = 88766181532.94E0_dp
     587          16 :          ipbv%a(15) = -79705131323.98_dp
     588             :       ELSE
     589           0 :          CPABORT("IPBV only for WATER")
     590             :       END IF
     591          48 :    END SUBROUTINE set_IPBV_ff
     592             : 
     593             : ! **************************************************************************************************
     594             : !> \brief Set up of the BMHFT force fields
     595             : !> \param at1 ...
     596             : !> \param at2 ...
     597             : !> \param ft ...
     598             : !> \author teo
     599             : ! **************************************************************************************************
     600          12 :    SUBROUTINE set_BMHFT_ff(at1, at2, ft)
     601             :       CHARACTER(LEN=*), INTENT(IN)                       :: at1, at2
     602             :       TYPE(ft_pot_type), POINTER                         :: ft
     603             : 
     604          12 :       ft%b = cp_unit_to_cp2k(3.1545_dp, "angstrom^-1")
     605          12 :       IF ((at1(1:2) == 'NA') .AND. (at2(1:2) == 'NA')) THEN
     606           4 :          ft%a = cp_unit_to_cp2k(424.097_dp, "eV")
     607           4 :          ft%c = cp_unit_to_cp2k(1.05_dp, "eV*angstrom^6")
     608           4 :          ft%d = cp_unit_to_cp2k(0.499_dp, "eV*angstrom^8")
     609           8 :       ELSEIF (((at1(1:2) == 'NA') .AND. (at2(1:2) == 'CL')) .OR. &
     610             :               ((at1(1:2) == 'CL') .AND. (at2(1:2) == 'NA'))) THEN
     611           4 :          ft%a = cp_unit_to_cp2k(1256.31_dp, "eV")
     612           4 :          ft%c = cp_unit_to_cp2k(7.00_dp, "eV*angstrom^6")
     613           4 :          ft%d = cp_unit_to_cp2k(8.676_dp, "eV*angstrom^8")
     614           4 :       ELSEIF ((at1(1:2) == 'CL') .AND. (at2(1:2) == 'CL')) THEN
     615           4 :          ft%a = cp_unit_to_cp2k(3488.998_dp, "eV")
     616           4 :          ft%c = cp_unit_to_cp2k(72.50_dp, "eV*angstrom^6")
     617           4 :          ft%d = cp_unit_to_cp2k(145.427_dp, "eV*angstrom^8")
     618             :       ELSE
     619           0 :          CPABORT("BMHFT only for NaCl")
     620             :       END IF
     621             : 
     622          12 :    END SUBROUTINE set_BMHFT_ff
     623             : 
     624             : ! **************************************************************************************************
     625             : !> \brief Set up of the BMHFTD force fields
     626             : !> \author Mathieu Salanne 05.2010
     627             : ! **************************************************************************************************
     628           0 :    SUBROUTINE set_BMHFTD_ff()
     629             : 
     630           0 :       CPABORT("No default parameters present for BMHFTD")
     631             : 
     632           0 :    END SUBROUTINE set_BMHFTD_ff
     633             : 
     634             : ! **************************************************************************************************
     635             : !> \brief Reads the EAM section
     636             : !> \param nonbonded ...
     637             : !> \param section ...
     638             : !> \param start ...
     639             : !> \param para_env ...
     640             : !> \param mm_section ...
     641             : !> \author teo
     642             : ! **************************************************************************************************
     643          12 :    SUBROUTINE read_eam_section(nonbonded, section, start, para_env, mm_section)
     644             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
     645             :       TYPE(section_vals_type), POINTER                   :: section
     646             :       INTEGER, INTENT(IN)                                :: start
     647             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     648             :       TYPE(section_vals_type), POINTER                   :: mm_section
     649             : 
     650             :       CHARACTER(LEN=default_string_length), &
     651          12 :          DIMENSION(:), POINTER                           :: atm_names
     652             :       INTEGER                                            :: isec, n_items
     653             : 
     654          12 :       CALL section_vals_get(section, n_repetition=n_items)
     655          32 :       DO isec = 1, n_items
     656          20 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
     657             : 
     658          40 :          nonbonded%pot(start + isec)%pot%type = ea_type
     659          20 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
     660          20 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
     661          20 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
     662          20 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
     663             :          CALL section_vals_val_get(section, "PARM_FILE_NAME", i_rep_section=isec, &
     664          20 :                                    c_val=nonbonded%pot(start + isec)%pot%set(1)%eam%eam_file_name)
     665          20 :          CALL read_eam_data(nonbonded%pot(start + isec)%pot%set(1)%eam, para_env, mm_section)
     666          32 :          nonbonded%pot(start + isec)%pot%rcutsq = nonbonded%pot(start + isec)%pot%set(1)%eam%acutal**2
     667             :       END DO
     668          12 :    END SUBROUTINE read_eam_section
     669             : 
     670             : ! **************************************************************************************************
     671             : !> \brief Reads the DEEPMD section
     672             : !> \param nonbonded ...
     673             : !> \param section ...
     674             : !> \param start ...
     675             : !> \author teo
     676             : ! **************************************************************************************************
     677           2 :    SUBROUTINE read_deepmd_section(nonbonded, section, start)
     678             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
     679             :       TYPE(section_vals_type), POINTER                   :: section
     680             :       INTEGER, INTENT(IN)                                :: start
     681             : 
     682             :       CHARACTER(LEN=default_string_length)               :: deepmd_file_name
     683             :       CHARACTER(LEN=default_string_length), &
     684           2 :          DIMENSION(:), POINTER                           :: atm_names
     685             :       INTEGER                                            :: isec, jsec, n_items
     686           2 :       INTEGER, DIMENSION(:), POINTER                     :: atm_deepmd_types
     687             : 
     688             :       n_items = 1
     689           2 :       isec = 1
     690           2 :       n_items = isec*n_items
     691           2 :       CALL section_vals_val_get(section, "ATOMS", c_vals=atm_names)
     692           2 :       CALL section_vals_val_get(section, "ATOMS_DEEPMD_TYPE", i_vals=atm_deepmd_types)
     693           2 :       CALL section_vals_val_get(section, "POT_FILE_NAME", c_val=deepmd_file_name)
     694             : 
     695           4 :       DO isec = 1, SIZE(atm_names)
     696           6 :          DO jsec = isec, SIZE(atm_names)
     697           4 :             nonbonded%pot(start + n_items)%pot%type = deepmd_type
     698           2 :             nonbonded%pot(start + n_items)%pot%at1 = atm_names(isec)
     699           2 :             nonbonded%pot(start + n_items)%pot%at2 = atm_names(jsec)
     700           2 :             CALL uppercase(nonbonded%pot(start + n_items)%pot%at1)
     701           2 :             CALL uppercase(nonbonded%pot(start + n_items)%pot%at2)
     702             : 
     703           2 :             nonbonded%pot(start + n_items)%pot%set(1)%deepmd%deepmd_file_name = discover_file(deepmd_file_name)
     704           2 :             nonbonded%pot(start + n_items)%pot%set(1)%deepmd%atom_deepmd_type = atm_deepmd_types(isec)
     705           2 :             nonbonded%pot(start + n_items)%pot%rcutsq = 0.0_dp
     706           4 :             n_items = n_items + 1
     707             :          END DO
     708             :       END DO
     709           2 :    END SUBROUTINE read_deepmd_section
     710             : 
     711             : ! **************************************************************************************************
     712             : !> \brief Reads the QUIP section
     713             : !> \param nonbonded ...
     714             : !> \param section ...
     715             : !> \param start ...
     716             : !> \author teo
     717             : ! **************************************************************************************************
     718           2 :    SUBROUTINE read_quip_section(nonbonded, section, start)
     719             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
     720             :       TYPE(section_vals_type), POINTER                   :: section
     721             :       INTEGER, INTENT(IN)                                :: start
     722             : 
     723             :       CHARACTER(LEN=default_string_length), &
     724           2 :          DIMENSION(:), POINTER                           :: args_str, atm_names
     725             :       INTEGER                                            :: is, isec, n_calc_args, n_items
     726             : 
     727           2 :       CALL section_vals_get(section, n_repetition=n_items)
     728           4 :       DO isec = 1, n_items
     729           2 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
     730             : 
     731           4 :          nonbonded%pot(start + isec)%pot%type = quip_type
     732           2 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
     733           2 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
     734           2 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
     735           2 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
     736             :          CALL section_vals_val_get(section, "PARM_FILE_NAME", i_rep_section=isec, &
     737           2 :                                    c_val=nonbonded%pot(start + isec)%pot%set(1)%quip%quip_file_name)
     738             :          CALL section_vals_val_get(section, "INIT_ARGS", i_rep_section=isec, &
     739           2 :                                    c_vals=args_str)
     740           2 :          nonbonded%pot(start + isec)%pot%set(1)%quip%init_args = ""
     741           6 :          DO is = 1, SIZE(args_str)
     742             :             nonbonded%pot(start + isec)%pot%set(1)%quip%init_args = &
     743             :                TRIM(nonbonded%pot(start + isec)%pot%set(1)%quip%init_args)// &
     744           6 :                " "//TRIM(args_str(is))
     745             :          END DO ! is
     746             :          CALL section_vals_val_get(section, "CALC_ARGS", i_rep_section=isec, &
     747           2 :                                    n_rep_val=n_calc_args)
     748           2 :          IF (n_calc_args > 0) THEN
     749             :             CALL section_vals_val_get(section, "CALC_ARGS", i_rep_section=isec, &
     750           2 :                                       c_vals=args_str)
     751           4 :             DO is = 1, SIZE(args_str)
     752             :                nonbonded%pot(start + isec)%pot%set(1)%quip%calc_args = &
     753             :                   TRIM(nonbonded%pot(start + isec)%pot%set(1)%quip%calc_args)// &
     754           4 :                   " "//TRIM(args_str(is))
     755             :             END DO ! is
     756             :          END IF
     757           6 :          nonbonded%pot(start + isec)%pot%rcutsq = 0.0_dp
     758             :       END DO
     759           2 :    END SUBROUTINE read_quip_section
     760             : 
     761             : ! **************************************************************************************************
     762             : !> \brief Reads the NEQUIP section
     763             : !> \param nonbonded ...
     764             : !> \param section ...
     765             : !> \param start ...
     766             : !> \author Gabriele Tocci
     767             : ! **************************************************************************************************
     768           4 :    SUBROUTINE read_nequip_section(nonbonded, section, start)
     769             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
     770             :       TYPE(section_vals_type), POINTER                   :: section
     771             :       INTEGER, INTENT(IN)                                :: start
     772             : 
     773             :       CHARACTER(LEN=default_string_length)               :: nequip_file_name, unit_cell, &
     774             :                                                             unit_coords, unit_energy, unit_forces
     775             :       CHARACTER(LEN=default_string_length), &
     776           4 :          DIMENSION(:), POINTER                           :: atm_names
     777             :       INTEGER                                            :: isec, jsec, n_items
     778             : 
     779             :       n_items = 1
     780           4 :       isec = 1
     781           4 :       n_items = isec*n_items
     782           4 :       CALL section_vals_val_get(section, "ATOMS", c_vals=atm_names)
     783           4 :       CALL section_vals_val_get(section, "PARM_FILE_NAME", c_val=nequip_file_name)
     784           4 :       CALL section_vals_val_get(section, "UNIT_COORDS", c_val=unit_coords)
     785           4 :       CALL section_vals_val_get(section, "UNIT_ENERGY", c_val=unit_energy)
     786           4 :       CALL section_vals_val_get(section, "UNIT_FORCES", c_val=unit_forces)
     787           4 :       CALL section_vals_val_get(section, "UNIT_CELL", c_val=unit_cell)
     788             : 
     789          12 :       DO isec = 1, SIZE(atm_names)
     790          24 :          DO jsec = isec, SIZE(atm_names)
     791          24 :             nonbonded%pot(start + n_items)%pot%type = nequip_type
     792          12 :             nonbonded%pot(start + n_items)%pot%at1 = atm_names(isec)
     793          12 :             nonbonded%pot(start + n_items)%pot%at2 = atm_names(jsec)
     794          12 :             CALL uppercase(nonbonded%pot(start + n_items)%pot%at1)
     795          12 :             CALL uppercase(nonbonded%pot(start + n_items)%pot%at2)
     796          12 :             nonbonded%pot(start + n_items)%pot%set(1)%nequip%nequip_file_name = discover_file(nequip_file_name)
     797          12 :             nonbonded%pot(start + n_items)%pot%set(1)%nequip%unit_coords = unit_coords
     798          12 :             nonbonded%pot(start + n_items)%pot%set(1)%nequip%unit_forces = unit_forces
     799          12 :             nonbonded%pot(start + n_items)%pot%set(1)%nequip%unit_energy = unit_energy
     800          12 :             nonbonded%pot(start + n_items)%pot%set(1)%nequip%unit_cell = unit_cell
     801          12 :             CALL read_nequip_data(nonbonded%pot(start + n_items)%pot%set(1)%nequip)
     802          12 :             CALL check_cp2k_atom_names_in_torch(atm_names, nonbonded%pot(start + n_items)%pot%set(1)%nequip%type_names_torch)
     803          12 :             nonbonded%pot(start + n_items)%pot%rcutsq = nonbonded%pot(start + n_items)%pot%set(1)%nequip%rcutsq
     804          20 :             n_items = n_items + 1
     805             :          END DO
     806             :       END DO
     807             : 
     808           4 :    END SUBROUTINE read_nequip_section
     809             : 
     810             : ! **************************************************************************************************
     811             : !> \brief Reads the ALLEGRO section
     812             : !> \param nonbonded ...
     813             : !> \param section ...
     814             : !> \param start ...
     815             : !> \author Gabriele Tocci
     816             : ! **************************************************************************************************
     817           4 :    SUBROUTINE read_allegro_section(nonbonded, section, start)
     818             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
     819             :       TYPE(section_vals_type), POINTER                   :: section
     820             :       INTEGER, INTENT(IN)                                :: start
     821             : 
     822             :       CHARACTER(LEN=default_string_length)               :: allegro_file_name, unit_cell, &
     823             :                                                             unit_coords, unit_energy, unit_forces
     824             :       CHARACTER(LEN=default_string_length), &
     825           4 :          DIMENSION(:), POINTER                           :: atm_names
     826             :       INTEGER                                            :: isec, jsec, n_items
     827             : 
     828             :       n_items = 1
     829           4 :       isec = 1
     830           4 :       n_items = isec*n_items
     831           4 :       CALL section_vals_val_get(section, "ATOMS", c_vals=atm_names)
     832           4 :       CALL section_vals_val_get(section, "PARM_FILE_NAME", c_val=allegro_file_name)
     833           4 :       CALL section_vals_val_get(section, "UNIT_COORDS", c_val=unit_coords)
     834           4 :       CALL section_vals_val_get(section, "UNIT_ENERGY", c_val=unit_energy)
     835           4 :       CALL section_vals_val_get(section, "UNIT_FORCES", c_val=unit_forces)
     836           4 :       CALL section_vals_val_get(section, "UNIT_CELL", c_val=unit_cell)
     837             : 
     838          10 :       DO isec = 1, SIZE(atm_names)
     839          18 :          DO jsec = isec, SIZE(atm_names)
     840          16 :             nonbonded%pot(start + n_items)%pot%type = allegro_type
     841           8 :             nonbonded%pot(start + n_items)%pot%at1 = atm_names(isec)
     842           8 :             nonbonded%pot(start + n_items)%pot%at2 = atm_names(jsec)
     843           8 :             CALL uppercase(nonbonded%pot(start + n_items)%pot%at1)
     844           8 :             CALL uppercase(nonbonded%pot(start + n_items)%pot%at2)
     845           8 :             nonbonded%pot(start + n_items)%pot%set(1)%allegro%allegro_file_name = discover_file(allegro_file_name)
     846           8 :             nonbonded%pot(start + n_items)%pot%set(1)%allegro%unit_coords = unit_coords
     847           8 :             nonbonded%pot(start + n_items)%pot%set(1)%allegro%unit_forces = unit_forces
     848           8 :             nonbonded%pot(start + n_items)%pot%set(1)%allegro%unit_energy = unit_energy
     849           8 :             nonbonded%pot(start + n_items)%pot%set(1)%allegro%unit_cell = unit_cell
     850           8 :             CALL read_allegro_data(nonbonded%pot(start + n_items)%pot%set(1)%allegro)
     851           8 :             CALL check_cp2k_atom_names_in_torch(atm_names, nonbonded%pot(start + n_items)%pot%set(1)%allegro%type_names_torch)
     852           8 :             nonbonded%pot(start + n_items)%pot%rcutsq = nonbonded%pot(start + n_items)%pot%set(1)%allegro%rcutsq
     853          14 :             n_items = n_items + 1
     854             :          END DO
     855             :       END DO
     856           4 :    END SUBROUTINE read_allegro_section
     857             : 
     858             : ! **************************************************************************************************
     859             : !> \brief Reads the LJ section
     860             : !> \param nonbonded ...
     861             : !> \param section ...
     862             : !> \param start ...
     863             : !> \author teo
     864             : ! **************************************************************************************************
     865        1006 :    SUBROUTINE read_lj_section(nonbonded, section, start)
     866             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
     867             :       TYPE(section_vals_type), POINTER                   :: section
     868             :       INTEGER, INTENT(IN)                                :: start
     869             : 
     870             :       CHARACTER(LEN=default_string_length), &
     871        1006 :          DIMENSION(:), POINTER                           :: atm_names
     872             :       INTEGER                                            :: isec, n_items, n_rep
     873             :       REAL(KIND=dp)                                      :: epsilon, rcut, sigma
     874             : 
     875        1006 :       CALL section_vals_get(section, n_repetition=n_items)
     876        4794 :       DO isec = 1, n_items
     877        3788 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
     878        3788 :          CALL section_vals_val_get(section, "EPSILON", i_rep_section=isec, r_val=epsilon)
     879        3788 :          CALL section_vals_val_get(section, "SIGMA", i_rep_section=isec, r_val=sigma)
     880        3788 :          CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
     881             : 
     882        7576 :          nonbonded%pot(start + isec)%pot%type = lj_charmm_type
     883        3788 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
     884        3788 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
     885        3788 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
     886        3788 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
     887        3788 :          nonbonded%pot(start + isec)%pot%set(1)%lj%epsilon = epsilon
     888        3788 :          nonbonded%pot(start + isec)%pot%set(1)%lj%sigma6 = sigma**6
     889        3788 :          nonbonded%pot(start + isec)%pot%set(1)%lj%sigma12 = sigma**12
     890        3788 :          nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut
     891             :          !
     892        3788 :          CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
     893        3788 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
     894           2 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
     895        3788 :          CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
     896        3788 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
     897       12372 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
     898             :       END DO
     899        1006 :    END SUBROUTINE read_lj_section
     900             : 
     901             : ! **************************************************************************************************
     902             : !> \brief Reads the WILLIAMS section
     903             : !> \param nonbonded ...
     904             : !> \param section ...
     905             : !> \param start ...
     906             : !> \author teo
     907             : ! **************************************************************************************************
     908         375 :    SUBROUTINE read_wl_section(nonbonded, section, start)
     909             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
     910             :       TYPE(section_vals_type), POINTER                   :: section
     911             :       INTEGER, INTENT(IN)                                :: start
     912             : 
     913             :       CHARACTER(LEN=default_string_length), &
     914         375 :          DIMENSION(:), POINTER                           :: atm_names
     915             :       INTEGER                                            :: isec, n_items, n_rep
     916             :       REAL(KIND=dp)                                      :: a, b, c, rcut
     917             : 
     918         375 :       CALL section_vals_get(section, n_repetition=n_items)
     919        1386 :       DO isec = 1, n_items
     920        1011 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
     921        1011 :          CALL section_vals_val_get(section, "A", i_rep_section=isec, r_val=a)
     922        1011 :          CALL section_vals_val_get(section, "B", i_rep_section=isec, r_val=b)
     923        1011 :          CALL section_vals_val_get(section, "C", i_rep_section=isec, r_val=c)
     924        1011 :          CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
     925             : 
     926        2022 :          nonbonded%pot(start + isec)%pot%type = wl_type
     927        1011 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
     928        1011 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
     929        1011 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
     930        1011 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
     931        1011 :          nonbonded%pot(start + isec)%pot%set(1)%willis%a = a
     932        1011 :          nonbonded%pot(start + isec)%pot%set(1)%willis%b = b
     933        1011 :          nonbonded%pot(start + isec)%pot%set(1)%willis%c = c
     934        1011 :          nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut
     935             :          !
     936        1011 :          CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
     937        1011 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
     938           0 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
     939        1011 :          CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
     940        1011 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
     941        3408 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
     942             :       END DO
     943         375 :    END SUBROUTINE read_wl_section
     944             : 
     945             : ! **************************************************************************************************
     946             : !> \brief Reads the GOODWIN section
     947             : !> \param nonbonded ...
     948             : !> \param section ...
     949             : !> \param start ...
     950             : !> \author teo
     951             : ! **************************************************************************************************
     952           0 :    SUBROUTINE read_gd_section(nonbonded, section, start)
     953             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
     954             :       TYPE(section_vals_type), POINTER                   :: section
     955             :       INTEGER, INTENT(IN)                                :: start
     956             : 
     957             :       CHARACTER(LEN=default_string_length), &
     958           0 :          DIMENSION(:), POINTER                           :: atm_names
     959             :       INTEGER                                            :: isec, m, mc, n_items, n_rep
     960             :       REAL(KIND=dp)                                      :: d, dc, rcut, vr0
     961             : 
     962           0 :       CALL section_vals_get(section, n_repetition=n_items)
     963           0 :       DO isec = 1, n_items
     964           0 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
     965           0 :          CALL section_vals_val_get(section, "VR0", i_rep_section=isec, r_val=vr0)
     966           0 :          CALL section_vals_val_get(section, "D", i_rep_section=isec, r_val=d)
     967           0 :          CALL section_vals_val_get(section, "DC", i_rep_section=isec, r_val=dc)
     968           0 :          CALL section_vals_val_get(section, "M", i_rep_section=isec, i_val=m)
     969           0 :          CALL section_vals_val_get(section, "MC", i_rep_section=isec, i_val=mc)
     970           0 :          CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
     971             : 
     972           0 :          nonbonded%pot(start + isec)%pot%type = gw_type
     973           0 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
     974           0 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
     975           0 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
     976           0 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
     977           0 :          nonbonded%pot(start + isec)%pot%set(1)%goodwin%vr0 = vr0
     978           0 :          nonbonded%pot(start + isec)%pot%set(1)%goodwin%d = d
     979           0 :          nonbonded%pot(start + isec)%pot%set(1)%goodwin%dc = dc
     980           0 :          nonbonded%pot(start + isec)%pot%set(1)%goodwin%m = m
     981           0 :          nonbonded%pot(start + isec)%pot%set(1)%goodwin%mc = mc
     982           0 :          nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut
     983             :          !
     984           0 :          CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
     985           0 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
     986           0 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
     987           0 :          CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
     988           0 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
     989           0 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
     990             :       END DO
     991           0 :    END SUBROUTINE read_gd_section
     992             : 
     993             : ! **************************************************************************************************
     994             : !> \brief Reads the IPBV section
     995             : !> \param nonbonded ...
     996             : !> \param section ...
     997             : !> \param start ...
     998             : !> \author teo
     999             : ! **************************************************************************************************
    1000          16 :    SUBROUTINE read_ipbv_section(nonbonded, section, start)
    1001             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
    1002             :       TYPE(section_vals_type), POINTER                   :: section
    1003             :       INTEGER, INTENT(IN)                                :: start
    1004             : 
    1005             :       CHARACTER(LEN=default_string_length), &
    1006          16 :          DIMENSION(:), POINTER                           :: atm_names
    1007             :       INTEGER                                            :: isec, n_items, n_rep
    1008             :       REAL(KIND=dp)                                      :: rcut
    1009             : 
    1010          16 :       CALL section_vals_get(section, n_repetition=n_items)
    1011          64 :       DO isec = 1, n_items
    1012          48 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    1013          96 :          nonbonded%pot(start + isec)%pot%type = ip_type
    1014          48 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
    1015          48 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
    1016          48 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
    1017          48 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
    1018             :          CALL set_IPBV_ff(nonbonded%pot(start + isec)%pot%at1, nonbonded%pot(start + isec)%pot%at2, &
    1019          48 :                           nonbonded%pot(start + isec)%pot%set(1)%ipbv)
    1020          48 :          CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
    1021          48 :          nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
    1022             :          !
    1023          48 :          CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
    1024          48 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
    1025           0 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
    1026          48 :          CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
    1027          48 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
    1028         112 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
    1029             :       END DO
    1030          16 :    END SUBROUTINE read_ipbv_section
    1031             : 
    1032             : ! **************************************************************************************************
    1033             : !> \brief Reads the BMHFT section
    1034             : !> \param nonbonded ...
    1035             : !> \param section ...
    1036             : !> \param start ...
    1037             : !> \author teo
    1038             : ! **************************************************************************************************
    1039           4 :    SUBROUTINE read_bmhft_section(nonbonded, section, start)
    1040             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
    1041             :       TYPE(section_vals_type), POINTER                   :: section
    1042             :       INTEGER, INTENT(IN)                                :: start
    1043             : 
    1044             :       CHARACTER(LEN=default_string_length), DIMENSION(2) :: map_atoms
    1045             :       CHARACTER(LEN=default_string_length), &
    1046           4 :          DIMENSION(:), POINTER                           :: atm_names
    1047             :       INTEGER                                            :: i, isec, n_items, n_rep
    1048             :       REAL(KIND=dp)                                      :: rcut
    1049             : 
    1050           4 :       CALL section_vals_get(section, n_repetition=n_items)
    1051          16 :       DO isec = 1, n_items
    1052          12 :          CALL cite_reference(Tosi1964a)
    1053          12 :          CALL cite_reference(Tosi1964b)
    1054          12 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    1055          24 :          nonbonded%pot(start + isec)%pot%type = ft_type
    1056          12 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
    1057          12 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
    1058          12 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
    1059          12 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
    1060             : 
    1061          12 :          CALL section_vals_val_get(section, "A", i_rep_section=isec, n_rep_val=i)
    1062          12 :          IF (i == 1) THEN
    1063             :             CALL section_vals_val_get(section, "A", i_rep_section=isec, &
    1064           0 :                                       r_val=nonbonded%pot(start + isec)%pot%set(1)%ft%a)
    1065             :             CALL section_vals_val_get(section, "B", i_rep_section=isec, &
    1066           0 :                                       r_val=nonbonded%pot(start + isec)%pot%set(1)%ft%b)
    1067             :             CALL section_vals_val_get(section, "C", i_rep_section=isec, &
    1068           0 :                                       r_val=nonbonded%pot(start + isec)%pot%set(1)%ft%c)
    1069             :             CALL section_vals_val_get(section, "D", i_rep_section=isec, &
    1070           0 :                                       r_val=nonbonded%pot(start + isec)%pot%set(1)%ft%d)
    1071             :          ELSE
    1072          12 :             CALL section_vals_val_get(section, "MAP_ATOMS", i_rep_section=isec, c_vals=atm_names)
    1073          36 :             map_atoms = atm_names
    1074          12 :             CALL uppercase(map_atoms(1))
    1075          12 :             CALL uppercase(map_atoms(2))
    1076          12 :             CALL set_BMHFT_ff(map_atoms(1), map_atoms(2), nonbonded%pot(start + isec)%pot%set(1)%ft)
    1077             :          END IF
    1078          12 :          CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
    1079          12 :          nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
    1080             :          !
    1081          12 :          CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
    1082          12 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
    1083           0 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
    1084          12 :          CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
    1085          12 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
    1086          40 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
    1087             :       END DO
    1088           4 :    END SUBROUTINE read_bmhft_section
    1089             : 
    1090             : ! **************************************************************************************************
    1091             : !> \brief Reads the BMHFTD section
    1092             : !> \param nonbonded ...
    1093             : !> \param section ...
    1094             : !> \param start ...
    1095             : !> \author Mathieu Salanne 05.2010
    1096             : ! **************************************************************************************************
    1097          18 :    SUBROUTINE read_bmhftd_section(nonbonded, section, start)
    1098             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
    1099             :       TYPE(section_vals_type), POINTER                   :: section
    1100             :       INTEGER, INTENT(IN)                                :: start
    1101             : 
    1102             :       CHARACTER(LEN=default_string_length), DIMENSION(2) :: map_atoms
    1103             :       CHARACTER(LEN=default_string_length), &
    1104          18 :          DIMENSION(:), POINTER                           :: atm_names
    1105             :       INTEGER                                            :: i, isec, n_items, n_rep
    1106             :       REAL(KIND=dp)                                      :: rcut
    1107          18 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: bd_vals
    1108             : 
    1109          18 :       NULLIFY (bd_vals)
    1110             : 
    1111          18 :       CALL section_vals_get(section, n_repetition=n_items)
    1112             : 
    1113          84 :       DO isec = 1, n_items
    1114          66 :          CALL cite_reference(Tosi1964a)
    1115          66 :          CALL cite_reference(Tosi1964b)
    1116          66 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    1117         132 :          nonbonded%pot(start + isec)%pot%type = ftd_type
    1118          66 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
    1119          66 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
    1120          66 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
    1121          66 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
    1122             : 
    1123          66 :          CALL section_vals_val_get(section, "A", i_rep_section=isec, n_rep_val=i)
    1124          66 :          IF (i == 1) THEN
    1125             :             CALL section_vals_val_get(section, "A", i_rep_section=isec, &
    1126          66 :                                       r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%a)
    1127             :             CALL section_vals_val_get(section, "B", i_rep_section=isec, &
    1128          66 :                                       r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%b)
    1129             :             CALL section_vals_val_get(section, "C", i_rep_section=isec, &
    1130          66 :                                       r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%c)
    1131             :             CALL section_vals_val_get(section, "D", i_rep_section=isec, &
    1132          66 :                                       r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%d)
    1133          66 :             CALL section_vals_val_get(section, "BD", i_rep_section=isec, r_vals=bd_vals)
    1134          66 :             IF (ASSOCIATED(bd_vals)) THEN
    1135          66 :                SELECT CASE (SIZE(bd_vals))
    1136             :                CASE (0)
    1137           0 :                   CPABORT("No values specified for parameter BD in section &BMHFTD")
    1138             :                CASE (1)
    1139         186 :                   nonbonded%pot(start + isec)%pot%set(1)%ftd%bd(1:2) = bd_vals(1)
    1140             :                CASE (2)
    1141          24 :                   nonbonded%pot(start + isec)%pot%set(1)%ftd%bd(1:2) = bd_vals(1:2)
    1142             :                CASE DEFAULT
    1143          66 :                   CPABORT("Too many values specified for parameter BD in section &BMHFTD")
    1144             :                END SELECT
    1145             :             ELSE
    1146           0 :                CPABORT("Parameter BD in section &BMHFTD was not specified")
    1147             :             END IF
    1148             :          ELSE
    1149           0 :             CALL section_vals_val_get(section, "MAP_ATOMS", i_rep_section=isec, c_vals=atm_names)
    1150           0 :             map_atoms = atm_names
    1151           0 :             CALL uppercase(map_atoms(1))
    1152           0 :             CALL uppercase(map_atoms(2))
    1153           0 :             CALL set_BMHFTD_ff()
    1154             :          END IF
    1155          66 :          CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
    1156          66 :          nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
    1157             :          !
    1158          66 :          CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
    1159          66 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
    1160           0 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
    1161          66 :          CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
    1162          66 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
    1163         216 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
    1164             :       END DO
    1165          18 :    END SUBROUTINE read_bmhftd_section
    1166             : 
    1167             : ! **************************************************************************************************
    1168             : !> \brief Reads the Buckingham 4 Ranges potential section
    1169             : !> \param nonbonded ...
    1170             : !> \param section ...
    1171             : !> \param start ...
    1172             : !> \par History
    1173             : !>      MK (11.11.2010): Automatic fit of the (default) polynomial coefficients
    1174             : !> \author MI,MK
    1175             : ! **************************************************************************************************
    1176         262 :    SUBROUTINE read_b4_section(nonbonded, section, start)
    1177             : 
    1178             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
    1179             :       TYPE(section_vals_type), POINTER                   :: section
    1180             :       INTEGER, INTENT(IN)                                :: start
    1181             : 
    1182             :       CHARACTER(LEN=default_string_length), &
    1183         262 :          DIMENSION(:), POINTER                           :: atm_names
    1184             :       INTEGER                                            :: i, ir, isec, n_items, n_rep, np1, np2
    1185             :       LOGICAL                                            :: explicit_poly1, explicit_poly2
    1186             :       REAL(KIND=dp)                                      :: a, b, c, eval_error, r1, r2, r3, rcut
    1187             :       REAL(KIND=dp), DIMENSION(10)                       :: v, x
    1188             :       REAL(KIND=dp), DIMENSION(10, 10)                   :: p, p_inv
    1189         262 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: coeff1, coeff2, list
    1190             : 
    1191         262 :       NULLIFY (coeff1)
    1192         262 :       NULLIFY (coeff2)
    1193             : 
    1194         262 :       CALL section_vals_get(section, n_repetition=n_items)
    1195             : 
    1196         524 :       DO isec = 1, n_items
    1197         262 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    1198         262 :          CALL section_vals_val_get(section, "A", i_rep_section=isec, r_val=a)
    1199         262 :          CALL section_vals_val_get(section, "B", i_rep_section=isec, r_val=b)
    1200         262 :          CALL section_vals_val_get(section, "C", i_rep_section=isec, r_val=c)
    1201         262 :          CALL section_vals_val_get(section, "R1", i_rep_section=isec, r_val=r1)
    1202         262 :          CALL section_vals_val_get(section, "R2", i_rep_section=isec, r_val=r2)
    1203         262 :          CALL section_vals_val_get(section, "R3", i_rep_section=isec, r_val=r3)
    1204         262 :          CALL section_vals_val_get(section, "POLY1", explicit=explicit_poly1, n_rep_val=n_rep)
    1205             :          ! Check if polynomial coefficients were specified for range 2 and 3 explicitly
    1206         262 :          IF (explicit_poly1) THEN
    1207          84 :             np1 = 0
    1208         168 :             DO ir = 1, n_rep
    1209          84 :                NULLIFY (list)
    1210          84 :                CALL section_vals_val_get(section, "POLY1", i_rep_val=ir, r_vals=list)
    1211         168 :                IF (ASSOCIATED(list)) THEN
    1212          84 :                   CALL reallocate(coeff1, 0, np1 + SIZE(list) - 1)
    1213         588 :                   DO i = 1, SIZE(list)
    1214         588 :                      coeff1(i + np1 - 1) = list(i)
    1215             :                   END DO
    1216          84 :                   np1 = np1 + SIZE(list)
    1217             :                END IF
    1218             :             END DO
    1219             :          END IF
    1220         262 :          CALL section_vals_val_get(section, "POLY2", explicit=explicit_poly2, n_rep_val=n_rep)
    1221         262 :          IF (explicit_poly2) THEN
    1222          84 :             np2 = 0
    1223         168 :             DO ir = 1, n_rep
    1224          84 :                NULLIFY (list)
    1225          84 :                CALL section_vals_val_get(section, "POLY2", i_rep_val=ir, r_vals=list)
    1226         168 :                IF (ASSOCIATED(list)) THEN
    1227          84 :                   CALL reallocate(coeff2, 0, np2 + SIZE(list) - 1)
    1228         420 :                   DO i = 1, SIZE(list)
    1229         420 :                      coeff2(i + np2 - 1) = list(i)
    1230             :                   END DO
    1231          84 :                   np2 = np2 + SIZE(list)
    1232             :                END IF
    1233             :             END DO
    1234             :          END IF
    1235             :          ! Default is a 5th/3rd-order polynomial fit
    1236         262 :          IF ((.NOT. explicit_poly1) .OR. (.NOT. explicit_poly2)) THEN
    1237             :             ! Build matrix p and vector v to calculate the polynomial coefficients
    1238             :             ! in the vector x from p*x = v
    1239         178 :             p(:, :) = 0.0_dp
    1240             :             ! Row 1: Match the 5th-order polynomial and the potential at r1
    1241         178 :             p(1, 1) = 1.0_dp
    1242        1068 :             DO i = 2, 6
    1243        1068 :                p(1, i) = p(1, i - 1)*r1
    1244             :             END DO
    1245             :             ! Row 2: Match the first derivatives of the 5th-order polynomial and the potential at r1
    1246        1068 :             DO i = 2, 6
    1247        1068 :                p(2, i) = REAL(i - 1, KIND=dp)*p(1, i - 1)
    1248             :             END DO
    1249             :             ! Row 3: Match the second derivatives of the 5th-order polynomial and the potential at r1
    1250         890 :             DO i = 3, 6
    1251         890 :                p(3, i) = REAL(i - 1, KIND=dp)*p(2, i - 1)
    1252             :             END DO
    1253             :             ! Row 4: Match the 5th-order and the 3rd-order polynomials at r2
    1254         178 :             p(4, 1) = 1.0_dp
    1255        1068 :             DO i = 2, 6
    1256        1068 :                p(4, i) = p(4, i - 1)*r2
    1257             :             END DO
    1258         178 :             p(4, 7) = -1.0_dp
    1259         712 :             DO i = 8, 10
    1260         712 :                p(4, i) = p(4, i - 1)*r2
    1261             :             END DO
    1262             :             ! Row 5: Match the first derivatives of the 5th-order and the 3rd-order polynomials at r2
    1263        1068 :             DO i = 2, 6
    1264        1068 :                p(5, i) = REAL(i - 1, KIND=dp)*p(4, i - 1)
    1265             :             END DO
    1266         712 :             DO i = 8, 10
    1267         712 :                p(5, i) = REAL(i - 7, KIND=dp)*p(4, i - 1)
    1268             :             END DO
    1269             :             ! Row 6: Match the second derivatives of the 5th-order and the 3rd-order polynomials at r2
    1270         890 :             DO i = 3, 6
    1271         890 :                p(6, i) = REAL(i - 1, KIND=dp)*p(5, i - 1)
    1272             :             END DO
    1273         534 :             DO i = 9, 10
    1274         534 :                p(6, i) = REAL(i - 7, KIND=dp)*p(5, i - 1)
    1275             :             END DO
    1276             :             ! Row 7: Minimum at r2, ie. the first derivative of the 3rd-order polynomial has to be zero at r2
    1277         712 :             DO i = 8, 10
    1278         712 :                p(7, i) = -p(5, i)
    1279             :             END DO
    1280             :             ! Row 8: Match the 3rd-order polynomial and the potential at r3
    1281         178 :             p(8, 7) = 1.0_dp
    1282         712 :             DO i = 8, 10
    1283         712 :                p(8, i) = p(8, i - 1)*r3
    1284             :             END DO
    1285             :             ! Row 9: Match the first derivatives of the 3rd-order polynomial and the potential at r3
    1286         712 :             DO i = 8, 10
    1287         712 :                p(9, i) = REAL(i - 7, KIND=dp)*p(8, i - 1)
    1288             :             END DO
    1289             :             ! Row 10: Match the second derivatives of the 3rd-order polynomial and the potential at r3
    1290         534 :             DO i = 9, 10
    1291         534 :                p(10, i) = REAL(i - 7, KIND=dp)*p(9, i - 1)
    1292             :             END DO
    1293             :             ! Build the vector v
    1294         178 :             v(1) = a*EXP(-b*r1)
    1295         178 :             v(2) = -b*v(1)
    1296         178 :             v(3) = -b*v(2)
    1297         890 :             v(4:7) = 0.0_dp
    1298         178 :             v(8) = -c/p(8, 10)**2 ! = -c/r3**6
    1299         178 :             v(9) = -6.0_dp*v(8)/r3
    1300         178 :             v(10) = -7.0_dp*v(9)/r3
    1301             :             ! Calculate p_inv the inverse of the matrix p
    1302         178 :             p_inv(:, :) = 0.0_dp
    1303         178 :             CALL invert_matrix(p, p_inv, eval_error)
    1304         178 :             IF (eval_error >= 1.0E-8_dp) &
    1305             :                CALL cp_warn(__LOCATION__, &
    1306             :                             "The polynomial fit for the BUCK4RANGES potential is only accurate to "// &
    1307           0 :                             TRIM(cp_to_string(eval_error)))
    1308             :             ! Get the 6 coefficients of the 5th-order polynomial -> x(1:6)
    1309             :             ! and the 4 coefficients of the 3rd-order polynomial -> x(7:10)
    1310       19758 :             x(:) = MATMUL(p_inv(:, :), v(:))
    1311             :          ELSE
    1312          84 :             x(:) = 0.0_dp
    1313             :          END IF
    1314             : 
    1315         262 :          CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
    1316             : 
    1317         524 :          nonbonded%pot(start + isec)%pot%type = b4_type
    1318         262 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
    1319         262 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
    1320         262 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
    1321         262 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
    1322         262 :          nonbonded%pot(start + isec)%pot%set(1)%buck4r%a = a
    1323         262 :          nonbonded%pot(start + isec)%pot%set(1)%buck4r%b = b
    1324         262 :          nonbonded%pot(start + isec)%pot%set(1)%buck4r%c = c
    1325         262 :          nonbonded%pot(start + isec)%pot%set(1)%buck4r%r1 = r1
    1326         262 :          nonbonded%pot(start + isec)%pot%set(1)%buck4r%r2 = r2
    1327         262 :          nonbonded%pot(start + isec)%pot%set(1)%buck4r%r3 = r3
    1328         262 :          IF ((.NOT. explicit_poly1) .OR. (.NOT. explicit_poly2)) THEN
    1329         178 :             nonbonded%pot(start + isec)%pot%set(1)%buck4r%npoly1 = 5
    1330        1246 :             nonbonded%pot(start + isec)%pot%set(1)%buck4r%poly1(0:5) = x(1:6)
    1331         178 :             nonbonded%pot(start + isec)%pot%set(1)%buck4r%npoly2 = 3
    1332         890 :             nonbonded%pot(start + isec)%pot%set(1)%buck4r%poly2(0:3) = x(7:10)
    1333             :          ELSE
    1334          84 :             nonbonded%pot(start + isec)%pot%set(1)%buck4r%npoly1 = np1 - 1
    1335          84 :             CPASSERT(np1 - 1 <= 10)
    1336        1092 :             nonbonded%pot(start + isec)%pot%set(1)%buck4r%poly1(0:np1 - 1) = coeff1(0:np1 - 1)
    1337          84 :             nonbonded%pot(start + isec)%pot%set(1)%buck4r%npoly2 = np2 - 1
    1338          84 :             CPASSERT(np2 - 1 <= 10)
    1339         756 :             nonbonded%pot(start + isec)%pot%set(1)%buck4r%poly2(0:np2 - 1) = coeff2(0:np2 - 1)
    1340             :          END IF
    1341         262 :          nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut
    1342             : 
    1343         262 :          IF (ASSOCIATED(coeff1)) THEN
    1344          84 :             DEALLOCATE (coeff1)
    1345             :          END IF
    1346         262 :          IF (ASSOCIATED(coeff2)) THEN
    1347          84 :             DEALLOCATE (coeff2)
    1348             :          END IF
    1349         262 :          CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
    1350         262 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
    1351           0 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
    1352         262 :          CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
    1353         262 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
    1354        1572 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
    1355             :       END DO
    1356             : 
    1357         262 :    END SUBROUTINE read_b4_section
    1358             : 
    1359             : ! **************************************************************************************************
    1360             : !> \brief Reads the GENPOT - generic potential section
    1361             : !> \param nonbonded ...
    1362             : !> \param section ...
    1363             : !> \param start ...
    1364             : !> \author Teodoro Laino - 10.2006
    1365             : ! **************************************************************************************************
    1366         576 :    SUBROUTINE read_gp_section(nonbonded, section, start)
    1367             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
    1368             :       TYPE(section_vals_type), POINTER                   :: section
    1369             :       INTEGER, INTENT(IN)                                :: start
    1370             : 
    1371             :       CHARACTER(LEN=default_string_length), &
    1372         576 :          DIMENSION(:), POINTER                           :: atm_names
    1373             :       INTEGER                                            :: isec, n_items, n_rep
    1374             :       REAL(KIND=dp)                                      :: rcut
    1375             : 
    1376         576 :       CALL section_vals_get(section, n_repetition=n_items)
    1377        3794 :       DO isec = 1, n_items
    1378        3218 :          NULLIFY (atm_names)
    1379        3218 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    1380        3218 :          CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
    1381        6436 :          nonbonded%pot(start + isec)%pot%type = gp_type
    1382        3218 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
    1383        3218 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
    1384        3218 :          nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut
    1385        3218 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
    1386        3218 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
    1387             :          ! Parse the genpot info
    1388             :          CALL get_generic_info(section, "FUNCTION", nonbonded%pot(start + isec)%pot%set(1)%gp%potential, &
    1389             :                                nonbonded%pot(start + isec)%pot%set(1)%gp%parameters, &
    1390             :                                nonbonded%pot(start + isec)%pot%set(1)%gp%values, &
    1391        3218 :                                size_variables=1, i_rep_sec=isec)
    1392        3218 :          nonbonded%pot(start + isec)%pot%set(1)%gp%variables = nonbonded%pot(start + isec)%pot%set(1)%gp%parameters(1)
    1393             :          !
    1394        3218 :          CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
    1395        3218 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
    1396          21 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
    1397        3218 :          CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
    1398        3218 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
    1399       10251 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
    1400             :       END DO
    1401         576 :    END SUBROUTINE read_gp_section
    1402             : 
    1403             : ! **************************************************************************************************
    1404             : !> \brief Reads the tersoff section
    1405             : !> \param nonbonded ...
    1406             : !> \param section ...
    1407             : !> \param start ...
    1408             : !> \param tersoff_section ...
    1409             : !> \author ikuo
    1410             : ! **************************************************************************************************
    1411          40 :    SUBROUTINE read_tersoff_section(nonbonded, section, start, tersoff_section)
    1412             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
    1413             :       TYPE(section_vals_type), POINTER                   :: section
    1414             :       INTEGER, INTENT(IN)                                :: start
    1415             :       TYPE(section_vals_type), POINTER                   :: tersoff_section
    1416             : 
    1417             :       CHARACTER(LEN=default_string_length), &
    1418          40 :          DIMENSION(:), POINTER                           :: atm_names
    1419             :       INTEGER                                            :: isec, n_items, n_rep
    1420             :       REAL(KIND=dp)                                      :: rcut, rcutsq
    1421             : 
    1422          40 :       CALL section_vals_get(section, n_repetition=n_items)
    1423          84 :       DO isec = 1, n_items
    1424          44 :          CALL cite_reference(Tersoff1988)
    1425          44 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    1426             : 
    1427          88 :          nonbonded%pot(start + isec)%pot%type = tersoff_type
    1428          44 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
    1429          44 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
    1430          44 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
    1431          44 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
    1432             : 
    1433             :          CALL section_vals_val_get(tersoff_section, "A", i_rep_section=isec, &
    1434          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%A)
    1435             :          CALL section_vals_val_get(tersoff_section, "B", i_rep_section=isec, &
    1436          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%B)
    1437             :          CALL section_vals_val_get(tersoff_section, "lambda1", i_rep_section=isec, &
    1438          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%lambda1)
    1439             :          CALL section_vals_val_get(tersoff_section, "lambda2", i_rep_section=isec, &
    1440          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%lambda2)
    1441             :          CALL section_vals_val_get(tersoff_section, "alpha", i_rep_section=isec, &
    1442          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%alpha)
    1443             :          CALL section_vals_val_get(tersoff_section, "beta", i_rep_section=isec, &
    1444          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%beta)
    1445             :          CALL section_vals_val_get(tersoff_section, "n", i_rep_section=isec, &
    1446          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%n)
    1447             :          CALL section_vals_val_get(tersoff_section, "c", i_rep_section=isec, &
    1448          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%c)
    1449             :          CALL section_vals_val_get(tersoff_section, "d", i_rep_section=isec, &
    1450          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%d)
    1451             :          CALL section_vals_val_get(tersoff_section, "h", i_rep_section=isec, &
    1452          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%h)
    1453             :          CALL section_vals_val_get(tersoff_section, "lambda3", i_rep_section=isec, &
    1454          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%lambda3)
    1455             :          CALL section_vals_val_get(tersoff_section, "bigR", i_rep_section=isec, &
    1456          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%bigR)
    1457             :          CALL section_vals_val_get(tersoff_section, "bigD", i_rep_section=isec, &
    1458          44 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%bigD)
    1459             : 
    1460             :          rcutsq = (nonbonded%pot(start + isec)%pot%set(1)%tersoff%bigR + &
    1461          44 :                    nonbonded%pot(start + isec)%pot%set(1)%tersoff%bigD)**2
    1462          44 :          nonbonded%pot(start + isec)%pot%set(1)%tersoff%rcutsq = rcutsq
    1463          44 :          nonbonded%pot(start + isec)%pot%rcutsq = rcutsq
    1464             : 
    1465             :          ! In case it is defined override the standard specification of RCUT
    1466          44 :          CALL section_vals_val_get(tersoff_section, "RCUT", i_rep_section=isec, n_rep_val=n_rep)
    1467          84 :          IF (n_rep == 1) THEN
    1468          26 :             CALL section_vals_val_get(tersoff_section, "RCUT", i_rep_section=isec, r_val=rcut)
    1469          26 :             nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
    1470             :          END IF
    1471             :       END DO
    1472          40 :    END SUBROUTINE read_tersoff_section
    1473             : 
    1474             : ! **************************************************************************************************
    1475             : !> \brief Reads the gal19 section
    1476             : !> \param nonbonded ...
    1477             : !> \param section ...
    1478             : !> \param start ...
    1479             : !> \param gal_section ...
    1480             : !> \author Clabaut Paul
    1481             : ! **************************************************************************************************
    1482           1 :    SUBROUTINE read_gal_section(nonbonded, section, start, gal_section)
    1483             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
    1484             :       TYPE(section_vals_type), POINTER                   :: section
    1485             :       INTEGER, INTENT(IN)                                :: start
    1486             :       TYPE(section_vals_type), POINTER                   :: gal_section
    1487             : 
    1488             :       CHARACTER(LEN=default_string_length), &
    1489           1 :          DIMENSION(:), POINTER                           :: atm_names
    1490             :       INTEGER                                            :: iatom, isec, n_items, n_rep, nval
    1491             :       LOGICAL                                            :: is_ok
    1492             :       REAL(KIND=dp)                                      :: rcut, rval
    1493           1 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: rvalues
    1494             :       TYPE(cp_sll_val_type), POINTER                     :: list
    1495             :       TYPE(section_vals_type), POINTER                   :: subsection
    1496             :       TYPE(val_type), POINTER                            :: val
    1497             : 
    1498           1 :       CALL section_vals_get(section, n_repetition=n_items)
    1499           2 :       DO isec = 1, n_items
    1500           1 :          CALL cite_reference(Clabaut2020)
    1501           1 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    1502             : 
    1503           2 :          nonbonded%pot(start + isec)%pot%type = gal_type
    1504           1 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
    1505           1 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
    1506           1 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
    1507           1 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
    1508             : 
    1509           1 :          CALL section_vals_val_get(section, "METALS", i_rep_section=isec, c_vals=atm_names)
    1510           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal%met1 = atm_names(1)
    1511           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal%met2 = atm_names(2)
    1512             : 
    1513             :          CALL section_vals_val_get(gal_section, "epsilon", i_rep_section=isec, &
    1514           1 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%epsilon)
    1515             :          CALL section_vals_val_get(gal_section, "bxy", i_rep_section=isec, &
    1516           1 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%bxy)
    1517             :          CALL section_vals_val_get(gal_section, "bz", i_rep_section=isec, &
    1518           1 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%bz)
    1519             : 
    1520           1 :          CALL section_vals_val_get(gal_section, "r", i_rep_section=isec, r_vals=rvalues)
    1521           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal%r1 = rvalues(1)
    1522           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal%r2 = rvalues(2)
    1523             : 
    1524             :          CALL section_vals_val_get(gal_section, "a1", i_rep_section=isec, &
    1525           1 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a1)
    1526             :          CALL section_vals_val_get(gal_section, "a2", i_rep_section=isec, &
    1527           1 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a2)
    1528             :          CALL section_vals_val_get(gal_section, "a3", i_rep_section=isec, &
    1529           1 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a3)
    1530             :          CALL section_vals_val_get(gal_section, "a4", i_rep_section=isec, &
    1531           1 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a4)
    1532             :          CALL section_vals_val_get(gal_section, "A", i_rep_section=isec, &
    1533           1 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a)
    1534             :          CALL section_vals_val_get(gal_section, "B", i_rep_section=isec, &
    1535           1 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%b)
    1536             :          CALL section_vals_val_get(gal_section, "C", i_rep_section=isec, &
    1537           1 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%c)
    1538           1 :          NULLIFY (list)
    1539           1 :          subsection => section_vals_get_subs_vals(section, "GCN", i_rep_section=isec)
    1540           1 :          CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=nval)
    1541           3 :          ALLOCATE (nonbonded%pot(start + isec)%pot%set(1)%gal%gcn(nval))
    1542           1 :          CALL section_vals_list_get(subsection, "_DEFAULT_KEYWORD_", list=list)
    1543         871 :          DO iatom = 1, nval
    1544             :             ! we use only the first default_string_length characters of each line
    1545         870 :             is_ok = cp_sll_val_next(list, val)
    1546         870 :             CALL val_get(val, r_val=rval)
    1547             :             ! assign values
    1548         871 :             nonbonded%pot(start + isec)%pot%set(1)%gal%gcn(iatom) = rval
    1549             :          END DO
    1550             : 
    1551             :          CALL section_vals_val_get(gal_section, "Fit_express", i_rep_section=isec, &
    1552           1 :                                    l_val=nonbonded%pot(start + isec)%pot%set(1)%gal%express)
    1553             : 
    1554             :          ! ! In case it is defined override the standard specification of RCUT
    1555           1 :          CALL section_vals_val_get(gal_section, "RCUT", i_rep_section=isec, n_rep_val=n_rep)
    1556           3 :          IF (n_rep == 1) THEN
    1557           1 :             CALL section_vals_val_get(gal_section, "RCUT", i_rep_section=isec, r_val=rcut)
    1558           1 :             nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
    1559           1 :             nonbonded%pot(start + isec)%pot%set(1)%gal%rcutsq = rcut**2
    1560             :          END IF
    1561             :       END DO
    1562           1 :    END SUBROUTINE read_gal_section
    1563             : 
    1564             : ! **************************************************************************************************
    1565             : !> \brief Reads the gal21 section
    1566             : !> \param nonbonded ...
    1567             : !> \param section ...
    1568             : !> \param start ...
    1569             : !> \param gal21_section ...
    1570             : !> \author Clabaut Paul
    1571             : ! **************************************************************************************************
    1572           1 :    SUBROUTINE read_gal21_section(nonbonded, section, start, gal21_section)
    1573             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
    1574             :       TYPE(section_vals_type), POINTER                   :: section
    1575             :       INTEGER, INTENT(IN)                                :: start
    1576             :       TYPE(section_vals_type), POINTER                   :: gal21_section
    1577             : 
    1578             :       CHARACTER(LEN=default_string_length), &
    1579           1 :          DIMENSION(:), POINTER                           :: atm_names
    1580             :       INTEGER                                            :: iatom, isec, n_items, n_rep, nval
    1581             :       LOGICAL                                            :: is_ok
    1582             :       REAL(KIND=dp)                                      :: rcut, rval
    1583           1 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: rvalues
    1584             :       TYPE(cp_sll_val_type), POINTER                     :: list
    1585             :       TYPE(section_vals_type), POINTER                   :: subsection
    1586             :       TYPE(val_type), POINTER                            :: val
    1587             : 
    1588           1 :       CALL section_vals_get(section, n_repetition=n_items)
    1589           2 :       DO isec = 1, n_items
    1590           1 :          CALL cite_reference(Clabaut2021)
    1591           1 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    1592             : 
    1593           2 :          nonbonded%pot(start + isec)%pot%type = gal21_type
    1594           1 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
    1595           1 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
    1596           1 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
    1597           1 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
    1598             : 
    1599           1 :          CALL section_vals_val_get(section, "METALS", i_rep_section=isec, c_vals=atm_names)
    1600           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%met1 = atm_names(1)
    1601           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%met2 = atm_names(2)
    1602             : 
    1603           1 :          CALL section_vals_val_get(gal21_section, "epsilon", i_rep_section=isec, r_vals=rvalues)
    1604           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%epsilon1 = rvalues(1)
    1605           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%epsilon2 = rvalues(2)
    1606           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%epsilon3 = rvalues(3)
    1607             : 
    1608           1 :          CALL section_vals_val_get(gal21_section, "bxy", i_rep_section=isec, r_vals=rvalues)
    1609           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%bxy1 = rvalues(1)
    1610           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%bxy2 = rvalues(2)
    1611             : 
    1612           1 :          CALL section_vals_val_get(gal21_section, "bz", i_rep_section=isec, r_vals=rvalues)
    1613           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%bz1 = rvalues(1)
    1614           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%bz2 = rvalues(2)
    1615             : 
    1616           1 :          CALL section_vals_val_get(gal21_section, "r", i_rep_section=isec, r_vals=rvalues)
    1617           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%r1 = rvalues(1)
    1618           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%r2 = rvalues(2)
    1619             : 
    1620           1 :          CALL section_vals_val_get(gal21_section, "a1", i_rep_section=isec, r_vals=rvalues)
    1621           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%a11 = rvalues(1)
    1622           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%a12 = rvalues(2)
    1623           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%a13 = rvalues(3)
    1624             : 
    1625           1 :          CALL section_vals_val_get(gal21_section, "a2", i_rep_section=isec, r_vals=rvalues)
    1626           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%a21 = rvalues(1)
    1627           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%a22 = rvalues(2)
    1628           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%a23 = rvalues(3)
    1629             : 
    1630           1 :          CALL section_vals_val_get(gal21_section, "a3", i_rep_section=isec, r_vals=rvalues)
    1631           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%a31 = rvalues(1)
    1632           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%a32 = rvalues(2)
    1633           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%a33 = rvalues(3)
    1634             : 
    1635           1 :          CALL section_vals_val_get(gal21_section, "a4", i_rep_section=isec, r_vals=rvalues)
    1636           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%a41 = rvalues(1)
    1637           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%a42 = rvalues(2)
    1638           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%a43 = rvalues(3)
    1639             : 
    1640           1 :          CALL section_vals_val_get(gal21_section, "A", i_rep_section=isec, r_vals=rvalues)
    1641           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%AO1 = rvalues(1)
    1642           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%AO2 = rvalues(2)
    1643             : 
    1644           1 :          CALL section_vals_val_get(gal21_section, "B", i_rep_section=isec, r_vals=rvalues)
    1645           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%BO1 = rvalues(1)
    1646           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%BO2 = rvalues(2)
    1647             : 
    1648             :          CALL section_vals_val_get(gal21_section, "C", i_rep_section=isec, &
    1649           1 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%gal21%c)
    1650             : 
    1651           1 :          CALL section_vals_val_get(gal21_section, "AH", i_rep_section=isec, r_vals=rvalues)
    1652           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%AH1 = rvalues(1)
    1653           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%AH2 = rvalues(2)
    1654             : 
    1655           1 :          CALL section_vals_val_get(gal21_section, "BH", i_rep_section=isec, r_vals=rvalues)
    1656           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%BH1 = rvalues(1)
    1657           1 :          nonbonded%pot(start + isec)%pot%set(1)%gal21%BH2 = rvalues(2)
    1658             : 
    1659           1 :          NULLIFY (list)
    1660           1 :          subsection => section_vals_get_subs_vals(section, "GCN", i_rep_section=isec)
    1661           1 :          CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=nval)
    1662           3 :          ALLOCATE (nonbonded%pot(start + isec)%pot%set(1)%gal21%gcn(nval))
    1663           1 :          CALL section_vals_list_get(subsection, "_DEFAULT_KEYWORD_", list=list)
    1664         871 :          DO iatom = 1, nval
    1665             :             ! we use only the first default_string_length characters of each line
    1666         870 :             is_ok = cp_sll_val_next(list, val)
    1667         870 :             CALL val_get(val, r_val=rval)
    1668             :             ! assign values
    1669         871 :             nonbonded%pot(start + isec)%pot%set(1)%gal21%gcn(iatom) = rval
    1670             :          END DO
    1671             : 
    1672             :          CALL section_vals_val_get(gal21_section, "Fit_express", i_rep_section=isec, &
    1673           1 :                                    l_val=nonbonded%pot(start + isec)%pot%set(1)%gal21%express)
    1674             : 
    1675             :          ! ! In case it is defined override the standard specification of RCUT
    1676           1 :          CALL section_vals_val_get(gal21_section, "RCUT", i_rep_section=isec, n_rep_val=n_rep)
    1677           3 :          IF (n_rep == 1) THEN
    1678           1 :             CALL section_vals_val_get(gal21_section, "RCUT", i_rep_section=isec, r_val=rcut)
    1679           1 :             nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
    1680           1 :             nonbonded%pot(start + isec)%pot%set(1)%gal21%rcutsq = rcut**2
    1681             :          END IF
    1682             :       END DO
    1683           1 :    END SUBROUTINE read_gal21_section
    1684             : 
    1685             : ! **************************************************************************************************
    1686             : !> \brief Reads the siepmann section
    1687             : !> \param nonbonded ...
    1688             : !> \param section ...
    1689             : !> \param start ...
    1690             : !> \param siepmann_section ...
    1691             : !> \author Dorothea Golze
    1692             : ! **************************************************************************************************
    1693           5 :    SUBROUTINE read_siepmann_section(nonbonded, section, start, siepmann_section)
    1694             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
    1695             :       TYPE(section_vals_type), POINTER                   :: section
    1696             :       INTEGER, INTENT(IN)                                :: start
    1697             :       TYPE(section_vals_type), POINTER                   :: siepmann_section
    1698             : 
    1699             :       CHARACTER(LEN=default_string_length), &
    1700           5 :          DIMENSION(:), POINTER                           :: atm_names
    1701             :       INTEGER                                            :: isec, n_items, n_rep
    1702             :       REAL(KIND=dp)                                      :: rcut
    1703             : 
    1704           5 :       CALL section_vals_get(section, n_repetition=n_items)
    1705          10 :       DO isec = 1, n_items
    1706           5 :          CALL cite_reference(Siepmann1995)
    1707           5 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    1708             : 
    1709          10 :          nonbonded%pot(start + isec)%pot%type = siepmann_type
    1710           5 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
    1711           5 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
    1712           5 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
    1713           5 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
    1714             : 
    1715             :          CALL section_vals_val_get(siepmann_section, "B", i_rep_section=isec, &
    1716           5 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%B)
    1717             :          CALL section_vals_val_get(siepmann_section, "D", i_rep_section=isec, &
    1718           5 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%D)
    1719             :          CALL section_vals_val_get(siepmann_section, "E", i_rep_section=isec, &
    1720           5 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%E)
    1721             :          CALL section_vals_val_get(siepmann_section, "F", i_rep_section=isec, &
    1722           5 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%F)
    1723             :          CALL section_vals_val_get(siepmann_section, "beta", i_rep_section=isec, &
    1724           5 :                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%beta)
    1725             :          CALL section_vals_val_get(siepmann_section, "ALLOW_OH_FORMATION", i_rep_section=isec, &
    1726           5 :                                    l_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%allow_oh_formation)
    1727             :          CALL section_vals_val_get(siepmann_section, "ALLOW_H3O_FORMATION", i_rep_section=isec, &
    1728           5 :                                    l_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%allow_h3o_formation)
    1729             :          CALL section_vals_val_get(siepmann_section, "ALLOW_O_FORMATION", i_rep_section=isec, &
    1730           5 :                                    l_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%allow_o_formation)
    1731             : 
    1732             :          ! ! In case it is defined override the standard specification of RCUT
    1733           5 :          CALL section_vals_val_get(siepmann_section, "RCUT", i_rep_section=isec, n_rep_val=n_rep)
    1734          10 :          IF (n_rep == 1) THEN
    1735           5 :             CALL section_vals_val_get(siepmann_section, "RCUT", i_rep_section=isec, r_val=rcut)
    1736           5 :             nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
    1737           5 :             nonbonded%pot(start + isec)%pot%set(1)%siepmann%rcutsq = rcut**2
    1738             :          END IF
    1739             :       END DO
    1740           5 :    END SUBROUTINE read_siepmann_section
    1741             : 
    1742             : ! **************************************************************************************************
    1743             : !> \brief Reads the Buckingham plus Morse potential section
    1744             : !> \param nonbonded ...
    1745             : !> \param section ...
    1746             : !> \param start ...
    1747             : !> \author MI
    1748             : ! **************************************************************************************************
    1749           8 :    SUBROUTINE read_bm_section(nonbonded, section, start)
    1750             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
    1751             :       TYPE(section_vals_type), POINTER                   :: section
    1752             :       INTEGER, INTENT(IN)                                :: start
    1753             : 
    1754             :       CHARACTER(LEN=default_string_length), &
    1755           8 :          DIMENSION(:), POINTER                           :: atm_names
    1756             :       INTEGER                                            :: isec, n_items, n_rep
    1757             :       REAL(KIND=dp)                                      :: a1, a2, b1, b2, beta, c, d, f0, r0, rcut
    1758             : 
    1759           8 :       CALL section_vals_get(section, n_repetition=n_items)
    1760          26 :       DO isec = 1, n_items
    1761          18 :          CALL cite_reference(Yamada2000)
    1762          18 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    1763          18 :          CALL section_vals_val_get(section, "F0", i_rep_section=isec, r_val=f0)
    1764          18 :          CALL section_vals_val_get(section, "A1", i_rep_section=isec, r_val=a1)
    1765          18 :          CALL section_vals_val_get(section, "A2", i_rep_section=isec, r_val=a2)
    1766          18 :          CALL section_vals_val_get(section, "B1", i_rep_section=isec, r_val=b1)
    1767          18 :          CALL section_vals_val_get(section, "B2", i_rep_section=isec, r_val=b2)
    1768          18 :          CALL section_vals_val_get(section, "C", i_rep_section=isec, r_val=c)
    1769          18 :          CALL section_vals_val_get(section, "D", i_rep_section=isec, r_val=d)
    1770          18 :          CALL section_vals_val_get(section, "R0", i_rep_section=isec, r_val=r0)
    1771          18 :          CALL section_vals_val_get(section, "Beta", i_rep_section=isec, r_val=beta)
    1772          18 :          CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
    1773             : 
    1774          36 :          nonbonded%pot(start + isec)%pot%type = bm_type
    1775          18 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
    1776          18 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
    1777          18 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
    1778          18 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
    1779          18 :          nonbonded%pot(start + isec)%pot%set(1)%buckmo%f0 = f0
    1780          18 :          nonbonded%pot(start + isec)%pot%set(1)%buckmo%a1 = a1
    1781          18 :          nonbonded%pot(start + isec)%pot%set(1)%buckmo%a2 = a2
    1782          18 :          nonbonded%pot(start + isec)%pot%set(1)%buckmo%b1 = b1
    1783          18 :          nonbonded%pot(start + isec)%pot%set(1)%buckmo%b2 = b2
    1784          18 :          nonbonded%pot(start + isec)%pot%set(1)%buckmo%c = c
    1785          18 :          nonbonded%pot(start + isec)%pot%set(1)%buckmo%d = d
    1786          18 :          nonbonded%pot(start + isec)%pot%set(1)%buckmo%r0 = r0
    1787          18 :          nonbonded%pot(start + isec)%pot%set(1)%buckmo%beta = beta
    1788          18 :          nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut
    1789             :          !
    1790          18 :          CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
    1791          18 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
    1792           0 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
    1793          18 :          CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
    1794          18 :          IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
    1795          62 :                                                    r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
    1796             :       END DO
    1797           8 :    END SUBROUTINE read_bm_section
    1798             : 
    1799             : ! **************************************************************************************************
    1800             : !> \brief Reads the TABPOT section
    1801             : !> \param nonbonded ...
    1802             : !> \param section ...
    1803             : !> \param start ...
    1804             : !> \param para_env ...
    1805             : !> \param mm_section ...
    1806             : !> \author Alex Mironenko, Da Teng
    1807             : ! **************************************************************************************************
    1808           8 :    SUBROUTINE read_tabpot_section(nonbonded, section, start, para_env, mm_section)
    1809             :       TYPE(pair_potential_p_type), POINTER               :: nonbonded
    1810             :       TYPE(section_vals_type), POINTER                   :: section
    1811             :       INTEGER, INTENT(IN)                                :: start
    1812             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1813             :       TYPE(section_vals_type), POINTER                   :: mm_section
    1814             : 
    1815             :       CHARACTER(LEN=default_string_length), &
    1816           8 :          DIMENSION(:), POINTER                           :: atm_names
    1817             :       INTEGER                                            :: isec, n_items
    1818             : 
    1819           8 :       CALL section_vals_get(section, n_repetition=n_items)
    1820          32 :       DO isec = 1, n_items
    1821          24 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    1822          48 :          nonbonded%pot(start + isec)%pot%type = tab_type
    1823          24 :          nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
    1824          24 :          nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
    1825          24 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
    1826          24 :          CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
    1827             :          CALL section_vals_val_get(section, "PARM_FILE_NAME", i_rep_section=isec, &
    1828          24 :                                    c_val=nonbonded%pot(start + isec)%pot%set(1)%tab%tabpot_file_name)
    1829          24 :          CALL read_tabpot_data(nonbonded%pot(start + isec)%pot%set(1)%tab, para_env, mm_section)
    1830          32 :          nonbonded%pot(start + isec)%pot%set(1)%tab%index = isec
    1831             :       END DO
    1832           8 :    END SUBROUTINE read_tabpot_section
    1833             : 
    1834             : ! **************************************************************************************************
    1835             : !> \brief Reads the CHARGE section
    1836             : !> \param charge_atm ...
    1837             : !> \param charge ...
    1838             : !> \param section ...
    1839             : !> \param start ...
    1840             : !> \author teo
    1841             : ! **************************************************************************************************
    1842        2109 :    SUBROUTINE read_chrg_section(charge_atm, charge, section, start)
    1843             :       CHARACTER(LEN=default_string_length), &
    1844             :          DIMENSION(:), POINTER                           :: charge_atm
    1845             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: charge
    1846             :       TYPE(section_vals_type), POINTER                   :: section
    1847             :       INTEGER, INTENT(IN)                                :: start
    1848             : 
    1849             :       CHARACTER(LEN=default_string_length)               :: atm_name
    1850             :       INTEGER                                            :: isec, n_items
    1851             : 
    1852        2109 :       CALL section_vals_get(section, n_repetition=n_items)
    1853        7266 :       DO isec = 1, n_items
    1854        5157 :          CALL section_vals_val_get(section, "ATOM", i_rep_section=isec, c_val=atm_name)
    1855        5157 :          charge_atm(start + isec) = atm_name
    1856        5157 :          CALL uppercase(charge_atm(start + isec))
    1857        7266 :          CALL section_vals_val_get(section, "CHARGE", i_rep_section=isec, r_val=charge(start + isec))
    1858             :       END DO
    1859        2109 :    END SUBROUTINE read_chrg_section
    1860             : 
    1861             : ! **************************************************************************************************
    1862             : !> \brief Reads the POLARIZABILITY section
    1863             : !> \param apol_atm ...
    1864             : !> \param apol ...
    1865             : !> \param damping_list ...
    1866             : !> \param section ...
    1867             : !> \param start ...
    1868             : !> \author Marcel Baer
    1869             : ! **************************************************************************************************
    1870          34 :    SUBROUTINE read_apol_section(apol_atm, apol, damping_list, section, &
    1871             :                                 start)
    1872             :       CHARACTER(LEN=default_string_length), &
    1873             :          DIMENSION(:), POINTER                           :: apol_atm
    1874             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: apol
    1875             :       TYPE(damping_info_type), DIMENSION(:), POINTER     :: damping_list
    1876             :       TYPE(section_vals_type), POINTER                   :: section
    1877             :       INTEGER, INTENT(IN)                                :: start
    1878             : 
    1879             :       CHARACTER(LEN=default_string_length)               :: atm_name
    1880             :       INTEGER                                            :: isec, isec_damp, n_damp, n_items, &
    1881             :                                                             start_damp, tmp_damp
    1882             :       TYPE(section_vals_type), POINTER                   :: tmp_section
    1883             : 
    1884          34 :       CALL section_vals_get(section, n_repetition=n_items)
    1885          34 :       NULLIFY (tmp_section)
    1886          34 :       n_damp = 0
    1887             : ! *** Counts number of DIPOLE%DAMPING sections ****
    1888         102 :       DO isec = 1, n_items
    1889             :          tmp_section => section_vals_get_subs_vals(section, "DAMPING", &
    1890          68 :                                                    i_rep_section=isec)
    1891          68 :          CALL section_vals_get(tmp_section, n_repetition=tmp_damp)
    1892         102 :          n_damp = n_damp + tmp_damp
    1893             : 
    1894             :       END DO
    1895             : 
    1896          34 :       IF (n_damp > 0) THEN
    1897          42 :          ALLOCATE (damping_list(1:n_damp))
    1898             :       END IF
    1899             : 
    1900             : ! *** Reads DIPOLE sections *****
    1901          34 :       start_damp = 0
    1902         102 :       DO isec = 1, n_items
    1903          68 :          CALL section_vals_val_get(section, "ATOM", i_rep_section=isec, c_val=atm_name)
    1904          68 :          apol_atm(start + isec) = atm_name
    1905          68 :          CALL uppercase(apol_atm(start + isec))
    1906          68 :          CALL section_vals_val_get(section, "APOL", i_rep_section=isec, r_val=apol(start + isec))
    1907             : 
    1908             :          tmp_section => section_vals_get_subs_vals(section, "DAMPING", &
    1909          68 :                                                    i_rep_section=isec)
    1910          68 :          CALL section_vals_get(tmp_section, n_repetition=tmp_damp)
    1911          80 :          DO isec_damp = 1, tmp_damp
    1912          12 :             damping_list(start_damp + isec_damp)%atm_name1 = apol_atm(start + isec)
    1913             :             CALL section_vals_val_get(tmp_section, "ATOM", i_rep_section=isec_damp, &
    1914          12 :                                       c_val=atm_name)
    1915          12 :             damping_list(start_damp + isec_damp)%atm_name2 = atm_name
    1916          12 :             CALL uppercase(damping_list(start_damp + isec_damp)%atm_name2)
    1917             :             CALL section_vals_val_get(tmp_section, "TYPE", i_rep_section=isec_damp, &
    1918          12 :                                       c_val=atm_name)
    1919          12 :             damping_list(start_damp + isec_damp)%dtype = atm_name
    1920          12 :             CALL uppercase(damping_list(start_damp + isec_damp)%dtype)
    1921             : 
    1922             :             CALL section_vals_val_get(tmp_section, "ORDER", i_rep_section=isec_damp, &
    1923          12 :                                       i_val=damping_list(start_damp + isec_damp)%order)
    1924             :             CALL section_vals_val_get(tmp_section, "BIJ", i_rep_section=isec_damp, &
    1925          12 :                                       r_val=damping_list(start_damp + isec_damp)%bij)
    1926             :             CALL section_vals_val_get(tmp_section, "CIJ", i_rep_section=isec_damp, &
    1927          80 :                                       r_val=damping_list(start_damp + isec_damp)%cij)
    1928             :          END DO
    1929         170 :          start_damp = start_damp + tmp_damp
    1930             : 
    1931             :       END DO
    1932             : 
    1933          34 :    END SUBROUTINE read_apol_section
    1934             : 
    1935             : ! **************************************************************************************************
    1936             : !> \brief Reads the QUADRUPOLE POLARIZABILITY section
    1937             : !> \param cpol_atm ...
    1938             : !> \param cpol ...
    1939             : !> \param section ...
    1940             : !> \param start ...
    1941             : !> \author Marcel Baer
    1942             : ! **************************************************************************************************
    1943           0 :    SUBROUTINE read_cpol_section(cpol_atm, cpol, section, start)
    1944             :       CHARACTER(LEN=default_string_length), &
    1945             :          DIMENSION(:), POINTER                           :: cpol_atm
    1946             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: cpol
    1947             :       TYPE(section_vals_type), POINTER                   :: section
    1948             :       INTEGER, INTENT(IN)                                :: start
    1949             : 
    1950             :       CHARACTER(LEN=default_string_length)               :: atm_name
    1951             :       INTEGER                                            :: isec, n_items
    1952             : 
    1953           0 :       CALL section_vals_get(section, n_repetition=n_items)
    1954           0 :       DO isec = 1, n_items
    1955           0 :          CALL section_vals_val_get(section, "ATOM", i_rep_section=isec, c_val=atm_name)
    1956           0 :          cpol_atm(start + isec) = atm_name
    1957           0 :          CALL uppercase(cpol_atm(start + isec))
    1958           0 :          CALL section_vals_val_get(section, "CPOL", i_rep_section=isec, r_val=cpol(start + isec))
    1959             :       END DO
    1960           0 :    END SUBROUTINE read_cpol_section
    1961             : 
    1962             : ! **************************************************************************************************
    1963             : !> \brief Reads the SHELL section
    1964             : !> \param shell_list ...
    1965             : !> \param section ...
    1966             : !> \param start ...
    1967             : !> \author Marcella Iannuzzi
    1968             : ! **************************************************************************************************
    1969         258 :    SUBROUTINE read_shell_section(shell_list, section, start)
    1970             : 
    1971             :       TYPE(shell_p_type), DIMENSION(:), POINTER          :: shell_list
    1972             :       TYPE(section_vals_type), POINTER                   :: section
    1973             :       INTEGER, INTENT(IN)                                :: start
    1974             : 
    1975             :       CHARACTER(LEN=default_string_length)               :: atm_name
    1976             :       INTEGER                                            :: i_rep, n_rep
    1977             :       REAL(dp)                                           :: ccharge, cutoff, k, maxdist, mfrac, &
    1978             :                                                             scharge
    1979             : 
    1980         258 :       CALL section_vals_get(section, n_repetition=n_rep)
    1981             : 
    1982         710 :       DO i_rep = 1, n_rep
    1983             :          CALL section_vals_val_get(section, "_SECTION_PARAMETERS_", &
    1984         452 :                                    c_val=atm_name, i_rep_section=i_rep)
    1985         452 :          CALL uppercase(atm_name)
    1986         452 :          shell_list(start + i_rep)%atm_name = atm_name
    1987         452 :          CALL section_vals_val_get(section, "CORE_CHARGE", i_rep_section=i_rep, r_val=ccharge)
    1988         452 :          shell_list(start + i_rep)%shell%charge_core = ccharge
    1989         452 :          CALL section_vals_val_get(section, "SHELL_CHARGE", i_rep_section=i_rep, r_val=scharge)
    1990         452 :          shell_list(start + i_rep)%shell%charge_shell = scharge
    1991         452 :          CALL section_vals_val_get(section, "MASS_FRACTION", i_rep_section=i_rep, r_val=mfrac)
    1992         452 :          shell_list(start + i_rep)%shell%massfrac = mfrac
    1993         452 :          CALL section_vals_val_get(section, "K2_SPRING", i_rep_section=i_rep, r_val=k)
    1994         452 :          IF (k < 0.0_dp) THEN
    1995             :             CALL cp_abort(__LOCATION__, &
    1996             :                           "An invalid value was specified for the force constant k2 of the core-shell "// &
    1997           0 :                           "spring potential")
    1998             :          END IF
    1999         452 :          shell_list(start + i_rep)%shell%k2_spring = k
    2000         452 :          CALL section_vals_val_get(section, "K4_SPRING", i_rep_section=i_rep, r_val=k)
    2001         452 :          IF (k < 0.0_dp) THEN
    2002             :             CALL cp_abort(__LOCATION__, &
    2003             :                           "An invalid value was specified for the force constant k4 of the core-shell "// &
    2004           0 :                           "spring potential")
    2005             :          END IF
    2006         452 :          shell_list(start + i_rep)%shell%k4_spring = k
    2007         452 :          CALL section_vals_val_get(section, "MAX_DISTANCE", i_rep_section=i_rep, r_val=maxdist)
    2008         452 :          shell_list(start + i_rep)%shell%max_dist = maxdist
    2009         452 :          CALL section_vals_val_get(section, "SHELL_CUTOFF", i_rep_section=i_rep, r_val=cutoff)
    2010        1614 :          shell_list(start + i_rep)%shell%shell_cutoff = cutoff
    2011             :       END DO
    2012             : 
    2013         258 :    END SUBROUTINE read_shell_section
    2014             : 
    2015             : ! **************************************************************************************************
    2016             : !> \brief Reads the BONDS section
    2017             : !> \param bond_kind ...
    2018             : !> \param bond_a ...
    2019             : !> \param bond_b ...
    2020             : !> \param bond_k ...
    2021             : !> \param bond_r0 ...
    2022             : !> \param bond_cs ...
    2023             : !> \param section ...
    2024             : !> \param start ...
    2025             : !> \author teo
    2026             : ! **************************************************************************************************
    2027         975 :    SUBROUTINE read_bonds_section(bond_kind, bond_a, bond_b, bond_k, bond_r0, bond_cs, section, start)
    2028             :       INTEGER, DIMENSION(:), POINTER                     :: bond_kind
    2029             :       CHARACTER(LEN=default_string_length), &
    2030             :          DIMENSION(:), POINTER                           :: bond_a, bond_b
    2031             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: bond_k
    2032             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: bond_r0, bond_cs
    2033             :       TYPE(section_vals_type), POINTER                   :: section
    2034             :       INTEGER, INTENT(IN)                                :: start
    2035             : 
    2036             :       CHARACTER(LEN=default_string_length), &
    2037         975 :          DIMENSION(:), POINTER                           :: atm_names
    2038             :       INTEGER                                            :: isec, k, n_items
    2039         975 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: Kvals
    2040             : 
    2041         975 :       NULLIFY (Kvals, atm_names)
    2042         975 :       CALL section_vals_get(section, n_repetition=n_items)
    2043        2826 :       DO isec = 1, n_items
    2044        1851 :          CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=bond_kind(start + isec))
    2045        1851 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    2046        1851 :          bond_a(start + isec) = atm_names(1)
    2047        1851 :          bond_b(start + isec) = atm_names(2)
    2048        1851 :          CALL uppercase(bond_a(start + isec))
    2049        1851 :          CALL uppercase(bond_b(start + isec))
    2050        1851 :          CALL section_vals_val_get(section, "K", i_rep_section=isec, r_vals=Kvals)
    2051        1851 :          CPASSERT(SIZE(Kvals) <= 3)
    2052        7404 :          bond_k(:, start + isec) = 0.0_dp
    2053        3742 :          DO k = 1, SIZE(Kvals)
    2054        3742 :             bond_k(k, start + isec) = Kvals(k)
    2055             :          END DO
    2056        1851 :          CALL section_vals_val_get(section, "R0", i_rep_section=isec, r_val=bond_r0(start + isec))
    2057        2826 :          CALL section_vals_val_get(section, "CS", i_rep_section=isec, r_val=bond_cs(start + isec))
    2058             :       END DO
    2059         975 :    END SUBROUTINE read_bonds_section
    2060             : 
    2061             : ! **************************************************************************************************
    2062             : !> \brief Reads the BENDS section
    2063             : !> \param bend_kind ...
    2064             : !> \param bend_a ...
    2065             : !> \param bend_b ...
    2066             : !> \param bend_c ...
    2067             : !> \param bend_k ...
    2068             : !> \param bend_theta0 ...
    2069             : !> \param bend_cb ...
    2070             : !> \param bend_r012 ...
    2071             : !> \param bend_r032 ...
    2072             : !> \param bend_kbs12 ...
    2073             : !> \param bend_kbs32 ...
    2074             : !> \param bend_kss ...
    2075             : !> \param bend_legendre ...
    2076             : !> \param section ...
    2077             : !> \param start ...
    2078             : !> \author teo
    2079             : ! **************************************************************************************************
    2080         939 :    SUBROUTINE read_bends_section(bend_kind, bend_a, bend_b, bend_c, bend_k, bend_theta0, bend_cb, &
    2081             :                                  bend_r012, bend_r032, bend_kbs12, bend_kbs32, bend_kss, bend_legendre, &
    2082             :                                  section, start)
    2083             :       INTEGER, DIMENSION(:), POINTER                     :: bend_kind
    2084             :       CHARACTER(LEN=default_string_length), &
    2085             :          DIMENSION(:), POINTER                           :: bend_a, bend_b, bend_c
    2086             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: bend_k, bend_theta0, bend_cb, bend_r012, &
    2087             :                                                             bend_r032, bend_kbs12, bend_kbs32, &
    2088             :                                                             bend_kss
    2089             :       TYPE(legendre_data_type), DIMENSION(:), POINTER    :: bend_legendre
    2090             :       TYPE(section_vals_type), POINTER                   :: section
    2091             :       INTEGER, INTENT(IN)                                :: start
    2092             : 
    2093             :       CHARACTER(LEN=default_string_length), &
    2094         939 :          DIMENSION(:), POINTER                           :: atm_names
    2095             :       INTEGER                                            :: isec, k, n_items, n_rep
    2096         939 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: Kvals, r_values
    2097             : 
    2098         939 :       NULLIFY (Kvals, atm_names)
    2099         939 :       CALL section_vals_get(section, n_repetition=n_items)
    2100        3058 :       bend_legendre%order = 0
    2101        3058 :       DO isec = 1, n_items
    2102        2119 :          CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=bend_kind(start + isec))
    2103        2119 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    2104        2119 :          bend_a(start + isec) = atm_names(1)
    2105        2119 :          bend_b(start + isec) = atm_names(2)
    2106        2119 :          bend_c(start + isec) = atm_names(3)
    2107        2119 :          CALL uppercase(bend_a(start + isec))
    2108        2119 :          CALL uppercase(bend_b(start + isec))
    2109        2119 :          CALL uppercase(bend_c(start + isec))
    2110        2119 :          CALL section_vals_val_get(section, "K", i_rep_section=isec, r_vals=Kvals)
    2111        2119 :          CPASSERT(SIZE(Kvals) == 1)
    2112        2119 :          bend_k(start + isec) = Kvals(1)
    2113        2119 :          CALL section_vals_val_get(section, "THETA0", i_rep_section=isec, r_val=bend_theta0(start + isec))
    2114        2119 :          CALL section_vals_val_get(section, "CB", i_rep_section=isec, r_val=bend_cb(start + isec))
    2115        2119 :          CALL section_vals_val_get(section, "R012", i_rep_section=isec, r_val=bend_r012(start + isec))
    2116        2119 :          CALL section_vals_val_get(section, "R032", i_rep_section=isec, r_val=bend_r032(start + isec))
    2117        2119 :          CALL section_vals_val_get(section, "KBS12", i_rep_section=isec, r_val=bend_kbs12(start + isec))
    2118        2119 :          CALL section_vals_val_get(section, "KBS32", i_rep_section=isec, r_val=bend_kbs32(start + isec))
    2119        2119 :          CALL section_vals_val_get(section, "KSS", i_rep_section=isec, r_val=bend_kss(start + isec))
    2120             :          ! get legendre based data
    2121        2119 :          CALL section_vals_val_get(section, "LEGENDRE", i_rep_section=isec, n_rep_val=n_rep)
    2122        5177 :          DO k = 1, n_rep
    2123        2119 :             CALL section_vals_val_get(section, "LEGENDRE", i_rep_val=k, r_vals=r_values, i_rep_section=isec)
    2124        2119 :             bend_legendre(start + isec)%order = SIZE(r_values)
    2125        2119 :             IF (ASSOCIATED(bend_legendre(start + isec)%coeffs)) THEN
    2126           0 :                DEALLOCATE (bend_legendre(start + isec)%coeffs)
    2127             :             END IF
    2128        6357 :             ALLOCATE (bend_legendre(start + isec)%coeffs(bend_legendre(start + isec)%order))
    2129       10675 :             bend_legendre(start + isec)%coeffs = r_values
    2130             :          END DO
    2131             :       END DO
    2132         939 :    END SUBROUTINE read_bends_section
    2133             : 
    2134             : ! **************************************************************************************************
    2135             : !> \brief ...
    2136             : !> \param ub_kind ...
    2137             : !> \param ub_a ...
    2138             : !> \param ub_b ...
    2139             : !> \param ub_c ...
    2140             : !> \param ub_k ...
    2141             : !> \param ub_r0 ...
    2142             : !> \param section ...
    2143             : !> \param start ...
    2144             : ! **************************************************************************************************
    2145         939 :    SUBROUTINE read_ubs_section(ub_kind, ub_a, ub_b, ub_c, ub_k, ub_r0, section, start)
    2146             :       INTEGER, DIMENSION(:), POINTER                     :: ub_kind
    2147             :       CHARACTER(LEN=default_string_length), &
    2148             :          DIMENSION(:), POINTER                           :: ub_a, ub_b, ub_c
    2149             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: ub_k
    2150             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: ub_r0
    2151             :       TYPE(section_vals_type), POINTER                   :: section
    2152             :       INTEGER, INTENT(IN)                                :: start
    2153             : 
    2154             :       CHARACTER(LEN=default_string_length), &
    2155         939 :          DIMENSION(:), POINTER                           :: atm_names
    2156             :       INTEGER                                            :: isec, k, n_items
    2157             :       LOGICAL                                            :: explicit
    2158         939 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: Kvals
    2159             :       TYPE(section_vals_type), POINTER                   :: subsection
    2160             : 
    2161         939 :       NULLIFY (atm_names)
    2162         939 :       CALL section_vals_get(section, n_repetition=n_items)
    2163        3058 :       DO isec = 1, n_items
    2164        2119 :          subsection => section_vals_get_subs_vals(section, "UB", i_rep_section=isec)
    2165        2119 :          CALL section_vals_get(subsection, explicit=explicit)
    2166        3058 :          IF (explicit) THEN
    2167           4 :             CALL section_vals_val_get(subsection, "KIND", i_val=ub_kind(start + isec))
    2168           4 :             CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    2169           4 :             ub_a(start + isec) = atm_names(1)
    2170           4 :             ub_b(start + isec) = atm_names(2)
    2171           4 :             ub_c(start + isec) = atm_names(3)
    2172           4 :             CALL uppercase(ub_a(start + isec))
    2173           4 :             CALL uppercase(ub_b(start + isec))
    2174           4 :             CALL uppercase(ub_c(start + isec))
    2175           4 :             CALL section_vals_val_get(subsection, "K", r_vals=Kvals)
    2176           4 :             CPASSERT(SIZE(Kvals) <= 3)
    2177          16 :             ub_k(:, start + isec) = 0.0_dp
    2178          12 :             DO k = 1, SIZE(Kvals)
    2179          12 :                ub_k(k, start + isec) = Kvals(k)
    2180             :             END DO
    2181           4 :             CALL section_vals_val_get(subsection, "R0", r_val=ub_r0(start + isec))
    2182             :          END IF
    2183             :       END DO
    2184         939 :    END SUBROUTINE read_ubs_section
    2185             : 
    2186             : ! **************************************************************************************************
    2187             : !> \brief Reads the TORSIONS section
    2188             : !> \param torsion_kind ...
    2189             : !> \param torsion_a ...
    2190             : !> \param torsion_b ...
    2191             : !> \param torsion_c ...
    2192             : !> \param torsion_d ...
    2193             : !> \param torsion_k ...
    2194             : !> \param torsion_phi0 ...
    2195             : !> \param torsion_m ...
    2196             : !> \param section ...
    2197             : !> \param start ...
    2198             : !> \author teo
    2199             : ! **************************************************************************************************
    2200           6 :    SUBROUTINE read_torsions_section(torsion_kind, torsion_a, torsion_b, torsion_c, torsion_d, torsion_k, &
    2201             :                                     torsion_phi0, torsion_m, section, start)
    2202             :       INTEGER, DIMENSION(:), POINTER                     :: torsion_kind
    2203             :       CHARACTER(LEN=default_string_length), &
    2204             :          DIMENSION(:), POINTER                           :: torsion_a, torsion_b, torsion_c, &
    2205             :                                                             torsion_d
    2206             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: torsion_k, torsion_phi0
    2207             :       INTEGER, DIMENSION(:), POINTER                     :: torsion_m
    2208             :       TYPE(section_vals_type), POINTER                   :: section
    2209             :       INTEGER, INTENT(IN)                                :: start
    2210             : 
    2211             :       CHARACTER(LEN=default_string_length), &
    2212           6 :          DIMENSION(:), POINTER                           :: atm_names
    2213             :       INTEGER                                            :: isec, n_items
    2214             : 
    2215           6 :       NULLIFY (atm_names)
    2216           6 :       CALL section_vals_get(section, n_repetition=n_items)
    2217          44 :       DO isec = 1, n_items
    2218          38 :          CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=torsion_kind(start + isec))
    2219          38 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    2220          38 :          torsion_a(start + isec) = atm_names(1)
    2221          38 :          torsion_b(start + isec) = atm_names(2)
    2222          38 :          torsion_c(start + isec) = atm_names(3)
    2223          38 :          torsion_d(start + isec) = atm_names(4)
    2224          38 :          CALL uppercase(torsion_a(start + isec))
    2225          38 :          CALL uppercase(torsion_b(start + isec))
    2226          38 :          CALL uppercase(torsion_c(start + isec))
    2227          38 :          CALL uppercase(torsion_d(start + isec))
    2228          38 :          CALL section_vals_val_get(section, "K", i_rep_section=isec, r_val=torsion_k(start + isec))
    2229          38 :          CALL section_vals_val_get(section, "PHI0", i_rep_section=isec, r_val=torsion_phi0(start + isec))
    2230          38 :          CALL section_vals_val_get(section, "M", i_rep_section=isec, i_val=torsion_m(start + isec))
    2231             :          ! Modify parameterisation for OPLS case
    2232          44 :          IF (torsion_kind(start + isec) .EQ. do_ff_opls) THEN
    2233          12 :             IF (torsion_phi0(start + isec) .NE. 0.0_dp) THEN
    2234             :                CALL cp_warn(__LOCATION__, "PHI0 parameter was non-zero "// &
    2235           0 :                             "for an OPLS-type TORSION.  It will be ignored.")
    2236             :             END IF
    2237          12 :             IF (MODULO(torsion_m(start + isec), 2) .EQ. 0) THEN
    2238             :                ! For even M, negate the cosine using a Pi phase factor
    2239           2 :                torsion_phi0(start + isec) = pi
    2240             :             END IF
    2241             :             ! the K parameter appears as K/2 in the OPLS parameterisation
    2242          12 :             torsion_k(start + isec) = torsion_k(start + isec)*0.5_dp
    2243             :          END IF
    2244             :       END DO
    2245           6 :    END SUBROUTINE read_torsions_section
    2246             : 
    2247             : ! **************************************************************************************************
    2248             : !> \brief Reads the IMPROPER section
    2249             : !> \param impr_kind ...
    2250             : !> \param impr_a ...
    2251             : !> \param impr_b ...
    2252             : !> \param impr_c ...
    2253             : !> \param impr_d ...
    2254             : !> \param impr_k ...
    2255             : !> \param impr_phi0 ...
    2256             : !> \param section ...
    2257             : !> \param start ...
    2258             : !> \author louis vanduyfhuys
    2259             : ! **************************************************************************************************
    2260           8 :    SUBROUTINE read_improper_section(impr_kind, impr_a, impr_b, impr_c, impr_d, impr_k, &
    2261             :                                     impr_phi0, section, start)
    2262             :       INTEGER, DIMENSION(:), POINTER                     :: impr_kind
    2263             :       CHARACTER(LEN=default_string_length), &
    2264             :          DIMENSION(:), POINTER                           :: impr_a, impr_b, impr_c, impr_d
    2265             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: impr_k, impr_phi0
    2266             :       TYPE(section_vals_type), POINTER                   :: section
    2267             :       INTEGER, INTENT(IN)                                :: start
    2268             : 
    2269             :       CHARACTER(LEN=default_string_length), &
    2270           8 :          DIMENSION(:), POINTER                           :: atm_names
    2271             :       INTEGER                                            :: isec, n_items
    2272             : 
    2273           8 :       NULLIFY (atm_names)
    2274           8 :       CALL section_vals_get(section, n_repetition=n_items)
    2275          16 :       DO isec = 1, n_items
    2276           8 :          CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=impr_kind(start + isec))
    2277           8 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    2278           8 :          impr_a(start + isec) = atm_names(1)
    2279           8 :          impr_b(start + isec) = atm_names(2)
    2280           8 :          impr_c(start + isec) = atm_names(3)
    2281           8 :          impr_d(start + isec) = atm_names(4)
    2282           8 :          CALL uppercase(impr_a(start + isec))
    2283           8 :          CALL uppercase(impr_b(start + isec))
    2284           8 :          CALL uppercase(impr_c(start + isec))
    2285           8 :          CALL uppercase(impr_d(start + isec))
    2286           8 :          CALL section_vals_val_get(section, "K", i_rep_section=isec, r_val=impr_k(start + isec))
    2287          16 :          CALL section_vals_val_get(section, "PHI0", i_rep_section=isec, r_val=impr_phi0(start + isec))
    2288             :       END DO
    2289           8 :    END SUBROUTINE read_improper_section
    2290             : 
    2291             : ! **************************************************************************************************
    2292             : !> \brief Reads the OPBEND section
    2293             : !> \param opbend_kind ...
    2294             : !> \param opbend_a ...
    2295             : !> \param opbend_b ...
    2296             : !> \param opbend_c ...
    2297             : !> \param opbend_d ...
    2298             : !> \param opbend_k ...
    2299             : !> \param opbend_phi0 ...
    2300             : !> \param section ...
    2301             : !> \param start ...
    2302             : !> \author louis vanduyfhuys
    2303             : ! **************************************************************************************************
    2304           2 :    SUBROUTINE read_opbend_section(opbend_kind, opbend_a, opbend_b, opbend_c, opbend_d, opbend_k, &
    2305             :                                   opbend_phi0, section, start)
    2306             :       INTEGER, DIMENSION(:), POINTER                     :: opbend_kind
    2307             :       CHARACTER(LEN=default_string_length), &
    2308             :          DIMENSION(:), POINTER                           :: opbend_a, opbend_b, opbend_c, opbend_d
    2309             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: opbend_k, opbend_phi0
    2310             :       TYPE(section_vals_type), POINTER                   :: section
    2311             :       INTEGER, INTENT(IN)                                :: start
    2312             : 
    2313             :       CHARACTER(LEN=default_string_length), &
    2314           2 :          DIMENSION(:), POINTER                           :: atm_names
    2315             :       INTEGER                                            :: isec, n_items
    2316             : 
    2317           2 :       NULLIFY (atm_names)
    2318           2 :       CALL section_vals_get(section, n_repetition=n_items)
    2319           4 :       DO isec = 1, n_items
    2320           2 :          CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=opbend_kind(start + isec))
    2321           2 :          CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
    2322           2 :          opbend_a(start + isec) = atm_names(1)
    2323           2 :          opbend_b(start + isec) = atm_names(2)
    2324           2 :          opbend_c(start + isec) = atm_names(3)
    2325           2 :          opbend_d(start + isec) = atm_names(4)
    2326           2 :          CALL uppercase(opbend_a(start + isec))
    2327           2 :          CALL uppercase(opbend_b(start + isec))
    2328           2 :          CALL uppercase(opbend_c(start + isec))
    2329           2 :          CALL uppercase(opbend_d(start + isec))
    2330           2 :          CALL section_vals_val_get(section, "K", i_rep_section=isec, r_val=opbend_k(start + isec))
    2331           4 :          CALL section_vals_val_get(section, "PHI0", i_rep_section=isec, r_val=opbend_phi0(start + isec))
    2332             :       END DO
    2333           2 :    END SUBROUTINE read_opbend_section
    2334             : 
    2335             : ! **************************************************************************************************
    2336             : !> \brief Reads the force_field input section
    2337             : !> \param ff_type ...
    2338             : !> \param para_env ...
    2339             : !> \param mm_section ...
    2340             : !> \par History
    2341             : !>      JGH (30.11.2001) : moved determination of setup variables to
    2342             : !>                         molecule_input
    2343             : !> \author CJM
    2344             : ! **************************************************************************************************
    2345        2639 :    SUBROUTINE read_force_field_section(ff_type, para_env, mm_section)
    2346             :       TYPE(force_field_type), INTENT(INOUT)              :: ff_type
    2347             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2348             :       TYPE(section_vals_type), POINTER                   :: mm_section
    2349             : 
    2350             :       TYPE(section_vals_type), POINTER                   :: ff_section
    2351             : 
    2352             :       NULLIFY (ff_section)
    2353        2639 :       ff_section => section_vals_get_subs_vals(mm_section, "FORCEFIELD")
    2354        2639 :       CALL read_force_field_section1(ff_section, mm_section, ff_type, para_env)
    2355        2639 :    END SUBROUTINE read_force_field_section
    2356             : 
    2357             : ! **************************************************************************************************
    2358             : !> \brief reads EAM potential from library
    2359             : !> \param eam ...
    2360             : !> \param para_env ...
    2361             : !> \param mm_section ...
    2362             : ! **************************************************************************************************
    2363          40 :    SUBROUTINE read_eam_data(eam, para_env, mm_section)
    2364             :       TYPE(eam_pot_type), POINTER                        :: eam
    2365             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2366             :       TYPE(section_vals_type), POINTER                   :: mm_section
    2367             : 
    2368             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'read_eam_data'
    2369             : 
    2370             :       INTEGER                                            :: handle, i, iw
    2371             :       TYPE(cp_logger_type), POINTER                      :: logger
    2372             :       TYPE(cp_parser_type)                               :: parser
    2373             : 
    2374          20 :       CALL timeset(routineN, handle)
    2375          20 :       NULLIFY (logger)
    2376          20 :       logger => cp_get_default_logger()
    2377             :       iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
    2378          20 :                                 extension=".mmLog")
    2379          20 :       IF (iw > 0) WRITE (iw, *) "Reading EAM data from: ", TRIM(eam%eam_file_name)
    2380          20 :       CALL parser_create(parser, TRIM(eam%eam_file_name), para_env=para_env)
    2381             : 
    2382          20 :       CALL parser_get_next_line(parser, 1)
    2383          20 :       IF (iw > 0) WRITE (iw, *) "Title: ", parser%input_line
    2384             : 
    2385          20 :       CALL parser_get_next_line(parser, 2)
    2386          20 :       READ (parser%input_line, *) eam%drar, eam%drhoar, eam%acutal, eam%npoints
    2387          20 :       eam%drar = cp_unit_to_cp2k(eam%drar, "angstrom")
    2388          20 :       eam%acutal = cp_unit_to_cp2k(eam%acutal, "angstrom")
    2389             :       ! Relocating arrays with the right size
    2390          20 :       CALL reallocate(eam%rho, 1, eam%npoints)
    2391          20 :       CALL reallocate(eam%rhop, 1, eam%npoints)
    2392          20 :       CALL reallocate(eam%rval, 1, eam%npoints)
    2393          20 :       CALL reallocate(eam%rhoval, 1, eam%npoints)
    2394          20 :       CALL reallocate(eam%phi, 1, eam%npoints)
    2395          20 :       CALL reallocate(eam%phip, 1, eam%npoints)
    2396          20 :       CALL reallocate(eam%frho, 1, eam%npoints)
    2397          20 :       CALL reallocate(eam%frhop, 1, eam%npoints)
    2398             :       ! Reading density and derivative of density (with respect to r)
    2399       64020 :       DO i = 1, eam%npoints
    2400       64000 :          CALL parser_get_next_line(parser, 1)
    2401       64000 :          READ (parser%input_line, *) eam%rho(i), eam%rhop(i)
    2402       64000 :          eam%rhop(i) = cp_unit_to_cp2k(eam%rhop(i), "angstrom^-1")
    2403       64000 :          eam%rval(i) = REAL(i - 1, KIND=dp)*eam%drar
    2404       64020 :          eam%rhoval(i) = REAL(i - 1, KIND=dp)*eam%drhoar
    2405             :       END DO
    2406             :       ! Reading pair potential PHI and its derivative (with respect to r)
    2407       64020 :       DO i = 1, eam%npoints
    2408       64000 :          CALL parser_get_next_line(parser, 1)
    2409       64000 :          READ (parser%input_line, *) eam%phi(i), eam%phip(i)
    2410       64000 :          eam%phi(i) = cp_unit_to_cp2k(eam%phi(i), "eV")
    2411       64020 :          eam%phip(i) = cp_unit_to_cp2k(eam%phip(i), "eV*angstrom^-1")
    2412             :       END DO
    2413             :       ! Reading embedded function and its derivative (with respect to density)
    2414       64020 :       DO i = 1, eam%npoints
    2415       64000 :          CALL parser_get_next_line(parser, 1)
    2416       64000 :          READ (parser%input_line, *) eam%frho(i), eam%frhop(i)
    2417       64000 :          eam%frho(i) = cp_unit_to_cp2k(eam%frho(i), "eV")
    2418       64020 :          eam%frhop(i) = cp_unit_to_cp2k(eam%frhop(i), "eV")
    2419             :       END DO
    2420             : 
    2421          20 :       IF (iw > 0) WRITE (iw, *) "Finished EAM data"
    2422          20 :       CALL parser_release(parser)
    2423          20 :       CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO")
    2424          20 :       CALL timestop(handle)
    2425             : 
    2426          60 :    END SUBROUTINE read_eam_data
    2427             : 
    2428             : ! **************************************************************************************************
    2429             : !> \brief reads NequIP potential from .pth file
    2430             : !> \param nequip ...
    2431             : !> \author Gabriele Tocci
    2432             : ! **************************************************************************************************
    2433          12 :    SUBROUTINE read_nequip_data(nequip)
    2434             :       TYPE(nequip_pot_type), POINTER                     :: nequip
    2435             : 
    2436             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'read_nequip_data'
    2437             :       CHARACTER(LEN=1), PARAMETER                        :: delimiter = ' '
    2438             : 
    2439          12 :       CHARACTER(LEN=100), ALLOCATABLE, DIMENSION(:)      :: tokenized_string
    2440             :       CHARACTER(LEN=default_path_length)                 :: allow_tf32_str, cutoff_str, types_str
    2441             :       INTEGER                                            :: handle
    2442             :       LOGICAL                                            :: allow_tf32, found
    2443             : 
    2444          12 :       CALL timeset(routineN, handle)
    2445             : 
    2446          12 :       INQUIRE (FILE=nequip%nequip_file_name, EXIST=found)
    2447          12 :       IF (.NOT. found) THEN
    2448             :          CALL cp_abort(__LOCATION__, &
    2449             :                        "Nequip model file <"//TRIM(nequip%nequip_file_name)// &
    2450           0 :                        "> not found.")
    2451             :       END IF
    2452             : 
    2453          12 :       nequip%nequip_version = torch_model_read_metadata(nequip%nequip_file_name, "nequip_version")
    2454          12 :       cutoff_str = torch_model_read_metadata(nequip%nequip_file_name, "r_max")
    2455          12 :       types_str = torch_model_read_metadata(nequip%nequip_file_name, "type_names")
    2456          12 :       CALL tokenize_string(TRIM(types_str), delimiter, tokenized_string)
    2457             : 
    2458          12 :       IF (ALLOCATED(nequip%type_names_torch)) THEN
    2459           0 :          DEALLOCATE (nequip%type_names_torch)
    2460             :       END IF
    2461          36 :       ALLOCATE (nequip%type_names_torch(SIZE(tokenized_string)))
    2462             : 
    2463          36 :       nequip%type_names_torch(:) = tokenized_string(:)
    2464             : 
    2465          12 :       READ (cutoff_str, *) nequip%rcutsq
    2466          12 :       nequip%rcutsq = cp_unit_to_cp2k(nequip%rcutsq, nequip%unit_coords)
    2467          12 :       nequip%rcutsq = nequip%rcutsq*nequip%rcutsq
    2468          12 :       nequip%unit_coords_val = cp_unit_to_cp2k(nequip%unit_coords_val, nequip%unit_coords)
    2469          12 :       nequip%unit_forces_val = cp_unit_to_cp2k(nequip%unit_forces_val, nequip%unit_forces)
    2470          12 :       nequip%unit_energy_val = cp_unit_to_cp2k(nequip%unit_energy_val, nequip%unit_energy)
    2471          12 :       nequip%unit_cell_val = cp_unit_to_cp2k(nequip%unit_cell_val, nequip%unit_cell)
    2472             : 
    2473          18 :       IF (torch_model_read_metadata(nequip%nequip_file_name, "default_dtype") == "float32" .AND. &
    2474             :           torch_model_read_metadata(nequip%nequip_file_name, "model_dtype") == "float32") THEN
    2475           6 :          nequip%do_nequip_sp = .TRUE.
    2476           6 :       ELSE IF (torch_model_read_metadata(nequip%nequip_file_name, "default_dtype") == "float64" .AND. &
    2477          12 :                torch_model_read_metadata(nequip%nequip_file_name, "model_dtype") == "float64") THEN
    2478           6 :          nequip%do_nequip_sp = .FALSE.
    2479             :       ELSE
    2480             :          CALL cp_abort(__LOCATION__, &
    2481             :                        "Both default_dtype and model_dtype should be either float32 or float64. Currently, default_dtype is <"// &
    2482             :                        torch_model_read_metadata(nequip%nequip_file_name, "default_dtype")//"> and model_dtype is <"// &
    2483           6 :                        torch_model_read_metadata(nequip%nequip_file_name, "model_dtype")//">.")
    2484             :       END IF
    2485             : 
    2486          12 :       allow_tf32_str = torch_model_read_metadata(nequip%nequip_file_name, "allow_tf32")
    2487          12 :       allow_tf32 = (TRIM(allow_tf32_str) == "1")
    2488          12 :       IF (TRIM(allow_tf32_str) /= "1" .AND. TRIM(allow_tf32_str) /= "0") THEN
    2489             :          CALL cp_abort(__LOCATION__, &
    2490             :                        "The value for allow_tf32 <"//TRIM(allow_tf32_str)// &
    2491           0 :                        "> is not supported. Check the .yaml and .pth files.")
    2492             :       END IF
    2493          12 :       CALL torch_allow_tf32(allow_tf32)
    2494             : 
    2495          12 :       CALL timestop(handle)
    2496          24 :    END SUBROUTINE read_nequip_data
    2497             : 
    2498             : ! **************************************************************************************************
    2499             : !> \brief reads ALLEGRO potential from .pth file
    2500             : !> \param allegro ...
    2501             : !> \author Gabriele Tocci
    2502             : ! **************************************************************************************************
    2503           8 :    SUBROUTINE read_allegro_data(allegro)
    2504             :       TYPE(allegro_pot_type), POINTER                    :: allegro
    2505             : 
    2506             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'read_allegro_data'
    2507             :       CHARACTER(LEN=1), PARAMETER                        :: delimiter = ' '
    2508             : 
    2509           8 :       CHARACTER(LEN=100), ALLOCATABLE, DIMENSION(:)      :: tokenized_string
    2510             :       CHARACTER(LEN=default_path_length)                 :: allow_tf32_str, cutoff_str, types_str
    2511             :       INTEGER                                            :: handle
    2512             :       LOGICAL                                            :: allow_tf32, found
    2513             : 
    2514           8 :       CALL timeset(routineN, handle)
    2515             : 
    2516           8 :       INQUIRE (FILE=allegro%allegro_file_name, EXIST=found)
    2517           8 :       IF (.NOT. found) THEN
    2518             :          CALL cp_abort(__LOCATION__, &
    2519             :                        "Allegro model file <"//TRIM(allegro%allegro_file_name)// &
    2520           0 :                        "> not found.")
    2521             :       END IF
    2522             : 
    2523           8 :       allegro%nequip_version = torch_model_read_metadata(allegro%allegro_file_name, "nequip_version")
    2524           8 :       IF (allegro%nequip_version == "") THEN
    2525             :          CALL cp_abort(__LOCATION__, &
    2526             :                        "Allegro model file <"//TRIM(allegro%allegro_file_name)// &
    2527           0 :                        "> has not been deployed; did you forget to run `nequip-deploy`?")
    2528             :       END IF
    2529           8 :       cutoff_str = torch_model_read_metadata(allegro%allegro_file_name, "r_max")
    2530           8 :       types_str = torch_model_read_metadata(allegro%allegro_file_name, "type_names")
    2531           8 :       CALL tokenize_string(TRIM(types_str), delimiter, tokenized_string)
    2532             : 
    2533           8 :       IF (ALLOCATED(allegro%type_names_torch)) THEN
    2534           0 :          DEALLOCATE (allegro%type_names_torch)
    2535             :       END IF
    2536          24 :       ALLOCATE (allegro%type_names_torch(SIZE(tokenized_string)))
    2537          28 :       allegro%type_names_torch(:) = tokenized_string(:)
    2538             : 
    2539           8 :       READ (cutoff_str, *) allegro%rcutsq
    2540           8 :       allegro%rcutsq = cp_unit_to_cp2k(allegro%rcutsq, allegro%unit_coords)
    2541           8 :       allegro%rcutsq = allegro%rcutsq*allegro%rcutsq
    2542           8 :       allegro%unit_coords_val = cp_unit_to_cp2k(allegro%unit_coords_val, allegro%unit_coords)
    2543           8 :       allegro%unit_forces_val = cp_unit_to_cp2k(allegro%unit_forces_val, allegro%unit_forces)
    2544           8 :       allegro%unit_energy_val = cp_unit_to_cp2k(allegro%unit_energy_val, allegro%unit_energy)
    2545           8 :       allegro%unit_cell_val = cp_unit_to_cp2k(allegro%unit_cell_val, allegro%unit_cell)
    2546             : 
    2547          10 :       IF (torch_model_read_metadata(allegro%allegro_file_name, "default_dtype") == "float32" .AND. &
    2548             :           torch_model_read_metadata(allegro%allegro_file_name, "model_dtype") == "float32") THEN
    2549           6 :          allegro%do_allegro_sp = .TRUE.
    2550           2 :       ELSE IF (torch_model_read_metadata(allegro%allegro_file_name, "default_dtype") == "float64" .AND. &
    2551           8 :                torch_model_read_metadata(allegro%allegro_file_name, "model_dtype") == "float64") THEN
    2552           2 :          allegro%do_allegro_sp = .FALSE.
    2553             :       ELSE
    2554             :          CALL cp_abort(__LOCATION__, &
    2555             :                        "Both default_dtype and model_dtype should be either float32 or float64. Currently, default_dtype is <"// &
    2556             :                        torch_model_read_metadata(allegro%allegro_file_name, "default_dtype")//"> and model_dtype is <"// &
    2557           2 :                        torch_model_read_metadata(allegro%allegro_file_name, "model_dtype")//">.")
    2558             :       END IF
    2559             : 
    2560           8 :       allow_tf32_str = torch_model_read_metadata(allegro%allegro_file_name, "allow_tf32")
    2561           8 :       allow_tf32 = (TRIM(allow_tf32_str) == "1")
    2562           8 :       IF (TRIM(allow_tf32_str) /= "1" .AND. TRIM(allow_tf32_str) /= "0") THEN
    2563             :          CALL cp_abort(__LOCATION__, &
    2564             :                        "The value for allow_tf32 <"//TRIM(allow_tf32_str)// &
    2565           0 :                        "> is not supported. Check the .yaml and .pth files.")
    2566             :       END IF
    2567           8 :       CALL torch_allow_tf32(allow_tf32)
    2568             : 
    2569           8 :       CALL timestop(handle)
    2570          16 :    END SUBROUTINE read_allegro_data
    2571             : 
    2572             : ! **************************************************************************************************
    2573             : !> \brief returns tokenized string of kinds from .pth file
    2574             : !> \param element ...
    2575             : !> \param delimiter ...
    2576             : !> \param tokenized_array ...
    2577             : !> \author Maria Bilichenko
    2578             : ! **************************************************************************************************
    2579          20 :    SUBROUTINE tokenize_string(element, delimiter, tokenized_array)
    2580             :       CHARACTER(LEN=*), INTENT(IN)                       :: element
    2581             :       CHARACTER(LEN=1), INTENT(IN)                       :: delimiter
    2582             :       CHARACTER(LEN=100), ALLOCATABLE, DIMENSION(:), &
    2583             :          INTENT(OUT)                                     :: tokenized_array
    2584             : 
    2585             :       CHARACTER(LEN=100)                                 :: temp_kinds
    2586             :       INTEGER                                            :: end_pos, i, num_elements, start
    2587          20 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: delim_positions
    2588             : 
    2589             :       ! Find positions of delimiter within element
    2590          60 :       ALLOCATE (delim_positions(LEN(element)))
    2591          90 :       delim_positions = .FALSE.
    2592             : 
    2593          90 :       DO i = 1, LEN(element)
    2594          90 :          IF (element(i:i) == delimiter) delim_positions(i) = .TRUE.
    2595             :       END DO
    2596             : 
    2597          90 :       num_elements = COUNT(delim_positions) + 1
    2598             : 
    2599          60 :       ALLOCATE (tokenized_array(num_elements))
    2600             : 
    2601          20 :       start = 1
    2602          64 :       DO i = 1, num_elements
    2603         210 :          IF (LEN(element) < 3 .AND. COUNT(delim_positions) == 0) THEN ! if there is 1 kind only and it has one or two
    2604             :             !characters (C or Cl) the end_pos will be the index of the last character (1 or 2)
    2605             :             end_pos = LEN(element)
    2606             :          ELSE ! else, the end_pos is determined by the index of the space - 1
    2607          42 :             end_pos = find_end_pos(start, delim_positions)
    2608             :          END IF
    2609          44 :          temp_kinds = element(start:end_pos)
    2610          44 :          IF (TRIM(temp_kinds) /= '') THEN
    2611          44 :             tokenized_array(i) = temp_kinds
    2612             :          END IF
    2613          64 :          start = end_pos + 2
    2614             :       END DO
    2615          20 :       DEALLOCATE (delim_positions)
    2616          20 :    END SUBROUTINE tokenize_string
    2617             : 
    2618             : ! **************************************************************************************************
    2619             : !> \brief finds the position of the atom by the spacing
    2620             : !> \param start ...
    2621             : !> \param delim_positions ...
    2622             : !> \return ...
    2623             : !> \author Maria Bilichenko
    2624             : ! **************************************************************************************************
    2625          42 :    INTEGER FUNCTION find_end_pos(start, delim_positions)
    2626             :       INTEGER, INTENT(IN)                                :: start
    2627             :       LOGICAL, DIMENSION(:), INTENT(IN)                  :: delim_positions
    2628             : 
    2629             :       INTEGER                                            :: end_pos, i
    2630             : 
    2631          42 :       end_pos = start
    2632          84 :       DO i = start, SIZE(delim_positions)
    2633          84 :          IF (delim_positions(i)) THEN
    2634          24 :             end_pos = i - 1
    2635          24 :             EXIT
    2636             :          END IF
    2637             :       END DO
    2638             : 
    2639          42 :       find_end_pos = end_pos
    2640          42 :    END FUNCTION find_end_pos
    2641             : 
    2642             : ! **************************************************************************************************
    2643             : !> \brief checks if all the ATOMS from *.inp file are available in *.pth file
    2644             : !> \param cp2k_inp_atom_types ...
    2645             : !> \param torch_atom_types ...
    2646             : !> \author Maria Bilichenko
    2647             : ! **************************************************************************************************
    2648          20 :    SUBROUTINE check_cp2k_atom_names_in_torch(cp2k_inp_atom_types, torch_atom_types)
    2649             :       CHARACTER(LEN=*), DIMENSION(:), INTENT(IN)         :: cp2k_inp_atom_types, torch_atom_types
    2650             : 
    2651             :       INTEGER                                            :: i, j
    2652             :       LOGICAL                                            :: found
    2653             : 
    2654          58 :       DO i = 1, SIZE(cp2k_inp_atom_types)
    2655          38 :          found = .FALSE.
    2656          62 :          DO j = 1, SIZE(torch_atom_types)
    2657          62 :             IF (TRIM(cp2k_inp_atom_types(i)) == TRIM(torch_atom_types(j))) THEN
    2658             :                found = .TRUE.
    2659             :                EXIT
    2660             :             END IF
    2661             :          END DO
    2662          58 :          IF (.NOT. found) THEN
    2663             :             CALL cp_abort(__LOCATION__, &
    2664             :                           "Atom "//TRIM(cp2k_inp_atom_types(i))// &
    2665           0 :                           " is defined in the CP2K input file but is missing in the torch model file")
    2666             :          END IF
    2667             :       END DO
    2668          20 :    END SUBROUTINE check_cp2k_atom_names_in_torch
    2669             : 
    2670             : ! **************************************************************************************************
    2671             : !> \brief reads TABPOT potential from file
    2672             : !> \param tab ...
    2673             : !> \param para_env ...
    2674             : !> \param mm_section ...
    2675             : !> \author Da Teng, Alex Mironenko
    2676             : ! **************************************************************************************************
    2677          48 :    SUBROUTINE read_tabpot_data(tab, para_env, mm_section)
    2678             :       TYPE(tab_pot_type), POINTER                        :: tab
    2679             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2680             :       TYPE(section_vals_type), POINTER                   :: mm_section
    2681             : 
    2682             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'read_tabpot_data'
    2683             : 
    2684             :       CHARACTER                                          :: d1, d2
    2685             :       INTEGER                                            :: d, handle, i, iw
    2686             :       TYPE(cp_logger_type), POINTER                      :: logger
    2687             :       TYPE(cp_parser_type)                               :: parser
    2688             : 
    2689          24 :       CALL timeset(routineN, handle)
    2690          24 :       NULLIFY (logger)
    2691          24 :       logger => cp_get_default_logger()
    2692             :       iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
    2693          24 :                                 extension=".mmLog")
    2694          24 :       IF (iw > 0) WRITE (iw, *) "Reading TABPOT data from: ", TRIM(tab%tabpot_file_name)
    2695          24 :       CALL parser_create(parser, TRIM(tab%tabpot_file_name), para_env=para_env)
    2696          24 :       CALL parser_get_next_line(parser, 1)
    2697          24 :       IF (iw > 0) WRITE (iw, *) "Title: ", parser%input_line
    2698          24 :       CALL parser_get_next_line(parser, 1)
    2699             : 
    2700             :       ! example format: N 1000 R 1.00 20.0
    2701             :       ! Assume the data is evenly spaced
    2702          24 :       READ (parser%input_line, *) d1, tab%npoints, d2, tab%dr, tab%rcut
    2703             : 
    2704             :       ! Relocating arrays with the right size
    2705          24 :       CALL reallocate(tab%r, 1, tab%npoints)
    2706          24 :       CALL reallocate(tab%e, 1, tab%npoints)
    2707          24 :       CALL reallocate(tab%f, 1, tab%npoints)
    2708             : 
    2709             :       ! Reading r, e, f
    2710       21912 :       DO i = 1, tab%npoints
    2711       21888 :          CALL parser_get_next_line(parser, 1)
    2712       21888 :          READ (parser%input_line, *) d, tab%r(i), tab%e(i), tab%f(i)
    2713       21888 :          tab%r(i) = cp_unit_to_cp2k(tab%r(i), "angstrom")
    2714       21888 :          tab%e(i) = cp_unit_to_cp2k(tab%e(i), "kcalmol")
    2715       21912 :          tab%f(i) = cp_unit_to_cp2k(tab%f(i), "kcalmol*angstrom^-1")
    2716             :       END DO
    2717             : 
    2718          24 :       tab%dr = tab%r(2) - tab%r(1)
    2719          24 :       tab%rcut = cp_unit_to_cp2k(tab%rcut, "angstrom")
    2720             : 
    2721          24 :       IF (iw > 0) WRITE (iw, *) "Finished TABPOT data"
    2722          24 :       CALL parser_release(parser)
    2723          24 :       CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO")
    2724          24 :       CALL timestop(handle)
    2725          72 :    END SUBROUTINE read_tabpot_data
    2726             : END MODULE force_fields_input

Generated by: LCOV version 1.15