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