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