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