Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \par History
10 : !> Subroutine input_torsions changed (DG) 05-Dec-2000
11 : !> Output formats changed (DG) 05-Dec-2000
12 : !> JGH (26-01-2002) : force field parameters stored in tables, not in
13 : !> matrices. Input changed to have parameters labeled by the position
14 : !> and not atom pairs (triples etc)
15 : !> Teo (11.2005) : Moved all information on force field pair_potential to
16 : !> a much lighter memory structure
17 : !> \author CJM
18 : ! **************************************************************************************************
19 : MODULE force_fields_util
20 :
21 : USE atomic_kind_types, ONLY: atomic_kind_type,&
22 : get_atomic_kind
23 : USE cell_types, ONLY: cell_type
24 : USE colvar_types, ONLY: dist_colvar_id
25 : USE cp_log_handling, ONLY: cp_get_default_logger,&
26 : cp_logger_get_default_io_unit,&
27 : cp_logger_type
28 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
29 : cp_print_key_unit_nr
30 : USE cp_units, ONLY: cp_unit_to_cp2k
31 : USE ewald_environment_types, ONLY: ewald_env_get,&
32 : ewald_environment_type
33 : USE fist_nonbond_env_types, ONLY: fist_nonbond_env_create,&
34 : fist_nonbond_env_set,&
35 : fist_nonbond_env_type
36 : USE force_field_kind_types, ONLY: &
37 : allocate_bend_kind_set, allocate_bond_kind_set, allocate_impr_kind_set, &
38 : allocate_opbend_kind_set, allocate_torsion_kind_set, allocate_ub_kind_set, bend_kind_type, &
39 : bond_kind_type, deallocate_bend_kind_set, deallocate_bond_kind_set, do_ff_undef, &
40 : impr_kind_dealloc_ref, impr_kind_type, opbend_kind_type, torsion_kind_dealloc_ref, &
41 : torsion_kind_type, ub_kind_dealloc_ref, ub_kind_type
42 : USE force_field_types, ONLY: amber_info_type,&
43 : charmm_info_type,&
44 : force_field_type,&
45 : gromos_info_type,&
46 : input_info_type
47 : USE force_fields_all, ONLY: &
48 : force_field_pack_bend, force_field_pack_bond, force_field_pack_charge, &
49 : force_field_pack_charges, force_field_pack_damp, force_field_pack_eicut, &
50 : force_field_pack_impr, force_field_pack_nonbond, force_field_pack_nonbond14, &
51 : force_field_pack_opbend, force_field_pack_pol, force_field_pack_radius, &
52 : force_field_pack_shell, force_field_pack_splines, force_field_pack_tors, &
53 : force_field_pack_ub, force_field_unique_bend, force_field_unique_bond, &
54 : force_field_unique_impr, force_field_unique_opbend, force_field_unique_tors, &
55 : force_field_unique_ub
56 : USE input_section_types, ONLY: section_vals_get,&
57 : section_vals_get_subs_vals,&
58 : section_vals_type,&
59 : section_vals_val_get
60 : USE kinds, ONLY: default_path_length,&
61 : default_string_length,&
62 : dp
63 : USE memory_utilities, ONLY: reallocate
64 : USE molecule_kind_types, ONLY: &
65 : atom_type, bend_type, bond_type, colvar_constraint_type, g3x3_constraint_type, &
66 : g4x6_constraint_type, get_molecule_kind, impr_type, molecule_kind_type, opbend_type, &
67 : set_molecule_kind, torsion_type, ub_type
68 : USE molecule_types, ONLY: get_molecule,&
69 : molecule_type
70 : USE pair_potential_types, ONLY: pair_potential_pp_type
71 : USE particle_types, ONLY: particle_type
72 : USE qmmm_types_low, ONLY: qmmm_env_mm_type
73 : USE shell_potential_types, ONLY: shell_kind_type
74 : USE string_utilities, ONLY: compress
75 : #include "./base/base_uses.f90"
76 :
77 : IMPLICIT NONE
78 :
79 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_fields_util'
80 :
81 : PRIVATE
82 : PUBLIC :: force_field_pack, &
83 : force_field_qeff_output, &
84 : clean_intra_force_kind, &
85 : get_generic_info
86 :
87 : CONTAINS
88 :
89 : ! **************************************************************************************************
90 : !> \brief Pack in all the information needed for the force_field
91 : !> \param particle_set ...
92 : !> \param atomic_kind_set ...
93 : !> \param molecule_kind_set ...
94 : !> \param molecule_set ...
95 : !> \param ewald_env ...
96 : !> \param fist_nonbond_env ...
97 : !> \param ff_type ...
98 : !> \param root_section ...
99 : !> \param qmmm ...
100 : !> \param qmmm_env ...
101 : !> \param mm_section ...
102 : !> \param subsys_section ...
103 : !> \param shell_particle_set ...
104 : !> \param core_particle_set ...
105 : !> \param cell ...
106 : ! **************************************************************************************************
107 10556 : SUBROUTINE force_field_pack(particle_set, atomic_kind_set, molecule_kind_set, &
108 : molecule_set, ewald_env, fist_nonbond_env, ff_type, root_section, qmmm, &
109 : qmmm_env, mm_section, subsys_section, shell_particle_set, core_particle_set, &
110 : cell)
111 :
112 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
113 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
114 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
115 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
116 : TYPE(ewald_environment_type), POINTER :: ewald_env
117 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
118 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
119 : TYPE(section_vals_type), POINTER :: root_section
120 : LOGICAL, INTENT(IN), OPTIONAL :: qmmm
121 : TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env
122 : TYPE(section_vals_type), POINTER :: mm_section, subsys_section
123 : TYPE(particle_type), DIMENSION(:), POINTER :: shell_particle_set, core_particle_set
124 : TYPE(cell_type), POINTER :: cell
125 :
126 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack'
127 :
128 : CHARACTER(LEN=default_string_length), &
129 2639 : DIMENSION(:), POINTER :: Ainfo
130 : INTEGER :: handle, iw, iw2, iw3, iw4, output_unit
131 : LOGICAL :: do_zbl, explicit, fatal, ignore_fatal, &
132 : my_qmmm
133 : REAL(KIND=dp) :: ewald_rcut, verlet_skin
134 2639 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
135 : TYPE(amber_info_type), POINTER :: amb_info
136 : TYPE(charmm_info_type), POINTER :: chm_info
137 : TYPE(cp_logger_type), POINTER :: logger
138 : TYPE(gromos_info_type), POINTER :: gro_info
139 : TYPE(input_info_type), POINTER :: inp_info
140 : TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond, potparm_nonbond14
141 : TYPE(section_vals_type), POINTER :: charges_section
142 :
143 2639 : CALL timeset(routineN, handle)
144 2639 : fatal = .FALSE.
145 2639 : ignore_fatal = ff_type%ignore_missing_critical
146 2639 : NULLIFY (logger, Ainfo, charges_section, charges)
147 2639 : logger => cp_get_default_logger()
148 : ! Error unit
149 2639 : output_unit = cp_logger_get_default_io_unit(logger)
150 :
151 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
152 2639 : extension=".mmLog")
153 : iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO/SPLINE_INFO", &
154 2639 : extension=".mmLog")
155 : iw3 = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO/SPLINE_DATA", &
156 2639 : extension=".mmLog")
157 : iw4 = cp_print_key_unit_nr(logger, mm_section, "PRINT%PROGRAM_RUN_INFO", &
158 2639 : extension=".mmLog")
159 2639 : NULLIFY (potparm_nonbond14, potparm_nonbond)
160 2639 : my_qmmm = .FALSE.
161 2639 : IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
162 2639 : inp_info => ff_type%inp_info
163 2639 : chm_info => ff_type%chm_info
164 2639 : gro_info => ff_type%gro_info
165 2639 : amb_info => ff_type%amb_info
166 : !-----------------------------------------------------------------------------
167 : !-----------------------------------------------------------------------------
168 : ! 1. Determine the number of unique bond kind and allocate bond_kind_set
169 : !-----------------------------------------------------------------------------
170 : CALL force_field_unique_bond(particle_set, molecule_kind_set, molecule_set, &
171 2639 : ff_type)
172 : !-----------------------------------------------------------------------------
173 : !-----------------------------------------------------------------------------
174 : ! 2. Determine the number of unique bend kind and allocate bend_kind_set
175 : !-----------------------------------------------------------------------------
176 : CALL force_field_unique_bend(particle_set, molecule_kind_set, molecule_set, &
177 2639 : ff_type)
178 : !-----------------------------------------------------------------------------
179 : !-----------------------------------------------------------------------------
180 : ! 3. Determine the number of unique Urey-Bradley kind and allocate ub_kind_set
181 : !-----------------------------------------------------------------------------
182 2639 : CALL force_field_unique_ub(particle_set, molecule_kind_set, molecule_set)
183 : !-----------------------------------------------------------------------------
184 : !-----------------------------------------------------------------------------
185 : ! 4. Determine the number of unique torsion kind and allocate torsion_kind_set
186 : !-----------------------------------------------------------------------------
187 : CALL force_field_unique_tors(particle_set, molecule_kind_set, molecule_set, &
188 2639 : ff_type)
189 : !-----------------------------------------------------------------------------
190 : !-----------------------------------------------------------------------------
191 : ! 5. Determine the number of unique impr kind and allocate impr_kind_set
192 : !-----------------------------------------------------------------------------
193 : CALL force_field_unique_impr(particle_set, molecule_kind_set, molecule_set, &
194 2639 : ff_type)
195 : !-----------------------------------------------------------------------------
196 : !-----------------------------------------------------------------------------
197 : ! 6. Determine the number of unique opbend kind and allocate opbend_kind_set
198 : !-----------------------------------------------------------------------------
199 : CALL force_field_unique_opbend(particle_set, molecule_kind_set, molecule_set, &
200 2639 : ff_type)
201 : !-----------------------------------------------------------------------------
202 : !-----------------------------------------------------------------------------
203 : ! 7. Bonds
204 : !-----------------------------------------------------------------------------
205 : CALL force_field_pack_bond(particle_set, molecule_kind_set, molecule_set, &
206 2639 : fatal, Ainfo, chm_info, inp_info, gro_info, amb_info)
207 : !-----------------------------------------------------------------------------
208 : !-----------------------------------------------------------------------------
209 : ! 8. Bends
210 : !-----------------------------------------------------------------------------
211 : CALL force_field_pack_bend(particle_set, molecule_kind_set, molecule_set, &
212 2639 : fatal, Ainfo, chm_info, inp_info, gro_info, amb_info)
213 : ! Give information and abort if any bond or angle parameter is missing..
214 2639 : CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw)
215 : !-----------------------------------------------------------------------------
216 : !-----------------------------------------------------------------------------
217 : ! 9. Urey-Bradley
218 : !-----------------------------------------------------------------------------
219 : CALL force_field_pack_ub(particle_set, molecule_kind_set, molecule_set, &
220 2639 : Ainfo, chm_info, inp_info, iw)
221 : !-----------------------------------------------------------------------------
222 : !-----------------------------------------------------------------------------
223 : ! 10. Torsions
224 : !-----------------------------------------------------------------------------
225 : CALL force_field_pack_tors(particle_set, molecule_kind_set, molecule_set, &
226 2639 : Ainfo, chm_info, inp_info, gro_info, amb_info, iw)
227 : !-----------------------------------------------------------------------------
228 : !-----------------------------------------------------------------------------
229 : ! 11. Impropers
230 : !-----------------------------------------------------------------------------
231 : CALL force_field_pack_impr(particle_set, molecule_kind_set, molecule_set, &
232 2639 : Ainfo, chm_info, inp_info, gro_info)
233 : !-----------------------------------------------------------------------------
234 : !-----------------------------------------------------------------------------
235 : ! 12. Out of plane bends
236 : !-----------------------------------------------------------------------------
237 : CALL force_field_pack_opbend(particle_set, molecule_kind_set, molecule_set, &
238 2639 : Ainfo, inp_info)
239 : ! Give information only if any Urey-Bradley, Torsion, improper or opbend is missing
240 : ! continue calculation..
241 2639 : CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw)
242 :
243 2639 : charges_section => section_vals_get_subs_vals(mm_section, "FORCEFIELD%CHARGES")
244 2639 : CALL section_vals_get(charges_section, explicit=explicit)
245 2639 : IF (.NOT. explicit) THEN
246 : !-----------------------------------------------------------------------------
247 : !-----------------------------------------------------------------------------
248 : ! 13.a Set up atomic_kind_set()%fist_potentail%[qeff] and shell
249 : ! potential parameters
250 : !-----------------------------------------------------------------------------
251 : CALL force_field_pack_charge(atomic_kind_set, qmmm_env, fatal, iw, iw4, &
252 2631 : Ainfo, my_qmmm, inp_info)
253 : ! Give information only if charge is missing and abort..
254 2631 : CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw)
255 : ELSE
256 : !-----------------------------------------------------------------------------
257 : !-----------------------------------------------------------------------------
258 : ! 13.b Setup a global array of classical charges - this avoids the packing and
259 : ! allows the usage of different charges for same atomic types
260 : !-----------------------------------------------------------------------------
261 : CALL force_field_pack_charges(charges, charges_section, particle_set, my_qmmm, &
262 8 : qmmm_env, inp_info, iw4)
263 : END IF
264 :
265 : !-----------------------------------------------------------------------------
266 : !-----------------------------------------------------------------------------
267 : ! 14. Set up the radius of the electrostatic multipole in Fist
268 : !-----------------------------------------------------------------------------
269 2639 : CALL force_field_pack_radius(atomic_kind_set, iw, subsys_section)
270 :
271 : !-----------------------------------------------------------------------------
272 : !-----------------------------------------------------------------------------
273 : ! 15. Set up the polarizable FF parameters
274 : !-----------------------------------------------------------------------------
275 2639 : CALL force_field_pack_pol(atomic_kind_set, iw, inp_info)
276 2639 : CALL force_field_pack_damp(atomic_kind_set, iw, inp_info)
277 :
278 : !-----------------------------------------------------------------------------
279 : !-----------------------------------------------------------------------------
280 : ! 16. Set up Shell potential
281 : !-----------------------------------------------------------------------------
282 : CALL force_field_pack_shell(particle_set, atomic_kind_set, &
283 : molecule_kind_set, molecule_set, root_section, subsys_section, &
284 2639 : shell_particle_set, core_particle_set, cell, iw, inp_info)
285 2639 : IF (ff_type%do_nonbonded) THEN
286 : !-----------------------------------------------------------------------------
287 : !-----------------------------------------------------------------------------
288 : ! 17. Set up potparm_nonbond14
289 : !-----------------------------------------------------------------------------
290 : ! Move the data from the info structures to potparm_nonbond
291 : CALL force_field_pack_nonbond14(atomic_kind_set, ff_type, qmmm_env, iw, Ainfo, &
292 2623 : chm_info, inp_info, gro_info, amb_info, potparm_nonbond14, ewald_env)
293 : ! Give information if any 1-4 is missing.. continue calculation..
294 2623 : CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw)
295 : ! Create the spline data
296 2623 : CALL section_vals_val_get(mm_section, "FORCEFIELD%ZBL_SCATTERING", l_val=do_zbl)
297 : CALL force_field_pack_splines(atomic_kind_set, ff_type, iw2, iw3, iw4, &
298 2623 : potparm_nonbond14, do_zbl, nonbonded_type="NONBONDED14")
299 : !-----------------------------------------------------------------------------
300 : !-----------------------------------------------------------------------------
301 : ! 18. Set up potparm_nonbond
302 : !-----------------------------------------------------------------------------
303 : ! Move the data from the info structures to potparm_nonbond
304 : CALL force_field_pack_nonbond(atomic_kind_set, ff_type, qmmm_env, &
305 : fatal, iw, Ainfo, chm_info, inp_info, gro_info, amb_info, &
306 2623 : potparm_nonbond, ewald_env)
307 : ! Give information and abort if any pair potential spline is missing..
308 2623 : CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw)
309 : ! Create the spline data
310 2623 : CALL section_vals_val_get(mm_section, "FORCEFIELD%ZBL_SCATTERING", l_val=do_zbl)
311 : CALL force_field_pack_splines(atomic_kind_set, ff_type, iw2, iw3, iw4, &
312 2623 : potparm_nonbond, do_zbl, nonbonded_type="NONBONDED")
313 : END IF
314 : !-----------------------------------------------------------------------------
315 : !-----------------------------------------------------------------------------
316 : ! 19. Create nonbond environment
317 : !-----------------------------------------------------------------------------
318 2639 : CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
319 : CALL section_vals_val_get(mm_section, "NEIGHBOR_LISTS%VERLET_SKIN", &
320 2639 : r_val=verlet_skin)
321 2639 : ALLOCATE (fist_nonbond_env)
322 : CALL fist_nonbond_env_create(fist_nonbond_env, atomic_kind_set, &
323 : potparm_nonbond14, potparm_nonbond, ff_type%do_nonbonded, &
324 : ff_type%do_electrostatics, verlet_skin, ewald_rcut, ff_type%ei_scale14, &
325 2639 : ff_type%vdw_scale14, ff_type%shift_cutoff)
326 2639 : CALL fist_nonbond_env_set(fist_nonbond_env, charges=charges)
327 : ! Compute the electrostatic interaction cutoffs.
328 : CALL force_field_pack_eicut(atomic_kind_set, ff_type, potparm_nonbond, &
329 2639 : ewald_env)
330 :
331 2639 : CALL cp_print_key_finished_output(iw4, logger, mm_section, "PRINT%PROGRAM_RUN_INFO")
332 2639 : CALL cp_print_key_finished_output(iw3, logger, mm_section, "PRINT%FF_INFO/SPLINE_DATA")
333 2639 : CALL cp_print_key_finished_output(iw2, logger, mm_section, "PRINT%FF_INFO/SPLINE_INFO")
334 2639 : CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO")
335 2639 : CALL timestop(handle)
336 :
337 2639 : END SUBROUTINE force_field_pack
338 :
339 : ! **************************************************************************************************
340 : !> \brief Store informations on possible missing ForceFields parameters
341 : !> \param fatal ...
342 : !> \param ignore_fatal ...
343 : !> \param array ...
344 : !> \param output_unit ...
345 : !> \param iw ...
346 : ! **************************************************************************************************
347 13155 : SUBROUTINE release_FF_missing_par(fatal, ignore_fatal, array, output_unit, iw)
348 : LOGICAL, INTENT(INOUT) :: fatal, ignore_fatal
349 : CHARACTER(LEN=default_string_length), &
350 : DIMENSION(:), POINTER :: array
351 : INTEGER, INTENT(IN) :: output_unit, iw
352 :
353 : INTEGER :: i
354 :
355 13155 : IF (ASSOCIATED(array)) THEN
356 3228 : IF (output_unit > 0) THEN
357 1971 : WRITE (output_unit, *)
358 : WRITE (output_unit, '(T2,"FORCEFIELD|",A)') &
359 1971 : " WARNING: A non Critical ForceField parameter is missing! CP2K GOES!", &
360 1971 : " Non critical parameters are:Urey-Bradley,Torsions,Impropers, Opbends and 1-4", &
361 3942 : " All missing parameters will not contribute to the potential energy!"
362 1971 : IF (fatal .OR. iw > 0) THEN
363 309 : WRITE (output_unit, *)
364 3029 : DO i = 1, SIZE(array)
365 4691 : WRITE (output_unit, '(A)') array(i)
366 : END DO
367 : END IF
368 1971 : IF (.NOT. fatal .AND. iw < 0) THEN
369 : WRITE (output_unit, '(T2,"FORCEFIELD|",A)') &
370 1662 : " Activate the print key FF_INFO to have a list of missing parameters"
371 1662 : WRITE (output_unit, *)
372 : END IF
373 : END IF
374 3228 : DEALLOCATE (array)
375 : END IF
376 13155 : IF (fatal) THEN
377 64 : IF (ignore_fatal) THEN
378 64 : IF (output_unit > 0) THEN
379 44 : WRITE (output_unit, *)
380 : WRITE (output_unit, '(T2,"FORCEFIELD|",A)') &
381 44 : " WARNING: Ignoring missing critical FF parameters! CP2K GOES!", &
382 44 : " Critical parameters are: Bonds, Bends and Charges", &
383 88 : " All missing parameters will not contribute to the potential energy!"
384 : END IF
385 : ELSE
386 0 : CPABORT("Missing critical ForceField parameters!")
387 : END IF
388 : END IF
389 13155 : END SUBROUTINE release_FF_missing_par
390 :
391 : ! **************************************************************************************************
392 : !> \brief Compute the total qeff charges for each molecule kind and total system
393 : !> \param particle_set ...
394 : !> \param molecule_kind_set ...
395 : !> \param molecule_set ...
396 : !> \param mm_section ...
397 : !> \param charges ...
398 : ! **************************************************************************************************
399 2639 : SUBROUTINE force_field_qeff_output(particle_set, molecule_kind_set, &
400 : molecule_set, mm_section, charges)
401 :
402 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
403 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
404 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
405 : TYPE(section_vals_type), POINTER :: mm_section
406 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
407 :
408 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_qeff_output'
409 :
410 : CHARACTER(LEN=default_string_length) :: atmname, molname
411 : INTEGER :: first, handle, iatom, imol, iw, j, jatom
412 : LOGICAL :: shell_active
413 : REAL(KIND=dp) :: mass, mass_mol, mass_sum, qeff, &
414 : qeff_mol, qeff_sum
415 2639 : TYPE(atom_type), DIMENSION(:), POINTER :: atom_list
416 : TYPE(atomic_kind_type), POINTER :: atomic_kind
417 : TYPE(cp_logger_type), POINTER :: logger
418 : TYPE(molecule_kind_type), POINTER :: molecule_kind
419 : TYPE(molecule_type), POINTER :: molecule
420 : TYPE(shell_kind_type), POINTER :: shell
421 :
422 2639 : CALL timeset(routineN, handle)
423 2639 : NULLIFY (logger)
424 2639 : logger => cp_get_default_logger()
425 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
426 2639 : extension=".mmLog")
427 :
428 2639 : qeff = 0.0_dp
429 2639 : qeff_mol = 0.0_dp
430 2639 : qeff_sum = 0.0_dp
431 2639 : mass_sum = 0.0_dp
432 :
433 : !-----------------------------------------------------------------------------
434 : !-----------------------------------------------------------------------------
435 : ! 1. Sum of [qeff,mass] for each molecule_kind
436 : !-----------------------------------------------------------------------------
437 74487 : DO imol = 1, SIZE(molecule_kind_set)
438 71848 : qeff_mol = 0.0_dp
439 71848 : mass_mol = 0.0_dp
440 71848 : molecule_kind => molecule_kind_set(imol)
441 :
442 71848 : j = molecule_kind_set(imol)%molecule_list(1)
443 71848 : molecule => molecule_set(j)
444 71848 : CALL get_molecule(molecule=molecule, first_atom=first)
445 :
446 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
447 71848 : name=molname, atom_list=atom_list)
448 266603 : DO iatom = 1, SIZE(atom_list)
449 194755 : atomic_kind => atom_list(iatom)%atomic_kind
450 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
451 194755 : name=atmname, qeff=qeff, mass=mass, shell_active=shell_active, shell=shell)
452 194755 : IF (shell_active) qeff = shell%charge_core + shell%charge_shell
453 194755 : IF (ASSOCIATED(charges)) THEN
454 30 : jatom = first - 1 + iatom
455 30 : qeff = charges(jatom)
456 : END IF
457 194755 : IF (iw > 0) WRITE (iw, *) " atom ", iatom, " ", TRIM(atmname), " charge = ", qeff
458 194755 : qeff_mol = qeff_mol + qeff
459 461358 : mass_mol = mass_mol + mass
460 : END DO
461 71848 : CALL set_molecule_kind(molecule_kind=molecule_kind, charge=qeff_mol, mass=mass_mol)
462 74487 : IF (iw > 0) WRITE (iw, *) " Mol Kind ", TRIM(molname), " charge = ", qeff_mol
463 : END DO
464 :
465 : !-----------------------------------------------------------------------------
466 : !-----------------------------------------------------------------------------
467 : ! 2. Sum of [qeff,mass] for particle_set
468 : !-----------------------------------------------------------------------------
469 660127 : DO iatom = 1, SIZE(particle_set)
470 657488 : atomic_kind => particle_set(iatom)%atomic_kind
471 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
472 657488 : name=atmname, qeff=qeff, mass=mass, shell_active=shell_active, shell=shell)
473 657488 : IF (shell_active) qeff = shell%charge_core + shell%charge_shell
474 657488 : IF (ASSOCIATED(charges)) THEN
475 36 : qeff = charges(iatom)
476 : END IF
477 666882 : IF (iw > 0) WRITE (iw, *) " atom ", iatom, " ", TRIM(atmname), &
478 18788 : " charge = ", qeff
479 657488 : qeff_sum = qeff_sum + qeff
480 1317615 : mass_sum = mass_sum + mass
481 : END DO
482 2639 : IF (iw > 0) WRITE (iw, '(A,F20.10)') " Total system charge = ", qeff_sum
483 2639 : IF (iw > 0) WRITE (iw, '(A,F20.10)') " Total system mass = ", mass_sum
484 :
485 : CALL cp_print_key_finished_output(iw, logger, mm_section, &
486 2639 : "PRINT%FF_INFO")
487 2639 : CALL timestop(handle)
488 2639 : END SUBROUTINE force_field_qeff_output
489 :
490 : ! **************************************************************************************************
491 : !> \brief Removes UNSET force field types
492 : !> \param molecule_kind_set ...
493 : !> \param mm_section ...
494 : ! **************************************************************************************************
495 2639 : SUBROUTINE clean_intra_force_kind(molecule_kind_set, mm_section)
496 :
497 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
498 : TYPE(section_vals_type), POINTER :: mm_section
499 :
500 : CHARACTER(len=*), PARAMETER :: routineN = 'clean_intra_force_kind'
501 :
502 : INTEGER :: atm2_a, atm2_b, atm2_c, atm_a, atm_b, atm_c, atm_d, counter, handle, i, ibend, &
503 : ibond, icolv, ig3x3, ig4x6, iimpr, ikind, iopbend, itorsion, iub, iw, j, k, nbend, nbond, &
504 : newkind, ng3x3, ng4x6, nimpr, nopbend, ntorsion, nub
505 2639 : INTEGER, POINTER :: bad1(:), bad2(:)
506 : LOGICAL :: unsetme, valid_kind
507 2639 : TYPE(bend_kind_type), DIMENSION(:), POINTER :: bend_kind_set, new_bend_kind_set
508 2639 : TYPE(bend_type), DIMENSION(:), POINTER :: bend_list, new_bend_list
509 2639 : TYPE(bond_kind_type), DIMENSION(:), POINTER :: bond_kind_set, new_bond_kind_set
510 2639 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list, new_bond_list
511 : TYPE(colvar_constraint_type), DIMENSION(:), &
512 2639 : POINTER :: colv_list
513 : TYPE(cp_logger_type), POINTER :: logger
514 2639 : TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
515 2639 : TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
516 2639 : TYPE(impr_kind_type), DIMENSION(:), POINTER :: impr_kind_set, new_impr_kind_set
517 2639 : TYPE(impr_type), DIMENSION(:), POINTER :: impr_list, new_impr_list
518 : TYPE(molecule_kind_type), POINTER :: molecule_kind
519 2639 : TYPE(opbend_kind_type), DIMENSION(:), POINTER :: new_opbend_kind_set, opbend_kind_set
520 2639 : TYPE(opbend_type), DIMENSION(:), POINTER :: new_opbend_list, opbend_list
521 2639 : TYPE(torsion_kind_type), DIMENSION(:), POINTER :: new_torsion_kind_set, torsion_kind_set
522 2639 : TYPE(torsion_type), DIMENSION(:), POINTER :: new_torsion_list, torsion_list
523 2639 : TYPE(ub_kind_type), DIMENSION(:), POINTER :: new_ub_kind_set, ub_kind_set
524 2639 : TYPE(ub_type), DIMENSION(:), POINTER :: new_ub_list, ub_list
525 :
526 2639 : CALL timeset(routineN, handle)
527 2639 : NULLIFY (logger)
528 2639 : logger => cp_get_default_logger()
529 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
530 2639 : extension=".mmLog")
531 : !-----------------------------------------------------------------------------
532 : !-----------------------------------------------------------------------------
533 : ! 1. Lets Tag the unwanted bonds due to the use of distance constraint
534 : !-----------------------------------------------------------------------------
535 74487 : DO ikind = 1, SIZE(molecule_kind_set)
536 71848 : molecule_kind => molecule_kind_set(ikind)
537 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
538 : colv_list=colv_list, &
539 : nbond=nbond, &
540 71848 : bond_list=bond_list)
541 74487 : IF (ASSOCIATED(colv_list)) THEN
542 804 : DO icolv = 1, SIZE(colv_list)
543 382 : IF ((colv_list(icolv)%type_id == dist_colvar_id) .AND. &
544 422 : ((.NOT. colv_list(icolv)%use_points) .OR. (SIZE(colv_list(icolv)%i_atoms) == 2))) THEN
545 284 : atm_a = colv_list(icolv)%i_atoms(1)
546 284 : atm_b = colv_list(icolv)%i_atoms(2)
547 3156 : DO ibond = 1, nbond
548 2872 : unsetme = .FALSE.
549 2872 : atm2_a = bond_list(ibond)%a
550 2872 : atm2_b = bond_list(ibond)%b
551 2872 : IF (atm2_a == atm_a .AND. atm2_b == atm_b) unsetme = .TRUE.
552 2872 : IF (atm2_a == atm_b .AND. atm2_b == atm_a) unsetme = .TRUE.
553 3254 : IF (unsetme) bond_list(ibond)%id_type = do_ff_undef
554 : END DO
555 : END IF
556 : END DO
557 : END IF
558 : END DO
559 : !-----------------------------------------------------------------------------
560 : !-----------------------------------------------------------------------------
561 : ! 2. Lets Tag the unwanted bends due to the use of distance constraint
562 : !-----------------------------------------------------------------------------
563 74487 : DO ikind = 1, SIZE(molecule_kind_set)
564 71848 : molecule_kind => molecule_kind_set(ikind)
565 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
566 : colv_list=colv_list, &
567 : nbend=nbend, &
568 71848 : bend_list=bend_list)
569 74487 : IF (ASSOCIATED(colv_list)) THEN
570 1576 : DO ibend = 1, nbend
571 1154 : unsetme = .FALSE.
572 1154 : atm_a = bend_list(ibend)%a
573 1154 : atm_b = bend_list(ibend)%b
574 1154 : atm_c = bend_list(ibend)%c
575 6774 : DO icolv = 1, SIZE(colv_list)
576 5656 : IF ((colv_list(icolv)%type_id == dist_colvar_id) .AND. &
577 1118 : ((.NOT. colv_list(icolv)%use_points) .OR. (SIZE(colv_list(icolv)%i_atoms) == 2))) THEN
578 5030 : atm2_a = colv_list(icolv)%i_atoms(1)
579 5030 : atm2_b = colv_list(icolv)%i_atoms(2)
580 : ! Check that the bonds we're constraining does not involve atoms defining
581 : ! a bend..
582 5030 : IF (((atm2_a == atm_a) .AND. (atm2_b == atm_c)) .OR. &
583 : ((atm2_a == atm_c) .AND. (atm2_b == atm_a))) THEN
584 : unsetme = .TRUE.
585 : EXIT
586 : END IF
587 : END IF
588 : END DO
589 1576 : IF (unsetme) bend_list(ibend)%id_type = do_ff_undef
590 : END DO
591 : END IF
592 : END DO
593 : !-----------------------------------------------------------------------------
594 : !-----------------------------------------------------------------------------
595 : ! 3. Lets Tag the unwanted bonds due to the use of 3x3
596 : !-----------------------------------------------------------------------------
597 74487 : DO ikind = 1, SIZE(molecule_kind_set)
598 71848 : molecule_kind => molecule_kind_set(ikind)
599 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
600 : ng3x3=ng3x3, &
601 : g3x3_list=g3x3_list, &
602 : nbond=nbond, &
603 71848 : bond_list=bond_list)
604 74625 : DO ig3x3 = 1, ng3x3
605 138 : atm_a = g3x3_list(ig3x3)%a
606 138 : atm_b = g3x3_list(ig3x3)%b
607 138 : atm_c = g3x3_list(ig3x3)%c
608 72302 : DO ibond = 1, nbond
609 316 : unsetme = .FALSE.
610 316 : atm2_a = bond_list(ibond)%a
611 316 : atm2_b = bond_list(ibond)%b
612 316 : IF (atm2_a == atm_a .AND. atm2_b == atm_b) unsetme = .TRUE.
613 316 : IF (atm2_a == atm_a .AND. atm2_b == atm_c) unsetme = .TRUE.
614 316 : IF (atm2_a == atm_c .AND. atm2_b == atm_c) unsetme = .TRUE.
615 454 : IF (unsetme) bond_list(ibond)%id_type = do_ff_undef
616 : END DO
617 : END DO
618 : END DO
619 : !-----------------------------------------------------------------------------
620 : !-----------------------------------------------------------------------------
621 : ! 4. Lets Tag the unwanted bends due to the use of 3x3
622 : !-----------------------------------------------------------------------------
623 74487 : DO ikind = 1, SIZE(molecule_kind_set)
624 71848 : molecule_kind => molecule_kind_set(ikind)
625 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
626 : ng3x3=ng3x3, &
627 : g3x3_list=g3x3_list, &
628 : nbend=nbend, &
629 71848 : bend_list=bend_list)
630 74625 : DO ig3x3 = 1, ng3x3
631 138 : atm_a = g3x3_list(ig3x3)%a
632 138 : atm_b = g3x3_list(ig3x3)%b
633 138 : atm_c = g3x3_list(ig3x3)%c
634 72108 : DO ibend = 1, nbend
635 122 : unsetme = .FALSE.
636 122 : atm2_a = bend_list(ibend)%a
637 122 : atm2_b = bend_list(ibend)%b
638 122 : atm2_c = bend_list(ibend)%c
639 122 : IF (atm2_a == atm_a .AND. atm2_b == atm_b .AND. atm2_c == atm_c) unsetme = .TRUE.
640 122 : IF (atm2_a == atm_a .AND. atm2_b == atm_c .AND. atm2_c == atm_b) unsetme = .TRUE.
641 122 : IF (atm2_a == atm_b .AND. atm2_b == atm_a .AND. atm2_c == atm_c) unsetme = .TRUE.
642 122 : IF (atm2_a == atm_b .AND. atm2_b == atm_c .AND. atm2_c == atm_a) unsetme = .TRUE.
643 260 : IF (unsetme) bend_list(ibend)%id_type = do_ff_undef
644 : END DO
645 : END DO
646 : END DO
647 : !-----------------------------------------------------------------------------
648 : !-----------------------------------------------------------------------------
649 : ! 5. Lets Tag the unwanted bonds due to the use of 4x6
650 : !-----------------------------------------------------------------------------
651 74487 : DO ikind = 1, SIZE(molecule_kind_set)
652 71848 : molecule_kind => molecule_kind_set(ikind)
653 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
654 : ng4x6=ng4x6, &
655 : g4x6_list=g4x6_list, &
656 : nbond=nbond, &
657 71848 : bond_list=bond_list)
658 74499 : DO ig4x6 = 1, ng4x6
659 12 : atm_a = g4x6_list(ig4x6)%a
660 12 : atm_b = g4x6_list(ig4x6)%b
661 12 : atm_c = g4x6_list(ig4x6)%c
662 12 : atm_d = g4x6_list(ig4x6)%d
663 71896 : DO ibond = 1, nbond
664 36 : unsetme = .FALSE.
665 36 : atm2_a = bond_list(ibond)%a
666 36 : atm2_b = bond_list(ibond)%b
667 36 : IF (atm2_a == atm_a .AND. atm2_b == atm_b) unsetme = .TRUE.
668 36 : IF (atm2_a == atm_a .AND. atm2_b == atm_c) unsetme = .TRUE.
669 36 : IF (atm2_a == atm_a .AND. atm2_b == atm_d) unsetme = .TRUE.
670 48 : IF (unsetme) bond_list(ibond)%id_type = do_ff_undef
671 : END DO
672 : END DO
673 : END DO
674 : !-----------------------------------------------------------------------------
675 : !-----------------------------------------------------------------------------
676 : ! 6. Lets Tag the unwanted bends due to the use of 4x6
677 : !-----------------------------------------------------------------------------
678 74487 : DO ikind = 1, SIZE(molecule_kind_set)
679 71848 : molecule_kind => molecule_kind_set(ikind)
680 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
681 : ng4x6=ng4x6, &
682 : g4x6_list=g4x6_list, &
683 : nbend=nbend, &
684 71848 : bend_list=bend_list)
685 74499 : DO ig4x6 = 1, ng4x6
686 12 : atm_a = g4x6_list(ig4x6)%a
687 12 : atm_b = g4x6_list(ig4x6)%b
688 12 : atm_c = g4x6_list(ig4x6)%c
689 12 : atm_d = g4x6_list(ig4x6)%d
690 71896 : DO ibend = 1, nbend
691 36 : unsetme = .FALSE.
692 36 : atm2_a = bend_list(ibend)%a
693 36 : atm2_b = bend_list(ibend)%b
694 36 : atm2_c = bend_list(ibend)%c
695 36 : IF (atm2_a == atm_b .AND. atm2_b == atm_a .AND. atm2_c == atm_c) unsetme = .TRUE.
696 36 : IF (atm2_a == atm_b .AND. atm2_b == atm_a .AND. atm2_c == atm_d) unsetme = .TRUE.
697 36 : IF (atm2_a == atm_c .AND. atm2_b == atm_a .AND. atm2_c == atm_d) unsetme = .TRUE.
698 48 : IF (unsetme) bend_list(ibend)%id_type = do_ff_undef
699 : END DO
700 : END DO
701 : END DO
702 : !-----------------------------------------------------------------------------
703 : !-----------------------------------------------------------------------------
704 : ! 7. Count the number of UNSET bond kinds there are
705 : !-----------------------------------------------------------------------------
706 74487 : DO ikind = 1, SIZE(molecule_kind_set)
707 71848 : molecule_kind => molecule_kind_set(ikind)
708 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
709 : nbond=nbond, &
710 : bond_kind_set=bond_kind_set, &
711 71848 : bond_list=bond_list)
712 74487 : IF (nbond > 0) THEN
713 29732 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old BOND Count: ", &
714 606 : SIZE(bond_list), SIZE(bond_kind_set)
715 31839 : IF (iw > 0) WRITE (iw, '(2I6)') (bond_list(ibond)%a, bond_list(ibond)%b, ibond=1, SIZE(bond_list))
716 29429 : NULLIFY (bad1, bad2)
717 88287 : ALLOCATE (bad1(SIZE(bond_kind_set)))
718 96336 : bad1(:) = 0
719 96336 : DO ibond = 1, SIZE(bond_kind_set)
720 66907 : unsetme = .FALSE.
721 66907 : IF (bond_kind_set(ibond)%id_type == do_ff_undef) unsetme = .TRUE.
722 66907 : valid_kind = .FALSE.
723 348762 : DO i = 1, SIZE(bond_list)
724 347348 : IF (bond_list(i)%id_type /= do_ff_undef .AND. &
725 1414 : bond_list(i)%bond_kind%kind_number == ibond) THEN
726 : valid_kind = .TRUE.
727 : EXIT
728 : END IF
729 : END DO
730 66907 : IF (.NOT. valid_kind) unsetme = .TRUE.
731 96336 : IF (unsetme) bad1(ibond) = 1
732 : END DO
733 96336 : IF (SUM(bad1) /= 0) THEN
734 2770 : counter = SIZE(bond_kind_set) - SUM(bad1)
735 804 : CALL allocate_bond_kind_set(new_bond_kind_set, counter)
736 804 : counter = 0
737 2770 : DO ibond = 1, SIZE(bond_kind_set)
738 2770 : IF (bad1(ibond) == 0) THEN
739 546 : counter = counter + 1
740 546 : new_bond_kind_set(counter) = bond_kind_set(ibond)
741 : END IF
742 : END DO
743 804 : counter = 0
744 2412 : ALLOCATE (bad2(SIZE(bond_list)))
745 4338 : bad2(:) = 0
746 4338 : DO ibond = 1, SIZE(bond_list)
747 3534 : unsetme = .FALSE.
748 3534 : IF (bond_list(ibond)%bond_kind%id_type == do_ff_undef) unsetme = .TRUE.
749 3534 : IF (bond_list(ibond)%id_type == do_ff_undef) unsetme = .TRUE.
750 4338 : IF (unsetme) bad2(ibond) = 1
751 : END DO
752 4338 : IF (SUM(bad2) /= 0) THEN
753 4338 : counter = SIZE(bond_list) - SUM(bad2)
754 2648 : ALLOCATE (new_bond_list(counter))
755 804 : counter = 0
756 4338 : DO ibond = 1, SIZE(bond_list)
757 4338 : IF (bad2(ibond) == 0) THEN
758 926 : counter = counter + 1
759 926 : new_bond_list(counter) = bond_list(ibond)
760 926 : newkind = bond_list(ibond)%bond_kind%kind_number
761 5958 : newkind = newkind - SUM(bad1(1:newkind))
762 926 : new_bond_list(counter)%bond_kind => new_bond_kind_set(newkind)
763 : END IF
764 : END DO
765 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
766 : nbond=SIZE(new_bond_list), &
767 : bond_kind_set=new_bond_kind_set, &
768 804 : bond_list=new_bond_list)
769 1350 : DO ibond = 1, SIZE(new_bond_kind_set)
770 1350 : new_bond_kind_set(ibond)%kind_number = ibond
771 : END DO
772 : END IF
773 804 : DEALLOCATE (bad2)
774 804 : CALL deallocate_bond_kind_set(bond_kind_set)
775 804 : DEALLOCATE (bond_list)
776 808 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New BOND Count: ", &
777 8 : SIZE(new_bond_list), SIZE(new_bond_kind_set)
778 808 : IF (iw > 0) WRITE (iw, '(2I6)') (new_bond_list(ibond)%a, new_bond_list(ibond)%b, &
779 8 : ibond=1, SIZE(new_bond_list))
780 : END IF
781 29429 : DEALLOCATE (bad1)
782 : END IF
783 : END DO
784 : !-----------------------------------------------------------------------------
785 : !-----------------------------------------------------------------------------
786 : ! 8. Count the number of UNSET bend kinds there are
787 : !-----------------------------------------------------------------------------
788 74487 : DO ikind = 1, SIZE(molecule_kind_set)
789 71848 : molecule_kind => molecule_kind_set(ikind)
790 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
791 : nbend=nbend, &
792 : bend_kind_set=bend_kind_set, &
793 71848 : bend_list=bend_list)
794 74487 : IF (nbend > 0) THEN
795 29402 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old BEND Count: ", &
796 606 : SIZE(bend_list), SIZE(bend_kind_set)
797 33071 : IF (iw > 0) WRITE (iw, '(3I6)') (bend_list(ibend)%a, bend_list(ibend)%b, &
798 4275 : bend_list(ibend)%c, ibend=1, SIZE(bend_list))
799 29099 : NULLIFY (bad1, bad2)
800 87297 : ALLOCATE (bad1(SIZE(bend_kind_set)))
801 124051 : bad1(:) = 0
802 124051 : DO ibend = 1, SIZE(bend_kind_set)
803 94952 : unsetme = .FALSE.
804 94952 : IF (bend_kind_set(ibend)%id_type == do_ff_undef) unsetme = .TRUE.
805 94952 : valid_kind = .FALSE.
806 1006023 : DO i = 1, SIZE(bend_list)
807 1004717 : IF (bend_list(i)%id_type /= do_ff_undef .AND. &
808 1306 : bend_list(i)%bend_kind%kind_number == ibend) THEN
809 : valid_kind = .TRUE.
810 : EXIT
811 : END IF
812 : END DO
813 94952 : IF (.NOT. valid_kind) unsetme = .TRUE.
814 124051 : IF (unsetme) bad1(ibend) = 1
815 : END DO
816 124051 : IF (SUM(bad1) /= 0) THEN
817 2786 : counter = SIZE(bend_kind_set) - SUM(bad1)
818 606 : CALL allocate_bend_kind_set(new_bend_kind_set, counter)
819 606 : counter = 0
820 2786 : DO ibend = 1, SIZE(bend_kind_set)
821 2786 : IF (bad1(ibend) == 0) THEN
822 868 : counter = counter + 1
823 868 : new_bend_kind_set(counter) = bend_kind_set(ibend)
824 : END IF
825 : END DO
826 606 : counter = 0
827 1818 : ALLOCATE (bad2(SIZE(bend_list)))
828 4322 : bad2(:) = 0
829 4322 : DO ibend = 1, SIZE(bend_list)
830 3716 : unsetme = .FALSE.
831 3716 : IF (bend_list(ibend)%bend_kind%id_type == do_ff_undef) unsetme = .TRUE.
832 3716 : IF (bend_list(ibend)%id_type == do_ff_undef) unsetme = .TRUE.
833 4322 : IF (unsetme) bad2(ibend) = 1
834 : END DO
835 4322 : IF (SUM(bad2) /= 0) THEN
836 4322 : counter = SIZE(bend_list) - SUM(bad2)
837 2900 : ALLOCATE (new_bend_list(counter))
838 606 : counter = 0
839 4322 : DO ibend = 1, SIZE(bend_list)
840 4322 : IF (bad2(ibend) == 0) THEN
841 1610 : counter = counter + 1
842 1610 : new_bend_list(counter) = bend_list(ibend)
843 1610 : newkind = bend_list(ibend)%bend_kind%kind_number
844 18122 : newkind = newkind - SUM(bad1(1:newkind))
845 1610 : new_bend_list(counter)%bend_kind => new_bend_kind_set(newkind)
846 : END IF
847 : END DO
848 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
849 : nbend=SIZE(new_bend_list), &
850 : bend_kind_set=new_bend_kind_set, &
851 606 : bend_list=new_bend_list)
852 1474 : DO ibend = 1, SIZE(new_bend_kind_set)
853 1474 : new_bend_kind_set(ibend)%kind_number = ibend
854 : END DO
855 : END IF
856 606 : DEALLOCATE (bad2)
857 606 : CALL deallocate_bend_kind_set(bend_kind_set)
858 606 : DEALLOCATE (bend_list)
859 610 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New BEND Count: ", &
860 8 : SIZE(new_bend_list), SIZE(new_bend_kind_set)
861 606 : IF (iw > 0) WRITE (iw, '(3I6)') (new_bend_list(ibend)%a, new_bend_list(ibend)%b, &
862 4 : new_bend_list(ibend)%c, ibend=1, SIZE(new_bend_list))
863 : END IF
864 29099 : DEALLOCATE (bad1)
865 : END IF
866 : END DO
867 :
868 : !-----------------------------------------------------------------------------
869 : !-----------------------------------------------------------------------------
870 : ! 9. Count the number of UNSET Urey-Bradley kinds there are
871 : !-----------------------------------------------------------------------------
872 74487 : DO ikind = 1, SIZE(molecule_kind_set)
873 71848 : molecule_kind => molecule_kind_set(ikind)
874 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
875 : nub=nub, &
876 : ub_kind_set=ub_kind_set, &
877 71848 : ub_list=ub_list)
878 74487 : IF (nub > 0) THEN
879 29388 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old UB Count: ", &
880 606 : SIZE(ub_list), SIZE(ub_kind_set)
881 33057 : IF (iw > 0) WRITE (iw, '(3I6)') (ub_list(iub)%a, ub_list(iub)%b, &
882 4275 : ub_list(iub)%c, iub=1, SIZE(ub_list))
883 29085 : NULLIFY (bad1, bad2)
884 87255 : ALLOCATE (bad1(SIZE(ub_kind_set)))
885 123879 : bad1(:) = 0
886 123879 : DO iub = 1, SIZE(ub_kind_set)
887 94794 : unsetme = .FALSE.
888 94794 : IF (ub_kind_set(iub)%id_type == do_ff_undef) unsetme = .TRUE.
889 94794 : valid_kind = .FALSE.
890 1927808 : DO i = 1, SIZE(ub_list)
891 1841384 : IF (ub_list(i)%id_type /= do_ff_undef .AND. &
892 86424 : ub_list(i)%ub_kind%kind_number == iub) THEN
893 : valid_kind = .TRUE.
894 : EXIT
895 : END IF
896 : END DO
897 94794 : IF (.NOT. valid_kind) unsetme = .TRUE.
898 123879 : IF (unsetme) bad1(iub) = 1
899 : END DO
900 123879 : IF (SUM(bad1) /= 0) THEN
901 123843 : counter = SIZE(ub_kind_set) - SUM(bad1)
902 29077 : CALL allocate_ub_kind_set(new_ub_kind_set, counter)
903 29077 : counter = 0
904 123843 : DO iub = 1, SIZE(ub_kind_set)
905 123843 : IF (bad1(iub) == 0) THEN
906 8342 : counter = counter + 1
907 8342 : new_ub_kind_set(counter) = ub_kind_set(iub)
908 : END IF
909 : END DO
910 29077 : counter = 0
911 87231 : ALLOCATE (bad2(SIZE(ub_list)))
912 167905 : bad2(:) = 0
913 167905 : DO iub = 1, SIZE(ub_list)
914 138828 : unsetme = .FALSE.
915 138828 : IF (ub_list(iub)%ub_kind%id_type == do_ff_undef) unsetme = .TRUE.
916 138828 : IF (ub_list(iub)%id_type == do_ff_undef) unsetme = .TRUE.
917 167905 : IF (unsetme) bad2(iub) = 1
918 : END DO
919 167905 : IF (SUM(bad2) /= 0) THEN
920 167905 : counter = SIZE(ub_list) - SUM(bad2)
921 79670 : ALLOCATE (new_ub_list(counter))
922 29077 : counter = 0
923 167905 : DO iub = 1, SIZE(ub_list)
924 167905 : IF (bad2(iub) == 0) THEN
925 19972 : counter = counter + 1
926 19972 : new_ub_list(counter) = ub_list(iub)
927 19972 : newkind = ub_list(iub)%ub_kind%kind_number
928 256656 : newkind = newkind - SUM(bad1(1:newkind))
929 19972 : new_ub_list(counter)%ub_kind => new_ub_kind_set(newkind)
930 : END IF
931 : END DO
932 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
933 : nub=SIZE(new_ub_list), &
934 : ub_kind_set=new_ub_kind_set, &
935 29077 : ub_list=new_ub_list)
936 37419 : DO iub = 1, SIZE(new_ub_kind_set)
937 37419 : new_ub_kind_set(iub)%kind_number = iub
938 : END DO
939 : END IF
940 29077 : DEALLOCATE (bad2)
941 29077 : CALL ub_kind_dealloc_ref(ub_kind_set)
942 29077 : DEALLOCATE (ub_list)
943 29380 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New UB Count: ", &
944 606 : SIZE(new_ub_list), SIZE(new_ub_kind_set)
945 29215 : IF (iw > 0) WRITE (iw, '(3I6)') (new_ub_list(iub)%a, new_ub_list(iub)%b, &
946 441 : new_ub_list(iub)%c, iub=1, SIZE(new_ub_list))
947 : END IF
948 29085 : DEALLOCATE (bad1)
949 : END IF
950 : END DO
951 :
952 : !-----------------------------------------------------------------------------
953 : !-----------------------------------------------------------------------------
954 : ! 10. Count the number of UNSET torsion kinds there are
955 : !-----------------------------------------------------------------------------
956 74487 : DO ikind = 1, SIZE(molecule_kind_set)
957 71848 : molecule_kind => molecule_kind_set(ikind)
958 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
959 : ntorsion=ntorsion, &
960 : torsion_kind_set=torsion_kind_set, &
961 71848 : torsion_list=torsion_list)
962 74487 : IF (ntorsion > 0) THEN
963 5801 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old TORSION Count: ", &
964 538 : SIZE(torsion_list), SIZE(torsion_kind_set)
965 10250 : IF (iw > 0) WRITE (iw, '(4I6)') (torsion_list(itorsion)%a, torsion_list(itorsion)%b, &
966 4987 : torsion_list(itorsion)%c, torsion_list(itorsion)%d, itorsion=1, SIZE(torsion_list))
967 5532 : NULLIFY (bad1, bad2)
968 16596 : ALLOCATE (bad1(SIZE(torsion_kind_set)))
969 98203 : bad1(:) = 0
970 98203 : DO itorsion = 1, SIZE(torsion_kind_set)
971 92671 : unsetme = .FALSE.
972 92671 : IF (torsion_kind_set(itorsion)%id_type == do_ff_undef) unsetme = .TRUE.
973 92671 : valid_kind = .FALSE.
974 2290223 : DO i = 1, SIZE(torsion_list)
975 2275981 : IF (torsion_list(i)%id_type /= do_ff_undef .AND. &
976 14242 : torsion_list(i)%torsion_kind%kind_number == itorsion) THEN
977 : valid_kind = .TRUE.
978 : EXIT
979 : END IF
980 : END DO
981 92671 : IF (.NOT. valid_kind) unsetme = .TRUE.
982 98203 : IF (unsetme) bad1(itorsion) = 1
983 : END DO
984 98203 : IF (SUM(bad1) /= 0) THEN
985 17930 : counter = SIZE(torsion_kind_set) - SUM(bad1)
986 1018 : CALL allocate_torsion_kind_set(new_torsion_kind_set, counter)
987 1018 : counter = 0
988 17930 : DO itorsion = 1, SIZE(torsion_kind_set)
989 17930 : IF (bad1(itorsion) == 0) THEN
990 2670 : counter = counter + 1
991 2670 : new_torsion_kind_set(counter) = torsion_kind_set(itorsion)
992 2670 : i = SIZE(torsion_kind_set(itorsion)%m)
993 2670 : j = SIZE(torsion_kind_set(itorsion)%k)
994 2670 : k = SIZE(torsion_kind_set(itorsion)%phi0)
995 8010 : ALLOCATE (new_torsion_kind_set(counter)%m(i))
996 8010 : ALLOCATE (new_torsion_kind_set(counter)%k(i))
997 5340 : ALLOCATE (new_torsion_kind_set(counter)%phi0(i))
998 11704 : new_torsion_kind_set(counter)%m = torsion_kind_set(itorsion)%m
999 11704 : new_torsion_kind_set(counter)%k = torsion_kind_set(itorsion)%k
1000 11704 : new_torsion_kind_set(counter)%phi0 = torsion_kind_set(itorsion)%phi0
1001 : END IF
1002 : END DO
1003 1018 : counter = 0
1004 3054 : ALLOCATE (bad2(SIZE(torsion_list)))
1005 42940 : bad2(:) = 0
1006 42940 : DO itorsion = 1, SIZE(torsion_list)
1007 41922 : unsetme = .FALSE.
1008 41922 : IF (torsion_list(itorsion)%torsion_kind%id_type == do_ff_undef) unsetme = .TRUE.
1009 41922 : IF (torsion_list(itorsion)%id_type == do_ff_undef) unsetme = .TRUE.
1010 42940 : IF (unsetme) bad2(itorsion) = 1
1011 : END DO
1012 42940 : IF (SUM(bad2) /= 0) THEN
1013 42940 : counter = SIZE(torsion_list) - SUM(bad2)
1014 8634 : ALLOCATE (new_torsion_list(counter))
1015 1018 : counter = 0
1016 42940 : DO itorsion = 1, SIZE(torsion_list)
1017 42940 : IF (bad2(itorsion) == 0) THEN
1018 6484 : counter = counter + 1
1019 6484 : new_torsion_list(counter) = torsion_list(itorsion)
1020 6484 : newkind = torsion_list(itorsion)%torsion_kind%kind_number
1021 128566 : newkind = newkind - SUM(bad1(1:newkind))
1022 6484 : new_torsion_list(counter)%torsion_kind => new_torsion_kind_set(newkind)
1023 : END IF
1024 : END DO
1025 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
1026 : ntorsion=SIZE(new_torsion_list), &
1027 : torsion_kind_set=new_torsion_kind_set, &
1028 1018 : torsion_list=new_torsion_list)
1029 3688 : DO itorsion = 1, SIZE(new_torsion_kind_set)
1030 3688 : new_torsion_kind_set(itorsion)%kind_number = itorsion
1031 : END DO
1032 : END IF
1033 1018 : DEALLOCATE (bad2)
1034 17930 : DO itorsion = 1, SIZE(torsion_kind_set)
1035 17930 : CALL torsion_kind_dealloc_ref(torsion_kind_set(itorsion))
1036 : END DO
1037 1018 : DEALLOCATE (torsion_kind_set)
1038 1018 : DEALLOCATE (torsion_list)
1039 1152 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New TORSION Count: ", &
1040 268 : SIZE(new_torsion_list), SIZE(new_torsion_kind_set)
1041 1165 : IF (iw > 0) WRITE (iw, '(4I6)') (new_torsion_list(itorsion)%a, new_torsion_list(itorsion)%b, &
1042 428 : new_torsion_list(itorsion)%c, new_torsion_list(itorsion)%d, itorsion=1, &
1043 415 : SIZE(new_torsion_list))
1044 : END IF
1045 5532 : DEALLOCATE (bad1)
1046 : END IF
1047 : END DO
1048 :
1049 : !-----------------------------------------------------------------------------
1050 : !-----------------------------------------------------------------------------
1051 : ! 11. Count the number of UNSET improper kinds there are
1052 : !-----------------------------------------------------------------------------
1053 74487 : DO ikind = 1, SIZE(molecule_kind_set)
1054 71848 : molecule_kind => molecule_kind_set(ikind)
1055 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
1056 : nimpr=nimpr, &
1057 : impr_kind_set=impr_kind_set, &
1058 71848 : impr_list=impr_list)
1059 74487 : IF (nimpr > 0) THEN
1060 1695 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old IMPROPER Count: ", &
1061 46 : SIZE(impr_list), SIZE(impr_kind_set)
1062 1672 : NULLIFY (bad1, bad2)
1063 5016 : ALLOCATE (bad1(SIZE(impr_kind_set)))
1064 6380 : bad1(:) = 0
1065 6380 : DO iimpr = 1, SIZE(impr_kind_set)
1066 4708 : unsetme = .FALSE.
1067 4708 : IF (impr_kind_set(iimpr)%id_type == do_ff_undef) unsetme = .TRUE.
1068 4708 : valid_kind = .FALSE.
1069 28820 : DO i = 1, SIZE(impr_list)
1070 27322 : IF (impr_list(i)%id_type /= do_ff_undef .AND. &
1071 1498 : impr_list(i)%impr_kind%kind_number == iimpr) THEN
1072 : valid_kind = .TRUE.
1073 : EXIT
1074 : END IF
1075 : END DO
1076 4708 : IF (.NOT. valid_kind) unsetme = .TRUE.
1077 6380 : IF (unsetme) bad1(iimpr) = 1
1078 : END DO
1079 6380 : IF (SUM(bad1) /= 0) THEN
1080 2012 : counter = SIZE(impr_kind_set) - SUM(bad1)
1081 390 : CALL allocate_impr_kind_set(new_impr_kind_set, counter)
1082 390 : counter = 0
1083 2012 : DO iimpr = 1, SIZE(impr_kind_set)
1084 2012 : IF (bad1(iimpr) == 0) THEN
1085 124 : counter = counter + 1
1086 124 : new_impr_kind_set(counter) = impr_kind_set(iimpr)
1087 : END IF
1088 : END DO
1089 390 : counter = 0
1090 1170 : ALLOCATE (bad2(SIZE(impr_list)))
1091 2420 : bad2(:) = 0
1092 2420 : DO iimpr = 1, SIZE(impr_list)
1093 2030 : unsetme = .FALSE.
1094 2030 : IF (impr_list(iimpr)%impr_kind%id_type == do_ff_undef) unsetme = .TRUE.
1095 2030 : IF (impr_list(iimpr)%id_type == do_ff_undef) unsetme = .TRUE.
1096 2420 : IF (unsetme) bad2(iimpr) = 1
1097 : END DO
1098 2420 : IF (SUM(bad2) /= 0) THEN
1099 2420 : counter = SIZE(impr_list) - SUM(bad2)
1100 1008 : ALLOCATE (new_impr_list(counter))
1101 390 : counter = 0
1102 2420 : DO iimpr = 1, SIZE(impr_list)
1103 2420 : IF (bad2(iimpr) == 0) THEN
1104 124 : counter = counter + 1
1105 124 : new_impr_list(counter) = impr_list(iimpr)
1106 124 : newkind = impr_list(iimpr)%impr_kind%kind_number
1107 324 : newkind = newkind - SUM(bad1(1:newkind))
1108 124 : new_impr_list(counter)%impr_kind => new_impr_kind_set(newkind)
1109 : END IF
1110 : END DO
1111 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
1112 : nimpr=SIZE(new_impr_list), &
1113 : impr_kind_set=new_impr_kind_set, &
1114 390 : impr_list=new_impr_list)
1115 514 : DO iimpr = 1, SIZE(new_impr_kind_set)
1116 514 : new_impr_kind_set(iimpr)%kind_number = iimpr
1117 : END DO
1118 : END IF
1119 390 : DEALLOCATE (bad2)
1120 2012 : DO iimpr = 1, SIZE(impr_kind_set)
1121 2012 : CALL impr_kind_dealloc_ref() !This Subroutine doesn't deallocate anything, maybe needs to be implemented
1122 : END DO
1123 390 : DEALLOCATE (impr_kind_set)
1124 390 : DEALLOCATE (impr_list)
1125 401 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New IMPROPER Count: ", &
1126 22 : SIZE(new_impr_list), SIZE(new_impr_kind_set)
1127 : END IF
1128 1672 : DEALLOCATE (bad1)
1129 : END IF
1130 : END DO
1131 :
1132 : !-----------------------------------------------------------------------------
1133 : !-----------------------------------------------------------------------------
1134 : ! 11. Count the number of UNSET opbends kinds there are
1135 : !-----------------------------------------------------------------------------
1136 74487 : DO ikind = 1, SIZE(molecule_kind_set)
1137 71848 : molecule_kind => molecule_kind_set(ikind)
1138 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
1139 : nopbend=nopbend, &
1140 : opbend_kind_set=opbend_kind_set, &
1141 71848 : opbend_list=opbend_list)
1142 74487 : IF (nopbend > 0) THEN
1143 1695 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old OPBEND Count: ", &
1144 46 : SIZE(opbend_list), SIZE(opbend_kind_set)
1145 1672 : NULLIFY (bad1, bad2)
1146 5016 : ALLOCATE (bad1(SIZE(opbend_kind_set)))
1147 6380 : bad1(:) = 0
1148 6380 : DO iopbend = 1, SIZE(opbend_kind_set)
1149 4708 : unsetme = .FALSE.
1150 4708 : IF (opbend_kind_set(iopbend)%id_type == do_ff_undef) unsetme = .TRUE.
1151 4708 : valid_kind = .FALSE.
1152 35848 : DO i = 1, SIZE(opbend_list)
1153 31142 : IF (opbend_list(i)%id_type /= do_ff_undef .AND. &
1154 4706 : opbend_list(i)%opbend_kind%kind_number == iopbend) THEN
1155 : valid_kind = .TRUE.
1156 : EXIT
1157 : END IF
1158 : END DO
1159 4708 : IF (.NOT. valid_kind) unsetme = .TRUE.
1160 6380 : IF (unsetme) bad1(iopbend) = 1
1161 : END DO
1162 6380 : IF (SUM(bad1) /= 0) THEN
1163 6376 : counter = SIZE(opbend_kind_set) - SUM(bad1)
1164 1670 : CALL allocate_opbend_kind_set(new_opbend_kind_set, counter)
1165 1670 : counter = 0
1166 6376 : DO iopbend = 1, SIZE(opbend_kind_set)
1167 6376 : IF (bad1(iopbend) == 0) THEN
1168 0 : counter = counter + 1
1169 0 : new_opbend_kind_set(counter) = opbend_kind_set(iopbend)
1170 : END IF
1171 : END DO
1172 1670 : counter = 0
1173 5010 : ALLOCATE (bad2(SIZE(opbend_list)))
1174 6980 : bad2(:) = 0
1175 6980 : DO iopbend = 1, SIZE(opbend_list)
1176 5310 : unsetme = .FALSE.
1177 5310 : IF (opbend_list(iopbend)%opbend_kind%id_type == do_ff_undef) unsetme = .TRUE.
1178 5310 : IF (opbend_list(iopbend)%id_type == do_ff_undef) unsetme = .TRUE.
1179 6980 : IF (unsetme) bad2(iopbend) = 1
1180 : END DO
1181 6980 : IF (SUM(bad2) /= 0) THEN
1182 6980 : counter = SIZE(opbend_list) - SUM(bad2)
1183 3340 : ALLOCATE (new_opbend_list(counter))
1184 1670 : counter = 0
1185 6980 : DO iopbend = 1, SIZE(opbend_list)
1186 6980 : IF (bad2(iopbend) == 0) THEN
1187 0 : counter = counter + 1
1188 0 : new_opbend_list(counter) = opbend_list(iopbend)
1189 0 : newkind = opbend_list(iopbend)%opbend_kind%kind_number
1190 0 : newkind = newkind - SUM(bad1(1:newkind))
1191 0 : new_opbend_list(counter)%opbend_kind => new_opbend_kind_set(newkind)
1192 : END IF
1193 : END DO
1194 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
1195 : nopbend=SIZE(new_opbend_list), &
1196 : opbend_kind_set=new_opbend_kind_set, &
1197 1670 : opbend_list=new_opbend_list)
1198 1670 : DO iopbend = 1, SIZE(new_opbend_kind_set)
1199 1670 : new_opbend_kind_set(iopbend)%kind_number = iopbend
1200 : END DO
1201 : END IF
1202 1670 : DEALLOCATE (bad2)
1203 1670 : DEALLOCATE (opbend_kind_set)
1204 1670 : DEALLOCATE (opbend_list)
1205 1692 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New OPBEND Count: ", &
1206 44 : SIZE(new_opbend_list), SIZE(new_opbend_kind_set)
1207 : END IF
1208 1672 : DEALLOCATE (bad1)
1209 : END IF
1210 : END DO
1211 :
1212 : !---------------------------------------------------------------------------
1213 : !---------------------------------------------------------------------------
1214 : ! 12. Count the number of UNSET NONBOND14 kinds there are
1215 : !- NEED TO REMOVE EXTRAS HERE - IKUO
1216 : !---------------------------------------------------------------------------
1217 : CALL cp_print_key_finished_output(iw, logger, mm_section, &
1218 2639 : "PRINT%FF_INFO")
1219 2639 : CALL timestop(handle)
1220 2639 : END SUBROUTINE clean_intra_force_kind
1221 :
1222 : ! **************************************************************************************************
1223 : !> \brief Reads from the input structure all information for generic functions
1224 : !> \param gen_section ...
1225 : !> \param func_name ...
1226 : !> \param xfunction ...
1227 : !> \param parameters ...
1228 : !> \param values ...
1229 : !> \param var_values ...
1230 : !> \param size_variables ...
1231 : !> \param i_rep_sec ...
1232 : !> \param input_variables ...
1233 : ! **************************************************************************************************
1234 4068 : SUBROUTINE get_generic_info(gen_section, func_name, xfunction, parameters, values, &
1235 4068 : var_values, size_variables, i_rep_sec, input_variables)
1236 : TYPE(section_vals_type), POINTER :: gen_section
1237 : CHARACTER(LEN=*), INTENT(IN) :: func_name
1238 : CHARACTER(LEN=default_path_length), INTENT(OUT) :: xfunction
1239 : CHARACTER(LEN=default_string_length), &
1240 : DIMENSION(:), POINTER :: parameters
1241 : REAL(KIND=dp), DIMENSION(:), POINTER :: values
1242 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: var_values
1243 : INTEGER, INTENT(IN), OPTIONAL :: size_variables, i_rep_sec
1244 : CHARACTER(LEN=*), DIMENSION(:), OPTIONAL :: input_variables
1245 :
1246 : CHARACTER(LEN=default_string_length), &
1247 4068 : DIMENSION(:), POINTER :: my_par, my_par_tmp, my_units, &
1248 4068 : my_units_tmp, my_var
1249 : INTEGER :: i, ind, irep, isize, j, mydim, n_par, &
1250 : n_units, n_val, nblank
1251 : LOGICAL :: check
1252 4068 : REAL(KIND=dp), DIMENSION(:), POINTER :: my_val, my_val_tmp
1253 :
1254 4068 : NULLIFY (my_var, my_par, my_val, my_par_tmp, my_val_tmp)
1255 4068 : NULLIFY (my_units)
1256 4068 : NULLIFY (my_units_tmp)
1257 4068 : IF (ASSOCIATED(parameters)) THEN
1258 322 : DEALLOCATE (parameters)
1259 : END IF
1260 4068 : IF (ASSOCIATED(values)) THEN
1261 322 : DEALLOCATE (values)
1262 : END IF
1263 4068 : irep = 1
1264 4068 : IF (PRESENT(i_rep_sec)) irep = i_rep_sec
1265 4068 : mydim = 0
1266 4068 : CALL section_vals_val_get(gen_section, TRIM(func_name), i_rep_section=irep, c_val=xfunction)
1267 4068 : CALL compress(xfunction, full=.TRUE.)
1268 4068 : IF (PRESENT(input_variables)) THEN
1269 1398 : ALLOCATE (my_var(SIZE(input_variables)))
1270 1864 : my_var = input_variables
1271 : ELSE
1272 3602 : CALL section_vals_val_get(gen_section, "VARIABLES", i_rep_section=irep, c_vals=my_var)
1273 : END IF
1274 4068 : IF (ASSOCIATED(my_var)) THEN
1275 4068 : mydim = SIZE(my_var)
1276 : END IF
1277 4068 : IF (PRESENT(size_variables)) THEN
1278 3218 : CPASSERT(mydim == size_variables)
1279 : END IF
1280 : ! Check for presence of Parameters
1281 4068 : CALL section_vals_val_get(gen_section, "PARAMETERS", i_rep_section=irep, n_rep_val=n_par)
1282 4068 : CALL section_vals_val_get(gen_section, "VALUES", i_rep_section=irep, n_rep_val=n_val)
1283 4068 : check = (n_par > 0) .EQV. (n_val > 0)
1284 4068 : CPASSERT(check)
1285 4068 : CALL section_vals_val_get(gen_section, "UNITS", i_rep_section=irep, n_rep_val=n_units)
1286 4068 : IF (n_par > 0) THEN
1287 : ! Parameters
1288 3454 : ALLOCATE (my_par(0))
1289 3454 : ALLOCATE (my_val(0))
1290 6908 : DO i = 1, n_par
1291 3454 : isize = SIZE(my_par)
1292 3454 : CALL section_vals_val_get(gen_section, "PARAMETERS", i_rep_section=irep, i_rep_val=i, c_vals=my_par_tmp)
1293 10268 : nblank = COUNT(my_par_tmp == "")
1294 3454 : CALL reallocate(my_par, 1, isize + SIZE(my_par_tmp) - nblank)
1295 3454 : ind = 0
1296 13722 : DO j = 1, SIZE(my_par_tmp)
1297 6814 : IF (my_par_tmp(j) == "") CYCLE
1298 6812 : ind = ind + 1
1299 10268 : my_par(isize + ind) = my_par_tmp(j)
1300 : END DO
1301 : END DO
1302 6910 : DO i = 1, n_val
1303 3456 : isize = SIZE(my_val)
1304 3456 : CALL section_vals_val_get(gen_section, "VALUES", i_rep_section=irep, i_rep_val=i, r_vals=my_val_tmp)
1305 3456 : CALL reallocate(my_val, 1, isize + SIZE(my_val_tmp))
1306 20534 : my_val(isize + 1:isize + SIZE(my_val_tmp)) = my_val_tmp
1307 : END DO
1308 3454 : CPASSERT(SIZE(my_par) == SIZE(my_val))
1309 : ! Optionally read the units for each parameter value
1310 3454 : ALLOCATE (my_units(0))
1311 3454 : IF (n_units > 0) THEN
1312 6 : DO i = 1, n_units
1313 4 : isize = SIZE(my_units)
1314 4 : CALL section_vals_val_get(gen_section, "UNITS", i_rep_section=irep, i_rep_val=i, c_vals=my_units_tmp)
1315 10 : nblank = COUNT(my_units_tmp == "")
1316 4 : CALL reallocate(my_units, 1, isize + SIZE(my_units_tmp) - nblank)
1317 4 : ind = 0
1318 12 : DO j = 1, SIZE(my_units_tmp)
1319 6 : IF (my_units_tmp(j) == "") CYCLE
1320 6 : ind = ind + 1
1321 10 : my_units(isize + ind) = my_units_tmp(j)
1322 : END DO
1323 : END DO
1324 2 : CPASSERT(SIZE(my_units) == SIZE(my_val))
1325 : END IF
1326 3454 : mydim = mydim + SIZE(my_val)
1327 3454 : IF (SIZE(my_val) == 0) THEN
1328 2 : DEALLOCATE (my_par)
1329 2 : DEALLOCATE (my_val)
1330 2 : DEALLOCATE (my_units)
1331 : END IF
1332 : END IF
1333 : ! Handle trivial case of a null function defined
1334 12204 : ALLOCATE (parameters(mydim))
1335 12204 : ALLOCATE (values(mydim))
1336 4068 : IF (mydim > 0) THEN
1337 14972 : parameters(1:SIZE(my_var)) = my_var
1338 9520 : values(1:SIZE(my_var)) = 0.0_dp
1339 4068 : IF (PRESENT(var_values)) THEN
1340 384 : CPASSERT(SIZE(var_values) == SIZE(my_var))
1341 2056 : values(1:SIZE(my_var)) = var_values
1342 : END IF
1343 4068 : IF (ASSOCIATED(my_val)) THEN
1344 10264 : DO i = 1, SIZE(my_val)
1345 6812 : parameters(SIZE(my_var) + i) = my_par(i)
1346 10264 : IF (n_units > 0) THEN
1347 6 : values(SIZE(my_var) + i) = cp_unit_to_cp2k(my_val(i), TRIM(ADJUSTL(my_units(i))))
1348 : ELSE
1349 6806 : values(SIZE(my_var) + i) = my_val(i)
1350 : END IF
1351 : END DO
1352 : END IF
1353 : END IF
1354 4068 : IF (ASSOCIATED(my_par)) THEN
1355 3452 : DEALLOCATE (my_par)
1356 : END IF
1357 4068 : IF (ASSOCIATED(my_val)) THEN
1358 3452 : DEALLOCATE (my_val)
1359 : END IF
1360 4068 : IF (ASSOCIATED(my_units)) THEN
1361 3452 : DEALLOCATE (my_units)
1362 : END IF
1363 4068 : IF (PRESENT(input_variables)) THEN
1364 466 : DEALLOCATE (my_var)
1365 : END IF
1366 12204 : END SUBROUTINE get_generic_info
1367 :
1368 : END MODULE force_fields_util
|
