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!> Splitting and cleaning the original force_field_pack - May 2007 9!> Teodoro Laino - Zurich University 10!> \author CJM 11! ************************************************************************************************** 12MODULE force_fields_all 13 14 USE atomic_kind_types, ONLY: atomic_kind_type,& 15 get_atomic_kind,& 16 get_atomic_kind_set,& 17 set_atomic_kind 18 USE atoms_input, ONLY: read_shell_coord_input 19 USE cell_types, ONLY: cell_type 20 USE cp_linked_list_input, ONLY: cp_sll_val_next,& 21 cp_sll_val_type 22 USE cp_log_handling, ONLY: cp_to_string 23 USE damping_dipole_types, ONLY: damping_p_create,& 24 damping_p_type,& 25 tang_toennies 26 USE ewald_environment_types, ONLY: ewald_env_get,& 27 ewald_env_set,& 28 ewald_environment_type 29 USE external_potential_types, ONLY: fist_potential_type,& 30 get_potential,& 31 set_potential 32 USE force_field_kind_types, ONLY: & 33 allocate_bend_kind_set, allocate_bond_kind_set, allocate_impr_kind_set, & 34 allocate_opbend_kind_set, allocate_torsion_kind_set, allocate_ub_kind_set, bend_kind_type, & 35 bond_kind_type, do_ff_amber, do_ff_charmm, do_ff_g87, do_ff_g96, do_ff_undef, & 36 impr_kind_type, opbend_kind_type, torsion_kind_type, ub_kind_type 37 USE force_field_types, ONLY: amber_info_type,& 38 charmm_info_type,& 39 force_field_type,& 40 gromos_info_type,& 41 input_info_type 42 USE input_constants, ONLY: do_qmmm_none 43 USE input_cp2k_binary_restarts, ONLY: read_binary_cs_coordinates 44 USE input_section_types, ONLY: section_vals_get,& 45 section_vals_get_subs_vals,& 46 section_vals_list_get,& 47 section_vals_type,& 48 section_vals_val_get 49 USE input_val_types, ONLY: val_get,& 50 val_type 51 USE kinds, ONLY: default_path_length,& 52 default_string_length,& 53 dp 54 USE mathconstants, ONLY: sqrthalf 55 USE memory_utilities, ONLY: reallocate 56 USE molecule_kind_types, ONLY: & 57 bend_type, bond_type, get_molecule_kind, impr_type, molecule_kind_type, opbend_type, & 58 set_molecule_kind, shell_type, torsion_type, ub_type 59 USE molecule_types, ONLY: get_molecule,& 60 molecule_type 61 USE pair_potential, ONLY: get_nonbond_storage,& 62 spline_nonbond_control 63 USE pair_potential_coulomb, ONLY: potential_coulomb 64 USE pair_potential_types, ONLY: & 65 ea_type, lj_charmm_type, lj_type, nn_type, nosh_nosh, nosh_sh, pair_potential_lj_create, & 66 pair_potential_pp_create, pair_potential_pp_type, pair_potential_single_add, & 67 pair_potential_single_clean, pair_potential_single_copy, pair_potential_single_type, & 68 quip_type, sh_sh, siepmann_type, tersoff_type 69 USE particle_types, ONLY: allocate_particle_set,& 70 particle_type 71 USE physcon, ONLY: bohr 72 USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm 73 USE qmmm_types_low, ONLY: qmmm_env_mm_type 74 USE shell_potential_types, ONLY: shell_create,& 75 shell_kind_type,& 76 shell_release 77 USE splines_types, ONLY: spline_data_p_release,& 78 spline_data_p_retain,& 79 spline_data_p_type,& 80 spline_env_release,& 81 spline_environment_type 82 USE string_utilities, ONLY: compress,& 83 integer_to_string,& 84 uppercase 85#include "./base/base_uses.f90" 86 87 IMPLICIT NONE 88 89 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_fields_all' 90 91 PRIVATE 92 PUBLIC :: force_field_unique_bond, & 93 force_field_unique_bend, & 94 force_field_unique_ub, & 95 force_field_unique_tors, & 96 force_field_unique_impr, & 97 force_field_unique_opbend, & 98 force_field_pack_bond, & 99 force_field_pack_bend, & 100 force_field_pack_ub, & 101 force_field_pack_tors, & 102 force_field_pack_impr, & 103 force_field_pack_opbend, & 104 force_field_pack_charge, & 105 force_field_pack_charges, & 106 force_field_pack_radius, & 107 force_field_pack_pol, & 108 force_field_pack_shell, & 109 force_field_pack_nonbond14, & 110 force_field_pack_nonbond, & 111 force_field_pack_splines, & 112 force_field_pack_eicut, & 113 force_field_pack_damp 114 115CONTAINS 116 117! ************************************************************************************************** 118!> \brief Determine the number of unique bond kind and allocate bond_kind_set 119!> \param particle_set ... 120!> \param molecule_kind_set ... 121!> \param molecule_set ... 122!> \param ff_type ... 123! ************************************************************************************************** 124 SUBROUTINE force_field_unique_bond(particle_set, & 125 molecule_kind_set, molecule_set, ff_type) 126 127 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 128 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 129 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 130 TYPE(force_field_type), INTENT(INOUT) :: ff_type 131 132 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_bond', & 133 routineP = moduleN//':'//routineN 134 135 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, name_atm_b, & 136 name_atm_b2 137 INTEGER :: atm_a, atm_b, counter, first, handle2, & 138 i, j, k, last, natom, nbond 139 INTEGER, DIMENSION(:), POINTER :: molecule_list 140 INTEGER, POINTER :: map_bond_kind(:) 141 LOGICAL :: found 142 TYPE(atomic_kind_type), POINTER :: atomic_kind 143 TYPE(bond_kind_type), DIMENSION(:), POINTER :: bond_kind_set 144 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list 145 TYPE(molecule_kind_type), POINTER :: molecule_kind 146 TYPE(molecule_type), POINTER :: molecule 147 148 CALL timeset(routineN, handle2) 149 DO i = 1, SIZE(molecule_kind_set) 150 molecule_kind => molecule_kind_set(i) 151 CALL get_molecule_kind(molecule_kind=molecule_kind, & 152 molecule_list=molecule_list, & 153 natom=natom, & 154 nbond=nbond, bond_list=bond_list) 155 molecule => molecule_set(molecule_list(1)) 156 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 157 IF (nbond > 0) THEN 158 ALLOCATE (map_bond_kind(nbond)) 159 counter = 0 160 IF ((ff_type%ff_type == do_ff_g96) .OR. (ff_type%ff_type == do_ff_g87)) THEN 161 DO j = 1, nbond 162 map_bond_kind(j) = j 163 END DO 164 counter = nbond 165 ELSE 166 DO j = 1, nbond 167 atm_a = bond_list(j)%a 168 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 169 CALL get_atomic_kind(atomic_kind=atomic_kind, & 170 name=name_atm_a) 171 atm_b = bond_list(j)%b 172 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 173 CALL get_atomic_kind(atomic_kind=atomic_kind, & 174 name=name_atm_b) 175 found = .FALSE. 176 DO k = 1, j - 1 177 atm_a = bond_list(k)%a 178 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 179 CALL get_atomic_kind(atomic_kind=atomic_kind, & 180 name=name_atm_a2) 181 atm_b = bond_list(k)%b 182 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 183 CALL get_atomic_kind(atomic_kind=atomic_kind, & 184 name=name_atm_b2) 185 IF ((((name_atm_a) == (name_atm_a2)) .AND. & 186 ((name_atm_b) == (name_atm_b2))) .OR. & 187 (((name_atm_a) == (name_atm_b2)) .AND. & 188 ((name_atm_b) == (name_atm_a2)))) THEN 189 found = .TRUE. 190 map_bond_kind(j) = map_bond_kind(k) 191 EXIT 192 END IF 193 END DO 194 IF (.NOT. found) THEN 195 counter = counter + 1 196 map_bond_kind(j) = counter 197 END IF 198 END DO 199 END IF 200 NULLIFY (bond_kind_set) 201 CALL allocate_bond_kind_set(bond_kind_set, counter) 202 DO j = 1, nbond 203 bond_list(j)%bond_kind => bond_kind_set(map_bond_kind(j)) 204 END DO 205 CALL set_molecule_kind(molecule_kind=molecule_kind, & 206 bond_kind_set=bond_kind_set, bond_list=bond_list) 207 DEALLOCATE (map_bond_kind) 208 END IF 209 END DO 210 CALL timestop(handle2) 211 212 END SUBROUTINE force_field_unique_bond 213 214! ************************************************************************************************** 215!> \brief Determine the number of unique bend kind and allocate bend_kind_set 216!> \param particle_set ... 217!> \param molecule_kind_set ... 218!> \param molecule_set ... 219!> \param ff_type ... 220! ************************************************************************************************** 221 SUBROUTINE force_field_unique_bend(particle_set, & 222 molecule_kind_set, molecule_set, ff_type) 223 224 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 225 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 226 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 227 TYPE(force_field_type), INTENT(INOUT) :: ff_type 228 229 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_bend', & 230 routineP = moduleN//':'//routineN 231 232 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, name_atm_b, & 233 name_atm_b2, name_atm_c, name_atm_c2 234 INTEGER :: atm_a, atm_b, atm_c, counter, first, & 235 handle2, i, j, k, last, natom, nbend 236 INTEGER, DIMENSION(:), POINTER :: molecule_list 237 INTEGER, POINTER :: map_bend_kind(:) 238 LOGICAL :: found 239 TYPE(atomic_kind_type), POINTER :: atomic_kind 240 TYPE(bend_kind_type), DIMENSION(:), POINTER :: bend_kind_set 241 TYPE(bend_type), DIMENSION(:), POINTER :: bend_list 242 TYPE(molecule_kind_type), POINTER :: molecule_kind 243 TYPE(molecule_type), POINTER :: molecule 244 245 CALL timeset(routineN, handle2) 246 DO i = 1, SIZE(molecule_kind_set) 247 molecule_kind => molecule_kind_set(i) 248 CALL get_molecule_kind(molecule_kind=molecule_kind, & 249 molecule_list=molecule_list, & 250 natom=natom, & 251 nbend=nbend, bend_list=bend_list) 252 molecule => molecule_set(molecule_list(1)) 253 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 254 IF (nbend > 0) THEN 255 ALLOCATE (map_bend_kind(nbend)) 256 counter = 0 257 IF ((ff_type%ff_type == do_ff_g96) .OR. (ff_type%ff_type == do_ff_g87)) THEN 258 DO j = 1, nbend 259 map_bend_kind(j) = j 260 END DO 261 counter = nbend 262 ELSE 263 DO j = 1, nbend 264 atm_a = bend_list(j)%a 265 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 266 CALL get_atomic_kind(atomic_kind=atomic_kind, & 267 name=name_atm_a) 268 atm_b = bend_list(j)%b 269 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 270 CALL get_atomic_kind(atomic_kind=atomic_kind, & 271 name=name_atm_b) 272 atm_c = bend_list(j)%c 273 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 274 CALL get_atomic_kind(atomic_kind=atomic_kind, & 275 name=name_atm_c) 276 found = .FALSE. 277 DO k = 1, j - 1 278 atm_a = bend_list(k)%a 279 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 280 CALL get_atomic_kind(atomic_kind=atomic_kind, & 281 name=name_atm_a2) 282 atm_b = bend_list(k)%b 283 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 284 CALL get_atomic_kind(atomic_kind=atomic_kind, & 285 name=name_atm_b2) 286 atm_c = bend_list(k)%c 287 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 288 CALL get_atomic_kind(atomic_kind=atomic_kind, & 289 name=name_atm_c2) 290 IF ((((name_atm_a) == (name_atm_a2)) .AND. & 291 ((name_atm_b) == (name_atm_b2)) .AND. & 292 ((name_atm_c) == (name_atm_c2))) .OR. & 293 (((name_atm_a) == (name_atm_c2)) .AND. & 294 ((name_atm_b) == (name_atm_b2)) .AND. & 295 ((name_atm_c) == (name_atm_a2)))) THEN 296 found = .TRUE. 297 map_bend_kind(j) = map_bend_kind(k) 298 EXIT 299 END IF 300 END DO 301 IF (.NOT. found) THEN 302 counter = counter + 1 303 map_bend_kind(j) = counter 304 END IF 305 END DO 306 END IF 307 NULLIFY (bend_kind_set) 308 CALL allocate_bend_kind_set(bend_kind_set, counter) 309 DO j = 1, nbend 310 bend_list(j)%bend_kind => bend_kind_set(map_bend_kind(j)) 311 END DO 312 CALL set_molecule_kind(molecule_kind=molecule_kind, & 313 bend_kind_set=bend_kind_set, bend_list=bend_list) 314 DEALLOCATE (map_bend_kind) 315 END IF 316 END DO 317 CALL timestop(handle2) 318 319 END SUBROUTINE force_field_unique_bend 320 321! ************************************************************************************************** 322!> \brief Determine the number of unique Urey-Bradley kind and allocate ub_kind_set 323!> \param particle_set ... 324!> \param molecule_kind_set ... 325!> \param molecule_set ... 326! ************************************************************************************************** 327 SUBROUTINE force_field_unique_ub(particle_set, & 328 molecule_kind_set, molecule_set) 329 330 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 331 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 332 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 333 334 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_ub', & 335 routineP = moduleN//':'//routineN 336 337 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, name_atm_b, & 338 name_atm_b2, name_atm_c, name_atm_c2 339 INTEGER :: atm_a, atm_b, atm_c, counter, first, & 340 handle2, i, j, k, last, natom, nub 341 INTEGER, DIMENSION(:), POINTER :: molecule_list 342 INTEGER, POINTER :: map_ub_kind(:) 343 LOGICAL :: found 344 TYPE(atomic_kind_type), POINTER :: atomic_kind 345 TYPE(molecule_kind_type), POINTER :: molecule_kind 346 TYPE(molecule_type), POINTER :: molecule 347 TYPE(ub_kind_type), DIMENSION(:), POINTER :: ub_kind_set 348 TYPE(ub_type), DIMENSION(:), POINTER :: ub_list 349 350 CALL timeset(routineN, handle2) 351 DO i = 1, SIZE(molecule_kind_set) 352 molecule_kind => molecule_kind_set(i) 353 CALL get_molecule_kind(molecule_kind=molecule_kind, & 354 molecule_list=molecule_list, & 355 natom=natom, & 356 nub=nub, ub_list=ub_list) 357 molecule => molecule_set(molecule_list(1)) 358 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 359 IF (nub > 0) THEN 360 ALLOCATE (map_ub_kind(nub)) 361 counter = 0 362 DO j = 1, nub 363 atm_a = ub_list(j)%a 364 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 365 CALL get_atomic_kind(atomic_kind=atomic_kind, & 366 name=name_atm_a) 367 atm_b = ub_list(j)%b 368 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 369 CALL get_atomic_kind(atomic_kind=atomic_kind, & 370 name=name_atm_b) 371 atm_c = ub_list(j)%c 372 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 373 CALL get_atomic_kind(atomic_kind=atomic_kind, & 374 name=name_atm_c) 375 found = .FALSE. 376 DO k = 1, j - 1 377 atm_a = ub_list(k)%a 378 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 379 CALL get_atomic_kind(atomic_kind=atomic_kind, & 380 name=name_atm_a2) 381 atm_b = ub_list(k)%b 382 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 383 CALL get_atomic_kind(atomic_kind=atomic_kind, & 384 name=name_atm_b2) 385 atm_c = ub_list(k)%c 386 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 387 CALL get_atomic_kind(atomic_kind=atomic_kind, & 388 name=name_atm_c2) 389 IF ((((name_atm_a) == (name_atm_a2)) .AND. & 390 ((name_atm_b) == (name_atm_b2)) .AND. & 391 ((name_atm_c) == (name_atm_c2))) .OR. & 392 (((name_atm_a) == (name_atm_c2)) .AND. & 393 ((name_atm_b) == (name_atm_b2)) .AND. & 394 ((name_atm_c) == (name_atm_a2)))) THEN 395 found = .TRUE. 396 map_ub_kind(j) = map_ub_kind(k) 397 EXIT 398 END IF 399 END DO 400 IF (.NOT. found) THEN 401 counter = counter + 1 402 map_ub_kind(j) = counter 403 END IF 404 END DO 405 CALL allocate_ub_kind_set(ub_kind_set, counter) 406 DO j = 1, nub 407 ub_list(j)%ub_kind => ub_kind_set(map_ub_kind(j)) 408 END DO 409 CALL set_molecule_kind(molecule_kind=molecule_kind, & 410 ub_kind_set=ub_kind_set, ub_list=ub_list) 411 DEALLOCATE (map_ub_kind) 412 END IF 413 END DO 414 CALL timestop(handle2) 415 416 END SUBROUTINE force_field_unique_ub 417 418! ************************************************************************************************** 419!> \brief Determine the number of unique torsion kind and allocate torsion_kind_set 420!> \param particle_set ... 421!> \param molecule_kind_set ... 422!> \param molecule_set ... 423!> \param ff_type ... 424! ************************************************************************************************** 425 SUBROUTINE force_field_unique_tors(particle_set, & 426 molecule_kind_set, molecule_set, ff_type) 427 428 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 429 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 430 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 431 TYPE(force_field_type), INTENT(INOUT) :: ff_type 432 433 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_tors', & 434 routineP = moduleN//':'//routineN 435 436 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, name_atm_b, & 437 name_atm_b2, name_atm_c, name_atm_c2, & 438 name_atm_d, name_atm_d2 439 INTEGER :: atm_a, atm_b, atm_c, atm_d, counter, & 440 first, handle2, i, j, k, last, natom, & 441 ntorsion 442 INTEGER, DIMENSION(:), POINTER :: molecule_list 443 INTEGER, POINTER :: map_torsion_kind(:) 444 LOGICAL :: found 445 TYPE(atomic_kind_type), POINTER :: atomic_kind 446 TYPE(molecule_kind_type), POINTER :: molecule_kind 447 TYPE(molecule_type), POINTER :: molecule 448 TYPE(torsion_kind_type), DIMENSION(:), POINTER :: torsion_kind_set 449 TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list 450 451 CALL timeset(routineN, handle2) 452 DO i = 1, SIZE(molecule_kind_set) 453 molecule_kind => molecule_kind_set(i) 454 CALL get_molecule_kind(molecule_kind=molecule_kind, & 455 molecule_list=molecule_list, & 456 natom=natom, & 457 ntorsion=ntorsion, torsion_list=torsion_list) 458 molecule => molecule_set(molecule_list(1)) 459 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 460 IF (ntorsion > 0) THEN 461 ALLOCATE (map_torsion_kind(ntorsion)) 462 counter = 0 463 IF ((ff_type%ff_type == do_ff_g96) .OR. (ff_type%ff_type == do_ff_g87)) THEN 464 DO j = 1, ntorsion 465 map_torsion_kind(j) = j 466 END DO 467 counter = ntorsion 468 ELSE 469 DO j = 1, ntorsion 470 atm_a = torsion_list(j)%a 471 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 472 CALL get_atomic_kind(atomic_kind=atomic_kind, & 473 name=name_atm_a) 474 atm_b = torsion_list(j)%b 475 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 476 CALL get_atomic_kind(atomic_kind=atomic_kind, & 477 name=name_atm_b) 478 atm_c = torsion_list(j)%c 479 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 480 CALL get_atomic_kind(atomic_kind=atomic_kind, & 481 name=name_atm_c) 482 atm_d = torsion_list(j)%d 483 atomic_kind => particle_set(atm_d + first - 1)%atomic_kind 484 CALL get_atomic_kind(atomic_kind=atomic_kind, & 485 name=name_atm_d) 486 found = .FALSE. 487 DO k = 1, j - 1 488 atm_a = torsion_list(k)%a 489 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 490 CALL get_atomic_kind(atomic_kind=atomic_kind, & 491 name=name_atm_a2) 492 atm_b = torsion_list(k)%b 493 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 494 CALL get_atomic_kind(atomic_kind=atomic_kind, & 495 name=name_atm_b2) 496 atm_c = torsion_list(k)%c 497 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 498 CALL get_atomic_kind(atomic_kind=atomic_kind, & 499 name=name_atm_c2) 500 atm_d = torsion_list(k)%d 501 atomic_kind => particle_set(atm_d + first - 1)%atomic_kind 502 CALL get_atomic_kind(atomic_kind=atomic_kind, & 503 name=name_atm_d2) 504 IF ((((name_atm_a) == (name_atm_a2)) .AND. & 505 ((name_atm_b) == (name_atm_b2)) .AND. & 506 ((name_atm_c) == (name_atm_c2)) .AND. & 507 ((name_atm_d) == (name_atm_d2))) .OR. & 508 (((name_atm_a) == (name_atm_d2)) .AND. & 509 ((name_atm_b) == (name_atm_c2)) .AND. & 510 ((name_atm_c) == (name_atm_b2)) .AND. & 511 ((name_atm_d) == (name_atm_a2)))) THEN 512 found = .TRUE. 513 map_torsion_kind(j) = map_torsion_kind(k) 514 EXIT 515 END IF 516 END DO 517 IF (.NOT. found) THEN 518 counter = counter + 1 519 map_torsion_kind(j) = counter 520 END IF 521 END DO 522 END IF 523 NULLIFY (torsion_kind_set) 524 CALL allocate_torsion_kind_set(torsion_kind_set, counter) 525 DO j = 1, ntorsion 526 torsion_list(j)%torsion_kind => torsion_kind_set(map_torsion_kind(j)) 527 END DO 528 CALL set_molecule_kind(molecule_kind=molecule_kind, & 529 torsion_kind_set=torsion_kind_set, torsion_list=torsion_list) 530 DEALLOCATE (map_torsion_kind) 531 END IF 532 END DO 533 CALL timestop(handle2) 534 535 END SUBROUTINE force_field_unique_tors 536 537! ************************************************************************************************** 538!> \brief Determine the number of unique impr kind and allocate impr_kind_set 539!> \param particle_set ... 540!> \param molecule_kind_set ... 541!> \param molecule_set ... 542!> \param ff_type ... 543! ************************************************************************************************** 544 SUBROUTINE force_field_unique_impr(particle_set, & 545 molecule_kind_set, molecule_set, ff_type) 546 547 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 548 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 549 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 550 TYPE(force_field_type), INTENT(INOUT) :: ff_type 551 552 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_impr', & 553 routineP = moduleN//':'//routineN 554 555 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, name_atm_b, & 556 name_atm_b2, name_atm_c, name_atm_c2, & 557 name_atm_d, name_atm_d2 558 INTEGER :: atm_a, atm_b, atm_c, atm_d, counter, & 559 first, handle2, i, j, k, last, natom, & 560 nimpr 561 INTEGER, DIMENSION(:), POINTER :: molecule_list 562 INTEGER, POINTER :: map_impr_kind(:) 563 LOGICAL :: found 564 TYPE(atomic_kind_type), POINTER :: atomic_kind 565 TYPE(impr_kind_type), DIMENSION(:), POINTER :: impr_kind_set 566 TYPE(impr_type), DIMENSION(:), POINTER :: impr_list 567 TYPE(molecule_kind_type), POINTER :: molecule_kind 568 TYPE(molecule_type), POINTER :: molecule 569 570 CALL timeset(routineN, handle2) 571 DO i = 1, SIZE(molecule_kind_set) 572 molecule_kind => molecule_kind_set(i) 573 CALL get_molecule_kind(molecule_kind=molecule_kind, & 574 molecule_list=molecule_list, & 575 natom=natom, & 576 nimpr=nimpr, impr_list=impr_list) 577 molecule => molecule_set(molecule_list(1)) 578 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 579 IF (nimpr > 0) THEN 580 ALLOCATE (map_impr_kind(nimpr)) 581 counter = 0 582 IF ((ff_type%ff_type == do_ff_g96) .OR. (ff_type%ff_type == do_ff_g87)) THEN 583 DO j = 1, nimpr 584 map_impr_kind(j) = j 585 END DO 586 counter = nimpr 587 ELSE 588 DO j = 1, nimpr 589 atm_a = impr_list(j)%a 590 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 591 CALL get_atomic_kind(atomic_kind=atomic_kind, & 592 name=name_atm_a) 593 atm_b = impr_list(j)%b 594 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 595 CALL get_atomic_kind(atomic_kind=atomic_kind, & 596 name=name_atm_b) 597 atm_c = impr_list(j)%c 598 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 599 CALL get_atomic_kind(atomic_kind=atomic_kind, & 600 name=name_atm_c) 601 atm_d = impr_list(j)%d 602 atomic_kind => particle_set(atm_d + first - 1)%atomic_kind 603 CALL get_atomic_kind(atomic_kind=atomic_kind, & 604 name=name_atm_d) 605 found = .FALSE. 606 DO k = 1, j - 1 607 atm_a = impr_list(k)%a 608 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 609 CALL get_atomic_kind(atomic_kind=atomic_kind, & 610 name=name_atm_a2) 611 atm_b = impr_list(k)%b 612 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 613 CALL get_atomic_kind(atomic_kind=atomic_kind, & 614 name=name_atm_b2) 615 atm_c = impr_list(k)%c 616 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 617 CALL get_atomic_kind(atomic_kind=atomic_kind, & 618 name=name_atm_c2) 619 atm_d = impr_list(k)%d 620 atomic_kind => particle_set(atm_d + first - 1)%atomic_kind 621 CALL get_atomic_kind(atomic_kind=atomic_kind, & 622 name=name_atm_d2) 623 IF ((((name_atm_a) == (name_atm_a2)) .AND. & 624 ((name_atm_b) == (name_atm_b2)) .AND. & 625 ((name_atm_c) == (name_atm_c2)) .AND. & 626 ((name_atm_d) == (name_atm_d2))) .OR. & 627 (((name_atm_a) == (name_atm_a2)) .AND. & 628 ((name_atm_b) == (name_atm_b2)) .AND. & 629 ((name_atm_c) == (name_atm_d2)) .AND. & 630 ((name_atm_d) == (name_atm_c2)))) THEN 631 found = .TRUE. 632 map_impr_kind(j) = map_impr_kind(k) 633 EXIT 634 END IF 635 END DO 636 IF (.NOT. found) THEN 637 counter = counter + 1 638 map_impr_kind(j) = counter 639 END IF 640 END DO 641 END IF 642 NULLIFY (impr_kind_set) 643 CALL allocate_impr_kind_set(impr_kind_set, counter) 644 DO j = 1, nimpr 645 impr_list(j)%impr_kind => impr_kind_set(map_impr_kind(j)) 646 END DO 647 CALL set_molecule_kind(molecule_kind=molecule_kind, & 648 impr_kind_set=impr_kind_set, impr_list=impr_list) 649 DEALLOCATE (map_impr_kind) 650 END IF 651 END DO 652 CALL timestop(handle2) 653 654 END SUBROUTINE force_field_unique_impr 655 656! ************************************************************************************************** 657!> \brief Determine the number of unique opbend kind and allocate opbend_kind_set 658!> based on the present impropers. With each improper, there also 659!> corresponds a opbend 660!> \param particle_set ... 661!> \param molecule_kind_set ... 662!> \param molecule_set ... 663!> \param ff_type ... 664! ************************************************************************************************** 665 SUBROUTINE force_field_unique_opbend(particle_set, & 666 molecule_kind_set, molecule_set, ff_type) 667 668 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 669 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 670 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 671 TYPE(force_field_type), INTENT(INOUT) :: ff_type 672 673 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_unique_opbend', & 674 routineP = moduleN//':'//routineN 675 676 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a2, name_atm_b, & 677 name_atm_b2, name_atm_c, name_atm_c2, & 678 name_atm_d, name_atm_d2 679 INTEGER :: atm_a, atm_b, atm_c, atm_d, counter, & 680 first, handle2, i, j, k, last, natom, & 681 nopbend 682 INTEGER, DIMENSION(:), POINTER :: molecule_list 683 INTEGER, POINTER :: map_opbend_kind(:) 684 LOGICAL :: found 685 TYPE(atomic_kind_type), POINTER :: atomic_kind 686 TYPE(molecule_kind_type), POINTER :: molecule_kind 687 TYPE(molecule_type), POINTER :: molecule 688 TYPE(opbend_kind_type), DIMENSION(:), POINTER :: opbend_kind_set 689 TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list 690 691 CALL timeset(routineN, handle2) 692 DO i = 1, SIZE(molecule_kind_set) 693 molecule_kind => molecule_kind_set(i) 694 CALL get_molecule_kind(molecule_kind=molecule_kind, & 695 molecule_list=molecule_list, & 696 natom=natom, & 697 nopbend=nopbend, opbend_list=opbend_list) 698 molecule => molecule_set(molecule_list(1)) 699 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 700 IF (nopbend > 0) THEN 701 ALLOCATE (map_opbend_kind(nopbend)) 702 counter = 0 703 IF ((ff_type%ff_type == do_ff_g96) .OR. (ff_type%ff_type == do_ff_g87)) THEN 704 DO j = 1, nopbend 705 map_opbend_kind(j) = j 706 END DO 707 counter = nopbend 708 ELSE 709 DO j = 1, nopbend 710 atm_a = opbend_list(j)%a 711 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 712 CALL get_atomic_kind(atomic_kind=atomic_kind, & 713 name=name_atm_a) 714 atm_b = opbend_list(j)%b 715 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 716 CALL get_atomic_kind(atomic_kind=atomic_kind, & 717 name=name_atm_b) 718 atm_c = opbend_list(j)%c 719 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 720 CALL get_atomic_kind(atomic_kind=atomic_kind, & 721 name=name_atm_c) 722 atm_d = opbend_list(j)%d 723 atomic_kind => particle_set(atm_d + first - 1)%atomic_kind 724 CALL get_atomic_kind(atomic_kind=atomic_kind, & 725 name=name_atm_d) 726 found = .FALSE. 727 DO k = 1, j - 1 728 atm_a = opbend_list(k)%a 729 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 730 CALL get_atomic_kind(atomic_kind=atomic_kind, & 731 name=name_atm_a2) 732 atm_b = opbend_list(k)%b 733 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 734 CALL get_atomic_kind(atomic_kind=atomic_kind, & 735 name=name_atm_b2) 736 atm_c = opbend_list(k)%c 737 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 738 CALL get_atomic_kind(atomic_kind=atomic_kind, & 739 name=name_atm_c2) 740 atm_d = opbend_list(k)%d 741 atomic_kind => particle_set(atm_d + first - 1)%atomic_kind 742 CALL get_atomic_kind(atomic_kind=atomic_kind, & 743 name=name_atm_d2) 744 IF ((((name_atm_a) == (name_atm_a2)) .AND. & 745 ((name_atm_b) == (name_atm_b2)) .AND. & 746 ((name_atm_c) == (name_atm_c2)) .AND. & 747 ((name_atm_d) == (name_atm_d2))) .OR. & 748 (((name_atm_a) == (name_atm_a2)) .AND. & 749 ((name_atm_b) == (name_atm_c2)) .AND. & 750 ((name_atm_c) == (name_atm_b2)) .AND. & 751 ((name_atm_d) == (name_atm_d2)))) THEN 752 found = .TRUE. 753 map_opbend_kind(j) = map_opbend_kind(k) 754 EXIT 755 END IF 756 END DO 757 IF (.NOT. found) THEN 758 counter = counter + 1 759 map_opbend_kind(j) = counter 760 END IF 761 END DO 762 END IF 763 NULLIFY (opbend_kind_set) 764 CALL allocate_opbend_kind_set(opbend_kind_set, counter) 765 DO j = 1, nopbend 766 opbend_list(j)%opbend_kind => opbend_kind_set(map_opbend_kind(j)) 767 END DO 768 CALL set_molecule_kind(molecule_kind=molecule_kind, & 769 opbend_kind_set=opbend_kind_set, opbend_list=opbend_list) 770 DEALLOCATE (map_opbend_kind) 771 END IF 772 END DO 773 CALL timestop(handle2) 774 775 END SUBROUTINE force_field_unique_opbend 776 777! ************************************************************************************************** 778!> \brief Pack in bonds information needed for the force_field 779!> \param particle_set ... 780!> \param molecule_kind_set ... 781!> \param molecule_set ... 782!> \param fatal ... 783!> \param Ainfo ... 784!> \param chm_info ... 785!> \param inp_info ... 786!> \param gro_info ... 787!> \param amb_info ... 788! ************************************************************************************************** 789 SUBROUTINE force_field_pack_bond(particle_set, molecule_kind_set, molecule_set, & 790 fatal, Ainfo, chm_info, inp_info, gro_info, amb_info) 791 792 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 793 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 794 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 795 LOGICAL :: fatal 796 CHARACTER(LEN=default_string_length), & 797 DIMENSION(:), POINTER :: Ainfo 798 TYPE(charmm_info_type), POINTER :: chm_info 799 TYPE(input_info_type), POINTER :: inp_info 800 TYPE(gromos_info_type), POINTER :: gro_info 801 TYPE(amber_info_type), POINTER :: amb_info 802 803 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_bond', & 804 routineP = moduleN//':'//routineN 805 806 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b 807 INTEGER :: atm_a, atm_b, first, handle2, i, itype, & 808 j, k, last, natom, nbond 809 INTEGER, DIMENSION(:), POINTER :: molecule_list 810 LOGICAL :: found, only_qm 811 TYPE(atomic_kind_type), POINTER :: atomic_kind 812 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list 813 TYPE(molecule_kind_type), POINTER :: molecule_kind 814 TYPE(molecule_type), POINTER :: molecule 815 816 CALL timeset(routineN, handle2) 817 818 DO i = 1, SIZE(molecule_kind_set) 819 molecule_kind => molecule_kind_set(i) 820 CALL get_molecule_kind(molecule_kind=molecule_kind, & 821 molecule_list=molecule_list, & 822 natom=natom, & 823 nbond=nbond, bond_list=bond_list) 824 molecule => molecule_set(molecule_list(1)) 825 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 826 DO j = 1, nbond 827 atm_a = bond_list(j)%a 828 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 829 CALL get_atomic_kind(atomic_kind=atomic_kind, & 830 name=name_atm_a) 831 atm_b = bond_list(j)%b 832 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 833 CALL get_atomic_kind(atomic_kind=atomic_kind, & 834 name=name_atm_b) 835 found = .FALSE. 836 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b) 837 CALL uppercase(name_atm_a) 838 CALL uppercase(name_atm_b) 839 840 ! loop over params from GROMOS 841 IF (ASSOCIATED(gro_info%bond_k)) THEN 842 k = SIZE(gro_info%bond_k) 843 itype = bond_list(j)%itype 844 IF (itype <= k) THEN 845 bond_list(j)%bond_kind%k(1) = gro_info%bond_k(itype) 846 bond_list(j)%bond_kind%r0 = gro_info%bond_r0(itype) 847 ELSE 848 itype = itype - k 849 bond_list(j)%bond_kind%k(1) = gro_info%solvent_k(itype) 850 bond_list(j)%bond_kind%r0 = gro_info%solvent_r0(itype) 851 END IF 852 bond_list(j)%bond_kind%id_type = gro_info%ff_gromos_type 853 bond_list(j)%id_type = gro_info%ff_gromos_type 854 found = .TRUE. 855 END IF 856 857 ! loop over params from CHARMM 858 IF (ASSOCIATED(chm_info%bond_a)) THEN 859 DO k = 1, SIZE(chm_info%bond_a) 860 IF ((((chm_info%bond_a(k)) == (name_atm_a)) .AND. & 861 ((chm_info%bond_b(k)) == (name_atm_b))) .OR. & 862 (((chm_info%bond_a(k)) == (name_atm_b)) .AND. & 863 ((chm_info%bond_b(k)) == (name_atm_a)))) THEN 864 bond_list(j)%bond_kind%id_type = do_ff_charmm 865 bond_list(j)%bond_kind%k(1) = chm_info%bond_k(k) 866 bond_list(j)%bond_kind%r0 = chm_info%bond_r0(k) 867 CALL issue_duplications(found, "Bond", name_atm_a, name_atm_b) 868 found = .TRUE. 869 EXIT 870 END IF 871 END DO 872 END IF 873 874 ! loop over params from AMBER 875 IF (ASSOCIATED(amb_info%bond_a)) THEN 876 DO k = 1, SIZE(amb_info%bond_a) 877 IF ((((amb_info%bond_a(k)) == (name_atm_a)) .AND. & 878 ((amb_info%bond_b(k)) == (name_atm_b))) .OR. & 879 (((amb_info%bond_a(k)) == (name_atm_b)) .AND. & 880 ((amb_info%bond_b(k)) == (name_atm_a)))) THEN 881 bond_list(j)%bond_kind%id_type = do_ff_amber 882 bond_list(j)%bond_kind%k(1) = amb_info%bond_k(k) 883 bond_list(j)%bond_kind%r0 = amb_info%bond_r0(k) 884 CALL issue_duplications(found, "Bond", name_atm_a, name_atm_b) 885 found = .TRUE. 886 EXIT 887 END IF 888 END DO 889 END IF 890 891 ! always have the input param last to overwrite all the other ones 892 IF (ASSOCIATED(inp_info%bond_a)) THEN 893 DO k = 1, SIZE(inp_info%bond_a) 894 IF ((((inp_info%bond_a(k)) == (name_atm_a)) .AND. & 895 ((inp_info%bond_b(k)) == (name_atm_b))) .OR. & 896 (((inp_info%bond_a(k)) == (name_atm_b)) .AND. & 897 ((inp_info%bond_b(k)) == (name_atm_a)))) THEN 898 bond_list(j)%bond_kind%id_type = inp_info%bond_kind(k) 899 bond_list(j)%bond_kind%k(:) = inp_info%bond_k(:, k) 900 bond_list(j)%bond_kind%r0 = inp_info%bond_r0(k) 901 bond_list(j)%bond_kind%cs = inp_info%bond_cs(k) 902 CALL issue_duplications(found, "Bond", name_atm_a, name_atm_b) 903 found = .TRUE. 904 EXIT 905 END IF 906 END DO 907 END IF 908 909 IF (.NOT. found) CALL store_FF_missing_par(atm1=TRIM(name_atm_a), & 910 atm2=TRIM(name_atm_b), & 911 fatal=fatal, & 912 type_name="Bond", & 913 array=Ainfo) 914 ! QM/MM modifications 915 IF (only_qm) THEN 916 bond_list(j)%id_type = do_ff_undef 917 bond_list(j)%bond_kind%id_type = do_ff_undef 918 END IF 919 END DO 920 CALL set_molecule_kind(molecule_kind=molecule_kind, & 921 bond_list=bond_list) 922 END DO 923 CALL timestop(handle2) 924 925 END SUBROUTINE force_field_pack_bond 926 927! ************************************************************************************************** 928!> \brief Pack in bends information needed for the force_field 929!> \param particle_set ... 930!> \param molecule_kind_set ... 931!> \param molecule_set ... 932!> \param fatal ... 933!> \param Ainfo ... 934!> \param chm_info ... 935!> \param inp_info ... 936!> \param gro_info ... 937!> \param amb_info ... 938! ************************************************************************************************** 939 SUBROUTINE force_field_pack_bend(particle_set, molecule_kind_set, molecule_set, & 940 fatal, Ainfo, chm_info, inp_info, gro_info, amb_info) 941 942 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 943 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 944 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 945 LOGICAL :: fatal 946 CHARACTER(LEN=default_string_length), & 947 DIMENSION(:), POINTER :: Ainfo 948 TYPE(charmm_info_type), POINTER :: chm_info 949 TYPE(input_info_type), POINTER :: inp_info 950 TYPE(gromos_info_type), POINTER :: gro_info 951 TYPE(amber_info_type), POINTER :: amb_info 952 953 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_bend', & 954 routineP = moduleN//':'//routineN 955 956 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, name_atm_c 957 INTEGER :: atm_a, atm_b, atm_c, first, handle2, i, & 958 itype, j, k, l, last, natom, nbend 959 INTEGER, DIMENSION(:), POINTER :: molecule_list 960 LOGICAL :: found, only_qm 961 TYPE(atomic_kind_type), POINTER :: atomic_kind 962 TYPE(bend_type), DIMENSION(:), POINTER :: bend_list 963 TYPE(molecule_kind_type), POINTER :: molecule_kind 964 TYPE(molecule_type), POINTER :: molecule 965 966 CALL timeset(routineN, handle2) 967 968 DO i = 1, SIZE(molecule_kind_set) 969 molecule_kind => molecule_kind_set(i) 970 CALL get_molecule_kind(molecule_kind=molecule_kind, & 971 molecule_list=molecule_list, & 972 natom=natom, & 973 nbend=nbend, bend_list=bend_list) 974 molecule => molecule_set(molecule_list(1)) 975 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 976 DO j = 1, nbend 977 atm_a = bend_list(j)%a 978 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 979 CALL get_atomic_kind(atomic_kind=atomic_kind, & 980 name=name_atm_a) 981 atm_b = bend_list(j)%b 982 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 983 CALL get_atomic_kind(atomic_kind=atomic_kind, & 984 name=name_atm_b) 985 atm_c = bend_list(j)%c 986 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 987 CALL get_atomic_kind(atomic_kind=atomic_kind, & 988 name=name_atm_c) 989 found = .FALSE. 990 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c) 991 CALL uppercase(name_atm_a) 992 CALL uppercase(name_atm_b) 993 CALL uppercase(name_atm_c) 994 995 ! loop over params from GROMOS 996 IF (ASSOCIATED(gro_info%bend_k)) THEN 997 k = SIZE(gro_info%bend_k) 998 itype = bend_list(j)%itype 999 IF (itype > 0) THEN 1000 bend_list(j)%bend_kind%k = gro_info%bend_k(itype) 1001 bend_list(j)%bend_kind%theta0 = gro_info%bend_theta0(itype) 1002 ELSE 1003 bend_list(j)%bend_kind%k = gro_info%bend_k(itype/k) 1004 bend_list(j)%bend_kind%theta0 = gro_info%bend_theta0(itype/k) 1005 END IF 1006 bend_list(j)%bend_kind%id_type = gro_info%ff_gromos_type 1007 bend_list(j)%id_type = gro_info%ff_gromos_type 1008 found = .TRUE. 1009 END IF 1010 1011 ! loop over params from CHARMM 1012 IF (ASSOCIATED(chm_info%bend_a)) THEN 1013 DO k = 1, SIZE(chm_info%bend_a) 1014 IF ((((chm_info%bend_a(k)) == (name_atm_a)) .AND. & 1015 ((chm_info%bend_b(k)) == (name_atm_b)) .AND. & 1016 ((chm_info%bend_c(k)) == (name_atm_c))) .OR. & 1017 (((chm_info%bend_a(k)) == (name_atm_c)) .AND. & 1018 ((chm_info%bend_b(k)) == (name_atm_b)) .AND. & 1019 ((chm_info%bend_c(k)) == (name_atm_a)))) THEN 1020 bend_list(j)%bend_kind%id_type = do_ff_charmm 1021 bend_list(j)%bend_kind%k = chm_info%bend_k(k) 1022 bend_list(j)%bend_kind%theta0 = chm_info%bend_theta0(k) 1023 CALL issue_duplications(found, "Bend", name_atm_a, name_atm_b, & 1024 name_atm_c) 1025 found = .TRUE. 1026 EXIT 1027 END IF 1028 END DO 1029 END IF 1030 1031 ! loop over params from AMBER 1032 IF (ASSOCIATED(amb_info%bend_a)) THEN 1033 DO k = 1, SIZE(amb_info%bend_a) 1034 IF ((((amb_info%bend_a(k)) == (name_atm_a)) .AND. & 1035 ((amb_info%bend_b(k)) == (name_atm_b)) .AND. & 1036 ((amb_info%bend_c(k)) == (name_atm_c))) .OR. & 1037 (((amb_info%bend_a(k)) == (name_atm_c)) .AND. & 1038 ((amb_info%bend_b(k)) == (name_atm_b)) .AND. & 1039 ((amb_info%bend_c(k)) == (name_atm_a)))) THEN 1040 bend_list(j)%bend_kind%id_type = do_ff_amber 1041 bend_list(j)%bend_kind%k = amb_info%bend_k(k) 1042 bend_list(j)%bend_kind%theta0 = amb_info%bend_theta0(k) 1043 CALL issue_duplications(found, "Bend", name_atm_a, name_atm_b, & 1044 name_atm_c) 1045 found = .TRUE. 1046 EXIT 1047 END IF 1048 END DO 1049 END IF 1050 1051 ! always have the input param last to overwrite all the other ones 1052 IF (ASSOCIATED(inp_info%bend_a)) THEN 1053 DO k = 1, SIZE(inp_info%bend_a) 1054 IF ((((inp_info%bend_a(k)) == (name_atm_a)) .AND. & 1055 ((inp_info%bend_b(k)) == (name_atm_b)) .AND. & 1056 ((inp_info%bend_c(k)) == (name_atm_c))) .OR. & 1057 (((inp_info%bend_a(k)) == (name_atm_c)) .AND. & 1058 ((inp_info%bend_b(k)) == (name_atm_b)) .AND. & 1059 ((inp_info%bend_c(k)) == (name_atm_a)))) THEN 1060 bend_list(j)%bend_kind%id_type = inp_info%bend_kind(k) 1061 bend_list(j)%bend_kind%k = inp_info%bend_k(k) 1062 bend_list(j)%bend_kind%theta0 = inp_info%bend_theta0(k) 1063 bend_list(j)%bend_kind%cb = inp_info%bend_cb(k) 1064 bend_list(j)%bend_kind%r012 = inp_info%bend_r012(k) 1065 bend_list(j)%bend_kind%r032 = inp_info%bend_r032(k) 1066 bend_list(j)%bend_kind%kbs12 = inp_info%bend_kbs12(k) 1067 bend_list(j)%bend_kind%kbs32 = inp_info%bend_kbs32(k) 1068 bend_list(j)%bend_kind%kss = inp_info%bend_kss(k) 1069 bend_list(j)%bend_kind%legendre%order = inp_info%bend_legendre(k)%order 1070 IF (bend_list(j)%bend_kind%legendre%order /= 0) THEN 1071 IF (ASSOCIATED(bend_list(j)%bend_kind%legendre%coeffs)) THEN 1072 DEALLOCATE (bend_list(j)%bend_kind%legendre%coeffs) 1073 END IF 1074 ALLOCATE (bend_list(j)%bend_kind%legendre%coeffs(bend_list(j)%bend_kind%legendre%order)) 1075 DO l = 1, bend_list(j)%bend_kind%legendre%order 1076 bend_list(j)%bend_kind%legendre%coeffs(l) = inp_info%bend_legendre(k)%coeffs(l) 1077 END DO 1078 END IF 1079 CALL issue_duplications(found, "Bend", name_atm_a, name_atm_b, & 1080 name_atm_c) 1081 found = .TRUE. 1082 EXIT 1083 END IF 1084 END DO 1085 END IF 1086 1087 IF (.NOT. found) CALL store_FF_missing_par(atm1=TRIM(name_atm_a), & 1088 atm2=TRIM(name_atm_b), & 1089 atm3=TRIM(name_atm_c), & 1090 fatal=fatal, & 1091 type_name="Angle", & 1092 array=Ainfo) 1093 ! QM/MM modifications 1094 IF (only_qm) THEN 1095 bend_list(j)%id_type = do_ff_undef 1096 bend_list(j)%bend_kind%id_type = do_ff_undef 1097 END IF 1098 END DO 1099 CALL set_molecule_kind(molecule_kind=molecule_kind, & 1100 bend_list=bend_list) 1101 END DO 1102 CALL timestop(handle2) 1103 1104 END SUBROUTINE force_field_pack_bend 1105 1106! ************************************************************************************************** 1107!> \brief Pack in Urey-Bradley information needed for the force_field 1108!> \param particle_set ... 1109!> \param molecule_kind_set ... 1110!> \param molecule_set ... 1111!> \param Ainfo ... 1112!> \param chm_info ... 1113!> \param inp_info ... 1114!> \param iw ... 1115! ************************************************************************************************** 1116 SUBROUTINE force_field_pack_ub(particle_set, molecule_kind_set, molecule_set, & 1117 Ainfo, chm_info, inp_info, iw) 1118 1119 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1120 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 1121 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 1122 CHARACTER(LEN=default_string_length), & 1123 DIMENSION(:), POINTER :: Ainfo 1124 TYPE(charmm_info_type), POINTER :: chm_info 1125 TYPE(input_info_type), POINTER :: inp_info 1126 INTEGER :: iw 1127 1128 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_ub', & 1129 routineP = moduleN//':'//routineN 1130 1131 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, name_atm_c 1132 INTEGER :: atm_a, atm_b, atm_c, first, handle2, i, & 1133 j, k, last, natom, nub 1134 INTEGER, DIMENSION(:), POINTER :: molecule_list 1135 LOGICAL :: found, only_qm 1136 TYPE(atomic_kind_type), POINTER :: atomic_kind 1137 TYPE(molecule_kind_type), POINTER :: molecule_kind 1138 TYPE(molecule_type), POINTER :: molecule 1139 TYPE(ub_type), DIMENSION(:), POINTER :: ub_list 1140 1141 CALL timeset(routineN, handle2) 1142 DO i = 1, SIZE(molecule_kind_set) 1143 molecule_kind => molecule_kind_set(i) 1144 CALL get_molecule_kind(molecule_kind=molecule_kind, & 1145 molecule_list=molecule_list, & 1146 natom=natom, & 1147 nub=nub, ub_list=ub_list) 1148 molecule => molecule_set(molecule_list(1)) 1149 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 1150 DO j = 1, nub 1151 atm_a = ub_list(j)%a 1152 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 1153 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1154 name=name_atm_a) 1155 atm_b = ub_list(j)%b 1156 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 1157 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1158 name=name_atm_b) 1159 atm_c = ub_list(j)%c 1160 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 1161 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1162 name=name_atm_c) 1163 found = .FALSE. 1164 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c) 1165 CALL uppercase(name_atm_a) 1166 CALL uppercase(name_atm_b) 1167 CALL uppercase(name_atm_c) 1168 1169 ! loop over params from GROMOS 1170 ! ikuo - None that I know... 1171 1172 ! loop over params from CHARMM 1173 IF (ASSOCIATED(chm_info%ub_a)) THEN 1174 DO k = 1, SIZE(chm_info%ub_a) 1175 IF ((((chm_info%ub_a(k)) == (name_atm_a)) .AND. & 1176 ((chm_info%ub_b(k)) == (name_atm_b)) .AND. & 1177 ((chm_info%ub_c(k)) == (name_atm_c))) .OR. & 1178 (((chm_info%ub_a(k)) == (name_atm_c)) .AND. & 1179 ((chm_info%ub_b(k)) == (name_atm_b)) .AND. & 1180 ((chm_info%ub_c(k)) == (name_atm_a)))) THEN 1181 ub_list(j)%ub_kind%id_type = do_ff_charmm 1182 ub_list(j)%ub_kind%k(1) = chm_info%ub_k(k) 1183 ub_list(j)%ub_kind%r0 = chm_info%ub_r0(k) 1184 IF (iw > 0) WRITE (iw, *) " Found UB ", TRIM(name_atm_a), " ", & 1185 TRIM(name_atm_b), " ", TRIM(name_atm_c) 1186 CALL issue_duplications(found, "Urey-Bradley", name_atm_a, & 1187 name_atm_b, name_atm_c) 1188 found = .TRUE. 1189 EXIT 1190 END IF 1191 END DO 1192 END IF 1193 1194 ! loop over params from AMBER 1195 ! teo - None that I know... 1196 1197 ! always have the input param last to overwrite all the other ones 1198 IF (ASSOCIATED(inp_info%ub_a)) THEN 1199 DO k = 1, SIZE(inp_info%ub_a) 1200 IF ((((inp_info%ub_a(k)) == (name_atm_a)) .AND. & 1201 ((inp_info%ub_b(k)) == (name_atm_b)) .AND. & 1202 ((inp_info%ub_c(k)) == (name_atm_c))) .OR. & 1203 (((inp_info%ub_a(k)) == (name_atm_c)) .AND. & 1204 ((inp_info%ub_b(k)) == (name_atm_b)) .AND. & 1205 ((inp_info%ub_c(k)) == (name_atm_a)))) THEN 1206 ub_list(j)%ub_kind%id_type = inp_info%ub_kind(k) 1207 ub_list(j)%ub_kind%k(:) = inp_info%ub_k(:, k) 1208 ub_list(j)%ub_kind%r0 = inp_info%ub_r0(k) 1209 CALL issue_duplications(found, "Urey-Bradley", name_atm_a, & 1210 name_atm_b, name_atm_c) 1211 found = .TRUE. 1212 EXIT 1213 END IF 1214 END DO 1215 END IF 1216 1217 IF (.NOT. found) THEN 1218 CALL store_FF_missing_par(atm1=TRIM(name_atm_a), & 1219 atm2=TRIM(name_atm_b), & 1220 atm3=TRIM(name_atm_c), & 1221 type_name="Urey-Bradley", & 1222 array=Ainfo) 1223 ub_list(j)%id_type = do_ff_undef 1224 ub_list(j)%ub_kind%id_type = do_ff_undef 1225 ub_list(j)%ub_kind%k = 0.0_dp 1226 ub_list(j)%ub_kind%r0 = 0.0_dp 1227 END IF 1228 1229 ! QM/MM modifications 1230 IF (only_qm) THEN 1231 ub_list(j)%id_type = do_ff_undef 1232 ub_list(j)%ub_kind%id_type = do_ff_undef 1233 END IF 1234 END DO 1235 CALL set_molecule_kind(molecule_kind=molecule_kind, & 1236 ub_list=ub_list) 1237 END DO 1238 CALL timestop(handle2) 1239 1240 END SUBROUTINE force_field_pack_ub 1241 1242! ************************************************************************************************** 1243!> \brief Pack in torsion information needed for the force_field 1244!> \param particle_set ... 1245!> \param molecule_kind_set ... 1246!> \param molecule_set ... 1247!> \param Ainfo ... 1248!> \param chm_info ... 1249!> \param inp_info ... 1250!> \param gro_info ... 1251!> \param amb_info ... 1252!> \param iw ... 1253! ************************************************************************************************** 1254 SUBROUTINE force_field_pack_tors(particle_set, molecule_kind_set, molecule_set, & 1255 Ainfo, chm_info, inp_info, gro_info, amb_info, iw) 1256 1257 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1258 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 1259 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 1260 CHARACTER(LEN=default_string_length), & 1261 DIMENSION(:), POINTER :: Ainfo 1262 TYPE(charmm_info_type), POINTER :: chm_info 1263 TYPE(input_info_type), POINTER :: inp_info 1264 TYPE(gromos_info_type), POINTER :: gro_info 1265 TYPE(amber_info_type), POINTER :: amb_info 1266 INTEGER, INTENT(IN) :: iw 1267 1268 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_tors', & 1269 routineP = moduleN//':'//routineN 1270 1271 CHARACTER(LEN=default_string_length) :: ldum, name_atm_a, name_atm_b, & 1272 name_atm_c, name_atm_d 1273 INTEGER :: atm_a, atm_b, atm_c, atm_d, first, & 1274 handle2, i, imul, itype, j, k, last, & 1275 natom, ntorsion 1276 INTEGER, DIMENSION(:), POINTER :: molecule_list 1277 LOGICAL :: found, only_qm 1278 TYPE(atomic_kind_type), POINTER :: atomic_kind 1279 TYPE(molecule_kind_type), POINTER :: molecule_kind 1280 TYPE(molecule_type), POINTER :: molecule 1281 TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list 1282 1283 CALL timeset(routineN, handle2) 1284 1285 DO i = 1, SIZE(molecule_kind_set) 1286 molecule_kind => molecule_kind_set(i) 1287 CALL get_molecule_kind(molecule_kind=molecule_kind, & 1288 molecule_list=molecule_list, & 1289 natom=natom, & 1290 ntorsion=ntorsion, torsion_list=torsion_list) 1291 molecule => molecule_set(molecule_list(1)) 1292 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 1293 DO j = 1, ntorsion 1294 IF (torsion_list(j)%torsion_kind%id_type == do_ff_undef) THEN 1295 atm_a = torsion_list(j)%a 1296 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 1297 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1298 name=name_atm_a) 1299 atm_b = torsion_list(j)%b 1300 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 1301 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1302 name=name_atm_b) 1303 atm_c = torsion_list(j)%c 1304 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 1305 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1306 name=name_atm_c) 1307 atm_d = torsion_list(j)%d 1308 atomic_kind => particle_set(atm_d + first - 1)%atomic_kind 1309 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1310 name=name_atm_d) 1311 found = .FALSE. 1312 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c, id4=name_atm_d) 1313 CALL uppercase(name_atm_a) 1314 CALL uppercase(name_atm_b) 1315 CALL uppercase(name_atm_c) 1316 CALL uppercase(name_atm_d) 1317 1318 ! loop over params from GROMOS 1319 IF (ASSOCIATED(gro_info%torsion_k)) THEN 1320 k = SIZE(gro_info%torsion_k) 1321 itype = torsion_list(j)%itype 1322 IF (itype > 0) THEN 1323 CALL reallocate(torsion_list(j)%torsion_kind%k, 1, 1) 1324 CALL reallocate(torsion_list(j)%torsion_kind%m, 1, 1) 1325 CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, 1) 1326 torsion_list(j)%torsion_kind%nmul = 1 1327 torsion_list(j)%torsion_kind%m(1) = gro_info%torsion_m(itype) 1328 torsion_list(j)%torsion_kind%k(1) = gro_info%torsion_k(itype) 1329 torsion_list(j)%torsion_kind%phi0(1) = gro_info%torsion_phi0(itype) 1330 ELSE 1331 CALL reallocate(torsion_list(j)%torsion_kind%k, 1, 1) 1332 CALL reallocate(torsion_list(j)%torsion_kind%m, 1, 1) 1333 CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, 1) 1334 torsion_list(j)%torsion_kind%nmul = 1 1335 torsion_list(j)%torsion_kind%m(1) = gro_info%torsion_m(itype/k) 1336 torsion_list(j)%torsion_kind%k(1) = gro_info%torsion_k(itype/k) 1337 torsion_list(j)%torsion_kind%phi0(1) = gro_info%torsion_phi0(itype/k) 1338 END IF 1339 torsion_list(j)%torsion_kind%id_type = gro_info%ff_gromos_type 1340 torsion_list(j)%id_type = gro_info%ff_gromos_type 1341 found = .TRUE. 1342 imul = torsion_list(j)%torsion_kind%nmul 1343 END IF 1344 1345 ! loop over params from CHARMM 1346 IF (ASSOCIATED(chm_info%torsion_a)) THEN 1347 DO k = 1, SIZE(chm_info%torsion_a) 1348 IF ((((chm_info%torsion_a(k)) == (name_atm_a)) .AND. & 1349 ((chm_info%torsion_b(k)) == (name_atm_b)) .AND. & 1350 ((chm_info%torsion_c(k)) == (name_atm_c)) .AND. & 1351 ((chm_info%torsion_d(k)) == (name_atm_d))) .OR. & 1352 (((chm_info%torsion_a(k)) == (name_atm_d)) .AND. & 1353 ((chm_info%torsion_b(k)) == (name_atm_c)) .AND. & 1354 ((chm_info%torsion_c(k)) == (name_atm_b)) .AND. & 1355 ((chm_info%torsion_d(k)) == (name_atm_a)))) THEN 1356 imul = torsion_list(j)%torsion_kind%nmul + 1 1357 CALL reallocate(torsion_list(j)%torsion_kind%k, 1, imul) 1358 CALL reallocate(torsion_list(j)%torsion_kind%m, 1, imul) 1359 CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, imul) 1360 torsion_list(j)%torsion_kind%id_type = do_ff_charmm 1361 torsion_list(j)%torsion_kind%k(imul) = chm_info%torsion_k(k) 1362 torsion_list(j)%torsion_kind%m(imul) = chm_info%torsion_m(k) 1363 torsion_list(j)%torsion_kind%phi0(imul) = chm_info%torsion_phi0(k) 1364 torsion_list(j)%torsion_kind%nmul = imul 1365 found = .TRUE. 1366 END IF 1367 END DO 1368 1369 IF (.NOT. found) THEN 1370 DO k = 1, SIZE(chm_info%torsion_a) 1371 IF ((((chm_info%torsion_a(k)) == ("X")) .AND. & 1372 ((chm_info%torsion_b(k)) == (name_atm_b)) .AND. & 1373 ((chm_info%torsion_c(k)) == (name_atm_c)) .AND. & 1374 ((chm_info%torsion_d(k)) == ("X"))) .OR. & 1375 (((chm_info%torsion_a(k)) == ("X")) .AND. & 1376 ((chm_info%torsion_b(k)) == (name_atm_c)) .AND. & 1377 ((chm_info%torsion_c(k)) == (name_atm_b)) .AND. & 1378 ((chm_info%torsion_d(k)) == ("X")))) THEN 1379 imul = torsion_list(j)%torsion_kind%nmul + 1 1380 CALL reallocate(torsion_list(j)%torsion_kind%k, 1, imul) 1381 CALL reallocate(torsion_list(j)%torsion_kind%m, 1, imul) 1382 CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, imul) 1383 torsion_list(j)%torsion_kind%id_type = do_ff_charmm 1384 torsion_list(j)%torsion_kind%k(imul) = chm_info%torsion_k(k) 1385 torsion_list(j)%torsion_kind%m(imul) = chm_info%torsion_m(k) 1386 torsion_list(j)%torsion_kind%phi0(imul) = chm_info%torsion_phi0(k) 1387 torsion_list(j)%torsion_kind%nmul = imul 1388 found = .TRUE. 1389 END IF 1390 END DO 1391 END IF 1392 END IF 1393 1394 ! loop over params from AMBER 1395 IF (ASSOCIATED(amb_info%torsion_a)) THEN 1396 DO k = 1, SIZE(amb_info%torsion_a) 1397 IF ((((amb_info%torsion_a(k)) == (name_atm_a)) .AND. & 1398 ((amb_info%torsion_b(k)) == (name_atm_b)) .AND. & 1399 ((amb_info%torsion_c(k)) == (name_atm_c)) .AND. & 1400 ((amb_info%torsion_d(k)) == (name_atm_d))) .OR. & 1401 (((amb_info%torsion_a(k)) == (name_atm_d)) .AND. & 1402 ((amb_info%torsion_b(k)) == (name_atm_c)) .AND. & 1403 ((amb_info%torsion_c(k)) == (name_atm_b)) .AND. & 1404 ((amb_info%torsion_d(k)) == (name_atm_a)))) THEN 1405 imul = torsion_list(j)%torsion_kind%nmul + 1 1406 CALL reallocate(torsion_list(j)%torsion_kind%k, 1, imul) 1407 CALL reallocate(torsion_list(j)%torsion_kind%m, 1, imul) 1408 CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, imul) 1409 torsion_list(j)%torsion_kind%id_type = do_ff_amber 1410 torsion_list(j)%torsion_kind%k(imul) = amb_info%torsion_k(k) 1411 torsion_list(j)%torsion_kind%m(imul) = amb_info%torsion_m(k) 1412 torsion_list(j)%torsion_kind%phi0(imul) = amb_info%torsion_phi0(k) 1413 torsion_list(j)%torsion_kind%nmul = imul 1414 found = .TRUE. 1415 END IF 1416 END DO 1417 1418 IF (.NOT. found) THEN 1419 DO k = 1, SIZE(amb_info%torsion_a) 1420 IF ((((amb_info%torsion_a(k)) == ("X")) .AND. & 1421 ((amb_info%torsion_b(k)) == (name_atm_b)) .AND. & 1422 ((amb_info%torsion_c(k)) == (name_atm_c)) .AND. & 1423 ((amb_info%torsion_d(k)) == ("X"))) .OR. & 1424 (((amb_info%torsion_a(k)) == ("X")) .AND. & 1425 ((amb_info%torsion_b(k)) == (name_atm_c)) .AND. & 1426 ((amb_info%torsion_c(k)) == (name_atm_b)) .AND. & 1427 ((amb_info%torsion_d(k)) == ("X")))) THEN 1428 imul = torsion_list(j)%torsion_kind%nmul + 1 1429 CALL reallocate(torsion_list(j)%torsion_kind%k, 1, imul) 1430 CALL reallocate(torsion_list(j)%torsion_kind%m, 1, imul) 1431 CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, imul) 1432 torsion_list(j)%torsion_kind%id_type = do_ff_amber 1433 torsion_list(j)%torsion_kind%k(imul) = amb_info%torsion_k(k) 1434 torsion_list(j)%torsion_kind%m(imul) = amb_info%torsion_m(k) 1435 torsion_list(j)%torsion_kind%phi0(imul) = amb_info%torsion_phi0(k) 1436 torsion_list(j)%torsion_kind%nmul = imul 1437 found = .TRUE. 1438 END IF 1439 END DO 1440 END IF 1441 END IF 1442 1443 ! always have the input param last to overwrite all the other ones 1444 IF (ASSOCIATED(inp_info%torsion_a)) THEN 1445 DO k = 1, SIZE(inp_info%torsion_a) 1446 IF ((((inp_info%torsion_a(k)) == (name_atm_a)) .AND. & 1447 ((inp_info%torsion_b(k)) == (name_atm_b)) .AND. & 1448 ((inp_info%torsion_c(k)) == (name_atm_c)) .AND. & 1449 ((inp_info%torsion_d(k)) == (name_atm_d))) .OR. & 1450 (((inp_info%torsion_a(k)) == (name_atm_d)) .AND. & 1451 ((inp_info%torsion_b(k)) == (name_atm_c)) .AND. & 1452 ((inp_info%torsion_c(k)) == (name_atm_b)) .AND. & 1453 ((inp_info%torsion_d(k)) == (name_atm_a)))) THEN 1454 imul = torsion_list(j)%torsion_kind%nmul + 1 1455 CALL reallocate(torsion_list(j)%torsion_kind%k, 1, imul) 1456 CALL reallocate(torsion_list(j)%torsion_kind%m, 1, imul) 1457 CALL reallocate(torsion_list(j)%torsion_kind%phi0, 1, imul) 1458 torsion_list(j)%torsion_kind%id_type = inp_info%torsion_kind(k) 1459 torsion_list(j)%torsion_kind%k(imul) = inp_info%torsion_k(k) 1460 torsion_list(j)%torsion_kind%m(imul) = inp_info%torsion_m(k) 1461 torsion_list(j)%torsion_kind%phi0(imul) = inp_info%torsion_phi0(k) 1462 torsion_list(j)%torsion_kind%nmul = imul 1463 found = .TRUE. 1464 END IF 1465 END DO 1466 END IF 1467 1468 IF (.NOT. found) THEN 1469 CALL store_FF_missing_par(atm1=TRIM(name_atm_a), & 1470 atm2=TRIM(name_atm_b), & 1471 atm3=TRIM(name_atm_c), & 1472 atm4=TRIM(name_atm_d), & 1473 type_name="Torsion", & 1474 array=Ainfo) 1475 torsion_list(j)%torsion_kind%id_type = do_ff_undef 1476 torsion_list(j)%id_type = do_ff_undef 1477 ELSE 1478 ldum = cp_to_string(imul) 1479 IF ((imul /= 1) .AND. (iw > 0)) & 1480 WRITE (iw, '(/,2("UTIL_INFO| ",A,/))') & 1481 "Multiple Torsion declarations: "//TRIM(name_atm_a)// & 1482 ","//TRIM(name_atm_b)//","//TRIM(name_atm_c)//" and "//TRIM(name_atm_d), & 1483 "Number of defined torsions: "//TRIM(ldum)//" ." 1484 END IF 1485 ! 1486 ! QM/MM modifications 1487 ! 1488 IF (only_qm) THEN 1489 IF (iw > 0) WRITE (iw, *) " Torsion PARAM between QM atoms ", j, " : ", & 1490 TRIM(name_atm_a), " ", & 1491 TRIM(name_atm_b), " ", & 1492 TRIM(name_atm_c), " ", & 1493 TRIM(name_atm_d), " ", & 1494 torsion_list(j)%a, & 1495 torsion_list(j)%b, & 1496 torsion_list(j)%c, & 1497 torsion_list(j)%d 1498 torsion_list(j)%torsion_kind%id_type = do_ff_undef 1499 torsion_list(j)%id_type = do_ff_undef 1500 END IF 1501 END IF 1502 END DO 1503 CALL set_molecule_kind(molecule_kind=molecule_kind, & 1504 torsion_list=torsion_list) 1505 END DO 1506 CALL timestop(handle2) 1507 1508 END SUBROUTINE force_field_pack_tors 1509 1510! ************************************************************************************************** 1511!> \brief Pack in impropers information needed for the force_field 1512!> \param particle_set ... 1513!> \param molecule_kind_set ... 1514!> \param molecule_set ... 1515!> \param Ainfo ... 1516!> \param chm_info ... 1517!> \param inp_info ... 1518!> \param gro_info ... 1519! ************************************************************************************************** 1520 SUBROUTINE force_field_pack_impr(particle_set, molecule_kind_set, molecule_set, & 1521 Ainfo, chm_info, inp_info, gro_info) 1522 1523 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1524 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 1525 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 1526 CHARACTER(LEN=default_string_length), & 1527 DIMENSION(:), POINTER :: Ainfo 1528 TYPE(charmm_info_type), POINTER :: chm_info 1529 TYPE(input_info_type), POINTER :: inp_info 1530 TYPE(gromos_info_type), POINTER :: gro_info 1531 1532 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_impr', & 1533 routineP = moduleN//':'//routineN 1534 1535 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, name_atm_c, & 1536 name_atm_d 1537 INTEGER :: atm_a, atm_b, atm_c, atm_d, first, & 1538 handle2, i, itype, j, k, last, natom, & 1539 nimpr 1540 INTEGER, DIMENSION(:), POINTER :: molecule_list 1541 LOGICAL :: found, only_qm 1542 TYPE(atomic_kind_type), POINTER :: atomic_kind 1543 TYPE(impr_type), DIMENSION(:), POINTER :: impr_list 1544 TYPE(molecule_kind_type), POINTER :: molecule_kind 1545 TYPE(molecule_type), POINTER :: molecule 1546 1547 CALL timeset(routineN, handle2) 1548 1549 DO i = 1, SIZE(molecule_kind_set) 1550 molecule_kind => molecule_kind_set(i) 1551 CALL get_molecule_kind(molecule_kind=molecule_kind, & 1552 molecule_list=molecule_list, & 1553 natom=natom, & 1554 nimpr=nimpr, impr_list=impr_list) 1555 molecule => molecule_set(molecule_list(1)) 1556 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 1557 DO j = 1, nimpr 1558 atm_a = impr_list(j)%a 1559 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 1560 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1561 name=name_atm_a) 1562 atm_b = impr_list(j)%b 1563 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 1564 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1565 name=name_atm_b) 1566 atm_c = impr_list(j)%c 1567 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 1568 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1569 name=name_atm_c) 1570 atm_d = impr_list(j)%d 1571 atomic_kind => particle_set(atm_d + first - 1)%atomic_kind 1572 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1573 name=name_atm_d) 1574 found = .FALSE. 1575 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c, id4=name_atm_d) 1576 CALL uppercase(name_atm_a) 1577 CALL uppercase(name_atm_b) 1578 CALL uppercase(name_atm_c) 1579 CALL uppercase(name_atm_d) 1580 1581 ! loop over params from GROMOS 1582 IF (ASSOCIATED(gro_info%impr_k)) THEN 1583 k = SIZE(gro_info%impr_k) 1584 itype = impr_list(j)%itype 1585 IF (itype > 0) THEN 1586 impr_list(j)%impr_kind%k = gro_info%impr_k(itype) 1587 impr_list(j)%impr_kind%phi0 = gro_info%impr_phi0(itype) 1588 ELSE 1589 impr_list(j)%impr_kind%k = gro_info%impr_k(itype) 1590 impr_list(j)%impr_kind%phi0 = gro_info%impr_phi0(itype) 1591 END IF 1592 found = .TRUE. 1593 impr_list(j)%impr_kind%id_type = gro_info%ff_gromos_type 1594 impr_list(j)%id_type = gro_info%ff_gromos_type 1595 END IF 1596 1597 ! loop over params from CHARMM 1598 IF (ASSOCIATED(chm_info%impr_a)) THEN 1599 DO k = 1, SIZE(chm_info%impr_a) 1600 IF ((((chm_info%impr_a(k)) == (name_atm_a)) .AND. & 1601 ((chm_info%impr_b(k)) == (name_atm_b)) .AND. & 1602 ((chm_info%impr_c(k)) == (name_atm_c)) .AND. & 1603 ((chm_info%impr_d(k)) == (name_atm_d))) .OR. & 1604 (((chm_info%impr_a(k)) == (name_atm_d)) .AND. & 1605 ((chm_info%impr_b(k)) == (name_atm_c)) .AND. & 1606 ((chm_info%impr_c(k)) == (name_atm_b)) .AND. & 1607 ((chm_info%impr_d(k)) == (name_atm_a)))) THEN 1608 impr_list(j)%impr_kind%id_type = do_ff_charmm 1609 impr_list(j)%impr_kind%k = chm_info%impr_k(k) 1610 impr_list(j)%impr_kind%phi0 = chm_info%impr_phi0(k) 1611 CALL issue_duplications(found, "Impropers", name_atm_a, name_atm_b, & 1612 name_atm_c, name_atm_d) 1613 found = .TRUE. 1614 EXIT 1615 END IF 1616 END DO 1617 IF (.NOT. found) THEN 1618 DO k = 1, SIZE(chm_info%impr_a) 1619 IF ((((chm_info%impr_a(k)) == (name_atm_a)) .AND. & 1620 ((chm_info%impr_b(k)) == ("X")) .AND. & 1621 ((chm_info%impr_c(k)) == ("X")) .AND. & 1622 ((chm_info%impr_d(k)) == (name_atm_d))) .OR. & 1623 (((chm_info%impr_a(k)) == (name_atm_d)) .AND. & 1624 ((chm_info%impr_b(k)) == ("X")) .AND. & 1625 ((chm_info%impr_c(k)) == ("X")) .AND. & 1626 ((chm_info%impr_d(k)) == (name_atm_a)))) THEN 1627 impr_list(j)%impr_kind%id_type = do_ff_charmm 1628 impr_list(j)%impr_kind%k = chm_info%impr_k(k) 1629 impr_list(j)%impr_kind%phi0 = chm_info%impr_phi0(k) 1630 CALL issue_duplications(found, "Impropers", name_atm_a, name_atm_b, & 1631 name_atm_c, name_atm_d) 1632 found = .TRUE. 1633 EXIT 1634 END IF 1635 END DO 1636 END IF 1637 END IF 1638 1639 ! loop over params from AMBER not needed since impropers in AMBER 1640 ! are treated like standard torsions 1641 1642 ! always have the input param last to overwrite all the other ones 1643 IF (ASSOCIATED(inp_info%impr_a)) THEN 1644 DO k = 1, SIZE(inp_info%impr_a) 1645 IF (((inp_info%impr_a(k)) == (name_atm_a)) .AND. & 1646 ((inp_info%impr_b(k)) == (name_atm_b)) .AND. & 1647 ((((inp_info%impr_c(k)) == (name_atm_c)) .AND. & 1648 ((inp_info%impr_d(k)) == (name_atm_d))) .OR. & 1649 (((inp_info%impr_c(k)) == (name_atm_d)) .AND. & 1650 ((inp_info%impr_d(k)) == (name_atm_c))))) THEN 1651 impr_list(j)%impr_kind%id_type = inp_info%impr_kind(k) 1652 impr_list(j)%impr_kind%k = inp_info%impr_k(k) 1653 IF (((inp_info%impr_c(k)) == (name_atm_c)) .AND. & 1654 ((inp_info%impr_d(k)) == (name_atm_d))) THEN 1655 impr_list(j)%impr_kind%phi0 = inp_info%impr_phi0(k) 1656 ELSE 1657 impr_list(j)%impr_kind%phi0 = -inp_info%impr_phi0(k) 1658 ! alternative solutions: 1659 ! - swap impr_list(j)%c with impr_list(j)%d and 1660 ! name_atom_c with name_atom_d and 1661 ! atm_c with atm_d 1662 ! - introduce impr_list(j)%impr_kind%sign. if one, the 1663 ! sign of phi is not changed in mol_force.f90. if minus 1664 ! one, the sign of phi is changed in mol_force.f90 1665 ! similar problems with parameters from charmm pot file 1666 ! above 1667 END IF 1668 CALL issue_duplications(found, "Impropers", name_atm_a, name_atm_b, & 1669 name_atm_c, name_atm_d) 1670 found = .TRUE. 1671 EXIT 1672 END IF 1673 END DO 1674 END IF 1675 1676 IF (.NOT. found) THEN 1677 CALL store_FF_missing_par(atm1=TRIM(name_atm_a), & 1678 atm2=TRIM(name_atm_b), & 1679 atm3=TRIM(name_atm_c), & 1680 atm4=TRIM(name_atm_d), & 1681 type_name="Improper", & 1682 array=Ainfo) 1683 impr_list(j)%impr_kind%k = 0.0_dp 1684 impr_list(j)%impr_kind%phi0 = 0.0_dp 1685 impr_list(j)%impr_kind%id_type = do_ff_undef 1686 impr_list(j)%id_type = do_ff_undef 1687 END IF 1688 ! 1689 ! QM/MM modifications 1690 ! 1691 IF (only_qm) THEN 1692 impr_list(j)%impr_kind%id_type = do_ff_undef 1693 impr_list(j)%id_type = do_ff_undef 1694 END IF 1695 1696 END DO 1697 CALL set_molecule_kind(molecule_kind=molecule_kind, impr_list=impr_list) 1698 END DO 1699 CALL timestop(handle2) 1700 1701 END SUBROUTINE force_field_pack_impr 1702 1703! ************************************************************************************************** 1704!> \brief Pack in opbend information needed for the force_field. 1705!> No loop over params for charmm, amber and gromos since these force 1706!> fields have no opbends 1707!> \param particle_set ... 1708!> \param molecule_kind_set ... 1709!> \param molecule_set ... 1710!> \param Ainfo ... 1711!> \param inp_info ... 1712!> \author Louis Vanduyfhuys 1713! ************************************************************************************************** 1714 SUBROUTINE force_field_pack_opbend(particle_set, molecule_kind_set, molecule_set, & 1715 Ainfo, inp_info) 1716 1717 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1718 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 1719 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 1720 CHARACTER(LEN=default_string_length), & 1721 DIMENSION(:), POINTER :: Ainfo 1722 TYPE(input_info_type), POINTER :: inp_info 1723 1724 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_opbend', & 1725 routineP = moduleN//':'//routineN 1726 1727 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_b, name_atm_c, & 1728 name_atm_d 1729 INTEGER :: atm_a, atm_b, atm_c, atm_d, first, & 1730 handle2, i, j, k, last, natom, nopbend 1731 INTEGER, DIMENSION(:), POINTER :: molecule_list 1732 LOGICAL :: found, only_qm 1733 TYPE(atomic_kind_type), POINTER :: atomic_kind 1734 TYPE(molecule_kind_type), POINTER :: molecule_kind 1735 TYPE(molecule_type), POINTER :: molecule 1736 TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list 1737 1738 CALL timeset(routineN, handle2) 1739 1740 DO i = 1, SIZE(molecule_kind_set) 1741 molecule_kind => molecule_kind_set(i) 1742 CALL get_molecule_kind(molecule_kind=molecule_kind, & 1743 molecule_list=molecule_list, & 1744 natom=natom, & 1745 nopbend=nopbend, opbend_list=opbend_list) 1746 molecule => molecule_set(molecule_list(1)) 1747 1748 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 1749 DO j = 1, nopbend 1750 atm_a = opbend_list(j)%a 1751 atomic_kind => particle_set(atm_a + first - 1)%atomic_kind 1752 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1753 name=name_atm_a) 1754 atm_b = opbend_list(j)%b 1755 atomic_kind => particle_set(atm_b + first - 1)%atomic_kind 1756 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1757 name=name_atm_b) 1758 atm_c = opbend_list(j)%c 1759 atomic_kind => particle_set(atm_c + first - 1)%atomic_kind 1760 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1761 name=name_atm_c) 1762 atm_d = opbend_list(j)%d 1763 atomic_kind => particle_set(atm_d + first - 1)%atomic_kind 1764 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1765 name=name_atm_d) 1766 found = .FALSE. 1767 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b, id3=name_atm_c, id4=name_atm_d) 1768 CALL uppercase(name_atm_a) 1769 CALL uppercase(name_atm_b) 1770 CALL uppercase(name_atm_c) 1771 CALL uppercase(name_atm_d) 1772 1773 ! always have the input param last to overwrite all the other ones 1774 IF (ASSOCIATED(inp_info%opbend_a)) THEN 1775 DO k = 1, SIZE(inp_info%opbend_a) 1776 IF (((inp_info%opbend_a(k)) == (name_atm_a)) .AND. & 1777 ((inp_info%opbend_d(k)) == (name_atm_d)) .AND. & 1778 ((((inp_info%opbend_c(k)) == (name_atm_c)) .AND. & 1779 ((inp_info%opbend_b(k)) == (name_atm_b))) .OR. & 1780 (((inp_info%opbend_c(k)) == (name_atm_b)) .AND. & 1781 ((inp_info%opbend_b(k)) == (name_atm_c))))) THEN 1782 opbend_list(j)%opbend_kind%id_type = inp_info%opbend_kind(k) 1783 opbend_list(j)%opbend_kind%k = inp_info%opbend_k(k) 1784 IF (((inp_info%opbend_c(k)) == (name_atm_c)) .AND. & 1785 ((inp_info%opbend_b(k)) == (name_atm_b))) THEN 1786 opbend_list(j)%opbend_kind%phi0 = inp_info%opbend_phi0(k) 1787 ELSE 1788 opbend_list(j)%opbend_kind%phi0 = -inp_info%opbend_phi0(k) 1789 ! alternative solutions: 1790 ! - swap opbend_list(j)%c with opbend_list(j)%b and 1791 ! name_atom_c with name_atom_b and 1792 ! atm_c with atm_b 1793 ! - introduce opbend_list(j)%opbend_kind%sign. if one, the 1794 ! sign of phi is not changed in mol_force.f90. if minus 1795 ! one, the sign of phi is changed in mol_force.f90 1796 1797 END IF 1798 CALL issue_duplications(found, "Out of plane bend", name_atm_a, name_atm_b, & 1799 name_atm_c, name_atm_d) 1800 found = .TRUE. 1801 EXIT 1802 END IF 1803 END DO 1804 END IF 1805 1806 IF (.NOT. found) THEN 1807 CALL store_FF_missing_par(atm1=TRIM(name_atm_a), & 1808 atm2=TRIM(name_atm_b), & 1809 atm3=TRIM(name_atm_c), & 1810 atm4=TRIM(name_atm_d), & 1811 type_name="Out of plane bend", & 1812 array=Ainfo) 1813 opbend_list(j)%opbend_kind%k = 0.0_dp 1814 opbend_list(j)%opbend_kind%phi0 = 0.0_dp 1815 opbend_list(j)%opbend_kind%id_type = do_ff_undef 1816 opbend_list(j)%id_type = do_ff_undef 1817 END IF 1818 ! 1819 ! QM/MM modifications 1820 ! 1821 IF (only_qm) THEN 1822 opbend_list(j)%opbend_kind%id_type = do_ff_undef 1823 opbend_list(j)%id_type = do_ff_undef 1824 END IF 1825 1826 END DO 1827 CALL set_molecule_kind(molecule_kind=molecule_kind, opbend_list=opbend_list) 1828 END DO 1829 CALL timestop(handle2) 1830 1831 END SUBROUTINE force_field_pack_opbend 1832 1833! ************************************************************************************************** 1834!> \brief Set up array of full charges 1835!> \param charges ... 1836!> \param charges_section ... 1837!> \param particle_set ... 1838!> \param my_qmmm ... 1839!> \param qmmm_env ... 1840!> \param inp_info ... 1841!> \param iw4 ... 1842!> \date 12.2010 1843!> \author Teodoro Laino (teodoro.laino@gmail.com) 1844! ************************************************************************************************** 1845 SUBROUTINE force_field_pack_charges(charges, charges_section, particle_set, & 1846 my_qmmm, qmmm_env, inp_info, iw4) 1847 1848 REAL(KIND=dp), DIMENSION(:), POINTER :: charges 1849 TYPE(section_vals_type), POINTER :: charges_section 1850 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1851 LOGICAL :: my_qmmm 1852 TYPE(qmmm_env_mm_type), POINTER :: qmmm_env 1853 TYPE(input_info_type), POINTER :: inp_info 1854 INTEGER :: iw4 1855 1856 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_charges', & 1857 routineP = moduleN//':'//routineN 1858 1859 CHARACTER(LEN=default_string_length) :: atmname 1860 INTEGER :: handle, iatom, ilink, j, nval 1861 LOGICAL :: found_p, is_link_atom, is_ok, & 1862 only_manybody, only_qm 1863 REAL(KIND=dp) :: charge, charge_tot, rval, scale_factor 1864 TYPE(atomic_kind_type), POINTER :: atomic_kind 1865 TYPE(cp_sll_val_type), POINTER :: list 1866 TYPE(fist_potential_type), POINTER :: fist_potential 1867 TYPE(val_type), POINTER :: val 1868 1869 CALL timeset(routineN, handle) 1870 charge_tot = 0.0_dp 1871 NULLIFY (list) 1872 1873 ! Not implemented for core-shell 1874 IF (ASSOCIATED(inp_info%shell_list)) THEN 1875 CPABORT("Array of charges not implemented for CORE-SHELL model!!") 1876 END IF 1877 1878 ! Allocate array to particle_set size 1879 CPASSERT(.NOT. (ASSOCIATED(charges))) 1880 ALLOCATE (charges(SIZE(particle_set))) 1881 1882 ! Fill with input values 1883 CALL section_vals_val_get(charges_section, "_DEFAULT_KEYWORD_", n_rep_val=nval) 1884 CPASSERT(nval == SIZE(charges)) 1885 CALL section_vals_list_get(charges_section, "_DEFAULT_KEYWORD_", list=list) 1886 DO iatom = 1, nval 1887 ! we use only the first default_string_length characters of each line 1888 is_ok = cp_sll_val_next(list, val) 1889 CALL val_get(val, r_val=rval) 1890 ! assign values 1891 charges(iatom) = rval 1892 1893 ! Perform a post-processing 1894 atomic_kind => particle_set(iatom)%atomic_kind 1895 CALL get_atomic_kind(atomic_kind=atomic_kind, & 1896 fist_potential=fist_potential, & 1897 name=atmname) 1898 CALL get_potential(potential=fist_potential, qeff=charge) 1899 1900 only_qm = qmmm_ff_precond_only_qm(id1=atmname, is_link=is_link_atom) 1901 CALL uppercase(atmname) 1902 IF (charge /= -HUGE(0.0_dp)) & 1903 CALL cp_warn(__LOCATION__, & 1904 "The charge for atom index ("//cp_to_string(iatom)//") and atom name ("// & 1905 TRIM(atmname)//") was already defined. The charge associated to this kind"// & 1906 " will be set to an uninitialized value and only the atom specific charge will be used! ") 1907 charge = -HUGE(0.0_dp) 1908 1909 ! Check if the potential really requires the charge definition.. 1910 IF (ASSOCIATED(inp_info%nonbonded)) THEN 1911 IF (ASSOCIATED(inp_info%nonbonded%pot)) THEN 1912 ! Let's find the nonbonded potential where this atom is involved 1913 only_manybody = .TRUE. 1914 found_p = .FALSE. 1915 DO j = 1, SIZE(inp_info%nonbonded%pot) 1916 IF (atmname == inp_info%nonbonded%pot(j)%pot%at1 .OR. & 1917 atmname == inp_info%nonbonded%pot(j)%pot%at2) THEN 1918 SELECT CASE (inp_info%nonbonded%pot(j)%pot%type(1)) 1919 CASE (ea_type, tersoff_type, siepmann_type) 1920 ! Charge is zero for EAM, TERSOFF and SIEPMANN type potential 1921 ! Do nothing.. 1922 CASE DEFAULT 1923 only_manybody = .FALSE. 1924 EXIT 1925 END SELECT 1926 found_p = .TRUE. 1927 END IF 1928 END DO 1929 IF (only_manybody .AND. found_p) THEN 1930 charges(iatom) = 0.0_dp 1931 END IF 1932 END IF 1933 END IF 1934 ! 1935 ! QM/MM modifications 1936 ! 1937 IF (only_qm .AND. my_qmmm) THEN 1938 IF (qmmm_env%qmmm_coupl_type /= do_qmmm_none) THEN 1939 scale_factor = 0.0_dp 1940 IF (is_link_atom) THEN 1941 ! 1942 ! Find the scaling factor... 1943 ! 1944 DO ilink = 1, SIZE(qmmm_env%mm_link_atoms) 1945 IF (iatom == qmmm_env%mm_link_atoms(ilink)) EXIT 1946 END DO 1947 CPASSERT(ilink <= SIZE(qmmm_env%mm_link_atoms)) 1948 scale_factor = qmmm_env%fist_scale_charge_link(ilink) 1949 END IF 1950 charges(iatom) = charges(iatom)*scale_factor 1951 END IF 1952 END IF 1953 END DO 1954 ! Sum up total charges for IO 1955 charge_tot = SUM(charges) 1956 ! Print Total Charge of the system 1957 IF (iw4 > 0) THEN 1958 WRITE (iw4, FMT='(T2,"CHARGE_INFO| Total Charge of the Classical System: ",T69,F12.6)') charge_tot 1959 END IF 1960 CALL timestop(handle) 1961 END SUBROUTINE force_field_pack_charges 1962 1963! ************************************************************************************************** 1964!> \brief Set up atomic_kind_set()%fist_potentail%[qeff] 1965!> and shell potential parameters 1966!> \param atomic_kind_set ... 1967!> \param qmmm_env ... 1968!> \param fatal ... 1969!> \param iw ... 1970!> \param iw4 ... 1971!> \param Ainfo ... 1972!> \param my_qmmm ... 1973!> \param inp_info ... 1974! ************************************************************************************************** 1975 SUBROUTINE force_field_pack_charge(atomic_kind_set, qmmm_env, fatal, iw, iw4, & 1976 Ainfo, my_qmmm, inp_info) 1977 1978 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1979 TYPE(qmmm_env_mm_type), POINTER :: qmmm_env 1980 LOGICAL :: fatal 1981 INTEGER :: iw, iw4 1982 CHARACTER(LEN=default_string_length), & 1983 DIMENSION(:), POINTER :: Ainfo 1984 LOGICAL :: my_qmmm 1985 TYPE(input_info_type), POINTER :: inp_info 1986 1987 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_charge', & 1988 routineP = moduleN//':'//routineN 1989 1990 CHARACTER(LEN=default_string_length) :: atmname 1991 INTEGER :: handle, i, ilink, j 1992 INTEGER, DIMENSION(:), POINTER :: my_atom_list 1993 LOGICAL :: found, found_p, is_link_atom, is_shell, & 1994 only_manybody, only_qm 1995 REAL(KIND=dp) :: charge, charge_tot, cs_charge, & 1996 scale_factor 1997 TYPE(atomic_kind_type), POINTER :: atomic_kind 1998 TYPE(fist_potential_type), POINTER :: fist_potential 1999 2000 CALL timeset(routineN, handle) 2001 charge_tot = 0.0_dp 2002 DO i = 1, SIZE(atomic_kind_set) 2003 atomic_kind => atomic_kind_set(i) 2004 CALL get_atomic_kind(atomic_kind=atomic_kind, & 2005 fist_potential=fist_potential, & 2006 atom_list=my_atom_list, & 2007 name=atmname) 2008 CALL get_potential(potential=fist_potential, qeff=charge) 2009 2010 is_shell = .FALSE. 2011 found = .FALSE. 2012 only_qm = qmmm_ff_precond_only_qm(id1=atmname, is_link=is_link_atom) 2013 CALL uppercase(atmname) 2014 IF (charge /= -HUGE(0.0_dp)) found = .TRUE. 2015 2016 ! Always have the input param last to overwrite all the other ones 2017 IF (ASSOCIATED(inp_info%charge_atm)) THEN 2018 DO j = 1, SIZE(inp_info%charge_atm) 2019 IF (iw > 0) WRITE (iw, *) "Charge Checking ::", TRIM(inp_info%charge_atm(j)), atmname 2020 IF ((inp_info%charge_atm(j)) == atmname) THEN 2021 charge = inp_info%charge(j) 2022 CALL issue_duplications(found, "Charge", atmname) 2023 found = .TRUE. 2024 END IF 2025 END DO 2026 END IF 2027 ! Check if the ATOM type has a core-shell associated.. In this case 2028 ! print a warning: the CHARGE will not be used if defined.. or we can 2029 ! even skip its definition.. 2030 IF (ASSOCIATED(inp_info%shell_list)) THEN 2031 DO j = 1, SIZE(inp_info%shell_list) 2032 IF ((inp_info%shell_list(j)%atm_name) == atmname) THEN 2033 is_shell = .TRUE. 2034 cs_charge = inp_info%shell_list(j)%shell%charge_core + & 2035 inp_info%shell_list(j)%shell%charge_shell 2036 charge = 0.0_dp 2037 IF (found) THEN 2038 IF (found) & 2039 CALL cp_warn(__LOCATION__, & 2040 "CORE-SHELL model defined for KIND ("//TRIM(atmname)//")"// & 2041 " ignoring charge definition! ") 2042 ELSE 2043 found = .TRUE. 2044 END IF 2045 END IF 2046 END DO 2047 END IF 2048 ! Check if the potential really requires the charge definition.. 2049 IF (ASSOCIATED(inp_info%nonbonded)) THEN 2050 IF (ASSOCIATED(inp_info%nonbonded%pot)) THEN 2051 ! Let's find the nonbonded potential where this atom is involved 2052 only_manybody = .TRUE. 2053 found_p = .FALSE. 2054 DO j = 1, SIZE(inp_info%nonbonded%pot) 2055 IF (atmname == inp_info%nonbonded%pot(j)%pot%at1 .OR. & 2056 atmname == inp_info%nonbonded%pot(j)%pot%at2) THEN 2057 SELECT CASE (inp_info%nonbonded%pot(j)%pot%type(1)) 2058 CASE (ea_type, tersoff_type, siepmann_type, quip_type) 2059 ! Charge is zero for EAM, TERSOFF and SIEPMANN type potential 2060 ! Do nothing.. 2061 CASE DEFAULT 2062 only_manybody = .FALSE. 2063 EXIT 2064 END SELECT 2065 found_p = .TRUE. 2066 END IF 2067 END DO 2068 IF (only_manybody .AND. found_p) THEN 2069 charge = 0.0_dp 2070 found = .TRUE. 2071 END IF 2072 END IF 2073 END IF 2074 IF (.NOT. found) THEN 2075 ! Set the charge to zero anyway in case the user decides to ignore 2076 ! missing critical parameters. 2077 charge = 0.0_dp 2078 CALL store_FF_missing_par(atm1=TRIM(atmname), & 2079 fatal=fatal, & 2080 type_name="Charge", & 2081 array=Ainfo) 2082 END IF 2083 ! 2084 ! QM/MM modifications 2085 ! 2086 IF (only_qm .AND. my_qmmm) THEN 2087 IF (qmmm_env%qmmm_coupl_type /= do_qmmm_none) THEN 2088 scale_factor = 0.0_dp 2089 IF (is_link_atom) THEN 2090 ! 2091 ! Find the scaling factor... 2092 ! 2093 DO ilink = 1, SIZE(qmmm_env%mm_link_atoms) 2094 IF (ANY(my_atom_list == qmmm_env%mm_link_atoms(ilink))) EXIT 2095 END DO 2096 CPASSERT(ilink <= SIZE(qmmm_env%mm_link_atoms)) 2097 scale_factor = qmmm_env%fist_scale_charge_link(ilink) 2098 END IF 2099 charge = charge*scale_factor 2100 END IF 2101 END IF 2102 2103 CALL set_potential(potential=fist_potential, qeff=charge) 2104 ! Sum up total charges for IO 2105 IF (found) THEN 2106 IF (is_shell) THEN 2107 charge_tot = charge_tot + atomic_kind%natom*cs_charge 2108 ELSE 2109 charge_tot = charge_tot + atomic_kind%natom*charge 2110 END IF 2111 END IF 2112 END DO 2113 ! Print Total Charge of the system 2114 IF (iw4 > 0) THEN 2115 WRITE (iw4, FMT='(T2,"CHARGE_INFO| Total Charge of the Classical System: ",T69,F12.6)') charge_tot 2116 END IF 2117 CALL timestop(handle) 2118 END SUBROUTINE force_field_pack_charge 2119 2120! ************************************************************************************************** 2121!> \brief Set up the radius of the electrostatic multipole in Fist 2122!> \param atomic_kind_set ... 2123!> \param iw ... 2124!> \param subsys_section ... 2125!> \author Toon.Verstraelen@gmail.com 2126! ************************************************************************************************** 2127 SUBROUTINE force_field_pack_radius(atomic_kind_set, iw, subsys_section) 2128 2129 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 2130 INTEGER, INTENT(IN) :: iw 2131 TYPE(section_vals_type), POINTER :: subsys_section 2132 2133 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_radius', & 2134 routineP = moduleN//':'//routineN 2135 2136 CHARACTER(LEN=default_string_length) :: inp_kind_name, kind_name 2137 INTEGER :: handle, i, i_rep, n_rep 2138 LOGICAL :: found 2139 REAL(KIND=dp) :: mm_radius 2140 TYPE(atomic_kind_type), POINTER :: atomic_kind 2141 TYPE(fist_potential_type), POINTER :: fist_potential 2142 TYPE(section_vals_type), POINTER :: kind_section 2143 2144 CALL timeset(routineN, handle) 2145 2146 kind_section => section_vals_get_subs_vals(subsys_section, "KIND") 2147 CALL section_vals_get(kind_section, n_repetition=n_rep) 2148 2149 DO i = 1, SIZE(atomic_kind_set) 2150 atomic_kind => atomic_kind_set(i) 2151 CALL get_atomic_kind(atomic_kind=atomic_kind, & 2152 fist_potential=fist_potential, name=kind_name) 2153 CALL uppercase(kind_name) 2154 found = .FALSE. 2155 2156 ! Try to find a matching KIND section in the SUBSYS section and read the 2157 ! MM_RADIUS field if it is present. In case the kind section is never 2158 ! encountered, the mm_radius remains zero. 2159 mm_radius = 0.0_dp 2160 DO i_rep = 1, n_rep 2161 CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", & 2162 c_val=inp_kind_name, i_rep_section=i_rep) 2163 CALL uppercase(inp_kind_name) 2164 IF (iw > 0) WRITE (iw, *) "Matching kinds for MM_RADIUS :: '", & 2165 TRIM(kind_name), "' with '", TRIM(inp_kind_name), "'" 2166 IF (TRIM(kind_name) == TRIM(inp_kind_name)) THEN 2167 CALL section_vals_val_get(kind_section, i_rep_section=i_rep, & 2168 keyword_name="MM_RADIUS", r_val=mm_radius) 2169 CALL issue_duplications(found, "MM_RADIUS", kind_name) 2170 found = .TRUE. 2171 END IF 2172 END DO 2173 2174 CALL set_potential(potential=fist_potential, mm_radius=mm_radius) 2175 END DO 2176 CALL timestop(handle) 2177 END SUBROUTINE force_field_pack_radius 2178 2179! ************************************************************************************************** 2180!> \brief Set up the polarizable FF parameters 2181!> \param atomic_kind_set ... 2182!> \param iw ... 2183!> \param inp_info ... 2184!> \author Toon.Verstraelen@gmail.com 2185! ************************************************************************************************** 2186 SUBROUTINE force_field_pack_pol(atomic_kind_set, iw, inp_info) 2187 2188 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 2189 INTEGER, INTENT(IN) :: iw 2190 TYPE(input_info_type), POINTER :: inp_info 2191 2192 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_pol', & 2193 routineP = moduleN//':'//routineN 2194 2195 CHARACTER(LEN=default_string_length) :: kind_name 2196 INTEGER :: handle, i, j 2197 LOGICAL :: found 2198 REAL(KIND=dp) :: apol, cpol 2199 TYPE(atomic_kind_type), POINTER :: atomic_kind 2200 TYPE(fist_potential_type), POINTER :: fist_potential 2201 2202 CALL timeset(routineN, handle) 2203 2204 DO i = 1, SIZE(atomic_kind_set) 2205 atomic_kind => atomic_kind_set(i) 2206 CALL get_atomic_kind(atomic_kind=atomic_kind, & 2207 fist_potential=fist_potential, & 2208 name=kind_name) 2209 CALL get_potential(potential=fist_potential, apol=apol, cpol=cpol) 2210 CALL uppercase(kind_name) 2211 found = .FALSE. 2212 2213 ! Always have the input param last to overwrite all the other ones 2214 IF (ASSOCIATED(inp_info%apol_atm)) THEN 2215 DO j = 1, SIZE(inp_info%apol_atm) 2216 IF (iw > 0) WRITE (iw, *) "Matching kinds for APOL :: '", & 2217 TRIM(kind_name), "' with '", TRIM(inp_info%apol_atm(j)), "'" 2218 IF ((inp_info%apol_atm(j)) == kind_name) THEN 2219 apol = inp_info%apol(j) 2220 CALL issue_duplications(found, "APOL", kind_name) 2221 found = .TRUE. 2222 END IF 2223 END DO 2224 END IF 2225 2226 IF (ASSOCIATED(inp_info%cpol_atm)) THEN 2227 DO j = 1, SIZE(inp_info%cpol_atm) 2228 IF (iw > 0) WRITE (iw, *) "Matching kinds for CPOL :: '", & 2229 TRIM(kind_name), "' with '", TRIM(inp_info%cpol_atm(j)), "'" 2230 IF ((inp_info%cpol_atm(j)) == kind_name) THEN 2231 cpol = inp_info%cpol(j) 2232 CALL issue_duplications(found, "CPOL", kind_name) 2233 found = .TRUE. 2234 END IF 2235 END DO 2236 END IF 2237 2238 CALL set_potential(potential=fist_potential, apol=apol, cpol=cpol) 2239 END DO 2240 CALL timestop(handle) 2241 END SUBROUTINE force_field_pack_pol 2242 2243! ************************************************************************************************** 2244!> \brief Set up damping parameters 2245!> \param atomic_kind_set ... 2246!> \param iw ... 2247!> \param inp_info ... 2248! ************************************************************************************************** 2249 SUBROUTINE force_field_pack_damp(atomic_kind_set, & 2250 iw, inp_info) 2251 2252 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 2253 INTEGER :: iw 2254 TYPE(input_info_type), POINTER :: inp_info 2255 2256 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_damp', & 2257 routineP = moduleN//':'//routineN 2258 2259 CHARACTER(len=default_string_length) :: atm_name1, atm_name2, my_atm_name1, & 2260 my_atm_name2 2261 INTEGER :: handle2, i, j, k, nkinds 2262 LOGICAL :: found 2263 TYPE(atomic_kind_type), POINTER :: atomic_kind, atomic_kind2 2264 TYPE(damping_p_type), POINTER :: damping 2265 2266 CALL timeset(routineN, handle2) 2267 NULLIFY (damping) 2268 nkinds = SIZE(atomic_kind_set) 2269 2270 DO j = 1, SIZE(atomic_kind_set) 2271 atomic_kind => atomic_kind_set(j) 2272 CALL get_atomic_kind(atomic_kind=atomic_kind, & 2273 name=atm_name1) 2274 CALL UPPERCASE(atm_name1) 2275 IF (ASSOCIATED(inp_info%damping_list)) THEN 2276 DO i = 1, SIZE(inp_info%damping_list) 2277 my_atm_name1 = inp_info%damping_list(i)%atm_name1 2278 my_atm_name2 = inp_info%damping_list(i)%atm_name2 2279 2280 IF (iw > 0) WRITE (iw, *) "Damping Checking ::", TRIM(my_atm_name1), & 2281 TRIM(atm_name1) 2282 IF (my_atm_name1 == atm_name1) THEN 2283 IF (.NOT. ASSOCIATED(damping)) THEN 2284 CALL damping_p_create(damping, nkinds) 2285 END IF 2286 2287 found = .FALSE. 2288 DO k = 1, SIZE(atomic_kind_set) 2289 atomic_kind2 => atomic_kind_set(k) 2290 CALL get_atomic_kind(atomic_kind=atomic_kind2, & 2291 name=atm_name2) 2292 CALL UPPERCASE(atm_name2) 2293 2294 IF (my_atm_name2 == atm_name2) THEN 2295 IF (damping%damp(k)%bij /= HUGE(0.0_dp)) found = .TRUE. 2296 CALL issue_duplications(found, "Damping", atm_name1) 2297 found = .TRUE. 2298 2299 SELECT CASE (TRIM(inp_info%damping_list(i)%dtype)) 2300 CASE ('TANG-TOENNIES') 2301 damping%damp(k)%itype = tang_toennies 2302 CASE DEFAULT 2303 CPABORT("Unknown damping type.") 2304 END SELECT 2305 damping%damp(k)%order = inp_info%damping_list(i)%order 2306 damping%damp(k)%bij = inp_info%damping_list(i)%bij 2307 damping%damp(k)%cij = inp_info%damping_list(i)%cij 2308 ENDIF 2309 2310 END DO 2311 IF (.NOT. found) & 2312 CALL cp_warn(__LOCATION__, & 2313 "Atom "//TRIM(my_atm_name2)// & 2314 " in damping parameters for atom "//TRIM(my_atm_name1)// & 2315 " not found! ") 2316 ENDIF 2317 END DO 2318 2319 ENDIF 2320 2321 CALL set_atomic_kind(atomic_kind=atomic_kind, damping=damping) 2322 NULLIFY (damping) 2323 END DO 2324 2325 CALL timestop(handle2) 2326 2327 END SUBROUTINE force_field_pack_damp 2328 2329! ************************************************************************************************** 2330!> \brief Set up shell potential parameters 2331!> \param particle_set ... 2332!> \param atomic_kind_set ... 2333!> \param molecule_kind_set ... 2334!> \param molecule_set ... 2335!> \param root_section ... 2336!> \param subsys_section ... 2337!> \param shell_particle_set ... 2338!> \param core_particle_set ... 2339!> \param cell ... 2340!> \param iw ... 2341!> \param inp_info ... 2342! ************************************************************************************************** 2343 SUBROUTINE force_field_pack_shell(particle_set, atomic_kind_set, & 2344 molecule_kind_set, molecule_set, root_section, subsys_section, & 2345 shell_particle_set, core_particle_set, cell, iw, inp_info) 2346 2347 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 2348 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 2349 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 2350 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 2351 TYPE(section_vals_type), POINTER :: root_section, subsys_section 2352 TYPE(particle_type), DIMENSION(:), POINTER :: shell_particle_set, core_particle_set 2353 TYPE(cell_type), POINTER :: cell 2354 INTEGER :: iw 2355 TYPE(input_info_type), POINTER :: inp_info 2356 2357 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_shell', & 2358 routineP = moduleN//':'//routineN 2359 2360 CHARACTER(LEN=default_string_length) :: atmname 2361 INTEGER :: counter, first, first_shell, handle2, i, & 2362 j, last, last_shell, n, natom, nmol, & 2363 nshell_tot 2364 INTEGER, DIMENSION(:), POINTER :: molecule_list, shell_list_tmp 2365 LOGICAL :: core_coord_read, found_shell, is_a_shell, is_link_atom, null_massfrac, only_qm, & 2366 save_mem, shell_adiabatic, shell_coord_read 2367 REAL(KIND=dp) :: atmmass 2368 TYPE(atomic_kind_type), POINTER :: atomic_kind 2369 TYPE(molecule_kind_type), POINTER :: molecule_kind 2370 TYPE(molecule_type), POINTER :: molecule 2371 TYPE(section_vals_type), POINTER :: global_section 2372 TYPE(shell_kind_type), POINTER :: shell 2373 TYPE(shell_type), DIMENSION(:), POINTER :: shell_list 2374 2375 CALL timeset(routineN, handle2) 2376 nshell_tot = 0 2377 n = 0 2378 first_shell = 1 2379 null_massfrac = .FALSE. 2380 core_coord_read = .FALSE. 2381 shell_coord_read = .FALSE. 2382 2383 NULLIFY (global_section) 2384 global_section => section_vals_get_subs_vals(root_section, "GLOBAL") 2385 CALL section_vals_val_get(global_section, "SAVE_MEM", l_val=save_mem) 2386 2387 DO i = 1, SIZE(atomic_kind_set) 2388 atomic_kind => atomic_kind_set(i) 2389 CALL get_atomic_kind(atomic_kind=atomic_kind, & 2390 name=atmname) 2391 2392 found_shell = .FALSE. 2393 only_qm = qmmm_ff_precond_only_qm(id1=atmname, is_link=is_link_atom) 2394 CALL uppercase(atmname) 2395 2396 ! The shell potential can be defined only from input 2397 IF (ASSOCIATED(inp_info%shell_list)) THEN 2398 DO j = 1, SIZE(inp_info%shell_list) 2399 IF (iw > 0) WRITE (iw, *) "Shell Checking ::", TRIM(inp_info%shell_list(j)%atm_name), atmname 2400 IF ((inp_info%shell_list(j)%atm_name) == atmname) THEN 2401 CALL get_atomic_kind(atomic_kind=atomic_kind, & 2402 shell=shell, mass=atmmass, natom=natom) 2403 IF (.NOT. ASSOCIATED(shell)) THEN 2404 CALL shell_create(shell) 2405 END IF 2406 nshell_tot = nshell_tot + natom 2407 shell%charge_core = inp_info%shell_list(j)%shell%charge_core 2408 shell%charge_shell = inp_info%shell_list(j)%shell%charge_shell 2409 shell%massfrac = inp_info%shell_list(j)%shell%massfrac 2410 IF (shell%massfrac < EPSILON(1.0_dp)) null_massfrac = .TRUE. 2411 shell%k2_spring = inp_info%shell_list(j)%shell%k2_spring 2412 shell%k4_spring = inp_info%shell_list(j)%shell%k4_spring 2413 shell%max_dist = inp_info%shell_list(j)%shell%max_dist 2414 shell%shell_cutoff = inp_info%shell_list(j)%shell%shell_cutoff 2415 shell%mass_shell = shell%massfrac*atmmass 2416 shell%mass_core = atmmass - shell%mass_shell 2417 CALL issue_duplications(found_shell, "Shell", atmname) 2418 found_shell = .TRUE. 2419 CALL set_atomic_kind(atomic_kind=atomic_kind, & 2420 shell=shell, shell_active=.TRUE.) 2421 CALL shell_release(shell) 2422 END IF 2423 END DO ! j shell kind 2424 END IF ! associated shell_list 2425 END DO ! i atomic kind 2426 2427 IF (iw > 0) WRITE (iw, *) "Total number of particles with a shell :: ", nshell_tot 2428 ! If shell-model is present: Create particle_set of shells (coord. vel. force) 2429 NULLIFY (shell_particle_set) 2430 NULLIFY (core_particle_set) 2431 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, shell_adiabatic=shell_adiabatic) 2432 IF (nshell_tot > 0) THEN 2433 IF (shell_adiabatic .AND. null_massfrac) THEN 2434 CPABORT("Shell-model adiabatic: at least one shell_kind has mass zero") 2435 END IF 2436 CALL allocate_particle_set(shell_particle_set, nshell_tot) 2437 CALL allocate_particle_set(core_particle_set, nshell_tot) 2438 counter = 0 2439 ! Initialise the shell (and core) coordinates with the particle (atomic) coordinates, 2440 ! count the shell and set pointers 2441 DO i = 1, SIZE(particle_set) 2442 NULLIFY (atomic_kind) 2443 NULLIFY (shell) 2444 atomic_kind => particle_set(i)%atomic_kind 2445 CALL get_atomic_kind(atomic_kind=atomic_kind, shell_active=is_a_shell) 2446 IF (is_a_shell) THEN 2447 counter = counter + 1 2448 particle_set(i)%shell_index = counter 2449 shell_particle_set(counter)%shell_index = counter 2450 shell_particle_set(counter)%atomic_kind => particle_set(i)%atomic_kind 2451 shell_particle_set(counter)%r(1:3) = particle_set(i)%r(1:3) 2452 shell_particle_set(counter)%atom_index = i 2453 core_particle_set(counter)%shell_index = counter 2454 core_particle_set(counter)%atomic_kind => particle_set(i)%atomic_kind 2455 core_particle_set(counter)%r(1:3) = particle_set(i)%r(1:3) 2456 core_particle_set(counter)%atom_index = i 2457 ELSE 2458 particle_set(i)%shell_index = 0 2459 END IF 2460 END DO 2461 CPASSERT(counter == nshell_tot) 2462 END IF 2463 2464 ! Read the shell (and core) coordinates from the restart file, if available 2465 CALL read_binary_cs_coordinates("SHELL", shell_particle_set, root_section, & 2466 subsys_section, shell_coord_read) 2467 CALL read_binary_cs_coordinates("CORE", core_particle_set, root_section, & 2468 subsys_section, core_coord_read) 2469 2470 IF (nshell_tot > 0) THEN 2471 ! Read the shell (and core) coordinates from the input, if no coordinates were found 2472 ! in the restart file 2473 IF (shell_adiabatic) THEN 2474 IF (.NOT. (core_coord_read .AND. shell_coord_read)) THEN 2475 CALL read_shell_coord_input(particle_set, shell_particle_set, cell, & 2476 subsys_section, core_particle_set, & 2477 save_mem=save_mem) 2478 END IF 2479 ELSE 2480 IF (.NOT. shell_coord_read) THEN 2481 CALL read_shell_coord_input(particle_set, shell_particle_set, cell, & 2482 subsys_section, save_mem=save_mem) 2483 END IF 2484 END IF 2485 ! Determine the number of shells per molecule kind 2486 n = 0 2487 DO i = 1, SIZE(molecule_kind_set) 2488 molecule_kind => molecule_kind_set(i) 2489 CALL get_molecule_kind(molecule_kind=molecule_kind, molecule_list=molecule_list, & 2490 natom=natom, nmolecule=nmol) 2491 molecule => molecule_set(molecule_list(1)) 2492 CALL get_molecule(molecule=molecule, first_atom=first, last_atom=last) 2493 ALLOCATE (shell_list_tmp(natom)) 2494 counter = 0 2495 DO j = first, last 2496 atomic_kind => particle_set(j)%atomic_kind 2497 CALL get_atomic_kind(atomic_kind=atomic_kind, shell_active=is_a_shell) 2498 IF (is_a_shell) THEN 2499 counter = counter + 1 2500 shell_list_tmp(counter) = j - first + 1 2501 first_shell = MIN(first_shell, MAX(1, particle_set(j)%shell_index)) 2502 END IF 2503 END DO ! j atom in molecule_kind i, molecule 1 of the molecule_list 2504 IF (counter /= 0) THEN 2505 ! Setup of fist_shell and last_shell for all molecules.. 2506 DO j = 1, SIZE(molecule_list) 2507 last_shell = first_shell + counter - 1 2508 molecule => molecule_set(molecule_list(j)) 2509 molecule%first_shell = first_shell 2510 molecule%last_shell = last_shell 2511 first_shell = last_shell + 1 2512 END DO 2513 ! Setup of shell_list 2514 CALL get_molecule_kind(molecule_kind=molecule_kind, shell_list=shell_list) 2515 IF (ASSOCIATED(shell_list)) THEN 2516 DEALLOCATE (shell_list) 2517 END IF 2518 ALLOCATE (shell_list(counter)) 2519 DO j = 1, counter 2520 shell_list(j)%a = shell_list_tmp(j) 2521 atomic_kind => particle_set(shell_list_tmp(j) + first - 1)%atomic_kind 2522 CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname, shell=shell) 2523 CALL uppercase(atmname) 2524 shell_list(j)%name = atmname 2525 shell_list(j)%shell_kind => shell 2526 END DO 2527 CALL set_molecule_kind(molecule_kind=molecule_kind, nshell=counter, shell_list=shell_list) 2528 END IF 2529 DEALLOCATE (shell_list_tmp) 2530 n = n + nmol*counter 2531 END DO ! i molecule kind 2532 END IF 2533 CPASSERT(first_shell - 1 == nshell_tot) 2534 CPASSERT(n == nshell_tot) 2535 CALL timestop(handle2) 2536 2537 END SUBROUTINE force_field_pack_shell 2538 2539! ************************************************************************************************** 2540!> \brief Assign input and potential info to potparm_nonbond14 2541!> \param atomic_kind_set ... 2542!> \param ff_type ... 2543!> \param qmmm_env ... 2544!> \param iw ... 2545!> \param Ainfo ... 2546!> \param chm_info ... 2547!> \param inp_info ... 2548!> \param gro_info ... 2549!> \param amb_info ... 2550!> \param potparm_nonbond14 ... 2551!> \param ewald_env ... 2552! ************************************************************************************************** 2553 SUBROUTINE force_field_pack_nonbond14(atomic_kind_set, ff_type, qmmm_env, iw, & 2554 Ainfo, chm_info, inp_info, gro_info, amb_info, potparm_nonbond14, ewald_env) 2555 2556 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 2557 TYPE(force_field_type), INTENT(INOUT) :: ff_type 2558 TYPE(qmmm_env_mm_type), POINTER :: qmmm_env 2559 INTEGER :: iw 2560 CHARACTER(LEN=default_string_length), & 2561 DIMENSION(:), POINTER :: Ainfo 2562 TYPE(charmm_info_type), POINTER :: chm_info 2563 TYPE(input_info_type), POINTER :: inp_info 2564 TYPE(gromos_info_type), POINTER :: gro_info 2565 TYPE(amber_info_type), POINTER :: amb_info 2566 TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond14 2567 TYPE(ewald_environment_type), POINTER :: ewald_env 2568 2569 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_nonbond14', & 2570 routineP = moduleN//':'//routineN 2571 2572 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a_local, & 2573 name_atm_b, name_atm_b_local 2574 INTEGER :: handle2, i, ii, j, jj, k, match_names 2575 LOGICAL :: found, found_a, found_b, only_qm, & 2576 use_qmmm_ff 2577 REAL(KIND=dp) :: epsilon0, epsilon_a, epsilon_b, & 2578 ewald_rcut, rmin, rmin2_a, rmin2_b 2579 TYPE(atomic_kind_type), POINTER :: atomic_kind 2580 TYPE(pair_potential_single_type), POINTER :: pot 2581 2582 use_qmmm_ff = qmmm_env%use_qmmm_ff 2583 NULLIFY (pot) 2584 CALL ewald_env_get(ewald_env, rcut=ewald_rcut) 2585 CALL timeset(routineN, handle2) 2586 CALL pair_potential_pp_create(potparm_nonbond14, SIZE(atomic_kind_set)) 2587 DO i = 1, SIZE(atomic_kind_set) 2588 atomic_kind => atomic_kind_set(i) 2589 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local) 2590 DO j = i, SIZE(atomic_kind_set) 2591 atomic_kind => atomic_kind_set(j) 2592 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local) 2593 found = .FALSE. 2594 found_a = .FALSE. 2595 found_b = .FALSE. 2596 name_atm_a = name_atm_a_local 2597 name_atm_b = name_atm_b_local 2598 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b) 2599 CALL uppercase(name_atm_a) 2600 CALL uppercase(name_atm_b) 2601 pot => potparm_nonbond14%pot(i, j)%pot 2602 2603 ! loop over params from GROMOS 2604 IF (ASSOCIATED(gro_info%nonbond_a_14)) THEN 2605 ii = 0 2606 jj = 0 2607 DO k = 1, SIZE(gro_info%nonbond_a_14) 2608 IF (TRIM(name_atm_a) == TRIM(gro_info%nonbond_a_14(k))) THEN 2609 ii = k 2610 found_a = .TRUE. 2611 EXIT 2612 END IF 2613 END DO 2614 DO k = 1, SIZE(gro_info%nonbond_a_14) 2615 IF (TRIM(name_atm_b) == TRIM(gro_info%nonbond_a_14(k))) THEN 2616 jj = k 2617 found_b = .TRUE. 2618 EXIT 2619 END IF 2620 END DO 2621 IF (ii /= 0 .AND. jj /= 0) THEN 2622 CALL pair_potential_lj_create(pot%set(1)%lj) 2623 pot%type = lj_type 2624 pot%at1 = name_atm_a 2625 pot%at2 = name_atm_b 2626 pot%set(1)%lj%epsilon = 1.0_dp 2627 pot%set(1)%lj%sigma6 = gro_info%nonbond_c6_14(ii, jj) 2628 pot%set(1)%lj%sigma12 = gro_info%nonbond_c12_14(ii, jj) 2629 pot%rcutsq = (10.0_dp*bohr)**2 2630 CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b) 2631 found = .TRUE. 2632 END IF 2633 END IF 2634 2635 ! loop over params from CHARMM 2636 ii = 0 2637 jj = 0 2638 IF (ASSOCIATED(chm_info%nonbond_a_14)) THEN 2639 DO k = 1, SIZE(chm_info%nonbond_a_14) 2640 IF ((name_atm_a) == (chm_info%nonbond_a_14(k))) THEN 2641 ii = k 2642 rmin2_a = chm_info%nonbond_rmin2_14(k) 2643 epsilon_a = chm_info%nonbond_eps_14(k) 2644 found_a = .TRUE. 2645 END IF 2646 END DO 2647 DO k = 1, SIZE(chm_info%nonbond_a_14) 2648 IF ((name_atm_b) == (chm_info%nonbond_a_14(k))) THEN 2649 jj = k 2650 rmin2_b = chm_info%nonbond_rmin2_14(k) 2651 epsilon_b = chm_info%nonbond_eps_14(k) 2652 found_b = .TRUE. 2653 END IF 2654 END DO 2655 END IF 2656 IF (ASSOCIATED(chm_info%nonbond_a)) THEN 2657 IF (.NOT. found_a) THEN 2658 DO k = 1, SIZE(chm_info%nonbond_a) 2659 IF ((name_atm_a) == (chm_info%nonbond_a(k))) THEN 2660 ii = k 2661 rmin2_a = chm_info%nonbond_rmin2(k) 2662 epsilon_a = chm_info%nonbond_eps(k) 2663 END IF 2664 END DO 2665 END IF 2666 IF (.NOT. found_b) THEN 2667 DO k = 1, SIZE(chm_info%nonbond_a) 2668 IF ((name_atm_b) == (chm_info%nonbond_a(k))) THEN 2669 jj = k 2670 rmin2_b = chm_info%nonbond_rmin2(k) 2671 epsilon_b = chm_info%nonbond_eps(k) 2672 END IF 2673 END DO 2674 END IF 2675 END IF 2676 IF (ii /= 0 .AND. jj /= 0) THEN 2677 rmin = rmin2_a + rmin2_b 2678 ! ABS to allow for mixing the two different sign conventions for epsilon 2679 epsilon0 = SQRT(ABS(epsilon_a*epsilon_b)) 2680 CALL pair_potential_lj_create(pot%set(1)%lj) 2681 pot%type = lj_charmm_type 2682 pot%at1 = name_atm_a 2683 pot%at2 = name_atm_b 2684 pot%set(1)%lj%epsilon = epsilon0 2685 pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6 2686 pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12 2687 pot%rcutsq = (10.0_dp*bohr)**2 2688 CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b) 2689 found = .TRUE. 2690 END IF 2691 2692 ! loop over params from AMBER 2693 IF (ASSOCIATED(amb_info%nonbond_a)) THEN 2694 ii = 0 2695 jj = 0 2696 IF (.NOT. found_a) THEN 2697 DO k = 1, SIZE(amb_info%nonbond_a) 2698 IF ((name_atm_a) == (amb_info%nonbond_a(k))) THEN 2699 ii = k 2700 rmin2_a = amb_info%nonbond_rmin2(k) 2701 epsilon_a = amb_info%nonbond_eps(k) 2702 END IF 2703 END DO 2704 END IF 2705 IF (.NOT. found_b) THEN 2706 DO k = 1, SIZE(amb_info%nonbond_a) 2707 IF ((name_atm_b) == (amb_info%nonbond_a(k))) THEN 2708 jj = k 2709 rmin2_b = amb_info%nonbond_rmin2(k) 2710 epsilon_b = amb_info%nonbond_eps(k) 2711 END IF 2712 END DO 2713 END IF 2714 IF (ii /= 0 .AND. jj /= 0) THEN 2715 rmin = rmin2_a + rmin2_b 2716 ! ABS to allow for mixing the two different sign conventions for epsilon 2717 epsilon0 = SQRT(ABS(epsilon_a*epsilon_b)) 2718 CALL pair_potential_lj_create(pot%set(1)%lj) 2719 pot%type = lj_charmm_type 2720 pot%at1 = name_atm_a 2721 pot%at2 = name_atm_b 2722 pot%set(1)%lj%epsilon = epsilon0 2723 pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6 2724 pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12 2725 pot%rcutsq = (10.0_dp*bohr)**2 2726 CALL issue_duplications(found, "Lennard-Jones", name_atm_a, & 2727 name_atm_b) 2728 found = .TRUE. 2729 END IF 2730 END IF 2731 2732 ! always have the input param last to overwrite all the other ones 2733 IF (ASSOCIATED(inp_info%nonbonded14)) THEN 2734 DO k = 1, SIZE(inp_info%nonbonded14%pot) 2735 IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & 2736 " with ", TRIM(inp_info%nonbonded14%pot(k)%pot%at1), & 2737 TRIM(inp_info%nonbonded14%pot(k)%pot%at2) 2738 IF ((((name_atm_a) == (inp_info%nonbonded14%pot(k)%pot%at1)) .AND. & 2739 ((name_atm_b) == (inp_info%nonbonded14%pot(k)%pot%at2))) .OR. & 2740 (((name_atm_b) == (inp_info%nonbonded14%pot(k)%pot%at1)) .AND. & 2741 ((name_atm_a) == (inp_info%nonbonded14%pot(k)%pot%at2)))) THEN 2742 IF (ff_type%multiple_potential) THEN 2743 CALL pair_potential_single_add(inp_info%nonbonded14%pot(k)%pot, pot) 2744 IF (found) & 2745 CALL cp_warn(__LOCATION__, & 2746 "Multiple ONFO declaration: "//TRIM(name_atm_a)// & 2747 " and "//TRIM(name_atm_b)//" ADDING! ") 2748 potparm_nonbond14%pot(i, j)%pot => pot 2749 potparm_nonbond14%pot(j, i)%pot => pot 2750 ELSE 2751 CALL pair_potential_single_copy(inp_info%nonbonded14%pot(k)%pot, pot) 2752 IF (found) & 2753 CALL cp_warn(__LOCATION__, & 2754 "Multiple ONFO declarations: "//TRIM(name_atm_a)// & 2755 " and "//TRIM(name_atm_b)//" OVERWRITING! ") 2756 END IF 2757 IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b) 2758 found = .TRUE. 2759 END IF 2760 END DO 2761 END IF 2762 ! at the very end we offer the possibility to overwrite the parameters for QM/MM 2763 ! nonbonded interactions 2764 IF (use_qmmm_ff) THEN 2765 match_names = 0 2766 IF ((name_atm_a) == (name_atm_a_local)) match_names = match_names + 1 2767 IF ((name_atm_b) == (name_atm_b_local)) match_names = match_names + 1 2768 IF (match_names == 1) THEN 2769 IF (ASSOCIATED(qmmm_env%inp_info%nonbonded14)) THEN 2770 DO k = 1, SIZE(qmmm_env%inp_info%nonbonded14%pot) 2771 IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & 2772 " with ", TRIM(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1), & 2773 TRIM(qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2) 2774 IF ((((name_atm_a) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)) .AND. & 2775 ((name_atm_b) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2))) .OR. & 2776 (((name_atm_b) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at1)) .AND. & 2777 ((name_atm_a) == (qmmm_env%inp_info%nonbonded14%pot(k)%pot%at2)))) THEN 2778 IF (qmmm_env%multiple_potential) THEN 2779 CALL pair_potential_single_add(qmmm_env%inp_info%nonbonded14%pot(k)%pot, pot) 2780 IF (found) & 2781 CALL cp_warn(__LOCATION__, & 2782 "Multiple ONFO declaration: "//TRIM(name_atm_a)// & 2783 " and "//TRIM(name_atm_b)//" ADDING QM/MM forcefield specifications! ") 2784 potparm_nonbond14%pot(i, j)%pot => pot 2785 potparm_nonbond14%pot(j, i)%pot => pot 2786 ELSE 2787 CALL pair_potential_single_copy(qmmm_env%inp_info%nonbonded14%pot(k)%pot, pot) 2788 IF (found) & 2789 CALL cp_warn(__LOCATION__, & 2790 "Multiple ONFO declaration: "//TRIM(name_atm_a)// & 2791 " and "//TRIM(name_atm_b)//" OVERWRITING QM/MM forcefield specifications! ") 2792 END IF 2793 IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), & 2794 " ", TRIM(name_atm_b) 2795 found = .TRUE. 2796 END IF 2797 END DO 2798 END IF 2799 END IF 2800 END IF 2801 2802 IF (.NOT. found) THEN 2803 CALL store_FF_missing_par(atm1=TRIM(name_atm_a), & 2804 atm2=TRIM(name_atm_b), & 2805 type_name="Spline_Bond_Env", & 2806 array=Ainfo) 2807 CALL pair_potential_single_clean(pot) 2808 pot%type = nn_type 2809 pot%at1 = name_atm_a 2810 pot%at2 = name_atm_b 2811 END IF 2812 ! If defined global RCUT let's use it 2813 IF (ff_type%rcut_nb > 0.0_dp) THEN 2814 pot%rcutsq = ff_type%rcut_nb*ff_type%rcut_nb 2815 END IF 2816 ! Cutoff is defined always as the maximum between the FF and Ewald 2817 pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut) 2818 IF (only_qm) THEN 2819 CALL pair_potential_single_clean(pot) 2820 END IF 2821 END DO ! atom kind j 2822 END DO ! atom kind i 2823 CALL timestop(handle2) 2824 2825 END SUBROUTINE force_field_pack_nonbond14 2826 2827! ************************************************************************************************** 2828!> \brief Assign input and potential info to potparm_nonbond 2829!> \param atomic_kind_set ... 2830!> \param ff_type ... 2831!> \param qmmm_env ... 2832!> \param fatal ... 2833!> \param iw ... 2834!> \param Ainfo ... 2835!> \param chm_info ... 2836!> \param inp_info ... 2837!> \param gro_info ... 2838!> \param amb_info ... 2839!> \param potparm_nonbond ... 2840!> \param ewald_env ... 2841! ************************************************************************************************** 2842 SUBROUTINE force_field_pack_nonbond(atomic_kind_set, ff_type, qmmm_env, fatal, & 2843 iw, Ainfo, chm_info, inp_info, gro_info, amb_info, potparm_nonbond, & 2844 ewald_env) 2845 2846 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 2847 TYPE(force_field_type), INTENT(INOUT) :: ff_type 2848 TYPE(qmmm_env_mm_type), POINTER :: qmmm_env 2849 LOGICAL :: fatal 2850 INTEGER :: iw 2851 CHARACTER(LEN=default_string_length), & 2852 DIMENSION(:), POINTER :: Ainfo 2853 TYPE(charmm_info_type), POINTER :: chm_info 2854 TYPE(input_info_type), POINTER :: inp_info 2855 TYPE(gromos_info_type), POINTER :: gro_info 2856 TYPE(amber_info_type), POINTER :: amb_info 2857 TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond 2858 TYPE(ewald_environment_type), POINTER :: ewald_env 2859 2860 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_nonbond', & 2861 routineP = moduleN//':'//routineN 2862 2863 CHARACTER(LEN=default_string_length) :: name_atm_a, name_atm_a_local, & 2864 name_atm_b, name_atm_b_local 2865 INTEGER :: handle2, i, ii, j, jj, k, match_names 2866 LOGICAL :: found, is_a_shell, is_b_shell, only_qm, & 2867 use_qmmm_ff 2868 REAL(KIND=dp) :: epsilon0, ewald_rcut, rmin 2869 TYPE(atomic_kind_type), POINTER :: atomic_kind 2870 TYPE(pair_potential_single_type), POINTER :: pot 2871 2872 CALL timeset(routineN, handle2) 2873 use_qmmm_ff = qmmm_env%use_qmmm_ff 2874 NULLIFY (pot) 2875 CALL ewald_env_get(ewald_env, rcut=ewald_rcut) 2876 CALL pair_potential_pp_create(potparm_nonbond, SIZE(atomic_kind_set)) 2877 DO i = 1, SIZE(atomic_kind_set) 2878 atomic_kind => atomic_kind_set(i) 2879 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_a_local, & 2880 shell_active=is_a_shell) 2881 DO j = i, SIZE(atomic_kind_set) 2882 atomic_kind => atomic_kind_set(j) 2883 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name_atm_b_local, & 2884 shell_active=is_b_shell) 2885 found = .FALSE. 2886 name_atm_a = name_atm_a_local 2887 name_atm_b = name_atm_b_local 2888 only_qm = qmmm_ff_precond_only_qm(id1=name_atm_a, id2=name_atm_b) 2889 CALL uppercase(name_atm_a) 2890 CALL uppercase(name_atm_b) 2891 pot => potparm_nonbond%pot(i, j)%pot 2892 2893 ! loop over params from GROMOS 2894 IF (ASSOCIATED(gro_info%nonbond_a)) THEN 2895 ii = 0 2896 jj = 0 2897 DO k = 1, SIZE(gro_info%nonbond_a) 2898 IF (TRIM(name_atm_a) == TRIM(gro_info%nonbond_a(k))) THEN 2899 ii = k 2900 EXIT 2901 END IF 2902 END DO 2903 DO k = 1, SIZE(gro_info%nonbond_a) 2904 IF (TRIM(name_atm_b) == TRIM(gro_info%nonbond_a(k))) THEN 2905 jj = k 2906 EXIT 2907 END IF 2908 END DO 2909 2910 IF (ii /= 0 .AND. jj /= 0) THEN 2911 CALL pair_potential_lj_create(pot%set(1)%lj) 2912 pot%type = lj_type 2913 pot%at1 = name_atm_a 2914 pot%at2 = name_atm_b 2915 pot%set(1)%lj%epsilon = 1.0_dp 2916 pot%set(1)%lj%sigma6 = gro_info%nonbond_c6(ii, jj) 2917 pot%set(1)%lj%sigma12 = gro_info%nonbond_c12(ii, jj) 2918 pot%rcutsq = (10.0_dp*bohr)**2 2919 CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b) 2920 found = .TRUE. 2921 END IF 2922 END IF 2923 2924 ! loop over params from CHARMM 2925 IF (ASSOCIATED(chm_info%nonbond_a)) THEN 2926 ii = 0 2927 jj = 0 2928 DO k = 1, SIZE(chm_info%nonbond_a) 2929 IF ((name_atm_a) == (chm_info%nonbond_a(k))) THEN 2930 ii = k 2931 END IF 2932 END DO 2933 DO k = 1, SIZE(chm_info%nonbond_a) 2934 IF ((name_atm_b) == (chm_info%nonbond_a(k))) THEN 2935 jj = k 2936 END IF 2937 END DO 2938 2939 IF (ii /= 0 .AND. jj /= 0) THEN 2940 rmin = chm_info%nonbond_rmin2(ii) + chm_info%nonbond_rmin2(jj) 2941 epsilon0 = SQRT(chm_info%nonbond_eps(ii)* & 2942 chm_info%nonbond_eps(jj)) 2943 CALL pair_potential_lj_create(pot%set(1)%lj) 2944 pot%type = lj_charmm_type 2945 pot%at1 = name_atm_a 2946 pot%at2 = name_atm_b 2947 pot%set(1)%lj%epsilon = epsilon0 2948 pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6 2949 pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12 2950 pot%rcutsq = (10.0_dp*bohr)**2 2951 CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b) 2952 found = .TRUE. 2953 END IF 2954 END IF 2955 2956 ! loop over params from AMBER 2957 IF (ASSOCIATED(amb_info%nonbond_a)) THEN 2958 ii = 0 2959 jj = 0 2960 DO k = 1, SIZE(amb_info%nonbond_a) 2961 IF ((name_atm_a) == (amb_info%nonbond_a(k))) THEN 2962 ii = k 2963 END IF 2964 END DO 2965 DO k = 1, SIZE(amb_info%nonbond_a) 2966 IF ((name_atm_b) == (amb_info%nonbond_a(k))) THEN 2967 jj = k 2968 END IF 2969 END DO 2970 2971 IF (ii /= 0 .AND. jj /= 0) THEN 2972 rmin = amb_info%nonbond_rmin2(ii) + amb_info%nonbond_rmin2(jj) 2973 epsilon0 = SQRT(amb_info%nonbond_eps(ii)*amb_info%nonbond_eps(jj)) 2974 CALL pair_potential_lj_create(pot%set(1)%lj) 2975 pot%type = lj_charmm_type 2976 pot%at1 = name_atm_a 2977 pot%at2 = name_atm_b 2978 pot%set(1)%lj%epsilon = epsilon0 2979 pot%set(1)%lj%sigma6 = 0.5_dp*rmin**6 2980 pot%set(1)%lj%sigma12 = 0.25_dp*rmin**12 2981 pot%rcutsq = (10.0_dp*bohr)**2 2982 CALL issue_duplications(found, "Lennard-Jones", name_atm_a, name_atm_b) 2983 found = .TRUE. 2984 END IF 2985 END IF 2986 2987 ! always have the input param last to overwrite all the other ones 2988 IF (ASSOCIATED(inp_info%nonbonded)) THEN 2989 DO k = 1, SIZE(inp_info%nonbonded%pot) 2990 IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1) == "*") .OR. & 2991 (TRIM(inp_info%nonbonded%pot(k)%pot%at2) == "*")) CYCLE 2992 2993 IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & 2994 " with ", TRIM(inp_info%nonbonded%pot(k)%pot%at1), & 2995 TRIM(inp_info%nonbonded%pot(k)%pot%at2) 2996 IF ((((name_atm_a) == (inp_info%nonbonded%pot(k)%pot%at1)) .AND. & 2997 ((name_atm_b) == (inp_info%nonbonded%pot(k)%pot%at2))) .OR. & 2998 (((name_atm_b) == (inp_info%nonbonded%pot(k)%pot%at1)) .AND. & 2999 ((name_atm_a) == (inp_info%nonbonded%pot(k)%pot%at2)))) THEN 3000 IF (ff_type%multiple_potential) THEN 3001 CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot, pot) 3002 IF (found) & 3003 CALL cp_warn(__LOCATION__, & 3004 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & 3005 " and "//TRIM(name_atm_b)//" ADDING! ") 3006 potparm_nonbond%pot(i, j)%pot => pot 3007 potparm_nonbond%pot(j, i)%pot => pot 3008 ELSE 3009 CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot, pot) 3010 IF (found) & 3011 CALL cp_warn(__LOCATION__, & 3012 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & 3013 " and "//TRIM(name_atm_b)//" OVERWRITING! ") 3014 END IF 3015 IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b) 3016 found = .TRUE. 3017 END IF 3018 END DO 3019 ! Check for wildcards for one of the two types (if not associated yet) 3020 IF (.NOT. found) THEN 3021 DO k = 1, SIZE(inp_info%nonbonded%pot) 3022 IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1) == "*") .EQV. & 3023 (TRIM(inp_info%nonbonded%pot(k)%pot%at2) == "*")) CYCLE 3024 3025 IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & 3026 " with ", TRIM(inp_info%nonbonded%pot(k)%pot%at1), & 3027 TRIM(inp_info%nonbonded%pot(k)%pot%at2) 3028 3029 IF ((name_atm_a == inp_info%nonbonded%pot(k)%pot%at1) .OR. & 3030 (name_atm_b == inp_info%nonbonded%pot(k)%pot%at2) .OR. & 3031 (name_atm_b == inp_info%nonbonded%pot(k)%pot%at1) .OR. & 3032 (name_atm_a == inp_info%nonbonded%pot(k)%pot%at2)) THEN 3033 IF (ff_type%multiple_potential) THEN 3034 CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot, pot) 3035 IF (found) & 3036 CALL cp_warn(__LOCATION__, & 3037 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & 3038 " and "//TRIM(name_atm_b)//" ADDING! ") 3039 potparm_nonbond%pot(i, j)%pot => pot 3040 potparm_nonbond%pot(j, i)%pot => pot 3041 ELSE 3042 CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot, pot) 3043 IF (found) & 3044 CALL cp_warn(__LOCATION__, & 3045 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & 3046 " and "//TRIM(name_atm_b)//" OVERWRITING! ") 3047 END IF 3048 IF (iw > 0) WRITE (iw, *) " FOUND (one WILDCARD)", TRIM(name_atm_a), " ", TRIM(name_atm_b) 3049 found = .TRUE. 3050 END IF 3051 END DO 3052 END IF 3053 ! Check for wildcards for both types (if not associated yet) 3054 IF (.NOT. found) THEN 3055 DO k = 1, SIZE(inp_info%nonbonded%pot) 3056 IF ((TRIM(inp_info%nonbonded%pot(k)%pot%at1) /= "*") .OR. & 3057 (TRIM(inp_info%nonbonded%pot(k)%pot%at2) /= "*")) CYCLE 3058 3059 IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & 3060 " with ", TRIM(inp_info%nonbonded%pot(k)%pot%at1), & 3061 TRIM(inp_info%nonbonded%pot(k)%pot%at2) 3062 3063 IF (ff_type%multiple_potential) THEN 3064 CALL pair_potential_single_add(inp_info%nonbonded%pot(k)%pot, pot) 3065 IF (found) & 3066 CALL cp_warn(__LOCATION__, & 3067 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & 3068 " and "//TRIM(name_atm_b)//" ADDING! ") 3069 potparm_nonbond%pot(i, j)%pot => pot 3070 potparm_nonbond%pot(j, i)%pot => pot 3071 ELSE 3072 CALL pair_potential_single_copy(inp_info%nonbonded%pot(k)%pot, pot) 3073 IF (found) & 3074 CALL cp_warn(__LOCATION__, & 3075 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & 3076 " and "//TRIM(name_atm_b)//" OVERWRITING! ") 3077 END IF 3078 IF (iw > 0) WRITE (iw, *) " FOUND (both WILDCARDS)", TRIM(name_atm_a), " ", TRIM(name_atm_b) 3079 found = .TRUE. 3080 END DO 3081 END IF 3082 END IF 3083 3084 ! at the very end we offer the possibility to overwrite the parameters for QM/MM 3085 ! nonbonded interactions 3086 IF (use_qmmm_ff) THEN 3087 match_names = 0 3088 IF ((name_atm_a) == (name_atm_a_local)) match_names = match_names + 1 3089 IF ((name_atm_b) == (name_atm_b_local)) match_names = match_names + 1 3090 IF (match_names == 1) THEN 3091 IF (ASSOCIATED(qmmm_env%inp_info%nonbonded)) THEN 3092 DO k = 1, SIZE(qmmm_env%inp_info%nonbonded%pot) 3093 IF (iw > 0) WRITE (iw, *) " TESTING ", TRIM(name_atm_a), TRIM(name_atm_b), & 3094 " with ", TRIM(qmmm_env%inp_info%nonbonded%pot(k)%pot%at1), & 3095 TRIM(qmmm_env%inp_info%nonbonded%pot(k)%pot%at2) 3096 IF ((((name_atm_a) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at1)) .AND. & 3097 ((name_atm_b) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at2))) .OR. & 3098 (((name_atm_b) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at1)) .AND. & 3099 ((name_atm_a) == (qmmm_env%inp_info%nonbonded%pot(k)%pot%at2)))) THEN 3100 IF (qmmm_env%multiple_potential) THEN 3101 CALL pair_potential_single_add(qmmm_env%inp_info%nonbonded%pot(k)%pot, pot) 3102 IF (found) & 3103 CALL cp_warn(__LOCATION__, & 3104 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & 3105 " and "//TRIM(name_atm_b)//" ADDING QM/MM forcefield specifications! ") 3106 potparm_nonbond%pot(i, j)%pot => pot 3107 potparm_nonbond%pot(j, i)%pot => pot 3108 ELSE 3109 CALL pair_potential_single_copy(qmmm_env%inp_info%nonbonded%pot(k)%pot, pot) 3110 IF (found) & 3111 CALL cp_warn(__LOCATION__, & 3112 "Multiple NONBONDED declaration: "//TRIM(name_atm_a)// & 3113 " and "//TRIM(name_atm_b)//" OVERWRITING QM/MM forcefield specifications! ") 3114 END IF 3115 IF (iw > 0) WRITE (iw, *) " FOUND ", TRIM(name_atm_a), " ", TRIM(name_atm_b) 3116 found = .TRUE. 3117 END IF 3118 END DO 3119 END IF 3120 END IF 3121 END IF 3122 IF (.NOT. found) THEN 3123 CALL store_FF_missing_par(atm1=TRIM(name_atm_a), & 3124 atm2=TRIM(name_atm_b), & 3125 type_name="Spline_Non_Bond_Env", & 3126 fatal=fatal, & 3127 array=Ainfo) 3128 END IF 3129 ! If defined global RCUT let's use it 3130 IF (ff_type%rcut_nb > 0.0_dp) THEN 3131 pot%rcutsq = ff_type%rcut_nb*ff_type%rcut_nb 3132 END IF 3133 ! Cutoff is defined always as the maximum between the FF and Ewald 3134 pot%rcutsq = MAX(pot%rcutsq, ewald_rcut*ewald_rcut) 3135 ! Set the shell type 3136 IF ((is_a_shell .AND. .NOT. is_b_shell) .OR. (is_b_shell .AND. .NOT. is_a_shell)) THEN 3137 pot%shell_type = nosh_sh 3138 ELSE IF (is_a_shell .AND. is_b_shell) THEN 3139 pot%shell_type = sh_sh 3140 ELSE 3141 pot%shell_type = nosh_nosh 3142 END IF 3143 IF (only_qm) THEN 3144 CALL pair_potential_single_clean(pot) 3145 END IF 3146 END DO ! jkind 3147 END DO ! ikind 3148 CALL timestop(handle2) 3149 END SUBROUTINE force_field_pack_nonbond 3150 3151! ************************************************************************************************** 3152!> \brief create the pair potential spline environment 3153!> \param atomic_kind_set ... 3154!> \param ff_type ... 3155!> \param iw2 ... 3156!> \param iw3 ... 3157!> \param iw4 ... 3158!> \param potparm ... 3159!> \param do_zbl ... 3160!> \param nonbonded_type ... 3161! ************************************************************************************************** 3162 SUBROUTINE force_field_pack_splines(atomic_kind_set, ff_type, iw2, iw3, iw4, & 3163 potparm, do_zbl, nonbonded_type) 3164 3165 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 3166 TYPE(force_field_type), INTENT(INOUT) :: ff_type 3167 INTEGER :: iw2, iw3, iw4 3168 TYPE(pair_potential_pp_type), POINTER :: potparm 3169 LOGICAL, INTENT(IN) :: do_zbl 3170 CHARACTER(LEN=*), INTENT(IN) :: nonbonded_type 3171 3172 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_splines', & 3173 routineP = moduleN//':'//routineN 3174 3175 INTEGER :: handle2, ikind, jkind, n 3176 TYPE(spline_data_p_type), DIMENSION(:), POINTER :: spl_p 3177 TYPE(spline_environment_type), POINTER :: spline_env 3178 3179 CALL timeset(routineN, handle2) 3180 ! Figure out which nonbonded interactions happen to be identical, and 3181 ! prepare storage for these, avoiding duplicates. 3182 NULLIFY (spline_env) 3183 CALL get_nonbond_storage(spline_env, potparm, atomic_kind_set, & 3184 do_zbl, shift_cutoff=ff_type%shift_cutoff) 3185 ! Effectively compute the spline data. 3186 CALL spline_nonbond_control(spline_env, potparm, & 3187 atomic_kind_set, eps_spline=ff_type%eps_spline, & 3188 max_energy=ff_type%max_energy, rlow_nb=ff_type%rlow_nb, & 3189 emax_spline=ff_type%emax_spline, npoints=ff_type%npoints, iw=iw2, iw2=iw3, iw3=iw4, & 3190 do_zbl=do_zbl, shift_cutoff=ff_type%shift_cutoff, & 3191 nonbonded_type=nonbonded_type) 3192 ! Let the pointers on potparm point to the splines generated in 3193 ! spline_nonbond_control. 3194 DO ikind = 1, SIZE(potparm%pot, 1) 3195 DO jkind = ikind, SIZE(potparm%pot, 2) 3196 n = spline_env%spltab(ikind, jkind) 3197 spl_p => spline_env%spl_pp(n)%spl_p 3198 CALL spline_data_p_retain(spl_p) 3199 CALL spline_data_p_release(potparm%pot(ikind, jkind)%pot%pair_spline_data) 3200 potparm%pot(ikind, jkind)%pot%pair_spline_data => spl_p 3201 END DO 3202 END DO 3203 CALL spline_env_release(spline_env) 3204 CALL timestop(handle2) 3205 3206 END SUBROUTINE force_field_pack_splines 3207 3208! ************************************************************************************************** 3209!> \brief Compute the electrostatic interaction cutoffs 3210!> \param atomic_kind_set ... 3211!> \param ff_type ... 3212!> \param potparm_nonbond ... 3213!> \param ewald_env ... 3214!> \author Toon.Verstraelen@gmail.com 3215! ************************************************************************************************** 3216 SUBROUTINE force_field_pack_eicut(atomic_kind_set, ff_type, & 3217 potparm_nonbond, ewald_env) 3218 3219 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 3220 TYPE(force_field_type), INTENT(IN) :: ff_type 3221 TYPE(pair_potential_pp_type), POINTER :: potparm_nonbond 3222 TYPE(ewald_environment_type), POINTER :: ewald_env 3223 3224 CHARACTER(len=*), PARAMETER :: routineN = 'force_field_pack_eicut', & 3225 routineP = moduleN//':'//routineN 3226 3227 INTEGER :: ewald_type, handle, i1, i2, nkinds 3228 REAL(KIND=dp) :: alpha, beta, mm_radius1, mm_radius2, & 3229 rcut2, rcut2_ewald, tmp 3230 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: interaction_cutoffs 3231 TYPE(atomic_kind_type), POINTER :: atomic_kind 3232 3233 CALL timeset(routineN, handle) 3234 3235 tmp = 0.0_dp 3236 nkinds = SIZE(atomic_kind_set) 3237 ! allocate the array with interaction cutoffs for the electrostatics, used 3238 ! to make the electrostatic interaction continuous at ewald_env%rcut 3239 ALLOCATE (interaction_cutoffs(3, nkinds, nkinds)) 3240 interaction_cutoffs = 0.0_dp 3241 3242 ! compute the interaction cutoff if SHIFT_CUTOFF is active 3243 IF (ff_type%shift_cutoff) THEN 3244 CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type, & 3245 rcut=rcut2_ewald) 3246 rcut2_ewald = rcut2_ewald*rcut2_ewald 3247 DO i1 = 1, nkinds 3248 atomic_kind => atomic_kind_set(i1) 3249 CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius1) 3250 DO i2 = 1, nkinds 3251 rcut2 = rcut2_ewald 3252 IF (ASSOCIATED(potparm_nonbond)) THEN 3253 rcut2 = MAX(potparm_nonbond%pot(i1, i2)%pot%rcutsq, rcut2_ewald) 3254 END IF 3255 IF (rcut2 > 0) THEN 3256 atomic_kind => atomic_kind_set(i2) 3257 CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius2) 3258 ! cutoff for core-core 3259 interaction_cutoffs(1, i1, i2) = potential_coulomb(rcut2, tmp, & 3260 1.0_dp, ewald_type, alpha, 0.0_dp, 0.0_dp) 3261 ! cutoff for core-shell, core-ion, shell-core or ion-core 3262 IF (mm_radius1 > 0.0_dp) THEN 3263 beta = sqrthalf/mm_radius1 3264 ELSE 3265 beta = 0.0_dp 3266 END IF 3267 interaction_cutoffs(2, i1, i2) = potential_coulomb(rcut2, tmp, & 3268 1.0_dp, ewald_type, alpha, beta, 0.0_dp) 3269 ! cutoff for shell-shell or ion-ion 3270 IF (mm_radius1 + mm_radius2 > 0.0_dp) THEN 3271 beta = sqrthalf/SQRT(mm_radius1*mm_radius1 + mm_radius2*mm_radius2) 3272 ELSE 3273 beta = 0.0_dp 3274 END IF 3275 interaction_cutoffs(3, i1, i2) = potential_coulomb(rcut2, tmp, & 3276 1.0_dp, ewald_type, alpha, beta, 0.0_dp) 3277 END IF 3278 END DO 3279 END DO 3280 END IF 3281 3282 CALL ewald_env_set(ewald_env, interaction_cutoffs=interaction_cutoffs) 3283 3284 CALL timestop(handle) 3285 END SUBROUTINE force_field_pack_eicut 3286 3287! ************************************************************************************************** 3288!> \brief Issues on screen a warning when repetitions are present in the 3289!> definition of the forcefield 3290!> \param found ... 3291!> \param tag_label ... 3292!> \param name_atm_a ... 3293!> \param name_atm_b ... 3294!> \param name_atm_c ... 3295!> \param name_atm_d ... 3296!> \author Teodoro Laino [tlaino] - University of Zurich 10.2008 3297! ************************************************************************************************** 3298 SUBROUTINE issue_duplications(found, tag_label, name_atm_a, name_atm_b, & 3299 name_atm_c, name_atm_d) 3300 3301 LOGICAL, INTENT(IN) :: found 3302 CHARACTER(LEN=*), INTENT(IN) :: tag_label, name_atm_a 3303 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: name_atm_b, name_atm_c, name_atm_d 3304 3305 CHARACTER(len=*), PARAMETER :: routineN = 'issue_duplications', & 3306 routineP = moduleN//':'//routineN 3307 3308 CHARACTER(LEN=default_string_length) :: item 3309 3310 item = "( "//TRIM(name_atm_a) 3311 IF (PRESENT(name_atm_b)) THEN 3312 item = TRIM(item)//" , "//TRIM(name_atm_b) 3313 END IF 3314 IF (PRESENT(name_atm_c)) THEN 3315 item = TRIM(item)//" , "//TRIM(name_atm_c) 3316 END IF 3317 IF (PRESENT(name_atm_d)) THEN 3318 item = TRIM(item)//" , "//TRIM(name_atm_d) 3319 END IF 3320 item = TRIM(item)//" )" 3321 IF (found) & 3322 CPWARN("Multiple "//TRIM(tag_label)//" declarations: "//TRIM(item)//" overwriting! ") 3323 3324 END SUBROUTINE issue_duplications 3325 3326! ************************************************************************************************** 3327!> \brief Store informations on possible missing ForceFields parameters 3328!> \param atm1 ... 3329!> \param atm2 ... 3330!> \param atm3 ... 3331!> \param atm4 ... 3332!> \param type_name ... 3333!> \param fatal ... 3334!> \param array ... 3335! ************************************************************************************************** 3336 SUBROUTINE store_FF_missing_par(atm1, atm2, atm3, atm4, type_name, fatal, array) 3337 CHARACTER(LEN=*), INTENT(IN) :: atm1 3338 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: atm2, atm3, atm4 3339 CHARACTER(LEN=*), INTENT(IN) :: type_name 3340 LOGICAL, INTENT(INOUT), OPTIONAL :: fatal 3341 CHARACTER(LEN=default_string_length), & 3342 DIMENSION(:), POINTER :: array 3343 3344 CHARACTER(len=*), PARAMETER :: routineN = 'store_FF_missing_par', & 3345 routineP = moduleN//':'//routineN 3346 3347 CHARACTER(LEN=10) :: sfmt 3348 CHARACTER(LEN=4) :: my_atm1, my_atm2, my_atm3, my_atm4 3349 CHARACTER(LEN=default_path_length) :: my_format 3350 INTEGER :: fmt, i, nsize 3351 LOGICAL :: found 3352 3353 nsize = 0 3354 fmt = 1 3355 my_format = '(T2,"FORCEFIELD| Missing ","'//TRIM(type_name)// & 3356 '",T40,"(",A4,")")' 3357 IF (PRESENT(atm2)) fmt = fmt + 1 3358 IF (PRESENT(atm3)) fmt = fmt + 1 3359 IF (PRESENT(atm4)) fmt = fmt + 1 3360 CALL integer_to_string(fmt - 1, sfmt) 3361 IF (fmt > 1) & 3362 my_format = '(T2,"FORCEFIELD| Missing ","'//TRIM(type_name)// & 3363 '",T40,"(",A4,'//TRIM(sfmt)//'(",",A4),")")' 3364 IF (PRESENT(fatal)) fatal = .TRUE. 3365 ! Check for previous already stored equal force fields 3366 IF (ASSOCIATED(array)) nsize = SIZE(array) 3367 found = .FALSE. 3368 IF (nsize >= 1) THEN 3369 DO i = 1, nsize 3370 SELECT CASE (type_name) 3371 CASE ("Bond") 3372 IF (INDEX(array(i) (21:39), "Bond") == 0) CYCLE 3373 my_atm1 = array(i) (41:44) 3374 my_atm2 = array(i) (46:49) 3375 CALL compress(my_atm1, .TRUE.) 3376 CALL compress(my_atm2, .TRUE.) 3377 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2)) .OR. & 3378 ((atm1 == my_atm2) .AND. (atm2 == my_atm1))) found = .TRUE. 3379 CASE ("Angle") 3380 IF (INDEX(array(i) (21:39), "Angle") == 0) CYCLE 3381 my_atm1 = array(i) (41:44) 3382 my_atm2 = array(i) (46:49) 3383 my_atm3 = array(i) (51:54) 3384 CALL compress(my_atm1, .TRUE.) 3385 CALL compress(my_atm2, .TRUE.) 3386 CALL compress(my_atm3, .TRUE.) 3387 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3)) .OR. & 3388 ((atm1 == my_atm3) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm1))) & 3389 found = .TRUE. 3390 CASE ("Urey-Bradley") 3391 IF (INDEX(array(i) (21:39), "Urey-Bradley") == 0) CYCLE 3392 my_atm1 = array(i) (41:44) 3393 my_atm2 = array(i) (46:49) 3394 my_atm3 = array(i) (51:54) 3395 CALL compress(my_atm1, .TRUE.) 3396 CALL compress(my_atm2, .TRUE.) 3397 CALL compress(my_atm3, .TRUE.) 3398 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3)) .OR. & 3399 ((atm1 == my_atm3) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm1))) & 3400 found = .TRUE. 3401 CASE ("Torsion") 3402 IF (INDEX(array(i) (21:39), "Torsion") == 0) CYCLE 3403 my_atm1 = array(i) (41:44) 3404 my_atm2 = array(i) (46:49) 3405 my_atm3 = array(i) (51:54) 3406 my_atm4 = array(i) (56:59) 3407 CALL compress(my_atm1, .TRUE.) 3408 CALL compress(my_atm2, .TRUE.) 3409 CALL compress(my_atm3, .TRUE.) 3410 CALL compress(my_atm4, .TRUE.) 3411 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm4)) .OR. & 3412 ((atm1 == my_atm4) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm1))) & 3413 found = .TRUE. 3414 CASE ("Improper") 3415 IF (INDEX(array(i) (21:39), "Improper") == 0) CYCLE 3416 my_atm1 = array(i) (41:44) 3417 my_atm2 = array(i) (46:49) 3418 my_atm3 = array(i) (51:54) 3419 my_atm4 = array(i) (56:59) 3420 CALL compress(my_atm1, .TRUE.) 3421 CALL compress(my_atm2, .TRUE.) 3422 CALL compress(my_atm3, .TRUE.) 3423 CALL compress(my_atm4, .TRUE.) 3424 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm4)) .OR. & 3425 ((atm1 == my_atm1) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm4)) .OR. & 3426 ((atm1 == my_atm1) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm4) .AND. (atm4 == my_atm3)) .OR. & 3427 ((atm1 == my_atm1) .AND. (atm2 == my_atm4) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm2)) .OR. & 3428 ((atm1 == my_atm1) .AND. (atm2 == my_atm4) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm3)) .OR. & 3429 ((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm4) .AND. (atm4 == my_atm3))) & 3430 found = .TRUE. 3431 3432 CASE ("Out of plane bend") 3433 IF (INDEX(array(i) (21:39), "Out of plane bend") == 0) CYCLE 3434 my_atm1 = array(i) (41:44) 3435 my_atm2 = array(i) (46:49) 3436 my_atm3 = array(i) (51:54) 3437 my_atm4 = array(i) (56:59) 3438 CALL compress(my_atm1, .TRUE.) 3439 CALL compress(my_atm2, .TRUE.) 3440 CALL compress(my_atm3, .TRUE.) 3441 CALL compress(my_atm4, .TRUE.) 3442 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2) .AND. (atm3 == my_atm3) .AND. (atm4 == my_atm4)) .OR. & 3443 ((atm1 == my_atm1) .AND. (atm2 == my_atm3) .AND. (atm3 == my_atm2) .AND. (atm4 == my_atm4))) & 3444 found = .TRUE. 3445 3446 CASE ("Charge") 3447 IF (INDEX(array(i) (21:39), "Charge") == 0) CYCLE 3448 my_atm1 = array(i) (41:44) 3449 CALL compress(my_atm1, .TRUE.) 3450 IF (atm1 == my_atm1) found = .TRUE. 3451 CASE ("Spline_Bond_Env", "Spline_Non_Bond_Env") 3452 IF (INDEX(array(i) (21:39), "Spline_") == 0) CYCLE 3453 fmt = 0 3454 my_atm1 = array(i) (41:44) 3455 my_atm2 = array(i) (46:49) 3456 CALL compress(my_atm1, .TRUE.) 3457 CALL compress(my_atm2, .TRUE.) 3458 IF (((atm1 == my_atm1) .AND. (atm2 == my_atm2)) .OR. & 3459 ((atm1 == my_atm2) .AND. (atm2 == my_atm1))) found = .TRUE. 3460 CASE DEFAULT 3461 ! Should never reach this point 3462 CPABORT("") 3463 END SELECT 3464 IF (found) EXIT 3465 END DO 3466 ENDIF 3467 IF (.NOT. found) THEN 3468 nsize = nsize + 1 3469 CALL reallocate(array, 1, nsize) 3470 SELECT CASE (fmt) 3471 CASE (1) 3472 WRITE (array(nsize), FMT=TRIM(my_format)) atm1 3473 CASE (2) 3474 WRITE (array(nsize), FMT=TRIM(my_format)) atm1, atm2 3475 CASE (3) 3476 WRITE (array(nsize), FMT=TRIM(my_format)) atm1, atm2, atm3 3477 CASE (4) 3478 WRITE (array(nsize), FMT=TRIM(my_format)) atm1, atm2, atm3, atm4 3479 END SELECT 3480 END IF 3481 3482 END SUBROUTINE store_FF_missing_par 3483 3484END MODULE force_fields_all 3485 3486