Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Creates and parses an input section to define a generic function
10 : !> \note
11 : !> This module handles the input definition and retrieval, while the task of actually parsing
12 : !> and evaluating a function is handled by common/fparser module. The general workflow is:
13 : !>
14 : !> 1. In the module that creates input sections, use the subroutine to add the subsection
15 : !> &GENERIC_FUNCTION with keywords FUNCTION, PARAMETERS, UNITS, VALUES and (optionally)
16 : !> VARIABLES. If a list of symbols has been specified (like input_variables=["X", "Y", "Z"])
17 : !> then they will be reserved as the variables and the VARIABLES keyword will not appear.
18 : !> If no such list is given then the VARIABLES keyword requesting input of symbol(s) for
19 : !> the variable(s) will appear, with no restriction on the number of symbols.
20 : !> If in addition to the function itself, the numerical derivatives are also needed, set
21 : !> do_deriv=.TRUE. which will create two additional keywords, DX and ERROR_LIMIT.
22 : !> NULLIFY(subsection)
23 : !> CALL create_function_section(subsection, "GENERIC_FUNCTION", "FUNCTION", &
24 : !> "This subsection defines a generic function.")
25 : !> ! Add extra keywords here if needed
26 : !> ! NULLIFY (keyword)
27 : !> ! CALL keyword_create(keyword, __LOCATION__, ...)
28 : !> ! CALL section_add_keyword(subsection, keyword)
29 : !> ! CALL keyword_release(keyword)
30 : !> ! ...
31 : !> ! Done adding extra keywords
32 : !> CALL section_add_subsection(section, subsection)
33 : !> CALL section_release(subsection)
34 : !>
35 : !> 2. Define a type that stores information about the function.
36 : !> TYPE my_func_type
37 : !> CHARACTER(LEN=default_path_length) :: function = ""
38 : !> CHARACTER(LEN=default_string_length), DIMENSION(:), ALLOCATABLE :: parameters
39 : !> REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: values
40 : !> ! To evaluate numerical derivatives also store DX and ERROR_LIMIT
41 : !> REAL(KIND=dp) :: dx = 0.0_dp, lerr = 0.0_dp
42 : !> END TYPE my_func_type
43 : !>
44 : !> 3. In the module that parses input, use the subroutine to get and save the definitions. The
45 : !> usage of input_variables should be consistent with above. Note that the variables and the
46 : !> parameters are all saved in my_func%parameters and my_func%values; with n_var variables
47 : !> and n_par parameters, both my_func%parameters and my_func%values will have a dimension of
48 : !> (n_var + n_par), variables first. To evaluate numerical derivatives, also pass
49 : !> dx=my_func%dx and error_limit=my_func%lerr.
50 : !> TYPE(my_func_type), POINTER :: my_func
51 : !> NULLIFY (my_func)
52 : !> CALL parse_function_section(gen_section, "FUNCTION", my_func%function, &
53 : !> my_func%parameters, my_func%values)
54 : !>
55 : !> 4. Use the subroutines from the fparser module to evaluate the function or its derivative.
56 : !> CALL initf(1)
57 : !> CALL parsef(1, TRIM(my_func%function), my_func%parameters)
58 : !> ! Substitute new values for the variables if needed
59 : !> ! my_func%values(1) = ...
60 : !> ! Done substituting
61 : !> result = evalf(1, my_func%values)
62 : !> CPASSERT(EvalErrType <= 0)
63 : !> ! To evaluate numerical derivative with respect to the k-th element of my_func%values:
64 : !> ! derivative = evalfd(1, k, my_func%values, my_func%dx, err)
65 : !> ! IF (ABS(err) > my_func%lerr) THEN
66 : !> ! CALL cp_warn(__LOCATION__, &
67 : !> ! "Error in computing numerical derivatives larger than allowed")
68 : !> ! END IF
69 : !> ! Done evaluating derivative
70 : !> CALL finalizef()
71 : !>
72 : !> \par History
73 : !> 10.2026 created
74 : !> \author Zilong He
75 : ! **************************************************************************************************
76 : MODULE generic_function_input
77 : USE cp_units, ONLY: cp_unit_to_cp2k
78 : USE fparser, ONLY: docf
79 : USE input_keyword_types, ONLY: keyword_create,&
80 : keyword_release,&
81 : keyword_type
82 : USE input_section_types, ONLY: section_add_keyword,&
83 : section_create,&
84 : section_type,&
85 : section_vals_type,&
86 : section_vals_val_get
87 : USE input_val_types, ONLY: char_t,&
88 : lchar_t,&
89 : real_t
90 : USE kinds, ONLY: default_path_length,&
91 : default_string_length,&
92 : dp
93 : USE memory_utilities, ONLY: reallocate
94 : USE string_utilities, ONLY: compress,&
95 : integer_to_string,&
96 : newline
97 : #include "./base/base_uses.f90"
98 :
99 : IMPLICIT NONE
100 : PRIVATE
101 :
102 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
103 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'generic_function_input'
104 :
105 : CHARACTER(len=*), PARAMETER, PRIVATE :: name_kw_dx = "DX", &
106 : name_kw_error = "ERROR_LIMIT", &
107 : name_kw_parameters = "PARAMETERS", &
108 : name_kw_units = "UNITS", &
109 : name_kw_values = "VALUES", &
110 : name_kw_variables = "VARIABLES"
111 :
112 : PUBLIC :: create_function_section, parse_function_section
113 :
114 : CONTAINS
115 :
116 : ! **************************************************************************************************
117 : !> \brief Create the input section for a generic function
118 : !> \param gen_section the section to be created for the generic function
119 : !> \param section_name the name of the subsection for the generic function
120 : !> \param func_name the name of the keyword for the function
121 : !> \param description a custom description for the subsection
122 : !> \param repeats whether the section can be repeated
123 : !> \param do_deriv whether keywords for evaluating numerical derivative are needed
124 : !> \param input_variables a list of symbols that will be parsed as variables
125 : !> \par History
126 : !> 10.2026 created from scratch matching argument names with parse_function_section
127 : ! **************************************************************************************************
128 40947 : SUBROUTINE create_function_section(gen_section, section_name, func_name, description, &
129 40947 : repeats, do_deriv, input_variables)
130 : TYPE(section_type), POINTER :: gen_section
131 : CHARACTER(LEN=*), INTENT(IN) :: section_name, func_name, description
132 : LOGICAL, INTENT(IN), OPTIONAL :: repeats, do_deriv
133 : CHARACTER(LEN=*), DIMENSION(:), INTENT(IN), &
134 : OPTIONAL :: input_variables
135 :
136 : CHARACTER(LEN=*), PARAMETER :: default_values_unit_str = "internal_cp2k"
137 :
138 40947 : CHARACTER(LEN=:), ALLOCATABLE :: input_variables_info, my_description
139 : CHARACTER(LEN=default_string_length) :: length, my_size_default, my_size_var
140 : INTEGER :: i, mydim
141 : LOGICAL :: has_kw_variables, my_do_deriv, my_repeats
142 : TYPE(keyword_type), POINTER :: keyword
143 :
144 40947 : my_repeats = .FALSE.
145 40947 : IF (PRESENT(repeats)) my_repeats = repeats
146 40947 : my_do_deriv = .FALSE.
147 40947 : IF (PRESENT(do_deriv)) my_do_deriv = do_deriv
148 40947 : has_kw_variables = .FALSE.
149 40947 : my_description = TRIM(description)
150 : mydim = 0
151 40947 : CALL integer_to_string(default_string_length, my_size_default)
152 40947 : CALL integer_to_string(default_path_length, length)
153 :
154 40947 : IF (PRESENT(input_variables)) THEN
155 24846 : has_kw_variables = .TRUE.
156 24846 : mydim = SIZE(input_variables)
157 24846 : CALL integer_to_string(mydim, my_size_var)
158 : input_variables_info = "The following "//TRIM(my_size_var)//" symbol(s) will be reserved "// &
159 24846 : "as the variable(s) to be substituted for internal evaluation: "
160 77481 : DO i = 1, SIZE(input_variables)
161 52635 : input_variables_info = input_variables_info//"`"//TRIM(input_variables(i))//"`"
162 77481 : IF (i < SIZE(input_variables)) THEN
163 27789 : input_variables_info = input_variables_info//", "
164 : ELSE
165 24846 : input_variables_info = input_variables_info//". "
166 : END IF
167 : END DO
168 24846 : input_variables_info = input_variables_info//"They cannot be used as parameters."
169 24846 : my_description = my_description//newline//newline//input_variables_info
170 : END IF
171 :
172 40947 : NULLIFY (keyword)
173 40947 : CPASSERT(.NOT. ASSOCIATED(gen_section))
174 :
175 : CALL section_create(gen_section, __LOCATION__, name=TRIM(section_name), &
176 : description=my_description, &
177 40947 : n_keywords=0, n_subsections=0, repeats=my_repeats)
178 :
179 : CALL keyword_create(keyword, __LOCATION__, name=TRIM(func_name), &
180 : description="Specifies a functional form. The maximum length "// &
181 : "is "//TRIM(ADJUSTL(length))//" characters. "//docf(), &
182 : usage=TRIM(func_name)//" A*X+B*(2/(1+EXP(-0.001*Y))-1)", &
183 40947 : type_of_var=lchar_t, n_var=1)
184 40947 : CALL section_add_keyword(gen_section, keyword)
185 40947 : CALL keyword_release(keyword)
186 :
187 : CALL keyword_create(keyword, __LOCATION__, name=name_kw_parameters, &
188 : description="Defines the symbols for the parameter(s) in the functional "// &
189 : "form, delimited by whitespaces. Symbols are case sensitive and each is "// &
190 : "at most "//TRIM(my_size_default)//" characters long.", &
191 : usage=name_kw_parameters//" A B", type_of_var=char_t, &
192 40947 : n_var=-1, repeats=.TRUE.)
193 40947 : CALL section_add_keyword(gen_section, keyword)
194 40947 : CALL keyword_release(keyword)
195 :
196 : CALL keyword_create(keyword, __LOCATION__, name=name_kw_units, &
197 : description="Optionally, defines CP2K unit string for each of the given "// &
198 : "parameter value(s). A valid CP2K unit string is one of the strings that "// &
199 : "is listed under the References/Units section of the CP2K manual, or a "// &
200 : "combination of multiple of these strings with `*`, `/` and `^`. It is "// &
201 : "assumed that the corresponding parameter value is specified in this unit "// &
202 : "which needs to be converted first to the internal unit during input parsing. "// &
203 : "Note that the meaningfulness of the unit is not checked. If omitted, all "// &
204 : "values will be treated as already in the internal unit and used literally.", &
205 : usage=name_kw_units//" eV*angstrom^-1 angstrom", type_of_var=char_t, &
206 40947 : n_var=-1, repeats=.TRUE.)
207 40947 : CALL section_add_keyword(gen_section, keyword)
208 40947 : CALL keyword_release(keyword)
209 :
210 : CALL keyword_create(keyword, __LOCATION__, name=name_kw_values, &
211 : description="Defines the value(s) of the parameter in the functional form.", &
212 : usage=name_kw_values//" 50 1.0E+2", type_of_var=real_t, &
213 40947 : n_var=-1, repeats=.TRUE., unit_str=default_values_unit_str)
214 40947 : CALL section_add_keyword(gen_section, keyword)
215 40947 : CALL keyword_release(keyword)
216 :
217 40947 : IF (.NOT. has_kw_variables) THEN
218 : CALL keyword_create(keyword, __LOCATION__, name=name_kw_variables, &
219 : description="Defines the symbols for the variable(s) in the functional "// &
220 : "form, delimited by whitespaces. Symbols are case sensitive and each is "// &
221 : "at most "//TRIM(my_size_default)//" characters long.", &
222 : usage=name_kw_variables//" X Y", type_of_var=char_t, &
223 16101 : n_var=-1)
224 16101 : CALL section_add_keyword(gen_section, keyword)
225 16101 : CALL keyword_release(keyword)
226 : END IF
227 :
228 40947 : IF (my_do_deriv) THEN
229 : CALL keyword_create(keyword, __LOCATION__, name=name_kw_dx, &
230 : description="Parameter used for computing the derivative with the Ridders' method.", &
231 11715 : usage=name_kw_dx//" <REAL>", default_r_val=0.1_dp, unit_str="bohr")
232 11715 : CALL section_add_keyword(gen_section, keyword)
233 11715 : CALL keyword_release(keyword)
234 :
235 : CALL keyword_create(keyword, __LOCATION__, name=name_kw_error, &
236 : description="Checks that the error in computing the derivative is not larger than "// &
237 : "the value set; in case error is larger a warning message is printed.", &
238 11715 : usage=name_kw_error//" <REAL>", default_r_val=1.0E-12_dp)
239 11715 : CALL section_add_keyword(gen_section, keyword)
240 11715 : CALL keyword_release(keyword)
241 : END IF
242 :
243 40947 : END SUBROUTINE create_function_section
244 :
245 : ! **************************************************************************************************
246 : !> \brief Reads from the input structure all information for generic functions
247 : !> \param gen_section the name of the section for the generic function
248 : !> \param func_name the name of the keyword for the function
249 : !> \param xfunction the function
250 : !> \param parameters the list of symbols for the parameters in the function
251 : !> \param values the list of values for the parameters in the function
252 : !> \param var_values the values of the variables in the function
253 : !> \param i_rep_sec the index of the section if it may be repeated
254 : !> \param input_variables the list of symbols for the variables in the function
255 : !> \param dx the dx for evaluating numerical derivatives of the function
256 : !> \param error_limit the error limit for evaluating numerical derivatives of the function
257 : !> \par History
258 : !> 10.2026 created by migrating the get_generic_info in force_fields_util.F to here
259 : ! **************************************************************************************************
260 4082 : SUBROUTINE parse_function_section(gen_section, func_name, xfunction, parameters, values, &
261 4082 : var_values, i_rep_sec, input_variables, dx, error_limit)
262 : TYPE(section_vals_type), POINTER :: gen_section
263 : CHARACTER(LEN=*), INTENT(IN) :: func_name
264 : CHARACTER(LEN=default_path_length), INTENT(OUT) :: xfunction
265 : CHARACTER(LEN=default_string_length), &
266 : ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: parameters
267 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
268 : INTENT(OUT) :: values
269 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: var_values
270 : INTEGER, INTENT(IN), OPTIONAL :: i_rep_sec
271 : CHARACTER(LEN=*), DIMENSION(:), INTENT(IN), &
272 : OPTIONAL :: input_variables
273 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: dx, error_limit
274 :
275 : CHARACTER(LEN=default_string_length), &
276 4082 : DIMENSION(:), POINTER :: my_par, my_par_tmp, my_units, &
277 4082 : my_units_tmp, my_var
278 : INTEGER :: i, ind, irep, isize, j, mydim, n_par, &
279 : n_units, n_val, nblank
280 : LOGICAL :: check
281 4082 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_val, my_val_tmp
282 :
283 4082 : NULLIFY (my_var, my_par, my_val, my_par_tmp, my_val_tmp)
284 4082 : NULLIFY (my_units)
285 4082 : NULLIFY (my_units_tmp)
286 4082 : irep = 1
287 3698 : IF (PRESENT(i_rep_sec)) irep = i_rep_sec
288 4082 : mydim = 0
289 4082 : CALL section_vals_val_get(gen_section, TRIM(func_name), i_rep_section=irep, c_val=xfunction)
290 4082 : CALL compress(xfunction, full=.TRUE.)
291 4082 : IF (PRESENT(dx)) THEN
292 850 : CALL section_vals_val_get(gen_section, name_kw_dx, i_rep_section=irep, r_val=dx)
293 : END IF
294 4082 : IF (PRESENT(error_limit)) THEN
295 850 : CALL section_vals_val_get(gen_section, name_kw_error, i_rep_section=irep, r_val=error_limit)
296 : END IF
297 4082 : IF (PRESENT(input_variables)) THEN
298 1410 : ALLOCATE (my_var(SIZE(input_variables)))
299 1876 : my_var = input_variables
300 : ELSE
301 3612 : CALL section_vals_val_get(gen_section, name_kw_variables, i_rep_section=irep, c_vals=my_var)
302 : END IF
303 4082 : IF (ASSOCIATED(my_var)) THEN
304 4082 : mydim = SIZE(my_var)
305 : END IF
306 : ! Check for presence of Parameters
307 4082 : CALL section_vals_val_get(gen_section, name_kw_parameters, i_rep_section=irep, n_rep_val=n_par)
308 4082 : CALL section_vals_val_get(gen_section, name_kw_values, i_rep_section=irep, n_rep_val=n_val)
309 4082 : check = (n_par > 0) .EQV. (n_val > 0)
310 4082 : CPASSERT(check)
311 4082 : CALL section_vals_val_get(gen_section, name_kw_units, i_rep_section=irep, n_rep_val=n_units)
312 4082 : IF (n_par > 0) THEN
313 : ! Parameters
314 3458 : ALLOCATE (my_par(0))
315 3458 : ALLOCATE (my_val(0))
316 6916 : DO i = 1, n_par
317 3458 : isize = SIZE(my_par)
318 3458 : CALL section_vals_val_get(gen_section, name_kw_parameters, i_rep_section=irep, i_rep_val=i, c_vals=my_par_tmp)
319 10278 : nblank = COUNT(my_par_tmp == "")
320 3458 : CALL reallocate(my_par, 1, isize + SIZE(my_par_tmp) - nblank)
321 3458 : ind = 0
322 13736 : DO j = 1, SIZE(my_par_tmp)
323 6820 : IF (my_par_tmp(j) == "") CYCLE
324 6818 : ind = ind + 1
325 10278 : my_par(isize + ind) = my_par_tmp(j)
326 : END DO
327 : END DO
328 6918 : DO i = 1, n_val
329 3460 : isize = SIZE(my_val)
330 3460 : CALL section_vals_val_get(gen_section, name_kw_values, i_rep_section=irep, i_rep_val=i, r_vals=my_val_tmp)
331 3460 : CALL reallocate(my_val, 1, isize + SIZE(my_val_tmp))
332 20554 : my_val(isize + 1:isize + SIZE(my_val_tmp)) = my_val_tmp
333 : END DO
334 3458 : CPASSERT(SIZE(my_par) == SIZE(my_val))
335 : ! Optionally read the units for each parameter value
336 3458 : ALLOCATE (my_units(0))
337 3458 : IF (n_units > 0) THEN
338 6 : DO i = 1, n_units
339 4 : isize = SIZE(my_units)
340 4 : CALL section_vals_val_get(gen_section, name_kw_units, i_rep_section=irep, i_rep_val=i, c_vals=my_units_tmp)
341 10 : nblank = COUNT(my_units_tmp == "")
342 4 : CALL reallocate(my_units, 1, isize + SIZE(my_units_tmp) - nblank)
343 4 : ind = 0
344 12 : DO j = 1, SIZE(my_units_tmp)
345 6 : IF (my_units_tmp(j) == "") CYCLE
346 6 : ind = ind + 1
347 10 : my_units(isize + ind) = my_units_tmp(j)
348 : END DO
349 : END DO
350 2 : CPASSERT(SIZE(my_units) == SIZE(my_val))
351 : END IF
352 3458 : mydim = mydim + SIZE(my_val)
353 3458 : IF (SIZE(my_val) == 0) THEN
354 2 : DEALLOCATE (my_par)
355 2 : DEALLOCATE (my_val)
356 2 : DEALLOCATE (my_units)
357 : END IF
358 : END IF
359 : ! Handle trivial case of a null function defined
360 12246 : ALLOCATE (parameters(mydim))
361 12246 : ALLOCATE (values(mydim))
362 4082 : IF (mydim > 0) THEN
363 9552 : parameters(1:SIZE(my_var)) = my_var
364 9552 : values(1:SIZE(my_var)) = 0.0_dp
365 4082 : IF (PRESENT(var_values)) THEN
366 384 : CPASSERT(SIZE(var_values) == SIZE(my_var))
367 1220 : values(1:SIZE(my_var)) = var_values
368 : END IF
369 4082 : IF (ASSOCIATED(my_val)) THEN
370 10274 : DO i = 1, SIZE(my_val)
371 6818 : parameters(SIZE(my_var) + i) = my_par(i)
372 10274 : IF (n_units > 0) THEN
373 6 : values(SIZE(my_var) + i) = cp_unit_to_cp2k(my_val(i), TRIM(ADJUSTL(my_units(i))))
374 : ELSE
375 6812 : values(SIZE(my_var) + i) = my_val(i)
376 : END IF
377 : END DO
378 : END IF
379 : END IF
380 4082 : IF (ASSOCIATED(my_par)) THEN
381 3456 : DEALLOCATE (my_par)
382 : END IF
383 4082 : IF (ASSOCIATED(my_val)) THEN
384 3456 : DEALLOCATE (my_val)
385 : END IF
386 4082 : IF (ASSOCIATED(my_units)) THEN
387 3456 : DEALLOCATE (my_units)
388 : END IF
389 4082 : IF (PRESENT(input_variables)) THEN
390 470 : DEALLOCATE (my_var)
391 : END IF
392 12246 : END SUBROUTINE parse_function_section
393 :
394 154962 : END MODULE generic_function_input
|
