(git:374b731)
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-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \par History
10!> Subroutine input_torsions changed (DG) 05-Dec-2000
11!> Output formats changed (DG) 05-Dec-2000
12!> JGH (26-01-2002) : force field parameters stored in tables, not in
13!> matrices. Input changed to have parameters labeled by the position
14!> and not atom pairs (triples etc)
15!> Teo (11.2005) : Moved all information on force field pair_potential to
16!> a much lighter memory structure
17!> Teo 09.2006 : Split all routines force_field I/O in a separate file
18!> \author CJM
19! **************************************************************************************************
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
779 n_items = 1
780 isec = 1
781 n_items = isec*n_items
782 CALL section_vals_val_get(section, "ATOMS", c_vals=atm_names)
783 CALL section_vals_val_get(section, "PARM_FILE_NAME", c_val=nequip_file_name)
784 CALL section_vals_val_get(section, "UNIT_COORDS", c_val=unit_coords)
785 CALL section_vals_val_get(section, "UNIT_ENERGY", c_val=unit_energy)
786 CALL section_vals_val_get(section, "UNIT_FORCES", c_val=unit_forces)
787 CALL section_vals_val_get(section, "UNIT_CELL", c_val=unit_cell)
788
789 DO isec = 1, SIZE(atm_names)
790 DO jsec = isec, SIZE(atm_names)
791 nonbonded%pot(start + n_items)%pot%type = nequip_type
792 nonbonded%pot(start + n_items)%pot%at1 = atm_names(isec)
793 nonbonded%pot(start + n_items)%pot%at2 = atm_names(jsec)
794 CALL uppercase(nonbonded%pot(start + n_items)%pot%at1)
795 CALL uppercase(nonbonded%pot(start + n_items)%pot%at2)
796 nonbonded%pot(start + n_items)%pot%set(1)%nequip%nequip_file_name = discover_file(nequip_file_name)
797 nonbonded%pot(start + n_items)%pot%set(1)%nequip%unit_coords = unit_coords
798 nonbonded%pot(start + n_items)%pot%set(1)%nequip%unit_forces = unit_forces
799 nonbonded%pot(start + n_items)%pot%set(1)%nequip%unit_energy = unit_energy
800 nonbonded%pot(start + n_items)%pot%set(1)%nequip%unit_cell = unit_cell
801 CALL read_nequip_data(nonbonded%pot(start + n_items)%pot%set(1)%nequip)
802 nonbonded%pot(start + n_items)%pot%rcutsq = nonbonded%pot(start + n_items)%pot%set(1)%nequip%rcutsq
803 n_items = n_items + 1
804 END DO
805 END DO
806
807 END SUBROUTINE read_nequip_section
808
809! **************************************************************************************************
810!> \brief Reads the ALLEGRO section
811!> \param nonbonded ...
812!> \param section ...
813!> \param start ...
814!> \author Gabriele Tocci
815! **************************************************************************************************
816 SUBROUTINE read_allegro_section(nonbonded, section, start)
817 TYPE(pair_potential_p_type), POINTER :: nonbonded
818 TYPE(section_vals_type), POINTER :: section
819 INTEGER, INTENT(IN) :: start
820
821 CHARACTER(LEN=default_string_length) :: allegro_file_name, unit_cell, &
822 unit_coords, unit_energy, unit_forces
823 CHARACTER(LEN=default_string_length), &
824 DIMENSION(:), POINTER :: atm_names
825 INTEGER :: isec, jsec, n_items
826
827 n_items = 1
828 isec = 1
829 n_items = isec*n_items
830 CALL section_vals_val_get(section, "ATOMS", c_vals=atm_names)
831 CALL section_vals_val_get(section, "PARM_FILE_NAME", c_val=allegro_file_name)
832 CALL section_vals_val_get(section, "UNIT_COORDS", c_val=unit_coords)
833 CALL section_vals_val_get(section, "UNIT_ENERGY", c_val=unit_energy)
834 CALL section_vals_val_get(section, "UNIT_FORCES", c_val=unit_forces)
835 CALL section_vals_val_get(section, "UNIT_CELL", c_val=unit_cell)
836
837 DO isec = 1, SIZE(atm_names)
838 DO jsec = isec, SIZE(atm_names)
839 nonbonded%pot(start + n_items)%pot%type = allegro_type
840 nonbonded%pot(start + n_items)%pot%at1 = atm_names(isec)
841 nonbonded%pot(start + n_items)%pot%at2 = atm_names(jsec)
842 CALL uppercase(nonbonded%pot(start + n_items)%pot%at1)
843 CALL uppercase(nonbonded%pot(start + n_items)%pot%at2)
844 nonbonded%pot(start + n_items)%pot%set(1)%allegro%allegro_file_name = discover_file(allegro_file_name)
845 nonbonded%pot(start + n_items)%pot%set(1)%allegro%unit_coords = unit_coords
846 nonbonded%pot(start + n_items)%pot%set(1)%allegro%unit_forces = unit_forces
847 nonbonded%pot(start + n_items)%pot%set(1)%allegro%unit_energy = unit_energy
848 nonbonded%pot(start + n_items)%pot%set(1)%allegro%unit_cell = unit_cell
849 CALL read_allegro_data(nonbonded%pot(start + n_items)%pot%set(1)%allegro)
850 nonbonded%pot(start + n_items)%pot%rcutsq = nonbonded%pot(start + n_items)%pot%set(1)%allegro%rcutsq
851 n_items = n_items + 1
852 END DO
853 END DO
854
855 END SUBROUTINE read_allegro_section
856
857! **************************************************************************************************
858!> \brief Reads the LJ section
859!> \param nonbonded ...
860!> \param section ...
861!> \param start ...
862!> \author teo
863! **************************************************************************************************
864 SUBROUTINE read_lj_section(nonbonded, section, start)
865 TYPE(pair_potential_p_type), POINTER :: nonbonded
866 TYPE(section_vals_type), POINTER :: section
867 INTEGER, INTENT(IN) :: start
868
869 CHARACTER(LEN=default_string_length), &
870 DIMENSION(:), POINTER :: atm_names
871 INTEGER :: isec, n_items, n_rep
872 REAL(kind=dp) :: epsilon, rcut, sigma
873
874 CALL section_vals_get(section, n_repetition=n_items)
875 DO isec = 1, n_items
876 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
877 CALL section_vals_val_get(section, "EPSILON", i_rep_section=isec, r_val=epsilon)
878 CALL section_vals_val_get(section, "SIGMA", i_rep_section=isec, r_val=sigma)
879 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
880
881 nonbonded%pot(start + isec)%pot%type = lj_charmm_type
882 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
883 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
884 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
885 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
886 nonbonded%pot(start + isec)%pot%set(1)%lj%epsilon = epsilon
887 nonbonded%pot(start + isec)%pot%set(1)%lj%sigma6 = sigma**6
888 nonbonded%pot(start + isec)%pot%set(1)%lj%sigma12 = sigma**12
889 nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut
890 !
891 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
892 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
893 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
894 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
895 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
896 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
897 END DO
898 END SUBROUTINE read_lj_section
899
900! **************************************************************************************************
901!> \brief Reads the WILLIAMS section
902!> \param nonbonded ...
903!> \param section ...
904!> \param start ...
905!> \author teo
906! **************************************************************************************************
907 SUBROUTINE read_wl_section(nonbonded, section, start)
908 TYPE(pair_potential_p_type), POINTER :: nonbonded
909 TYPE(section_vals_type), POINTER :: section
910 INTEGER, INTENT(IN) :: start
911
912 CHARACTER(LEN=default_string_length), &
913 DIMENSION(:), POINTER :: atm_names
914 INTEGER :: isec, n_items, n_rep
915 REAL(kind=dp) :: a, b, c, rcut
916
917 CALL section_vals_get(section, n_repetition=n_items)
918 DO isec = 1, n_items
919 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
920 CALL section_vals_val_get(section, "A", i_rep_section=isec, r_val=a)
921 CALL section_vals_val_get(section, "B", i_rep_section=isec, r_val=b)
922 CALL section_vals_val_get(section, "C", i_rep_section=isec, r_val=c)
923 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
924
925 nonbonded%pot(start + isec)%pot%type = wl_type
926 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
927 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
928 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
929 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
930 nonbonded%pot(start + isec)%pot%set(1)%willis%a = a
931 nonbonded%pot(start + isec)%pot%set(1)%willis%b = b
932 nonbonded%pot(start + isec)%pot%set(1)%willis%c = c
933 nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut
934 !
935 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
936 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
937 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
938 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
939 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
940 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
941 END DO
942 END SUBROUTINE read_wl_section
943
944! **************************************************************************************************
945!> \brief Reads the GOODWIN section
946!> \param nonbonded ...
947!> \param section ...
948!> \param start ...
949!> \author teo
950! **************************************************************************************************
951 SUBROUTINE read_gd_section(nonbonded, section, start)
952 TYPE(pair_potential_p_type), POINTER :: nonbonded
953 TYPE(section_vals_type), POINTER :: section
954 INTEGER, INTENT(IN) :: start
955
956 CHARACTER(LEN=default_string_length), &
957 DIMENSION(:), POINTER :: atm_names
958 INTEGER :: isec, m, mc, n_items, n_rep
959 REAL(kind=dp) :: d, dc, rcut, vr0
960
961 CALL section_vals_get(section, n_repetition=n_items)
962 DO isec = 1, n_items
963 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
964 CALL section_vals_val_get(section, "VR0", i_rep_section=isec, r_val=vr0)
965 CALL section_vals_val_get(section, "D", i_rep_section=isec, r_val=d)
966 CALL section_vals_val_get(section, "DC", i_rep_section=isec, r_val=dc)
967 CALL section_vals_val_get(section, "M", i_rep_section=isec, i_val=m)
968 CALL section_vals_val_get(section, "MC", i_rep_section=isec, i_val=mc)
969 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
970
971 nonbonded%pot(start + isec)%pot%type = gw_type
972 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
973 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
974 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
975 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
976 nonbonded%pot(start + isec)%pot%set(1)%goodwin%vr0 = vr0
977 nonbonded%pot(start + isec)%pot%set(1)%goodwin%d = d
978 nonbonded%pot(start + isec)%pot%set(1)%goodwin%dc = dc
979 nonbonded%pot(start + isec)%pot%set(1)%goodwin%m = m
980 nonbonded%pot(start + isec)%pot%set(1)%goodwin%mc = mc
981 nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut
982 !
983 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
984 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
985 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
986 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
987 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
988 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
989 END DO
990 END SUBROUTINE read_gd_section
991
992! **************************************************************************************************
993!> \brief Reads the IPBV section
994!> \param nonbonded ...
995!> \param section ...
996!> \param start ...
997!> \author teo
998! **************************************************************************************************
999 SUBROUTINE read_ipbv_section(nonbonded, section, start)
1000 TYPE(pair_potential_p_type), POINTER :: nonbonded
1001 TYPE(section_vals_type), POINTER :: section
1002 INTEGER, INTENT(IN) :: start
1003
1004 CHARACTER(LEN=default_string_length), &
1005 DIMENSION(:), POINTER :: atm_names
1006 INTEGER :: isec, n_items, n_rep
1007 REAL(kind=dp) :: rcut
1008
1009 CALL section_vals_get(section, n_repetition=n_items)
1010 DO isec = 1, n_items
1011 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
1012 nonbonded%pot(start + isec)%pot%type = ip_type
1013 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
1014 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
1015 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
1016 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
1017 CALL set_ipbv_ff(nonbonded%pot(start + isec)%pot%at1, nonbonded%pot(start + isec)%pot%at2, &
1018 nonbonded%pot(start + isec)%pot%set(1)%ipbv)
1019 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
1020 nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
1021 !
1022 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
1023 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
1024 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
1025 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
1026 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
1027 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
1028 END DO
1029 END SUBROUTINE read_ipbv_section
1030
1031! **************************************************************************************************
1032!> \brief Reads the BMHFT section
1033!> \param nonbonded ...
1034!> \param section ...
1035!> \param start ...
1036!> \author teo
1037! **************************************************************************************************
1038 SUBROUTINE read_bmhft_section(nonbonded, section, start)
1039 TYPE(pair_potential_p_type), POINTER :: nonbonded
1040 TYPE(section_vals_type), POINTER :: section
1041 INTEGER, INTENT(IN) :: start
1042
1043 CHARACTER(LEN=default_string_length), DIMENSION(2) :: map_atoms
1044 CHARACTER(LEN=default_string_length), &
1045 DIMENSION(:), POINTER :: atm_names
1046 INTEGER :: i, isec, n_items, n_rep
1047 REAL(kind=dp) :: rcut
1048
1049 CALL section_vals_get(section, n_repetition=n_items)
1050 DO isec = 1, n_items
1051 CALL cite_reference(tosi1964a)
1052 CALL cite_reference(tosi1964b)
1053 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
1054 nonbonded%pot(start + isec)%pot%type = ft_type
1055 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
1056 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
1057 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
1058 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
1059
1060 CALL section_vals_val_get(section, "A", i_rep_section=isec, n_rep_val=i)
1061 IF (i == 1) THEN
1062 CALL section_vals_val_get(section, "A", i_rep_section=isec, &
1063 r_val=nonbonded%pot(start + isec)%pot%set(1)%ft%a)
1064 CALL section_vals_val_get(section, "B", i_rep_section=isec, &
1065 r_val=nonbonded%pot(start + isec)%pot%set(1)%ft%b)
1066 CALL section_vals_val_get(section, "C", i_rep_section=isec, &
1067 r_val=nonbonded%pot(start + isec)%pot%set(1)%ft%c)
1068 CALL section_vals_val_get(section, "D", i_rep_section=isec, &
1069 r_val=nonbonded%pot(start + isec)%pot%set(1)%ft%d)
1070 ELSE
1071 CALL section_vals_val_get(section, "MAP_ATOMS", i_rep_section=isec, c_vals=atm_names)
1072 map_atoms = atm_names
1073 CALL uppercase(map_atoms(1))
1074 CALL uppercase(map_atoms(2))
1075 CALL set_bmhft_ff(map_atoms(1), map_atoms(2), nonbonded%pot(start + isec)%pot%set(1)%ft)
1076 END IF
1077 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
1078 nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
1079 !
1080 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
1081 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
1082 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
1083 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
1084 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
1085 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
1086 END DO
1087 END SUBROUTINE read_bmhft_section
1088
1089! **************************************************************************************************
1090!> \brief Reads the BMHFTD section
1091!> \param nonbonded ...
1092!> \param section ...
1093!> \param start ...
1094!> \author Mathieu Salanne 05.2010
1095! **************************************************************************************************
1096 SUBROUTINE read_bmhftd_section(nonbonded, section, start)
1097 TYPE(pair_potential_p_type), POINTER :: nonbonded
1098 TYPE(section_vals_type), POINTER :: section
1099 INTEGER, INTENT(IN) :: start
1100
1101 CHARACTER(LEN=default_string_length), DIMENSION(2) :: map_atoms
1102 CHARACTER(LEN=default_string_length), &
1103 DIMENSION(:), POINTER :: atm_names
1104 INTEGER :: i, isec, n_items, n_rep
1105 REAL(kind=dp) :: rcut
1106 REAL(kind=dp), DIMENSION(:), POINTER :: bd_vals
1107
1108 NULLIFY (bd_vals)
1109
1110 CALL section_vals_get(section, n_repetition=n_items)
1111
1112 DO isec = 1, n_items
1113 CALL cite_reference(tosi1964a)
1114 CALL cite_reference(tosi1964b)
1115 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
1116 nonbonded%pot(start + isec)%pot%type = ftd_type
1117 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
1118 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
1119 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
1120 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
1121
1122 CALL section_vals_val_get(section, "A", i_rep_section=isec, n_rep_val=i)
1123 IF (i == 1) THEN
1124 CALL section_vals_val_get(section, "A", i_rep_section=isec, &
1125 r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%a)
1126 CALL section_vals_val_get(section, "B", i_rep_section=isec, &
1127 r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%b)
1128 CALL section_vals_val_get(section, "C", i_rep_section=isec, &
1129 r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%c)
1130 CALL section_vals_val_get(section, "D", i_rep_section=isec, &
1131 r_val=nonbonded%pot(start + isec)%pot%set(1)%ftd%d)
1132 CALL section_vals_val_get(section, "BD", i_rep_section=isec, r_vals=bd_vals)
1133 IF (ASSOCIATED(bd_vals)) THEN
1134 SELECT CASE (SIZE(bd_vals))
1135 CASE (0)
1136 cpabort("No values specified for parameter BD in section &BMHFTD")
1137 CASE (1)
1138 nonbonded%pot(start + isec)%pot%set(1)%ftd%bd(1:2) = bd_vals(1)
1139 CASE (2)
1140 nonbonded%pot(start + isec)%pot%set(1)%ftd%bd(1:2) = bd_vals(1:2)
1141 CASE DEFAULT
1142 cpabort("Too many values specified for parameter BD in section &BMHFTD")
1143 END SELECT
1144 ELSE
1145 cpabort("Parameter BD in section &BMHFTD was not specified")
1146 END IF
1147 ELSE
1148 CALL section_vals_val_get(section, "MAP_ATOMS", i_rep_section=isec, c_vals=atm_names)
1149 map_atoms = atm_names
1150 CALL uppercase(map_atoms(1))
1151 CALL uppercase(map_atoms(2))
1152 CALL set_bmhftd_ff()
1153 END IF
1154 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
1155 nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
1156 !
1157 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
1158 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
1159 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
1160 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
1161 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
1162 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
1163 END DO
1164 END SUBROUTINE read_bmhftd_section
1165
1166! **************************************************************************************************
1167!> \brief Reads the Buckingham 4 Ranges potential section
1168!> \param nonbonded ...
1169!> \param section ...
1170!> \param start ...
1171!> \par History
1172!> MK (11.11.2010): Automatic fit of the (default) polynomial coefficients
1173!> \author MI,MK
1174! **************************************************************************************************
1175 SUBROUTINE read_b4_section(nonbonded, section, start)
1176
1177 TYPE(pair_potential_p_type), POINTER :: nonbonded
1178 TYPE(section_vals_type), POINTER :: section
1179 INTEGER, INTENT(IN) :: start
1180
1181 CHARACTER(LEN=default_string_length), &
1182 DIMENSION(:), POINTER :: atm_names
1183 INTEGER :: i, ir, isec, n_items, n_rep, np1, np2
1184 LOGICAL :: explicit_poly1, explicit_poly2
1185 REAL(kind=dp) :: a, b, c, eval_error, r1, r2, r3, rcut
1186 REAL(kind=dp), DIMENSION(10) :: v, x
1187 REAL(kind=dp), DIMENSION(10, 10) :: p, p_inv
1188 REAL(kind=dp), DIMENSION(:), POINTER :: coeff1, coeff2, list
1189
1190 NULLIFY (coeff1)
1191 NULLIFY (coeff2)
1192
1193 CALL section_vals_get(section, n_repetition=n_items)
1194
1195 DO isec = 1, n_items
1196 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
1197 CALL section_vals_val_get(section, "A", i_rep_section=isec, r_val=a)
1198 CALL section_vals_val_get(section, "B", i_rep_section=isec, r_val=b)
1199 CALL section_vals_val_get(section, "C", i_rep_section=isec, r_val=c)
1200 CALL section_vals_val_get(section, "R1", i_rep_section=isec, r_val=r1)
1201 CALL section_vals_val_get(section, "R2", i_rep_section=isec, r_val=r2)
1202 CALL section_vals_val_get(section, "R3", i_rep_section=isec, r_val=r3)
1203 CALL section_vals_val_get(section, "POLY1", explicit=explicit_poly1, n_rep_val=n_rep)
1204 ! Check if polynomial coefficients were specified for range 2 and 3 explicitly
1205 IF (explicit_poly1) THEN
1206 np1 = 0
1207 DO ir = 1, n_rep
1208 NULLIFY (list)
1209 CALL section_vals_val_get(section, "POLY1", i_rep_val=ir, r_vals=list)
1210 IF (ASSOCIATED(list)) THEN
1211 CALL reallocate(coeff1, 0, np1 + SIZE(list) - 1)
1212 DO i = 1, SIZE(list)
1213 coeff1(i + np1 - 1) = list(i)
1214 END DO
1215 np1 = np1 + SIZE(list)
1216 END IF
1217 END DO
1218 END IF
1219 CALL section_vals_val_get(section, "POLY2", explicit=explicit_poly2, n_rep_val=n_rep)
1220 IF (explicit_poly2) THEN
1221 np2 = 0
1222 DO ir = 1, n_rep
1223 NULLIFY (list)
1224 CALL section_vals_val_get(section, "POLY2", i_rep_val=ir, r_vals=list)
1225 IF (ASSOCIATED(list)) THEN
1226 CALL reallocate(coeff2, 0, np2 + SIZE(list) - 1)
1227 DO i = 1, SIZE(list)
1228 coeff2(i + np2 - 1) = list(i)
1229 END DO
1230 np2 = np2 + SIZE(list)
1231 END IF
1232 END DO
1233 END IF
1234 ! Default is a 5th/3rd-order polynomial fit
1235 IF ((.NOT. explicit_poly1) .OR. (.NOT. explicit_poly2)) THEN
1236 ! Build matrix p and vector v to calculate the polynomial coefficients
1237 ! in the vector x from p*x = v
1238 p(:, :) = 0.0_dp
1239 ! Row 1: Match the 5th-order polynomial and the potential at r1
1240 p(1, 1) = 1.0_dp
1241 DO i = 2, 6
1242 p(1, i) = p(1, i - 1)*r1
1243 END DO
1244 ! Row 2: Match the first derivatives of the 5th-order polynomial and the potential at r1
1245 DO i = 2, 6
1246 p(2, i) = real(i - 1, kind=dp)*p(1, i - 1)
1247 END DO
1248 ! Row 3: Match the second derivatives of the 5th-order polynomial and the potential at r1
1249 DO i = 3, 6
1250 p(3, i) = real(i - 1, kind=dp)*p(2, i - 1)
1251 END DO
1252 ! Row 4: Match the 5th-order and the 3rd-order polynomials at r2
1253 p(4, 1) = 1.0_dp
1254 DO i = 2, 6
1255 p(4, i) = p(4, i - 1)*r2
1256 END DO
1257 p(4, 7) = -1.0_dp
1258 DO i = 8, 10
1259 p(4, i) = p(4, i - 1)*r2
1260 END DO
1261 ! Row 5: Match the first derivatives of the 5th-order and the 3rd-order polynomials at r2
1262 DO i = 2, 6
1263 p(5, i) = real(i - 1, kind=dp)*p(4, i - 1)
1264 END DO
1265 DO i = 8, 10
1266 p(5, i) = real(i - 7, kind=dp)*p(4, i - 1)
1267 END DO
1268 ! Row 6: Match the second derivatives of the 5th-order and the 3rd-order polynomials at r2
1269 DO i = 3, 6
1270 p(6, i) = real(i - 1, kind=dp)*p(5, i - 1)
1271 END DO
1272 DO i = 9, 10
1273 p(6, i) = real(i - 7, kind=dp)*p(5, i - 1)
1274 END DO
1275 ! Row 7: Minimum at r2, ie. the first derivative of the 3rd-order polynomial has to be zero at r2
1276 DO i = 8, 10
1277 p(7, i) = -p(5, i)
1278 END DO
1279 ! Row 8: Match the 3rd-order polynomial and the potential at r3
1280 p(8, 7) = 1.0_dp
1281 DO i = 8, 10
1282 p(8, i) = p(8, i - 1)*r3
1283 END DO
1284 ! Row 9: Match the first derivatives of the 3rd-order polynomial and the potential at r3
1285 DO i = 8, 10
1286 p(9, i) = real(i - 7, kind=dp)*p(8, i - 1)
1287 END DO
1288 ! Row 10: Match the second derivatives of the 3rd-order polynomial and the potential at r3
1289 DO i = 9, 10
1290 p(10, i) = real(i - 7, kind=dp)*p(9, i - 1)
1291 END DO
1292 ! Build the vector v
1293 v(1) = a*exp(-b*r1)
1294 v(2) = -b*v(1)
1295 v(3) = -b*v(2)
1296 v(4:7) = 0.0_dp
1297 v(8) = -c/p(8, 10)**2 ! = -c/r3**6
1298 v(9) = -6.0_dp*v(8)/r3
1299 v(10) = -7.0_dp*v(9)/r3
1300 ! Calculate p_inv the inverse of the matrix p
1301 p_inv(:, :) = 0.0_dp
1302 CALL invert_matrix(p, p_inv, eval_error)
1303 IF (eval_error >= 1.0e-8_dp) &
1304 CALL cp_warn(__location__, &
1305 "The polynomial fit for the BUCK4RANGES potential is only accurate to "// &
1306 trim(cp_to_string(eval_error)))
1307 ! Get the 6 coefficients of the 5th-order polynomial -> x(1:6)
1308 ! and the 4 coefficients of the 3rd-order polynomial -> x(7:10)
1309 x(:) = matmul(p_inv(:, :), v(:))
1310 ELSE
1311 x(:) = 0.0_dp
1312 END IF
1313
1314 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
1315
1316 nonbonded%pot(start + isec)%pot%type = b4_type
1317 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
1318 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
1319 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
1320 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
1321 nonbonded%pot(start + isec)%pot%set(1)%buck4r%a = a
1322 nonbonded%pot(start + isec)%pot%set(1)%buck4r%b = b
1323 nonbonded%pot(start + isec)%pot%set(1)%buck4r%c = c
1324 nonbonded%pot(start + isec)%pot%set(1)%buck4r%r1 = r1
1325 nonbonded%pot(start + isec)%pot%set(1)%buck4r%r2 = r2
1326 nonbonded%pot(start + isec)%pot%set(1)%buck4r%r3 = r3
1327 IF ((.NOT. explicit_poly1) .OR. (.NOT. explicit_poly2)) THEN
1328 nonbonded%pot(start + isec)%pot%set(1)%buck4r%npoly1 = 5
1329 nonbonded%pot(start + isec)%pot%set(1)%buck4r%poly1(0:5) = x(1:6)
1330 nonbonded%pot(start + isec)%pot%set(1)%buck4r%npoly2 = 3
1331 nonbonded%pot(start + isec)%pot%set(1)%buck4r%poly2(0:3) = x(7:10)
1332 ELSE
1333 nonbonded%pot(start + isec)%pot%set(1)%buck4r%npoly1 = np1 - 1
1334 cpassert(np1 - 1 <= 10)
1335 nonbonded%pot(start + isec)%pot%set(1)%buck4r%poly1(0:np1 - 1) = coeff1(0:np1 - 1)
1336 nonbonded%pot(start + isec)%pot%set(1)%buck4r%npoly2 = np2 - 1
1337 cpassert(np2 - 1 <= 10)
1338 nonbonded%pot(start + isec)%pot%set(1)%buck4r%poly2(0:np2 - 1) = coeff2(0:np2 - 1)
1339 END IF
1340 nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut
1341
1342 IF (ASSOCIATED(coeff1)) THEN
1343 DEALLOCATE (coeff1)
1344 END IF
1345 IF (ASSOCIATED(coeff2)) THEN
1346 DEALLOCATE (coeff2)
1347 END IF
1348 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
1349 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
1350 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
1351 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
1352 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
1353 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
1354 END DO
1355
1356 END SUBROUTINE read_b4_section
1357
1358! **************************************************************************************************
1359!> \brief Reads the GENPOT - generic potential section
1360!> \param nonbonded ...
1361!> \param section ...
1362!> \param start ...
1363!> \author Teodoro Laino - 10.2006
1364! **************************************************************************************************
1365 SUBROUTINE read_gp_section(nonbonded, section, start)
1366 TYPE(pair_potential_p_type), POINTER :: nonbonded
1367 TYPE(section_vals_type), POINTER :: section
1368 INTEGER, INTENT(IN) :: start
1369
1370 CHARACTER(LEN=default_string_length), &
1371 DIMENSION(:), POINTER :: atm_names
1372 INTEGER :: isec, n_items, n_rep
1373 REAL(kind=dp) :: rcut
1374
1375 CALL section_vals_get(section, n_repetition=n_items)
1376 DO isec = 1, n_items
1377 NULLIFY (atm_names)
1378 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
1379 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
1380 nonbonded%pot(start + isec)%pot%type = gp_type
1381 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
1382 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
1383 nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut
1384 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
1385 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
1386 ! Parse the genpot info
1387 CALL get_generic_info(section, "FUNCTION", nonbonded%pot(start + isec)%pot%set(1)%gp%potential, &
1388 nonbonded%pot(start + isec)%pot%set(1)%gp%parameters, &
1389 nonbonded%pot(start + isec)%pot%set(1)%gp%values, &
1390 size_variables=1, i_rep_sec=isec)
1391 nonbonded%pot(start + isec)%pot%set(1)%gp%variables = nonbonded%pot(start + isec)%pot%set(1)%gp%parameters(1)
1392 !
1393 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
1394 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
1395 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
1396 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
1397 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
1398 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
1399 END DO
1400 END SUBROUTINE read_gp_section
1401
1402! **************************************************************************************************
1403!> \brief Reads the tersoff section
1404!> \param nonbonded ...
1405!> \param section ...
1406!> \param start ...
1407!> \param tersoff_section ...
1408!> \author ikuo
1409! **************************************************************************************************
1410 SUBROUTINE read_tersoff_section(nonbonded, section, start, tersoff_section)
1411 TYPE(pair_potential_p_type), POINTER :: nonbonded
1412 TYPE(section_vals_type), POINTER :: section
1413 INTEGER, INTENT(IN) :: start
1414 TYPE(section_vals_type), POINTER :: tersoff_section
1415
1416 CHARACTER(LEN=default_string_length), &
1417 DIMENSION(:), POINTER :: atm_names
1418 INTEGER :: isec, n_items, n_rep
1419 REAL(kind=dp) :: rcut, rcutsq
1420
1421 CALL section_vals_get(section, n_repetition=n_items)
1422 DO isec = 1, n_items
1423 CALL cite_reference(tersoff1988)
1424 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
1425
1426 nonbonded%pot(start + isec)%pot%type = tersoff_type
1427 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
1428 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
1429 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
1430 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
1431
1432 CALL section_vals_val_get(tersoff_section, "A", i_rep_section=isec, &
1433 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%A)
1434 CALL section_vals_val_get(tersoff_section, "B", i_rep_section=isec, &
1435 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%B)
1436 CALL section_vals_val_get(tersoff_section, "lambda1", i_rep_section=isec, &
1437 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%lambda1)
1438 CALL section_vals_val_get(tersoff_section, "lambda2", i_rep_section=isec, &
1439 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%lambda2)
1440 CALL section_vals_val_get(tersoff_section, "alpha", i_rep_section=isec, &
1441 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%alpha)
1442 CALL section_vals_val_get(tersoff_section, "beta", i_rep_section=isec, &
1443 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%beta)
1444 CALL section_vals_val_get(tersoff_section, "n", i_rep_section=isec, &
1445 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%n)
1446 CALL section_vals_val_get(tersoff_section, "c", i_rep_section=isec, &
1447 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%c)
1448 CALL section_vals_val_get(tersoff_section, "d", i_rep_section=isec, &
1449 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%d)
1450 CALL section_vals_val_get(tersoff_section, "h", i_rep_section=isec, &
1451 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%h)
1452 CALL section_vals_val_get(tersoff_section, "lambda3", i_rep_section=isec, &
1453 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%lambda3)
1454 CALL section_vals_val_get(tersoff_section, "bigR", i_rep_section=isec, &
1455 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%bigR)
1456 CALL section_vals_val_get(tersoff_section, "bigD", i_rep_section=isec, &
1457 r_val=nonbonded%pot(start + isec)%pot%set(1)%tersoff%bigD)
1458
1459 rcutsq = (nonbonded%pot(start + isec)%pot%set(1)%tersoff%bigR + &
1460 nonbonded%pot(start + isec)%pot%set(1)%tersoff%bigD)**2
1461 nonbonded%pot(start + isec)%pot%set(1)%tersoff%rcutsq = rcutsq
1462 nonbonded%pot(start + isec)%pot%rcutsq = rcutsq
1463
1464 ! In case it is defined override the standard specification of RCUT
1465 CALL section_vals_val_get(tersoff_section, "RCUT", i_rep_section=isec, n_rep_val=n_rep)
1466 IF (n_rep == 1) THEN
1467 CALL section_vals_val_get(tersoff_section, "RCUT", i_rep_section=isec, r_val=rcut)
1468 nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
1469 END IF
1470 END DO
1471 END SUBROUTINE read_tersoff_section
1472
1473! **************************************************************************************************
1474!> \brief Reads the gal19 section
1475!> \param nonbonded ...
1476!> \param section ...
1477!> \param start ...
1478!> \param gal_section ...
1479!> \author Clabaut Paul
1480! **************************************************************************************************
1481 SUBROUTINE read_gal_section(nonbonded, section, start, gal_section)
1482 TYPE(pair_potential_p_type), POINTER :: nonbonded
1483 TYPE(section_vals_type), POINTER :: section
1484 INTEGER, INTENT(IN) :: start
1485 TYPE(section_vals_type), POINTER :: gal_section
1486
1487 CHARACTER(LEN=default_string_length), &
1488 DIMENSION(:), POINTER :: atm_names
1489 INTEGER :: iatom, isec, n_items, n_rep, nval
1490 LOGICAL :: is_ok
1491 REAL(kind=dp) :: rcut, rval
1492 REAL(kind=dp), DIMENSION(:), POINTER :: rvalues
1493 TYPE(cp_sll_val_type), POINTER :: list
1494 TYPE(section_vals_type), POINTER :: subsection
1495 TYPE(val_type), POINTER :: val
1496
1497 CALL section_vals_get(section, n_repetition=n_items)
1498 DO isec = 1, n_items
1499 CALL cite_reference(clabaut2020)
1500 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
1501
1502 nonbonded%pot(start + isec)%pot%type = gal_type
1503 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
1504 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
1505 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
1506 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
1507
1508 CALL section_vals_val_get(section, "METALS", i_rep_section=isec, c_vals=atm_names)
1509 nonbonded%pot(start + isec)%pot%set(1)%gal%met1 = atm_names(1)
1510 nonbonded%pot(start + isec)%pot%set(1)%gal%met2 = atm_names(2)
1511
1512 CALL section_vals_val_get(gal_section, "epsilon", i_rep_section=isec, &
1513 r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%epsilon)
1514 CALL section_vals_val_get(gal_section, "bxy", i_rep_section=isec, &
1515 r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%bxy)
1516 CALL section_vals_val_get(gal_section, "bz", i_rep_section=isec, &
1517 r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%bz)
1518
1519 CALL section_vals_val_get(gal_section, "r", i_rep_section=isec, r_vals=rvalues)
1520 nonbonded%pot(start + isec)%pot%set(1)%gal%r1 = rvalues(1)
1521 nonbonded%pot(start + isec)%pot%set(1)%gal%r2 = rvalues(2)
1522
1523 CALL section_vals_val_get(gal_section, "a1", i_rep_section=isec, &
1524 r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a1)
1525 CALL section_vals_val_get(gal_section, "a2", i_rep_section=isec, &
1526 r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a2)
1527 CALL section_vals_val_get(gal_section, "a3", i_rep_section=isec, &
1528 r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a3)
1529 CALL section_vals_val_get(gal_section, "a4", i_rep_section=isec, &
1530 r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a4)
1531 CALL section_vals_val_get(gal_section, "A", i_rep_section=isec, &
1532 r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%a)
1533 CALL section_vals_val_get(gal_section, "B", i_rep_section=isec, &
1534 r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%b)
1535 CALL section_vals_val_get(gal_section, "C", i_rep_section=isec, &
1536 r_val=nonbonded%pot(start + isec)%pot%set(1)%gal%c)
1537 NULLIFY (list)
1538 subsection => section_vals_get_subs_vals(section, "GCN", i_rep_section=isec)
1539 CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=nval)
1540 ALLOCATE (nonbonded%pot(start + isec)%pot%set(1)%gal%gcn(nval))
1541 CALL section_vals_list_get(subsection, "_DEFAULT_KEYWORD_", list=list)
1542 DO iatom = 1, nval
1543 ! we use only the first default_string_length characters of each line
1544 is_ok = cp_sll_val_next(list, val)
1545 CALL val_get(val, r_val=rval)
1546 ! assign values
1547 nonbonded%pot(start + isec)%pot%set(1)%gal%gcn(iatom) = rval
1548 END DO
1549
1550 CALL section_vals_val_get(gal_section, "Fit_express", i_rep_section=isec, &
1551 l_val=nonbonded%pot(start + isec)%pot%set(1)%gal%express)
1552
1553 ! ! In case it is defined override the standard specification of RCUT
1554 CALL section_vals_val_get(gal_section, "RCUT", i_rep_section=isec, n_rep_val=n_rep)
1555 IF (n_rep == 1) THEN
1556 CALL section_vals_val_get(gal_section, "RCUT", i_rep_section=isec, r_val=rcut)
1557 nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
1558 nonbonded%pot(start + isec)%pot%set(1)%gal%rcutsq = rcut**2
1559 END IF
1560 END DO
1561 END SUBROUTINE read_gal_section
1562
1563! **************************************************************************************************
1564!> \brief Reads the gal21 section
1565!> \param nonbonded ...
1566!> \param section ...
1567!> \param start ...
1568!> \param gal21_section ...
1569!> \author Clabaut Paul
1570! **************************************************************************************************
1571 SUBROUTINE read_gal21_section(nonbonded, section, start, gal21_section)
1572 TYPE(pair_potential_p_type), POINTER :: nonbonded
1573 TYPE(section_vals_type), POINTER :: section
1574 INTEGER, INTENT(IN) :: start
1575 TYPE(section_vals_type), POINTER :: gal21_section
1576
1577 CHARACTER(LEN=default_string_length), &
1578 DIMENSION(:), POINTER :: atm_names
1579 INTEGER :: iatom, isec, n_items, n_rep, nval
1580 LOGICAL :: is_ok
1581 REAL(kind=dp) :: rcut, rval
1582 REAL(kind=dp), DIMENSION(:), POINTER :: rvalues
1583 TYPE(cp_sll_val_type), POINTER :: list
1584 TYPE(section_vals_type), POINTER :: subsection
1585 TYPE(val_type), POINTER :: val
1586
1587 CALL section_vals_get(section, n_repetition=n_items)
1588 DO isec = 1, n_items
1589 CALL cite_reference(clabaut2021)
1590 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
1591
1592 nonbonded%pot(start + isec)%pot%type = gal21_type
1593 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
1594 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
1595 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
1596 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
1597
1598 CALL section_vals_val_get(section, "METALS", i_rep_section=isec, c_vals=atm_names)
1599 nonbonded%pot(start + isec)%pot%set(1)%gal21%met1 = atm_names(1)
1600 nonbonded%pot(start + isec)%pot%set(1)%gal21%met2 = atm_names(2)
1601
1602 CALL section_vals_val_get(gal21_section, "epsilon", i_rep_section=isec, r_vals=rvalues)
1603 nonbonded%pot(start + isec)%pot%set(1)%gal21%epsilon1 = rvalues(1)
1604 nonbonded%pot(start + isec)%pot%set(1)%gal21%epsilon2 = rvalues(2)
1605 nonbonded%pot(start + isec)%pot%set(1)%gal21%epsilon3 = rvalues(3)
1606
1607 CALL section_vals_val_get(gal21_section, "bxy", i_rep_section=isec, r_vals=rvalues)
1608 nonbonded%pot(start + isec)%pot%set(1)%gal21%bxy1 = rvalues(1)
1609 nonbonded%pot(start + isec)%pot%set(1)%gal21%bxy2 = rvalues(2)
1610
1611 CALL section_vals_val_get(gal21_section, "bz", i_rep_section=isec, r_vals=rvalues)
1612 nonbonded%pot(start + isec)%pot%set(1)%gal21%bz1 = rvalues(1)
1613 nonbonded%pot(start + isec)%pot%set(1)%gal21%bz2 = rvalues(2)
1614
1615 CALL section_vals_val_get(gal21_section, "r", i_rep_section=isec, r_vals=rvalues)
1616 nonbonded%pot(start + isec)%pot%set(1)%gal21%r1 = rvalues(1)
1617 nonbonded%pot(start + isec)%pot%set(1)%gal21%r2 = rvalues(2)
1618
1619 CALL section_vals_val_get(gal21_section, "a1", i_rep_section=isec, r_vals=rvalues)
1620 nonbonded%pot(start + isec)%pot%set(1)%gal21%a11 = rvalues(1)
1621 nonbonded%pot(start + isec)%pot%set(1)%gal21%a12 = rvalues(2)
1622 nonbonded%pot(start + isec)%pot%set(1)%gal21%a13 = rvalues(3)
1623
1624 CALL section_vals_val_get(gal21_section, "a2", i_rep_section=isec, r_vals=rvalues)
1625 nonbonded%pot(start + isec)%pot%set(1)%gal21%a21 = rvalues(1)
1626 nonbonded%pot(start + isec)%pot%set(1)%gal21%a22 = rvalues(2)
1627 nonbonded%pot(start + isec)%pot%set(1)%gal21%a23 = rvalues(3)
1628
1629 CALL section_vals_val_get(gal21_section, "a3", i_rep_section=isec, r_vals=rvalues)
1630 nonbonded%pot(start + isec)%pot%set(1)%gal21%a31 = rvalues(1)
1631 nonbonded%pot(start + isec)%pot%set(1)%gal21%a32 = rvalues(2)
1632 nonbonded%pot(start + isec)%pot%set(1)%gal21%a33 = rvalues(3)
1633
1634 CALL section_vals_val_get(gal21_section, "a4", i_rep_section=isec, r_vals=rvalues)
1635 nonbonded%pot(start + isec)%pot%set(1)%gal21%a41 = rvalues(1)
1636 nonbonded%pot(start + isec)%pot%set(1)%gal21%a42 = rvalues(2)
1637 nonbonded%pot(start + isec)%pot%set(1)%gal21%a43 = rvalues(3)
1638
1639 CALL section_vals_val_get(gal21_section, "A", i_rep_section=isec, r_vals=rvalues)
1640 nonbonded%pot(start + isec)%pot%set(1)%gal21%AO1 = rvalues(1)
1641 nonbonded%pot(start + isec)%pot%set(1)%gal21%AO2 = rvalues(2)
1642
1643 CALL section_vals_val_get(gal21_section, "B", i_rep_section=isec, r_vals=rvalues)
1644 nonbonded%pot(start + isec)%pot%set(1)%gal21%BO1 = rvalues(1)
1645 nonbonded%pot(start + isec)%pot%set(1)%gal21%BO2 = rvalues(2)
1646
1647 CALL section_vals_val_get(gal21_section, "C", i_rep_section=isec, &
1648 r_val=nonbonded%pot(start + isec)%pot%set(1)%gal21%c)
1649
1650 CALL section_vals_val_get(gal21_section, "AH", i_rep_section=isec, r_vals=rvalues)
1651 nonbonded%pot(start + isec)%pot%set(1)%gal21%AH1 = rvalues(1)
1652 nonbonded%pot(start + isec)%pot%set(1)%gal21%AH2 = rvalues(2)
1653
1654 CALL section_vals_val_get(gal21_section, "BH", i_rep_section=isec, r_vals=rvalues)
1655 nonbonded%pot(start + isec)%pot%set(1)%gal21%BH1 = rvalues(1)
1656 nonbonded%pot(start + isec)%pot%set(1)%gal21%BH2 = rvalues(2)
1657
1658 NULLIFY (list)
1659 subsection => section_vals_get_subs_vals(section, "GCN", i_rep_section=isec)
1660 CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=nval)
1661 ALLOCATE (nonbonded%pot(start + isec)%pot%set(1)%gal21%gcn(nval))
1662 CALL section_vals_list_get(subsection, "_DEFAULT_KEYWORD_", list=list)
1663 DO iatom = 1, nval
1664 ! we use only the first default_string_length characters of each line
1665 is_ok = cp_sll_val_next(list, val)
1666 CALL val_get(val, r_val=rval)
1667 ! assign values
1668 nonbonded%pot(start + isec)%pot%set(1)%gal21%gcn(iatom) = rval
1669 END DO
1670
1671 CALL section_vals_val_get(gal21_section, "Fit_express", i_rep_section=isec, &
1672 l_val=nonbonded%pot(start + isec)%pot%set(1)%gal21%express)
1673
1674 ! ! In case it is defined override the standard specification of RCUT
1675 CALL section_vals_val_get(gal21_section, "RCUT", i_rep_section=isec, n_rep_val=n_rep)
1676 IF (n_rep == 1) THEN
1677 CALL section_vals_val_get(gal21_section, "RCUT", i_rep_section=isec, r_val=rcut)
1678 nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
1679 nonbonded%pot(start + isec)%pot%set(1)%gal21%rcutsq = rcut**2
1680 END IF
1681 END DO
1682 END SUBROUTINE read_gal21_section
1683
1684! **************************************************************************************************
1685!> \brief Reads the siepmann section
1686!> \param nonbonded ...
1687!> \param section ...
1688!> \param start ...
1689!> \param siepmann_section ...
1690!> \author Dorothea Golze
1691! **************************************************************************************************
1692 SUBROUTINE read_siepmann_section(nonbonded, section, start, siepmann_section)
1693 TYPE(pair_potential_p_type), POINTER :: nonbonded
1694 TYPE(section_vals_type), POINTER :: section
1695 INTEGER, INTENT(IN) :: start
1696 TYPE(section_vals_type), POINTER :: siepmann_section
1697
1698 CHARACTER(LEN=default_string_length), &
1699 DIMENSION(:), POINTER :: atm_names
1700 INTEGER :: isec, n_items, n_rep
1701 REAL(kind=dp) :: rcut
1702
1703 CALL section_vals_get(section, n_repetition=n_items)
1704 DO isec = 1, n_items
1705 CALL cite_reference(siepmann1995)
1706 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
1707
1708 nonbonded%pot(start + isec)%pot%type = siepmann_type
1709 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
1710 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
1711 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
1712 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
1713
1714 CALL section_vals_val_get(siepmann_section, "B", i_rep_section=isec, &
1715 r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%B)
1716 CALL section_vals_val_get(siepmann_section, "D", i_rep_section=isec, &
1717 r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%D)
1718 CALL section_vals_val_get(siepmann_section, "E", i_rep_section=isec, &
1719 r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%E)
1720 CALL section_vals_val_get(siepmann_section, "F", i_rep_section=isec, &
1721 r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%F)
1722 CALL section_vals_val_get(siepmann_section, "beta", i_rep_section=isec, &
1723 r_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%beta)
1724 CALL section_vals_val_get(siepmann_section, "ALLOW_OH_FORMATION", i_rep_section=isec, &
1725 l_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%allow_oh_formation)
1726 CALL section_vals_val_get(siepmann_section, "ALLOW_H3O_FORMATION", i_rep_section=isec, &
1727 l_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%allow_h3o_formation)
1728 CALL section_vals_val_get(siepmann_section, "ALLOW_O_FORMATION", i_rep_section=isec, &
1729 l_val=nonbonded%pot(start + isec)%pot%set(1)%siepmann%allow_o_formation)
1730
1731 ! ! In case it is defined override the standard specification of RCUT
1732 CALL section_vals_val_get(siepmann_section, "RCUT", i_rep_section=isec, n_rep_val=n_rep)
1733 IF (n_rep == 1) THEN
1734 CALL section_vals_val_get(siepmann_section, "RCUT", i_rep_section=isec, r_val=rcut)
1735 nonbonded%pot(start + isec)%pot%rcutsq = rcut**2
1736 nonbonded%pot(start + isec)%pot%set(1)%siepmann%rcutsq = rcut**2
1737 END IF
1738 END DO
1739 END SUBROUTINE read_siepmann_section
1740
1741! **************************************************************************************************
1742!> \brief Reads the Buckingham plus Morse potential section
1743!> \param nonbonded ...
1744!> \param section ...
1745!> \param start ...
1746!> \author MI
1747! **************************************************************************************************
1748 SUBROUTINE read_bm_section(nonbonded, section, start)
1749 TYPE(pair_potential_p_type), POINTER :: nonbonded
1750 TYPE(section_vals_type), POINTER :: section
1751 INTEGER, INTENT(IN) :: start
1752
1753 CHARACTER(LEN=default_string_length), &
1754 DIMENSION(:), POINTER :: atm_names
1755 INTEGER :: isec, n_items, n_rep
1756 REAL(kind=dp) :: a1, a2, b1, b2, beta, c, d, f0, r0, rcut
1757
1758 CALL section_vals_get(section, n_repetition=n_items)
1759 DO isec = 1, n_items
1760 CALL cite_reference(yamada2000)
1761 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
1762 CALL section_vals_val_get(section, "F0", i_rep_section=isec, r_val=f0)
1763 CALL section_vals_val_get(section, "A1", i_rep_section=isec, r_val=a1)
1764 CALL section_vals_val_get(section, "A2", i_rep_section=isec, r_val=a2)
1765 CALL section_vals_val_get(section, "B1", i_rep_section=isec, r_val=b1)
1766 CALL section_vals_val_get(section, "B2", i_rep_section=isec, r_val=b2)
1767 CALL section_vals_val_get(section, "C", i_rep_section=isec, r_val=c)
1768 CALL section_vals_val_get(section, "D", i_rep_section=isec, r_val=d)
1769 CALL section_vals_val_get(section, "R0", i_rep_section=isec, r_val=r0)
1770 CALL section_vals_val_get(section, "Beta", i_rep_section=isec, r_val=beta)
1771 CALL section_vals_val_get(section, "RCUT", i_rep_section=isec, r_val=rcut)
1772
1773 nonbonded%pot(start + isec)%pot%type = bm_type
1774 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
1775 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
1776 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
1777 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
1778 nonbonded%pot(start + isec)%pot%set(1)%buckmo%f0 = f0
1779 nonbonded%pot(start + isec)%pot%set(1)%buckmo%a1 = a1
1780 nonbonded%pot(start + isec)%pot%set(1)%buckmo%a2 = a2
1781 nonbonded%pot(start + isec)%pot%set(1)%buckmo%b1 = b1
1782 nonbonded%pot(start + isec)%pot%set(1)%buckmo%b2 = b2
1783 nonbonded%pot(start + isec)%pot%set(1)%buckmo%c = c
1784 nonbonded%pot(start + isec)%pot%set(1)%buckmo%d = d
1785 nonbonded%pot(start + isec)%pot%set(1)%buckmo%r0 = r0
1786 nonbonded%pot(start + isec)%pot%set(1)%buckmo%beta = beta
1787 nonbonded%pot(start + isec)%pot%rcutsq = rcut*rcut
1788 !
1789 CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, n_rep_val=n_rep)
1790 IF (n_rep == 1) CALL section_vals_val_get(section, "RMIN", i_rep_section=isec, &
1791 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmin)
1792 CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, n_rep_val=n_rep)
1793 IF (n_rep == 1) CALL section_vals_val_get(section, "RMAX", i_rep_section=isec, &
1794 r_val=nonbonded%pot(start + isec)%pot%set(1)%rmax)
1795 END DO
1796 END SUBROUTINE read_bm_section
1797
1798! **************************************************************************************************
1799!> \brief Reads the TABPOT section
1800!> \param nonbonded ...
1801!> \param section ...
1802!> \param start ...
1803!> \param para_env ...
1804!> \param mm_section ...
1805!> \author Alex Mironenko, Da Teng
1806! **************************************************************************************************
1807 SUBROUTINE read_tabpot_section(nonbonded, section, start, para_env, mm_section)
1808 TYPE(pair_potential_p_type), POINTER :: nonbonded
1809 TYPE(section_vals_type), POINTER :: section
1810 INTEGER, INTENT(IN) :: start
1811 TYPE(mp_para_env_type), POINTER :: para_env
1812 TYPE(section_vals_type), POINTER :: mm_section
1813
1814 CHARACTER(LEN=default_string_length), &
1815 DIMENSION(:), POINTER :: atm_names
1816 INTEGER :: isec, n_items
1817
1818 CALL section_vals_get(section, n_repetition=n_items)
1819 DO isec = 1, n_items
1820 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
1821 nonbonded%pot(start + isec)%pot%type = tab_type
1822 nonbonded%pot(start + isec)%pot%at1 = atm_names(1)
1823 nonbonded%pot(start + isec)%pot%at2 = atm_names(2)
1824 CALL uppercase(nonbonded%pot(start + isec)%pot%at1)
1825 CALL uppercase(nonbonded%pot(start + isec)%pot%at2)
1826 CALL section_vals_val_get(section, "PARM_FILE_NAME", i_rep_section=isec, &
1827 c_val=nonbonded%pot(start + isec)%pot%set(1)%tab%tabpot_file_name)
1828 CALL read_tabpot_data(nonbonded%pot(start + isec)%pot%set(1)%tab, para_env, mm_section)
1829 nonbonded%pot(start + isec)%pot%set(1)%tab%index = isec
1830 END DO
1831 END SUBROUTINE read_tabpot_section
1832
1833! **************************************************************************************************
1834!> \brief Reads the CHARGE section
1835!> \param charge_atm ...
1836!> \param charge ...
1837!> \param section ...
1838!> \param start ...
1839!> \author teo
1840! **************************************************************************************************
1841 SUBROUTINE read_chrg_section(charge_atm, charge, section, start)
1842 CHARACTER(LEN=default_string_length), &
1843 DIMENSION(:), POINTER :: charge_atm
1844 REAL(kind=dp), DIMENSION(:), POINTER :: charge
1845 TYPE(section_vals_type), POINTER :: section
1846 INTEGER, INTENT(IN) :: start
1847
1848 CHARACTER(LEN=default_string_length) :: atm_name
1849 INTEGER :: isec, n_items
1850
1851 CALL section_vals_get(section, n_repetition=n_items)
1852 DO isec = 1, n_items
1853 CALL section_vals_val_get(section, "ATOM", i_rep_section=isec, c_val=atm_name)
1854 charge_atm(start + isec) = atm_name
1855 CALL uppercase(charge_atm(start + isec))
1856 CALL section_vals_val_get(section, "CHARGE", i_rep_section=isec, r_val=charge(start + isec))
1857 END DO
1858 END SUBROUTINE read_chrg_section
1859
1860! **************************************************************************************************
1861!> \brief Reads the POLARIZABILITY section
1862!> \param apol_atm ...
1863!> \param apol ...
1864!> \param damping_list ...
1865!> \param section ...
1866!> \param start ...
1867!> \author Marcel Baer
1868! **************************************************************************************************
1869 SUBROUTINE read_apol_section(apol_atm, apol, damping_list, section, &
1870 start)
1871 CHARACTER(LEN=default_string_length), &
1872 DIMENSION(:), POINTER :: apol_atm
1873 REAL(kind=dp), DIMENSION(:), POINTER :: apol
1874 TYPE(damping_info_type), DIMENSION(:), POINTER :: damping_list
1875 TYPE(section_vals_type), POINTER :: section
1876 INTEGER, INTENT(IN) :: start
1877
1878 CHARACTER(LEN=default_string_length) :: atm_name
1879 INTEGER :: isec, isec_damp, n_damp, n_items, &
1880 start_damp, tmp_damp
1881 TYPE(section_vals_type), POINTER :: tmp_section
1882
1883 CALL section_vals_get(section, n_repetition=n_items)
1884 NULLIFY (tmp_section)
1885 n_damp = 0
1886! *** Counts number of DIPOLE%DAMPING sections ****
1887 DO isec = 1, n_items
1888 tmp_section => section_vals_get_subs_vals(section, "DAMPING", &
1889 i_rep_section=isec)
1890 CALL section_vals_get(tmp_section, n_repetition=tmp_damp)
1891 n_damp = n_damp + tmp_damp
1892
1893 END DO
1894
1895 IF (n_damp > 0) THEN
1896 ALLOCATE (damping_list(1:n_damp))
1897 END IF
1898
1899! *** Reads DIPOLE sections *****
1900 start_damp = 0
1901 DO isec = 1, n_items
1902 CALL section_vals_val_get(section, "ATOM", i_rep_section=isec, c_val=atm_name)
1903 apol_atm(start + isec) = atm_name
1904 CALL uppercase(apol_atm(start + isec))
1905 CALL section_vals_val_get(section, "APOL", i_rep_section=isec, r_val=apol(start + isec))
1906
1907 tmp_section => section_vals_get_subs_vals(section, "DAMPING", &
1908 i_rep_section=isec)
1909 CALL section_vals_get(tmp_section, n_repetition=tmp_damp)
1910 DO isec_damp = 1, tmp_damp
1911 damping_list(start_damp + isec_damp)%atm_name1 = apol_atm(start + isec)
1912 CALL section_vals_val_get(tmp_section, "ATOM", i_rep_section=isec_damp, &
1913 c_val=atm_name)
1914 damping_list(start_damp + isec_damp)%atm_name2 = atm_name
1915 CALL uppercase(damping_list(start_damp + isec_damp)%atm_name2)
1916 CALL section_vals_val_get(tmp_section, "TYPE", i_rep_section=isec_damp, &
1917 c_val=atm_name)
1918 damping_list(start_damp + isec_damp)%dtype = atm_name
1919 CALL uppercase(damping_list(start_damp + isec_damp)%dtype)
1920
1921 CALL section_vals_val_get(tmp_section, "ORDER", i_rep_section=isec_damp, &
1922 i_val=damping_list(start_damp + isec_damp)%order)
1923 CALL section_vals_val_get(tmp_section, "BIJ", i_rep_section=isec_damp, &
1924 r_val=damping_list(start_damp + isec_damp)%bij)
1925 CALL section_vals_val_get(tmp_section, "CIJ", i_rep_section=isec_damp, &
1926 r_val=damping_list(start_damp + isec_damp)%cij)
1927 END DO
1928 start_damp = start_damp + tmp_damp
1929
1930 END DO
1931
1932 END SUBROUTINE read_apol_section
1933
1934! **************************************************************************************************
1935!> \brief Reads the QUADRUPOLE POLARIZABILITY section
1936!> \param cpol_atm ...
1937!> \param cpol ...
1938!> \param section ...
1939!> \param start ...
1940!> \author Marcel Baer
1941! **************************************************************************************************
1942 SUBROUTINE read_cpol_section(cpol_atm, cpol, section, start)
1943 CHARACTER(LEN=default_string_length), &
1944 DIMENSION(:), POINTER :: cpol_atm
1945 REAL(kind=dp), DIMENSION(:), POINTER :: cpol
1946 TYPE(section_vals_type), POINTER :: section
1947 INTEGER, INTENT(IN) :: start
1948
1949 CHARACTER(LEN=default_string_length) :: atm_name
1950 INTEGER :: isec, n_items
1951
1952 CALL section_vals_get(section, n_repetition=n_items)
1953 DO isec = 1, n_items
1954 CALL section_vals_val_get(section, "ATOM", i_rep_section=isec, c_val=atm_name)
1955 cpol_atm(start + isec) = atm_name
1956 CALL uppercase(cpol_atm(start + isec))
1957 CALL section_vals_val_get(section, "CPOL", i_rep_section=isec, r_val=cpol(start + isec))
1958 END DO
1959 END SUBROUTINE read_cpol_section
1960
1961! **************************************************************************************************
1962!> \brief Reads the SHELL section
1963!> \param shell_list ...
1964!> \param section ...
1965!> \param start ...
1966!> \author Marcella Iannuzzi
1967! **************************************************************************************************
1968 SUBROUTINE read_shell_section(shell_list, section, start)
1969
1970 TYPE(shell_p_type), DIMENSION(:), POINTER :: shell_list
1971 TYPE(section_vals_type), POINTER :: section
1972 INTEGER, INTENT(IN) :: start
1973
1974 CHARACTER(LEN=default_string_length) :: atm_name
1975 INTEGER :: i_rep, n_rep
1976 REAL(dp) :: ccharge, cutoff, k, maxdist, mfrac, &
1977 scharge
1978
1979 CALL section_vals_get(section, n_repetition=n_rep)
1980
1981 DO i_rep = 1, n_rep
1982 CALL section_vals_val_get(section, "_SECTION_PARAMETERS_", &
1983 c_val=atm_name, i_rep_section=i_rep)
1984 CALL uppercase(atm_name)
1985 shell_list(start + i_rep)%atm_name = atm_name
1986 CALL section_vals_val_get(section, "CORE_CHARGE", i_rep_section=i_rep, r_val=ccharge)
1987 shell_list(start + i_rep)%shell%charge_core = ccharge
1988 CALL section_vals_val_get(section, "SHELL_CHARGE", i_rep_section=i_rep, r_val=scharge)
1989 shell_list(start + i_rep)%shell%charge_shell = scharge
1990 CALL section_vals_val_get(section, "MASS_FRACTION", i_rep_section=i_rep, r_val=mfrac)
1991 shell_list(start + i_rep)%shell%massfrac = mfrac
1992 CALL section_vals_val_get(section, "K2_SPRING", i_rep_section=i_rep, r_val=k)
1993 IF (k < 0.0_dp) THEN
1994 CALL cp_abort(__location__, &
1995 "An invalid value was specified for the force constant k2 of the core-shell "// &
1996 "spring potential")
1997 END IF
1998 shell_list(start + i_rep)%shell%k2_spring = k
1999 CALL section_vals_val_get(section, "K4_SPRING", i_rep_section=i_rep, r_val=k)
2000 IF (k < 0.0_dp) THEN
2001 CALL cp_abort(__location__, &
2002 "An invalid value was specified for the force constant k4 of the core-shell "// &
2003 "spring potential")
2004 END IF
2005 shell_list(start + i_rep)%shell%k4_spring = k
2006 CALL section_vals_val_get(section, "MAX_DISTANCE", i_rep_section=i_rep, r_val=maxdist)
2007 shell_list(start + i_rep)%shell%max_dist = maxdist
2008 CALL section_vals_val_get(section, "SHELL_CUTOFF", i_rep_section=i_rep, r_val=cutoff)
2009 shell_list(start + i_rep)%shell%shell_cutoff = cutoff
2010 END DO
2011
2012 END SUBROUTINE read_shell_section
2013
2014! **************************************************************************************************
2015!> \brief Reads the BONDS section
2016!> \param bond_kind ...
2017!> \param bond_a ...
2018!> \param bond_b ...
2019!> \param bond_k ...
2020!> \param bond_r0 ...
2021!> \param bond_cs ...
2022!> \param section ...
2023!> \param start ...
2024!> \author teo
2025! **************************************************************************************************
2026 SUBROUTINE read_bonds_section(bond_kind, bond_a, bond_b, bond_k, bond_r0, bond_cs, section, start)
2027 INTEGER, DIMENSION(:), POINTER :: bond_kind
2028 CHARACTER(LEN=default_string_length), &
2029 DIMENSION(:), POINTER :: bond_a, bond_b
2030 REAL(kind=dp), DIMENSION(:, :), POINTER :: bond_k
2031 REAL(kind=dp), DIMENSION(:), POINTER :: bond_r0, bond_cs
2032 TYPE(section_vals_type), POINTER :: section
2033 INTEGER, INTENT(IN) :: start
2034
2035 CHARACTER(LEN=default_string_length), &
2036 DIMENSION(:), POINTER :: atm_names
2037 INTEGER :: isec, k, n_items
2038 REAL(kind=dp), DIMENSION(:), POINTER :: kvals
2039
2040 NULLIFY (kvals, atm_names)
2041 CALL section_vals_get(section, n_repetition=n_items)
2042 DO isec = 1, n_items
2043 CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=bond_kind(start + isec))
2044 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
2045 bond_a(start + isec) = atm_names(1)
2046 bond_b(start + isec) = atm_names(2)
2047 CALL uppercase(bond_a(start + isec))
2048 CALL uppercase(bond_b(start + isec))
2049 CALL section_vals_val_get(section, "K", i_rep_section=isec, r_vals=kvals)
2050 cpassert(SIZE(kvals) <= 3)
2051 bond_k(:, start + isec) = 0.0_dp
2052 DO k = 1, SIZE(kvals)
2053 bond_k(k, start + isec) = kvals(k)
2054 END DO
2055 CALL section_vals_val_get(section, "R0", i_rep_section=isec, r_val=bond_r0(start + isec))
2056 CALL section_vals_val_get(section, "CS", i_rep_section=isec, r_val=bond_cs(start + isec))
2057 END DO
2058 END SUBROUTINE read_bonds_section
2059
2060! **************************************************************************************************
2061!> \brief Reads the BENDS section
2062!> \param bend_kind ...
2063!> \param bend_a ...
2064!> \param bend_b ...
2065!> \param bend_c ...
2066!> \param bend_k ...
2067!> \param bend_theta0 ...
2068!> \param bend_cb ...
2069!> \param bend_r012 ...
2070!> \param bend_r032 ...
2071!> \param bend_kbs12 ...
2072!> \param bend_kbs32 ...
2073!> \param bend_kss ...
2074!> \param bend_legendre ...
2075!> \param section ...
2076!> \param start ...
2077!> \author teo
2078! **************************************************************************************************
2079 SUBROUTINE read_bends_section(bend_kind, bend_a, bend_b, bend_c, bend_k, bend_theta0, bend_cb, &
2080 bend_r012, bend_r032, bend_kbs12, bend_kbs32, bend_kss, bend_legendre, &
2081 section, start)
2082 INTEGER, DIMENSION(:), POINTER :: bend_kind
2083 CHARACTER(LEN=default_string_length), &
2084 DIMENSION(:), POINTER :: bend_a, bend_b, bend_c
2085 REAL(kind=dp), DIMENSION(:), POINTER :: bend_k, bend_theta0, bend_cb, bend_r012, &
2086 bend_r032, bend_kbs12, bend_kbs32, &
2087 bend_kss
2088 TYPE(legendre_data_type), DIMENSION(:), POINTER :: bend_legendre
2089 TYPE(section_vals_type), POINTER :: section
2090 INTEGER, INTENT(IN) :: start
2091
2092 CHARACTER(LEN=default_string_length), &
2093 DIMENSION(:), POINTER :: atm_names
2094 INTEGER :: isec, k, n_items, n_rep
2095 REAL(kind=dp), DIMENSION(:), POINTER :: kvals, r_values
2096
2097 NULLIFY (kvals, atm_names)
2098 CALL section_vals_get(section, n_repetition=n_items)
2099 bend_legendre%order = 0
2100 DO isec = 1, n_items
2101 CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=bend_kind(start + isec))
2102 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
2103 bend_a(start + isec) = atm_names(1)
2104 bend_b(start + isec) = atm_names(2)
2105 bend_c(start + isec) = atm_names(3)
2106 CALL uppercase(bend_a(start + isec))
2107 CALL uppercase(bend_b(start + isec))
2108 CALL uppercase(bend_c(start + isec))
2109 CALL section_vals_val_get(section, "K", i_rep_section=isec, r_vals=kvals)
2110 cpassert(SIZE(kvals) == 1)
2111 bend_k(start + isec) = kvals(1)
2112 CALL section_vals_val_get(section, "THETA0", i_rep_section=isec, r_val=bend_theta0(start + isec))
2113 CALL section_vals_val_get(section, "CB", i_rep_section=isec, r_val=bend_cb(start + isec))
2114 CALL section_vals_val_get(section, "R012", i_rep_section=isec, r_val=bend_r012(start + isec))
2115 CALL section_vals_val_get(section, "R032", i_rep_section=isec, r_val=bend_r032(start + isec))
2116 CALL section_vals_val_get(section, "KBS12", i_rep_section=isec, r_val=bend_kbs12(start + isec))
2117 CALL section_vals_val_get(section, "KBS32", i_rep_section=isec, r_val=bend_kbs32(start + isec))
2118 CALL section_vals_val_get(section, "KSS", i_rep_section=isec, r_val=bend_kss(start + isec))
2119 ! get legendre based data
2120 CALL section_vals_val_get(section, "LEGENDRE", i_rep_section=isec, n_rep_val=n_rep)
2121 DO k = 1, n_rep
2122 CALL section_vals_val_get(section, "LEGENDRE", i_rep_val=k, r_vals=r_values, i_rep_section=isec)
2123 bend_legendre(start + isec)%order = SIZE(r_values)
2124 IF (ASSOCIATED(bend_legendre(start + isec)%coeffs)) THEN
2125 DEALLOCATE (bend_legendre(start + isec)%coeffs)
2126 END IF
2127 ALLOCATE (bend_legendre(start + isec)%coeffs(bend_legendre(start + isec)%order))
2128 bend_legendre(start + isec)%coeffs = r_values
2129 END DO
2130 END DO
2131 END SUBROUTINE read_bends_section
2132
2133! **************************************************************************************************
2134!> \brief ...
2135!> \param ub_kind ...
2136!> \param ub_a ...
2137!> \param ub_b ...
2138!> \param ub_c ...
2139!> \param ub_k ...
2140!> \param ub_r0 ...
2141!> \param section ...
2142!> \param start ...
2143! **************************************************************************************************
2144 SUBROUTINE read_ubs_section(ub_kind, ub_a, ub_b, ub_c, ub_k, ub_r0, section, start)
2145 INTEGER, DIMENSION(:), POINTER :: ub_kind
2146 CHARACTER(LEN=default_string_length), &
2147 DIMENSION(:), POINTER :: ub_a, ub_b, ub_c
2148 REAL(kind=dp), DIMENSION(:, :), POINTER :: ub_k
2149 REAL(kind=dp), DIMENSION(:), POINTER :: ub_r0
2150 TYPE(section_vals_type), POINTER :: section
2151 INTEGER, INTENT(IN) :: start
2152
2153 CHARACTER(LEN=default_string_length), &
2154 DIMENSION(:), POINTER :: atm_names
2155 INTEGER :: isec, k, n_items
2156 LOGICAL :: explicit
2157 REAL(kind=dp), DIMENSION(:), POINTER :: kvals
2158 TYPE(section_vals_type), POINTER :: subsection
2159
2160 NULLIFY (atm_names)
2161 CALL section_vals_get(section, n_repetition=n_items)
2162 DO isec = 1, n_items
2163 subsection => section_vals_get_subs_vals(section, "UB", i_rep_section=isec)
2164 CALL section_vals_get(subsection, explicit=explicit)
2165 IF (explicit) THEN
2166 CALL section_vals_val_get(subsection, "KIND", i_val=ub_kind(start + isec))
2167 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
2168 ub_a(start + isec) = atm_names(1)
2169 ub_b(start + isec) = atm_names(2)
2170 ub_c(start + isec) = atm_names(3)
2171 CALL uppercase(ub_a(start + isec))
2172 CALL uppercase(ub_b(start + isec))
2173 CALL uppercase(ub_c(start + isec))
2174 CALL section_vals_val_get(subsection, "K", r_vals=kvals)
2175 cpassert(SIZE(kvals) <= 3)
2176 ub_k(:, start + isec) = 0.0_dp
2177 DO k = 1, SIZE(kvals)
2178 ub_k(k, start + isec) = kvals(k)
2179 END DO
2180 CALL section_vals_val_get(subsection, "R0", r_val=ub_r0(start + isec))
2181 END IF
2182 END DO
2183 END SUBROUTINE read_ubs_section
2184
2185! **************************************************************************************************
2186!> \brief Reads the TORSIONS section
2187!> \param torsion_kind ...
2188!> \param torsion_a ...
2189!> \param torsion_b ...
2190!> \param torsion_c ...
2191!> \param torsion_d ...
2192!> \param torsion_k ...
2193!> \param torsion_phi0 ...
2194!> \param torsion_m ...
2195!> \param section ...
2196!> \param start ...
2197!> \author teo
2198! **************************************************************************************************
2199 SUBROUTINE read_torsions_section(torsion_kind, torsion_a, torsion_b, torsion_c, torsion_d, torsion_k, &
2200 torsion_phi0, torsion_m, section, start)
2201 INTEGER, DIMENSION(:), POINTER :: torsion_kind
2202 CHARACTER(LEN=default_string_length), &
2203 DIMENSION(:), POINTER :: torsion_a, torsion_b, torsion_c, &
2204 torsion_d
2205 REAL(kind=dp), DIMENSION(:), POINTER :: torsion_k, torsion_phi0
2206 INTEGER, DIMENSION(:), POINTER :: torsion_m
2207 TYPE(section_vals_type), POINTER :: section
2208 INTEGER, INTENT(IN) :: start
2209
2210 CHARACTER(LEN=default_string_length), &
2211 DIMENSION(:), POINTER :: atm_names
2212 INTEGER :: isec, n_items
2213
2214 NULLIFY (atm_names)
2215 CALL section_vals_get(section, n_repetition=n_items)
2216 DO isec = 1, n_items
2217 CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=torsion_kind(start + isec))
2218 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
2219 torsion_a(start + isec) = atm_names(1)
2220 torsion_b(start + isec) = atm_names(2)
2221 torsion_c(start + isec) = atm_names(3)
2222 torsion_d(start + isec) = atm_names(4)
2223 CALL uppercase(torsion_a(start + isec))
2224 CALL uppercase(torsion_b(start + isec))
2225 CALL uppercase(torsion_c(start + isec))
2226 CALL uppercase(torsion_d(start + isec))
2227 CALL section_vals_val_get(section, "K", i_rep_section=isec, r_val=torsion_k(start + isec))
2228 CALL section_vals_val_get(section, "PHI0", i_rep_section=isec, r_val=torsion_phi0(start + isec))
2229 CALL section_vals_val_get(section, "M", i_rep_section=isec, i_val=torsion_m(start + isec))
2230 ! Modify parameterisation for OPLS case
2231 IF (torsion_kind(start + isec) .EQ. do_ff_opls) THEN
2232 IF (torsion_phi0(start + isec) .NE. 0.0_dp) THEN
2233 CALL cp_warn(__location__, "PHI0 parameter was non-zero "// &
2234 "for an OPLS-type TORSION. It will be ignored.")
2235 END IF
2236 IF (modulo(torsion_m(start + isec), 2) .EQ. 0) THEN
2237 ! For even M, negate the cosine using a Pi phase factor
2238 torsion_phi0(start + isec) = pi
2239 END IF
2240 ! the K parameter appears as K/2 in the OPLS parameterisation
2241 torsion_k(start + isec) = torsion_k(start + isec)*0.5_dp
2242 END IF
2243 END DO
2244 END SUBROUTINE read_torsions_section
2245
2246! **************************************************************************************************
2247!> \brief Reads the IMPROPER section
2248!> \param impr_kind ...
2249!> \param impr_a ...
2250!> \param impr_b ...
2251!> \param impr_c ...
2252!> \param impr_d ...
2253!> \param impr_k ...
2254!> \param impr_phi0 ...
2255!> \param section ...
2256!> \param start ...
2257!> \author louis vanduyfhuys
2258! **************************************************************************************************
2259 SUBROUTINE read_improper_section(impr_kind, impr_a, impr_b, impr_c, impr_d, impr_k, &
2260 impr_phi0, section, start)
2261 INTEGER, DIMENSION(:), POINTER :: impr_kind
2262 CHARACTER(LEN=default_string_length), &
2263 DIMENSION(:), POINTER :: impr_a, impr_b, impr_c, impr_d
2264 REAL(kind=dp), DIMENSION(:), POINTER :: impr_k, impr_phi0
2265 TYPE(section_vals_type), POINTER :: section
2266 INTEGER, INTENT(IN) :: start
2267
2268 CHARACTER(LEN=default_string_length), &
2269 DIMENSION(:), POINTER :: atm_names
2270 INTEGER :: isec, n_items
2271
2272 NULLIFY (atm_names)
2273 CALL section_vals_get(section, n_repetition=n_items)
2274 DO isec = 1, n_items
2275 CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=impr_kind(start + isec))
2276 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
2277 impr_a(start + isec) = atm_names(1)
2278 impr_b(start + isec) = atm_names(2)
2279 impr_c(start + isec) = atm_names(3)
2280 impr_d(start + isec) = atm_names(4)
2281 CALL uppercase(impr_a(start + isec))
2282 CALL uppercase(impr_b(start + isec))
2283 CALL uppercase(impr_c(start + isec))
2284 CALL uppercase(impr_d(start + isec))
2285 CALL section_vals_val_get(section, "K", i_rep_section=isec, r_val=impr_k(start + isec))
2286 CALL section_vals_val_get(section, "PHI0", i_rep_section=isec, r_val=impr_phi0(start + isec))
2287 END DO
2288 END SUBROUTINE read_improper_section
2289
2290! **************************************************************************************************
2291!> \brief Reads the OPBEND section
2292!> \param opbend_kind ...
2293!> \param opbend_a ...
2294!> \param opbend_b ...
2295!> \param opbend_c ...
2296!> \param opbend_d ...
2297!> \param opbend_k ...
2298!> \param opbend_phi0 ...
2299!> \param section ...
2300!> \param start ...
2301!> \author louis vanduyfhuys
2302! **************************************************************************************************
2303 SUBROUTINE read_opbend_section(opbend_kind, opbend_a, opbend_b, opbend_c, opbend_d, opbend_k, &
2304 opbend_phi0, section, start)
2305 INTEGER, DIMENSION(:), POINTER :: opbend_kind
2306 CHARACTER(LEN=default_string_length), &
2307 DIMENSION(:), POINTER :: opbend_a, opbend_b, opbend_c, opbend_d
2308 REAL(kind=dp), DIMENSION(:), POINTER :: opbend_k, opbend_phi0
2309 TYPE(section_vals_type), POINTER :: section
2310 INTEGER, INTENT(IN) :: start
2311
2312 CHARACTER(LEN=default_string_length), &
2313 DIMENSION(:), POINTER :: atm_names
2314 INTEGER :: isec, n_items
2315
2316 NULLIFY (atm_names)
2317 CALL section_vals_get(section, n_repetition=n_items)
2318 DO isec = 1, n_items
2319 CALL section_vals_val_get(section, "KIND", i_rep_section=isec, i_val=opbend_kind(start + isec))
2320 CALL section_vals_val_get(section, "ATOMS", i_rep_section=isec, c_vals=atm_names)
2321 opbend_a(start + isec) = atm_names(1)
2322 opbend_b(start + isec) = atm_names(2)
2323 opbend_c(start + isec) = atm_names(3)
2324 opbend_d(start + isec) = atm_names(4)
2325 CALL uppercase(opbend_a(start + isec))
2326 CALL uppercase(opbend_b(start + isec))
2327 CALL uppercase(opbend_c(start + isec))
2328 CALL uppercase(opbend_d(start + isec))
2329 CALL section_vals_val_get(section, "K", i_rep_section=isec, r_val=opbend_k(start + isec))
2330 CALL section_vals_val_get(section, "PHI0", i_rep_section=isec, r_val=opbend_phi0(start + isec))
2331 END DO
2332 END SUBROUTINE read_opbend_section
2333
2334! **************************************************************************************************
2335!> \brief Reads the force_field input section
2336!> \param ff_type ...
2337!> \param para_env ...
2338!> \param mm_section ...
2339!> \par History
2340!> JGH (30.11.2001) : moved determination of setup variables to
2341!> molecule_input
2342!> \author CJM
2343! **************************************************************************************************
2344 SUBROUTINE read_force_field_section(ff_type, para_env, mm_section)
2345 TYPE(force_field_type), INTENT(INOUT) :: ff_type
2346 TYPE(mp_para_env_type), POINTER :: para_env
2347 TYPE(section_vals_type), POINTER :: mm_section
2348
2349 TYPE(section_vals_type), POINTER :: ff_section
2350
2351 NULLIFY (ff_section)
2352 ff_section => section_vals_get_subs_vals(mm_section, "FORCEFIELD")
2353 CALL read_force_field_section1(ff_section, mm_section, ff_type, para_env)
2354 END SUBROUTINE read_force_field_section
2355
2356! **************************************************************************************************
2357!> \brief reads EAM potential from library
2358!> \param eam ...
2359!> \param para_env ...
2360!> \param mm_section ...
2361! **************************************************************************************************
2362 SUBROUTINE read_eam_data(eam, para_env, mm_section)
2363 TYPE(eam_pot_type), POINTER :: eam
2364 TYPE(mp_para_env_type), POINTER :: para_env
2365 TYPE(section_vals_type), POINTER :: mm_section
2366
2367 CHARACTER(len=*), PARAMETER :: routinen = 'read_eam_data'
2368
2369 INTEGER :: handle, i, iw
2370 TYPE(cp_logger_type), POINTER :: logger
2371 TYPE(cp_parser_type) :: parser
2372
2373 CALL timeset(routinen, handle)
2374 NULLIFY (logger)
2375 logger => cp_get_default_logger()
2376 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
2377 extension=".mmLog")
2378 IF (iw > 0) WRITE (iw, *) "Reading EAM data from: ", trim(eam%eam_file_name)
2379 CALL parser_create(parser, trim(eam%eam_file_name), para_env=para_env)
2380
2381 CALL parser_get_next_line(parser, 1)
2382 IF (iw > 0) WRITE (iw, *) "Title: ", parser%input_line
2383
2384 CALL parser_get_next_line(parser, 2)
2385 READ (parser%input_line, *) eam%drar, eam%drhoar, eam%acutal, eam%npoints
2386 eam%drar = cp_unit_to_cp2k(eam%drar, "angstrom")
2387 eam%acutal = cp_unit_to_cp2k(eam%acutal, "angstrom")
2388 ! Relocating arrays with the right size
2389 CALL reallocate(eam%rho, 1, eam%npoints)
2390 CALL reallocate(eam%rhop, 1, eam%npoints)
2391 CALL reallocate(eam%rval, 1, eam%npoints)
2392 CALL reallocate(eam%rhoval, 1, eam%npoints)
2393 CALL reallocate(eam%phi, 1, eam%npoints)
2394 CALL reallocate(eam%phip, 1, eam%npoints)
2395 CALL reallocate(eam%frho, 1, eam%npoints)
2396 CALL reallocate(eam%frhop, 1, eam%npoints)
2397 ! Reading density and derivative of density (with respect to r)
2398 DO i = 1, eam%npoints
2399 CALL parser_get_next_line(parser, 1)
2400 READ (parser%input_line, *) eam%rho(i), eam%rhop(i)
2401 eam%rhop(i) = cp_unit_to_cp2k(eam%rhop(i), "angstrom^-1")
2402 eam%rval(i) = real(i - 1, kind=dp)*eam%drar
2403 eam%rhoval(i) = real(i - 1, kind=dp)*eam%drhoar
2404 END DO
2405 ! Reading pair potential PHI and its derivative (with respect to r)
2406 DO i = 1, eam%npoints
2407 CALL parser_get_next_line(parser, 1)
2408 READ (parser%input_line, *) eam%phi(i), eam%phip(i)
2409 eam%phi(i) = cp_unit_to_cp2k(eam%phi(i), "eV")
2410 eam%phip(i) = cp_unit_to_cp2k(eam%phip(i), "eV*angstrom^-1")
2411 END DO
2412 ! Reading embedded function and its derivative (with respect to density)
2413 DO i = 1, eam%npoints
2414 CALL parser_get_next_line(parser, 1)
2415 READ (parser%input_line, *) eam%frho(i), eam%frhop(i)
2416 eam%frho(i) = cp_unit_to_cp2k(eam%frho(i), "eV")
2417 eam%frhop(i) = cp_unit_to_cp2k(eam%frhop(i), "eV")
2418 END DO
2419
2420 IF (iw > 0) WRITE (iw, *) "Finished EAM data"
2421 CALL parser_release(parser)
2422 CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO")
2423 CALL timestop(handle)
2424
2425 END SUBROUTINE read_eam_data
2426
2427! **************************************************************************************************
2428!> \brief reads NequIP potential from .pth file
2429!> \param nequip ...
2430!> \author Gabriele Tocci
2431! **************************************************************************************************
2432 SUBROUTINE read_nequip_data(nequip)
2433 TYPE(nequip_pot_type), POINTER :: nequip
2434
2435 CHARACTER(len=*), PARAMETER :: routinen = 'read_nequip_data'
2436
2437 CHARACTER(LEN=default_path_length) :: allow_tf32_str, config_str, cutoff_str
2438 INTEGER :: handle
2439 LOGICAL :: allow_tf32, found
2440
2441 CALL timeset(routinen, handle)
2442
2443 INQUIRE (file=nequip%nequip_file_name, exist=found)
2444 IF (.NOT. found) THEN
2445 CALL cp_abort(__location__, &
2446 "Nequip model file <"//trim(nequip%nequip_file_name)// &
2447 "> not found.")
2448 END IF
2449
2450 nequip%nequip_version = torch_model_read_metadata(nequip%nequip_file_name, "nequip_version")
2451 cutoff_str = torch_model_read_metadata(nequip%nequip_file_name, "r_max")
2452 READ (cutoff_str, *) nequip%rcutsq
2453 nequip%rcutsq = cp_unit_to_cp2k(nequip%rcutsq, nequip%unit_coords)
2454 nequip%rcutsq = nequip%rcutsq*nequip%rcutsq
2455 nequip%unit_coords_val = cp_unit_to_cp2k(nequip%unit_coords_val, nequip%unit_coords)
2456 nequip%unit_forces_val = cp_unit_to_cp2k(nequip%unit_forces_val, nequip%unit_forces)
2457 nequip%unit_energy_val = cp_unit_to_cp2k(nequip%unit_energy_val, nequip%unit_energy)
2458 nequip%unit_cell_val = cp_unit_to_cp2k(nequip%unit_cell_val, nequip%unit_cell)
2459
2460 ! look in config which contains all the .yaml file options to see if we use float32 or float64
2461 config_str = torch_model_read_metadata(nequip%nequip_file_name, "config")
2462 CALL read_default_dtype(config_str, nequip%do_nequip_sp)
2463
2464 allow_tf32_str = torch_model_read_metadata(nequip%nequip_file_name, "allow_tf32")
2465 allow_tf32 = (trim(allow_tf32_str) == "1")
2466 IF (trim(allow_tf32_str) /= "1" .AND. trim(allow_tf32_str) /= "0") THEN
2467 CALL cp_abort(__location__, &
2468 "The value for allow_tf32 <"//trim(allow_tf32_str)// &
2469 "> is not supported. Check the .yaml and .pth files.")
2470 END IF
2471 CALL torch_allow_tf32(allow_tf32)
2472
2473 CALL timestop(handle)
2474 END SUBROUTINE read_nequip_data
2475
2476! **************************************************************************************************
2477!> \brief reads ALLEGRO potential from .pth file
2478!> \param allegro ...
2479!> \author Gabriele Tocci
2480! **************************************************************************************************
2481 SUBROUTINE read_allegro_data(allegro)
2482 TYPE(allegro_pot_type), POINTER :: allegro
2483
2484 CHARACTER(len=*), PARAMETER :: routinen = 'read_allegro_data'
2485
2486 CHARACTER(LEN=default_path_length) :: allow_tf32_str, config_str, cutoff_str
2487 INTEGER :: handle
2488 LOGICAL :: allow_tf32, found
2489
2490 CALL timeset(routinen, handle)
2491
2492 INQUIRE (file=allegro%allegro_file_name, exist=found)
2493 IF (.NOT. found) THEN
2494 CALL cp_abort(__location__, &
2495 "Allegro model file <"//trim(allegro%allegro_file_name)// &
2496 "> not found.")
2497 END IF
2498
2499 allegro%nequip_version = torch_model_read_metadata(allegro%allegro_file_name, "nequip_version")
2500 IF (allegro%nequip_version == "") THEN
2501 CALL cp_abort(__location__, &
2502 "Allegro model file <"//trim(allegro%allegro_file_name)// &
2503 "> has not been deployed; did you forget to run `nequip-deploy`?")
2504 END IF
2505 cutoff_str = torch_model_read_metadata(allegro%allegro_file_name, "r_max")
2506 READ (cutoff_str, *) allegro%rcutsq
2507 allegro%rcutsq = cp_unit_to_cp2k(allegro%rcutsq, allegro%unit_coords)
2508 allegro%rcutsq = allegro%rcutsq*allegro%rcutsq
2509 allegro%unit_coords_val = cp_unit_to_cp2k(allegro%unit_coords_val, allegro%unit_coords)
2510 allegro%unit_forces_val = cp_unit_to_cp2k(allegro%unit_forces_val, allegro%unit_forces)
2511 allegro%unit_energy_val = cp_unit_to_cp2k(allegro%unit_energy_val, allegro%unit_energy)
2512 allegro%unit_cell_val = cp_unit_to_cp2k(allegro%unit_cell_val, allegro%unit_cell)
2513
2514 ! look in config which contains all the .yaml file options to see if we use float32 or float64
2515 config_str = torch_model_read_metadata(allegro%allegro_file_name, "config")
2516 CALL read_default_dtype(config_str, allegro%do_allegro_sp)
2517
2518 allow_tf32_str = torch_model_read_metadata(allegro%allegro_file_name, "allow_tf32")
2519 allow_tf32 = (trim(allow_tf32_str) == "1")
2520 IF (trim(allow_tf32_str) /= "1" .AND. trim(allow_tf32_str) /= "0") THEN
2521 CALL cp_abort(__location__, &
2522 "The value for allow_tf32 <"//trim(allow_tf32_str)// &
2523 "> is not supported. Check the .yaml and .pth files.")
2524 END IF
2525 CALL torch_allow_tf32(allow_tf32)
2526
2527 CALL timestop(handle)
2528 END SUBROUTINE read_allegro_data
2529
2530! **************************************************************************************************
2531!> \brief reads the default_dtype used in the Allegro/NequIP model by parsing the config file
2532!> \param config_str ...
2533!> \param do_model_sp ...
2534!> \author Gabriele Tocci
2535! **************************************************************************************************
2536 SUBROUTINE read_default_dtype(config_str, do_model_sp)
2537
2538 CHARACTER(LEN=default_path_length) :: config_str
2539 LOGICAL :: do_model_sp
2540
2541 CHARACTER(len=*), PARAMETER :: routinen = 'read_default_dtype'
2542
2543 INTEGER :: handle, i, idx, len_config
2544
2545 CALL timeset(routinen, handle)
2546
2547 len_config = len_trim(config_str)
2548 idx = index(config_str, "default_dtype:")
2549 IF (idx /= 0) THEN
2550 i = idx + 14 ! skip over "default_dtype:"
2551 DO WHILE (i <= len_config .AND. config_str(i:i) == " ")
2552 i = i + 1 ! skip over any whitespace
2553 END DO
2554
2555 IF (i > len_config) THEN
2556 CALL cp_abort(__location__, &
2557 "No default_dtype found, check the Nequip/Allegro .yaml or .pth files."// &
2558 " Default_dtype should be either <float32> or <float64>.")
2559 ELSE IF (config_str(i:i + 6) == "float32") THEN
2560 do_model_sp = .true.
2561 ELSE IF (config_str(i:i + 6) == "float64") THEN
2562 do_model_sp = .false.
2563 ELSE
2564 CALL cp_abort(__location__, &
2565 "The default_dtype should be either <float32> or <float64>."// &
2566 " Check the NequIP/Allegro .yaml and .pth files.")
2567 END IF
2568 END IF
2569
2570 CALL timestop(handle)
2571 END SUBROUTINE read_default_dtype
2572
2573! **************************************************************************************************
2574!> \brief reads TABPOT potential from file
2575!> \param tab ...
2576!> \param para_env ...
2577!> \param mm_section ...
2578!> \author Da Teng, Alex Mironenko
2579! **************************************************************************************************
2580 SUBROUTINE read_tabpot_data(tab, para_env, mm_section)
2581 TYPE(tab_pot_type), POINTER :: tab
2582 TYPE(mp_para_env_type), POINTER :: para_env
2583 TYPE(section_vals_type), POINTER :: mm_section
2584
2585 CHARACTER(len=*), PARAMETER :: routinen = 'read_tabpot_data'
2586
2587 CHARACTER :: d1, d2
2588 INTEGER :: d, handle, i, iw
2589 TYPE(cp_logger_type), POINTER :: logger
2590 TYPE(cp_parser_type) :: parser
2591
2592 CALL timeset(routinen, handle)
2593 NULLIFY (logger)
2594 logger => cp_get_default_logger()
2595 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
2596 extension=".mmLog")
2597 IF (iw > 0) WRITE (iw, *) "Reading TABPOT data from: ", trim(tab%tabpot_file_name)
2598 CALL parser_create(parser, trim(tab%tabpot_file_name), para_env=para_env)
2599 CALL parser_get_next_line(parser, 1)
2600 IF (iw > 0) WRITE (iw, *) "Title: ", parser%input_line
2601 CALL parser_get_next_line(parser, 1)
2602
2603 ! example format: N 1000 R 1.00 20.0
2604 ! Assume the data is evenly spaced
2605 READ (parser%input_line, *) d1, tab%npoints, d2, tab%dr, tab%rcut
2606
2607 ! Relocating arrays with the right size
2608 CALL reallocate(tab%r, 1, tab%npoints)
2609 CALL reallocate(tab%e, 1, tab%npoints)
2610 CALL reallocate(tab%f, 1, tab%npoints)
2611
2612 ! Reading r, e, f
2613 DO i = 1, tab%npoints
2614 CALL parser_get_next_line(parser, 1)
2615 READ (parser%input_line, *) d, tab%r(i), tab%e(i), tab%f(i)
2616 tab%r(i) = cp_unit_to_cp2k(tab%r(i), "angstrom")
2617 tab%e(i) = cp_unit_to_cp2k(tab%e(i), "kcalmol")
2618 tab%f(i) = cp_unit_to_cp2k(tab%f(i), "kcalmol*angstrom^-1")
2619 END DO
2620
2621 tab%dr = tab%r(2) - tab%r(1)
2622 tab%rcut = cp_unit_to_cp2k(tab%rcut, "angstrom")
2623
2624 IF (iw > 0) WRITE (iw, *) "Finished TABPOT data"
2625 CALL parser_release(parser)
2626 CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO")
2627 CALL timestop(handle)
2628 END SUBROUTINE read_tabpot_data
2629
2630END 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:117
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:153
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
mathlib::invert_matrix
Definition mathlib.F:69
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:28
bibliography::tosi1964b
integer, save, public tosi1964b
Definition bibliography.F:43
bibliography::tersoff1988
integer, save, public tersoff1988
Definition bibliography.F:43
bibliography::tosi1964a
integer, save, public tosi1964a
Definition bibliography.F:43
bibliography::siepmann1995
integer, save, public siepmann1995
Definition bibliography.F:43
bibliography::yamada2000
integer, save, public yamada2000
Definition bibliography.F:43
bibliography::clabaut2021
integer, save, public clabaut2021
Definition bibliography.F:43
bibliography::clabaut2020
integer, save, public clabaut2020
Definition bibliography.F:43
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:510
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:803
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:1013
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:952
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:2345
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:1366
force_fields_input::read_wl_section
subroutine, public read_wl_section(nonbonded, section, start)
Reads the WILLIAMS section.
Definition force_fields_input.F:908
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:1842
force_fields_input::read_lj_section
subroutine, public read_lj_section(nonbonded, section, start)
Reads the LJ section.
Definition force_fields_input.F:865
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:1149
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:724
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:697
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:1040
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:1109
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:3361
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:1104
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:1028
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:126
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:203
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:407
pair_potential_types::tab_pot_type
Definition pair_potential_types.F:346
shell_potential_types::shell_p_type
Definition shell_potential_types.F:40