Forcefield: Add table potential type TABPOT

cp2k · Oct 7, 2022 · 6df1cbb · 6df1cbb
1 parent 829ad49
commit 6df1cbb
Show file tree

Hide file tree

Showing 16 changed files with 8,417 additions and 32 deletions.
diff --git a/src/force_fields_input.F b/src/force_fields_input.F
@@ -67,7 +67,7 @@ MODULE force_fields_input
         ft_type, ftd_type, gal21_type, gal_type, gp_type, gw_type, ip_type, ipbv_pot_type, &
         lj_charmm_type, no_potential_single_allocation, pair_potential_p_type, &
         pair_potential_reallocate, potential_single_allocation, quip_type, siepmann_type, &
-        tersoff_type, wl_type
+        tab_pot_type, tab_type, tersoff_type, wl_type
    USE shell_potential_types,           ONLY: shell_p_create,&
                                               shell_p_type
    USE string_utilities,                ONLY: uppercase
@@ -101,7 +101,8 @@ SUBROUTINE read_force_field_section1(ff_section, mm_section, ff_type, para_env)
       TYPE(cp_para_env_type), POINTER                    :: para_env
 
       INTEGER :: nb4, nbends, nbm, nbmhft, nbmhftd, nbonds, nchg, neam, ngal, ngal21, ngd, ngp, &
-         nimpr, nipbv, nlj, nopbend, nquip, nshell, nsiepmann, ntersoff, ntors, ntot, nubs, nwl
+         nimpr, nipbv, nlj, nopbend, nquip, nshell, nsiepmann, ntab, ntersoff, ntors, ntot, nubs, &
+         nwl
       LOGICAL                                            :: explicit, unique_spline
       REAL(KIND=dp)                                      :: min_eps_spline_allowed
       TYPE(input_info_type), POINTER                     :: inp_info
@@ -277,6 +278,14 @@ SUBROUTINE read_force_field_section1(ff_section, mm_section, ff_type, para_env)
             CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + nquip, quip=.TRUE.)
             CALL read_quip_section(inp_info%nonbonded, tmp_section2, ntot)
          END IF
+
+         tmp_section2 => section_vals_get_subs_vals(tmp_section, "TABPOT")
+         CALL section_vals_get(tmp_section2, explicit=explicit, n_repetition=ntab)
+         ntot = nlj + nwl + neam + ngd + nipbv + nbmhft + nbmhftd + nb4 + nbm + ngp + ntersoff + ngal + ngal21 + nsiepmann + nquip
+         IF (explicit) THEN
+            CALL pair_potential_reallocate(inp_info%nonbonded, 1, ntot + ntab, tab=.TRUE.)
+            CALL read_tabpot_section(inp_info%nonbonded, tmp_section2, ntot, para_env, mm_section)
+         END IF
       END IF
 
       tmp_section => section_vals_get_subs_vals(ff_section, "NONBONDED14")
@@ -1603,6 +1612,44 @@ SUBROUTINE read_bm_section(nonbonded, section, start)
       END DO
    END SUBROUTINE read_bm_section
 
+! **************************************************************************************************
+!> \brief Reads the TABPOT section
+!> \param nonbonded ...
+!> \param section ...
+!> \param start ...
+!> \param para_env ...
+!> \param mm_section ...
+!> \author Alex Mironenko, Da Teng
+! **************************************************************************************************
+   SUBROUTINE read_tabpot_section(nonbonded, section, start, para_env, mm_section)
+      TYPE(pair_potential_p_type), POINTER               :: nonbonded
+      TYPE(section_vals_type), POINTER                   :: section
+      INTEGER, INTENT(IN)                                :: start
+      TYPE(cp_para_env_type), POINTER                    :: para_env
+      TYPE(section_vals_type), POINTER                   :: mm_section
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'read_tabpot_section', &
+         routineP = moduleN//':'//routineN
+
+      CHARACTER(LEN=default_string_length), &
+         DIMENSION(:), POINTER                           :: atm_names
+      INTEGER                                            :: isec, n_items
+
+      CALL section_vals_get(section, n_repetition=n_items)
+      DO isec = 1, n_items
+         CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
+         nonbonded%pot(start + isec)%pot%type = tab_type
+         nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
+         nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
+         CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
+         CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
+         CALL section_vals_val_get(section, "PARM_FILE_NAME", i_rep_section=isec, &
+                                   c_val=nonbonded%pot(start + isec)%pot%set(1)%tab%tabpot_file_name)
+         CALL read_tabpot_data(nonbonded%pot(start + isec)%pot%set(1)%tab, para_env, mm_section)
+         nonbonded%pot(start + isec)%pot%set(1)%tab%index = isec
+      END DO
+   END SUBROUTINE read_tabpot_section
+
 ! **************************************************************************************************
 !> \brief Reads the CHARGE section
 !> \param charge_atm ...
