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!> Teodor Laino 09.2006 - Major rewriting with linear scaling routines 10! ************************************************************************************************** 11MODULE topology_generate_util 12 USE atomic_kind_types, ONLY: atomic_kind_type,& 13 deallocate_atomic_kind_set,& 14 set_atomic_kind 15 USE cell_types, ONLY: pbc 16 USE cp_log_handling, ONLY: cp_get_default_logger,& 17 cp_logger_get_default_io_unit,& 18 cp_logger_type,& 19 cp_to_string 20 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 21 cp_print_key_unit_nr 22 USE cp_para_types, ONLY: cp_para_env_type 23 USE cp_units, ONLY: cp_unit_to_cp2k 24 USE fist_neighbor_list_types, ONLY: fist_neighbor_deallocate,& 25 fist_neighbor_type 26 USE fist_neighbor_lists, ONLY: build_fist_neighbor_lists 27 USE input_constants, ONLY: do_add,& 28 do_bondparm_covalent,& 29 do_bondparm_vdw,& 30 do_conn_off,& 31 do_conn_user,& 32 do_remove 33 USE input_section_types, ONLY: section_vals_get,& 34 section_vals_get_subs_vals,& 35 section_vals_type,& 36 section_vals_val_get 37 USE kinds, ONLY: default_string_length,& 38 dp 39 USE mathlib, ONLY: matvec_3x3 40 USE memory_utilities, ONLY: reallocate 41 USE particle_types, ONLY: allocate_particle_set,& 42 deallocate_particle_set,& 43 particle_type 44 USE periodic_table, ONLY: get_ptable_info 45 USE qmmm_types_low, ONLY: qmmm_env_mm_type 46 USE string_table, ONLY: id2str,& 47 s2s,& 48 str2id 49 USE string_utilities, ONLY: integer_to_string,& 50 uppercase 51 USE topology_types, ONLY: atom_info_type,& 52 connectivity_info_type,& 53 topology_parameters_type 54 USE topology_util, ONLY: array1_list_type,& 55 array2_list_type,& 56 find_molecule,& 57 give_back_molecule,& 58 reorder_list_array,& 59 reorder_structure 60 USE util, ONLY: find_boundary,& 61 sort 62#include "./base/base_uses.f90" 63 64 IMPLICIT NONE 65 66 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_generate_util' 67 68 PRIVATE 69 LOGICAL, PARAMETER :: debug_this_module = .FALSE. 70 71 PUBLIC :: topology_generate_bend, & 72 topology_generate_bond, & 73 topology_generate_dihe, & 74 topology_generate_impr, & 75 topology_generate_onfo, & 76 topology_generate_ub, & 77 topology_generate_molecule, & 78 topology_generate_molname 79 80CONTAINS 81 82! ************************************************************************************************** 83!> \brief Generates molnames: useful when the connectivity on file does not 84!> provide them 85!> \param conn_info ... 86!> \param natom ... 87!> \param natom_prev ... 88!> \param nbond_prev ... 89!> \param id_molname ... 90!> \author Teodoro Laino [tlaino] - University of Zurich 10.2008 91! ************************************************************************************************** 92 SUBROUTINE topology_generate_molname(conn_info, natom, natom_prev, nbond_prev, & 93 id_molname) 94 TYPE(connectivity_info_type), POINTER :: conn_info 95 INTEGER, INTENT(IN) :: natom, natom_prev, nbond_prev 96 INTEGER, DIMENSION(:), INTENT(INOUT) :: id_molname 97 98 CHARACTER(len=*), PARAMETER :: routineN = 'topology_generate_molname', & 99 routineP = moduleN//':'//routineN 100 CHARACTER(LEN=default_string_length), PARAMETER :: basename = "MOL" 101 102 CHARACTER(LEN=default_string_length) :: molname 103 INTEGER :: i, N, nmol 104 LOGICAL :: check 105 TYPE(array1_list_type), ALLOCATABLE, DIMENSION(:) :: atom_bond_list 106 107! convert a simple list of bonds to a list of bonds per atom 108! (each bond is present in the forward and backward direction 109 110 ALLOCATE (atom_bond_list(natom)) 111 DO I = 1, natom 112 ALLOCATE (atom_bond_list(I)%array1(0)) 113 ENDDO 114 N = 0 115 IF (ASSOCIATED(conn_info%bond_a)) N = SIZE(conn_info%bond_a) - nbond_prev 116 CALL reorder_structure(atom_bond_list, conn_info%bond_a(nbond_prev + 1:) - natom_prev, & 117 conn_info%bond_b(nbond_prev + 1:) - natom_prev, N) 118 119 nmol = 0 120 check = ALL(id_molname == str2id(s2s("__UNDEF__"))) .OR. ALL(id_molname /= str2id(s2s("__UNDEF__"))) 121 CPASSERT(check) 122 DO i = 1, natom 123 IF (id2str(id_molname(i)) == "__UNDEF__") THEN 124 molname = TRIM(basename)//ADJUSTL(cp_to_string(nmol)) 125 CALL generate_molname_low(i, atom_bond_list, molname, id_molname) 126 nmol = nmol + 1 127 END IF 128 END DO 129 DO I = 1, natom 130 DEALLOCATE (atom_bond_list(I)%array1) 131 ENDDO 132 DEALLOCATE (atom_bond_list) 133 134 END SUBROUTINE topology_generate_molname 135 136! ************************************************************************************************** 137!> \brief Generates molnames: useful when the connectivity on file does not 138!> provide them 139!> \param i ... 140!> \param atom_bond_list ... 141!> \param molname ... 142!> \param id_molname ... 143!> \author Teodoro Laino [tlaino] - University of Zurich 10.2008 144! ************************************************************************************************** 145 RECURSIVE SUBROUTINE generate_molname_low(i, atom_bond_list, molname, id_molname) 146 INTEGER, INTENT(IN) :: i 147 TYPE(array1_list_type), DIMENSION(:) :: atom_bond_list 148 CHARACTER(LEN=default_string_length), INTENT(IN) :: molname 149 INTEGER, DIMENSION(:), INTENT(INOUT) :: id_molname 150 151 CHARACTER(len=*), PARAMETER :: routineN = 'generate_molname_low', & 152 routineP = moduleN//':'//routineN 153 154 INTEGER :: j, k 155 156 IF (debug_this_module) THEN 157 WRITE (*, *) "Entered with :", i 158 WRITE (*, *) TRIM(molname)//": entering with i:", i, " full series to test:: ", atom_bond_list(i)%array1 159 IF ((TRIM(id2str(id_molname(i))) /= "__UNDEF__") .AND. & 160 (TRIM(id2str(id_molname(i))) /= TRIM(molname))) THEN 161 WRITE (*, *) "Atom (", i, ") has already a molecular name assigned ! ("//TRIM(id2str(id_molname(i)))//")." 162 WRITE (*, *) "New molecular name would be: ("//TRIM(molname)//")." 163 WRITE (*, *) "Detecting something wrong in the molecular setup!" 164 CPABORT("") 165 END IF 166 END IF 167 id_molname(i) = str2id(molname) 168 DO j = 1, SIZE(atom_bond_list(i)%array1) 169 k = atom_bond_list(i)%array1(j) 170 IF (debug_this_module) WRITE (*, *) "entering with i:", i, "testing :", k 171 IF (k == -1) CYCLE 172 atom_bond_list(i)%array1(j) = -1 173 WHERE (atom_bond_list(k)%array1 == i) atom_bond_list(k)%array1 = -1 174 CALL generate_molname_low(k, atom_bond_list, molname, id_molname) 175 END DO 176 END SUBROUTINE generate_molname_low 177 178! ************************************************************************************************** 179!> \brief Use information from bond list to generate molecule. (ie clustering) 180!> \param topology ... 181!> \param qmmm ... 182!> \param qmmm_env ... 183!> \param subsys_section ... 184! ************************************************************************************************** 185 SUBROUTINE topology_generate_molecule(topology, qmmm, qmmm_env, subsys_section) 186 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 187 LOGICAL, INTENT(in), OPTIONAL :: qmmm 188 TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env 189 TYPE(section_vals_type), POINTER :: subsys_section 190 191 CHARACTER(len=*), PARAMETER :: routineN = 'topology_generate_molecule', & 192 routineP = moduleN//':'//routineN 193 INTEGER, PARAMETER :: nblock = 100 194 195 INTEGER :: atom_in_kind, atom_in_mol, first, handle, handle2, i, iatm, iatom, iend, ifirst, & 196 ilast, inum, istart, itype, iw, j, jump1, jump2, last, max_mol_num, mol_num, mol_res, & 197 mol_typ, myind, N, natom, nlocl, ntype, resid 198 INTEGER, DIMENSION(:), POINTER :: qm_atom_index, wrk1, wrk2 199 LOGICAL :: do_again, found, my_qmmm 200 TYPE(array1_list_type), ALLOCATABLE, DIMENSION(:) :: atom_bond_list 201 TYPE(atom_info_type), POINTER :: atom_info 202 TYPE(connectivity_info_type), POINTER :: conn_info 203 TYPE(cp_logger_type), POINTER :: logger 204 205 NULLIFY (logger) 206 logger => cp_get_default_logger() 207 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", & 208 extension=".subsysLog") 209 CALL timeset(routineN, handle) 210 NULLIFY (qm_atom_index) 211 NULLIFY (wrk1) 212 NULLIFY (wrk2) 213 214 atom_info => topology%atom_info 215 conn_info => topology%conn_info 216 ! 217 ! QM/MM coordinate_control 218 ! 219 my_qmmm = .FALSE. 220 IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm 221 222 natom = topology%natoms 223 IF (ASSOCIATED(atom_info%map_mol_typ)) DEALLOCATE (atom_info%map_mol_typ) 224 ALLOCATE (atom_info%map_mol_typ(natom)) 225 226 IF (ASSOCIATED(atom_info%map_mol_num)) DEALLOCATE (atom_info%map_mol_num) 227 ALLOCATE (atom_info%map_mol_num(natom)) 228 229 IF (ASSOCIATED(atom_info%map_mol_res)) DEALLOCATE (atom_info%map_mol_res) 230 ALLOCATE (atom_info%map_mol_res(natom)) 231 232 ! Initialisation 233 atom_info%map_mol_typ(:) = 0 234 atom_info%map_mol_num(:) = -1 235 atom_info%map_mol_res(:) = 1 236 237 ! Parse the atom list to find the different molecule types and residues 238 ntype = 1 239 atom_info%map_mol_typ(1) = 1 240 resid = 1 241 CALL reallocate(wrk1, 1, nblock) 242 wrk1(1) = atom_info%id_molname(1) 243 DO iatom = 2, natom 244 IF (topology%conn_type == do_conn_off) THEN 245 ! No connectivity: each atom becomes a molecule of its own molecule kind 246 ntype = ntype + 1 247 atom_info%map_mol_typ(iatom) = ntype 248 ELSE IF (topology%conn_type == do_conn_user) THEN 249 ! User-defined connectivity: 5th column of COORD section or molecule 250 ! or residue name in the case of PDB files 251 IF (atom_info%id_molname(iatom) == atom_info%id_molname(iatom - 1)) THEN 252 atom_info%map_mol_typ(iatom) = atom_info%map_mol_typ(iatom - 1) 253 IF (atom_info%id_resname(iatom) == atom_info%id_resname(iatom - 1)) THEN 254 atom_info%map_mol_res(iatom) = atom_info%map_mol_res(iatom - 1) 255 ELSE 256 resid = resid + 1 257 atom_info%map_mol_res(iatom) = resid 258 END IF 259 ELSE 260 ! Check if the type is already known 261 found = .FALSE. 262 DO itype = 1, ntype 263 IF (atom_info%id_molname(iatom) == wrk1(itype)) THEN 264 atom_info%map_mol_typ(iatom) = itype 265 found = .TRUE. 266 EXIT 267 END IF 268 END DO 269 IF (.NOT. found) THEN 270 ntype = ntype + 1 271 atom_info%map_mol_typ(iatom) = ntype 272 IF (ntype > SIZE(wrk1)) CALL reallocate(wrk1, 1, 2*SIZE(wrk1)) 273 wrk1(ntype) = atom_info%id_molname(iatom) 274 END IF 275 resid = resid + 1 276 atom_info%map_mol_res(iatom) = resid 277 END IF 278 ELSE 279 IF (atom_info%id_molname(iatom - 1) == atom_info%id_molname(iatom)) THEN 280 atom_info%map_mol_typ(iatom) = ntype 281 ELSE 282 ntype = ntype + 1 283 atom_info%map_mol_typ(iatom) = ntype 284 END IF 285 END IF 286 END DO 287 DEALLOCATE (wrk1) 288 289 IF (iw > 0) WRITE (iw, '(/,T2,A)') "Start of molecule generation" 290 291 ! convert a simple list of bonds to a list of bonds per atom 292 ! (each bond is present in the forward and backward direction 293 ALLOCATE (atom_bond_list(natom)) 294 DO I = 1, natom 295 ALLOCATE (atom_bond_list(I)%array1(0)) 296 ENDDO 297 N = 0 298 IF (ASSOCIATED(conn_info%bond_a)) N = SIZE(conn_info%bond_a) 299 CALL reorder_structure(atom_bond_list, conn_info%bond_a, conn_info%bond_b, N) 300 CALL find_molecule(atom_bond_list, atom_info%map_mol_num, atom_info%id_molname) 301 DO I = 1, natom 302 DEALLOCATE (atom_bond_list(I)%array1) 303 ENDDO 304 DEALLOCATE (atom_bond_list) 305 IF (iw > 0) WRITE (iw, '(/,T2,A)') "End of molecule generation" 306 307 ! Modify according map_mol_typ the array map_mol_num 308 IF (iw > 0) WRITE (iw, '(/,T2,A)') "Checking for non-continuous generated molecules" 309 ! Check molecule number 310 ALLOCATE (wrk1(natom)) 311 ALLOCATE (wrk2(natom)) 312 wrk1 = atom_info%map_mol_num 313 314 IF (debug_this_module) THEN 315 DO i = 1, natom 316 WRITE (*, '(2I10)') i, atom_info%map_mol_num(i) 317 END DO 318 END IF 319 320 CALL sort(wrk1, natom, wrk2) 321 istart = 1 322 mol_typ = wrk1(istart) 323 DO i = 2, natom 324 IF (mol_typ /= wrk1(i)) THEN 325 iend = i - 1 326 first = MINVAL(wrk2(istart:iend)) 327 last = MAXVAL(wrk2(istart:iend)) 328 nlocl = last - first + 1 329 IF (iend - istart + 1 /= nlocl) THEN 330 IF (debug_this_module) WRITE (*, *) iend, istart, iend - istart + 1, first, last, nlocl 331 CALL cp_abort(__LOCATION__, & 332 "CP2K requires molecules to be contiguous and we have detected a non contiguous one!! "// & 333 "In particular a molecule defined from index ("//cp_to_string(first)//") to ("// & 334 cp_to_string(last)//") contains other molecules, not connected! "// & 335 "Too late at this stage everything should be already ordered! "// & 336 "If you have not yet employed the REORDERING keyword, please do so."// & 337 "It may help to fix this issue.") 338 END IF 339 istart = i 340 mol_typ = wrk1(istart) 341 END IF 342 END DO 343 iend = i - 1 344 first = MINVAL(wrk2(istart:iend)) 345 last = MAXVAL(wrk2(istart:iend)) 346 nlocl = last - first + 1 347 IF (iend - istart + 1 /= nlocl) THEN 348 IF (debug_this_module) WRITE (*, *) iend, istart, iend - istart + 1, first, last, nlocl 349 CALL cp_abort(__LOCATION__, & 350 "CP2K requires molecules to be contiguous and we have detected a non contiguous one!! "// & 351 "In particular a molecule defined from index ("//cp_to_string(first)//") to ("// & 352 cp_to_string(last)//") contains other molecules, not connected! "// & 353 "Too late at this stage everything should be already ordered! "// & 354 "If you have not yet employed the REORDERING keyword, please do so."// & 355 "It may help to fix this issue.") 356 END IF 357 DEALLOCATE (wrk1) 358 DEALLOCATE (wrk2) 359 IF (iw > 0) WRITE (iw, '(/,T2,A)') "End of check" 360 361 IF (iw > 0) WRITE (UNIT=iw, FMT="(/,T2,A)") "Start of renumbering molecules" 362 IF (topology%conn_type == do_conn_user) THEN 363 mol_num = 1 364 atom_info%map_mol_num(1) = 1 365 DO iatom = 2, natom 366 IF (atom_info%id_molname(iatom) /= atom_info%id_molname(iatom - 1)) THEN 367 mol_num = 1 368 ELSE IF (atom_info%map_mol_res(iatom) /= atom_info%map_mol_res(iatom - 1)) THEN 369 mol_num = mol_num + 1 370 END IF 371 atom_info%map_mol_num(iatom) = mol_num 372 END DO 373 ELSE 374 mol_typ = atom_info%map_mol_typ(1) 375 mol_num = atom_info%map_mol_num(1) 376 DO i = 2, natom 377 IF (atom_info%map_mol_typ(i) /= mol_typ) THEN 378 myind = atom_info%map_mol_num(i) - mol_num + 1 379 CPASSERT(myind /= atom_info%map_mol_num(i - 1)) 380 mol_typ = atom_info%map_mol_typ(i) 381 mol_num = atom_info%map_mol_num(i) 382 END IF 383 atom_info%map_mol_num(i) = atom_info%map_mol_num(i) - mol_num + 1 384 END DO 385 END IF 386 IF (iw > 0) WRITE (UNIT=iw, FMT="(/,T2,A)") "End of renumbering molecules" 387 388 ! Optionally, use the residues as molecules 389 CALL timeset(routineN//"_PARA_RES", handle2) 390 IF (iw > 0) WRITE (UNIT=iw, FMT="(/,T2,A,L2)") "Starting PARA_RES: ", topology%para_res 391 IF (topology%para_res) THEN 392 IF (topology%conn_type == do_conn_user) THEN 393 atom_info%id_molname(:) = atom_info%id_resname(:) 394 ntype = 1 395 atom_info%map_mol_typ(1) = 1 396 mol_num = 1 397 atom_info%map_mol_num(1) = 1 398 DO iatom = 2, natom 399 IF (atom_info%id_molname(iatom) /= atom_info%id_molname(iatom - 1)) THEN 400 ntype = ntype + 1 401 mol_num = 1 402 ELSE IF (atom_info%map_mol_res(iatom) /= atom_info%map_mol_res(iatom - 1)) THEN 403 mol_num = mol_num + 1 404 END IF 405 atom_info%map_mol_typ(iatom) = ntype 406 atom_info%map_mol_num(iatom) = mol_num 407 END DO 408 ELSE 409 mol_res = 1 410 mol_typ = atom_info%map_mol_typ(1) 411 mol_num = atom_info%map_mol_num(1) 412 atom_info%map_mol_res(1) = mol_res 413 DO i = 2, natom 414 IF ((atom_info%resid(i - 1) /= atom_info%resid(i)) .OR. & 415 (atom_info%id_resname(i - 1) /= atom_info%id_resname(i))) THEN 416 mol_res = mol_res + 1 417 END IF 418 IF ((atom_info%map_mol_typ(i) /= mol_typ) .OR. & 419 (atom_info%map_mol_num(i) /= mol_num)) THEN 420 mol_typ = atom_info%map_mol_typ(i) 421 mol_num = atom_info%map_mol_num(i) 422 mol_res = 1 423 END IF 424 atom_info%map_mol_res(i) = mol_res 425 END DO 426 END IF 427 END IF 428 IF (iw > 0) WRITE (UNIT=iw, FMT="(/,T2,A)") "End of PARA_RES" 429 CALL timestop(handle2) 430 431 IF (iw > 0) THEN 432 DO iatom = 1, natom 433 WRITE (iw, '(4(1X,A,":",I0),2(1X,A,1X,A))') "iatom", iatom, & 434 "map_mol_typ", atom_info%map_mol_typ(iatom), & 435 "map_mol_num", atom_info%map_mol_num(iatom), & 436 "map_mol_res", atom_info%map_mol_res(iatom), & 437 "mol_name:", TRIM(id2str(atom_info%id_molname(iatom))), & 438 "res_name:", TRIM(id2str(atom_info%id_resname(iatom))) 439 END DO 440 END IF 441 442 IF (my_qmmm) THEN 443 do_again = .FALSE. 444 IF (iw > 0) WRITE (iw, *) "MAP_MOL_NUM ", atom_info%map_mol_num 445 IF (iw > 0) WRITE (iw, *) "MAP_MOL_TYP ", atom_info%map_mol_typ 446 IF (iw > 0) WRITE (iw, *) "MAP_MOL_RES ", atom_info%map_mol_res 447 ALLOCATE (qm_atom_index(SIZE(qmmm_env%qm_atom_index))) 448 qm_atom_index = qmmm_env%qm_atom_index 449 CPASSERT(ALL(qm_atom_index /= 0)) 450 DO myind = 1, SIZE(qm_atom_index) 451 IF (qm_atom_index(myind) == 0) CYCLE 452 CALL find_boundary(atom_info%map_mol_typ, natom, ifirst, ilast, & 453 atom_info%map_mol_typ(qm_atom_index(myind))) 454 CALL find_boundary(atom_info%map_mol_typ, atom_info%map_mol_num, natom, ifirst, ilast, & 455 atom_info%map_mol_typ(qm_atom_index(myind)), atom_info%map_mol_num(qm_atom_index(myind))) 456 IF (iw > 0) WRITE (iw, *) "qm fragment:: ifirst, ilast", ifirst, ilast 457 CPASSERT(((ifirst /= 0) .OR. (ilast /= natom))) 458 DO iatm = ifirst, ilast 459 atom_info%id_molname(iatm) = str2id(s2s("_QM_"// & 460 TRIM(id2str(atom_info%id_molname(iatm))))) 461 IF (iw > 0) WRITE (iw, *) "QM Molecule name :: ", id2str(atom_info%id_molname(iatm)) 462 WHERE (qm_atom_index == iatm) qm_atom_index = 0 463 END DO 464 DO iatm = 1, ifirst - 1 465 IF (ANY(qm_atom_index == iatm)) do_again = .TRUE. 466 END DO 467 DO iatm = ilast + 1, natom 468 IF (ANY(qm_atom_index == iatm)) do_again = .TRUE. 469 END DO 470 IF (iw > 0) WRITE (iw, *) " Another QM fragment? :: ", do_again 471 IF (ifirst /= 1) THEN 472 jump1 = atom_info%map_mol_typ(ifirst) - atom_info%map_mol_typ(ifirst - 1) 473 CPASSERT(jump1 <= 1 .AND. jump1 >= 0) 474 jump1 = ABS(jump1 - 1) 475 ELSE 476 jump1 = 0 477 END IF 478 IF (ilast /= natom) THEN 479 jump2 = atom_info%map_mol_typ(ilast + 1) - atom_info%map_mol_typ(ilast) 480 CPASSERT(jump2 <= 1 .AND. jump2 >= 0) 481 jump2 = ABS(jump2 - 1) 482 ELSE 483 jump2 = 0 484 END IF 485 486 ! Changing mol_type consistently 487 DO iatm = ifirst, natom 488 atom_info%map_mol_typ(iatm) = atom_info%map_mol_typ(iatm) + jump1 489 END DO 490 DO iatm = ilast + 1, natom 491 atom_info%map_mol_typ(iatm) = atom_info%map_mol_typ(iatm) + jump2 492 END DO 493 IF (jump1 == 1) THEN 494 DO iatm = ifirst, ilast 495 atom_info%map_mol_num(iatm) = 1 496 END DO 497 END IF 498 499 IF (jump2 == 1) THEN 500 CALL find_boundary(atom_info%map_mol_typ, natom, first, last, atom_info%map_mol_typ(ilast + 1)) 501 CALL find_boundary(atom_info%map_mol_typ, atom_info%map_mol_num, natom, ifirst, ilast, & 502 atom_info%map_mol_typ(ilast + 1), atom_info%map_mol_num(ilast + 1)) 503 atom_in_mol = ilast - ifirst + 1 504 inum = 1 505 DO iatm = first, last, atom_in_mol 506 atom_info%map_mol_num(iatm:iatm + atom_in_mol - 1) = inum 507 inum = inum + 1 508 END DO 509 END IF 510 511 IF (.NOT. do_again) EXIT 512 END DO 513 DEALLOCATE (qm_atom_index) 514 515 IF (iw > 0) THEN 516 WRITE (iw, *) "After the QM/MM Setup:" 517 DO iatom = 1, natom 518 WRITE (iw, *) " iatom,map_mol_typ,map_mol_num ", iatom, & 519 atom_info%map_mol_typ(iatom), atom_info%map_mol_num(iatom) 520 END DO 521 END IF 522 END IF 523 ! 524 ! Further check : see if the number of atoms belonging to same molecule kinds 525 ! are equal 526 IF (iw > 0) THEN 527 WRITE (iw, *) "SUMMARY:: Number of molecule kinds found:", ntype 528 ntype = MAXVAL(atom_info%map_mol_typ) 529 DO i = 1, ntype 530 atom_in_kind = COUNT(atom_info%map_mol_typ == i) 531 WRITE (iw, *) "Molecule kind:", i, " contains", atom_in_kind, " atoms" 532 IF (atom_in_kind <= 1) CYCLE 533 CALL find_boundary(atom_info%map_mol_typ, natom, first, last, i) 534 WRITE (iw, *) "Boundary atoms:", first, last 535 CPASSERT(last - first + 1 == atom_in_kind) 536 max_mol_num = MAXVAL(atom_info%map_mol_num(first:last)) 537 WRITE (iw, *) "Number of molecules of kind", i, "is ::", max_mol_num 538 atom_in_mol = atom_in_kind/max_mol_num 539 WRITE (iw, *) "Number of atoms per each molecule:", atom_in_mol 540 WRITE (iw, *) "MAP_MOL_TYP::", atom_info%map_mol_typ(first:last) 541 WRITE (iw, *) "MAP_MOL_NUM::", atom_info%map_mol_num(first:last) 542 WRITE (iw, *) "MAP_MOL_RES::", atom_info%map_mol_res(first:last) 543 ! 544 DO j = 1, max_mol_num 545 IF (COUNT(atom_info%map_mol_num(first:last) == j) /= atom_in_mol) THEN 546 WRITE (iw, *) "molecule type:", i, "molecule num:", j, " has ", & 547 COUNT(atom_info%map_mol_num(first:last) == j), & 548 " atoms instead of ", atom_in_mol, " ." 549 CALL cp_abort(__LOCATION__, & 550 "Two molecules of the same kind have "// & 551 "been created with different numbers of atoms!") 552 END IF 553 END DO 554 END DO 555 END IF 556 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 557 "PRINT%TOPOLOGY_INFO/UTIL_INFO") 558 CALL timestop(handle) 559 END SUBROUTINE topology_generate_molecule 560 561! ************************************************************************************************** 562!> \brief Use info from periodic table and assumptions to generate bonds 563!> \param topology ... 564!> \param para_env ... 565!> \param subsys_section ... 566!> \author Teodoro Laino 09.2006 567! ************************************************************************************************** 568 SUBROUTINE topology_generate_bond(topology, para_env, subsys_section) 569 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 570 TYPE(cp_para_env_type), POINTER :: para_env 571 TYPE(section_vals_type), POINTER :: subsys_section 572 573 CHARACTER(len=*), PARAMETER :: routineN = 'topology_generate_bond', & 574 routineP = moduleN//':'//routineN 575 576 CHARACTER(LEN=2) :: upper_sym_1 577 INTEGER :: cbond, handle, handle2, i, iatm1, iatm2, iatom, ibond, idim, iw, j, jatom, k, & 578 n_bonds, n_heavy_bonds, n_hydr_bonds, n_rep, natom, npairs, output_unit 579 INTEGER, ALLOCATABLE, DIMENSION(:) :: bond_a, bond_b, list, map_nb 580 INTEGER, DIMENSION(:), POINTER :: isolated_atoms, tmp_v 581 LOGICAL :: connectivity_ok, explicit, print_info 582 LOGICAL, ALLOCATABLE, DIMENSION(:) :: h_list 583 REAL(KIND=dp) :: bondparm_factor, cell_v(3), dr(3), & 584 ksign, my_maxrad, r2, r2_min, rbond, & 585 rbond2, tmp 586 REAL(KIND=dp), DIMENSION(1, 1) :: r_max, r_minsq 587 REAL(KIND=dp), DIMENSION(:), POINTER :: radius 588 REAL(KIND=dp), DIMENSION(:, :), POINTER :: pbc_coord 589 TYPE(array2_list_type), DIMENSION(:), POINTER :: bond_list 590 TYPE(atom_info_type), POINTER :: atom_info 591 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 592 TYPE(atomic_kind_type), POINTER :: atomic_kind 593 TYPE(connectivity_info_type), POINTER :: conn_info 594 TYPE(cp_logger_type), POINTER :: logger 595 TYPE(fist_neighbor_type), POINTER :: nonbonded 596 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 597 TYPE(section_vals_type), POINTER :: bond_section, generate_section, & 598 isolated_section 599 600 NULLIFY (logger, particle_set, atomic_kind_set, nonbonded, bond_section, generate_section) 601 NULLIFY (isolated_atoms, tmp_v) 602 CALL timeset(routineN, handle) 603 logger => cp_get_default_logger() 604 output_unit = cp_logger_get_default_io_unit(logger) 605 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", & 606 extension=".subsysLog") 607 ! Get atoms that one considers isolated (like ions in solution) 608 ALLOCATE (isolated_atoms(0)) 609 generate_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE") 610 isolated_section => section_vals_get_subs_vals(generate_section, "ISOLATED_ATOMS") 611 CALL section_vals_get(isolated_section, explicit=explicit) 612 IF (explicit) THEN 613 CALL section_vals_val_get(isolated_section, "LIST", n_rep_val=n_rep) 614 DO i = 1, n_rep 615 CALL section_vals_val_get(isolated_section, "LIST", i_vals=tmp_v, i_rep_val=i) 616 CALL reallocate(isolated_atoms, 1, SIZE(isolated_atoms) + SIZE(tmp_v)) 617 isolated_atoms(SIZE(isolated_atoms) - SIZE(tmp_v) + 1:SIZE(isolated_atoms)) = tmp_v 618 END DO 619 END IF 620 atom_info => topology%atom_info 621 conn_info => topology%conn_info 622 bondparm_factor = topology%bondparm_factor 623 cbond = 0 624 natom = topology%natoms 625 NULLIFY (radius) 626 ! Allocate temporary arrays 627 ALLOCATE (radius(natom)) 628 ALLOCATE (list(natom)) 629 ALLOCATE (h_list(natom)) 630 ALLOCATE (pbc_coord(3, natom)) 631 h_list = .FALSE. 632 CALL timeset(TRIM(routineN)//"_1", handle2) 633 DO iatom = 1, natom 634 list(iatom) = iatom 635 upper_sym_1 = id2str(atom_info%id_element(iatom)) 636 IF (topology%bondparm_type == do_bondparm_covalent) THEN 637 CALL get_ptable_info(symbol=upper_sym_1, covalent_radius=radius(iatom)) 638 ELSE IF (topology%bondparm_type == do_bondparm_vdw) THEN 639 CALL get_ptable_info(symbol=upper_sym_1, vdw_radius=radius(iatom)) 640 ELSE 641 CPABORT("Illegal bondparm_type") 642 END IF 643 IF (upper_sym_1 == "H ") h_list(iatom) = .TRUE. 644 ! isolated atoms? put the radius to 0.0_dp 645 IF (ANY(isolated_atoms == iatom)) radius(iatom) = 0.0_dp 646 radius(iatom) = cp_unit_to_cp2k(radius(iatom), "angstrom") 647 IF (iw > 0) WRITE (iw, '(T2,"GENERATE|",5X,A,T50,A5,T60,A,T69,F12.6)') & 648 "In topology_generate_bond :: iatom = ", upper_sym_1, & 649 "radius:", radius(iatom) 650 END DO 651 CALL timestop(handle2) 652 CALL timeset(TRIM(routineN)//"_2", handle2) 653 ! Initialize fake particle_set and atomic_kinds to generate the bond list 654 ! using the neighboring list routine 655 ALLOCATE (atomic_kind_set(1)) 656 CALL allocate_particle_set(particle_set, natom) 657 ! 658 my_maxrad = MAXVAL(radius)*2.0_dp 659 atomic_kind => atomic_kind_set(1) 660 CALL set_atomic_kind(atomic_kind=atomic_kind, kind_number=1, & 661 name="XXX", element_symbol="XXX", mass=0.0_dp, atom_list=list) 662 CALL section_vals_val_get(subsys_section, "TOPOLOGY%GENERATE%BONDLENGTH_MAX", r_val=tmp) 663 r_max = tmp 664 IF (my_maxrad*bondparm_factor > r_max(1, 1) .AND. (.NOT. topology%molname_generated)) THEN 665 IF (output_unit > 0) THEN 666 WRITE (output_unit, '(T2,"GENERATE|",A)') & 667 " ERROR in connectivity generation!", & 668 " The THRESHOLD to select possible bonds is larger than the max. bondlength", & 669 " used to build the neighbors lists. Increase the BONDLENGTH_MAX parameter" 670 WRITE (output_unit, '(T2,"GENERATE|",2(A,F11.6),A)') & 671 " Present THRESHOLD (", my_maxrad*bondparm_factor, " ). "// & 672 " Present BONDLENGTH_MAX (", r_max(1, 1), " )" 673 END IF 674 CPABORT("") 675 END IF 676 DO i = 1, natom 677 particle_set(i)%atomic_kind => atomic_kind_set(1) 678 particle_set(i)%r(1) = atom_info%r(1, i) 679 particle_set(i)%r(2) = atom_info%r(2, i) 680 particle_set(i)%r(3) = atom_info%r(3, i) 681 pbc_coord(:, i) = pbc(atom_info%r(:, i), topology%cell) 682 END DO 683 CALL section_vals_val_get(subsys_section, "TOPOLOGY%GENERATE%BONDLENGTH_MIN", r_val=tmp) 684 r_minsq = tmp*tmp 685 CALL timestop(handle2) 686 CALL timeset(TRIM(routineN)//"_3", handle2) 687 CALL build_fist_neighbor_lists(atomic_kind_set, particle_set, & 688 cell=topology%cell, r_max=r_max, r_minsq=r_minsq, & 689 ei_scale14=1.0_dp, vdw_scale14=1.0_dp, nonbonded=nonbonded, & 690 para_env=para_env, build_from_scratch=.TRUE., geo_check=.TRUE., & 691 mm_section=generate_section) 692 IF (iw > 0) THEN 693 WRITE (iw, '(T2,"GENERATE| Number of prescreened bonds (neighbors):",T71,I10)') & 694 nonbonded%neighbor_kind_pairs(1)%npairs 695 END IF 696 npairs = 0 697 DO i = 1, SIZE(nonbonded%neighbor_kind_pairs) 698 npairs = npairs + nonbonded%neighbor_kind_pairs(i)%npairs 699 END DO 700 ALLOCATE (bond_a(npairs)) 701 ALLOCATE (bond_b(npairs)) 702 ALLOCATE (map_nb(npairs)) 703 idim = 0 704 DO j = 1, SIZE(nonbonded%neighbor_kind_pairs) 705 DO i = 1, nonbonded%neighbor_kind_pairs(j)%npairs 706 idim = idim + 1 707 bond_a(idim) = nonbonded%neighbor_kind_pairs(j)%list(1, i) 708 bond_b(idim) = nonbonded%neighbor_kind_pairs(j)%list(2, i) 709 map_nb(idim) = j 710 END DO 711 END DO 712 CALL timestop(handle2) 713 CALL timeset(TRIM(routineN)//"_4", handle2) 714 ! We have a list of neighbors let's order the list w.r.t. the particle number 715 ALLOCATE (bond_list(natom)) 716 DO I = 1, natom 717 ALLOCATE (bond_list(I)%array1(0)) 718 ALLOCATE (bond_list(I)%array2(0)) 719 ENDDO 720 CALL reorder_structure(bond_list, bond_a, bond_b, map_nb, SIZE(bond_a)) 721 DEALLOCATE (bond_a) 722 DEALLOCATE (bond_b) 723 DEALLOCATE (map_nb) 724 ! Find the Real bonds in the system 725 ! Let's start with heavy atoms.. hydrogens will be treated only later on... 726 ! Heavy atoms loop 727 CALL reallocate(conn_info%bond_a, 1, 1) 728 CALL reallocate(conn_info%bond_b, 1, 1) 729 connectivity_ok = .FALSE. 730 ! No need to check consistency between provided molecule name and 731 ! generated connectivity since we overrided the molecule definition. 732 IF (topology%create_molecules) THEN 733 atom_info%id_molname = str2id(s2s("TO_DEFINE_LATER")) 734 ! A real name assignment will then be performed in the reorder module.. 735 END IF 736 ! It may happen that the connectivity created is fault for the missing 737 ! of one bond.. this external loop ensures that everything was created 738 ! fits exactly with the definition of molecules.. 739 DO WHILE (.NOT. connectivity_ok) 740 n_heavy_bonds = 0 741 n_bonds = 0 742 DO iatm1 = 1, natom 743 IF (h_list(iatm1)) CYCLE 744 DO j = 1, SIZE(bond_list(iatm1)%array1) 745 iatm2 = bond_list(iatm1)%array1(j) 746 IF (atom_info%id_molname(iatm1) /= atom_info%id_molname(iatm2)) CYCLE 747 IF (h_list(iatm2) .OR. (iatm2 <= iatm1)) CYCLE 748 k = bond_list(iatm1)%array2(j) 749 ksign = SIGN(1.0_dp, REAL(k, KIND=dp)) 750 k = ABS(k) 751 CALL matvec_3x3(cell_v, topology%cell%hmat, & 752 REAL(nonbonded%neighbor_kind_pairs(k)%cell_vector, KIND=dp)) 753 dr = pbc_coord(:, iatm1) - pbc_coord(:, iatm2) - ksign*cell_v 754 r2 = DOT_PRODUCT(dr, dr) 755 IF (r2 <= r_minsq(1, 1)) THEN 756 CALL cp_abort(__LOCATION__, & 757 "bond distance between atoms less then the smallest distance provided "// & 758 "in input "//cp_to_string(tmp)//" [bohr]") 759 END IF 760 ! Screen isolated atoms 761 IF ((ANY(isolated_atoms == iatm1)) .OR. (ANY(isolated_atoms == iatm2))) CYCLE 762 763 ! Screen neighbors 764 IF (topology%bondparm_type == do_bondparm_covalent) THEN 765 rbond = radius(iatm1) + radius(iatm2) 766 ELSE IF (topology%bondparm_type == do_bondparm_vdw) THEN 767 rbond = MAX(radius(iatm1), radius(iatm2)) 768 END IF 769 rbond2 = rbond*rbond 770 rbond2 = rbond2*(bondparm_factor)**2 771 !Test the distance to the sum of the covalent radius 772 IF (r2 <= rbond2) THEN 773 n_heavy_bonds = n_heavy_bonds + 1 774 CALL add_bonds_list(conn_info, iatm1, iatm2, n_heavy_bonds) 775 END IF 776 END DO 777 END DO 778 n_hydr_bonds = 0 779 n_bonds = n_heavy_bonds 780 ! Now check bonds formed by hydrogens... 781 ! The hydrogen valence is 1 so we can choose the closest atom.. 782 IF (output_unit > 0) WRITE (output_unit, *) 783 DO iatm1 = 1, natom 784 IF (.NOT. h_list(iatm1)) CYCLE 785 r2_min = HUGE(0.0_dp) 786 ibond = -1 787 print_info = .TRUE. 788 DO j = 1, SIZE(bond_list(iatm1)%array1) 789 iatm2 = bond_list(iatm1)%array1(j) 790 print_info = .FALSE. 791 IF (atom_info%id_molname(iatm1) /= atom_info%id_molname(iatm2)) CYCLE 792 IF (h_list(iatm2) .AND. (iatm2 <= iatm1)) CYCLE 793 ! Screen isolated atoms 794 IF ((ANY(isolated_atoms == iatm1)) .OR. (ANY(isolated_atoms == iatm2))) CYCLE 795 796 k = bond_list(iatm1)%array2(j) 797 ksign = SIGN(1.0_dp, REAL(k, KIND=dp)) 798 k = ABS(k) 799 CALL matvec_3x3(cell_v, topology%cell%hmat, & 800 REAL(nonbonded%neighbor_kind_pairs(k)%cell_vector, KIND=dp)) 801 dr = pbc_coord(:, iatm1) - pbc_coord(:, iatm2) - ksign*cell_v 802 r2 = DOT_PRODUCT(dr, dr) 803 IF (r2 <= r_minsq(1, 1)) THEN 804 CALL cp_abort(__LOCATION__, & 805 "bond distance between atoms less then the smallest distance provided "// & 806 "in input "//cp_to_string(tmp)//" [bohr]") 807 END IF 808 IF (r2 <= r2_min) THEN 809 r2_min = r2 810 ibond = iatm2 811 END IF 812 END DO 813 IF (ibond == -1) THEN 814 IF (output_unit > 0 .AND. print_info) THEN 815 WRITE (output_unit, '(T2,"GENERATE|",1X,A,I10,A)') & 816 "WARNING:: No connections detected for Hydrogen - Atom Nr:", iatm1, " !" 817 END IF 818 ELSE 819 n_hydr_bonds = n_hydr_bonds + 1 820 n_bonds = n_bonds + 1 821 CALL add_bonds_list(conn_info, MIN(iatm1, ibond), MAX(iatm1, ibond), n_bonds) 822 END IF 823 END DO 824 IF (output_unit > 0) THEN 825 WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') & 826 " Preliminary Number of Bonds generated:", n_bonds 827 END IF 828 ! External defined bonds (useful for complex connectivity) 829 bond_section => section_vals_get_subs_vals(generate_section, "BOND") 830 CALL connectivity_external_control(section=bond_section, & 831 Iarray1=conn_info%bond_a, & 832 Iarray2=conn_info%bond_b, & 833 nvar=n_bonds, & 834 topology=topology, & 835 output_unit=output_unit) 836 ! Resize arrays to their proper size.. 837 CALL reallocate(conn_info%bond_a, 1, n_bonds) 838 CALL reallocate(conn_info%bond_b, 1, n_bonds) 839 IF (topology%create_molecules) THEN 840 ! Since we create molecule names we're not sure that all atoms are contiguous 841 ! so we need to reorder them on the basis of the generated name 842 IF (.NOT. topology%reorder_atom) THEN 843 topology%reorder_atom = .TRUE. 844 IF (output_unit > 0) WRITE (output_unit, '(T2,"GENERATE|",A)') & 845 " Molecules names have been generated. Now reordering particle set in order to have ", & 846 " atoms belonging to the same molecule in a sequential order." 847 END IF 848 connectivity_ok = .TRUE. 849 ELSE 850 ! Check created connectivity and possibly give the OK to proceed 851 connectivity_ok = check_generate_mol(conn_info%bond_a, conn_info%bond_b, & 852 atom_info, bondparm_factor, output_unit) 853 END IF 854 IF (my_maxrad*bondparm_factor > r_max(1, 1) .AND. (.NOT. topology%molname_generated)) THEN 855 IF (output_unit > 0) THEN 856 WRITE (output_unit, '(T2,"GENERATE|",A)') & 857 " ERROR in connectivity generation!", & 858 " The THRESHOLD to select possible bonds is bigger than the MAX bondlength", & 859 " used to build the neighbors lists. Increase the BONDLENGTH_MAX patameter" 860 WRITE (output_unit, '(T2,"GENERATE|",2(A,F11.6),A)') & 861 " Present THRESHOLD (", my_maxrad*bondparm_factor, " ). "// & 862 " Present BONDLENGTH_MAX (", r_max(1, 1), " )" 863 END IF 864 CPABORT("") 865 END IF 866 END DO 867 IF (connectivity_ok .AND. (output_unit > 0)) THEN 868 WRITE (output_unit, '(T2,"GENERATE|",A)') & 869 " Achieved consistency in connectivity generation." 870 END IF 871 CALL fist_neighbor_deallocate(nonbonded) 872 CALL timestop(handle2) 873 CALL timeset(TRIM(routineN)//"_6", handle2) 874 ! Deallocate temporary working arrays 875 DO I = 1, natom 876 DEALLOCATE (bond_list(I)%array1) 877 DEALLOCATE (bond_list(I)%array2) 878 ENDDO 879 DEALLOCATE (bond_list) 880 DEALLOCATE (pbc_coord) 881 DEALLOCATE (radius) 882 DEALLOCATE (list) 883 CALL deallocate_particle_set(particle_set) 884 CALL deallocate_atomic_kind_set(atomic_kind_set) 885 ! 886 CALL timestop(handle2) 887 IF (output_unit > 0 .AND. n_bonds > 0) THEN 888 WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of Bonds generated:", & 889 n_bonds 890 END IF 891 CALL timeset(TRIM(routineN)//"_7", handle2) 892 ! If PARA_RES then activate RESIDUES 893 CALL reallocate(conn_info%c_bond_a, 1, 0) 894 CALL reallocate(conn_info%c_bond_b, 1, 0) 895 IF (topology%para_res) THEN 896 DO ibond = 1, SIZE(conn_info%bond_a) 897 iatom = conn_info%bond_a(ibond) 898 jatom = conn_info%bond_b(ibond) 899 IF ((atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) .OR. & 900 (atom_info%resid(iatom) /= atom_info%resid(jatom)) .OR. & 901 (atom_info%id_resname(iatom) /= atom_info%id_resname(jatom))) THEN 902 IF (iw > 0) WRITE (iw, *) " PARA_RES, bond between molecules atom ", & 903 iatom, jatom 904 cbond = cbond + 1 905 CALL reallocate(conn_info%c_bond_a, 1, cbond) 906 CALL reallocate(conn_info%c_bond_b, 1, cbond) 907 conn_info%c_bond_a(cbond) = iatom 908 conn_info%c_bond_b(cbond) = jatom 909 ELSE 910 IF (atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) THEN 911 CPABORT("Bonds between different molecule types?") 912 END IF 913 END IF 914 END DO 915 END IF 916 CALL timestop(handle2) 917 DEALLOCATE (isolated_atoms) 918 CALL timestop(handle) 919 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 920 "PRINT%TOPOLOGY_INFO/GENERATE_INFO") 921 END SUBROUTINE topology_generate_bond 922 923! ************************************************************************************************** 924!> \brief Performs a check on the generated connectivity 925!> \param bond_a ... 926!> \param bond_b ... 927!> \param atom_info ... 928!> \param bondparm_factor ... 929!> \param output_unit ... 930!> \return ... 931!> \author Teodoro Laino 09.2006 932! ************************************************************************************************** 933 FUNCTION check_generate_mol(bond_a, bond_b, atom_info, bondparm_factor, output_unit) & 934 RESULT(conn_ok) 935 INTEGER, DIMENSION(:), POINTER :: bond_a, bond_b 936 TYPE(atom_info_type), POINTER :: atom_info 937 REAL(KIND=dp), INTENT(INOUT) :: bondparm_factor 938 INTEGER, INTENT(IN) :: output_unit 939 LOGICAL :: conn_ok 940 941 CHARACTER(len=*), PARAMETER :: routineN = 'check_generate_mol', & 942 routineP = moduleN//':'//routineN 943 944 CHARACTER(LEN=10) :: ctmp1, ctmp2, ctmp3 945 INTEGER :: handle, i, idim, itype, j, mol_natom, & 946 natom, nsize 947 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: mol_info_tmp 948 INTEGER, DIMENSION(:), POINTER :: mol_map, mol_map_o, wrk 949 INTEGER, DIMENSION(:, :), POINTER :: mol_info 950 LOGICAL, DIMENSION(:), POINTER :: icheck 951 TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list 952 953 CALL timeset(routineN, handle) 954 conn_ok = .TRUE. 955 natom = SIZE(atom_info%id_atmname) 956 ALLOCATE (bond_list(natom)) 957 DO I = 1, natom 958 ALLOCATE (bond_list(I)%array1(0)) 959 ENDDO 960 CALL reorder_structure(bond_list, bond_a, bond_b, SIZE(bond_a)) 961 ALLOCATE (mol_map(natom)) 962 ALLOCATE (mol_map_o(natom)) 963 ALLOCATE (wrk(natom)) 964 965 DO i = 1, natom 966 mol_map(i) = atom_info%id_molname(i) 967 END DO 968 mol_map_o = mol_map 969 970 CALL sort(mol_map, natom, wrk) 971 ! 972 ! mol(i,1) : stores id of the molecule 973 ! mol(i,2) : stores the total number of atoms forming that kind of molecule 974 ! mol(i,3) : contains the number of molecules generated for that kind 975 ! mol(i,4) : contains the number of atoms forming one molecule of that kind 976 ! Connectivity will be considered correct only if for each i: 977 ! 978 ! mol(i,2) = mol(i,3)*mol(i,4) 979 ! 980 ! If not, very probably, a bond is missing increase bondparm by 10% and let's 981 ! check if the newest connectivity is bug free.. 982 ! 983 984 ALLOCATE (mol_info_tmp(natom, 2)) 985 986 itype = mol_map(1) 987 nsize = 1 988 idim = 1 989 mol_info_tmp(1, 1) = itype 990 DO i = 2, natom 991 IF (mol_map(i) /= itype) THEN 992 nsize = nsize + 1 993 itype = mol_map(i) 994 mol_info_tmp(nsize, 1) = itype 995 mol_info_tmp(nsize - 1, 2) = idim 996 idim = 1 997 ELSE 998 idim = idim + 1 999 END IF 1000 END DO 1001 mol_info_tmp(nsize, 2) = idim 1002 1003 ALLOCATE (mol_info(nsize, 4)) 1004 mol_info(1:nsize, 1:2) = mol_info_tmp(1:nsize, 1:2) 1005 DEALLOCATE (mol_info_tmp) 1006 1007 DO i = 1, nsize 1008 mol_info(i, 3) = 0 1009 mol_info(i, 4) = 0 1010 END DO 1011 ! 1012 ALLOCATE (icheck(natom)) 1013 icheck = .FALSE. 1014 DO i = 1, natom 1015 IF (icheck(i)) CYCLE 1016 itype = mol_map_o(i) 1017 mol_natom = 0 1018 CALL give_back_molecule(icheck, bond_list, i, mol_natom, mol_map_o, mol_map_o(i)) 1019 DO j = 1, SIZE(mol_info) 1020 IF (itype == mol_info(j, 1)) EXIT 1021 END DO 1022 mol_info(j, 3) = mol_info(j, 3) + 1 1023 IF (mol_info(j, 4) == 0) mol_info(j, 4) = mol_natom 1024 IF (mol_info(j, 4) /= mol_natom) THEN 1025 ! Two same molecules have been found with different number 1026 ! of atoms. This usually indicates a missing bond in the 1027 ! generated connectivity 1028 ! Set connectivity to .false. EXIT and increase bondparm_factor by 1.05 1029 conn_ok = .FALSE. 1030 bondparm_factor = bondparm_factor*1.05_dp 1031 IF (output_unit < 0) EXIT 1032 WRITE (output_unit, '(/,T2,"GENERATE|",A)') " WARNING in connectivity generation!" 1033 WRITE (output_unit, '(T2,"GENERATE|",A)') & 1034 ' Two molecules/residues named ('//TRIM(id2str(itype))//') have different '// & 1035 ' number of atoms.' 1036 CALL integer_to_string(i, ctmp1) 1037 CALL integer_to_string(mol_natom, ctmp2) 1038 CALL integer_to_string(mol_info(j, 4), ctmp3) 1039 WRITE (output_unit, '(T2,"GENERATE|",A)') ' Molecule starting at position ('// & 1040 TRIM(ctmp1)//') has Nr. <'//TRIM(ctmp2)// & 1041 '> of atoms.', ' while the other same molecules have Nr. <'// & 1042 TRIM(ctmp3)//'> of atoms!' 1043 WRITE (output_unit, '(T2,"GENERATE|",A)') & 1044 ' Increasing bondparm_factor by 1.05.. An error was found in the generated', & 1045 ' connectivity. Retry...' 1046 WRITE (output_unit, '(T2,"GENERATE|",A,F11.6,A,/)') & 1047 " Present value of BONDPARM_FACTOR (", bondparm_factor, " )." 1048 EXIT 1049 END IF 1050 END DO 1051 1052 DEALLOCATE (icheck) 1053 DEALLOCATE (mol_info) 1054 DEALLOCATE (mol_map) 1055 DEALLOCATE (mol_map_o) 1056 DEALLOCATE (wrk) 1057 DO I = 1, natom 1058 DEALLOCATE (bond_list(I)%array1) 1059 ENDDO 1060 DEALLOCATE (bond_list) 1061 CALL timestop(handle) 1062 END FUNCTION check_generate_mol 1063 1064! ************************************************************************************************** 1065!> \brief Add/Remove a bond to the generated list 1066!> Particularly useful for system with complex connectivity 1067!> \param section ... 1068!> \param Iarray1 ... 1069!> \param Iarray2 ... 1070!> \param Iarray3 ... 1071!> \param Iarray4 ... 1072!> \param nvar ... 1073!> \param topology ... 1074!> \param output_unit ... 1075!> \param is_impr ... 1076!> \author Teodoro Laino 09.2006 1077! ************************************************************************************************** 1078 SUBROUTINE connectivity_external_control(section, Iarray1, Iarray2, Iarray3, Iarray4, nvar, & 1079 topology, output_unit, is_impr) 1080 TYPE(section_vals_type), POINTER :: section 1081 INTEGER, DIMENSION(:), POINTER :: Iarray1, Iarray2 1082 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: Iarray3, Iarray4 1083 INTEGER, INTENT(INOUT) :: nvar 1084 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 1085 INTEGER, INTENT(IN) :: output_unit 1086 LOGICAL, INTENT(IN), OPTIONAL :: is_impr 1087 1088 CHARACTER(len=*), PARAMETER :: routineN = 'connectivity_external_control', & 1089 routineP = moduleN//':'//routineN 1090 1091 CHARACTER(LEN=8) :: fmt 1092 INTEGER :: do_action, do_it, i, j, k, n_rep, & 1093 n_rep_val, natom, new_size, nsize 1094 INTEGER, DIMENSION(:), POINTER :: atlist, Ilist1, Ilist2, Ilist3, Ilist4 1095 LOGICAL :: explicit, ip3, ip4 1096 1097 natom = topology%natoms 1098 ! Preliminary sort of arrays 1099 ip3 = PRESENT(Iarray3) 1100 ip4 = PRESENT(Iarray4) 1101 nsize = 2 1102 IF (ip3) nsize = nsize + 1 1103 IF (ip3 .AND. ip4) nsize = nsize + 1 1104 ! Put the lists always in the canonical order 1105 CALL reorder_list_array(Iarray1, Iarray2, Iarray3, Iarray4, nsize, nvar) 1106 ! Go on with external control 1107 CALL section_vals_get(section, explicit=explicit, n_repetition=n_rep) 1108 IF (explicit) THEN 1109 NULLIFY (Ilist1, Ilist2, Ilist3, Ilist4, atlist) 1110 ALLOCATE (Ilist1(nvar)) 1111 ALLOCATE (Ilist2(nvar)) 1112 Ilist1 = Iarray1(1:nvar) 1113 Ilist2 = Iarray2(1:nvar) 1114 SELECT CASE (nsize) 1115 CASE (2) !do nothing 1116 CASE (3) 1117 ALLOCATE (Ilist3(nvar)) 1118 Ilist3 = Iarray3(1:nvar) 1119 CASE (4) 1120 ALLOCATE (Ilist3(nvar)) 1121 ALLOCATE (Ilist4(nvar)) 1122 Ilist3 = Iarray3(1:nvar) 1123 Ilist4 = Iarray4(1:nvar) 1124 CASE DEFAULT 1125 ! Should never reach this point 1126 CPABORT("") 1127 END SELECT 1128 CALL list_canonical_order(Ilist1, Ilist2, Ilist3, Ilist4, nsize, is_impr) 1129 ! 1130 DO i = 1, n_rep 1131 CALL section_vals_val_get(section, "ATOMS", i_rep_section=i, n_rep_val=n_rep_val) 1132 CALL section_vals_val_get(section, "_SECTION_PARAMETERS_", i_rep_section=i, & 1133 i_val=do_action) 1134 ! 1135 DO j = 1, n_rep_val 1136 CALL section_vals_val_get(section, "ATOMS", i_rep_section=i, i_rep_val=j, & 1137 i_vals=atlist) 1138 CPASSERT(SIZE(atlist) == nsize) 1139 CALL integer_to_string(nsize - 1, fmt) 1140 CALL check_element_list(do_it, do_action, atlist, Ilist1, Ilist2, Ilist3, Ilist4, & 1141 is_impr) 1142 IF (do_action == do_add) THEN 1143 ! Add to the element to the list 1144 IF (do_it > 0) THEN 1145 nvar = nvar + 1 1146 IF (output_unit > 0) THEN 1147 WRITE (output_unit, '(T2,"ADD|",1X,A,I6,'//TRIM(fmt)//'(A,I6),A,T64,A,I6)') & 1148 "element (", & 1149 atlist(1), (",", atlist(k), k=2, nsize), ") added.", " NEW size::", nvar 1150 END IF 1151 IF (nvar > SIZE(Iarray1)) THEN 1152 new_size = INT(5 + 1.2*nvar) 1153 CALL reallocate(Iarray1, 1, new_size) 1154 CALL reallocate(Iarray2, 1, new_size) 1155 SELECT CASE (nsize) 1156 CASE (3) 1157 CALL reallocate(Iarray3, 1, new_size) 1158 CASE (4) 1159 CALL reallocate(Iarray3, 1, new_size) 1160 CALL reallocate(Iarray4, 1, new_size) 1161 END SELECT 1162 END IF 1163 ! Using Ilist instead of atlist the canonical order is preserved.. 1164 Iarray1(do_it + 1:nvar) = Iarray1(do_it:nvar - 1) 1165 Iarray2(do_it + 1:nvar) = Iarray2(do_it:nvar - 1) 1166 Iarray1(do_it) = Ilist1(do_it) 1167 Iarray2(do_it) = Ilist2(do_it) 1168 SELECT CASE (nsize) 1169 CASE (3) 1170 Iarray3(do_it + 1:nvar) = Iarray3(do_it:nvar - 1) 1171 Iarray3(do_it) = Ilist3(do_it) 1172 CASE (4) 1173 Iarray3(do_it + 1:nvar) = Iarray3(do_it:nvar - 1) 1174 Iarray4(do_it + 1:nvar) = Iarray4(do_it:nvar - 1) 1175 Iarray3(do_it) = Ilist3(do_it) 1176 Iarray4(do_it) = Ilist4(do_it) 1177 END SELECT 1178 ELSE 1179 IF (output_unit > 0) THEN 1180 WRITE (output_unit, '(T2,"ADD|",1X,A,I6,'//TRIM(fmt)//'(A,I6),A,T80,A)') & 1181 "element (", & 1182 atlist(1), (",", atlist(k), k=2, nsize), ") already found.", "X" 1183 END IF 1184 END IF 1185 ELSE 1186 ! Remove element from the list 1187 IF (do_it > 0) THEN 1188 nvar = nvar - 1 1189 IF (output_unit > 0) THEN 1190 WRITE (output_unit, '(T2,"RMV|",1X,A,I6,'//TRIM(fmt)//'(A,I6),A,T64,A,I6)') & 1191 "element (", & 1192 atlist(1), (",", atlist(k), k=2, nsize), ") removed.", " NEW size::", nvar 1193 END IF 1194 Iarray1(do_it:nvar) = Iarray1(do_it + 1:nvar + 1) 1195 Iarray2(do_it:nvar) = Iarray2(do_it + 1:nvar + 1) 1196 Iarray1(nvar + 1) = -HUGE(0) 1197 Iarray2(nvar + 1) = -HUGE(0) 1198 SELECT CASE (nsize) 1199 CASE (3) 1200 Iarray3(do_it:nvar) = Iarray3(do_it + 1:nvar + 1) 1201 Iarray3(nvar + 1) = -HUGE(0) 1202 CASE (4) 1203 Iarray3(do_it:nvar) = Iarray3(do_it + 1:nvar + 1) 1204 Iarray4(do_it:nvar) = Iarray4(do_it + 1:nvar + 1) 1205 Iarray3(nvar + 1) = -HUGE(0) 1206 Iarray4(nvar + 1) = -HUGE(0) 1207 END SELECT 1208 ELSE 1209 IF (output_unit > 0) THEN 1210 WRITE (output_unit, '(T2,"RMV|",1X,A,I6,'//TRIM(fmt)//'(A,I6),A,T80,A)') & 1211 "element (", & 1212 atlist(1), (",", atlist(k), k=2, nsize), ") not found.", "X" 1213 END IF 1214 END IF 1215 END IF 1216 1217 END DO 1218 END DO 1219 DEALLOCATE (Ilist1) 1220 DEALLOCATE (Ilist2) 1221 SELECT CASE (nsize) 1222 CASE (2) ! do nothing 1223 CASE (3) 1224 DEALLOCATE (Ilist3) 1225 CASE (4) 1226 DEALLOCATE (Ilist3) 1227 DEALLOCATE (Ilist4) 1228 CASE DEFAULT 1229 ! Should never reach this point 1230 CPABORT("") 1231 END SELECT 1232 END IF 1233 END SUBROUTINE connectivity_external_control 1234 1235! ************************************************************************************************** 1236!> \brief Orders list in the canonical order: the extrema of the list are such 1237!> that the first extrema is always smaller or equal to the last extrema. 1238!> \param Ilist1 ... 1239!> \param Ilist2 ... 1240!> \param Ilist3 ... 1241!> \param Ilist4 ... 1242!> \param nsize ... 1243!> \param is_impr ... 1244!> \author Teodoro Laino 09.2006 1245! ************************************************************************************************** 1246 SUBROUTINE list_canonical_order(Ilist1, Ilist2, Ilist3, Ilist4, nsize, is_impr) 1247 INTEGER, DIMENSION(:), POINTER :: Ilist1, Ilist2 1248 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: Ilist3, Ilist4 1249 INTEGER, INTENT(IN) :: nsize 1250 LOGICAL, INTENT(IN), OPTIONAL :: is_impr 1251 1252 INTEGER :: i, ss(3), tmp1, tmp2, tmp3, tt(3) 1253 LOGICAL :: do_impr 1254 1255 do_impr = .FALSE. 1256 IF (PRESENT(is_impr)) do_impr = is_impr 1257 SELECT CASE (nsize) 1258 CASE (2) 1259 DO i = 1, SIZE(Ilist1) 1260 tmp1 = Ilist1(i) 1261 tmp2 = Ilist2(i) 1262 Ilist1(i) = MIN(tmp1, tmp2) 1263 Ilist2(i) = MAX(tmp1, tmp2) 1264 END DO 1265 CASE (3) 1266 DO i = 1, SIZE(Ilist1) 1267 tmp1 = Ilist1(i) 1268 tmp2 = Ilist3(i) 1269 Ilist1(i) = MIN(tmp1, tmp2) 1270 Ilist3(i) = MAX(tmp1, tmp2) 1271 END DO 1272 CASE (4) 1273 DO i = 1, SIZE(Ilist1) 1274 IF (.NOT. do_impr) THEN 1275 tmp1 = Ilist1(i) 1276 tmp2 = Ilist4(i) 1277 Ilist1(i) = MIN(tmp1, tmp2) 1278 IF (Ilist1(i) == tmp2) THEN 1279 tmp3 = Ilist3(i) 1280 Ilist3(i) = Ilist2(i) 1281 Ilist2(i) = tmp3 1282 END IF 1283 Ilist4(i) = MAX(tmp1, tmp2) 1284 ELSE 1285 tt(1) = Ilist2(i) 1286 tt(2) = Ilist3(i) 1287 tt(3) = Ilist4(i) 1288 CALL sort(tt, 3, ss) 1289 Ilist2(i) = tt(1) 1290 Ilist3(i) = tt(2) 1291 Ilist4(i) = tt(3) 1292 END IF 1293 END DO 1294 END SELECT 1295 1296 END SUBROUTINE list_canonical_order 1297 1298! ************************************************************************************************** 1299!> \brief finds an element in the ordered list 1300!> \param do_it ... 1301!> \param do_action ... 1302!> \param atlist ... 1303!> \param Ilist1 ... 1304!> \param Ilist2 ... 1305!> \param Ilist3 ... 1306!> \param Ilist4 ... 1307!> \param is_impr ... 1308!> \author Teodoro Laino 09.2006 1309! ************************************************************************************************** 1310 SUBROUTINE check_element_list(do_it, do_action, atlist, Ilist1, Ilist2, Ilist3, Ilist4, & 1311 is_impr) 1312 INTEGER, INTENT(OUT) :: do_it 1313 INTEGER, INTENT(IN) :: do_action 1314 INTEGER, DIMENSION(:), POINTER :: atlist, Ilist1, Ilist2 1315 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: Ilist3, Ilist4 1316 LOGICAL, INTENT(IN), OPTIONAL :: is_impr 1317 1318 INTEGER :: i, iend, istart, ndim, new_size, nsize, & 1319 ss(3), tmp1, tmp2, tmp3, tt(3) 1320 INTEGER, DIMENSION(4) :: tmp 1321 LOGICAL :: do_impr, found 1322 1323 do_impr = .FALSE. 1324 IF (PRESENT(is_impr)) do_impr = is_impr 1325 found = .FALSE. 1326 nsize = SIZE(atlist) 1327 ndim = SIZE(Ilist1) 1328 DO i = 1, nsize 1329 tmp(i) = atlist(i) 1330 END DO 1331 SELECT CASE (nsize) 1332 CASE (2) 1333 tmp1 = tmp(1) 1334 tmp2 = tmp(2) 1335 tmp(1) = MIN(tmp1, tmp2) 1336 tmp(2) = MAX(tmp1, tmp2) 1337 CASE (3) 1338 tmp1 = tmp(1) 1339 tmp2 = tmp(3) 1340 tmp(1) = MIN(tmp1, tmp2) 1341 tmp(3) = MAX(tmp1, tmp2) 1342 CASE (4) 1343 IF (.NOT. do_impr) THEN 1344 tmp1 = tmp(1) 1345 tmp2 = tmp(4) 1346 tmp(1) = MIN(tmp1, tmp2) 1347 IF (tmp(1) == tmp2) THEN 1348 tmp3 = tmp(3) 1349 tmp(3) = tmp(2) 1350 tmp(2) = tmp3 1351 END IF 1352 tmp(4) = MAX(tmp1, tmp2) 1353 ELSE 1354 tt(1) = tmp(2) 1355 tt(2) = tmp(3) 1356 tt(3) = tmp(4) 1357 CALL sort(tt, 3, ss) 1358 tmp(2) = tt(1) 1359 tmp(3) = tt(2) 1360 tmp(4) = tt(3) 1361 END IF 1362 END SELECT 1363 ! boundary to search 1364 DO istart = 1, ndim 1365 IF (Ilist1(istart) >= tmp(1)) EXIT 1366 END DO 1367 ! if nothing there stay within bounds 1368 IF (istart <= ndim) THEN 1369 IF (Ilist1(istart) > tmp(1) .AND. (istart /= 1)) istart = istart - 1 1370 ENDIF 1371 DO iend = istart, ndim 1372 IF (Ilist1(iend) /= tmp(1)) EXIT 1373 END DO 1374 IF (iend == ndim + 1) iend = ndim 1375 ! Final search in array 1376 SELECT CASE (nsize) 1377 CASE (2) 1378 DO i = istart, iend 1379 IF ((Ilist1(i) > tmp(1)) .OR. (Ilist2(i) > tmp(2))) EXIT 1380 IF ((Ilist1(i) == tmp(1)) .AND. (Ilist2(i) == tmp(2))) THEN 1381 found = .TRUE. 1382 EXIT 1383 END IF 1384 END DO 1385 CASE (3) 1386 DO i = istart, iend 1387 IF ((Ilist1(i) > tmp(1)) .OR. (Ilist2(i) > tmp(2)) .OR. (Ilist3(i) > tmp(3))) EXIT 1388 IF ((Ilist1(i) == tmp(1)) .AND. (Ilist2(i) == tmp(2)) .AND. (Ilist3(i) == tmp(3))) THEN 1389 found = .TRUE. 1390 EXIT 1391 END IF 1392 END DO 1393 CASE (4) 1394 DO i = istart, iend 1395 IF ((Ilist1(i) > tmp(1)) .OR. (Ilist2(i) > tmp(2)) .OR. (Ilist3(i) > tmp(3)) .OR. (Ilist4(i) > tmp(4))) EXIT 1396 IF ((Ilist1(i) == tmp(1)) .AND. (Ilist2(i) == tmp(2)) & 1397 .AND. (Ilist3(i) == tmp(3)) .AND. (Ilist4(i) == tmp(4))) THEN 1398 found = .TRUE. 1399 EXIT 1400 END IF 1401 END DO 1402 END SELECT 1403 SELECT CASE (do_action) 1404 CASE (do_add) 1405 IF (found) THEN 1406 do_it = -i 1407 ! Nothing to modify. Element already present 1408 ! in this case ABS(do_it) gives the exact location of the element 1409 ! in the list 1410 ELSE 1411 ! Let's add the element in the right place of the list.. so that we can keep the 1412 ! canonical order 1413 ! in this case do_it gives the index of the list with indexes bigger than 1414 ! the one we're searching for 1415 ! At the end do_it gives the exact location of the element in the canonical list 1416 do_it = i 1417 new_size = ndim + 1 1418 CALL reallocate(Ilist1, 1, new_size) 1419 CALL reallocate(Ilist2, 1, new_size) 1420 Ilist1(i + 1:new_size) = Ilist1(i:ndim) 1421 Ilist2(i + 1:new_size) = Ilist2(i:ndim) 1422 Ilist1(i) = tmp(1) 1423 Ilist2(i) = tmp(2) 1424 SELECT CASE (nsize) 1425 CASE (3) 1426 CALL reallocate(Ilist3, 1, new_size) 1427 Ilist3(i + 1:new_size) = Ilist3(i:ndim) 1428 Ilist3(i) = tmp(3) 1429 CASE (4) 1430 CALL reallocate(Ilist3, 1, new_size) 1431 CALL reallocate(Ilist4, 1, new_size) 1432 Ilist3(i + 1:new_size) = Ilist3(i:ndim) 1433 Ilist4(i + 1:new_size) = Ilist4(i:ndim) 1434 Ilist3(i) = tmp(3) 1435 Ilist4(i) = tmp(4) 1436 END SELECT 1437 END IF 1438 CASE (do_remove) 1439 IF (found) THEN 1440 do_it = i 1441 ! Let's delete the element in position do_it 1442 new_size = ndim - 1 1443 Ilist1(i:new_size) = Ilist1(i + 1:ndim) 1444 Ilist2(i:new_size) = Ilist2(i + 1:ndim) 1445 CALL reallocate(Ilist1, 1, new_size) 1446 CALL reallocate(Ilist2, 1, new_size) 1447 SELECT CASE (nsize) 1448 CASE (3) 1449 Ilist3(i:new_size) = Ilist3(i + 1:ndim) 1450 CALL reallocate(Ilist3, 1, new_size) 1451 CASE (4) 1452 Ilist3(i:new_size) = Ilist3(i + 1:ndim) 1453 Ilist4(i:new_size) = Ilist4(i + 1:ndim) 1454 CALL reallocate(Ilist3, 1, new_size) 1455 CALL reallocate(Ilist4, 1, new_size) 1456 END SELECT 1457 ELSE 1458 do_it = -i 1459 ! Nothing to modify. Element not present in the list 1460 ! in this case ABS(do_it) gives the exact location of the element 1461 ! in the list 1462 END IF 1463 END SELECT 1464 END SUBROUTINE check_element_list 1465 1466! ************************************************************************************************** 1467!> \brief Adds a bond to the generated bond list 1468!> \param conn_info ... 1469!> \param atm1 ... 1470!> \param atm2 ... 1471!> \param n_bonds ... 1472!> \author Teodoro Laino 09.2006 1473! ************************************************************************************************** 1474 SUBROUTINE add_bonds_list(conn_info, atm1, atm2, n_bonds) 1475 TYPE(connectivity_info_type), POINTER :: conn_info 1476 INTEGER, INTENT(IN) :: atm1, atm2, n_bonds 1477 1478 INTEGER :: new_size, old_size 1479 1480 old_size = SIZE(conn_info%bond_a) 1481 IF (n_bonds > old_size) THEN 1482 new_size = INT(5 + 1.2*old_size) 1483 CALL reallocate(conn_info%bond_a, 1, new_size) 1484 CALL reallocate(conn_info%bond_b, 1, new_size) 1485 END IF 1486 conn_info%bond_a(n_bonds) = atm1 1487 conn_info%bond_b(n_bonds) = atm2 1488 END SUBROUTINE add_bonds_list 1489 1490! ************************************************************************************************** 1491!> \brief Using a list of bonds, generate a list of bends 1492!> \param topology ... 1493!> \param subsys_section ... 1494!> \author Teodoro Laino 09.2006 1495! ************************************************************************************************** 1496 SUBROUTINE topology_generate_bend(topology, subsys_section) 1497 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 1498 TYPE(section_vals_type), POINTER :: subsys_section 1499 1500 CHARACTER(len=*), PARAMETER :: routineN = 'topology_generate_bend', & 1501 routineP = moduleN//':'//routineN 1502 1503 INTEGER :: handle, handle2, i, iw, natom, nbond, & 1504 nsize, ntheta, output_unit 1505 TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list 1506 TYPE(connectivity_info_type), POINTER :: conn_info 1507 TYPE(cp_logger_type), POINTER :: logger 1508 TYPE(section_vals_type), POINTER :: bend_section 1509 1510 NULLIFY (logger) 1511 logger => cp_get_default_logger() 1512 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", & 1513 extension=".subsysLog") 1514 CALL timeset(routineN, handle) 1515 output_unit = cp_logger_get_default_io_unit(logger) 1516 conn_info => topology%conn_info 1517 nbond = 0 1518 ntheta = 0 1519 natom = topology%natoms 1520 ! This call is for connectivity off 1521 IF (ASSOCIATED(conn_info%bond_a)) THEN 1522 nbond = SIZE(conn_info%bond_a) 1523 ELSE 1524 CALL reallocate(conn_info%bond_a, 1, nbond) 1525 CALL reallocate(conn_info%bond_b, 1, nbond) 1526 END IF 1527 IF (nbond /= 0) THEN 1528 nsize = INT(5 + 1.2*ntheta) 1529 CALL reallocate(conn_info%theta_a, 1, nsize) 1530 CALL reallocate(conn_info%theta_b, 1, nsize) 1531 CALL reallocate(conn_info%theta_c, 1, nsize) 1532 ! Get list of bonds to pre-process theta 1533 ALLOCATE (bond_list(natom)) 1534 DO I = 1, natom 1535 ALLOCATE (bond_list(I)%array1(0)) 1536 ENDDO 1537 CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond) 1538 ! All the dirty job is handled by this routine.. for bends it_levl is equal 3 1539 CALL timeset(routineN//"_1", handle2) 1540 CALL match_iterative_path(Iarray1=bond_list, & 1541 Iarray2=bond_list, & 1542 max_levl=3, & 1543 nvar=ntheta, & 1544 Oarray1=conn_info%theta_a, & 1545 Oarray2=conn_info%theta_b, & 1546 Oarray3=conn_info%theta_c) 1547 CALL timestop(handle2) 1548 DO I = 1, natom 1549 DEALLOCATE (bond_list(I)%array1) 1550 ENDDO 1551 DEALLOCATE (bond_list) 1552 IF (output_unit > 0) THEN 1553 WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Preliminary Number of Bends generated:", & 1554 ntheta 1555 END IF 1556 ! External defined bends (useful for complex connectivity) 1557 bend_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE%ANGLE") 1558 CALL connectivity_external_control(section=bend_section, & 1559 Iarray1=conn_info%theta_a, & 1560 Iarray2=conn_info%theta_b, & 1561 Iarray3=conn_info%theta_c, & 1562 nvar=ntheta, & 1563 topology=topology, & 1564 output_unit=output_unit) 1565 END IF 1566 ! Resize arrays to their proper size.. 1567 CALL reallocate(conn_info%theta_a, 1, ntheta) 1568 CALL reallocate(conn_info%theta_b, 1, ntheta) 1569 CALL reallocate(conn_info%theta_c, 1, ntheta) 1570 IF (output_unit > 0 .AND. ntheta > 0) THEN 1571 WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of Bends generated:", & 1572 ntheta 1573 END IF 1574 CALL timestop(handle) 1575 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 1576 "PRINT%TOPOLOGY_INFO/GENERATE_INFO") 1577 END SUBROUTINE topology_generate_bend 1578 1579! 1580 1581! ************************************************************************************************** 1582!> \brief Routine matching iteratively along a graph 1583!> \param Iarray1 ... 1584!> \param Iarray2 ... 1585!> \param Iarray3 ... 1586!> \param max_levl ... 1587!> \param Oarray1 ... 1588!> \param Oarray2 ... 1589!> \param Oarray3 ... 1590!> \param Oarray4 ... 1591!> \param Ilist ... 1592!> \param it_levl ... 1593!> \param nvar ... 1594!> \author Teodoro Laino 09.2006 1595! ************************************************************************************************** 1596 RECURSIVE SUBROUTINE match_iterative_path(Iarray1, Iarray2, Iarray3, & 1597 max_levl, Oarray1, Oarray2, Oarray3, Oarray4, Ilist, it_levl, nvar) 1598 TYPE(array1_list_type), DIMENSION(:), POINTER :: Iarray1 1599 TYPE(array1_list_type), DIMENSION(:), OPTIONAL, & 1600 POINTER :: Iarray2, Iarray3 1601 INTEGER, INTENT(IN) :: max_levl 1602 INTEGER, DIMENSION(:), POINTER :: Oarray1, Oarray2 1603 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: Oarray3, Oarray4 1604 INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL :: Ilist 1605 INTEGER, INTENT(IN), OPTIONAL :: it_levl 1606 INTEGER, INTENT(INOUT) :: nvar 1607 1608 CHARACTER(len=*), PARAMETER :: routineN = 'match_iterative_path', & 1609 routineP = moduleN//':'//routineN 1610 1611 INTEGER :: i, ind, j, my_levl, natom 1612 INTEGER, ALLOCATABLE, DIMENSION(:) :: my_list 1613 LOGICAL :: check 1614 TYPE(array1_list_type), DIMENSION(:), POINTER :: wrk 1615 1616 check = max_levl >= 2 .AND. max_levl <= 4 1617 CPASSERT(check) 1618 IF (.NOT. PRESENT(Ilist)) THEN 1619 SELECT CASE (max_levl) 1620 CASE (2) 1621 CPASSERT(.NOT. PRESENT(Iarray2)) 1622 CPASSERT(.NOT. PRESENT(Iarray3)) 1623 CPASSERT(.NOT. PRESENT(Oarray3)) 1624 CPASSERT(.NOT. PRESENT(Oarray4)) 1625 CASE (3) 1626 CPASSERT(PRESENT(Iarray2)) 1627 CPASSERT(.NOT. PRESENT(Iarray3)) 1628 CPASSERT(PRESENT(Oarray3)) 1629 CPASSERT(.NOT. PRESENT(Oarray4)) 1630 CASE (4) 1631 CPASSERT(PRESENT(Iarray2)) 1632 CPASSERT(PRESENT(Iarray3)) 1633 CPASSERT(PRESENT(Oarray3)) 1634 CPASSERT(PRESENT(Oarray4)) 1635 END SELECT 1636 END IF 1637 natom = SIZE(Iarray1) 1638 IF (.NOT. PRESENT(Ilist)) THEN 1639 ! Start a new loop.. Only the first time the routine is called 1640 ALLOCATE (my_list(max_levl)) 1641 DO i = 1, natom 1642 my_levl = 1 1643 my_list = -1 1644 my_list(my_levl) = i 1645 CALL match_iterative_path(Iarray1=Iarray1, & 1646 Iarray2=Iarray2, & 1647 Iarray3=Iarray3, & 1648 it_levl=my_levl + 1, & 1649 max_levl=max_levl, & 1650 Oarray1=Oarray1, & 1651 Oarray2=Oarray2, & 1652 Oarray3=Oarray3, & 1653 Oarray4=Oarray4, & 1654 nvar=nvar, & 1655 Ilist=my_list) 1656 END DO 1657 DEALLOCATE (my_list) 1658 ELSE 1659 SELECT CASE (it_levl) 1660 CASE (2) 1661 wrk => Iarray1 1662 CASE (3) 1663 wrk => Iarray2 1664 CASE (4) 1665 wrk => Iarray3 1666 END SELECT 1667 i = Ilist(it_levl - 1) 1668 DO j = 1, SIZE(Iarray1(i)%array1) 1669 ind = wrk(i)%array1(j) 1670 IF (ANY(Ilist == ind)) CYCLE 1671 IF (it_levl < max_levl) THEN 1672 Ilist(it_levl) = ind 1673 CALL match_iterative_path(Iarray1=Iarray1, & 1674 Iarray2=Iarray2, & 1675 Iarray3=Iarray3, & 1676 it_levl=it_levl + 1, & 1677 max_levl=max_levl, & 1678 Oarray1=Oarray1, & 1679 Oarray2=Oarray2, & 1680 Oarray3=Oarray3, & 1681 Oarray4=Oarray4, & 1682 nvar=nvar, & 1683 Ilist=Ilist) 1684 Ilist(it_levl) = -1 1685 ELSEIF (it_levl == max_levl) THEN 1686 IF (Ilist(1) > ind) CYCLE 1687 Ilist(it_levl) = ind 1688 nvar = nvar + 1 1689 SELECT CASE (it_levl) 1690 CASE (2) 1691 IF (nvar > SIZE(Oarray1)) THEN 1692 CALL reallocate(Oarray1, 1, INT(5 + 1.2*nvar)) 1693 CALL reallocate(Oarray2, 1, INT(5 + 1.2*nvar)) 1694 END IF 1695 Oarray1(nvar) = Ilist(1) 1696 Oarray2(nvar) = Ilist(2) 1697 CASE (3) 1698 IF (nvar > SIZE(Oarray1)) THEN 1699 CALL reallocate(Oarray1, 1, INT(5 + 1.2*nvar)) 1700 CALL reallocate(Oarray2, 1, INT(5 + 1.2*nvar)) 1701 CALL reallocate(Oarray3, 1, INT(5 + 1.2*nvar)) 1702 END IF 1703 Oarray1(nvar) = Ilist(1) 1704 Oarray2(nvar) = Ilist(2) 1705 Oarray3(nvar) = Ilist(3) 1706 CASE (4) 1707 IF (nvar > SIZE(Oarray1)) THEN 1708 CALL reallocate(Oarray1, 1, INT(5 + 1.2*nvar)) 1709 CALL reallocate(Oarray2, 1, INT(5 + 1.2*nvar)) 1710 CALL reallocate(Oarray3, 1, INT(5 + 1.2*nvar)) 1711 CALL reallocate(Oarray4, 1, INT(5 + 1.2*nvar)) 1712 END IF 1713 Oarray1(nvar) = Ilist(1) 1714 Oarray2(nvar) = Ilist(2) 1715 Oarray3(nvar) = Ilist(3) 1716 Oarray4(nvar) = Ilist(4) 1717 CASE DEFAULT 1718 !should never reach this point 1719 CPABORT("") 1720 END SELECT 1721 Ilist(it_levl) = -1 1722 ELSE 1723 !should never reach this point 1724 CPABORT("") 1725 END IF 1726 END DO 1727 END IF 1728 END SUBROUTINE match_iterative_path 1729 1730! 1731 1732! ************************************************************************************************** 1733!> \brief The list of Urey-Bradley is equal to the list of bends 1734!> \param topology ... 1735!> \param subsys_section ... 1736! ************************************************************************************************** 1737 SUBROUTINE topology_generate_ub(topology, subsys_section) 1738 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 1739 TYPE(section_vals_type), POINTER :: subsys_section 1740 1741 CHARACTER(len=*), PARAMETER :: routineN = 'topology_generate_ub', & 1742 routineP = moduleN//':'//routineN 1743 1744 INTEGER :: handle, itheta, iw, ntheta, output_unit 1745 TYPE(connectivity_info_type), POINTER :: conn_info 1746 TYPE(cp_logger_type), POINTER :: logger 1747 1748 NULLIFY (logger) 1749 logger => cp_get_default_logger() 1750 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", & 1751 extension=".subsysLog") 1752 output_unit = cp_logger_get_default_io_unit(logger) 1753 CALL timeset(routineN, handle) 1754 conn_info => topology%conn_info 1755 ntheta = SIZE(conn_info%theta_a) 1756 CALL reallocate(conn_info%ub_a, 1, ntheta) 1757 CALL reallocate(conn_info%ub_b, 1, ntheta) 1758 CALL reallocate(conn_info%ub_c, 1, ntheta) 1759 1760 DO itheta = 1, ntheta 1761 conn_info%ub_a(itheta) = conn_info%theta_a(itheta) 1762 conn_info%ub_b(itheta) = conn_info%theta_b(itheta) 1763 conn_info%ub_c(itheta) = conn_info%theta_c(itheta) 1764 END DO 1765 IF (output_unit > 0 .AND. ntheta > 0) THEN 1766 WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of UB generated:", & 1767 ntheta 1768 END IF 1769 CALL timestop(handle) 1770 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 1771 "PRINT%TOPOLOGY_INFO/GENERATE_INFO") 1772 1773 END SUBROUTINE topology_generate_ub 1774 1775! ************************************************************************************************** 1776!> \brief Generate a list of torsions from bonds 1777!> \param topology ... 1778!> \param subsys_section ... 1779!> \author Teodoro Laino 09.2006 1780! ************************************************************************************************** 1781 SUBROUTINE topology_generate_dihe(topology, subsys_section) 1782 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 1783 TYPE(section_vals_type), POINTER :: subsys_section 1784 1785 CHARACTER(len=*), PARAMETER :: routineN = 'topology_generate_dihe', & 1786 routineP = moduleN//':'//routineN 1787 1788 INTEGER :: handle, i, iw, natom, nbond, nphi, & 1789 nsize, output_unit 1790 TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list 1791 TYPE(connectivity_info_type), POINTER :: conn_info 1792 TYPE(cp_logger_type), POINTER :: logger 1793 TYPE(section_vals_type), POINTER :: torsion_section 1794 1795 NULLIFY (logger) 1796 logger => cp_get_default_logger() 1797 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", & 1798 extension=".subsysLog") 1799 output_unit = cp_logger_get_default_io_unit(logger) 1800 CALL timeset(routineN, handle) 1801 conn_info => topology%conn_info 1802 nphi = 0 1803 nbond = SIZE(conn_info%bond_a) 1804 IF (nbond /= 0) THEN 1805 nsize = INT(5 + 1.2*nphi) 1806 CALL reallocate(conn_info%phi_a, 1, nsize) 1807 CALL reallocate(conn_info%phi_b, 1, nsize) 1808 CALL reallocate(conn_info%phi_c, 1, nsize) 1809 CALL reallocate(conn_info%phi_d, 1, nsize) 1810 ! Get list of bonds to pre-process phi 1811 natom = topology%natoms 1812 ALLOCATE (bond_list(natom)) 1813 DO I = 1, natom 1814 ALLOCATE (bond_list(I)%array1(0)) 1815 ENDDO 1816 CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond) 1817 ! All the dirty job is handled by this routine.. for torsions it_levl is equal 4 1818 CALL match_iterative_path(Iarray1=bond_list, & 1819 Iarray2=bond_list, & 1820 Iarray3=bond_list, & 1821 max_levl=4, & 1822 nvar=nphi, & 1823 Oarray1=conn_info%phi_a, & 1824 Oarray2=conn_info%phi_b, & 1825 Oarray3=conn_info%phi_c, & 1826 Oarray4=conn_info%phi_d) 1827 DO I = 1, natom 1828 DEALLOCATE (bond_list(I)%array1) 1829 ENDDO 1830 DEALLOCATE (bond_list) 1831 IF (output_unit > 0) THEN 1832 WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Preliminary Number of Torsions generated:", & 1833 nphi 1834 END IF 1835 ! External defined torsions (useful for complex connectivity) 1836 torsion_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE%TORSION") 1837 CALL connectivity_external_control(section=torsion_section, & 1838 Iarray1=conn_info%phi_a, & 1839 Iarray2=conn_info%phi_b, & 1840 Iarray3=conn_info%phi_c, & 1841 Iarray4=conn_info%phi_d, & 1842 nvar=nphi, & 1843 topology=topology, & 1844 output_unit=output_unit) 1845 END IF 1846 ! Resize arrays to their proper size.. 1847 CALL reallocate(conn_info%phi_a, 1, nphi) 1848 CALL reallocate(conn_info%phi_b, 1, nphi) 1849 CALL reallocate(conn_info%phi_c, 1, nphi) 1850 CALL reallocate(conn_info%phi_d, 1, nphi) 1851 IF (output_unit > 0 .AND. nphi > 0) THEN 1852 WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of Torsions generated:", & 1853 nphi 1854 END IF 1855 CALL timestop(handle) 1856 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 1857 "PRINT%TOPOLOGY_INFO/GENERATE_INFO") 1858 1859 END SUBROUTINE topology_generate_dihe 1860 1861! ************************************************************************************************** 1862!> \brief Using a list of bends, generate a list of impr 1863!> \param topology ... 1864!> \param subsys_section ... 1865!> \author Teodoro Laino 09.2006 1866! ************************************************************************************************** 1867 SUBROUTINE topology_generate_impr(topology, subsys_section) 1868 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 1869 TYPE(section_vals_type), POINTER :: subsys_section 1870 1871 CHARACTER(len=*), PARAMETER :: routineN = 'topology_generate_impr', & 1872 routineP = moduleN//':'//routineN 1873 1874 CHARACTER(LEN=2) :: atm_symbol 1875 INTEGER :: handle, i, ind, iw, j, natom, nbond, & 1876 nimpr, nsize, output_unit 1877 LOGICAL :: accept_impr 1878 TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list 1879 TYPE(atom_info_type), POINTER :: atom_info 1880 TYPE(connectivity_info_type), POINTER :: conn_info 1881 TYPE(cp_logger_type), POINTER :: logger 1882 TYPE(section_vals_type), POINTER :: impr_section 1883 1884 NULLIFY (logger) 1885 logger => cp_get_default_logger() 1886 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", & 1887 extension=".subsysLog") 1888 output_unit = cp_logger_get_default_io_unit(logger) 1889 CALL timeset(routineN, handle) 1890 atom_info => topology%atom_info 1891 conn_info => topology%conn_info 1892 natom = topology%natoms 1893 nimpr = 0 1894 nbond = SIZE(conn_info%bond_a) 1895 IF (nbond /= 0) THEN 1896 nsize = INT(5 + 1.2*nimpr) 1897 CALL reallocate(conn_info%impr_a, 1, nsize) 1898 CALL reallocate(conn_info%impr_b, 1, nsize) 1899 CALL reallocate(conn_info%impr_c, 1, nsize) 1900 CALL reallocate(conn_info%impr_d, 1, nsize) 1901 ! Get list of bonds to pre-process phi 1902 ALLOCATE (bond_list(natom)) 1903 DO I = 1, natom 1904 ALLOCATE (bond_list(I)%array1(0)) 1905 ENDDO 1906 CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond) 1907 DO I = 1, natom 1908 ! Count all atoms with three bonds 1909 IF (SIZE(bond_list(I)%array1) == 3) THEN 1910 ! Problematic cases:: 1911 ! Nitrogen 1912 accept_impr = .TRUE. 1913 atm_symbol = id2str(atom_info%id_element(i)) 1914 CALL uppercase(atm_symbol) 1915 IF (atm_symbol == "N ") THEN 1916 accept_impr = .FALSE. 1917 ! Impropers on Nitrogen only when there is another atom close to it 1918 ! with other 3 bonds 1919 DO j = 1, 3 1920 ind = bond_list(I)%array1(j) 1921 IF (SIZE(bond_list(ind)%array1) == 3) accept_impr = .TRUE. 1922 END DO 1923 END IF 1924 IF (.NOT. accept_impr) CYCLE 1925 nimpr = nimpr + 1 1926 IF (nimpr > SIZE(conn_info%impr_a)) THEN 1927 nsize = INT(5 + 1.2*nimpr) 1928 CALL reallocate(conn_info%impr_a, 1, nsize) 1929 CALL reallocate(conn_info%impr_b, 1, nsize) 1930 CALL reallocate(conn_info%impr_c, 1, nsize) 1931 CALL reallocate(conn_info%impr_d, 1, nsize) 1932 END IF 1933 conn_info%impr_a(nimpr) = i 1934 conn_info%impr_b(nimpr) = bond_list(I)%array1(1) 1935 conn_info%impr_c(nimpr) = bond_list(I)%array1(2) 1936 conn_info%impr_d(nimpr) = bond_list(I)%array1(3) 1937 END IF 1938 END DO 1939 DO I = 1, natom 1940 DEALLOCATE (bond_list(I)%array1) 1941 ENDDO 1942 DEALLOCATE (bond_list) 1943 ! External defined impropers (useful for complex connectivity) 1944 impr_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE%IMPROPER") 1945 CALL connectivity_external_control(section=impr_section, & 1946 Iarray1=conn_info%impr_a, & 1947 Iarray2=conn_info%impr_b, & 1948 Iarray3=conn_info%impr_c, & 1949 Iarray4=conn_info%impr_d, & 1950 nvar=nimpr, & 1951 topology=topology, & 1952 output_unit=output_unit, & 1953 is_impr=.TRUE.) 1954 END IF 1955 ! Resize arrays to their proper size.. 1956 CALL reallocate(conn_info%impr_a, 1, nimpr) 1957 CALL reallocate(conn_info%impr_b, 1, nimpr) 1958 CALL reallocate(conn_info%impr_c, 1, nimpr) 1959 CALL reallocate(conn_info%impr_d, 1, nimpr) 1960 IF (output_unit > 0 .AND. nimpr > 0) THEN 1961 WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of Impropers generated:", & 1962 nimpr 1963 END IF 1964 CALL timestop(handle) 1965 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 1966 "PRINT%TOPOLOGY_INFO/GENERATE_INFO") 1967 1968 END SUBROUTINE topology_generate_impr 1969 1970! ************************************************************************************************** 1971!> \brief Using a list of torsion, generate a list of onfo 1972!> \param topology ... 1973!> \param subsys_section ... 1974! ************************************************************************************************** 1975 SUBROUTINE topology_generate_onfo(topology, subsys_section) 1976 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 1977 TYPE(section_vals_type), POINTER :: subsys_section 1978 1979 CHARACTER(len=*), PARAMETER :: routineN = 'topology_generate_onfo', & 1980 routineP = moduleN//':'//routineN 1981 1982 INTEGER :: atom_a, atom_b, handle, i, ionfo, iw, & 1983 natom, nbond, nphi, ntheta, output_unit 1984 TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list, phi_list, theta_list 1985 TYPE(connectivity_info_type), POINTER :: conn_info 1986 TYPE(cp_logger_type), POINTER :: logger 1987 1988 NULLIFY (logger) 1989 logger => cp_get_default_logger() 1990 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GENERATE_INFO", & 1991 extension=".subsysLog") 1992 output_unit = cp_logger_get_default_io_unit(logger) 1993 CALL timeset(routineN, handle) 1994 1995 conn_info => topology%conn_info 1996 natom = topology%natoms 1997 1998 ! Get list of bonds (sic). Get a list of bonded neighbors for every atom. 1999 ALLOCATE (bond_list(natom)) 2000 DO i = 1, natom 2001 ALLOCATE (bond_list(i)%array1(0)) 2002 ENDDO 2003 nbond = SIZE(conn_info%bond_a) 2004 CALL reorder_structure(bond_list, conn_info%bond_a, conn_info%bond_b, nbond) 2005 2006 ! Get a list of next nearest neighbors for every atom. 2007 ALLOCATE (theta_list(natom)) 2008 DO i = 1, natom 2009 ALLOCATE (theta_list(i)%array1(0)) 2010 ENDDO 2011 ntheta = SIZE(conn_info%theta_a) 2012 CALL reorder_structure(theta_list, conn_info%theta_a, conn_info%theta_c, ntheta) 2013 2014 ! Get a list of next next nearest neighbors for every atom. 2015 ALLOCATE (phi_list(natom)) 2016 DO i = 1, natom 2017 ALLOCATE (phi_list(i)%array1(0)) 2018 ENDDO 2019 nphi = SIZE(conn_info%phi_a) 2020 CALL reorder_structure(phi_list, conn_info%phi_a, conn_info%phi_d, nphi) 2021 2022 ! Allocate enough (possible too much) 2023 CALL reallocate(conn_info%onfo_a, 1, nphi) 2024 CALL reallocate(conn_info%onfo_b, 1, nphi) 2025 2026 ionfo = 0 2027 DO atom_a = 1, natom 2028 DO i = 1, SIZE(phi_list(atom_a)%array1) 2029 atom_b = phi_list(atom_a)%array1(i) 2030 ! Avoid trivial duplicates. 2031 IF (atom_a > atom_b) CYCLE 2032 ! Avoid onfo's in 4-rings. 2033 IF (ANY(atom_b == bond_list(atom_a)%array1)) CYCLE 2034 ! Avoid onfo's in 5-rings. 2035 IF (ANY(atom_b == theta_list(atom_a)%array1)) CYCLE 2036 ! Avoid onfo's in 6-rings. 2037 IF (ANY(atom_b == phi_list(atom_a)%array1(:i - 1))) CYCLE 2038 ionfo = ionfo + 1 2039 conn_info%onfo_a(ionfo) = atom_a 2040 conn_info%onfo_b(ionfo) = atom_b 2041 END DO 2042 END DO 2043 2044 ! Reallocate such that just enough memory is used. 2045 CALL reallocate(conn_info%onfo_a, 1, ionfo) 2046 CALL reallocate(conn_info%onfo_b, 1, ionfo) 2047 2048 ! Deallocate bond_list 2049 DO i = 1, natom 2050 DEALLOCATE (bond_list(i)%array1) 2051 ENDDO 2052 DEALLOCATE (bond_list) 2053 ! Deallocate theta_list 2054 DO i = 1, natom 2055 DEALLOCATE (theta_list(i)%array1) 2056 ENDDO 2057 DEALLOCATE (theta_list) 2058 ! Deallocate phi_list 2059 DO i = 1, natom 2060 DEALLOCATE (phi_list(i)%array1) 2061 ENDDO 2062 DEALLOCATE (phi_list) 2063 2064 ! Final output 2065 IF (output_unit > 0 .AND. ionfo > 0) THEN 2066 WRITE (output_unit, '(T2,"GENERATE|",1X,A,T71,I10)') " Number of 1-4 interactions generated:", & 2067 ionfo 2068 END IF 2069 CALL timestop(handle) 2070 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 2071 "PRINT%TOPOLOGY_INFO/GENERATE_INFO") 2072 2073 END SUBROUTINE topology_generate_onfo 2074 2075END MODULE topology_generate_util 2076