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