@@ -2197,4 +2244,62 @@ SUBROUTINE read_eam_data(eam, para_env, mm_section)
 
    END SUBROUTINE read_eam_data
 
+! **************************************************************************************************
+!> \brief reads TABPOT potential from file
+!> \param tab ...
+!> \param para_env ...
+!> \param mm_section ...
+!> \author Da Teng, Alex Mironenko
+! **************************************************************************************************
+   SUBROUTINE read_tabpot_data(tab, para_env, mm_section)
+      TYPE(tab_pot_type), POINTER                        :: tab
+      TYPE(cp_para_env_type), POINTER                    :: para_env
+      TYPE(section_vals_type), POINTER                   :: mm_section
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'read_tabpot_data', &
+         routineP = moduleN//':'//routineN
+
+      CHARACTER                                          :: d1, d2
+      INTEGER                                            :: d, handle, i, iw
+      TYPE(cp_logger_type), POINTER                      :: logger
+      TYPE(cp_parser_type), POINTER                      :: parser
+
+      CALL timeset(routineN, handle)
+      NULLIFY (parser, logger)
+      logger => cp_get_default_logger()
+      iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
+                                extension=".mmLog")
+      IF (iw > 0) WRITE (iw, *) "Reading TABPOT data from: ", TRIM(tab%tabpot_file_name)
+      CALL parser_create(parser, TRIM(tab%tabpot_file_name), para_env=para_env)
+      CALL parser_get_next_line(parser, 1)
+      IF (iw > 0) WRITE (iw, *) "Title: ", parser%input_line
+      CALL parser_get_next_line(parser, 1)
+
+      ! example format: N 1000 R 1.00 20.0
+      ! Assume the data is evenly spaced
+      READ (parser%input_line, *) d1, tab%npoints, d2, tab%dr, tab%rcut
+
+      ! Relocating arrays with the right size
+      CALL reallocate(tab%r, 1, tab%npoints)
+      CALL reallocate(tab%e, 1, tab%npoints)
+      CALL reallocate(tab%f, 1, tab%npoints)
+
+      ! Reading r, e, f
+      DO i = 1, tab%npoints
+         CALL parser_get_next_line(parser, 1)
+         READ (parser%input_line, *) d, tab%r(i), tab%e(i), tab%f(i)
+         tab%r(i) = cp_unit_to_cp2k(tab%r(i), "angstrom")
+         tab%e(i) = cp_unit_to_cp2k(tab%e(i), "kcalmol")
+         tab%f(i) = cp_unit_to_cp2k(tab%f(i), "kcalmol*angstrom^-1")
+      END DO
+
+      tab%dr = tab%r(2) - tab%r(1)
+      tab%rcut = cp_unit_to_cp2k(tab%rcut, "angstrom")
+
+      IF (iw > 0) WRITE (iw, *) "Finished TABPOT data"
+      CALL parser_release(parser)
+      CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO")
+      CALL timestop(handle)
+   END SUBROUTINE read_tabpot_data
+
 END MODULE force_fields_input
diff --git a/src/input_cp2k_mm.F b/src/input_cp2k_mm.F
@@ -60,7 +60,7 @@ MODULE input_cp2k_mm
    PUBLIC :: create_mm_section, create_dipoles_section
    PUBLIC :: create_NONBONDED14_section, create_LJ_section, create_Williams_section, &
              create_Goodwin_section, &
-             create_GENPOT_section, create_neighbor_lists_section
+             create_GENPOT_section, create_TABPOT_section, create_neighbor_lists_section
    PUBLIC :: create_CHARGE_section
 !***
 CONTAINS
@@ -1214,6 +1214,10 @@ SUBROUTINE create_NONBONDED_section(section)
       CALL section_add_subsection(section, subsection)
       CALL section_release(subsection)
 
+      CALL create_TABPOT_section(subsection)
+      CALL section_add_subsection(section, subsection)
+      CALL section_release(subsection)
+
    END SUBROUTINE create_NONBONDED_section
 
 ! **************************************************************************************************
@@ -2674,6 +2678,47 @@ SUBROUTINE create_Gal21_section(section)
 
    END SUBROUTINE create_Gal21_section
 
