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