(git:374b731)
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-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! **************************************************************************************************
20
23 USE cell_types, ONLY: cell_type
30 USE cp_units, ONLY: cp_unit_to_cp2k
36 USE force_field_kind_types, ONLY: &
47 USE force_fields_all, ONLY: &
60 USE kinds, ONLY: default_path_length,&
62 dp
64 USE molecule_kind_types, ONLY: &
68 USE molecule_types, ONLY: get_molecule,&
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, &
86
87CONTAINS
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
1368END 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...
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)
...
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,...
unit conversion facility
Definition cp_units.F:30
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
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, quip_data, nequip_data, allegro_data, deepmd_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_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_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,...
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_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.
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_ub(particle_set, molecule_kind_set, molecule_set)
Determine the number of unique Urey-Bradley kind and allocate ub_kind_set.
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...
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_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.
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.
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_potentail%[qeff] and shell potential parameters.
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.
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.
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.
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_eicut(atomic_kind_set, ff_type, potparm_nonbond, ewald_env)
Compute the electrostatic interaction cutoffs.
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.
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 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.
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
integer, parameter, public default_path_length
Definition kinds.F:58
Utility routines for the memory handling.
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.
Utilities for string manipulations.
subroutine, public compress(string, full)
Eliminate multiple space characters in a string. If full is .TRUE., then all spaces are eliminated.
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
type of a logger, at the moment it contains just a print level starting at which level it should be l...