+! **************************************************************************************************
+!> \brief This section specifies the input parameters for TABPOT  potential type
+!> \param section the section to create
+!> \author teo, Alex Mironenko, Da Teng
+! **************************************************************************************************
+   SUBROUTINE create_TABPOT_section(section)
+
+      TYPE(section_type), POINTER                        :: section
+
+      CHARACTER(len=*), PARAMETER :: routineN = 'create_TABPOT_section', &
+         routineP = moduleN//':'//routineN
+
+      TYPE(keyword_type), POINTER                        :: keyword
+
+      CPASSERT(.NOT. ASSOCIATED(section))
+
+      CALL section_create(section, __LOCATION__, name="TABPOT", &
+                          description="This section specifies the input parameters for TABPOT potential type.", &
+                          n_keywords=1, n_subsections=0, repeats=.TRUE.)
+
+      NULLIFY (keyword)
+      CALL keyword_create(keyword, __LOCATION__, name="ATOMS", &
+                          description="Defines the atomic kind involved", &
+                          usage="ATOMS {KIND1} {KIND2}", type_of_var=char_t, &
+                          n_var=2)
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
+      CALL keyword_create(keyword, __LOCATION__, name="PARM_FILE_NAME", &
+                          variants=(/"PARMFILE"/), &
+                          description="Specifies the filename that contains the tabulated NONBONDED potential. "// &
+                          "File structure: the third line of the potential file contains a title. "// &
+                          "The 4th line contains: 'N', number of data points, 'R', lower bound of distance, distance cutoff. "// &
+                          "Follow "// &
+                          "in order npoints lines for index, distance [A], energy [kcal/mol], and force [kcal/mol/A]", &
+                          usage="PARM_FILE_NAME {FILENAME}", default_lc_val="")
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
+   END SUBROUTINE create_TABPOT_section
+
 ! **************************************************************************************************
 !> \brief This section specifies the input parameters for the subsection GCN of GAL19 and GAL21
 !>        potential type

diff --git a/src/pair_potential.F b/src/pair_potential.F
@@ -32,7 +32,7 @@ MODULE pair_potential
         b4_type, bm_type, compare_pot, ea_type, ft_type, ftd_type, gal21_type, gal_type, gp_type, &
         gw_type, ip_type, list_pot, lj_charmm_type, lj_type, multi_type, nn_type, &
         pair_potential_pp_type, pair_potential_single_type, potential_single_allocation, &
-        quip_type, siepmann_type, tersoff_type, wl_type
+        quip_type, siepmann_type, tab_type, tersoff_type, wl_type
    USE pair_potential_util,             ONLY: ener_pot,&
                                               ener_zbl,&
                                               zbl_matching_polinomial
@@ -186,7 +186,7 @@ SUBROUTINE spline_nonbond_control(spline_env, potparm, atomic_kind_set, &
             DO k = 1, SIZE(pot%type)
                SELECT CASE (pot%type(k))
                CASE (lj_type, lj_charmm_type, wl_type, gw_type, ft_type, ftd_type, ip_type, &
-                     b4_type, bm_type, gp_type, ea_type, quip_type)
+                     b4_type, bm_type, gp_type, ea_type, quip_type, tab_type)
                   pot%no_pp = .FALSE.
                CASE (tersoff_type)
                   pot%no_mb = .FALSE.
@@ -628,6 +628,8 @@ SUBROUTINE get_nonbond_storage(spline_env, potparm, atomic_kind_set, do_zbl, &
             nvar = 30 + nvar
          CASE (nn_type)
             nvar = nvar
+         CASE (tab_type)
+            nvar = 4 + nvar
          CASE DEFAULT
             CPABORT("")
          END SELECT
@@ -781,6 +783,11 @@ SUBROUTINE get_nonbond_storage(spline_env, potparm, atomic_kind_set, do_zbl, &
                         pot_par(nk, 28) = potparm%pot(i, j)%pot%set(1)%gal21%AH2
                         pot_par(nk, 29) = potparm%pot(i, j)%pot%set(1)%gal21%BH1
                         pot_par(nk, 30) = potparm%pot(i, j)%pot%set(1)%gal21%BH2
+                     CASE (tab_type)
+                        pot_par(nk, 1) = potparm%pot(i, j)%pot%set(1)%tab%dr
+                        pot_par(nk, 2) = potparm%pot(i, j)%pot%set(1)%tab%rcut
+                        pot_par(nk, 3) = potparm%pot(i, j)%pot%set(1)%tab%npoints
+                        pot_par(nk, 4) = potparm%pot(i, j)%pot%set(1)%tab%index
                      CASE (nn_type)
                         ! no checks
                      CASE DEFAULT