LCOV - code coverage report
Current view: top level - src - topology_generate_util.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 1002 1066 94.0 %
Date: 2024-03-28 07:31:50 Functions: 15 15 100.0 %

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

Generated by: LCOV version 1.15