Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \par History
10 : !> Subroutine input_torsions changed (DG) 05-Dec-2000
11 : !> Output formats changed (DG) 05-Dec-2000
12 : !> JGH (26-01-2002) : force field parameters stored in tables, not in
13 : !> matrices. Input changed to have parameters labeled by the position
14 : !> and not atom pairs (triples etc)
15 : !> Teo (11.2005) : Moved all information on force field pair_potential to
16 : !> a much lighter memory structure
17 : !> \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 ewald_environment_types, ONLY: ewald_env_get,&
31 : ewald_environment_type
32 : USE fist_nonbond_env_types, ONLY: fist_nonbond_env_create,&
33 : fist_nonbond_env_set,&
34 : fist_nonbond_env_type
35 : USE force_field_kind_types, ONLY: &
36 : allocate_bend_kind_set, allocate_bond_kind_set, allocate_impr_kind_set, &
37 : allocate_opbend_kind_set, allocate_torsion_kind_set, allocate_ub_kind_set, bend_kind_type, &
38 : bond_kind_type, deallocate_bend_kind_set, deallocate_bond_kind_set, do_ff_undef, &
39 : impr_kind_dealloc_ref, impr_kind_type, opbend_kind_type, torsion_kind_dealloc_ref, &
40 : torsion_kind_type, ub_kind_dealloc_ref, ub_kind_type
41 : USE force_field_types, ONLY: amber_info_type,&
42 : charmm_info_type,&
43 : force_field_type,&
44 : gromos_info_type,&
45 : input_info_type
46 : USE force_fields_all, ONLY: &
47 : force_field_pack_bend, force_field_pack_bond, force_field_pack_charge, &
48 : force_field_pack_charges, force_field_pack_damp, force_field_pack_eicut, &
49 : force_field_pack_impr, force_field_pack_nonbond, force_field_pack_nonbond14, &
50 : force_field_pack_opbend, force_field_pack_pol, force_field_pack_radius, &
51 : force_field_pack_shell, force_field_pack_splines, force_field_pack_tors, &
52 : force_field_pack_ub, force_field_unique_bend, force_field_unique_bond, &
53 : force_field_unique_impr, force_field_unique_opbend, force_field_unique_tors, &
54 : force_field_unique_ub
55 : USE input_section_types, ONLY: section_vals_get,&
56 : section_vals_get_subs_vals,&
57 : section_vals_type,&
58 : section_vals_val_get
59 : USE kinds, ONLY: default_string_length,&
60 : dp
61 : USE molecule_kind_types, ONLY: &
62 : atom_type, bend_type, bond_type, colvar_constraint_type, g3x3_constraint_type, &
63 : g4x6_constraint_type, get_molecule_kind, impr_type, molecule_kind_type, opbend_type, &
64 : set_molecule_kind, torsion_type, ub_type
65 : USE molecule_types, ONLY: get_molecule,&
66 : molecule_type
67 : USE pair_potential_types, ONLY: pair_potential_pp_type
68 : USE particle_types, ONLY: particle_type
69 : USE qmmm_types_low, ONLY: qmmm_env_mm_type
70 : USE shell_potential_types, ONLY: shell_kind_type
71 : #include "./base/base_uses.f90"
72 :
73 : IMPLICIT NONE
74 :
75 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_fields_util'
76 :
77 : PRIVATE
78 : PUBLIC :: force_field_pack, &
79 : force_field_qeff_output, &
80 : clean_intra_force_kind
81 :
82 : CONTAINS
83 :
84 : ! **************************************************************************************************
85 : !> \brief Pack in all the information needed for the force_field
86 : !> \param particle_set ...
87 : !> \param atomic_kind_set ...
88 : !> \param molecule_kind_set ...
89 : !> \param molecule_set ...
90 : !> \param ewald_env ...
91 : !> \param fist_nonbond_env ...
92 : !> \param ff_type ...
93 : !> \param root_section ...
94 : !> \param qmmm ...
95 : !> \param qmmm_env ...
96 : !> \param mm_section ...
97 : !> \param subsys_section ...
98 : !> \param shell_particle_set ...
99 : !> \param core_particle_set ...
100 : !> \param cell ...
101 : ! **************************************************************************************************
102 10796 : SUBROUTINE force_field_pack(particle_set, atomic_kind_set, molecule_kind_set, &
103 : molecule_set, ewald_env, fist_nonbond_env, ff_type, root_section, qmmm, &
104 : qmmm_env, mm_section, subsys_section, shell_particle_set, core_particle_set, &
105 : cell)
106 :
107 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
108 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
109 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
110 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
111 : TYPE(ewald_environment_type), POINTER :: ewald_env
112 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
113 : TYPE(force_field_type), INTENT(INOUT) :: ff_type
114 : TYPE(section_vals_type), POINTER :: root_section
115 : LOGICAL, INTENT(IN), OPTIONAL :: qmmm
116 : TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env
117 : TYPE(section_vals_type), POINTER :: mm_section, subsys_section
118 : TYPE(particle_type), DIMENSION(:), POINTER :: shell_particle_set, core_particle_set
119 : TYPE(cell_type), POINTER :: cell
120 :
121 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack'
122 :
123 : CHARACTER(LEN=default_string_length), &
124 2699 : DIMENSION(:), POINTER :: Ainfo
125 : INTEGER :: handle, iw, iw2, iw3, iw4, output_unit
126 : LOGICAL :: do_zbl, explicit, fatal, ignore_fatal, &
127 : my_qmmm
128 : REAL(KIND=dp) :: ewald_rcut, verlet_skin
129 2699 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
130 : TYPE(amber_info_type), POINTER :: amb_info
131 : TYPE(charmm_info_type), POINTER :: chm_info
132 : TYPE(cp_logger_type), POINTER :: logger
133 : TYPE(gromos_info_type), POINTER :: gro_info
134 : TYPE(input_info_type), POINTER :: inp_info
135 : TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond, potparm_nonbond14
136 : TYPE(section_vals_type), POINTER :: charges_section
137 :
138 2699 : CALL timeset(routineN, handle)
139 2699 : fatal = .FALSE.
140 2699 : ignore_fatal = ff_type%ignore_missing_critical
141 2699 : NULLIFY (logger, Ainfo, charges_section, charges)
142 2699 : logger => cp_get_default_logger()
143 : ! Error unit
144 2699 : output_unit = cp_logger_get_default_io_unit(logger)
145 :
146 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
147 2699 : extension=".mmLog")
148 : iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO/SPLINE_INFO", &
149 2699 : extension=".mmLog")
150 : iw3 = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO/SPLINE_DATA", &
151 2699 : extension=".mmLog")
152 : iw4 = cp_print_key_unit_nr(logger, mm_section, "PRINT%PROGRAM_RUN_INFO", &
153 2699 : extension=".mmLog")
154 2699 : NULLIFY (potparm_nonbond14, potparm_nonbond)
155 2699 : my_qmmm = .FALSE.
156 2699 : IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
157 2699 : inp_info => ff_type%inp_info
158 2699 : chm_info => ff_type%chm_info
159 2699 : gro_info => ff_type%gro_info
160 2699 : amb_info => ff_type%amb_info
161 : !-----------------------------------------------------------------------------
162 : ! 1. Determine the number of unique bond kind and allocate bond_kind_set
163 : !-----------------------------------------------------------------------------
164 : CALL force_field_unique_bond(particle_set, molecule_kind_set, molecule_set, &
165 2699 : ff_type, iw)
166 : !-----------------------------------------------------------------------------
167 : ! 2. Determine the number of unique bend kind and allocate bend_kind_set
168 : !-----------------------------------------------------------------------------
169 : CALL force_field_unique_bend(particle_set, molecule_kind_set, molecule_set, &
170 2699 : ff_type, iw)
171 : !-----------------------------------------------------------------------------
172 : ! 3. Determine the number of unique Urey-Bradley kind and allocate ub_kind_set
173 : !-----------------------------------------------------------------------------
174 2699 : CALL force_field_unique_ub(particle_set, molecule_kind_set, molecule_set, iw)
175 : !-----------------------------------------------------------------------------
176 : ! 4. Determine the number of unique torsion kind and allocate torsion_kind_set
177 : !-----------------------------------------------------------------------------
178 : CALL force_field_unique_tors(particle_set, molecule_kind_set, molecule_set, &
179 2699 : ff_type, iw)
180 : !-----------------------------------------------------------------------------
181 : ! 5. Determine the number of unique impr kind and allocate impr_kind_set
182 : !-----------------------------------------------------------------------------
183 : CALL force_field_unique_impr(particle_set, molecule_kind_set, molecule_set, &
184 2699 : ff_type, iw)
185 : !-----------------------------------------------------------------------------
186 : ! 6. Determine the number of unique opbend kind and allocate opbend_kind_set
187 : !-----------------------------------------------------------------------------
188 : CALL force_field_unique_opbend(particle_set, molecule_kind_set, molecule_set, &
189 2699 : ff_type, iw)
190 : !-----------------------------------------------------------------------------
191 : ! 7. Bonds
192 : !-----------------------------------------------------------------------------
193 : CALL force_field_pack_bond(particle_set, molecule_kind_set, molecule_set, &
194 : fatal, Ainfo, chm_info, inp_info, gro_info, &
195 2699 : amb_info, iw)
196 : !-----------------------------------------------------------------------------
197 : ! 8. Bends
198 : !-----------------------------------------------------------------------------
199 : CALL force_field_pack_bend(particle_set, molecule_kind_set, molecule_set, &
200 : fatal, Ainfo, chm_info, inp_info, gro_info, &
201 2699 : amb_info, iw)
202 : ! Give information and abort if any bond or angle parameter is missing..
203 2699 : CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw)
204 : !-----------------------------------------------------------------------------
205 : ! 9. Urey-Bradley
206 : !-----------------------------------------------------------------------------
207 : CALL force_field_pack_ub(particle_set, molecule_kind_set, molecule_set, &
208 2699 : Ainfo, chm_info, inp_info, iw)
209 : !-----------------------------------------------------------------------------
210 : ! 10. Torsions
211 : !-----------------------------------------------------------------------------
212 : CALL force_field_pack_tors(particle_set, molecule_kind_set, molecule_set, &
213 2699 : Ainfo, chm_info, inp_info, gro_info, amb_info, iw)
214 : !-----------------------------------------------------------------------------
215 : ! 11. Impropers
216 : !-----------------------------------------------------------------------------
217 : CALL force_field_pack_impr(particle_set, molecule_kind_set, molecule_set, &
218 2699 : Ainfo, chm_info, inp_info, gro_info, iw)
219 : !-----------------------------------------------------------------------------
220 : ! 12. Out of plane bends
221 : !-----------------------------------------------------------------------------
222 : CALL force_field_pack_opbend(particle_set, molecule_kind_set, molecule_set, &
223 2699 : Ainfo, inp_info, iw)
224 : ! Give information only if any Urey-Bradley, Torsion, improper or opbend is missing
225 : ! continue calculation
226 2699 : CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw)
227 :
228 2699 : charges_section => section_vals_get_subs_vals(mm_section, "FORCEFIELD%CHARGES")
229 2699 : CALL section_vals_get(charges_section, explicit=explicit)
230 2699 : IF (.NOT. explicit) THEN
231 : !-----------------------------------------------------------------------------
232 : ! 13.a Set up atomic_kind_set()%fist_potentail%[qeff] and shell
233 : ! potential parameters
234 : !-----------------------------------------------------------------------------
235 : CALL force_field_pack_charge(atomic_kind_set, qmmm_env, fatal, iw, iw4, &
236 2691 : Ainfo, my_qmmm, inp_info)
237 : ! Give information only if charge is missing and abort..
238 2691 : CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw)
239 : ELSE
240 : !-----------------------------------------------------------------------------
241 : ! 13.b Setup a global array of classical charges - this avoids the packing and
242 : ! allows the usage of different charges for same atomic types
243 : !-----------------------------------------------------------------------------
244 : CALL force_field_pack_charges(charges, charges_section, particle_set, my_qmmm, &
245 8 : qmmm_env, inp_info, iw4)
246 : END IF
247 : !-----------------------------------------------------------------------------
248 : ! 14. Set up the radius of the electrostatic multipole in Fist
249 : !-----------------------------------------------------------------------------
250 2699 : CALL force_field_pack_radius(atomic_kind_set, iw, subsys_section)
251 : !-----------------------------------------------------------------------------
252 : ! 15. Set up the polarizable FF parameters
253 : !-----------------------------------------------------------------------------
254 2699 : CALL force_field_pack_pol(atomic_kind_set, iw, inp_info)
255 2699 : CALL force_field_pack_damp(atomic_kind_set, iw, inp_info)
256 : !-----------------------------------------------------------------------------
257 : ! 16. Set up Shell potential
258 : !-----------------------------------------------------------------------------
259 : CALL force_field_pack_shell(particle_set, atomic_kind_set, &
260 : molecule_kind_set, molecule_set, root_section, subsys_section, &
261 2699 : shell_particle_set, core_particle_set, cell, iw, inp_info)
262 2699 : IF (ff_type%do_nonbonded) THEN
263 : !-----------------------------------------------------------------------------
264 : ! 17. Set up potparm_nonbond14
265 : !-----------------------------------------------------------------------------
266 : ! Move the data from the info structures to potparm_nonbond
267 : CALL force_field_pack_nonbond14(atomic_kind_set, ff_type, qmmm_env, iw, Ainfo, &
268 2683 : chm_info, inp_info, gro_info, amb_info, potparm_nonbond14, ewald_env)
269 : ! Give information if any 1-4 is missing.. continue calculation..
270 2683 : CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw)
271 : ! Create the spline data
272 2683 : CALL section_vals_val_get(mm_section, "FORCEFIELD%ZBL_SCATTERING", l_val=do_zbl)
273 : CALL force_field_pack_splines(atomic_kind_set, ff_type, iw2, iw3, iw4, &
274 2683 : potparm_nonbond14, do_zbl, nonbonded_type="NONBONDED14")
275 : !-----------------------------------------------------------------------------
276 : ! 18. Set up potparm_nonbond
277 : !-----------------------------------------------------------------------------
278 : ! Move the data from the info structures to potparm_nonbond
279 : CALL force_field_pack_nonbond(atomic_kind_set, ff_type, qmmm_env, &
280 : fatal, iw, Ainfo, chm_info, inp_info, gro_info, amb_info, &
281 2683 : potparm_nonbond, ewald_env)
282 : ! Give information and abort if any pair potential spline is missing..
283 2683 : CALL release_FF_missing_par(fatal, ignore_fatal, AInfo, output_unit, iw)
284 : ! Create the spline data
285 2683 : CALL section_vals_val_get(mm_section, "FORCEFIELD%ZBL_SCATTERING", l_val=do_zbl)
286 : CALL force_field_pack_splines(atomic_kind_set, ff_type, iw2, iw3, iw4, &
287 2683 : potparm_nonbond, do_zbl, nonbonded_type="NONBONDED")
288 : END IF
289 : !-----------------------------------------------------------------------------
290 : ! 19. Create nonbond environment
291 : !-----------------------------------------------------------------------------
292 2699 : CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
293 : CALL section_vals_val_get(mm_section, "NEIGHBOR_LISTS%VERLET_SKIN", &
294 2699 : r_val=verlet_skin)
295 2699 : ALLOCATE (fist_nonbond_env)
296 : CALL fist_nonbond_env_create(fist_nonbond_env, atomic_kind_set, &
297 : potparm_nonbond14, potparm_nonbond, ff_type%do_nonbonded, &
298 : ff_type%do_electrostatics, verlet_skin, ewald_rcut, ff_type%ei_scale14, &
299 2699 : ff_type%vdw_scale14, ff_type%shift_cutoff)
300 2699 : CALL fist_nonbond_env_set(fist_nonbond_env, charges=charges)
301 : ! Compute the electrostatic interaction cutoffs.
302 : CALL force_field_pack_eicut(atomic_kind_set, ff_type, potparm_nonbond, &
303 2699 : ewald_env, iw)
304 :
305 2699 : CALL cp_print_key_finished_output(iw4, logger, mm_section, "PRINT%PROGRAM_RUN_INFO")
306 2699 : CALL cp_print_key_finished_output(iw3, logger, mm_section, "PRINT%FF_INFO/SPLINE_DATA")
307 2699 : CALL cp_print_key_finished_output(iw2, logger, mm_section, "PRINT%FF_INFO/SPLINE_INFO")
308 2699 : CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO")
309 :
310 2699 : CALL timestop(handle)
311 :
312 2699 : END SUBROUTINE force_field_pack
313 :
314 : ! **************************************************************************************************
315 : !> \brief Store informations on possible missing ForceFields parameters
316 : !> \param fatal ...
317 : !> \param ignore_fatal ...
318 : !> \param array ...
319 : !> \param output_unit ...
320 : !> \param iw ...
321 : ! **************************************************************************************************
322 13455 : SUBROUTINE release_FF_missing_par(fatal, ignore_fatal, array, output_unit, iw)
323 : LOGICAL, INTENT(INOUT) :: fatal, ignore_fatal
324 : CHARACTER(LEN=default_string_length), &
325 : DIMENSION(:), POINTER :: array
326 : INTEGER, INTENT(IN) :: output_unit, iw
327 :
328 : INTEGER :: i
329 :
330 13455 : IF (ASSOCIATED(array)) THEN
331 3286 : IF (output_unit > 0) THEN
332 2002 : WRITE (output_unit, *)
333 : WRITE (output_unit, '(T2,"FORCEFIELD|",A)') &
334 2002 : " WARNING: A non Critical ForceField parameter is missing! CP2K GOES!", &
335 2002 : " Non critical parameters are:Urey-Bradley,Torsions,Impropers, Opbends and 1-4", &
336 4004 : " All missing parameters will not contribute to the potential energy!"
337 2002 : IF (fatal .OR. iw > 0) THEN
338 311 : WRITE (output_unit, *)
339 3033 : DO i = 1, SIZE(array)
340 4724 : WRITE (output_unit, '(A)') array(i)
341 : END DO
342 : END IF
343 2002 : IF (.NOT. fatal .AND. iw < 0) THEN
344 : WRITE (output_unit, '(T2,"FORCEFIELD|",A)') &
345 1691 : " Activate the print key FF_INFO to have a list of missing parameters"
346 1691 : WRITE (output_unit, *)
347 : END IF
348 : END IF
349 3286 : DEALLOCATE (array)
350 : END IF
351 13455 : IF (fatal) THEN
352 62 : IF (ignore_fatal) THEN
353 62 : IF (output_unit > 0) THEN
354 43 : WRITE (output_unit, *)
355 : WRITE (output_unit, '(T2,"FORCEFIELD|",A)') &
356 43 : " WARNING: Ignoring missing critical FF parameters! CP2K GOES!", &
357 43 : " Critical parameters are: Bonds, Bends and Charges", &
358 86 : " All missing parameters will not contribute to the potential energy!"
359 : END IF
360 : ELSE
361 0 : CPABORT("Missing critical ForceField parameters!")
362 : END IF
363 : END IF
364 13455 : END SUBROUTINE release_FF_missing_par
365 :
366 : ! **************************************************************************************************
367 : !> \brief Compute the total qeff charges for each molecule kind and total system
368 : !> \param particle_set ...
369 : !> \param molecule_kind_set ...
370 : !> \param molecule_set ...
371 : !> \param mm_section ...
372 : !> \param charges ...
373 : ! **************************************************************************************************
374 2699 : SUBROUTINE force_field_qeff_output(particle_set, molecule_kind_set, &
375 : molecule_set, mm_section, charges)
376 :
377 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
378 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
379 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
380 : TYPE(section_vals_type), POINTER :: mm_section
381 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
382 :
383 : CHARACTER(len=*), PARAMETER :: routineN = 'force_field_qeff_output'
384 :
385 : CHARACTER(LEN=default_string_length) :: atmname, molname
386 : INTEGER :: first, handle, iatom, imol, iw, j, jatom
387 : LOGICAL :: shell_active
388 : REAL(KIND=dp) :: mass, mass_mol, mass_sum, qeff, &
389 : qeff_mol, qeff_sum
390 2699 : TYPE(atom_type), DIMENSION(:), POINTER :: atom_list
391 : TYPE(atomic_kind_type), POINTER :: atomic_kind
392 : TYPE(cp_logger_type), POINTER :: logger
393 : TYPE(molecule_kind_type), POINTER :: molecule_kind
394 : TYPE(molecule_type), POINTER :: molecule
395 : TYPE(shell_kind_type), POINTER :: shell
396 :
397 2699 : CALL timeset(routineN, handle)
398 :
399 2699 : NULLIFY (logger)
400 2699 : logger => cp_get_default_logger()
401 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
402 2699 : extension=".mmLog")
403 :
404 : qeff = 0.0_dp
405 : qeff_mol = 0.0_dp
406 2699 : qeff_sum = 0.0_dp
407 2699 : mass_sum = 0.0_dp
408 : !-----------------------------------------------------------------------------
409 : ! 1. Sum of [qeff,mass] for each molecule_kind
410 : !-----------------------------------------------------------------------------
411 75167 : DO imol = 1, SIZE(molecule_kind_set)
412 72468 : qeff_mol = 0.0_dp
413 72468 : mass_mol = 0.0_dp
414 72468 : molecule_kind => molecule_kind_set(imol)
415 :
416 72468 : j = molecule_kind_set(imol)%molecule_list(1)
417 72468 : molecule => molecule_set(j)
418 72468 : CALL get_molecule(molecule=molecule, first_atom=first)
419 :
420 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
421 72468 : name=molname, atom_list=atom_list)
422 256935 : DO iatom = 1, SIZE(atom_list)
423 184467 : atomic_kind => atom_list(iatom)%atomic_kind
424 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
425 184467 : name=atmname, qeff=qeff, mass=mass, shell_active=shell_active, shell=shell)
426 184467 : IF (shell_active) qeff = shell%charge_core + shell%charge_shell
427 184467 : IF (ASSOCIATED(charges)) THEN
428 30 : jatom = first - 1 + iatom
429 30 : qeff = charges(jatom)
430 : END IF
431 184467 : IF (iw > 0) WRITE (iw, *) " atom ", iatom, " ", TRIM(atmname), " charge = ", qeff
432 184467 : qeff_mol = qeff_mol + qeff
433 441402 : mass_mol = mass_mol + mass
434 : END DO
435 72468 : CALL set_molecule_kind(molecule_kind=molecule_kind, charge=qeff_mol, mass=mass_mol)
436 75167 : IF (iw > 0) WRITE (iw, *) " Mol Kind ", TRIM(molname), " charge = ", qeff_mol
437 : END DO
438 : !-----------------------------------------------------------------------------
439 : ! 2. Sum of [qeff,mass] for particle_set
440 : !-----------------------------------------------------------------------------
441 668721 : DO iatom = 1, SIZE(particle_set)
442 666022 : atomic_kind => particle_set(iatom)%atomic_kind
443 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
444 666022 : name=atmname, qeff=qeff, mass=mass, shell_active=shell_active, shell=shell)
445 666022 : IF (shell_active) qeff = shell%charge_core + shell%charge_shell
446 666022 : IF (ASSOCIATED(charges)) THEN
447 36 : qeff = charges(iatom)
448 : END IF
449 675366 : IF (iw > 0) WRITE (iw, *) " atom ", iatom, " ", TRIM(atmname), &
450 18688 : " charge = ", qeff
451 666022 : qeff_sum = qeff_sum + qeff
452 1334743 : mass_sum = mass_sum + mass
453 : END DO
454 2699 : IF (iw > 0) WRITE (iw, '(A,F20.10)') " Total system charge = ", qeff_sum
455 2699 : IF (iw > 0) WRITE (iw, '(A,F20.10)') " Total system mass = ", mass_sum
456 :
457 : CALL cp_print_key_finished_output(iw, logger, mm_section, &
458 2699 : "PRINT%FF_INFO")
459 :
460 2699 : CALL timestop(handle)
461 :
462 2699 : END SUBROUTINE force_field_qeff_output
463 :
464 : ! **************************************************************************************************
465 : !> \brief Removes UNSET force field types
466 : !> \param molecule_kind_set ...
467 : !> \param mm_section ...
468 : ! **************************************************************************************************
469 2699 : SUBROUTINE clean_intra_force_kind(molecule_kind_set, mm_section)
470 :
471 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
472 : TYPE(section_vals_type), POINTER :: mm_section
473 :
474 : CHARACTER(len=*), PARAMETER :: routineN = 'clean_intra_force_kind'
475 :
476 : INTEGER :: atm2_a, atm2_b, atm2_c, atm_a, atm_b, atm_c, atm_d, counter, handle, i, ibend, &
477 : ibond, icolv, ig3x3, ig4x6, iimpr, ikind, iopbend, itorsion, iub, iw, j, k, nbend, nbond, &
478 : newkind, ng3x3, ng4x6, nimpr, nopbend, ntorsion, nub
479 2699 : INTEGER, POINTER :: bad1(:), bad2(:)
480 : LOGICAL :: unsetme, valid_kind
481 2699 : TYPE(bend_kind_type), DIMENSION(:), POINTER :: bend_kind_set, new_bend_kind_set
482 2699 : TYPE(bend_type), DIMENSION(:), POINTER :: bend_list, new_bend_list
483 2699 : TYPE(bond_kind_type), DIMENSION(:), POINTER :: bond_kind_set, new_bond_kind_set
484 2699 : TYPE(bond_type), DIMENSION(:), POINTER :: bond_list, new_bond_list
485 : TYPE(colvar_constraint_type), DIMENSION(:), &
486 2699 : POINTER :: colv_list
487 : TYPE(cp_logger_type), POINTER :: logger
488 2699 : TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
489 2699 : TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
490 2699 : TYPE(impr_kind_type), DIMENSION(:), POINTER :: impr_kind_set, new_impr_kind_set
491 2699 : TYPE(impr_type), DIMENSION(:), POINTER :: impr_list, new_impr_list
492 : TYPE(molecule_kind_type), POINTER :: molecule_kind
493 2699 : TYPE(opbend_kind_type), DIMENSION(:), POINTER :: new_opbend_kind_set, opbend_kind_set
494 2699 : TYPE(opbend_type), DIMENSION(:), POINTER :: new_opbend_list, opbend_list
495 2699 : TYPE(torsion_kind_type), DIMENSION(:), POINTER :: new_torsion_kind_set, torsion_kind_set
496 2699 : TYPE(torsion_type), DIMENSION(:), POINTER :: new_torsion_list, torsion_list
497 2699 : TYPE(ub_kind_type), DIMENSION(:), POINTER :: new_ub_kind_set, ub_kind_set
498 2699 : TYPE(ub_type), DIMENSION(:), POINTER :: new_ub_list, ub_list
499 :
500 2699 : CALL timeset(routineN, handle)
501 2699 : NULLIFY (logger)
502 2699 : logger => cp_get_default_logger()
503 : iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
504 2699 : extension=".mmLog")
505 : !-----------------------------------------------------------------------------
506 : ! 1. Lets Tag the unwanted bonds due to the use of distance constraint
507 : !-----------------------------------------------------------------------------
508 75167 : DO ikind = 1, SIZE(molecule_kind_set)
509 72468 : molecule_kind => molecule_kind_set(ikind)
510 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
511 : colv_list=colv_list, &
512 : nbond=nbond, &
513 72468 : bond_list=bond_list)
514 75167 : IF (ASSOCIATED(colv_list)) THEN
515 834 : DO icolv = 1, SIZE(colv_list)
516 382 : IF ((colv_list(icolv)%type_id == dist_colvar_id) .AND. &
517 452 : ((.NOT. colv_list(icolv)%use_points) .OR. (SIZE(colv_list(icolv)%i_atoms) == 2))) THEN
518 284 : atm_a = colv_list(icolv)%i_atoms(1)
519 284 : atm_b = colv_list(icolv)%i_atoms(2)
520 3156 : DO ibond = 1, nbond
521 2872 : unsetme = .FALSE.
522 2872 : atm2_a = bond_list(ibond)%a
523 2872 : atm2_b = bond_list(ibond)%b
524 2872 : IF (atm2_a == atm_a .AND. atm2_b == atm_b) unsetme = .TRUE.
525 2872 : IF (atm2_a == atm_b .AND. atm2_b == atm_a) unsetme = .TRUE.
526 3254 : IF (unsetme) bond_list(ibond)%id_type = do_ff_undef
527 : END DO
528 : END IF
529 : END DO
530 : END IF
531 : END DO
532 : !-----------------------------------------------------------------------------
533 : ! 2. Lets Tag the unwanted bends due to the use of distance constraint
534 : !-----------------------------------------------------------------------------
535 75167 : DO ikind = 1, SIZE(molecule_kind_set)
536 72468 : molecule_kind => molecule_kind_set(ikind)
537 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
538 : colv_list=colv_list, &
539 : nbend=nbend, &
540 72468 : bend_list=bend_list)
541 75167 : IF (ASSOCIATED(colv_list)) THEN
542 1606 : DO ibend = 1, nbend
543 1154 : unsetme = .FALSE.
544 1154 : atm_a = bend_list(ibend)%a
545 1154 : atm_b = bend_list(ibend)%b
546 1154 : atm_c = bend_list(ibend)%c
547 6774 : DO icolv = 1, SIZE(colv_list)
548 5656 : IF ((colv_list(icolv)%type_id == dist_colvar_id) .AND. &
549 1118 : ((.NOT. colv_list(icolv)%use_points) .OR. (SIZE(colv_list(icolv)%i_atoms) == 2))) THEN
550 5030 : atm2_a = colv_list(icolv)%i_atoms(1)
551 5030 : atm2_b = colv_list(icolv)%i_atoms(2)
552 : ! Check that the bonds we're constraining does not involve atoms defining
553 : ! a bend..
554 5030 : IF (((atm2_a == atm_a) .AND. (atm2_b == atm_c)) .OR. &
555 : ((atm2_a == atm_c) .AND. (atm2_b == atm_a))) THEN
556 : unsetme = .TRUE.
557 : EXIT
558 : END IF
559 : END IF
560 : END DO
561 1606 : IF (unsetme) bend_list(ibend)%id_type = do_ff_undef
562 : END DO
563 : END IF
564 : END DO
565 : !-----------------------------------------------------------------------------
566 : ! 3. Lets Tag the unwanted bonds due to the use of 3x3
567 : !-----------------------------------------------------------------------------
568 75167 : DO ikind = 1, SIZE(molecule_kind_set)
569 72468 : molecule_kind => molecule_kind_set(ikind)
570 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
571 : ng3x3=ng3x3, &
572 : g3x3_list=g3x3_list, &
573 : nbond=nbond, &
574 72468 : bond_list=bond_list)
575 75303 : DO ig3x3 = 1, ng3x3
576 136 : atm_a = g3x3_list(ig3x3)%a
577 136 : atm_b = g3x3_list(ig3x3)%b
578 136 : atm_c = g3x3_list(ig3x3)%c
579 72916 : DO ibond = 1, nbond
580 312 : unsetme = .FALSE.
581 312 : atm2_a = bond_list(ibond)%a
582 312 : atm2_b = bond_list(ibond)%b
583 312 : IF (atm2_a == atm_a .AND. atm2_b == atm_b) unsetme = .TRUE.
584 312 : IF (atm2_a == atm_a .AND. atm2_b == atm_c) unsetme = .TRUE.
585 312 : IF (atm2_a == atm_c .AND. atm2_b == atm_c) unsetme = .TRUE.
586 448 : IF (unsetme) bond_list(ibond)%id_type = do_ff_undef
587 : END DO
588 : END DO
589 : END DO
590 : !-----------------------------------------------------------------------------
591 : ! 4. Lets Tag the unwanted bends due to the use of 3x3
592 : !-----------------------------------------------------------------------------
593 75167 : DO ikind = 1, SIZE(molecule_kind_set)
594 72468 : molecule_kind => molecule_kind_set(ikind)
595 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
596 : ng3x3=ng3x3, &
597 : g3x3_list=g3x3_list, &
598 : nbend=nbend, &
599 72468 : bend_list=bend_list)
600 75303 : DO ig3x3 = 1, ng3x3
601 136 : atm_a = g3x3_list(ig3x3)%a
602 136 : atm_b = g3x3_list(ig3x3)%b
603 136 : atm_c = g3x3_list(ig3x3)%c
604 72724 : DO ibend = 1, nbend
605 120 : unsetme = .FALSE.
606 120 : atm2_a = bend_list(ibend)%a
607 120 : atm2_b = bend_list(ibend)%b
608 120 : atm2_c = bend_list(ibend)%c
609 120 : IF (atm2_a == atm_a .AND. atm2_b == atm_b .AND. atm2_c == atm_c) unsetme = .TRUE.
610 120 : IF (atm2_a == atm_a .AND. atm2_b == atm_c .AND. atm2_c == atm_b) unsetme = .TRUE.
611 120 : IF (atm2_a == atm_b .AND. atm2_b == atm_a .AND. atm2_c == atm_c) unsetme = .TRUE.
612 120 : IF (atm2_a == atm_b .AND. atm2_b == atm_c .AND. atm2_c == atm_a) unsetme = .TRUE.
613 256 : IF (unsetme) bend_list(ibend)%id_type = do_ff_undef
614 : END DO
615 : END DO
616 : END DO
617 : !-----------------------------------------------------------------------------
618 : ! 5. Lets Tag the unwanted bonds due to the use of 4x6
619 : !-----------------------------------------------------------------------------
620 75167 : DO ikind = 1, SIZE(molecule_kind_set)
621 72468 : molecule_kind => molecule_kind_set(ikind)
622 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
623 : ng4x6=ng4x6, &
624 : g4x6_list=g4x6_list, &
625 : nbond=nbond, &
626 72468 : bond_list=bond_list)
627 75179 : DO ig4x6 = 1, ng4x6
628 12 : atm_a = g4x6_list(ig4x6)%a
629 12 : atm_b = g4x6_list(ig4x6)%b
630 12 : atm_c = g4x6_list(ig4x6)%c
631 12 : atm_d = g4x6_list(ig4x6)%d
632 72516 : DO ibond = 1, nbond
633 36 : unsetme = .FALSE.
634 36 : atm2_a = bond_list(ibond)%a
635 36 : atm2_b = bond_list(ibond)%b
636 36 : IF (atm2_a == atm_a .AND. atm2_b == atm_b) unsetme = .TRUE.
637 36 : IF (atm2_a == atm_a .AND. atm2_b == atm_c) unsetme = .TRUE.
638 36 : IF (atm2_a == atm_a .AND. atm2_b == atm_d) unsetme = .TRUE.
639 48 : IF (unsetme) bond_list(ibond)%id_type = do_ff_undef
640 : END DO
641 : END DO
642 : END DO
643 : !-----------------------------------------------------------------------------
644 : ! 6. Lets Tag the unwanted bends due to the use of 4x6
645 : !-----------------------------------------------------------------------------
646 75167 : DO ikind = 1, SIZE(molecule_kind_set)
647 72468 : molecule_kind => molecule_kind_set(ikind)
648 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
649 : ng4x6=ng4x6, &
650 : g4x6_list=g4x6_list, &
651 : nbend=nbend, &
652 72468 : bend_list=bend_list)
653 75179 : DO ig4x6 = 1, ng4x6
654 12 : atm_a = g4x6_list(ig4x6)%a
655 12 : atm_b = g4x6_list(ig4x6)%b
656 12 : atm_c = g4x6_list(ig4x6)%c
657 12 : atm_d = g4x6_list(ig4x6)%d
658 72516 : DO ibend = 1, nbend
659 36 : unsetme = .FALSE.
660 36 : atm2_a = bend_list(ibend)%a
661 36 : atm2_b = bend_list(ibend)%b
662 36 : atm2_c = bend_list(ibend)%c
663 36 : IF (atm2_a == atm_b .AND. atm2_b == atm_a .AND. atm2_c == atm_c) unsetme = .TRUE.
664 36 : IF (atm2_a == atm_b .AND. atm2_b == atm_a .AND. atm2_c == atm_d) unsetme = .TRUE.
665 36 : IF (atm2_a == atm_c .AND. atm2_b == atm_a .AND. atm2_c == atm_d) unsetme = .TRUE.
666 48 : IF (unsetme) bend_list(ibend)%id_type = do_ff_undef
667 : END DO
668 : END DO
669 : END DO
670 : !-----------------------------------------------------------------------------
671 : ! 7. Count the number of UNSET bond kinds there are
672 : !-----------------------------------------------------------------------------
673 75167 : DO ikind = 1, SIZE(molecule_kind_set)
674 72468 : molecule_kind => molecule_kind_set(ikind)
675 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
676 : nbond=nbond, &
677 : bond_kind_set=bond_kind_set, &
678 72468 : bond_list=bond_list)
679 75167 : IF (nbond > 0) THEN
680 29760 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old BOND Count: ", &
681 606 : SIZE(bond_list), SIZE(bond_kind_set)
682 31867 : IF (iw > 0) WRITE (iw, '(2I6)') (bond_list(ibond)%a, bond_list(ibond)%b, ibond=1, SIZE(bond_list))
683 29457 : NULLIFY (bad1, bad2)
684 88371 : ALLOCATE (bad1(SIZE(bond_kind_set)))
685 96530 : bad1(:) = 0
686 96530 : DO ibond = 1, SIZE(bond_kind_set)
687 67073 : unsetme = .FALSE.
688 67073 : IF (bond_kind_set(ibond)%id_type == do_ff_undef) unsetme = .TRUE.
689 67073 : valid_kind = .FALSE.
690 350648 : DO i = 1, SIZE(bond_list)
691 349234 : IF (bond_list(i)%id_type /= do_ff_undef .AND. &
692 1414 : bond_list(i)%bond_kind%kind_number == ibond) THEN
693 : valid_kind = .TRUE.
694 : EXIT
695 : END IF
696 : END DO
697 67073 : IF (.NOT. valid_kind) unsetme = .TRUE.
698 96530 : IF (unsetme) bad1(ibond) = 1
699 : END DO
700 96530 : IF (SUM(bad1) /= 0) THEN
701 2770 : counter = SIZE(bond_kind_set) - SUM(bad1)
702 804 : CALL allocate_bond_kind_set(new_bond_kind_set, counter)
703 804 : counter = 0
704 2770 : DO ibond = 1, SIZE(bond_kind_set)
705 2770 : IF (bad1(ibond) == 0) THEN
706 546 : counter = counter + 1
707 546 : new_bond_kind_set(counter) = bond_kind_set(ibond)
708 : END IF
709 : END DO
710 804 : counter = 0
711 2412 : ALLOCATE (bad2(SIZE(bond_list)))
712 4338 : bad2(:) = 0
713 4338 : DO ibond = 1, SIZE(bond_list)
714 3534 : unsetme = .FALSE.
715 3534 : IF (bond_list(ibond)%bond_kind%id_type == do_ff_undef) unsetme = .TRUE.
716 3534 : IF (bond_list(ibond)%id_type == do_ff_undef) unsetme = .TRUE.
717 4338 : IF (unsetme) bad2(ibond) = 1
718 : END DO
719 4338 : IF (SUM(bad2) /= 0) THEN
720 4338 : counter = SIZE(bond_list) - SUM(bad2)
721 2648 : ALLOCATE (new_bond_list(counter))
722 804 : counter = 0
723 4338 : DO ibond = 1, SIZE(bond_list)
724 4338 : IF (bad2(ibond) == 0) THEN
725 926 : counter = counter + 1
726 926 : new_bond_list(counter) = bond_list(ibond)
727 926 : newkind = bond_list(ibond)%bond_kind%kind_number
728 5958 : newkind = newkind - SUM(bad1(1:newkind))
729 926 : new_bond_list(counter)%bond_kind => new_bond_kind_set(newkind)
730 : END IF
731 : END DO
732 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
733 : nbond=SIZE(new_bond_list), &
734 : bond_kind_set=new_bond_kind_set, &
735 804 : bond_list=new_bond_list)
736 1350 : DO ibond = 1, SIZE(new_bond_kind_set)
737 1350 : new_bond_kind_set(ibond)%kind_number = ibond
738 : END DO
739 : END IF
740 804 : DEALLOCATE (bad2)
741 804 : CALL deallocate_bond_kind_set(bond_kind_set)
742 804 : DEALLOCATE (bond_list)
743 808 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New BOND Count: ", &
744 8 : SIZE(new_bond_list), SIZE(new_bond_kind_set)
745 808 : IF (iw > 0) WRITE (iw, '(2I6)') (new_bond_list(ibond)%a, new_bond_list(ibond)%b, &
746 8 : ibond=1, SIZE(new_bond_list))
747 : END IF
748 29457 : DEALLOCATE (bad1)
749 : END IF
750 : END DO
751 : !-----------------------------------------------------------------------------
752 : ! 8. Count the number of UNSET bend kinds there are
753 : !-----------------------------------------------------------------------------
754 75167 : DO ikind = 1, SIZE(molecule_kind_set)
755 72468 : molecule_kind => molecule_kind_set(ikind)
756 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
757 : nbend=nbend, &
758 : bend_kind_set=bend_kind_set, &
759 72468 : bend_list=bend_list)
760 75167 : IF (nbend > 0) THEN
761 29430 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old BEND Count: ", &
762 606 : SIZE(bend_list), SIZE(bend_kind_set)
763 33099 : IF (iw > 0) WRITE (iw, '(3I6)') (bend_list(ibend)%a, bend_list(ibend)%b, &
764 4275 : bend_list(ibend)%c, ibend=1, SIZE(bend_list))
765 29127 : NULLIFY (bad1, bad2)
766 87381 : ALLOCATE (bad1(SIZE(bend_kind_set)))
767 124357 : bad1(:) = 0
768 124357 : DO ibend = 1, SIZE(bend_kind_set)
769 95230 : unsetme = .FALSE.
770 95230 : IF (bend_kind_set(ibend)%id_type == do_ff_undef) unsetme = .TRUE.
771 95230 : valid_kind = .FALSE.
772 1012011 : DO i = 1, SIZE(bend_list)
773 1010705 : IF (bend_list(i)%id_type /= do_ff_undef .AND. &
774 1306 : bend_list(i)%bend_kind%kind_number == ibend) THEN
775 : valid_kind = .TRUE.
776 : EXIT
777 : END IF
778 : END DO
779 95230 : IF (.NOT. valid_kind) unsetme = .TRUE.
780 124357 : IF (unsetme) bad1(ibend) = 1
781 : END DO
782 124357 : IF (SUM(bad1) /= 0) THEN
783 2786 : counter = SIZE(bend_kind_set) - SUM(bad1)
784 606 : CALL allocate_bend_kind_set(new_bend_kind_set, counter)
785 606 : counter = 0
786 2786 : DO ibend = 1, SIZE(bend_kind_set)
787 2786 : IF (bad1(ibend) == 0) THEN
788 868 : counter = counter + 1
789 868 : new_bend_kind_set(counter) = bend_kind_set(ibend)
790 : END IF
791 : END DO
792 606 : counter = 0
793 1818 : ALLOCATE (bad2(SIZE(bend_list)))
794 4322 : bad2(:) = 0
795 4322 : DO ibend = 1, SIZE(bend_list)
796 3716 : unsetme = .FALSE.
797 3716 : IF (bend_list(ibend)%bend_kind%id_type == do_ff_undef) unsetme = .TRUE.
798 3716 : IF (bend_list(ibend)%id_type == do_ff_undef) unsetme = .TRUE.
799 4322 : IF (unsetme) bad2(ibend) = 1
800 : END DO
801 4322 : IF (SUM(bad2) /= 0) THEN
802 4322 : counter = SIZE(bend_list) - SUM(bad2)
803 2900 : ALLOCATE (new_bend_list(counter))
804 606 : counter = 0
805 4322 : DO ibend = 1, SIZE(bend_list)
806 4322 : IF (bad2(ibend) == 0) THEN
807 1610 : counter = counter + 1
808 1610 : new_bend_list(counter) = bend_list(ibend)
809 1610 : newkind = bend_list(ibend)%bend_kind%kind_number
810 18122 : newkind = newkind - SUM(bad1(1:newkind))
811 1610 : new_bend_list(counter)%bend_kind => new_bend_kind_set(newkind)
812 : END IF
813 : END DO
814 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
815 : nbend=SIZE(new_bend_list), &
816 : bend_kind_set=new_bend_kind_set, &
817 606 : bend_list=new_bend_list)
818 1474 : DO ibend = 1, SIZE(new_bend_kind_set)
819 1474 : new_bend_kind_set(ibend)%kind_number = ibend
820 : END DO
821 : END IF
822 606 : DEALLOCATE (bad2)
823 606 : CALL deallocate_bend_kind_set(bend_kind_set)
824 606 : DEALLOCATE (bend_list)
825 610 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New BEND Count: ", &
826 8 : SIZE(new_bend_list), SIZE(new_bend_kind_set)
827 606 : IF (iw > 0) WRITE (iw, '(3I6)') (new_bend_list(ibend)%a, new_bend_list(ibend)%b, &
828 4 : new_bend_list(ibend)%c, ibend=1, SIZE(new_bend_list))
829 : END IF
830 29127 : DEALLOCATE (bad1)
831 : END IF
832 : END DO
833 : !-----------------------------------------------------------------------------
834 : ! 9. Count the number of UNSET Urey-Bradley kinds there are
835 : !-----------------------------------------------------------------------------
836 75167 : DO ikind = 1, SIZE(molecule_kind_set)
837 72468 : molecule_kind => molecule_kind_set(ikind)
838 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
839 : nub=nub, &
840 : ub_kind_set=ub_kind_set, &
841 72468 : ub_list=ub_list)
842 75167 : IF (nub > 0) THEN
843 29416 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old UB Count: ", &
844 606 : SIZE(ub_list), SIZE(ub_kind_set)
845 33085 : IF (iw > 0) WRITE (iw, '(3I6)') (ub_list(iub)%a, ub_list(iub)%b, &
846 4275 : ub_list(iub)%c, iub=1, SIZE(ub_list))
847 29113 : NULLIFY (bad1, bad2)
848 87339 : ALLOCATE (bad1(SIZE(ub_kind_set)))
849 124185 : bad1(:) = 0
850 124185 : DO iub = 1, SIZE(ub_kind_set)
851 95072 : unsetme = .FALSE.
852 95072 : IF (ub_kind_set(iub)%id_type == do_ff_undef) unsetme = .TRUE.
853 95072 : valid_kind = .FALSE.
854 1955458 : DO i = 1, SIZE(ub_list)
855 1868756 : IF (ub_list(i)%id_type /= do_ff_undef .AND. &
856 86702 : ub_list(i)%ub_kind%kind_number == iub) THEN
857 : valid_kind = .TRUE.
858 : EXIT
859 : END IF
860 : END DO
861 95072 : IF (.NOT. valid_kind) unsetme = .TRUE.
862 124185 : IF (unsetme) bad1(iub) = 1
863 : END DO
864 124185 : IF (SUM(bad1) /= 0) THEN
865 124149 : counter = SIZE(ub_kind_set) - SUM(bad1)
866 29105 : CALL allocate_ub_kind_set(new_ub_kind_set, counter)
867 29105 : counter = 0
868 124149 : DO iub = 1, SIZE(ub_kind_set)
869 124149 : IF (bad1(iub) == 0) THEN
870 8342 : counter = counter + 1
871 8342 : new_ub_kind_set(counter) = ub_kind_set(iub)
872 : END IF
873 : END DO
874 29105 : counter = 0
875 87315 : ALLOCATE (bad2(SIZE(ub_list)))
876 170197 : bad2(:) = 0
877 170197 : DO iub = 1, SIZE(ub_list)
878 141092 : unsetme = .FALSE.
879 141092 : IF (ub_list(iub)%ub_kind%id_type == do_ff_undef) unsetme = .TRUE.
880 141092 : IF (ub_list(iub)%id_type == do_ff_undef) unsetme = .TRUE.
881 170197 : IF (unsetme) bad2(iub) = 1
882 : END DO
883 170197 : IF (SUM(bad2) /= 0) THEN
884 170197 : counter = SIZE(ub_list) - SUM(bad2)
885 79726 : ALLOCATE (new_ub_list(counter))
886 29105 : counter = 0
887 170197 : DO iub = 1, SIZE(ub_list)
888 170197 : IF (bad2(iub) == 0) THEN
889 19972 : counter = counter + 1
890 19972 : new_ub_list(counter) = ub_list(iub)
891 19972 : newkind = ub_list(iub)%ub_kind%kind_number
892 256656 : newkind = newkind - SUM(bad1(1:newkind))
893 19972 : new_ub_list(counter)%ub_kind => new_ub_kind_set(newkind)
894 : END IF
895 : END DO
896 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
897 : nub=SIZE(new_ub_list), &
898 : ub_kind_set=new_ub_kind_set, &
899 29105 : ub_list=new_ub_list)
900 37447 : DO iub = 1, SIZE(new_ub_kind_set)
901 37447 : new_ub_kind_set(iub)%kind_number = iub
902 : END DO
903 : END IF
904 29105 : DEALLOCATE (bad2)
905 29105 : CALL ub_kind_dealloc_ref(ub_kind_set)
906 29105 : DEALLOCATE (ub_list)
907 29408 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New UB Count: ", &
908 606 : SIZE(new_ub_list), SIZE(new_ub_kind_set)
909 29243 : IF (iw > 0) WRITE (iw, '(3I6)') (new_ub_list(iub)%a, new_ub_list(iub)%b, &
910 441 : new_ub_list(iub)%c, iub=1, SIZE(new_ub_list))
911 : END IF
912 29113 : DEALLOCATE (bad1)
913 : END IF
914 : END DO
915 : !-----------------------------------------------------------------------------
916 : ! 10. Count the number of UNSET torsion kinds there are
917 : !-----------------------------------------------------------------------------
918 75167 : DO ikind = 1, SIZE(molecule_kind_set)
919 72468 : molecule_kind => molecule_kind_set(ikind)
920 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
921 : ntorsion=ntorsion, &
922 : torsion_kind_set=torsion_kind_set, &
923 72468 : torsion_list=torsion_list)
924 75167 : IF (ntorsion > 0) THEN
925 5813 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old TORSION Count: ", &
926 538 : SIZE(torsion_list), SIZE(torsion_kind_set)
927 10262 : IF (iw > 0) WRITE (iw, '(4I6)') (torsion_list(itorsion)%a, torsion_list(itorsion)%b, &
928 4987 : torsion_list(itorsion)%c, torsion_list(itorsion)%d, itorsion=1, SIZE(torsion_list))
929 5544 : NULLIFY (bad1, bad2)
930 16632 : ALLOCATE (bad1(SIZE(torsion_kind_set)))
931 98627 : bad1(:) = 0
932 98627 : DO itorsion = 1, SIZE(torsion_kind_set)
933 93083 : unsetme = .FALSE.
934 93083 : IF (torsion_kind_set(itorsion)%id_type == do_ff_undef) unsetme = .TRUE.
935 93083 : valid_kind = .FALSE.
936 2320901 : DO i = 1, SIZE(torsion_list)
937 2306549 : IF (torsion_list(i)%id_type /= do_ff_undef .AND. &
938 14352 : torsion_list(i)%torsion_kind%kind_number == itorsion) THEN
939 : valid_kind = .TRUE.
940 : EXIT
941 : END IF
942 : END DO
943 93083 : IF (.NOT. valid_kind) unsetme = .TRUE.
944 98627 : IF (unsetme) bad1(itorsion) = 1
945 : END DO
946 98627 : IF (SUM(bad1) /= 0) THEN
947 18312 : counter = SIZE(torsion_kind_set) - SUM(bad1)
948 1026 : CALL allocate_torsion_kind_set(new_torsion_kind_set, counter)
949 1026 : counter = 0
950 18312 : DO itorsion = 1, SIZE(torsion_kind_set)
951 18312 : IF (bad1(itorsion) == 0) THEN
952 2934 : counter = counter + 1
953 2934 : new_torsion_kind_set(counter) = torsion_kind_set(itorsion)
954 2934 : i = SIZE(torsion_kind_set(itorsion)%m)
955 2934 : j = SIZE(torsion_kind_set(itorsion)%k)
956 2934 : k = SIZE(torsion_kind_set(itorsion)%phi0)
957 8802 : ALLOCATE (new_torsion_kind_set(counter)%m(i))
958 8802 : ALLOCATE (new_torsion_kind_set(counter)%k(i))
959 5868 : ALLOCATE (new_torsion_kind_set(counter)%phi0(i))
960 12772 : new_torsion_kind_set(counter)%m = torsion_kind_set(itorsion)%m
961 12772 : new_torsion_kind_set(counter)%k = torsion_kind_set(itorsion)%k
962 12772 : new_torsion_kind_set(counter)%phi0 = torsion_kind_set(itorsion)%phi0
963 : END IF
964 : END DO
965 1026 : counter = 0
966 3078 : ALLOCATE (bad2(SIZE(torsion_list)))
967 46332 : bad2(:) = 0
968 46332 : DO itorsion = 1, SIZE(torsion_list)
969 45306 : unsetme = .FALSE.
970 45306 : IF (torsion_list(itorsion)%torsion_kind%id_type == do_ff_undef) unsetme = .TRUE.
971 45306 : IF (torsion_list(itorsion)%id_type == do_ff_undef) unsetme = .TRUE.
972 46332 : IF (unsetme) bad2(itorsion) = 1
973 : END DO
974 46332 : IF (SUM(bad2) /= 0) THEN
975 46332 : counter = SIZE(torsion_list) - SUM(bad2)
976 9424 : ALLOCATE (new_torsion_list(counter))
977 1026 : counter = 0
978 46332 : DO itorsion = 1, SIZE(torsion_list)
979 46332 : IF (bad2(itorsion) == 0) THEN
980 7252 : counter = counter + 1
981 7252 : new_torsion_list(counter) = torsion_list(itorsion)
982 7252 : newkind = torsion_list(itorsion)%torsion_kind%kind_number
983 146518 : newkind = newkind - SUM(bad1(1:newkind))
984 7252 : new_torsion_list(counter)%torsion_kind => new_torsion_kind_set(newkind)
985 : END IF
986 : END DO
987 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
988 : ntorsion=SIZE(new_torsion_list), &
989 : torsion_kind_set=new_torsion_kind_set, &
990 1026 : torsion_list=new_torsion_list)
991 3960 : DO itorsion = 1, SIZE(new_torsion_kind_set)
992 3960 : new_torsion_kind_set(itorsion)%kind_number = itorsion
993 : END DO
994 : END IF
995 1026 : DEALLOCATE (bad2)
996 18312 : DO itorsion = 1, SIZE(torsion_kind_set)
997 18312 : CALL torsion_kind_dealloc_ref(torsion_kind_set(itorsion))
998 : END DO
999 1026 : DEALLOCATE (torsion_kind_set)
1000 1026 : DEALLOCATE (torsion_list)
1001 1160 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New TORSION Count: ", &
1002 268 : SIZE(new_torsion_list), SIZE(new_torsion_kind_set)
1003 1173 : IF (iw > 0) WRITE (iw, '(4I6)') (new_torsion_list(itorsion)%a, new_torsion_list(itorsion)%b, &
1004 428 : new_torsion_list(itorsion)%c, new_torsion_list(itorsion)%d, itorsion=1, &
1005 415 : SIZE(new_torsion_list))
1006 : END IF
1007 5544 : DEALLOCATE (bad1)
1008 : END IF
1009 : END DO
1010 : !-----------------------------------------------------------------------------
1011 : ! 11. Count the number of UNSET improper kinds there are
1012 : !-----------------------------------------------------------------------------
1013 75167 : DO ikind = 1, SIZE(molecule_kind_set)
1014 72468 : molecule_kind => molecule_kind_set(ikind)
1015 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
1016 : nimpr=nimpr, &
1017 : impr_kind_set=impr_kind_set, &
1018 72468 : impr_list=impr_list)
1019 75167 : IF (nimpr > 0) THEN
1020 1701 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old IMPROPER Count: ", &
1021 46 : SIZE(impr_list), SIZE(impr_kind_set)
1022 1678 : NULLIFY (bad1, bad2)
1023 5034 : ALLOCATE (bad1(SIZE(impr_kind_set)))
1024 6506 : bad1(:) = 0
1025 6506 : DO iimpr = 1, SIZE(impr_kind_set)
1026 4828 : unsetme = .FALSE.
1027 4828 : IF (impr_kind_set(iimpr)%id_type == do_ff_undef) unsetme = .TRUE.
1028 4828 : valid_kind = .FALSE.
1029 32300 : DO i = 1, SIZE(impr_list)
1030 30682 : IF (impr_list(i)%id_type /= do_ff_undef .AND. &
1031 1618 : impr_list(i)%impr_kind%kind_number == iimpr) THEN
1032 : valid_kind = .TRUE.
1033 : EXIT
1034 : END IF
1035 : END DO
1036 4828 : IF (.NOT. valid_kind) unsetme = .TRUE.
1037 6506 : IF (unsetme) bad1(iimpr) = 1
1038 : END DO
1039 6506 : IF (SUM(bad1) /= 0) THEN
1040 2138 : counter = SIZE(impr_kind_set) - SUM(bad1)
1041 396 : CALL allocate_impr_kind_set(new_impr_kind_set, counter)
1042 396 : counter = 0
1043 2138 : DO iimpr = 1, SIZE(impr_kind_set)
1044 2138 : IF (bad1(iimpr) == 0) THEN
1045 124 : counter = counter + 1
1046 124 : new_impr_kind_set(counter) = impr_kind_set(iimpr)
1047 : END IF
1048 : END DO
1049 : counter = 0
1050 1188 : ALLOCATE (bad2(SIZE(impr_list)))
1051 2594 : bad2(:) = 0
1052 2594 : DO iimpr = 1, SIZE(impr_list)
1053 2198 : unsetme = .FALSE.
1054 2198 : IF (impr_list(iimpr)%impr_kind%id_type == do_ff_undef) unsetme = .TRUE.
1055 2198 : IF (impr_list(iimpr)%id_type == do_ff_undef) unsetme = .TRUE.
1056 2594 : IF (unsetme) bad2(iimpr) = 1
1057 : END DO
1058 2594 : IF (SUM(bad2) /= 0) THEN
1059 2594 : counter = SIZE(impr_list) - SUM(bad2)
1060 1020 : ALLOCATE (new_impr_list(counter))
1061 396 : counter = 0
1062 2594 : DO iimpr = 1, SIZE(impr_list)
1063 2594 : IF (bad2(iimpr) == 0) THEN
1064 124 : counter = counter + 1
1065 124 : new_impr_list(counter) = impr_list(iimpr)
1066 124 : newkind = impr_list(iimpr)%impr_kind%kind_number
1067 324 : newkind = newkind - SUM(bad1(1:newkind))
1068 124 : new_impr_list(counter)%impr_kind => new_impr_kind_set(newkind)
1069 : END IF
1070 : END DO
1071 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
1072 : nimpr=SIZE(new_impr_list), &
1073 : impr_kind_set=new_impr_kind_set, &
1074 396 : impr_list=new_impr_list)
1075 520 : DO iimpr = 1, SIZE(new_impr_kind_set)
1076 520 : new_impr_kind_set(iimpr)%kind_number = iimpr
1077 : END DO
1078 : END IF
1079 396 : DEALLOCATE (bad2)
1080 2138 : DO iimpr = 1, SIZE(impr_kind_set)
1081 2138 : CALL impr_kind_dealloc_ref() !This Subroutine doesn't deallocate anything, maybe needs to be implemented
1082 : END DO
1083 396 : DEALLOCATE (impr_kind_set)
1084 396 : DEALLOCATE (impr_list)
1085 407 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New IMPROPER Count: ", &
1086 22 : SIZE(new_impr_list), SIZE(new_impr_kind_set)
1087 : END IF
1088 1678 : DEALLOCATE (bad1)
1089 : END IF
1090 : END DO
1091 : !-----------------------------------------------------------------------------
1092 : ! 11. Count the number of UNSET opbends kinds there are
1093 : !-----------------------------------------------------------------------------
1094 75167 : DO ikind = 1, SIZE(molecule_kind_set)
1095 72468 : molecule_kind => molecule_kind_set(ikind)
1096 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
1097 : nopbend=nopbend, &
1098 : opbend_kind_set=opbend_kind_set, &
1099 72468 : opbend_list=opbend_list)
1100 75167 : IF (nopbend > 0) THEN
1101 1701 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old OPBEND Count: ", &
1102 46 : SIZE(opbend_list), SIZE(opbend_kind_set)
1103 1678 : NULLIFY (bad1, bad2)
1104 5034 : ALLOCATE (bad1(SIZE(opbend_kind_set)))
1105 6506 : bad1(:) = 0
1106 6506 : DO iopbend = 1, SIZE(opbend_kind_set)
1107 4828 : unsetme = .FALSE.
1108 4828 : IF (opbend_kind_set(iopbend)%id_type == do_ff_undef) unsetme = .TRUE.
1109 4828 : valid_kind = .FALSE.
1110 39328 : DO i = 1, SIZE(opbend_list)
1111 34502 : IF (opbend_list(i)%id_type /= do_ff_undef .AND. &
1112 4826 : opbend_list(i)%opbend_kind%kind_number == iopbend) THEN
1113 : valid_kind = .TRUE.
1114 : EXIT
1115 : END IF
1116 : END DO
1117 4828 : IF (.NOT. valid_kind) unsetme = .TRUE.
1118 6506 : IF (unsetme) bad1(iopbend) = 1
1119 : END DO
1120 6506 : IF (SUM(bad1) /= 0) THEN
1121 6502 : counter = SIZE(opbend_kind_set) - SUM(bad1)
1122 1676 : CALL allocate_opbend_kind_set(new_opbend_kind_set, counter)
1123 1676 : counter = 0
1124 6502 : DO iopbend = 1, SIZE(opbend_kind_set)
1125 6502 : IF (bad1(iopbend) == 0) THEN
1126 0 : counter = counter + 1
1127 0 : new_opbend_kind_set(counter) = opbend_kind_set(iopbend)
1128 : END IF
1129 : END DO
1130 : counter = 0
1131 5028 : ALLOCATE (bad2(SIZE(opbend_list)))
1132 7154 : bad2(:) = 0
1133 7154 : DO iopbend = 1, SIZE(opbend_list)
1134 5478 : unsetme = .FALSE.
1135 5478 : IF (opbend_list(iopbend)%opbend_kind%id_type == do_ff_undef) unsetme = .TRUE.
1136 5478 : IF (opbend_list(iopbend)%id_type == do_ff_undef) unsetme = .TRUE.
1137 7154 : IF (unsetme) bad2(iopbend) = 1
1138 : END DO
1139 7154 : IF (SUM(bad2) /= 0) THEN
1140 7154 : counter = SIZE(opbend_list) - SUM(bad2)
1141 3352 : ALLOCATE (new_opbend_list(counter))
1142 1676 : counter = 0
1143 7154 : DO iopbend = 1, SIZE(opbend_list)
1144 7154 : IF (bad2(iopbend) == 0) THEN
1145 0 : counter = counter + 1
1146 0 : new_opbend_list(counter) = opbend_list(iopbend)
1147 0 : newkind = opbend_list(iopbend)%opbend_kind%kind_number
1148 0 : newkind = newkind - SUM(bad1(1:newkind))
1149 0 : new_opbend_list(counter)%opbend_kind => new_opbend_kind_set(newkind)
1150 : END IF
1151 : END DO
1152 : CALL set_molecule_kind(molecule_kind=molecule_kind, &
1153 : nopbend=SIZE(new_opbend_list), &
1154 : opbend_kind_set=new_opbend_kind_set, &
1155 1676 : opbend_list=new_opbend_list)
1156 1676 : DO iopbend = 1, SIZE(new_opbend_kind_set)
1157 1676 : new_opbend_kind_set(iopbend)%kind_number = iopbend
1158 : END DO
1159 : END IF
1160 1676 : DEALLOCATE (bad2)
1161 1676 : DEALLOCATE (opbend_kind_set)
1162 1676 : DEALLOCATE (opbend_list)
1163 1698 : IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New OPBEND Count: ", &
1164 44 : SIZE(new_opbend_list), SIZE(new_opbend_kind_set)
1165 : END IF
1166 1678 : DEALLOCATE (bad1)
1167 : END IF
1168 : END DO
1169 : !---------------------------------------------------------------------------
1170 : ! 12. Count the number of UNSET NONBOND14 kinds there are
1171 : !- NEED TO REMOVE EXTRAS HERE - IKUO
1172 : !---------------------------------------------------------------------------
1173 : CALL cp_print_key_finished_output(iw, logger, mm_section, &
1174 2699 : "PRINT%FF_INFO")
1175 :
1176 2699 : CALL timestop(handle)
1177 :
1178 2699 : END SUBROUTINE clean_intra_force_kind
1179 :
1180 : END MODULE force_fields_util
|
