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