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