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