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