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