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