1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Collection of subroutine needed for topology related things 8!> \par History 9!> jgh (23-05-2004) Last atom of molecule information added 10! ************************************************************************************************** 11 12MODULE topology_connectivity_util 13 USE cp_log_handling, ONLY: cp_get_default_logger,& 14 cp_logger_get_default_io_unit,& 15 cp_logger_type 16 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 17 cp_print_key_unit_nr 18 USE force_field_kind_types, ONLY: do_ff_charmm,& 19 do_ff_harmonic 20 USE input_constants, ONLY: do_conn_g87,& 21 do_conn_g96,& 22 do_conn_user 23 USE input_section_types, ONLY: section_vals_type,& 24 section_vals_val_get 25 USE kinds, ONLY: default_string_length 26 USE memory_utilities, ONLY: reallocate 27 USE molecule_kind_types, ONLY: & 28 allocate_molecule_kind_set, atom_type, bend_type, bond_type, get_molecule_kind, impr_type, & 29 molecule_kind_type, opbend_type, set_molecule_kind, torsion_type, ub_type 30 USE molecule_types, ONLY: allocate_molecule_set,& 31 get_molecule,& 32 local_states_type,& 33 molecule_type,& 34 set_molecule,& 35 set_molecule_set 36 USE string_table, ONLY: id2str 37 USE topology_types, ONLY: atom_info_type,& 38 connectivity_info_type,& 39 topology_parameters_type 40 USE util, ONLY: sort 41#include "./base/base_uses.f90" 42 43 IMPLICIT NONE 44 45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_connectivity_util' 46 47 PRIVATE 48 PUBLIC :: topology_connectivity_pack, topology_conn_multiple 49 50CONTAINS 51 52! ************************************************************************************************** 53!> \brief topology connectivity pack 54!> \param molecule_kind_set ... 55!> \param molecule_set ... 56!> \param topology ... 57!> \param subsys_section ... 58!> \par History 11/2009 (Louis Vanduyhuys): added Out of Plane bends based on 59!> impropers in topology 60! ************************************************************************************************** 61 SUBROUTINE topology_connectivity_pack(molecule_kind_set, molecule_set, & 62 topology, subsys_section) 63 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 64 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 65 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 66 TYPE(section_vals_type), POINTER :: subsys_section 67 68 CHARACTER(len=*), PARAMETER :: routineN = 'topology_connectivity_pack', & 69 routineP = moduleN//':'//routineN 70 71 CHARACTER(LEN=default_string_length) :: name 72 INTEGER :: counter, first, handle, handle2, i, ibend, ibond, idim, iimpr, ikind, imol, & 73 inter_bends, inter_bonds, inter_imprs, inter_torsions, inter_ubs, intra_bends, & 74 intra_bonds, intra_imprs, intra_torsions, intra_ubs, inum, ires, istart_mol, istart_typ, & 75 itorsion, ityp, iub, iw, j, j1, j2, j3, j4, jind, last, min_index, natom, nelectron, & 76 nsgf, nval_tot1, nval_tot2, nvar1, nvar2, output_unit, stat 77 INTEGER, DIMENSION(:), POINTER :: c_var_a, c_var_b, c_var_c, c_var_d, c_var_type, & 78 first_list, last_list, map_atom_mol, map_atom_type, map_cvar_mol, map_cvars, map_var_mol, & 79 map_vars, molecule_list 80 INTEGER, DIMENSION(:, :), POINTER :: bnd_ctype, bnd_type 81 LOGICAL :: found, found_last 82 TYPE(atom_info_type), POINTER :: atom_info 83 TYPE(atom_type), DIMENSION(:), POINTER :: atom_list 84 TYPE(bend_type), DIMENSION(:), POINTER :: bend_list 85 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list 86 TYPE(connectivity_info_type), POINTER :: conn_info 87 TYPE(cp_logger_type), POINTER :: logger 88 TYPE(impr_type), DIMENSION(:), POINTER :: impr_list 89 TYPE(local_states_type), DIMENSION(:), POINTER :: lmi 90 TYPE(molecule_kind_type), POINTER :: molecule_kind 91 TYPE(molecule_type), POINTER :: molecule 92 TYPE(opbend_type), DIMENSION(:), POINTER :: opbend_list 93 TYPE(torsion_type), DIMENSION(:), POINTER :: torsion_list 94 TYPE(ub_type), DIMENSION(:), POINTER :: ub_list 95 96 NULLIFY (logger) 97 logger => cp_get_default_logger() 98 output_unit = cp_logger_get_default_io_unit(logger) 99 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", & 100 extension=".subsysLog") 101 CALL timeset(routineN, handle) 102 103 atom_info => topology%atom_info 104 conn_info => topology%conn_info 105 ALLOCATE (map_atom_mol(topology%natoms)) 106 ALLOCATE (map_atom_type(topology%natoms)) 107 !----------------------------------------------------------------------------- 108 !----------------------------------------------------------------------------- 109 ! 1. Set the topology%[nmol_type,nmol,nmol_conn] 110 !----------------------------------------------------------------------------- 111 CALL timeset(routineN//"_1", handle2) 112 natom = topology%natoms 113 topology%nmol = 1 114 topology%nmol_type = 1 115 topology%nmol_conn = 0 116 map_atom_mol = -1 117 map_atom_type = -1 118 map_atom_mol(1) = 1 119 map_atom_type(1) = 1 120 DO i = 1, natom - 1 121 IF ((atom_info%map_mol_typ(i + 1) /= atom_info%map_mol_typ(i)) .OR. & 122 ((atom_info%map_mol_res(i + 1) /= atom_info%map_mol_res(i)) .AND. & 123 (.NOT. (topology%conn_type == do_conn_user)))) THEN 124 topology%nmol_type = topology%nmol_type + 1 125 END IF 126 map_atom_type(i + 1) = topology%nmol_type 127 IF ((atom_info%map_mol_typ(i + 1) /= atom_info%map_mol_typ(i)) .OR. & 128 (atom_info%map_mol_num(i + 1) /= atom_info%map_mol_num(i)) .OR. & 129 (atom_info%map_mol_res(i + 1) /= atom_info%map_mol_res(i))) THEN 130 topology%nmol = topology%nmol + 1 131 END IF 132 map_atom_mol(i + 1) = topology%nmol 133 IF ((atom_info%map_mol_typ(i + 1) == atom_info%map_mol_typ(i)) .AND. & 134 (atom_info%map_mol_num(i + 1) == atom_info%map_mol_num(i)) .AND. & 135 (atom_info%map_mol_res(i + 1) /= atom_info%map_mol_res(i))) THEN 136 topology%nmol_conn = topology%nmol_conn + 1 137 END IF 138 END DO 139 IF (iw > 0) WRITE (iw, *) "topology%nmol ::", topology%nmol 140 IF (iw > 0) WRITE (iw, *) "topology%nmol_type ::", topology%nmol_type 141 IF (iw > 0) WRITE (iw, *) "topology%nmol_conn ::", topology%nmol_conn 142 CALL timestop(handle2) 143 !----------------------------------------------------------------------------- 144 !----------------------------------------------------------------------------- 145 ! 1.1 Clean the temporary arrays to avoid quadratic loops around 146 ! after this fix all topology_pack will be linear scaling 147 !----------------------------------------------------------------------------- 148 CALL timeset(routineN//"_1.1", handle2) 149 istart_mol = map_atom_mol(1) 150 istart_typ = map_atom_type(1) 151 DO i = 2, natom 152 IF ((map_atom_mol(i) /= istart_mol) .AND. (map_atom_type(i) == istart_typ)) THEN 153 map_atom_mol(i) = -map_atom_mol(i) 154 ELSE IF (map_atom_type(i) /= istart_typ) THEN 155 istart_mol = map_atom_mol(i) 156 istart_typ = map_atom_type(i) 157 END IF 158 END DO 159 CALL timestop(handle2) 160 !----------------------------------------------------------------------------- 161 !----------------------------------------------------------------------------- 162 ! 2. Allocate the molecule_kind_set 163 !----------------------------------------------------------------------------- 164 CALL timeset(routineN//"_2", handle2) 165 IF (topology%nmol_type <= 0) THEN 166 CPABORT("No molecule kind defined") 167 ELSE 168 NULLIFY (molecule_kind_set) 169 i = topology%nmol_type 170 CALL allocate_molecule_kind_set(molecule_kind_set, i) 171 IF (iw > 0) WRITE (iw, *) " Allocated molecule_kind_set, Dimenstion of ", & 172 SIZE(molecule_kind_set) 173 END IF 174 CALL timestop(handle2) 175 !----------------------------------------------------------------------------- 176 !----------------------------------------------------------------------------- 177 ! 3. Allocate the molecule_set 178 !----------------------------------------------------------------------------- 179 CALL timeset(routineN//"_3", handle2) 180 IF (topology%nmol <= 0) THEN 181 CPABORT("No molecule defined") 182 ELSE 183 NULLIFY (molecule_set) 184 i = topology%nmol 185 CALL allocate_molecule_set(molecule_set, i) 186 IF (iw > 0) WRITE (iw, *) " Allocated molecule_set, dimension of ", & 187 topology%nmol 188 END IF 189 CALL timestop(handle2) 190 !----------------------------------------------------------------------------- 191 !----------------------------------------------------------------------------- 192 ! 4. Set the molecule_kind_set%[kind_number,name,nsgf,nelectron] 193 !----------------------------------------------------------------------------- 194 CALL timeset(routineN//"_4", handle2) 195 natom = topology%natoms 196 ikind = -1 197 DO i = 1, natom 198 IF (ikind /= map_atom_type(i)) THEN 199 ikind = map_atom_type(i) 200 molecule_kind => molecule_kind_set(ikind) 201 nsgf = 0 202 nelectron = 0 203 name = TRIM(id2str(atom_info%id_molname(i))) 204 CALL set_molecule_kind(molecule_kind=molecule_kind, & 205 kind_number=ikind, & 206 molname_generated=topology%molname_generated, & 207 name=TRIM(name), & 208 nsgf=nsgf, & 209 nelectron=nelectron) 210 END IF 211 END DO 212 CALL timestop(handle2) 213 !----------------------------------------------------------------------------- 214 !----------------------------------------------------------------------------- 215 ! 5. Set the molecule_list for molecule_kind in molecule_kind_set 216 !----------------------------------------------------------------------------- 217 CALL timeset(routineN//"_5", handle2) 218 natom = topology%natoms 219 ikind = map_atom_type(1) 220 imol = ABS(map_atom_mol(1)) 221 counter = 0 222 DO i = 1, natom - 1 223 IF (ikind /= map_atom_type(i + 1)) THEN 224 found = .TRUE. 225 found_last = .FALSE. 226 imol = ABS(map_atom_mol(i)) 227 ELSEIF (ikind == topology%nmol_type) THEN 228 found = .TRUE. 229 found_last = .TRUE. 230 imol = ABS(map_atom_mol(natom)) 231 ELSE 232 found = .FALSE. 233 found_last = .FALSE. 234 END IF 235 236 IF (found) THEN 237 ALLOCATE (molecule_list(imol - counter)) 238 DO j = 1, SIZE(molecule_list) 239 molecule_list(j) = j + counter 240 END DO 241 molecule_kind => molecule_kind_set(ikind) 242 CALL set_molecule_kind(molecule_kind=molecule_kind, & 243 molecule_list=molecule_list) 244 IF (iw > 0) WRITE (iw, *) " molecule_list", ikind, molecule_list(:) 245 IF (found_last) EXIT 246 counter = imol 247 ikind = map_atom_type(i + 1) 248 END IF 249 END DO 250 ! Treat separately the case in which the last atom is also a molecule 251 IF (i == natom) THEN 252 imol = ABS(map_atom_mol(natom)) 253 ! Last atom is also a molecule by itself 254 ALLOCATE (molecule_list(imol - counter)) 255 DO j = 1, SIZE(molecule_list) 256 molecule_list(j) = j + counter 257 END DO 258 molecule_kind => molecule_kind_set(ikind) 259 CALL set_molecule_kind(molecule_kind=molecule_kind, & 260 molecule_list=molecule_list) 261 IF (iw > 0) WRITE (iw, *) " molecule_list", ikind, molecule_list(:) 262 END IF 263 CALL timestop(handle2) 264 !----------------------------------------------------------------------------- 265 !----------------------------------------------------------------------------- 266 ! 6. Set the molecule_set(imol)%molecule_kind via set_molecule 267 !----------------------------------------------------------------------------- 268 CALL timeset(routineN//"_6", handle2) 269 DO ikind = 1, SIZE(molecule_kind_set) 270 molecule_kind => molecule_kind_set(ikind) 271 CALL get_molecule_kind(molecule_kind=molecule_kind, & 272 molecule_list=molecule_list) 273 DO i = 1, SIZE(molecule_list) 274 molecule => molecule_set(molecule_list(i)) 275 CALL set_molecule(molecule, molecule_kind=molecule_kind) 276 END DO 277 END DO 278 CALL timestop(handle2) 279 !----------------------------------------------------------------------------- 280 !----------------------------------------------------------------------------- 281 ! 7. Set the molecule_set(imol)%[first_atom,last_atom] via set_molecule_set 282 !----------------------------------------------------------------------------- 283 ALLOCATE (first_list(SIZE(molecule_set))) 284 ALLOCATE (last_list(SIZE(molecule_set))) 285 CALL timeset(routineN//"_7", handle2) 286 first_list(:) = 0 287 last_list(:) = 0 288 ityp = atom_info%map_mol_typ(1) 289 inum = atom_info%map_mol_num(1) 290 ires = atom_info%map_mol_res(1) 291 imol = 1 292 first_list(1) = 1 293 DO j = 2, natom 294 IF ((atom_info%map_mol_typ(j) /= ityp) .OR. & 295 (atom_info%map_mol_num(j) /= inum) .OR. & 296 (atom_info%map_mol_res(j) /= ires)) THEN 297 ityp = atom_info%map_mol_typ(j) 298 inum = atom_info%map_mol_num(j) 299 ires = atom_info%map_mol_res(j) 300 imol = imol + 1 301 first_list(imol) = j 302 END IF 303 END DO 304 CPASSERT(imol == topology%nmol) 305 DO ikind = 1, topology%nmol - 1 306 last_list(ikind) = first_list(ikind + 1) - 1 307 END DO 308 last_list(topology%nmol) = topology%natoms 309 CALL set_molecule_set(molecule_set, first_list, last_list) 310 CALL timestop(handle2) 311 !----------------------------------------------------------------------------- 312 !----------------------------------------------------------------------------- 313 ! 8. Set and NULLIFY the molecule_set(imol)%lmi via set_molecule_set 314 !----------------------------------------------------------------------------- 315 CALL timeset(routineN//"_8", handle2) 316 DO imol = 1, SIZE(molecule_set) 317 molecule => molecule_set(imol) 318 NULLIFY (lmi) 319 CALL set_molecule(molecule, lmi=lmi) 320 END DO 321 CALL timestop(handle2) 322 !----------------------------------------------------------------------------- 323 !----------------------------------------------------------------------------- 324 ! 9. Set the atom_list for molecule_kind in molecule_kind_set (PART 1) 325 !----------------------------------------------------------------------------- 326 CALL timeset(routineN//"_9", handle2) 327 counter = 0 328 DO imol = 1, SIZE(molecule_set) 329 molecule => molecule_set(imol) 330 molecule_kind => molecule_set(imol)%molecule_kind 331 CALL get_molecule_kind(molecule_kind=molecule_kind, & 332 kind_number=i) 333 IF (counter /= i) THEN 334 counter = i 335 CALL get_molecule(molecule=molecule, & 336 first_atom=first, last_atom=last) 337 natom = 0 338 IF (first /= 0 .AND. last /= 0) natom = last - first + 1 339 ALLOCATE (atom_list(natom)) 340 DO i = 1, natom 341 !Atomic kind information will be filled in (PART 2) 342 NULLIFY (atom_list(i)%atomic_kind) 343 atom_list(i)%id_name = atom_info%id_atmname(i + first - 1) 344 IF (iw > 0) WRITE (iw, '(5X,A,3I5,1X,A5)') "atom_list ", & 345 imol, counter, i, TRIM(id2str(atom_list(i)%id_name)) 346 END DO 347 CALL set_molecule_kind(molecule_kind=molecule_kind, atom_list=atom_list) 348 END IF 349 END DO 350 CALL timestop(handle2) 351 !----------------------------------------------------------------------------- 352 !----------------------------------------------------------------------------- 353 ! 10. Set the molecule_kind%[nbond,bond_list] via set_molecule_kind 354 !----------------------------------------------------------------------------- 355 CALL timeset(routineN//"_10", handle2) 356 ! First map bonds on molecules 357 nvar1 = 0 358 nvar2 = 0 359 NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype) 360 IF (ASSOCIATED(conn_info%bond_a)) nvar1 = SIZE(conn_info%bond_a) 361 IF (ASSOCIATED(conn_info%c_bond_a)) nvar2 = SIZE(conn_info%c_bond_a) 362 nval_tot1 = nvar1 363 nval_tot2 = 0 364 ALLOCATE (map_var_mol(nvar1)) 365 ALLOCATE (map_cvar_mol(nvar2)) 366 map_var_mol = -1 367 map_cvar_mol = -1 368 DO i = 1, nvar1 369 j1 = map_atom_mol(conn_info%bond_a(i)) 370 j2 = map_atom_mol(conn_info%bond_b(i)) 371 IF (j1 == j2) THEN 372 IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%bond_a(i)) 373 END IF 374 END DO 375 DO i = 1, nvar2 376 min_index = MIN(conn_info%c_bond_a(i), conn_info%c_bond_b(i)) 377 j1 = map_atom_mol(min_index) 378 IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index) 379 END DO 380 CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1) 381 CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2) 382 DO i = 1, topology%nmol_type 383 intra_bonds = 0 384 inter_bonds = 0 385 IF (ALL(bnd_type(:, i) > 0)) THEN 386 intra_bonds = bnd_type(2, i) - bnd_type(1, i) + 1 387 END IF 388 IF (ALL(bnd_ctype(:, i) > 0)) THEN 389 inter_bonds = bnd_ctype(2, i) - bnd_ctype(1, i) + 1 390 END IF 391 ibond = intra_bonds + inter_bonds 392 IF (iw > 0) THEN 393 WRITE (iw, *) " Total number bonds for molecule type ", i, " :", ibond 394 WRITE (iw, *) " intra (bonds inside molecules) :: ", intra_bonds 395 WRITE (iw, *) " inter (bonds between molecules) :: ", inter_bonds 396 END IF 397 molecule_kind => molecule_kind_set(i) 398 nval_tot2 = nval_tot2 + ibond*SIZE(molecule_kind%molecule_list) 399 400 ALLOCATE (bond_list(ibond)) 401 ibond = 0 402 DO j = bnd_type(1, i), bnd_type(2, i) 403 IF (j == 0) CYCLE 404 ibond = ibond + 1 405 jind = map_vars(j) 406 first = first_list(map_atom_mol(conn_info%bond_a(jind))) 407 bond_list(ibond)%a = conn_info%bond_a(jind) - first + 1 408 bond_list(ibond)%b = conn_info%bond_b(jind) - first + 1 409 ! Set by default id_type to charmm and modify when handling the forcefield 410 bond_list(ibond)%id_type = do_ff_charmm 411 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN 412 bond_list(ibond)%itype = conn_info%bond_type(jind) 413 END IF 414 !point this to the right bond_kind_type if using force field 415 NULLIFY (bond_list(ibond)%bond_kind) 416 IF (iw > 0) THEN 417 WRITE (iw, '(7X,A,I3,1X,A,I5,I5,1X,A,I5,I5)') "molecule_kind", & 418 i, " intra bond", & 419 conn_info%bond_a(jind), & 420 conn_info%bond_b(jind), & 421 "offset number at", & 422 conn_info%bond_a(jind) - first + 1, & 423 conn_info%bond_b(jind) - first + 1 424 END IF 425 END DO 426 DO j = bnd_ctype(1, i), bnd_ctype(2, i) 427 IF (j == 0) CYCLE 428 ibond = ibond + 1 429 jind = map_cvars(j) 430 min_index = MIN(conn_info%c_bond_a(jind), conn_info%c_bond_b(jind)) 431 first = first_list(map_atom_mol(min_index)) 432 bond_list(ibond)%a = conn_info%c_bond_a(jind) - first + 1 433 bond_list(ibond)%b = conn_info%c_bond_b(jind) - first + 1 434 ! Set by default id_type to charmm and modify when handling the forcefield 435 bond_list(ibond)%id_type = do_ff_charmm 436 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN 437 bond_list(ibond)%itype = conn_info%c_bond_type(jind) 438 END IF 439 !point this to the right bond_kind_type if using force field 440 NULLIFY (bond_list(ibond)%bond_kind) 441 IF (iw > 0) THEN 442 WRITE (iw, '(7X,A,I3,1X,A,I5,I5,1X,A,I5,I5)') "molecule_kind", & 443 i, " inter bond", & 444 conn_info%c_bond_a(jind), & 445 conn_info%c_bond_b(jind), & 446 "offset number at", & 447 conn_info%c_bond_a(jind) - first + 1, & 448 conn_info%c_bond_b(jind) - first + 1 449 END IF 450 END DO 451 CALL set_molecule_kind(molecule_kind=molecule_kind, & 452 nbond=SIZE(bond_list), bond_list=bond_list) 453 END DO 454 IF ((nval_tot1 /= nval_tot2) .AND. (output_unit > 0)) THEN 455 WRITE (output_unit, '(/)') 456 WRITE (output_unit, '(T5,A)') "ERROR| Mismatching found between the total number of atoms" 457 WRITE (output_unit, '(T5,A)') "ERROR| and the number of atoms computed multiplying the Nr." 458 WRITE (output_unit, '(T5,A)') "ERROR| of molecules by the number of atoms building that" 459 WRITE (output_unit, '(T5,A)') "ERROR| kind of molecule." 460 WRITE (output_unit, '(T5,A)') "ERROR| This happens when the connectivity is wrongly built" 461 WRITE (output_unit, '(T5,A)') "ERROR| One example could be two same kind of molecules have" 462 WRITE (output_unit, '(T5,A)') "ERROR| a different number of atoms. Check the connectivity!" 463 END IF 464 CPASSERT(nval_tot1 == nval_tot2) 465 DEALLOCATE (map_var_mol) 466 DEALLOCATE (map_cvar_mol) 467 DEALLOCATE (map_vars) 468 DEALLOCATE (map_cvars) 469 DEALLOCATE (bnd_type) 470 DEALLOCATE (bnd_ctype) 471 CALL timestop(handle2) 472 !----------------------------------------------------------------------------- 473 !----------------------------------------------------------------------------- 474 ! 11. Set the molecule_kind%[nbend,bend_list] via set_molecule_kind 475 !----------------------------------------------------------------------------- 476 ! Allocate c_var_a, c_var_b, c_var_c, c_var_type 477 CALL timeset(routineN//"_11_pre", handle2) 478 idim = 0 479 ALLOCATE (c_var_a(idim)) 480 ALLOCATE (c_var_b(idim)) 481 ALLOCATE (c_var_c(idim)) 482 found = ASSOCIATED(conn_info%theta_type) 483 IF (found) THEN 484 ALLOCATE (c_var_type(idim)) 485 END IF 486 IF (ASSOCIATED(conn_info%c_bond_a) .AND. ASSOCIATED(conn_info%theta_a)) THEN 487 DO j = 1, SIZE(conn_info%theta_a) 488 j1 = map_atom_mol(conn_info%theta_a(j)) 489 j2 = map_atom_mol(conn_info%theta_b(j)) 490 j3 = map_atom_mol(conn_info%theta_c(j)) 491 IF (j1 /= j2 .OR. j2 /= j3) THEN 492 idim = idim + 1 493 END IF 494 END DO 495 CALL reallocate(c_var_a, 1, idim) 496 CALL reallocate(c_var_b, 1, idim) 497 CALL reallocate(c_var_c, 1, idim) 498 IF (found) THEN 499 CALL reallocate(c_var_type, 1, idim) 500 END IF 501 idim = 0 502 DO j = 1, SIZE(conn_info%theta_a) 503 j1 = map_atom_mol(conn_info%theta_a(j)) 504 j2 = map_atom_mol(conn_info%theta_b(j)) 505 j3 = map_atom_mol(conn_info%theta_c(j)) 506 IF (j1 /= j2 .OR. j2 /= j3) THEN 507 idim = idim + 1 508 c_var_a(idim) = conn_info%theta_a(j) 509 c_var_b(idim) = conn_info%theta_b(j) 510 c_var_c(idim) = conn_info%theta_c(j) 511 IF (found) THEN 512 c_var_type(idim) = conn_info%theta_type(j) 513 END IF 514 END IF 515 END DO 516 END IF 517 CALL timestop(handle2) 518 CALL timeset(routineN//"_11", handle2) 519 ! map bends on molecules 520 nvar1 = 0 521 nvar2 = 0 522 NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype) 523 IF (ASSOCIATED(conn_info%theta_a)) nvar1 = SIZE(conn_info%theta_a) 524 IF (ASSOCIATED(c_var_a)) nvar2 = SIZE(c_var_a) 525 nval_tot1 = nvar1 526 nval_tot2 = 0 527 ALLOCATE (map_var_mol(nvar1)) 528 ALLOCATE (map_cvar_mol(nvar2)) 529 map_var_mol = -1 530 map_cvar_mol = -1 531 DO i = 1, nvar1 532 j1 = map_atom_mol(conn_info%theta_a(i)) 533 j2 = map_atom_mol(conn_info%theta_b(i)) 534 j3 = map_atom_mol(conn_info%theta_c(i)) 535 IF (j1 == j2 .AND. j2 == j3) THEN 536 IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%theta_a(i)) 537 END IF 538 END DO 539 DO i = 1, nvar2 540 min_index = MIN(c_var_a(i), c_var_b(i), c_var_c(i)) 541 j1 = map_atom_mol(min_index) 542 IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index) 543 END DO 544 CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1) 545 CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2) 546 DO i = 1, topology%nmol_type 547 intra_bends = 0 548 inter_bends = 0 549 IF (ALL(bnd_type(:, i) > 0)) THEN 550 intra_bends = bnd_type(2, i) - bnd_type(1, i) + 1 551 END IF 552 IF (ALL(bnd_ctype(:, i) > 0)) THEN 553 inter_bends = bnd_ctype(2, i) - bnd_ctype(1, i) + 1 554 END IF 555 ibend = intra_bends + inter_bends 556 IF (iw > 0) THEN 557 WRITE (iw, *) " Total number of angles for molecule type ", i, " :", ibend 558 WRITE (iw, *) " intra (angles inside molecules) :: ", intra_bends 559 WRITE (iw, *) " inter (angles between molecules) :: ", inter_bends 560 END IF 561 molecule_kind => molecule_kind_set(i) 562 nval_tot2 = nval_tot2 + ibend*SIZE(molecule_kind%molecule_list) 563 ALLOCATE (bend_list(ibend)) 564 ibend = 0 565 DO j = bnd_type(1, i), bnd_type(2, i) 566 IF (j == 0) CYCLE 567 ibend = ibend + 1 568 jind = map_vars(j) 569 first = first_list(map_atom_mol(conn_info%theta_a(jind))) 570 bend_list(ibend)%a = conn_info%theta_a(jind) - first + 1 571 bend_list(ibend)%b = conn_info%theta_b(jind) - first + 1 572 bend_list(ibend)%c = conn_info%theta_c(jind) - first + 1 573 ! Set by default id_type to charmm and modify when handling the forcefield 574 bend_list(ibend)%id_type = do_ff_charmm 575 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN 576 bend_list(ibend)%itype = conn_info%theta_type(jind) 577 END IF 578 !point this to the right bend_kind_type if using force field 579 NULLIFY (bend_list(ibend)%bend_kind) 580 IF (iw > 0) THEN 581 WRITE (iw, '(7X,A,I3,1X,A,I5,I5,I5,1X,A,I5,I5,I5)') & 582 "molecule_kind", ikind, "intra bend", & 583 conn_info%theta_a(jind), & 584 conn_info%theta_b(jind), & 585 conn_info%theta_c(jind), & 586 "offset number at", & 587 conn_info%theta_a(jind) - first + 1, & 588 conn_info%theta_b(jind) - first + 1, & 589 conn_info%theta_c(jind) - first + 1 590 END IF 591 END DO 592 DO j = bnd_ctype(1, i), bnd_ctype(2, i) 593 IF (j == 0) CYCLE 594 ibend = ibend + 1 595 jind = map_cvars(j) 596 min_index = MIN(c_var_a(jind), c_var_b(jind), c_var_c(jind)) 597 first = first_list(map_atom_mol(min_index)) 598 bend_list(ibend)%a = c_var_a(jind) - first + 1 599 bend_list(ibend)%b = c_var_b(jind) - first + 1 600 bend_list(ibend)%c = c_var_c(jind) - first + 1 601 ! Set by default id_type to charmm and modify when handling the forcefield 602 bend_list(ibend)%id_type = do_ff_charmm 603 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN 604 bend_list(ibend)%itype = c_var_type(jind) 605 END IF 606 !point this to the right bend_kind_type if using force field 607 NULLIFY (bend_list(ibend)%bend_kind) 608 IF (iw > 0) THEN 609 WRITE (iw, '(7X,A,I3,1X,A,I5,I5,I5,1X,A,I5,I5,I5)') & 610 "molecule_kind", ikind, "inter bend", & 611 c_var_a(jind), & 612 c_var_b(jind), & 613 c_var_c(jind), & 614 "offset number at", & 615 c_var_a(jind) - first + 1, & 616 c_var_b(jind) - first + 1, & 617 c_var_c(jind) - first + 1 618 END IF 619 END DO 620 CALL set_molecule_kind(molecule_kind=molecule_kind, & 621 nbend=SIZE(bend_list), bend_list=bend_list) 622 END DO 623 CPASSERT(nval_tot1 == nval_tot2) 624 DEALLOCATE (map_var_mol) 625 DEALLOCATE (map_cvar_mol) 626 DEALLOCATE (map_vars) 627 DEALLOCATE (map_cvars) 628 DEALLOCATE (bnd_type) 629 DEALLOCATE (bnd_ctype) 630 DEALLOCATE (c_var_a) 631 DEALLOCATE (c_var_b) 632 DEALLOCATE (c_var_c) 633 IF (found) THEN 634 DEALLOCATE (c_var_type) 635 END IF 636 CALL timestop(handle2) 637 !----------------------------------------------------------------------------- 638 !----------------------------------------------------------------------------- 639 ! 12. Set the molecule_kind%[nub,ub_list] via set_molecule_kind 640 !----------------------------------------------------------------------------- 641 CALL timeset(routineN//"_12_pre", handle2) 642 idim = 0 643 ALLOCATE (c_var_a(idim)) 644 ALLOCATE (c_var_b(idim)) 645 ALLOCATE (c_var_c(idim)) 646 IF (ASSOCIATED(conn_info%c_bond_a) .AND. ASSOCIATED(conn_info%ub_a)) THEN 647 DO j = 1, SIZE(conn_info%ub_a) 648 j1 = map_atom_mol(conn_info%ub_a(j)) 649 j2 = map_atom_mol(conn_info%ub_b(j)) 650 j3 = map_atom_mol(conn_info%ub_c(j)) 651 IF (j1 /= j2 .OR. j2 /= j3) THEN 652 idim = idim + 1 653 END IF 654 END DO 655 CALL reallocate(c_var_a, 1, idim) 656 CALL reallocate(c_var_b, 1, idim) 657 CALL reallocate(c_var_c, 1, idim) 658 idim = 0 659 DO j = 1, SIZE(conn_info%ub_a) 660 j1 = map_atom_mol(conn_info%ub_a(j)) 661 j2 = map_atom_mol(conn_info%ub_b(j)) 662 j3 = map_atom_mol(conn_info%ub_c(j)) 663 IF (j1 /= j2 .OR. j2 /= j3) THEN 664 idim = idim + 1 665 c_var_a(idim) = conn_info%ub_a(j) 666 c_var_b(idim) = conn_info%ub_b(j) 667 c_var_c(idim) = conn_info%ub_c(j) 668 END IF 669 END DO 670 END IF 671 CALL timestop(handle2) 672 CALL timeset(routineN//"_12", handle2) 673 ! map UBs on molecules 674 nvar1 = 0 675 nvar2 = 0 676 NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype) 677 IF (ASSOCIATED(conn_info%ub_a)) nvar1 = SIZE(conn_info%ub_a) 678 IF (ASSOCIATED(c_var_a)) nvar2 = SIZE(c_var_a) 679 nval_tot1 = nvar1 680 nval_tot2 = 0 681 ALLOCATE (map_var_mol(nvar1)) 682 ALLOCATE (map_cvar_mol(nvar2)) 683 map_var_mol = -1 684 map_cvar_mol = -1 685 DO i = 1, nvar1 686 j1 = map_atom_mol(conn_info%ub_a(i)) 687 j2 = map_atom_mol(conn_info%ub_b(i)) 688 j3 = map_atom_mol(conn_info%ub_c(i)) 689 IF (j1 == j2 .AND. j2 == j3) THEN 690 IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%ub_a(i)) 691 END IF 692 END DO 693 DO i = 1, nvar2 694 min_index = MIN(c_var_a(i), c_var_b(i), c_var_c(i)) 695 j1 = map_atom_mol(min_index) 696 IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index) 697 END DO 698 CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1) 699 CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2) 700 DO i = 1, topology%nmol_type 701 intra_ubs = 0 702 inter_ubs = 0 703 IF (ALL(bnd_type(:, i) > 0)) THEN 704 intra_ubs = bnd_type(2, i) - bnd_type(1, i) + 1 705 END IF 706 IF (ALL(bnd_ctype(:, i) > 0)) THEN 707 inter_ubs = bnd_ctype(2, i) - bnd_ctype(1, i) + 1 708 END IF 709 iub = intra_ubs + inter_ubs 710 IF (iw > 0) THEN 711 WRITE (iw, *) " Total number of Urey-Bradley for molecule type ", i, " :", iub 712 WRITE (iw, *) " intra (UB inside molecules) :: ", intra_ubs 713 WRITE (iw, *) " inter (UB between molecules) :: ", inter_ubs 714 END IF 715 molecule_kind => molecule_kind_set(i) 716 nval_tot2 = nval_tot2 + iub*SIZE(molecule_kind%molecule_list) 717 ALLOCATE (ub_list(iub)) 718 iub = 0 719 DO j = bnd_type(1, i), bnd_type(2, i) 720 IF (j == 0) CYCLE 721 iub = iub + 1 722 jind = map_vars(j) 723 first = first_list(map_atom_mol(conn_info%ub_a(jind))) 724 ub_list(iub)%a = conn_info%ub_a(jind) - first + 1 725 ub_list(iub)%b = conn_info%ub_b(jind) - first + 1 726 ub_list(iub)%c = conn_info%ub_c(jind) - first + 1 727 ub_list(iub)%id_type = do_ff_charmm 728 !point this to the right ub_kind_type if using force field 729 NULLIFY (ub_list(iub)%ub_kind) 730 IF (iw > 0) THEN 731 WRITE (iw, '(7X,A,I3,1X,A,I5,I5,I5,1X,A,I5,I5,I5)') & 732 "molecule_kind", i, "intra UB", & 733 conn_info%ub_a(jind), & 734 conn_info%ub_b(jind), & 735 conn_info%ub_c(jind), & 736 "offset number at", & 737 conn_info%ub_a(jind) - first + 1, & 738 conn_info%ub_b(jind) - first + 1, & 739 conn_info%ub_c(jind) - first + 1 740 END IF 741 END DO 742 DO j = bnd_ctype(1, i), bnd_ctype(2, i) 743 IF (j == 0) CYCLE 744 iub = iub + 1 745 jind = map_cvars(j) 746 min_index = MIN(c_var_a(jind), c_var_b(jind), c_var_c(jind)) 747 first = first_list(map_atom_mol(min_index)) 748 ub_list(iub)%a = c_var_a(jind) - first + 1 749 ub_list(iub)%b = c_var_b(jind) - first + 1 750 ub_list(iub)%c = c_var_c(jind) - first + 1 751 ub_list(iub)%id_type = do_ff_charmm 752 !point this to the right ub_kind_type if using force field 753 NULLIFY (ub_list(iub)%ub_kind) 754 IF (iw > 0) THEN 755 WRITE (iw, '(7X,A,I3,1X,A,I5,I5,I5,1X,A,I5,I5,I5)') & 756 "molecule_kind", i, "inter UB", & 757 c_var_a(jind), & 758 c_var_b(jind), & 759 c_var_c(jind), & 760 "offset number at", & 761 c_var_a(jind) - first + 1, & 762 c_var_b(jind) - first + 1, & 763 c_var_c(jind) - first + 1 764 END IF 765 END DO 766 CALL set_molecule_kind(molecule_kind=molecule_kind, & 767 nub=SIZE(ub_list), ub_list=ub_list) 768 END DO 769 CPASSERT(nval_tot1 == nval_tot2) 770 DEALLOCATE (map_var_mol) 771 DEALLOCATE (map_cvar_mol) 772 DEALLOCATE (map_vars) 773 DEALLOCATE (map_cvars) 774 DEALLOCATE (bnd_type) 775 DEALLOCATE (bnd_ctype) 776 DEALLOCATE (c_var_a) 777 DEALLOCATE (c_var_b) 778 DEALLOCATE (c_var_c) 779 CALL timestop(handle2) 780 !----------------------------------------------------------------------------- 781 !----------------------------------------------------------------------------- 782 ! 13. Set the molecule_kind%[ntorsion,torsion_list] via set_molecule_kind 783 !----------------------------------------------------------------------------- 784 ! Allocate c_var_a, c_var_b, c_var_c, c_var_d, c_var_type 785 CALL timeset(routineN//"_13_pre", handle2) 786 idim = 0 787 ALLOCATE (c_var_a(idim)) 788 ALLOCATE (c_var_b(idim)) 789 ALLOCATE (c_var_c(idim)) 790 ALLOCATE (c_var_d(idim)) 791 found = ASSOCIATED(conn_info%phi_type) 792 IF (found) THEN 793 ALLOCATE (c_var_type(idim)) 794 END IF 795 IF (ASSOCIATED(conn_info%c_bond_a) .AND. ASSOCIATED(conn_info%phi_a)) THEN 796 DO j = 1, SIZE(conn_info%phi_a) 797 j1 = map_atom_mol(conn_info%phi_a(j)) 798 j2 = map_atom_mol(conn_info%phi_b(j)) 799 j3 = map_atom_mol(conn_info%phi_c(j)) 800 j4 = map_atom_mol(conn_info%phi_d(j)) 801 IF (j1 /= j2 .OR. j2 /= j3 .OR. j3 /= j4) THEN 802 idim = idim + 1 803 END IF 804 END DO 805 CALL reallocate(c_var_a, 1, idim) 806 CALL reallocate(c_var_b, 1, idim) 807 CALL reallocate(c_var_c, 1, idim) 808 CALL reallocate(c_var_d, 1, idim) 809 IF (found) THEN 810 CALL reallocate(c_var_type, 1, idim) 811 END IF 812 idim = 0 813 DO j = 1, SIZE(conn_info%phi_a) 814 j1 = map_atom_mol(conn_info%phi_a(j)) 815 j2 = map_atom_mol(conn_info%phi_b(j)) 816 j3 = map_atom_mol(conn_info%phi_c(j)) 817 j4 = map_atom_mol(conn_info%phi_d(j)) 818 IF (j1 /= j2 .OR. j2 /= j3 .OR. j3 /= j4) THEN 819 idim = idim + 1 820 c_var_a(idim) = conn_info%phi_a(j) 821 c_var_b(idim) = conn_info%phi_b(j) 822 c_var_c(idim) = conn_info%phi_c(j) 823 c_var_d(idim) = conn_info%phi_d(j) 824 IF (found) THEN 825 c_var_type(idim) = conn_info%phi_type(j) 826 END IF 827 END IF 828 END DO 829 END IF 830 CALL timestop(handle2) 831 CALL timeset(routineN//"_13", handle2) 832 ! map torsions on molecules 833 nvar1 = 0 834 nvar2 = 0 835 NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype) 836 IF (ASSOCIATED(conn_info%phi_a)) nvar1 = SIZE(conn_info%phi_a) 837 IF (ASSOCIATED(c_var_a)) nvar2 = SIZE(c_var_a) 838 nval_tot1 = nvar1 839 nval_tot2 = 0 840 ALLOCATE (map_var_mol(nvar1)) 841 ALLOCATE (map_cvar_mol(nvar2)) 842 map_var_mol = -1 843 map_cvar_mol = -1 844 DO i = 1, nvar1 845 j1 = map_atom_mol(conn_info%phi_a(i)) 846 j2 = map_atom_mol(conn_info%phi_b(i)) 847 j3 = map_atom_mol(conn_info%phi_c(i)) 848 j4 = map_atom_mol(conn_info%phi_d(i)) 849 IF (j1 == j2 .AND. j2 == j3 .AND. j3 == j4) THEN 850 IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%phi_a(i)) 851 END IF 852 END DO 853 DO i = 1, nvar2 854 min_index = MIN(c_var_a(i), c_var_b(i), c_var_c(i), c_var_d(i)) 855 j1 = map_atom_mol(min_index) 856 IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index) 857 END DO 858 CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1) 859 CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2) 860 DO i = 1, topology%nmol_type 861 intra_torsions = 0 862 inter_torsions = 0 863 IF (ALL(bnd_type(:, i) > 0)) THEN 864 intra_torsions = bnd_type(2, i) - bnd_type(1, i) + 1 865 END IF 866 IF (ALL(bnd_ctype(:, i) > 0)) THEN 867 inter_torsions = bnd_ctype(2, i) - bnd_ctype(1, i) + 1 868 END IF 869 itorsion = intra_torsions + inter_torsions 870 IF (iw > 0) THEN 871 WRITE (iw, *) " Total number of torsions for molecule type ", i, " :", itorsion 872 WRITE (iw, *) " intra (torsions inside molecules) :: ", intra_torsions 873 WRITE (iw, *) " inter (torsions between molecules) :: ", inter_torsions 874 END IF 875 molecule_kind => molecule_kind_set(i) 876 nval_tot2 = nval_tot2 + itorsion*SIZE(molecule_kind%molecule_list) 877 ALLOCATE (torsion_list(itorsion)) 878 itorsion = 0 879 DO j = bnd_type(1, i), bnd_type(2, i) 880 IF (j == 0) CYCLE 881 itorsion = itorsion + 1 882 jind = map_vars(j) 883 first = first_list(map_atom_mol(conn_info%phi_a(jind))) 884 torsion_list(itorsion)%a = conn_info%phi_a(jind) - first + 1 885 torsion_list(itorsion)%b = conn_info%phi_b(jind) - first + 1 886 torsion_list(itorsion)%c = conn_info%phi_c(jind) - first + 1 887 torsion_list(itorsion)%d = conn_info%phi_d(jind) - first + 1 888 ! Set by default id_type to charmm and modify when handling the forcefield 889 torsion_list(itorsion)%id_type = do_ff_charmm 890 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN 891 torsion_list(itorsion)%itype = conn_info%phi_type(jind) 892 END IF 893 !point this to the right torsion_kind_type if using force field 894 NULLIFY (torsion_list(itorsion)%torsion_kind) 895 IF (iw > 0) THEN 896 WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') & 897 "molecule_kind", i, "intra TOR", & 898 conn_info%phi_a(jind), & 899 conn_info%phi_b(jind), & 900 conn_info%phi_c(jind), & 901 conn_info%phi_d(jind), & 902 "offset number at", & 903 conn_info%phi_a(jind) - first + 1, & 904 conn_info%phi_b(jind) - first + 1, & 905 conn_info%phi_c(jind) - first + 1, & 906 conn_info%phi_d(jind) - first + 1 907 END IF 908 END DO 909 DO j = bnd_ctype(1, i), bnd_ctype(2, i) 910 IF (j == 0) CYCLE 911 itorsion = itorsion + 1 912 jind = map_cvars(j) 913 min_index = MIN(c_var_a(jind), c_var_b(jind), c_var_c(jind), c_var_d(jind)) 914 first = first_list(map_atom_mol(min_index)) 915 torsion_list(itorsion)%a = c_var_a(jind) - first + 1 916 torsion_list(itorsion)%b = c_var_b(jind) - first + 1 917 torsion_list(itorsion)%c = c_var_c(jind) - first + 1 918 torsion_list(itorsion)%d = c_var_d(jind) - first + 1 919 ! Set by default id_type to charmm and modify when handling the forcefield 920 torsion_list(itorsion)%id_type = do_ff_charmm 921 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN 922 torsion_list(itorsion)%itype = c_var_type(jind) 923 END IF 924 !point this to the right torsion_kind_type if using force field 925 NULLIFY (torsion_list(itorsion)%torsion_kind) 926 IF (iw > 0) THEN 927 WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') & 928 "molecule_kind", i, "inter TOR", & 929 c_var_a(jind), & 930 c_var_b(jind), & 931 c_var_c(jind), & 932 c_var_d(jind), & 933 "offset number at", & 934 c_var_a(jind) - first + 1, & 935 c_var_b(jind) - first + 1, & 936 c_var_c(jind) - first + 1, & 937 c_var_d(jind) - first + 1 938 END IF 939 END DO 940 CALL set_molecule_kind(molecule_kind=molecule_kind, & 941 ntorsion=SIZE(torsion_list), torsion_list=torsion_list) 942 END DO 943 CPASSERT(nval_tot1 == nval_tot2) 944 DEALLOCATE (map_var_mol) 945 DEALLOCATE (map_cvar_mol) 946 DEALLOCATE (map_vars) 947 DEALLOCATE (map_cvars) 948 DEALLOCATE (bnd_type) 949 DEALLOCATE (bnd_ctype) 950 DEALLOCATE (c_var_a) 951 DEALLOCATE (c_var_b) 952 DEALLOCATE (c_var_c) 953 DEALLOCATE (c_var_d) 954 IF (found) THEN 955 DEALLOCATE (c_var_type) 956 END IF 957 CALL timestop(handle2) 958 !----------------------------------------------------------------------------- 959 !----------------------------------------------------------------------------- 960 ! 14. Set the molecule_kind%[nimpr,impr_list] via set_molecule_kind 961 ! Also set the molecule_kind%[nopbend,opbend_list] 962 !----------------------------------------------------------------------------- 963 ! Allocate c_var_a, c_var_b, c_var_c, c_var_d, c_var_type 964 CALL timeset(routineN//"_14_pre", handle2) 965 idim = 0 966 ALLOCATE (c_var_a(idim)) 967 ALLOCATE (c_var_b(idim)) 968 ALLOCATE (c_var_c(idim)) 969 ALLOCATE (c_var_d(idim)) 970 found = ASSOCIATED(conn_info%impr_type) 971 IF (found) THEN 972 ALLOCATE (c_var_type(idim)) 973 END IF 974 IF (ASSOCIATED(conn_info%c_bond_a) .AND. ASSOCIATED(conn_info%impr_a)) THEN 975 DO j = 1, SIZE(conn_info%impr_a) 976 j1 = map_atom_mol(conn_info%impr_a(j)) 977 j2 = map_atom_mol(conn_info%impr_b(j)) 978 j3 = map_atom_mol(conn_info%impr_c(j)) 979 j4 = map_atom_mol(conn_info%impr_d(j)) 980 IF (j1 /= j2 .OR. j2 /= j3 .OR. j3 /= j4) THEN 981 idim = idim + 1 982 END IF 983 END DO 984 CALL reallocate(c_var_a, 1, idim) 985 CALL reallocate(c_var_b, 1, idim) 986 CALL reallocate(c_var_c, 1, idim) 987 CALL reallocate(c_var_d, 1, idim) 988 IF (found) THEN 989 CALL reallocate(c_var_type, 1, idim) 990 END IF 991 idim = 0 992 DO j = 1, SIZE(conn_info%impr_a) 993 j1 = map_atom_mol(conn_info%impr_a(j)) 994 j2 = map_atom_mol(conn_info%impr_b(j)) 995 j3 = map_atom_mol(conn_info%impr_c(j)) 996 j4 = map_atom_mol(conn_info%impr_d(j)) 997 IF (j1 /= j2 .OR. j2 /= j3 .OR. j3 /= j4) THEN 998 idim = idim + 1 999 c_var_a(idim) = conn_info%impr_a(j) 1000 c_var_b(idim) = conn_info%impr_b(j) 1001 c_var_c(idim) = conn_info%impr_c(j) 1002 c_var_d(idim) = conn_info%impr_d(j) 1003 IF (found) THEN 1004 c_var_type(idim) = conn_info%impr_type(j) 1005 END IF 1006 END IF 1007 END DO 1008 END IF 1009 CALL timestop(handle2) 1010 CALL timeset(routineN//"_14", handle2) 1011 ! map imprs on molecules 1012 nvar1 = 0 1013 nvar2 = 0 1014 NULLIFY (map_vars, map_cvars, bnd_type, bnd_ctype) 1015 IF (ASSOCIATED(conn_info%impr_a)) nvar1 = SIZE(conn_info%impr_a) 1016 IF (ASSOCIATED(c_var_a)) nvar2 = SIZE(c_var_a) 1017 nval_tot1 = nvar1 1018 nval_tot2 = 0 1019 ALLOCATE (map_var_mol(nvar1)) 1020 ALLOCATE (map_cvar_mol(nvar2)) 1021 map_var_mol = -1 1022 map_cvar_mol = -1 1023 DO i = 1, nvar1 1024 j1 = map_atom_mol(conn_info%impr_a(i)) 1025 j2 = map_atom_mol(conn_info%impr_b(i)) 1026 j3 = map_atom_mol(conn_info%impr_c(i)) 1027 j4 = map_atom_mol(conn_info%impr_d(i)) 1028 IF (j1 == j2 .AND. j2 == j3 .AND. j3 == j4) THEN 1029 IF (j1 > 0) map_var_mol(i) = map_atom_type(conn_info%impr_a(i)) 1030 END IF 1031 END DO 1032 DO i = 1, nvar2 1033 min_index = MIN(c_var_a(i), c_var_b(i), c_var_c(i), c_var_d(i)) 1034 j1 = map_atom_mol(min_index) 1035 IF (j1 > 0) map_cvar_mol(i) = map_atom_type(min_index) 1036 END DO 1037 CALL find_bnd_typ(topology%nmol_type, map_vars, map_var_mol, bnd_type, nvar1) 1038 CALL find_bnd_typ(topology%nmol_type, map_cvars, map_cvar_mol, bnd_ctype, nvar2) 1039 DO i = 1, topology%nmol_type 1040 intra_imprs = 0 1041 inter_imprs = 0 1042 IF (ALL(bnd_type(:, i) > 0)) THEN 1043 intra_imprs = bnd_type(2, i) - bnd_type(1, i) + 1 1044 END IF 1045 IF (ALL(bnd_ctype(:, i) > 0)) THEN 1046 inter_imprs = bnd_ctype(2, i) - bnd_ctype(1, i) + 1 1047 END IF 1048 iimpr = intra_imprs + inter_imprs 1049 IF (iw > 0) THEN 1050 WRITE (iw, *) " Total number of imprs for molecule type ", i, " :", iimpr 1051 WRITE (iw, *) " intra (imprs inside molecules) :: ", intra_imprs 1052 WRITE (iw, *) " inter (imprs between molecules) :: ", inter_imprs 1053 WRITE (iw, *) " Total number of opbends for molecule type ", i, " :", iimpr 1054 WRITE (iw, *) " intra (opbends inside molecules) :: ", intra_imprs 1055 WRITE (iw, *) " inter (opbends between molecules) :: ", inter_imprs 1056 END IF 1057 molecule_kind => molecule_kind_set(i) 1058 nval_tot2 = nval_tot2 + iimpr*SIZE(molecule_kind%molecule_list) 1059 ALLOCATE (impr_list(iimpr), STAT=stat) 1060 ALLOCATE (opbend_list(iimpr), STAT=stat) 1061 CPASSERT(stat == 0) 1062 iimpr = 0 1063 DO j = bnd_type(1, i), bnd_type(2, i) 1064 IF (j == 0) CYCLE 1065 iimpr = iimpr + 1 1066 jind = map_vars(j) 1067 first = first_list(map_atom_mol(conn_info%impr_a(jind))) 1068 impr_list(iimpr)%a = conn_info%impr_a(jind) - first + 1 1069 impr_list(iimpr)%b = conn_info%impr_b(jind) - first + 1 1070 impr_list(iimpr)%c = conn_info%impr_c(jind) - first + 1 1071 impr_list(iimpr)%d = conn_info%impr_d(jind) - first + 1 1072 ! Atom sequence for improper is A B C D in which A is central atom, 1073 ! B is deviating atom and C & D are secondairy atoms. Atom sequence for 1074 ! opbend is B D C A in which A is central atom, B is deviating. Hence 1075 ! to create an opbend out of an improper, B and D need to be interchanged. 1076 opbend_list(iimpr)%a = conn_info%impr_b(jind) - first + 1 1077 opbend_list(iimpr)%b = conn_info%impr_d(jind) - first + 1 1078 opbend_list(iimpr)%c = conn_info%impr_c(jind) - first + 1 1079 opbend_list(iimpr)%d = conn_info%impr_a(jind) - first + 1 1080 ! Set by default id_type of improper to charmm and modify when handling the forcefield 1081 impr_list(iimpr)%id_type = do_ff_charmm 1082 ! Set by default id_type of opbend to harmonic and modify when handling the forcefield 1083 opbend_list(iimpr)%id_type = do_ff_harmonic 1084 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN 1085 impr_list(iimpr)%itype = conn_info%impr_type(jind) 1086 END IF 1087 !point this to the right impr_kind_type if using force field 1088 NULLIFY (impr_list(iimpr)%impr_kind) 1089 NULLIFY (opbend_list(iimpr)%opbend_kind) 1090 IF (iw > 0) THEN 1091 WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') & 1092 "molecule_kind", i, "intra IMPR", & 1093 conn_info%impr_a(jind), & 1094 conn_info%impr_b(jind), & 1095 conn_info%impr_c(jind), & 1096 conn_info%impr_d(jind), & 1097 "offset number at", & 1098 conn_info%impr_a(jind) - first + 1, & 1099 conn_info%impr_b(jind) - first + 1, & 1100 conn_info%impr_c(jind) - first + 1, & 1101 conn_info%impr_d(jind) - first + 1 1102 WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') & 1103 "molecule_kind", i, "intra OPBEND", & 1104 conn_info%impr_b(jind), & 1105 conn_info%impr_d(jind), & 1106 conn_info%impr_c(jind), & 1107 conn_info%impr_a(jind), & 1108 "offset number at", & 1109 conn_info%impr_b(jind) - first + 1, & 1110 conn_info%impr_d(jind) - first + 1, & 1111 conn_info%impr_c(jind) - first + 1, & 1112 conn_info%impr_a(jind) - first + 1 1113 END IF 1114 END DO 1115 DO j = bnd_ctype(1, i), bnd_ctype(2, i) 1116 IF (j == 0) CYCLE 1117 iimpr = iimpr + 1 1118 jind = map_cvars(j) 1119 min_index = MIN(c_var_a(jind), c_var_b(jind), c_var_c(jind), c_var_d(jind)) 1120 first = first_list(map_atom_mol(min_index)) 1121 impr_list(iimpr)%a = c_var_a(jind) - first + 1 1122 impr_list(iimpr)%b = c_var_b(jind) - first + 1 1123 impr_list(iimpr)%c = c_var_c(jind) - first + 1 1124 impr_list(iimpr)%d = c_var_d(jind) - first + 1 1125 opbend_list(iimpr)%a = c_var_b(jind) - first + 1 1126 opbend_list(iimpr)%b = c_var_d(jind) - first + 1 1127 opbend_list(iimpr)%c = c_var_c(jind) - first + 1 1128 opbend_list(iimpr)%d = c_var_a(jind) - first + 1 1129 ! Set by default id_type of improper to charmm and modify when handling the forcefield 1130 impr_list(iimpr)%id_type = do_ff_charmm 1131 ! Set by default id_type of opbend to harmonic and modify when handling the forcefield 1132 opbend_list(iimpr)%id_type = do_ff_harmonic 1133 IF ((topology%conn_type == do_conn_g96) .OR. (topology%conn_type == do_conn_g87)) THEN 1134 impr_list(iimpr)%itype = c_var_type(jind) 1135 END IF 1136 !point this to the right impr_kind_type and opbend_kind_type if using force field 1137 NULLIFY (impr_list(iimpr)%impr_kind) 1138 NULLIFY (opbend_list(iimpr)%opbend_kind) 1139 IF (iw > 0) THEN 1140 WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') & 1141 "molecule_kind", i, "inter IMPR", & 1142 c_var_a(jind), & 1143 c_var_b(jind), & 1144 c_var_c(jind), & 1145 c_var_d(jind), & 1146 "offset number at", & 1147 c_var_a(jind) - first + 1, & 1148 c_var_b(jind) - first + 1, & 1149 c_var_c(jind) - first + 1, & 1150 c_var_d(jind) - first + 1 1151 WRITE (iw, '(7X,A,I3,1X,A,I4,I4,I4,I4,1X,A,I4,I4,I4,I4)') & 1152 "molecule_kind", i, "inter OPBEND", & 1153 c_var_b(jind), & 1154 c_var_d(jind), & 1155 c_var_c(jind), & 1156 c_var_a(jind), & 1157 "offset number at", & 1158 c_var_b(jind) - first + 1, & 1159 c_var_d(jind) - first + 1, & 1160 c_var_c(jind) - first + 1, & 1161 c_var_a(jind) - first + 1 1162 END IF 1163 END DO 1164 CALL set_molecule_kind(molecule_kind=molecule_kind, & 1165 nimpr=SIZE(impr_list), impr_list=impr_list) 1166 CALL set_molecule_kind(molecule_kind=molecule_kind, & 1167 nopbend=SIZE(opbend_list), opbend_list=opbend_list) 1168 END DO 1169 CPASSERT(nval_tot1 == nval_tot2) 1170 DEALLOCATE (map_var_mol) 1171 DEALLOCATE (map_cvar_mol) 1172 DEALLOCATE (map_vars) 1173 DEALLOCATE (map_cvars) 1174 DEALLOCATE (bnd_type) 1175 DEALLOCATE (bnd_ctype) 1176 DEALLOCATE (c_var_a) 1177 DEALLOCATE (c_var_b) 1178 DEALLOCATE (c_var_c) 1179 DEALLOCATE (c_var_d) 1180 IF (found) THEN 1181 DEALLOCATE (c_var_type) 1182 END IF 1183 CALL timestop(handle2) 1184 !----------------------------------------------------------------------------- 1185 !----------------------------------------------------------------------------- 1186 ! Finally deallocate some stuff. 1187 !----------------------------------------------------------------------------- 1188 DEALLOCATE (first_list) 1189 DEALLOCATE (last_list) 1190 DEALLOCATE (map_atom_mol) 1191 DEALLOCATE (map_atom_type) 1192 CALL timestop(handle) 1193 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 1194 "PRINT%TOPOLOGY_INFO/UTIL_INFO") 1195 END SUBROUTINE topology_connectivity_pack 1196 1197! ************************************************************************************************** 1198!> \brief used to achieve linear scaling in the connectivity_pack 1199!> \param nmol_type ... 1200!> \param map_vars ... 1201!> \param map_var_mol ... 1202!> \param bnd_type ... 1203!> \param nvar1 ... 1204!> \par History 1205!> Teodoro Laino 1206! ************************************************************************************************** 1207 SUBROUTINE find_bnd_typ(nmol_type, map_vars, map_var_mol, bnd_type, nvar1) 1208 INTEGER, INTENT(IN) :: nmol_type 1209 INTEGER, DIMENSION(:), POINTER :: map_vars, map_var_mol 1210 INTEGER, DIMENSION(:, :), POINTER :: bnd_type 1211 INTEGER, INTENT(IN) :: nvar1 1212 1213 CHARACTER(len=*), PARAMETER :: routineN = 'find_bnd_typ', routineP = moduleN//':'//routineN 1214 1215 INTEGER :: i, ibond, j 1216 1217 ALLOCATE (map_vars(nvar1)) 1218 CALL sort(map_var_mol, nvar1, map_vars) 1219 ALLOCATE (bnd_type(2, nmol_type)) 1220 bnd_type = 0 1221 IF (nvar1 == 0) RETURN 1222 DO j = 1, nvar1 1223 IF (map_var_mol(j) /= -1) EXIT 1224 END DO 1225 IF (j == nvar1 + 1) RETURN 1226 i = map_var_mol(j) 1227 bnd_type(1, i) = j 1228 DO ibond = j, nvar1 1229 IF (map_var_mol(ibond) /= i) THEN 1230 bnd_type(2, i) = ibond - 1 1231 i = map_var_mol(ibond) 1232 bnd_type(1, i) = ibond 1233 END IF 1234 END DO 1235 bnd_type(2, i) = nvar1 1236 1237 END SUBROUTINE find_bnd_typ 1238 1239! ************************************************************************************************** 1240!> \brief Handles the multiple unit cell option for the connectivity 1241!> \param topology ... 1242!> \param subsys_section ... 1243!> \author Teodoro Laino [tlaino] - 06.2009 1244! ************************************************************************************************** 1245 SUBROUTINE topology_conn_multiple(topology, subsys_section) 1246 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 1247 TYPE(section_vals_type), POINTER :: subsys_section 1248 1249 CHARACTER(len=*), PARAMETER :: routineN = 'topology_conn_multiple', & 1250 routineP = moduleN//':'//routineN 1251 1252 INTEGER :: a, fac, i, ind, j, k, m, natoms_orig, & 1253 nbond, nbond_c, nimpr, nonfo, nphi, & 1254 ntheta, nub 1255 INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell 1256 TYPE(connectivity_info_type), POINTER :: conn_info 1257 1258 NULLIFY (multiple_unit_cell) 1259 CALL section_vals_val_get(subsys_section, "TOPOLOGY%MULTIPLE_UNIT_CELL", & 1260 i_vals=multiple_unit_cell) 1261 IF (ANY(multiple_unit_cell /= 1)) THEN 1262 fac = PRODUCT(multiple_unit_cell) 1263 conn_info => topology%conn_info 1264 1265 nbond = 0 1266 IF (ASSOCIATED(conn_info%bond_a)) THEN 1267 nbond = SIZE(conn_info%bond_a) 1268 CALL reallocate(conn_info%bond_a, 1, nbond*fac) 1269 CALL reallocate(conn_info%bond_b, 1, nbond*fac) 1270 END IF 1271 1272 ntheta = 0 1273 IF (ASSOCIATED(conn_info%theta_a)) THEN 1274 ntheta = SIZE(conn_info%theta_a) 1275 CALL reallocate(conn_info%theta_a, 1, ntheta*fac) 1276 CALL reallocate(conn_info%theta_b, 1, ntheta*fac) 1277 CALL reallocate(conn_info%theta_c, 1, ntheta*fac) 1278 END IF 1279 1280 nphi = 0 1281 IF (ASSOCIATED(conn_info%phi_a)) THEN 1282 nphi = SIZE(conn_info%phi_a) 1283 CALL reallocate(conn_info%phi_a, 1, nphi*fac) 1284 CALL reallocate(conn_info%phi_b, 1, nphi*fac) 1285 CALL reallocate(conn_info%phi_c, 1, nphi*fac) 1286 CALL reallocate(conn_info%phi_d, 1, nphi*fac) 1287 END IF 1288 1289 nimpr = 0 1290 IF (ASSOCIATED(conn_info%impr_a)) THEN 1291 nimpr = SIZE(conn_info%impr_a) 1292 CALL reallocate(conn_info%impr_a, 1, nimpr*fac) 1293 CALL reallocate(conn_info%impr_b, 1, nimpr*fac) 1294 CALL reallocate(conn_info%impr_c, 1, nimpr*fac) 1295 CALL reallocate(conn_info%impr_d, 1, nimpr*fac) 1296 END IF 1297 1298 nbond_c = 0 1299 IF (ASSOCIATED(conn_info%c_bond_a)) THEN 1300 nbond_c = SIZE(conn_info%c_bond_a) 1301 CALL reallocate(conn_info%c_bond_a, 1, nbond_c*fac) 1302 CALL reallocate(conn_info%c_bond_b, 1, nbond_c*fac) 1303 END IF 1304 1305 nub = 0 1306 IF (ASSOCIATED(conn_info%ub_a)) THEN 1307 nub = SIZE(conn_info%ub_a) 1308 CALL reallocate(conn_info%ub_a, 1, nub*fac) 1309 CALL reallocate(conn_info%ub_b, 1, nub*fac) 1310 CALL reallocate(conn_info%ub_c, 1, nub*fac) 1311 END IF 1312 1313 nonfo = 0 1314 IF (ASSOCIATED(conn_info%onfo_a)) THEN 1315 nonfo = SIZE(conn_info%onfo_a) 1316 CALL reallocate(conn_info%onfo_a, 1, nonfo*fac) 1317 CALL reallocate(conn_info%onfo_b, 1, nonfo*fac) 1318 END IF 1319 1320 natoms_orig = topology%natoms/fac 1321 ind = 0 1322 DO k = 1, multiple_unit_cell(3) 1323 DO j = 1, multiple_unit_cell(2) 1324 DO i = 1, multiple_unit_cell(1) 1325 ind = ind + 1 1326 IF (ind == 1) CYCLE 1327 a = (ind - 1)*natoms_orig 1328 1329 ! Bonds 1330 IF (nbond > 0) THEN 1331 m = (ind - 1)*nbond 1332 conn_info%bond_a(m + 1:m + nbond) = conn_info%bond_a(1:nbond) + a 1333 conn_info%bond_b(m + 1:m + nbond) = conn_info%bond_b(1:nbond) + a 1334 END IF 1335 ! Theta 1336 IF (ntheta > 0) THEN 1337 m = (ind - 1)*ntheta 1338 conn_info%theta_a(m + 1:m + ntheta) = conn_info%theta_a(1:ntheta) + a 1339 conn_info%theta_b(m + 1:m + ntheta) = conn_info%theta_b(1:ntheta) + a 1340 conn_info%theta_c(m + 1:m + ntheta) = conn_info%theta_c(1:ntheta) + a 1341 END IF 1342 ! Phi 1343 IF (nphi > 0) THEN 1344 m = (ind - 1)*nphi 1345 conn_info%phi_a(m + 1:m + nphi) = conn_info%phi_a(1:nphi) + a 1346 conn_info%phi_b(m + 1:m + nphi) = conn_info%phi_b(1:nphi) + a 1347 conn_info%phi_c(m + 1:m + nphi) = conn_info%phi_c(1:nphi) + a 1348 conn_info%phi_d(m + 1:m + nphi) = conn_info%phi_d(1:nphi) + a 1349 END IF 1350 ! Impropers 1351 IF (nimpr > 0) THEN 1352 m = (ind - 1)*nimpr 1353 conn_info%impr_a(m + 1:m + nimpr) = conn_info%impr_a(1:nimpr) + a 1354 conn_info%impr_b(m + 1:m + nimpr) = conn_info%impr_b(1:nimpr) + a 1355 conn_info%impr_c(m + 1:m + nimpr) = conn_info%impr_c(1:nimpr) + a 1356 conn_info%impr_d(m + 1:m + nimpr) = conn_info%impr_d(1:nimpr) + a 1357 END IF 1358 ! Para_res 1359 IF (nbond_c > 0) THEN 1360 m = (ind - 1)*nbond_c 1361 conn_info%c_bond_a(m + 1:m + nbond_c) = conn_info%c_bond_a(1:nbond_c) + a 1362 conn_info%c_bond_b(m + 1:m + nbond_c) = conn_info%c_bond_b(1:nbond_c) + a 1363 END IF 1364 ! Urey-Bradley 1365 IF (nub > 0) THEN 1366 m = (ind - 1)*nub 1367 conn_info%ub_a(m + 1:m + nub) = conn_info%ub_a(1:nub) + a 1368 conn_info%ub_b(m + 1:m + nub) = conn_info%ub_b(1:nub) + a 1369 conn_info%ub_c(m + 1:m + nub) = conn_info%ub_c(1:nub) + a 1370 END IF 1371 ! ONFO 1372 IF (nonfo > 0) THEN 1373 m = (ind - 1)*nonfo 1374 conn_info%onfo_a(m + 1:m + nonfo) = conn_info%onfo_a(1:nonfo) + a 1375 conn_info%onfo_b(m + 1:m + nonfo) = conn_info%onfo_b(1:nonfo) + a 1376 END IF 1377 END DO 1378 END DO 1379 END DO 1380 END IF 1381 1382 END SUBROUTINE topology_conn_multiple 1383 1384END MODULE topology_connectivity_util 1385