(git:0973b94)
Loading...
Searching...
No Matches
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-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! **************************************************************************************************
20
23 USE cell_types, ONLY: cell_type
35 USE force_field_kind_types, ONLY: &
46 USE force_fields_all, ONLY: &
59 USE kinds, ONLY: default_string_length,&
60 dp
61 USE molecule_kind_types, ONLY: &
65 USE molecule_types, ONLY: get_molecule,&
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, &
81
82CONTAINS
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 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 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 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 CALL timeset(routinen, handle)
139 fatal = .false.
140 ignore_fatal = ff_type%ignore_missing_critical
141 NULLIFY (logger, ainfo, charges_section, charges)
142 logger => cp_get_default_logger()
143 ! Error unit
144 output_unit = cp_logger_get_default_io_unit(logger)
145
146 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
147 extension=".mmLog")
148 iw2 = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO/SPLINE_INFO", &
149 extension=".mmLog")
150 iw3 = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO/SPLINE_DATA", &
151 extension=".mmLog")
152 iw4 = cp_print_key_unit_nr(logger, mm_section, "PRINT%PROGRAM_RUN_INFO", &
153 extension=".mmLog")
154 NULLIFY (potparm_nonbond14, potparm_nonbond)
155 my_qmmm = .false.
156 IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
157 inp_info => ff_type%inp_info
158 chm_info => ff_type%chm_info
159 gro_info => ff_type%gro_info
160 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 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 ff_type, iw)
171 !-----------------------------------------------------------------------------
172 ! 3. Determine the number of unique Urey-Bradley kind and allocate ub_kind_set
173 !-----------------------------------------------------------------------------
174 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 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 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 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 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 amb_info, iw)
202 ! Give information and abort if any bond or angle parameter is missing..
203 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 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 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 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 ainfo, inp_info, iw)
224 ! Give information only if any Urey-Bradley, Torsion, improper or opbend is missing
225 ! continue calculation
226 CALL release_ff_missing_par(fatal, ignore_fatal, ainfo, output_unit, iw)
227
228 charges_section => section_vals_get_subs_vals(mm_section, "FORCEFIELD%CHARGES")
229 CALL section_vals_get(charges_section, explicit=explicit)
230 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 ainfo, my_qmmm, inp_info)
237 ! Give information only if charge is missing and abort..
238 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 qmmm_env, inp_info, iw4)
246 END IF
247 !-----------------------------------------------------------------------------
248 ! 14. Set up the radius of the electrostatic multipole in Fist
249 !-----------------------------------------------------------------------------
250 CALL force_field_pack_radius(atomic_kind_set, iw, subsys_section)
251 !-----------------------------------------------------------------------------
252 ! 15. Set up the polarizable FF parameters
253 !-----------------------------------------------------------------------------
254 CALL force_field_pack_pol(atomic_kind_set, iw, inp_info)
255 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 shell_particle_set, core_particle_set, cell, iw, inp_info)
262 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 chm_info, inp_info, gro_info, amb_info, potparm_nonbond14, ewald_env)
269 ! Give information if any 1-4 is missing.. continue calculation..
270 CALL release_ff_missing_par(fatal, ignore_fatal, ainfo, output_unit, iw)
271 ! Create the spline data
272 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 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 potparm_nonbond, ewald_env)
282 ! Give information and abort if any pair potential spline is missing..
283 CALL release_ff_missing_par(fatal, ignore_fatal, ainfo, output_unit, iw)
284 ! Create the spline data
285 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 potparm_nonbond, do_zbl, nonbonded_type="NONBONDED")
288 END IF
289 !-----------------------------------------------------------------------------
290 ! 19. Create nonbond environment
291 !-----------------------------------------------------------------------------
292 CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
293 CALL section_vals_val_get(mm_section, "NEIGHBOR_LISTS%VERLET_SKIN", &
294 r_val=verlet_skin)
295 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 ff_type%vdw_scale14, ff_type%shift_cutoff)
300 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 ewald_env, iw)
304
305 CALL cp_print_key_finished_output(iw4, logger, mm_section, "PRINT%PROGRAM_RUN_INFO")
306 CALL cp_print_key_finished_output(iw3, logger, mm_section, "PRINT%FF_INFO/SPLINE_DATA")
307 CALL cp_print_key_finished_output(iw2, logger, mm_section, "PRINT%FF_INFO/SPLINE_INFO")
308 CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO")
309
310 CALL timestop(handle)
311
312 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 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 IF (ASSOCIATED(array)) THEN
331 IF (output_unit > 0) THEN
332 WRITE (output_unit, *)
333 WRITE (output_unit, '(T2,"FORCEFIELD|",A)') &
334 " WARNING: A non Critical ForceField parameter is missing! CP2K GOES!", &
335 " Non critical parameters are:Urey-Bradley,Torsions,Impropers, Opbends and 1-4", &
336 " All missing parameters will not contribute to the potential energy!"
337 IF (fatal .OR. iw > 0) THEN
338 WRITE (output_unit, *)
339 DO i = 1, SIZE(array)
340 WRITE (output_unit, '(A)') array(i)
341 END DO
342 END IF
343 IF (.NOT. fatal .AND. iw < 0) THEN
344 WRITE (output_unit, '(T2,"FORCEFIELD|",A)') &
345 " Activate the print key FF_INFO to have a list of missing parameters"
346 WRITE (output_unit, *)
347 END IF
348 END IF
349 DEALLOCATE (array)
350 END IF
351 IF (fatal) THEN
352 IF (ignore_fatal) THEN
353 IF (output_unit > 0) THEN
354 WRITE (output_unit, *)
355 WRITE (output_unit, '(T2,"FORCEFIELD|",A)') &
356 " WARNING: Ignoring missing critical FF parameters! CP2K GOES!", &
357 " Critical parameters are: Bonds, Bends and Charges", &
358 " All missing parameters will not contribute to the potential energy!"
359 END IF
360 ELSE
361 cpabort("Missing critical ForceField parameters!")
362 END IF
363 END IF
364 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 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 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 CALL timeset(routinen, handle)
398
399 NULLIFY (logger)
400 logger => cp_get_default_logger()
401 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
402 extension=".mmLog")
403
404 qeff = 0.0_dp
405 qeff_mol = 0.0_dp
406 qeff_sum = 0.0_dp
407 mass_sum = 0.0_dp
408 !-----------------------------------------------------------------------------
409 ! 1. Sum of [qeff,mass] for each molecule_kind
410 !-----------------------------------------------------------------------------
411 DO imol = 1, SIZE(molecule_kind_set)
412 qeff_mol = 0.0_dp
413 mass_mol = 0.0_dp
414 molecule_kind => molecule_kind_set(imol)
415
416 j = molecule_kind_set(imol)%molecule_list(1)
417 molecule => molecule_set(j)
418 CALL get_molecule(molecule=molecule, first_atom=first)
419
420 CALL get_molecule_kind(molecule_kind=molecule_kind, &
421 name=molname, atom_list=atom_list)
422 DO iatom = 1, SIZE(atom_list)
423 atomic_kind => atom_list(iatom)%atomic_kind
424 CALL get_atomic_kind(atomic_kind=atomic_kind, &
425 name=atmname, qeff=qeff, mass=mass, shell_active=shell_active, shell=shell)
426 IF (shell_active) qeff = shell%charge_core + shell%charge_shell
427 IF (ASSOCIATED(charges)) THEN
428 jatom = first - 1 + iatom
429 qeff = charges(jatom)
430 END IF
431 IF (iw > 0) WRITE (iw, *) " atom ", iatom, " ", trim(atmname), " charge = ", qeff
432 qeff_mol = qeff_mol + qeff
433 mass_mol = mass_mol + mass
434 END DO
435 CALL set_molecule_kind(molecule_kind=molecule_kind, charge=qeff_mol, mass=mass_mol)
436 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 DO iatom = 1, SIZE(particle_set)
442 atomic_kind => particle_set(iatom)%atomic_kind
443 CALL get_atomic_kind(atomic_kind=atomic_kind, &
444 name=atmname, qeff=qeff, mass=mass, shell_active=shell_active, shell=shell)
445 IF (shell_active) qeff = shell%charge_core + shell%charge_shell
446 IF (ASSOCIATED(charges)) THEN
447 qeff = charges(iatom)
448 END IF
449 IF (iw > 0) WRITE (iw, *) " atom ", iatom, " ", trim(atmname), &
450 " charge = ", qeff
451 qeff_sum = qeff_sum + qeff
452 mass_sum = mass_sum + mass
453 END DO
454 IF (iw > 0) WRITE (iw, '(A,F20.10)') " Total system charge = ", qeff_sum
455 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 "PRINT%FF_INFO")
459
460 CALL timestop(handle)
461
462 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 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 INTEGER, POINTER :: bad1(:), bad2(:)
480 LOGICAL :: unsetme, valid_kind
481 TYPE(bend_kind_type), DIMENSION(:), POINTER :: bend_kind_set, new_bend_kind_set
482 TYPE(bend_type), DIMENSION(:), POINTER :: bend_list, new_bend_list
483 TYPE(bond_kind_type), DIMENSION(:), POINTER :: bond_kind_set, new_bond_kind_set
484 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list, new_bond_list
485 TYPE(colvar_constraint_type), DIMENSION(:), &
486 POINTER :: colv_list
487 TYPE(cp_logger_type), POINTER :: logger
488 TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list
489 TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list
490 TYPE(impr_kind_type), DIMENSION(:), POINTER :: impr_kind_set, new_impr_kind_set
491 TYPE(impr_type), DIMENSION(:), POINTER :: impr_list, new_impr_list
492 TYPE(molecule_kind_type), POINTER :: molecule_kind
493 TYPE(opbend_kind_type), DIMENSION(:), POINTER :: new_opbend_kind_set, opbend_kind_set
494 TYPE(opbend_type), DIMENSION(:), POINTER :: new_opbend_list, opbend_list
495 TYPE(torsion_kind_type), DIMENSION(:), POINTER :: new_torsion_kind_set, torsion_kind_set
496 TYPE(torsion_type), DIMENSION(:), POINTER :: new_torsion_list, torsion_list
497 TYPE(ub_kind_type), DIMENSION(:), POINTER :: new_ub_kind_set, ub_kind_set
498 TYPE(ub_type), DIMENSION(:), POINTER :: new_ub_list, ub_list
499
500 CALL timeset(routinen, handle)
501 NULLIFY (logger)
502 logger => cp_get_default_logger()
503 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", &
504 extension=".mmLog")
505 !-----------------------------------------------------------------------------
506 ! 1. Lets Tag the unwanted bonds due to the use of distance constraint
507 !-----------------------------------------------------------------------------
508 DO ikind = 1, SIZE(molecule_kind_set)
509 molecule_kind => molecule_kind_set(ikind)
510 CALL get_molecule_kind(molecule_kind=molecule_kind, &
511 colv_list=colv_list, &
512 nbond=nbond, &
513 bond_list=bond_list)
514 IF (ASSOCIATED(colv_list)) THEN
515 DO icolv = 1, SIZE(colv_list)
516 IF ((colv_list(icolv)%type_id == dist_colvar_id) .AND. &
517 ((.NOT. colv_list(icolv)%use_points) .OR. (SIZE(colv_list(icolv)%i_atoms) == 2))) THEN
518 atm_a = colv_list(icolv)%i_atoms(1)
519 atm_b = colv_list(icolv)%i_atoms(2)
520 DO ibond = 1, nbond
521 unsetme = .false.
522 atm2_a = bond_list(ibond)%a
523 atm2_b = bond_list(ibond)%b
524 IF (atm2_a == atm_a .AND. atm2_b == atm_b) unsetme = .true.
525 IF (atm2_a == atm_b .AND. atm2_b == atm_a) unsetme = .true.
526 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 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 nbend=nbend, &
540 bend_list=bend_list)
541 IF (ASSOCIATED(colv_list)) THEN
542 DO ibend = 1, nbend
543 unsetme = .false.
544 atm_a = bend_list(ibend)%a
545 atm_b = bend_list(ibend)%b
546 atm_c = bend_list(ibend)%c
547 DO icolv = 1, SIZE(colv_list)
548 IF ((colv_list(icolv)%type_id == dist_colvar_id) .AND. &
549 ((.NOT. colv_list(icolv)%use_points) .OR. (SIZE(colv_list(icolv)%i_atoms) == 2))) THEN
550 atm2_a = colv_list(icolv)%i_atoms(1)
551 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 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 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 DO ikind = 1, SIZE(molecule_kind_set)
569 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 bond_list=bond_list)
575 DO ig3x3 = 1, ng3x3
576 atm_a = g3x3_list(ig3x3)%a
577 atm_b = g3x3_list(ig3x3)%b
578 atm_c = g3x3_list(ig3x3)%c
579 DO ibond = 1, nbond
580 unsetme = .false.
581 atm2_a = bond_list(ibond)%a
582 atm2_b = bond_list(ibond)%b
583 IF (atm2_a == atm_a .AND. atm2_b == atm_b) unsetme = .true.
584 IF (atm2_a == atm_a .AND. atm2_b == atm_c) unsetme = .true.
585 IF (atm2_a == atm_c .AND. atm2_b == atm_c) unsetme = .true.
586 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 DO ikind = 1, SIZE(molecule_kind_set)
594 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 bend_list=bend_list)
600 DO ig3x3 = 1, ng3x3
601 atm_a = g3x3_list(ig3x3)%a
602 atm_b = g3x3_list(ig3x3)%b
603 atm_c = g3x3_list(ig3x3)%c
604 DO ibend = 1, nbend
605 unsetme = .false.
606 atm2_a = bend_list(ibend)%a
607 atm2_b = bend_list(ibend)%b
608 atm2_c = bend_list(ibend)%c
609 IF (atm2_a == atm_a .AND. atm2_b == atm_b .AND. atm2_c == atm_c) unsetme = .true.
610 IF (atm2_a == atm_a .AND. atm2_b == atm_c .AND. atm2_c == atm_b) unsetme = .true.
611 IF (atm2_a == atm_b .AND. atm2_b == atm_a .AND. atm2_c == atm_c) unsetme = .true.
612 IF (atm2_a == atm_b .AND. atm2_b == atm_c .AND. atm2_c == atm_a) unsetme = .true.
613 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 DO ikind = 1, SIZE(molecule_kind_set)
621 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 bond_list=bond_list)
627 DO ig4x6 = 1, ng4x6
628 atm_a = g4x6_list(ig4x6)%a
629 atm_b = g4x6_list(ig4x6)%b
630 atm_c = g4x6_list(ig4x6)%c
631 atm_d = g4x6_list(ig4x6)%d
632 DO ibond = 1, nbond
633 unsetme = .false.
634 atm2_a = bond_list(ibond)%a
635 atm2_b = bond_list(ibond)%b
636 IF (atm2_a == atm_a .AND. atm2_b == atm_b) unsetme = .true.
637 IF (atm2_a == atm_a .AND. atm2_b == atm_c) unsetme = .true.
638 IF (atm2_a == atm_a .AND. atm2_b == atm_d) unsetme = .true.
639 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 DO ikind = 1, SIZE(molecule_kind_set)
647 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 bend_list=bend_list)
653 DO ig4x6 = 1, ng4x6
654 atm_a = g4x6_list(ig4x6)%a
655 atm_b = g4x6_list(ig4x6)%b
656 atm_c = g4x6_list(ig4x6)%c
657 atm_d = g4x6_list(ig4x6)%d
658 DO ibend = 1, nbend
659 unsetme = .false.
660 atm2_a = bend_list(ibend)%a
661 atm2_b = bend_list(ibend)%b
662 atm2_c = bend_list(ibend)%c
663 IF (atm2_a == atm_b .AND. atm2_b == atm_a .AND. atm2_c == atm_c) unsetme = .true.
664 IF (atm2_a == atm_b .AND. atm2_b == atm_a .AND. atm2_c == atm_d) unsetme = .true.
665 IF (atm2_a == atm_c .AND. atm2_b == atm_a .AND. atm2_c == atm_d) unsetme = .true.
666 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 DO ikind = 1, SIZE(molecule_kind_set)
674 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 bond_list=bond_list)
679 IF (nbond > 0) THEN
680 IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old BOND Count: ", &
681 SIZE(bond_list), SIZE(bond_kind_set)
682 IF (iw > 0) WRITE (iw, '(2I6)') (bond_list(ibond)%a, bond_list(ibond)%b, ibond=1, SIZE(bond_list))
683 NULLIFY (bad1, bad2)
684 ALLOCATE (bad1(SIZE(bond_kind_set)))
685 bad1(:) = 0
686 DO ibond = 1, SIZE(bond_kind_set)
687 unsetme = .false.
688 IF (bond_kind_set(ibond)%id_type == do_ff_undef) unsetme = .true.
689 valid_kind = .false.
690 DO i = 1, SIZE(bond_list)
691 IF (bond_list(i)%id_type /= do_ff_undef .AND. &
692 bond_list(i)%bond_kind%kind_number == ibond) THEN
693 valid_kind = .true.
694 EXIT
695 END IF
696 END DO
697 IF (.NOT. valid_kind) unsetme = .true.
698 IF (unsetme) bad1(ibond) = 1
699 END DO
700 IF (sum(bad1) /= 0) THEN
701 counter = SIZE(bond_kind_set) - sum(bad1)
702 CALL allocate_bond_kind_set(new_bond_kind_set, counter)
703 counter = 0
704 DO ibond = 1, SIZE(bond_kind_set)
705 IF (bad1(ibond) == 0) THEN
706 counter = counter + 1
707 new_bond_kind_set(counter) = bond_kind_set(ibond)
708 END IF
709 END DO
710 counter = 0
711 ALLOCATE (bad2(SIZE(bond_list)))
712 bad2(:) = 0
713 DO ibond = 1, SIZE(bond_list)
714 unsetme = .false.
715 IF (bond_list(ibond)%bond_kind%id_type == do_ff_undef) unsetme = .true.
716 IF (bond_list(ibond)%id_type == do_ff_undef) unsetme = .true.
717 IF (unsetme) bad2(ibond) = 1
718 END DO
719 IF (sum(bad2) /= 0) THEN
720 counter = SIZE(bond_list) - sum(bad2)
721 ALLOCATE (new_bond_list(counter))
722 counter = 0
723 DO ibond = 1, SIZE(bond_list)
724 IF (bad2(ibond) == 0) THEN
725 counter = counter + 1
726 new_bond_list(counter) = bond_list(ibond)
727 newkind = bond_list(ibond)%bond_kind%kind_number
728 newkind = newkind - sum(bad1(1:newkind))
729 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 bond_list=new_bond_list)
736 DO ibond = 1, SIZE(new_bond_kind_set)
737 new_bond_kind_set(ibond)%kind_number = ibond
738 END DO
739 END IF
740 DEALLOCATE (bad2)
741 CALL deallocate_bond_kind_set(bond_kind_set)
742 DEALLOCATE (bond_list)
743 IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New BOND Count: ", &
744 SIZE(new_bond_list), SIZE(new_bond_kind_set)
745 IF (iw > 0) WRITE (iw, '(2I6)') (new_bond_list(ibond)%a, new_bond_list(ibond)%b, &
746 ibond=1, SIZE(new_bond_list))
747 END IF
748 DEALLOCATE (bad1)
749 END IF
750 END DO
751 !-----------------------------------------------------------------------------
752 ! 8. Count the number of UNSET bend kinds there are
753 !-----------------------------------------------------------------------------
754 DO ikind = 1, SIZE(molecule_kind_set)
755 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 bend_list=bend_list)
760 IF (nbend > 0) THEN
761 IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old BEND Count: ", &
762 SIZE(bend_list), SIZE(bend_kind_set)
763 IF (iw > 0) WRITE (iw, '(3I6)') (bend_list(ibend)%a, bend_list(ibend)%b, &
764 bend_list(ibend)%c, ibend=1, SIZE(bend_list))
765 NULLIFY (bad1, bad2)
766 ALLOCATE (bad1(SIZE(bend_kind_set)))
767 bad1(:) = 0
768 DO ibend = 1, SIZE(bend_kind_set)
769 unsetme = .false.
770 IF (bend_kind_set(ibend)%id_type == do_ff_undef) unsetme = .true.
771 valid_kind = .false.
772 DO i = 1, SIZE(bend_list)
773 IF (bend_list(i)%id_type /= do_ff_undef .AND. &
774 bend_list(i)%bend_kind%kind_number == ibend) THEN
775 valid_kind = .true.
776 EXIT
777 END IF
778 END DO
779 IF (.NOT. valid_kind) unsetme = .true.
780 IF (unsetme) bad1(ibend) = 1
781 END DO
782 IF (sum(bad1) /= 0) THEN
783 counter = SIZE(bend_kind_set) - sum(bad1)
784 CALL allocate_bend_kind_set(new_bend_kind_set, counter)
785 counter = 0
786 DO ibend = 1, SIZE(bend_kind_set)
787 IF (bad1(ibend) == 0) THEN
788 counter = counter + 1
789 new_bend_kind_set(counter) = bend_kind_set(ibend)
790 END IF
791 END DO
792 counter = 0
793 ALLOCATE (bad2(SIZE(bend_list)))
794 bad2(:) = 0
795 DO ibend = 1, SIZE(bend_list)
796 unsetme = .false.
797 IF (bend_list(ibend)%bend_kind%id_type == do_ff_undef) unsetme = .true.
798 IF (bend_list(ibend)%id_type == do_ff_undef) unsetme = .true.
799 IF (unsetme) bad2(ibend) = 1
800 END DO
801 IF (sum(bad2) /= 0) THEN
802 counter = SIZE(bend_list) - sum(bad2)
803 ALLOCATE (new_bend_list(counter))
804 counter = 0
805 DO ibend = 1, SIZE(bend_list)
806 IF (bad2(ibend) == 0) THEN
807 counter = counter + 1
808 new_bend_list(counter) = bend_list(ibend)
809 newkind = bend_list(ibend)%bend_kind%kind_number
810 newkind = newkind - sum(bad1(1:newkind))
811 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 bend_list=new_bend_list)
818 DO ibend = 1, SIZE(new_bend_kind_set)
819 new_bend_kind_set(ibend)%kind_number = ibend
820 END DO
821 END IF
822 DEALLOCATE (bad2)
823 CALL deallocate_bend_kind_set(bend_kind_set)
824 DEALLOCATE (bend_list)
825 IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New BEND Count: ", &
826 SIZE(new_bend_list), SIZE(new_bend_kind_set)
827 IF (iw > 0) WRITE (iw, '(3I6)') (new_bend_list(ibend)%a, new_bend_list(ibend)%b, &
828 new_bend_list(ibend)%c, ibend=1, SIZE(new_bend_list))
829 END IF
830 DEALLOCATE (bad1)
831 END IF
832 END DO
833 !-----------------------------------------------------------------------------
834 ! 9. Count the number of UNSET Urey-Bradley kinds there are
835 !-----------------------------------------------------------------------------
836 DO ikind = 1, SIZE(molecule_kind_set)
837 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 ub_list=ub_list)
842 IF (nub > 0) THEN
843 IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old UB Count: ", &
844 SIZE(ub_list), SIZE(ub_kind_set)
845 IF (iw > 0) WRITE (iw, '(3I6)') (ub_list(iub)%a, ub_list(iub)%b, &
846 ub_list(iub)%c, iub=1, SIZE(ub_list))
847 NULLIFY (bad1, bad2)
848 ALLOCATE (bad1(SIZE(ub_kind_set)))
849 bad1(:) = 0
850 DO iub = 1, SIZE(ub_kind_set)
851 unsetme = .false.
852 IF (ub_kind_set(iub)%id_type == do_ff_undef) unsetme = .true.
853 valid_kind = .false.
854 DO i = 1, SIZE(ub_list)
855 IF (ub_list(i)%id_type /= do_ff_undef .AND. &
856 ub_list(i)%ub_kind%kind_number == iub) THEN
857 valid_kind = .true.
858 EXIT
859 END IF
860 END DO
861 IF (.NOT. valid_kind) unsetme = .true.
862 IF (unsetme) bad1(iub) = 1
863 END DO
864 IF (sum(bad1) /= 0) THEN
865 counter = SIZE(ub_kind_set) - sum(bad1)
866 CALL allocate_ub_kind_set(new_ub_kind_set, counter)
867 counter = 0
868 DO iub = 1, SIZE(ub_kind_set)
869 IF (bad1(iub) == 0) THEN
870 counter = counter + 1
871 new_ub_kind_set(counter) = ub_kind_set(iub)
872 END IF
873 END DO
874 counter = 0
875 ALLOCATE (bad2(SIZE(ub_list)))
876 bad2(:) = 0
877 DO iub = 1, SIZE(ub_list)
878 unsetme = .false.
879 IF (ub_list(iub)%ub_kind%id_type == do_ff_undef) unsetme = .true.
880 IF (ub_list(iub)%id_type == do_ff_undef) unsetme = .true.
881 IF (unsetme) bad2(iub) = 1
882 END DO
883 IF (sum(bad2) /= 0) THEN
884 counter = SIZE(ub_list) - sum(bad2)
885 ALLOCATE (new_ub_list(counter))
886 counter = 0
887 DO iub = 1, SIZE(ub_list)
888 IF (bad2(iub) == 0) THEN
889 counter = counter + 1
890 new_ub_list(counter) = ub_list(iub)
891 newkind = ub_list(iub)%ub_kind%kind_number
892 newkind = newkind - sum(bad1(1:newkind))
893 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 ub_list=new_ub_list)
900 DO iub = 1, SIZE(new_ub_kind_set)
901 new_ub_kind_set(iub)%kind_number = iub
902 END DO
903 END IF
904 DEALLOCATE (bad2)
905 CALL ub_kind_dealloc_ref(ub_kind_set)
906 DEALLOCATE (ub_list)
907 IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New UB Count: ", &
908 SIZE(new_ub_list), SIZE(new_ub_kind_set)
909 IF (iw > 0) WRITE (iw, '(3I6)') (new_ub_list(iub)%a, new_ub_list(iub)%b, &
910 new_ub_list(iub)%c, iub=1, SIZE(new_ub_list))
911 END IF
912 DEALLOCATE (bad1)
913 END IF
914 END DO
915 !-----------------------------------------------------------------------------
916 ! 10. Count the number of UNSET torsion kinds there are
917 !-----------------------------------------------------------------------------
918 DO ikind = 1, SIZE(molecule_kind_set)
919 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 torsion_list=torsion_list)
924 IF (ntorsion > 0) THEN
925 IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old TORSION Count: ", &
926 SIZE(torsion_list), SIZE(torsion_kind_set)
927 IF (iw > 0) WRITE (iw, '(4I6)') (torsion_list(itorsion)%a, torsion_list(itorsion)%b, &
928 torsion_list(itorsion)%c, torsion_list(itorsion)%d, itorsion=1, SIZE(torsion_list))
929 NULLIFY (bad1, bad2)
930 ALLOCATE (bad1(SIZE(torsion_kind_set)))
931 bad1(:) = 0
932 DO itorsion = 1, SIZE(torsion_kind_set)
933 unsetme = .false.
934 IF (torsion_kind_set(itorsion)%id_type == do_ff_undef) unsetme = .true.
935 valid_kind = .false.
936 DO i = 1, SIZE(torsion_list)
937 IF (torsion_list(i)%id_type /= do_ff_undef .AND. &
938 torsion_list(i)%torsion_kind%kind_number == itorsion) THEN
939 valid_kind = .true.
940 EXIT
941 END IF
942 END DO
943 IF (.NOT. valid_kind) unsetme = .true.
944 IF (unsetme) bad1(itorsion) = 1
945 END DO
946 IF (sum(bad1) /= 0) THEN
947 counter = SIZE(torsion_kind_set) - sum(bad1)
948 CALL allocate_torsion_kind_set(new_torsion_kind_set, counter)
949 counter = 0
950 DO itorsion = 1, SIZE(torsion_kind_set)
951 IF (bad1(itorsion) == 0) THEN
952 counter = counter + 1
953 new_torsion_kind_set(counter) = torsion_kind_set(itorsion)
954 i = SIZE(torsion_kind_set(itorsion)%m)
955 j = SIZE(torsion_kind_set(itorsion)%k)
956 k = SIZE(torsion_kind_set(itorsion)%phi0)
957 ALLOCATE (new_torsion_kind_set(counter)%m(i))
958 ALLOCATE (new_torsion_kind_set(counter)%k(i))
959 ALLOCATE (new_torsion_kind_set(counter)%phi0(i))
960 new_torsion_kind_set(counter)%m = torsion_kind_set(itorsion)%m
961 new_torsion_kind_set(counter)%k = torsion_kind_set(itorsion)%k
962 new_torsion_kind_set(counter)%phi0 = torsion_kind_set(itorsion)%phi0
963 END IF
964 END DO
965 counter = 0
966 ALLOCATE (bad2(SIZE(torsion_list)))
967 bad2(:) = 0
968 DO itorsion = 1, SIZE(torsion_list)
969 unsetme = .false.
970 IF (torsion_list(itorsion)%torsion_kind%id_type == do_ff_undef) unsetme = .true.
971 IF (torsion_list(itorsion)%id_type == do_ff_undef) unsetme = .true.
972 IF (unsetme) bad2(itorsion) = 1
973 END DO
974 IF (sum(bad2) /= 0) THEN
975 counter = SIZE(torsion_list) - sum(bad2)
976 ALLOCATE (new_torsion_list(counter))
977 counter = 0
978 DO itorsion = 1, SIZE(torsion_list)
979 IF (bad2(itorsion) == 0) THEN
980 counter = counter + 1
981 new_torsion_list(counter) = torsion_list(itorsion)
982 newkind = torsion_list(itorsion)%torsion_kind%kind_number
983 newkind = newkind - sum(bad1(1:newkind))
984 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 torsion_list=new_torsion_list)
991 DO itorsion = 1, SIZE(new_torsion_kind_set)
992 new_torsion_kind_set(itorsion)%kind_number = itorsion
993 END DO
994 END IF
995 DEALLOCATE (bad2)
996 DO itorsion = 1, SIZE(torsion_kind_set)
997 CALL torsion_kind_dealloc_ref(torsion_kind_set(itorsion))
998 END DO
999 DEALLOCATE (torsion_kind_set)
1000 DEALLOCATE (torsion_list)
1001 IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New TORSION Count: ", &
1002 SIZE(new_torsion_list), SIZE(new_torsion_kind_set)
1003 IF (iw > 0) WRITE (iw, '(4I6)') (new_torsion_list(itorsion)%a, new_torsion_list(itorsion)%b, &
1004 new_torsion_list(itorsion)%c, new_torsion_list(itorsion)%d, itorsion=1, &
1005 SIZE(new_torsion_list))
1006 END IF
1007 DEALLOCATE (bad1)
1008 END IF
1009 END DO
1010 !-----------------------------------------------------------------------------
1011 ! 11. Count the number of UNSET improper kinds there are
1012 !-----------------------------------------------------------------------------
1013 DO ikind = 1, SIZE(molecule_kind_set)
1014 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 impr_list=impr_list)
1019 IF (nimpr > 0) THEN
1020 IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old IMPROPER Count: ", &
1021 SIZE(impr_list), SIZE(impr_kind_set)
1022 NULLIFY (bad1, bad2)
1023 ALLOCATE (bad1(SIZE(impr_kind_set)))
1024 bad1(:) = 0
1025 DO iimpr = 1, SIZE(impr_kind_set)
1026 unsetme = .false.
1027 IF (impr_kind_set(iimpr)%id_type == do_ff_undef) unsetme = .true.
1028 valid_kind = .false.
1029 DO i = 1, SIZE(impr_list)
1030 IF (impr_list(i)%id_type /= do_ff_undef .AND. &
1031 impr_list(i)%impr_kind%kind_number == iimpr) THEN
1032 valid_kind = .true.
1033 EXIT
1034 END IF
1035 END DO
1036 IF (.NOT. valid_kind) unsetme = .true.
1037 IF (unsetme) bad1(iimpr) = 1
1038 END DO
1039 IF (sum(bad1) /= 0) THEN
1040 counter = SIZE(impr_kind_set) - sum(bad1)
1041 CALL allocate_impr_kind_set(new_impr_kind_set, counter)
1042 counter = 0
1043 DO iimpr = 1, SIZE(impr_kind_set)
1044 IF (bad1(iimpr) == 0) THEN
1045 counter = counter + 1
1046 new_impr_kind_set(counter) = impr_kind_set(iimpr)
1047 END IF
1048 END DO
1049 counter = 0
1050 ALLOCATE (bad2(SIZE(impr_list)))
1051 bad2(:) = 0
1052 DO iimpr = 1, SIZE(impr_list)
1053 unsetme = .false.
1054 IF (impr_list(iimpr)%impr_kind%id_type == do_ff_undef) unsetme = .true.
1055 IF (impr_list(iimpr)%id_type == do_ff_undef) unsetme = .true.
1056 IF (unsetme) bad2(iimpr) = 1
1057 END DO
1058 IF (sum(bad2) /= 0) THEN
1059 counter = SIZE(impr_list) - sum(bad2)
1060 ALLOCATE (new_impr_list(counter))
1061 counter = 0
1062 DO iimpr = 1, SIZE(impr_list)
1063 IF (bad2(iimpr) == 0) THEN
1064 counter = counter + 1
1065 new_impr_list(counter) = impr_list(iimpr)
1066 newkind = impr_list(iimpr)%impr_kind%kind_number
1067 newkind = newkind - sum(bad1(1:newkind))
1068 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 impr_list=new_impr_list)
1075 DO iimpr = 1, SIZE(new_impr_kind_set)
1076 new_impr_kind_set(iimpr)%kind_number = iimpr
1077 END DO
1078 END IF
1079 DEALLOCATE (bad2)
1080 DO iimpr = 1, SIZE(impr_kind_set)
1081 CALL impr_kind_dealloc_ref() !This Subroutine doesn't deallocate anything, maybe needs to be implemented
1082 END DO
1083 DEALLOCATE (impr_kind_set)
1084 DEALLOCATE (impr_list)
1085 IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New IMPROPER Count: ", &
1086 SIZE(new_impr_list), SIZE(new_impr_kind_set)
1087 END IF
1088 DEALLOCATE (bad1)
1089 END IF
1090 END DO
1091 !-----------------------------------------------------------------------------
1092 ! 11. Count the number of UNSET opbends kinds there are
1093 !-----------------------------------------------------------------------------
1094 DO ikind = 1, SIZE(molecule_kind_set)
1095 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 opbend_list=opbend_list)
1100 IF (nopbend > 0) THEN
1101 IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") Old OPBEND Count: ", &
1102 SIZE(opbend_list), SIZE(opbend_kind_set)
1103 NULLIFY (bad1, bad2)
1104 ALLOCATE (bad1(SIZE(opbend_kind_set)))
1105 bad1(:) = 0
1106 DO iopbend = 1, SIZE(opbend_kind_set)
1107 unsetme = .false.
1108 IF (opbend_kind_set(iopbend)%id_type == do_ff_undef) unsetme = .true.
1109 valid_kind = .false.
1110 DO i = 1, SIZE(opbend_list)
1111 IF (opbend_list(i)%id_type /= do_ff_undef .AND. &
1112 opbend_list(i)%opbend_kind%kind_number == iopbend) THEN
1113 valid_kind = .true.
1114 EXIT
1115 END IF
1116 END DO
1117 IF (.NOT. valid_kind) unsetme = .true.
1118 IF (unsetme) bad1(iopbend) = 1
1119 END DO
1120 IF (sum(bad1) /= 0) THEN
1121 counter = SIZE(opbend_kind_set) - sum(bad1)
1122 CALL allocate_opbend_kind_set(new_opbend_kind_set, counter)
1123 counter = 0
1124 DO iopbend = 1, SIZE(opbend_kind_set)
1125 IF (bad1(iopbend) == 0) THEN
1126 counter = counter + 1
1127 new_opbend_kind_set(counter) = opbend_kind_set(iopbend)
1128 END IF
1129 END DO
1130 counter = 0
1131 ALLOCATE (bad2(SIZE(opbend_list)))
1132 bad2(:) = 0
1133 DO iopbend = 1, SIZE(opbend_list)
1134 unsetme = .false.
1135 IF (opbend_list(iopbend)%opbend_kind%id_type == do_ff_undef) unsetme = .true.
1136 IF (opbend_list(iopbend)%id_type == do_ff_undef) unsetme = .true.
1137 IF (unsetme) bad2(iopbend) = 1
1138 END DO
1139 IF (sum(bad2) /= 0) THEN
1140 counter = SIZE(opbend_list) - sum(bad2)
1141 ALLOCATE (new_opbend_list(counter))
1142 counter = 0
1143 DO iopbend = 1, SIZE(opbend_list)
1144 IF (bad2(iopbend) == 0) THEN
1145 counter = counter + 1
1146 new_opbend_list(counter) = opbend_list(iopbend)
1147 newkind = opbend_list(iopbend)%opbend_kind%kind_number
1148 newkind = newkind - sum(bad1(1:newkind))
1149 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 opbend_list=new_opbend_list)
1156 DO iopbend = 1, SIZE(new_opbend_kind_set)
1157 new_opbend_kind_set(iopbend)%kind_number = iopbend
1158 END DO
1159 END IF
1160 DEALLOCATE (bad2)
1161 DEALLOCATE (opbend_kind_set)
1162 DEALLOCATE (opbend_list)
1163 IF (iw > 0) WRITE (iw, *) " Mol(", ikind, ") New OPBEND Count: ", &
1164 SIZE(new_opbend_list), SIZE(new_opbend_kind_set)
1165 END IF
1166 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 "PRINT%FF_INFO")
1175
1176 CALL timestop(handle)
1177
1178 END SUBROUTINE clean_intra_force_kind
1179
1180END MODULE force_fields_util
Define the atomic kind types and their sub types.
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.
Handles all functions related to the CELL.
Definition cell_types.F:15
Initialize the collective variables types.
integer, parameter, public dist_colvar_id
various routines to log and control the output. The idea is that decisions about where to log should ...
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...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
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,...
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, file_access)
...
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.
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
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, nequip_data, deepmd_data, ace_data, charges)
sets a fist_nonbond_env
Define all structure types related to force field kinds.
pure subroutine, public deallocate_bend_kind_set(bend_kind_set)
Deallocate a bend kind set.
integer, parameter, public do_ff_undef
pure subroutine, public allocate_impr_kind_set(impr_kind_set, nkind)
Allocate and initialize a impr kind set.
pure subroutine, public torsion_kind_dealloc_ref(torsion_kind)
Deallocate a torsion kind element.
pure subroutine, public allocate_torsion_kind_set(torsion_kind_set, nkind)
Allocate and initialize a torsion kind set.
pure subroutine, public allocate_bond_kind_set(bond_kind_set, nkind)
Allocate and initialize a bond kind set.
pure subroutine, public allocate_opbend_kind_set(opbend_kind_set, nkind)
Allocate and initialize a opbend kind set.
pure subroutine, public ub_kind_dealloc_ref(ub_kind_set)
Deallocate a ub kind set.
pure subroutine, public allocate_bend_kind_set(bend_kind_set, nkind)
Allocate and initialize a bend kind set.
pure subroutine, public deallocate_bond_kind_set(bond_kind_set)
Deallocate a bond kind set.
pure subroutine, public allocate_ub_kind_set(ub_kind_set, nkind)
Allocate and initialize a ub kind set.
pure subroutine, public impr_kind_dealloc_ref()
Deallocate a impr kind element.
Define all structures types related to force_fields.
subroutine, public force_field_unique_ub(particle_set, molecule_kind_set, molecule_set, iw)
Determine the number of unique Urey-Bradley kind and allocate ub_kind_set.
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.
subroutine, public force_field_unique_bond(particle_set, molecule_kind_set, molecule_set, ff_type, iw)
Determine the number of unique bond kind and allocate bond_kind_set.
subroutine, public force_field_unique_opbend(particle_set, molecule_kind_set, molecule_set, ff_type, iw)
Determine the number of unique opbend kind and allocate opbend_kind_set based on the present improper...
subroutine, public force_field_pack_impr(particle_set, molecule_kind_set, molecule_set, ainfo, chm_info, inp_info, gro_info, iw)
Pack in impropers information needed for the force_field.
subroutine, public force_field_pack_radius(atomic_kind_set, iw, subsys_section)
Set up the radius of the electrostatic multipole in Fist.
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.
subroutine, public force_field_pack_opbend(particle_set, molecule_kind_set, molecule_set, ainfo, inp_info, iw)
Pack in opbend information needed for the force_field. No loop over params for charmm,...
subroutine, public force_field_pack_charges(charges, charges_section, particle_set, my_qmmm, qmmm_env, inp_info, iw4)
Set up array of full charges.
subroutine, public force_field_pack_pol(atomic_kind_set, iw, inp_info)
Set up the polarizable FF parameters.
subroutine, public force_field_unique_bend(particle_set, molecule_kind_set, molecule_set, ff_type, iw)
Determine the number of unique bend kind and allocate bend_kind_set.
subroutine, public force_field_unique_impr(particle_set, molecule_kind_set, molecule_set, ff_type, iw)
Determine the number of unique impr kind and allocate impr_kind_set.
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
subroutine, public force_field_unique_tors(particle_set, molecule_kind_set, molecule_set, ff_type, iw)
Determine the number of unique torsion kind and allocate torsion_kind_set.
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.
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_potential%[qeff] and shell potential parameters.
subroutine, public force_field_pack_bend(particle_set, molecule_kind_set, molecule_set, fatal, ainfo, chm_info, inp_info, gro_info, amb_info, iw)
Pack in bends information needed for the force_field.
subroutine, public force_field_pack_eicut(atomic_kind_set, ff_type, potparm_nonbond, ewald_env, iw)
Compute the electrostatic interaction cutoffs.
subroutine, public force_field_pack_bond(particle_set, molecule_kind_set, molecule_set, fatal, ainfo, chm_info, inp_info, gro_info, amb_info, iw)
Pack in bonds information needed for the force_field.
subroutine, public force_field_pack_damp(atomic_kind_set, iw, inp_info)
Set up damping parameters.
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.
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.
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.
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.
subroutine, public clean_intra_force_kind(molecule_kind_set, mm_section)
Removes UNSET force field types.
objects that represent the structure of input sections and the data contained in an input section
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
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
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
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Define the molecule kind structure types and the corresponding functionality.
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.
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.
Define the data structure for the molecule information.
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.
Define the data structure for the particle information.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
type of a logger, at the moment it contains just a print level starting at which level it should be l...