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