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