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