LCOV - code coverage report
Current view: top level - src - topology_util.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 574 592 97.0 %
Date: 2024-03-28 07:31:50 Functions: 15 17 88.2 %

          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             : !>     jgh (23-05-2004) Last atom of molecule information added
      12             : ! **************************************************************************************************
      13             : MODULE topology_util
      14             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      15             :                                               cp_logger_get_default_io_unit,&
      16             :                                               cp_logger_type,&
      17             :                                               cp_to_string
      18             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      19             :                                               cp_print_key_unit_nr
      20             :    USE graphcon,                        ONLY: graph_type,&
      21             :                                               hash_molecule,&
      22             :                                               reorder_graph,&
      23             :                                               vertex
      24             :    USE input_section_types,             ONLY: section_vals_get,&
      25             :                                               section_vals_get_subs_vals,&
      26             :                                               section_vals_type,&
      27             :                                               section_vals_val_get,&
      28             :                                               section_vals_val_set
      29             :    USE kinds,                           ONLY: default_string_length,&
      30             :                                               dp
      31             :    USE mm_mapping_library,              ONLY: amber_map,&
      32             :                                               charmm_map,&
      33             :                                               gromos_map
      34             :    USE periodic_table,                  ONLY: get_ptable_info
      35             :    USE qmmm_types_low,                  ONLY: qmmm_env_mm_type
      36             :    USE string_table,                    ONLY: id2str,&
      37             :                                               str2id
      38             :    USE string_utilities,                ONLY: uppercase
      39             :    USE topology_types,                  ONLY: atom_info_type,&
      40             :                                               connectivity_info_type,&
      41             :                                               topology_parameters_type
      42             :    USE util,                            ONLY: sort
      43             : #include "./base/base_uses.f90"
      44             : 
      45             :    IMPLICIT NONE
      46             : 
      47             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_util'
      48             : 
      49             : ! **************************************************************************************************
      50             :    TYPE array1_list_type
      51             :       INTEGER, POINTER, DIMENSION(:) :: array1 => NULL()
      52             :    END TYPE array1_list_type
      53             : 
      54             : ! **************************************************************************************************
      55             :    TYPE array2_list_type
      56             :       INTEGER, POINTER, DIMENSION(:) :: array1 => NULL()
      57             :       INTEGER, POINTER, DIMENSION(:) :: array2 => NULL()
      58             :    END TYPE array2_list_type
      59             : 
      60             :    PRIVATE
      61             :    PUBLIC :: topology_set_atm_mass, &
      62             :              topology_reorder_atoms, &
      63             :              topology_molecules_check, &
      64             :              check_subsys_element, &
      65             :              reorder_structure, &
      66             :              find_molecule, &
      67             :              array1_list_type, &
      68             :              array2_list_type, &
      69             :              give_back_molecule, &
      70             :              reorder_list_array, &
      71             :              tag_molecule
      72             : 
      73             :    INTERFACE reorder_structure
      74             :       MODULE PROCEDURE reorder_structure1d, reorder_structure2d
      75             :    END INTERFACE
      76             : 
      77             : CONTAINS
      78             : 
      79             : ! **************************************************************************************************
      80             : !> \brief ...
      81             : !> \param topology ...
      82             : !> \param qmmm ...
      83             : !> \param qmmm_env_mm ...
      84             : !> \param subsys_section ...
      85             : !> \param force_env_section ...
      86             : !> \par History
      87             : !>      Teodoro Laino 09.2006 - Rewritten with a graph matching algorithm
      88             : ! **************************************************************************************************
      89         314 :    SUBROUTINE topology_reorder_atoms(topology, qmmm, qmmm_env_mm, subsys_section, force_env_section)
      90             :       TYPE(topology_parameters_type), INTENT(INOUT)      :: topology
      91             :       LOGICAL, INTENT(in), OPTIONAL                      :: qmmm
      92             :       TYPE(qmmm_env_mm_type), OPTIONAL, POINTER          :: qmmm_env_mm
      93             :       TYPE(section_vals_type), POINTER                   :: subsys_section, force_env_section
      94             : 
      95             :       CHARACTER(len=*), PARAMETER :: routineN = 'topology_reorder_atoms'
      96             : 
      97             :       CHARACTER(LEN=default_string_length)               :: mol_id
      98         314 :       CHARACTER(LEN=default_string_length), POINTER      :: molname(:), telement(:), &
      99         314 :                                                             tlabel_atmname(:), tlabel_molname(:), &
     100         314 :                                                             tlabel_resname(:)
     101             :       INTEGER :: handle, i, iatm, iindex, ikind, imol, imol_ref, iref, iund, iw, j, k, location, &
     102             :          max_mol_num, mm_index, n, n_rep, n_var, natom, natom_loc, nkind, nlinks, old_hash, &
     103             :          old_mol, output_unit, qm_index, unique_mol
     104         314 :       INTEGER, DIMENSION(:), POINTER                     :: mm_link_atoms, qm_atom_index
     105         314 :       INTEGER, POINTER :: atm_map1(:), atm_map2(:), atm_map3(:), map_atom_type(:), &
     106         314 :          map_mol_hash(:), mm_indexes_n(:), mm_indexes_v(:), mol_bnd(:, :), mol_hash(:), &
     107         314 :          mol_num(:), new_position(:), order(:), tmp_n(:), tmp_v(:), wrk(:)
     108             :       LOGICAL                                            :: explicit, matches, my_qmmm
     109         314 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: tr
     110         314 :       REAL(KIND=dp), POINTER                             :: tatm_charge(:), tatm_mass(:)
     111         314 :       TYPE(array1_list_type), ALLOCATABLE, DIMENSION(:)  :: atom_bond_list
     112             :       TYPE(atom_info_type), POINTER                      :: atom_info
     113             :       TYPE(connectivity_info_type), POINTER              :: conn_info
     114             :       TYPE(cp_logger_type), POINTER                      :: logger
     115         314 :       TYPE(graph_type), DIMENSION(:), POINTER            :: reference_set
     116             :       TYPE(section_vals_type), POINTER                   :: generate_section, isolated_section, &
     117             :                                                             qm_kinds, qmmm_link_section, &
     118             :                                                             qmmm_section
     119         314 :       TYPE(vertex), DIMENSION(:), POINTER                :: reference, unordered
     120             : 
     121         314 :       NULLIFY (logger, generate_section, isolated_section, tmp_v, tmp_n)
     122         628 :       logger => cp_get_default_logger()
     123             :       iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
     124         314 :                                 extension=".subsysLog")
     125         314 :       CALL timeset(routineN, handle)
     126         314 :       output_unit = cp_logger_get_default_io_unit(logger)
     127         314 :       IF (output_unit > 0) WRITE (output_unit, '(T2,"REORDER |  ")')
     128             : 
     129         314 :       my_qmmm = .FALSE.
     130         314 :       IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env_mm)) my_qmmm = qmmm
     131             : 
     132         314 :       atom_info => topology%atom_info
     133         314 :       conn_info => topology%conn_info
     134         314 :       natom = topology%natoms
     135             : 
     136         314 :       NULLIFY (new_position, reference_set)
     137         314 :       NULLIFY (tlabel_atmname, telement, mol_num, tlabel_molname, tlabel_resname)
     138         314 :       NULLIFY (tr, tatm_charge, tatm_mass, atm_map1, atm_map2, atm_map3)
     139             :       ! This routine can be called only at a very high level where these structures are still
     140             :       ! not even taken into account...
     141         314 :       CPASSERT(.NOT. ASSOCIATED(atom_info%map_mol_num))
     142         314 :       CPASSERT(.NOT. ASSOCIATED(atom_info%map_mol_typ))
     143         314 :       CPASSERT(.NOT. ASSOCIATED(atom_info%map_mol_res))
     144             :       !ALLOCATE all the temporary arrays needed for this routine
     145         942 :       ALLOCATE (new_position(natom))
     146         942 :       ALLOCATE (mol_num(natom))
     147         942 :       ALLOCATE (molname(natom))
     148         942 :       ALLOCATE (tlabel_atmname(natom))
     149         942 :       ALLOCATE (tlabel_molname(natom))
     150         942 :       ALLOCATE (tlabel_resname(natom))
     151         942 :       ALLOCATE (tr(3, natom))
     152         942 :       ALLOCATE (tatm_charge(natom))
     153         942 :       ALLOCATE (tatm_mass(natom))
     154         942 :       ALLOCATE (telement(natom))
     155         942 :       ALLOCATE (atm_map1(natom))
     156         942 :       ALLOCATE (atm_map2(natom))
     157             : 
     158             :       ! The only information we have at this level is the connectivity of the system.
     159             :       ! 0. Build a list of mapping atom types
     160         942 :       ALLOCATE (map_atom_type(natom))
     161             :       ! 1. Build a list of bonds
     162        9500 :       ALLOCATE (atom_bond_list(natom))
     163        8872 :       DO I = 1, natom
     164        8558 :          map_atom_type(I) = atom_info%id_atmname(i)
     165        8872 :          ALLOCATE (atom_bond_list(I)%array1(0))
     166             :       END DO
     167         314 :       N = 0
     168         314 :       IF (ASSOCIATED(conn_info%bond_a)) N = SIZE(conn_info%bond_a)
     169         314 :       CALL reorder_structure(atom_bond_list, conn_info%bond_a, conn_info%bond_b, N)
     170         942 :       ALLOCATE (atom_info%map_mol_num(natom))
     171        8872 :       atom_info%map_mol_num = -1
     172         314 :       CALL find_molecule(atom_bond_list, atom_info%map_mol_num, atom_info%id_molname)
     173        8872 :       max_mol_num = MAXVAL(atom_info%map_mol_num)
     174             :       ! In atom_info%map_mol_num have already been mapped molecules
     175         942 :       ALLOCATE (mol_bnd(2, max_mol_num))
     176         942 :       ALLOCATE (mol_hash(max_mol_num))
     177         942 :       ALLOCATE (map_mol_hash(max_mol_num))
     178             :       ! 2. Sort the map_mol_num array.. atm_map1 contains now the mapped index
     179             :       !    of the reordered array
     180         314 :       CALL sort(atom_info%map_mol_num, natom, atm_map1)
     181         314 :       old_mol = 0
     182         314 :       iindex = 0
     183         314 :       imol = 0
     184        8872 :       DO i = 1, natom
     185        8558 :          IF (old_mol .NE. atom_info%map_mol_num(I)) THEN
     186        1420 :             old_mol = atom_info%map_mol_num(I)
     187        1420 :             iindex = 0
     188        1420 :             IF (imol > 0) THEN
     189        1106 :                mol_bnd(2, imol) = i - 1
     190             :             END IF
     191        1420 :             imol = imol + 1
     192        1420 :             mol_bnd(1, imol) = i
     193             :          END IF
     194        8558 :          iindex = iindex + 1
     195        8872 :          atm_map2(atm_map1(i)) = iindex
     196             :       END DO
     197         314 :       mol_bnd(2, imol) = natom
     198             :       ! Indexes of the two molecules to check
     199         314 :       iref = 1
     200         314 :       iund = max_mol_num/2 + 1
     201             :       ! Allocate reference and unordered
     202         314 :       NULLIFY (reference, unordered)
     203             :       ! This is the real matching of graphs
     204        1734 :       DO j = 1, max_mol_num
     205             :          CALL setup_graph(j, reference, map_atom_type, &
     206        1420 :                           atom_bond_list, mol_bnd, atm_map1, atm_map2)
     207             : 
     208        4260 :          ALLOCATE (order(SIZE(reference)))
     209        1420 :          CALL hash_molecule(reference, order, mol_hash(j))
     210             : 
     211        1420 :          DEALLOCATE (order)
     212        9978 :          DO I = 1, SIZE(reference)
     213        9978 :             DEALLOCATE (reference(I)%bonds)
     214             :          END DO
     215        1734 :          DEALLOCATE (reference)
     216             :       END DO
     217             :       ! Reorder molecules hashes
     218         314 :       CALL sort(mol_hash, max_mol_num, map_mol_hash)
     219             :       ! Now find unique molecules and reorder atoms too (if molecules match)
     220         314 :       old_hash = -1
     221         314 :       unique_mol = 0
     222         314 :       natom_loc = 0
     223         314 :       IF (output_unit > 0) THEN
     224             :          WRITE (output_unit, '(T2,"REORDER |  ",A)') &
     225         157 :             "Reordering Molecules. The Reordering of molecules will override all", &
     226         157 :             "information regarding molecule names and residue names.", &
     227         314 :             "New ones will be provided on the basis of the connectivity!"
     228             :       END IF
     229        1734 :       DO j = 1, max_mol_num
     230        1734 :          IF (mol_hash(j) .NE. old_hash) THEN
     231         326 :             unique_mol = unique_mol + 1
     232         326 :             old_hash = mol_hash(j)
     233             :             CALL setup_graph_set(reference_set, unique_mol, map_mol_hash(j), &
     234             :                                  map_atom_type, atom_bond_list, mol_bnd, atm_map1, atm_map2, &
     235         326 :                                  atm_map3)
     236             :             ! Reorder Last added reference
     237         326 :             mol_id = TRIM(ADJUSTL(cp_to_string(unique_mol)))
     238        5636 :             DO i = 1, SIZE(atm_map3)
     239        5310 :                natom_loc = natom_loc + 1
     240        5310 :                new_position(natom_loc) = atm_map3(i)
     241        5310 :                molname(natom_loc) = mol_id
     242        5636 :                mol_num(natom_loc) = unique_mol
     243             :             END DO
     244         326 :             DEALLOCATE (atm_map3)
     245             :          ELSE
     246             :             CALL setup_graph(map_mol_hash(j), unordered, map_atom_type, &
     247        1094 :                              atom_bond_list, mol_bnd, atm_map1, atm_map2, atm_map3)
     248        1150 :             DO imol_ref = 1, unique_mol
     249             :                !
     250        1150 :                reference => reference_set(imol_ref)%graph
     251        3450 :                ALLOCATE (order(SIZE(reference)))
     252        1150 :                CALL reorder_graph(reference, unordered, order, matches)
     253        1150 :                IF (matches) EXIT
     254        2300 :                DEALLOCATE (order)
     255             :             END DO
     256        1094 :             IF (matches) THEN
     257             :                ! Reorder according the correct index
     258        3282 :                ALLOCATE (wrk(SIZE(order)))
     259        1094 :                CALL sort(order, SIZE(order), wrk)
     260        4342 :                DO i = 1, SIZE(order)
     261        3248 :                   natom_loc = natom_loc + 1
     262        3248 :                   new_position(natom_loc) = atm_map3(wrk(i))
     263        3248 :                   molname(natom_loc) = mol_id
     264        4342 :                   mol_num(natom_loc) = unique_mol
     265             :                END DO
     266             :                !
     267        1094 :                DEALLOCATE (order)
     268        1094 :                DEALLOCATE (wrk)
     269             :             ELSE
     270           0 :                unique_mol = unique_mol + 1
     271             :                CALL setup_graph_set(reference_set, unique_mol, map_mol_hash(j), &
     272             :                                     map_atom_type, atom_bond_list, mol_bnd, atm_map1, atm_map2, &
     273           0 :                                     atm_map3)
     274             :                ! Reorder Last added reference
     275           0 :                mol_id = TRIM(ADJUSTL(cp_to_string(unique_mol)))
     276           0 :                DO i = 1, SIZE(atm_map3)
     277           0 :                   natom_loc = natom_loc + 1
     278           0 :                   new_position(natom_loc) = atm_map3(i)
     279           0 :                   molname(natom_loc) = mol_id
     280           0 :                   mol_num(natom_loc) = unique_mol
     281             :                END DO
     282           0 :                DEALLOCATE (atm_map3)
     283             :             END IF
     284        4342 :             DO I = 1, SIZE(unordered)
     285        4342 :                DEALLOCATE (unordered(I)%bonds)
     286             :             END DO
     287        1094 :             DEALLOCATE (unordered)
     288        1094 :             DEALLOCATE (atm_map3)
     289             :          END IF
     290             :       END DO
     291         314 :       IF (output_unit > 0) THEN
     292         157 :          WRITE (output_unit, '(T2,"REORDER |  ",A,I7,A)') "Number of unique molecules found:", unique_mol, "."
     293             :       END IF
     294         314 :       CPASSERT(natom_loc == natom)
     295         314 :       DEALLOCATE (map_atom_type)
     296         314 :       DEALLOCATE (atm_map1)
     297         314 :       DEALLOCATE (atm_map2)
     298         314 :       DEALLOCATE (mol_bnd)
     299         314 :       DEALLOCATE (mol_hash)
     300         314 :       DEALLOCATE (map_mol_hash)
     301             :       ! Deallocate working arrays
     302        8872 :       DO I = 1, natom
     303        8872 :          DEALLOCATE (atom_bond_list(I)%array1)
     304             :       END DO
     305         314 :       DEALLOCATE (atom_bond_list)
     306         314 :       DEALLOCATE (atom_info%map_mol_num)
     307             :       ! Deallocate reference_set
     308         640 :       DO i = 1, SIZE(reference_set)
     309        5636 :          DO j = 1, SIZE(reference_set(i)%graph)
     310        5636 :             DEALLOCATE (reference_set(i)%graph(j)%bonds)
     311             :          END DO
     312         640 :          DEALLOCATE (reference_set(i)%graph)
     313             :       END DO
     314         314 :       DEALLOCATE (reference_set)
     315             :       !Lets swap the atoms now
     316        8872 :       DO iatm = 1, natom
     317        8558 :          location = new_position(iatm)
     318        8558 :          tlabel_molname(iatm) = id2str(atom_info%id_molname(location))
     319        8558 :          tlabel_resname(iatm) = id2str(atom_info%id_resname(location))
     320        8558 :          tlabel_atmname(iatm) = id2str(atom_info%id_atmname(location))
     321        8558 :          telement(iatm) = id2str(atom_info%id_element(location))
     322        8558 :          tr(1, iatm) = atom_info%r(1, location)
     323        8558 :          tr(2, iatm) = atom_info%r(2, location)
     324        8558 :          tr(3, iatm) = atom_info%r(3, location)
     325        8558 :          tatm_charge(iatm) = atom_info%atm_charge(location)
     326        8872 :          tatm_mass(iatm) = atom_info%atm_mass(location)
     327             :       END DO
     328         314 :       IF (topology%create_molecules) THEN
     329        8858 :          DO iatm = 1, natom
     330        8546 :             tlabel_molname(iatm) = "MOL"//TRIM(molname(iatm))
     331        8858 :             tlabel_resname(iatm) = "R"//TRIM(molname(iatm))
     332             :          END DO
     333         312 :          topology%create_molecules = .FALSE.
     334             :       END IF
     335        8872 :       DO iatm = 1, natom
     336        8558 :          atom_info%id_molname(iatm) = str2id(tlabel_molname(iatm))
     337        8558 :          atom_info%id_resname(iatm) = str2id(tlabel_resname(iatm))
     338        8558 :          atom_info%id_atmname(iatm) = str2id(tlabel_atmname(iatm))
     339        8558 :          atom_info%id_element(iatm) = str2id(telement(iatm))
     340        8558 :          atom_info%resid(iatm) = mol_num(iatm)
     341        8558 :          atom_info%r(1, iatm) = tr(1, iatm)
     342        8558 :          atom_info%r(2, iatm) = tr(2, iatm)
     343        8558 :          atom_info%r(3, iatm) = tr(3, iatm)
     344        8558 :          atom_info%atm_charge(iatm) = tatm_charge(iatm)
     345        8872 :          atom_info%atm_mass(iatm) = tatm_mass(iatm)
     346             :       END DO
     347             : 
     348             :       ! Let's reorder all the list provided in the input..
     349         942 :       ALLOCATE (wrk(SIZE(new_position)))
     350         314 :       CALL sort(new_position, SIZE(new_position), wrk)
     351             : 
     352             :       ! NOTE: In the future the whole list of possible integers should be updated at this level..
     353             :       ! Let's update all the integer lists in the qmmm_env_mm section and in the input sections
     354             :       !       from where qmmm_env_qm will be read.
     355         314 :       IF (my_qmmm) THEN
     356             :          ! Update the qmmm_env_mm
     357           2 :          CPASSERT(SIZE(qmmm_env_mm%qm_atom_index) /= 0)
     358           8 :          CPASSERT(ALL(qmmm_env_mm%qm_atom_index /= 0))
     359           6 :          ALLOCATE (qm_atom_index(SIZE(qmmm_env_mm%qm_atom_index)))
     360           8 :          DO iatm = 1, SIZE(qmmm_env_mm%qm_atom_index)
     361           8 :             qm_atom_index(iatm) = wrk(qmmm_env_mm%qm_atom_index(iatm))
     362             :          END DO
     363          16 :          qmmm_env_mm%qm_atom_index = qm_atom_index
     364           8 :          CPASSERT(ALL(qmmm_env_mm%qm_atom_index /= 0))
     365           2 :          DEALLOCATE (qm_atom_index)
     366             :          ! Update the qmmm_section: MM_INDEX of the QM_KIND
     367           2 :          qmmm_section => section_vals_get_subs_vals(force_env_section, "QMMM")
     368           2 :          qm_kinds => section_vals_get_subs_vals(qmmm_section, "QM_KIND")
     369           2 :          CALL section_vals_get(qm_kinds, n_repetition=nkind)
     370           6 :          DO ikind = 1, nkind
     371           4 :             CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
     372          10 :             DO k = 1, n_var
     373             :                CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
     374           4 :                                          i_vals=mm_indexes_v)
     375          12 :                ALLOCATE (mm_indexes_n(SIZE(mm_indexes_v)))
     376           8 :                DO j = 1, SIZE(mm_indexes_v)
     377           8 :                   mm_indexes_n(j) = wrk(mm_indexes_v(j))
     378             :                END DO
     379             :                CALL section_vals_val_set(qm_kinds, "MM_INDEX", i_rep_section=ikind, &
     380           8 :                                          i_vals_ptr=mm_indexes_n, i_rep_val=k)
     381             :             END DO
     382             :          END DO
     383             :          ! Handle the link atoms
     384           4 :          IF (qmmm_env_mm%qmmm_link) THEN
     385             :             ! Update the qmmm_env_mm
     386           2 :             CPASSERT(SIZE(qmmm_env_mm%mm_link_atoms) > 0)
     387           6 :             ALLOCATE (mm_link_atoms(SIZE(qmmm_env_mm%mm_link_atoms)))
     388           4 :             DO iatm = 1, SIZE(qmmm_env_mm%mm_link_atoms)
     389           4 :                mm_link_atoms(iatm) = wrk(qmmm_env_mm%mm_link_atoms(iatm))
     390             :             END DO
     391           8 :             qmmm_env_mm%mm_link_atoms = mm_link_atoms
     392           4 :             CPASSERT(ALL(qmmm_env_mm%mm_link_atoms /= 0))
     393           2 :             DEALLOCATE (mm_link_atoms)
     394             :             ! Update the qmmm_section: MM_INDEX,QM_INDEX of the LINK atom list
     395           2 :             qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
     396           2 :             CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
     397           2 :             CPASSERT(nlinks /= 0)
     398           6 :             DO ikind = 1, nlinks
     399           2 :                CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
     400           2 :                CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
     401           2 :                mm_index = wrk(mm_index)
     402           2 :                qm_index = wrk(qm_index)
     403           2 :                CALL section_vals_val_set(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
     404           6 :                CALL section_vals_val_set(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
     405             :             END DO
     406             :          END IF
     407             :       END IF
     408             :       !
     409             :       !LIST of ISOLATED atoms
     410             :       !
     411         314 :       generate_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE")
     412         314 :       isolated_section => section_vals_get_subs_vals(generate_section, "ISOLATED_ATOMS")
     413         314 :       CALL section_vals_get(isolated_section, explicit=explicit)
     414         314 :       IF (explicit) THEN
     415           4 :          CALL section_vals_val_get(isolated_section, "LIST", n_rep_val=n_rep)
     416          10 :          DO i = 1, n_rep
     417           6 :             CALL section_vals_val_get(isolated_section, "LIST", i_vals=tmp_v, i_rep_val=i)
     418          18 :             ALLOCATE (tmp_n(SIZE(tmp_v)))
     419          50 :             DO j = 1, SIZE(tmp_v)
     420          50 :                tmp_n(j) = wrk(tmp_v(j))
     421             :             END DO
     422          10 :             CALL section_vals_val_set(isolated_section, "LIST", i_vals_ptr=tmp_n, i_rep_val=i)
     423             :          END DO
     424             :       END IF
     425         314 :       DEALLOCATE (wrk)
     426             :       !DEALLOCATE all the temporary arrays needed for this routine
     427         314 :       DEALLOCATE (new_position)
     428         314 :       DEALLOCATE (tlabel_atmname)
     429         314 :       DEALLOCATE (tlabel_molname)
     430         314 :       DEALLOCATE (tlabel_resname)
     431         314 :       DEALLOCATE (telement)
     432         314 :       DEALLOCATE (tr)
     433         314 :       DEALLOCATE (tatm_charge)
     434         314 :       DEALLOCATE (tatm_mass)
     435         314 :       DEALLOCATE (molname)
     436         314 :       DEALLOCATE (mol_num)
     437             : 
     438             :       ! DEALLOCATE the bond structures in the connectivity
     439         314 :       DEALLOCATE (conn_info%bond_a)
     440         314 :       DEALLOCATE (conn_info%bond_b)
     441         314 :       IF (output_unit > 0) WRITE (output_unit, '(T2,"REORDER |  ")')
     442         314 :       CALL timestop(handle)
     443             :       CALL cp_print_key_finished_output(iw, logger, subsys_section, &
     444         314 :                                         "PRINT%TOPOLOGY_INFO/UTIL_INFO")
     445        1256 :    END SUBROUTINE topology_reorder_atoms
     446             : 
     447             : ! **************************************************************************************************
     448             : !> \brief Set up a SET of graph kind
     449             : !> \param graph_set ...
     450             : !> \param idim ...
     451             : !> \param ind ...
     452             : !> \param array2 ...
     453             : !> \param atom_bond_list ...
     454             : !> \param map_mol ...
     455             : !> \param atm_map1 ...
     456             : !> \param atm_map2 ...
     457             : !> \param atm_map3 ...
     458             : !> \author Teodoro Laino 10.2006
     459             : ! **************************************************************************************************
     460         326 :    SUBROUTINE setup_graph_set(graph_set, idim, ind, array2, atom_bond_list, map_mol, &
     461             :                               atm_map1, atm_map2, atm_map3)
     462             :       TYPE(graph_type), DIMENSION(:), POINTER            :: graph_set
     463             :       INTEGER, INTENT(IN)                                :: idim, ind
     464             :       INTEGER, DIMENSION(:), POINTER                     :: array2
     465             :       TYPE(array1_list_type), DIMENSION(:), INTENT(IN)   :: atom_bond_list
     466             :       INTEGER, DIMENSION(:, :), POINTER                  :: map_mol
     467             :       INTEGER, DIMENSION(:), POINTER                     :: atm_map1, atm_map2, atm_map3
     468             : 
     469             :       INTEGER                                            :: ldim
     470         326 :       TYPE(graph_type), DIMENSION(:), POINTER            :: tmp_graph_set
     471             : 
     472         326 :       ldim = 0
     473         326 :       NULLIFY (tmp_graph_set)
     474         326 :       IF (ASSOCIATED(graph_set)) THEN
     475          12 :          ldim = SIZE(graph_set)
     476          12 :          CPASSERT(ldim + 1 == idim)
     477             :          MARK_USED(idim)
     478             :          NULLIFY (tmp_graph_set)
     479          12 :          CALL allocate_graph_set(graph_set, tmp_graph_set)
     480             :       END IF
     481         326 :       CALL allocate_graph_set(tmp_graph_set, graph_set, ldim, ldim + 1)
     482             :       CALL setup_graph(ind, graph_set(ldim + 1)%graph, array2, &
     483         326 :                        atom_bond_list, map_mol, atm_map1, atm_map2, atm_map3)
     484             : 
     485         326 :    END SUBROUTINE setup_graph_set
     486             : 
     487             : ! **************************************************************************************************
     488             : !> \brief Allocate a new graph_set deallocating an old one..
     489             : !> \param graph_set_in ...
     490             : !> \param graph_set_out ...
     491             : !> \param ldim_in ...
     492             : !> \param ldim_out ...
     493             : !> \author Teodoro Laino 10.2006
     494             : ! **************************************************************************************************
     495         338 :    SUBROUTINE allocate_graph_set(graph_set_in, graph_set_out, ldim_in, ldim_out)
     496             :       TYPE(graph_type), DIMENSION(:), POINTER            :: graph_set_in, graph_set_out
     497             :       INTEGER, INTENT(IN), OPTIONAL                      :: ldim_in, ldim_out
     498             : 
     499             :       INTEGER                                            :: b_dim, i, j, mydim_in, mydim_out, v_dim
     500             : 
     501         338 :       CPASSERT(.NOT. ASSOCIATED(graph_set_out))
     502         338 :       mydim_in = 0
     503         338 :       mydim_out = 0
     504         338 :       IF (ASSOCIATED(graph_set_in)) THEN
     505          24 :          mydim_in = SIZE(graph_set_in)
     506          24 :          mydim_out = SIZE(graph_set_in)
     507             :       END IF
     508         338 :       IF (PRESENT(ldim_in)) mydim_in = ldim_in
     509         338 :       IF (PRESENT(ldim_out)) mydim_out = ldim_out
     510        1014 :       ALLOCATE (graph_set_out(mydim_out))
     511         696 :       DO i = 1, mydim_out
     512         696 :          NULLIFY (graph_set_out(i)%graph)
     513             :       END DO
     514             :       ! Copy graph structure into the temporary array
     515         370 :       DO i = 1, mydim_in
     516          32 :          v_dim = SIZE(graph_set_in(i)%graph)
     517          96 :          ALLOCATE (graph_set_out(i)%graph(v_dim))
     518         268 :          DO j = 1, v_dim
     519         236 :             graph_set_out(i)%graph(j)%kind = graph_set_in(i)%graph(j)%kind
     520         236 :             b_dim = SIZE(graph_set_in(i)%graph(j)%bonds)
     521         700 :             ALLOCATE (graph_set_out(i)%graph(j)%bonds(b_dim))
     522        1304 :             graph_set_out(i)%graph(j)%bonds = graph_set_in(i)%graph(j)%bonds
     523         268 :             DEALLOCATE (graph_set_in(i)%graph(j)%bonds)
     524             :          END DO
     525         370 :          DEALLOCATE (graph_set_in(i)%graph)
     526             :       END DO
     527         338 :       IF (ASSOCIATED(graph_set_in)) THEN
     528          24 :          DEALLOCATE (graph_set_in)
     529             :       END IF
     530             : 
     531         338 :    END SUBROUTINE allocate_graph_set
     532             : 
     533             : ! **************************************************************************************************
     534             : !> \brief Set up a graph kind
     535             : !> \param ind ...
     536             : !> \param graph ...
     537             : !> \param array2 ...
     538             : !> \param atom_bond_list ...
     539             : !> \param map_mol ...
     540             : !> \param atm_map1 ...
     541             : !> \param atm_map2 ...
     542             : !> \param atm_map3 ...
     543             : !> \author Teodoro Laino 09.2006
     544             : ! **************************************************************************************************
     545        2840 :    SUBROUTINE setup_graph(ind, graph, array2, atom_bond_list, map_mol, &
     546             :                           atm_map1, atm_map2, atm_map3)
     547             :       INTEGER, INTENT(IN)                                :: ind
     548             :       TYPE(vertex), DIMENSION(:), POINTER                :: graph
     549             :       INTEGER, DIMENSION(:), POINTER                     :: array2
     550             :       TYPE(array1_list_type), DIMENSION(:), INTENT(IN)   :: atom_bond_list
     551             :       INTEGER, DIMENSION(:, :), POINTER                  :: map_mol
     552             :       INTEGER, DIMENSION(:), POINTER                     :: atm_map1, atm_map2
     553             :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: atm_map3
     554             : 
     555             :       INTEGER                                            :: i, idim, ifirst, ilast, j, nbonds, &
     556             :                                                             nelement
     557             : 
     558        2840 :       IF (PRESENT(atm_map3)) THEN
     559        1420 :          CPASSERT(.NOT. ASSOCIATED(atm_map3))
     560             :       END IF
     561        2840 :       CPASSERT(.NOT. ASSOCIATED(graph))
     562             :       ! Setup reference graph
     563        2840 :       idim = 0
     564        2840 :       ifirst = map_mol(1, ind)
     565        2840 :       ilast = map_mol(2, ind)
     566        2840 :       nelement = ilast - ifirst + 1
     567        8520 :       ALLOCATE (graph(nelement))
     568        2840 :       IF (PRESENT(atm_map3)) THEN
     569        4260 :          ALLOCATE (atm_map3(nelement))
     570             :       END IF
     571       19956 :       DO i = ifirst, ilast
     572       17116 :          idim = idim + 1
     573       17116 :          graph(idim)%kind = array2(atm_map1(i))
     574       17116 :          nbonds = SIZE(atom_bond_list(atm_map1(i))%array1)
     575       51260 :          ALLOCATE (graph(idim)%bonds(nbonds))
     576       49220 :          DO j = 1, nbonds
     577       49220 :             graph(idim)%bonds(j) = atm_map2(atom_bond_list(atm_map1(i))%array1(j))
     578             :          END DO
     579       19956 :          IF (PRESENT(atm_map3)) THEN
     580        8558 :             atm_map3(idim) = atm_map1(i)
     581             :          END IF
     582             :       END DO
     583             : 
     584        2840 :    END SUBROUTINE setup_graph
     585             : 
     586             : ! **************************************************************************************************
     587             : !> \brief Order arrays of lists
     588             : !> \param Ilist1 ...
     589             : !> \param Ilist2 ...
     590             : !> \param Ilist3 ...
     591             : !> \param Ilist4 ...
     592             : !> \param nsize ...
     593             : !> \param ndim ...
     594             : !> \author Teodoro Laino 09.2006
     595             : ! **************************************************************************************************
     596      259080 :    RECURSIVE SUBROUTINE reorder_list_array(Ilist1, Ilist2, Ilist3, Ilist4, nsize, ndim)
     597             :       INTEGER, DIMENSION(:), POINTER                     :: Ilist1
     598             :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: Ilist2, Ilist3, Ilist4
     599             :       INTEGER, INTENT(IN)                                :: nsize, ndim
     600             : 
     601             :       INTEGER                                            :: i, iend, istart, ldim
     602      259080 :       INTEGER, DIMENSION(:), POINTER                     :: tmp, wrk
     603             : 
     604           0 :       CPASSERT(nsize > 0)
     605      754274 :       ALLOCATE (wrk(Ndim))
     606      259080 :       CALL sort(Ilist1, Ndim, wrk)
     607      259080 :       IF (nsize /= 1) THEN
     608      323576 :          ALLOCATE (tmp(Ndim))
     609     1051786 :          tmp = Ilist2(1:Ndim)
     610      525893 :          DO i = 1, Ndim
     611      525893 :             Ilist2(i) = tmp(wrk(i))
     612             :          END DO
     613       19604 :          SELECT CASE (nsize)
     614             :          CASE (3)
     615      292716 :             tmp = Ilist3(1:Ndim)
     616      146358 :             DO i = 1, Ndim
     617      146358 :                Ilist3(i) = tmp(wrk(i))
     618             :             END DO
     619             :          CASE (4)
     620       96172 :             tmp = Ilist3(1:Ndim)
     621       48086 :             DO i = 1, Ndim
     622       48086 :                Ilist3(i) = tmp(wrk(i))
     623             :             END DO
     624       96172 :             tmp = Ilist4(1:Ndim)
     625      160823 :             DO i = 1, Ndim
     626       48086 :                Ilist4(i) = tmp(wrk(i))
     627             :             END DO
     628             :          END SELECT
     629      112737 :          DEALLOCATE (tmp)
     630      112737 :          istart = 1
     631      525893 :          DO i = 1, Ndim
     632      525893 :             IF (Ilist1(i) /= Ilist1(istart)) THEN
     633      133628 :                iend = i - 1
     634      133628 :                ldim = iend - istart + 1
     635             :                CALL reorder_list_array_low(Ilist2, Ilist3, Ilist4, nsize, &
     636      309942 :                                            ldim, istart, iend)
     637      133628 :                istart = i
     638             :             END IF
     639             :          END DO
     640             :          ! Last term to sort
     641      112737 :          iend = Ndim
     642      112737 :          ldim = iend - istart + 1
     643             :          CALL reorder_list_array_low(Ilist2, Ilist3, Ilist4, nsize, &
     644      311343 :                                      ldim, istart, iend)
     645             :       END IF
     646      259080 :       DEALLOCATE (wrk)
     647      259080 :    END SUBROUTINE reorder_list_array
     648             : 
     649             : ! **************************************************************************************************
     650             : !> \brief Low level routine for ordering arrays of lists
     651             : !> \param Ilist2 ...
     652             : !> \param Ilist3 ...
     653             : !> \param Ilist4 ...
     654             : !> \param nsize ...
     655             : !> \param ldim ...
     656             : !> \param istart ...
     657             : !> \param iend ...
     658             : !> \author Teodoro Laino 09.2006
     659             : ! **************************************************************************************************
     660      246365 :    RECURSIVE SUBROUTINE reorder_list_array_low(Ilist2, Ilist3, Ilist4, nsize, &
     661             :                                                ldim, istart, iend)
     662             :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: Ilist2, Ilist3, Ilist4
     663             :       INTEGER, INTENT(IN)                                :: nsize, ldim, istart, iend
     664             : 
     665      246365 :       INTEGER, DIMENSION(:), POINTER                     :: tmp_2, tmp_3, tmp_4
     666             : 
     667      392708 :       SELECT CASE (nsize)
     668             :       CASE (2)
     669      430698 :          ALLOCATE (tmp_2(ldim))
     670      776582 :          tmp_2(:) = Ilist2(istart:iend)
     671      146343 :          CALL reorder_list_array(tmp_2, nsize=nsize - 1, ndim=ldim)
     672      776582 :          Ilist2(istart:iend) = tmp_2(:)
     673      146343 :          DEALLOCATE (tmp_2)
     674             :       CASE (3)
     675      243508 :          ALLOCATE (tmp_2(ldim))
     676      243508 :          ALLOCATE (tmp_3(ldim))
     677      417976 :          tmp_2(:) = Ilist2(istart:iend)
     678      417976 :          tmp_3(:) = Ilist3(istart:iend)
     679       82234 :          CALL reorder_list_array(tmp_2, tmp_3, nsize=nsize - 1, ndim=ldim)
     680      417976 :          Ilist2(istart:iend) = tmp_2(:)
     681      417976 :          Ilist3(istart:iend) = tmp_3(:)
     682       82234 :          DEALLOCATE (tmp_2)
     683       82234 :          DEALLOCATE (tmp_3)
     684             :       CASE (4)
     685       50254 :          ALLOCATE (tmp_2(ldim))
     686       50254 :          ALLOCATE (tmp_3(ldim))
     687       50254 :          ALLOCATE (tmp_4(ldim))
     688      124484 :          tmp_2(:) = Ilist2(istart:iend)
     689      124484 :          tmp_3(:) = Ilist3(istart:iend)
     690      124484 :          tmp_4(:) = Ilist4(istart:iend)
     691       17788 :          CALL reorder_list_array(tmp_2, tmp_3, tmp_4, nsize=nsize - 1, ndim=ldim)
     692      124484 :          Ilist2(istart:iend) = tmp_2(:)
     693      124484 :          Ilist3(istart:iend) = tmp_3(:)
     694      124484 :          Ilist4(istart:iend) = tmp_4(:)
     695       17788 :          DEALLOCATE (tmp_2)
     696       17788 :          DEALLOCATE (tmp_3)
     697       17788 :          DEALLOCATE (tmp_4)
     698             :       END SELECT
     699             : 
     700      246365 :    END SUBROUTINE reorder_list_array_low
     701             : 
     702             : ! **************************************************************************************************
     703             : !> \brief ...
     704             : !> \param icheck ...
     705             : !> \param bond_list ...
     706             : !> \param i ...
     707             : !> \param mol_natom ...
     708             : !> \param mol_map ...
     709             : !> \param my_mol ...
     710             : !> \author Teodoro Laino 09.2006
     711             : ! **************************************************************************************************
     712      205158 :    RECURSIVE SUBROUTINE give_back_molecule(icheck, bond_list, i, mol_natom, mol_map, my_mol)
     713             :       LOGICAL, DIMENSION(:), POINTER                     :: icheck
     714             :       TYPE(array1_list_type), DIMENSION(:), POINTER      :: bond_list
     715             :       INTEGER, INTENT(IN)                                :: i
     716             :       INTEGER, INTENT(INOUT)                             :: mol_natom
     717             :       INTEGER, DIMENSION(:), POINTER                     :: mol_map
     718             :       INTEGER, INTENT(IN)                                :: my_mol
     719             : 
     720             :       INTEGER                                            :: j, k
     721             : 
     722      205158 :       IF (mol_map(i) == my_mol) THEN
     723      205150 :          icheck(i) = .TRUE.
     724      419526 :          DO j = 1, SIZE(bond_list(i)%array1)
     725      214376 :             k = bond_list(i)%array1(j)
     726      214376 :             IF (icheck(k)) CYCLE
     727      106188 :             mol_natom = mol_natom + 1
     728      419526 :             CALL give_back_molecule(icheck, bond_list, k, mol_natom, mol_map, my_mol)
     729             :          END DO
     730             :       ELSE
     731             :          ! Do nothing means only that bonds were found between molecules
     732             :          ! as we defined them.. This could easily be a bond detected but not
     733             :          ! physically present.. so just skip this part and go on..
     734             :       END IF
     735      205158 :    END SUBROUTINE give_back_molecule
     736             : 
     737             : ! **************************************************************************************************
     738             : !> \brief gives back a mapping of molecules.. icheck needs to be initialized with -1
     739             : !> \param icheck ...
     740             : !> \param bond_list ...
     741             : !> \param i ...
     742             : !> \param my_mol ...
     743             : !> \author Teodoro Laino 04.2007 - Zurich University
     744             : ! **************************************************************************************************
     745        6706 :    RECURSIVE SUBROUTINE tag_molecule(icheck, bond_list, i, my_mol)
     746             :       INTEGER, DIMENSION(:), POINTER                     :: icheck
     747             :       TYPE(array1_list_type), DIMENSION(:), POINTER      :: bond_list
     748             :       INTEGER, INTENT(IN)                                :: i, my_mol
     749             : 
     750             :       INTEGER                                            :: j, k
     751             : 
     752        6706 :       icheck(i) = my_mol
     753       17898 :       DO j = 1, SIZE(bond_list(i)%array1)
     754       11192 :          k = bond_list(i)%array1(j)
     755       11192 :          IF (k <= i) CYCLE
     756       17898 :          CALL tag_molecule(icheck, bond_list, k, my_mol)
     757             :       END DO
     758             : 
     759        6706 :    END SUBROUTINE tag_molecule
     760             : 
     761             : ! **************************************************************************************************
     762             : !> \brief Given two lists of corresponding element a complex type is built in
     763             : !>      which each element of the type has a list of mapping elements
     764             : !> \param work ...
     765             : !> \param list1 ...
     766             : !> \param list2 ...
     767             : !> \param N ...
     768             : !> \author Teodoro Laino 08.2006
     769             : ! **************************************************************************************************
     770       66485 :    SUBROUTINE reorder_structure1d(work, list1, list2, N)
     771             :       TYPE(array1_list_type), DIMENSION(:), &
     772             :          INTENT(INOUT)                                   :: work
     773             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: list1, list2
     774             :       INTEGER, INTENT(IN)                                :: N
     775             : 
     776             :       INTEGER                                            :: I, index1, index2, Nsize
     777       66485 :       INTEGER, DIMENSION(:), POINTER                     :: wrk_tmp
     778             : 
     779     3367771 :       DO I = 1, N
     780     3301286 :          index1 = list1(I)
     781     3301286 :          index2 = list2(I)
     782             : 
     783     3301286 :          wrk_tmp => work(index1)%array1
     784     3301286 :          Nsize = SIZE(wrk_tmp)
     785     9903858 :          ALLOCATE (work(index1)%array1(Nsize + 1))
     786     6453143 :          work(index1)%array1(1:Nsize) = wrk_tmp
     787     3301286 :          work(index1)%array1(Nsize + 1) = index2
     788     3301286 :          DEALLOCATE (wrk_tmp)
     789             : 
     790     3301286 :          wrk_tmp => work(index2)%array1
     791     3301286 :          Nsize = SIZE(wrk_tmp)
     792     9903858 :          ALLOCATE (work(index2)%array1(Nsize + 1))
     793     4518671 :          work(index2)%array1(1:Nsize) = wrk_tmp
     794     3301286 :          work(index2)%array1(Nsize + 1) = index1
     795     3367771 :          DEALLOCATE (wrk_tmp)
     796             :       END DO
     797             : 
     798       66485 :    END SUBROUTINE reorder_structure1d
     799             : 
     800             : ! **************************************************************************************************
     801             : !> \brief Given two lists of corresponding element a complex type is built in
     802             : !>      which each element of the type has a list of mapping elements
     803             : !> \param work ...
     804             : !> \param list1 ...
     805             : !> \param list2 ...
     806             : !> \param list3 ...
     807             : !> \param N ...
     808             : !> \author Teodoro Laino 09.2006
     809             : ! **************************************************************************************************
     810        7265 :    SUBROUTINE reorder_structure2d(work, list1, list2, list3, N)
     811             :       TYPE(array2_list_type), DIMENSION(:), &
     812             :          INTENT(INOUT)                                   :: work
     813             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: list1, list2, list3
     814             :       INTEGER, INTENT(IN)                                :: N
     815             : 
     816             :       INTEGER                                            :: I, index1, index2, index3, Nsize
     817        7265 :       INTEGER, DIMENSION(:), POINTER                     :: wrk_tmp
     818             : 
     819     1138238 :       DO I = 1, N
     820     1130973 :          index1 = list1(I)
     821     1130973 :          index2 = list2(I)
     822     1130973 :          index3 = list3(I)
     823             : 
     824     1130973 :          wrk_tmp => work(index1)%array1
     825     1130973 :          Nsize = SIZE(wrk_tmp)
     826     3392919 :          ALLOCATE (work(index1)%array1(Nsize + 1))
     827    23831301 :          work(index1)%array1(1:Nsize) = wrk_tmp
     828     1130973 :          work(index1)%array1(Nsize + 1) = index2
     829     1130973 :          DEALLOCATE (wrk_tmp)
     830             : 
     831     1130973 :          wrk_tmp => work(index2)%array1
     832     1130973 :          Nsize = SIZE(wrk_tmp)
     833     3392919 :          ALLOCATE (work(index2)%array1(Nsize + 1))
     834    15715757 :          work(index2)%array1(1:Nsize) = wrk_tmp
     835     1130973 :          work(index2)%array1(Nsize + 1) = index1
     836     1130973 :          DEALLOCATE (wrk_tmp)
     837             : 
     838     1130973 :          wrk_tmp => work(index1)%array2
     839     1130973 :          Nsize = SIZE(wrk_tmp)
     840     3392919 :          ALLOCATE (work(index1)%array2(Nsize + 1))
     841    23831301 :          work(index1)%array2(1:Nsize) = wrk_tmp
     842     1130973 :          work(index1)%array2(Nsize + 1) = index3
     843     1130973 :          DEALLOCATE (wrk_tmp)
     844             : 
     845     1130973 :          wrk_tmp => work(index2)%array2
     846     1130973 :          Nsize = SIZE(wrk_tmp)
     847     3392919 :          ALLOCATE (work(index2)%array2(Nsize + 1))
     848    15715757 :          work(index2)%array2(1:Nsize) = wrk_tmp
     849     1130973 :          work(index2)%array2(Nsize + 1) = -index3
     850     1138238 :          DEALLOCATE (wrk_tmp)
     851             :       END DO
     852             : 
     853        7265 :    END SUBROUTINE reorder_structure2d
     854             : 
     855             : ! **************************************************************************************************
     856             : !> \brief each atom will be assigned a molecule number based on bonded fragments
     857             : !>      The array mol_info should be initialized with -1 before calling the
     858             : !>      find_molecule routine
     859             : !> \param atom_bond_list ...
     860             : !> \param mol_info ...
     861             : !> \param mol_name ...
     862             : !> \author Joost 05.2006
     863             : ! **************************************************************************************************
     864        9684 :    SUBROUTINE find_molecule(atom_bond_list, mol_info, mol_name)
     865             :       TYPE(array1_list_type), DIMENSION(:), INTENT(IN)   :: atom_bond_list
     866             :       INTEGER, DIMENSION(:), POINTER                     :: mol_info, mol_name
     867             : 
     868             :       INTEGER                                            :: I, my_mol_name, N, nmol
     869             : 
     870        9684 :       N = SIZE(atom_bond_list)
     871        9684 :       nmol = 0
     872      768707 :       DO I = 1, N
     873      768707 :          IF (mol_info(I) == -1) THEN
     874      326834 :             nmol = nmol + 1
     875      326834 :             my_mol_name = mol_name(I)
     876             :             CALL spread_mol(atom_bond_list, mol_info, i, nmol, my_mol_name, &
     877      326834 :                             mol_name)
     878             :          END IF
     879             :       END DO
     880        9684 :    END SUBROUTINE find_molecule
     881             : 
     882             : ! **************************************************************************************************
     883             : !> \brief spreads the molnumber over the bonded list
     884             : !> \param atom_bond_list ...
     885             : !> \param mol_info ...
     886             : !> \param iatom ...
     887             : !> \param imol ...
     888             : !> \param my_mol_name ...
     889             : !> \param mol_name ...
     890             : !> \author Joost 05.2006
     891             : ! **************************************************************************************************
     892      759023 :    RECURSIVE SUBROUTINE spread_mol(atom_bond_list, mol_info, iatom, imol, &
     893             :                                    my_mol_name, mol_name)
     894             :       TYPE(array1_list_type), DIMENSION(:), INTENT(IN)   :: atom_bond_list
     895             :       INTEGER, DIMENSION(:), POINTER                     :: mol_info
     896             :       INTEGER, INTENT(IN)                                :: iatom, imol, my_mol_name
     897             :       INTEGER, DIMENSION(:), POINTER                     :: mol_name
     898             : 
     899             :       INTEGER                                            :: atom_b, i
     900             : 
     901      759023 :       mol_info(iatom) = imol
     902     1791349 :       DO I = 1, SIZE(atom_bond_list(iatom)%array1)
     903     1032326 :          atom_b = atom_bond_list(iatom)%array1(I)
     904             :          ! In this way we're really sure that all atoms belong to the same
     905             :          ! molecule. This should correct possible errors in the generation of
     906             :          ! the bond list..
     907     1791349 :          IF (mol_info(atom_b) == -1 .AND. my_mol_name == mol_name(atom_b)) THEN
     908      432189 :             CALL spread_mol(atom_bond_list, mol_info, atom_b, imol, my_mol_name, mol_name)
     909      432189 :             IF (mol_info(atom_b) /= imol) CPABORT("internal error")
     910             :          END IF
     911             :       END DO
     912      759023 :    END SUBROUTINE spread_mol
     913             : 
     914             : ! **************************************************************************************************
     915             : !> \brief Use info from periodic table and set atm_mass
     916             : !> \param topology ...
     917             : !> \param subsys_section ...
     918             : ! **************************************************************************************************
     919        8893 :    SUBROUTINE topology_set_atm_mass(topology, subsys_section)
     920             :       TYPE(topology_parameters_type), INTENT(INOUT)      :: topology
     921             :       TYPE(section_vals_type), POINTER                   :: subsys_section
     922             : 
     923             :       CHARACTER(len=*), PARAMETER :: routineN = 'topology_set_atm_mass'
     924             : 
     925             :       CHARACTER(LEN=2)                                   :: upper_sym_1
     926             :       CHARACTER(LEN=default_string_length)               :: atmname_upper
     927             :       CHARACTER(LEN=default_string_length), &
     928        8893 :          DIMENSION(:), POINTER                           :: keyword
     929             :       INTEGER                                            :: handle, i, i_rep, iatom, ielem_found, &
     930             :                                                             iw, n_rep, natom
     931             :       LOGICAL                                            :: user_defined
     932        8893 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: mass
     933             :       TYPE(atom_info_type), POINTER                      :: atom_info
     934             :       TYPE(cp_logger_type), POINTER                      :: logger
     935             :       TYPE(section_vals_type), POINTER                   :: kind_section
     936             : 
     937        8893 :       NULLIFY (logger)
     938       17786 :       logger => cp_get_default_logger()
     939             :       iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
     940        8893 :                                 extension=".subsysLog")
     941        8893 :       CALL timeset(routineN, handle)
     942             : 
     943        8893 :       atom_info => topology%atom_info
     944        8893 :       natom = topology%natoms
     945             : 
     946             :       ! Available external info
     947        8893 :       kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
     948        8893 :       CALL section_vals_get(kind_section, n_repetition=n_rep)
     949       23657 :       ALLOCATE (keyword(n_rep))
     950       23657 :       ALLOCATE (mass(n_rep))
     951       20704 :       mass = HUGE(0.0_dp)
     952       20704 :       DO i_rep = 1, n_rep
     953             :          CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
     954       11811 :                                    c_val=keyword(i_rep), i_rep_section=i_rep)
     955       11811 :          CALL uppercase(keyword(i_rep))
     956             :          CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
     957       11811 :                                    keyword_name="MASS", n_rep_val=i)
     958       11811 :          IF (i > 0) CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
     959       20800 :                                               keyword_name="MASS", r_val=mass(i_rep))
     960             :       END DO
     961             :       !
     962      283931 :       DO iatom = 1, natom
     963             :          !If we reach this point then we've definitely identified the element..
     964             :          !Let's look if an external mass has been defined..
     965      275038 :          user_defined = .FALSE.
     966      386454 :          DO i = 1, SIZE(keyword)
     967      112010 :             atmname_upper = id2str(atom_info%id_atmname(iatom))
     968      112010 :             CALL uppercase(atmname_upper)
     969      386454 :             IF (TRIM(atmname_upper) == TRIM(keyword(i)) .AND. mass(i) /= HUGE(0.0_dp)) THEN
     970         594 :                atom_info%atm_mass(iatom) = mass(i)
     971             :                user_defined = .TRUE.
     972             :                EXIT
     973             :             END IF
     974             :          END DO
     975             :          ! If name didn't match let's try with the element
     976      274444 :          IF (.NOT. user_defined) THEN
     977      274444 :             upper_sym_1 = id2str(atom_info%id_element(iatom))
     978      274444 :             CALL get_ptable_info(symbol=upper_sym_1, ielement=ielem_found, amass=atom_info%atm_mass(iatom))
     979             :          END IF
     980      276442 :          IF (iw > 0) WRITE (iw, '(7X,A,A5,A,F12.5)') "In topology_set_atm_mass :: element = ", &
     981       11701 :             id2str(atom_info%id_element(iatom)), " a_mass ", atom_info%atm_mass(iatom)
     982             :       END DO
     983        8893 :       DEALLOCATE (keyword)
     984        8893 :       DEALLOCATE (mass)
     985             : 
     986        8893 :       CALL timestop(handle)
     987             :       CALL cp_print_key_finished_output(iw, logger, subsys_section, &
     988        8893 :                                         "PRINT%TOPOLOGY_INFO/UTIL_INFO")
     989             : 
     990       26679 :    END SUBROUTINE topology_set_atm_mass
     991             : 
     992             : ! **************************************************************************************************
     993             : !> \brief Check and verify that all molecules of the same kind are bonded the same
     994             : !> \param topology ...
     995             : !> \param subsys_section ...
     996             : ! **************************************************************************************************
     997        8842 :    SUBROUTINE topology_molecules_check(topology, subsys_section)
     998             :       TYPE(topology_parameters_type), INTENT(INOUT)      :: topology
     999             :       TYPE(section_vals_type), POINTER                   :: subsys_section
    1000             : 
    1001             :       CHARACTER(len=*), PARAMETER :: routineN = 'topology_molecules_check'
    1002             : 
    1003             :       INTEGER                                            :: counter, first, first_loc, handle, i, &
    1004             :                                                             iatom, iw, k, loc_counter, mol_num, &
    1005             :                                                             mol_typ, n, natom
    1006             :       LOGICAL                                            :: icheck_num, icheck_typ
    1007        8842 :       TYPE(array1_list_type), ALLOCATABLE, DIMENSION(:)  :: atom_bond_list
    1008             :       TYPE(atom_info_type), POINTER                      :: atom_info
    1009             :       TYPE(connectivity_info_type), POINTER              :: conn_info
    1010             :       TYPE(cp_logger_type), POINTER                      :: logger
    1011             : 
    1012        8842 :       NULLIFY (logger)
    1013        8842 :       logger => cp_get_default_logger()
    1014             :       iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
    1015        8842 :                                 extension=".subsysLog")
    1016        8842 :       CALL timeset(routineN, handle)
    1017             : 
    1018        8842 :       atom_info => topology%atom_info
    1019        8842 :       conn_info => topology%conn_info
    1020        8842 :       natom = topology%natoms
    1021             : 
    1022        8868 :       IF (iw > 0) WRITE (iw, '(A)') "Start of Molecule_Check", &
    1023          52 :          "  Checking consistency between the generated molecules"
    1024             : 
    1025      773327 :       ALLOCATE (atom_bond_list(natom))
    1026      755643 :       DO I = 1, natom
    1027      755643 :          ALLOCATE (atom_bond_list(I)%array1(0))
    1028             :       END DO
    1029        8842 :       N = 0
    1030        8842 :       IF (ASSOCIATED(conn_info%bond_a)) N = SIZE(conn_info%bond_a)
    1031        8842 :       CALL reorder_structure(atom_bond_list, conn_info%bond_a, conn_info%bond_b, N)
    1032             : 
    1033        8842 :       mol_typ = atom_info%map_mol_typ(1)
    1034        8842 :       mol_num = atom_info%map_mol_num(1)
    1035        8842 :       counter = 1
    1036        8842 :       loc_counter = 1
    1037        8842 :       first = 1
    1038        8842 :       first_loc = 1
    1039      746801 :       DO iatom = 2, natom
    1040      737959 :          icheck_num = (atom_info%map_mol_num(iatom) == mol_num)
    1041      737959 :          icheck_typ = (atom_info%map_mol_typ(iatom) == mol_typ)
    1042      737959 :          IF ((icheck_typ .AND. (.NOT. icheck_num)) .OR. (.NOT. icheck_typ)) THEN
    1043             :             !-----------------------------------------------------------------------------
    1044             :             !-----------------------------------------------------------------------------
    1045             :             ! 1. Check each molecule have the same number of atoms
    1046             :             !-----------------------------------------------------------------------------
    1047      284986 :             IF (counter /= loc_counter) THEN
    1048             :                CALL cp_abort(__LOCATION__, &
    1049             :                              "different number of atoms for same molecule kind"// &
    1050             :                              " molecule type  = "//cp_to_string(mol_typ)// &
    1051             :                              " molecule number= "//cp_to_string(mol_num)// &
    1052             :                              " expected number of atoms="//cp_to_string(counter)//" found="// &
    1053           0 :                              cp_to_string(loc_counter))
    1054             :             END IF
    1055             :          END IF
    1056      737959 :          IF (.NOT. icheck_typ) THEN
    1057      116749 :             first = iatom
    1058      116749 :             first_loc = iatom
    1059      116749 :             counter = 1
    1060      116749 :             loc_counter = 1
    1061      116749 :             mol_typ = atom_info%map_mol_typ(iatom)
    1062             :          END IF
    1063      737959 :          IF (icheck_num) THEN
    1064      569030 :             IF (icheck_typ) loc_counter = loc_counter + 1
    1065             :             !-----------------------------------------------------------------------------
    1066             :             !-----------------------------------------------------------------------------
    1067             :             ! 2. Check that each molecule has the same atom sequences
    1068             :             !-----------------------------------------------------------------------------
    1069      569030 :             IF (atom_info%id_atmname(iatom) /= &
    1070             :                 atom_info%id_atmname(first + loc_counter - 1)) THEN
    1071             :                CALL cp_abort(__LOCATION__, &
    1072             :                              "different atom name for same molecule kind"// &
    1073             :                              " atom number    = "//cp_to_string(iatom)// &
    1074             :                              " molecule type  = "//cp_to_string(mol_typ)// &
    1075             :                              " molecule number= "//cp_to_string(mol_num)// &
    1076             :                              " expected atom name="//TRIM(id2str(atom_info%id_atmname(first + loc_counter - 1)))// &
    1077           0 :                              " found="//TRIM(id2str(atom_info%id_atmname(iatom))))
    1078             :             END IF
    1079             :             !-----------------------------------------------------------------------------
    1080             :             !-----------------------------------------------------------------------------
    1081             :             ! 3. Check that each molecule have the same bond sequences
    1082             :             !-----------------------------------------------------------------------------
    1083      569030 :             IF (SIZE(atom_bond_list(iatom)%array1) /= SIZE(atom_bond_list(first + loc_counter - 1)%array1)) THEN
    1084             :                CALL cp_abort(__LOCATION__, &
    1085             :                              "different number of bonds for same molecule kind"// &
    1086             :                              " molecule type  = "//cp_to_string(mol_typ)// &
    1087             :                              " molecule number= "//cp_to_string(mol_num)// &
    1088             :                              " expected bonds="// &
    1089             :                              cp_to_string(SIZE(atom_bond_list(first + loc_counter - 1)%array1))//" - "// &
    1090             :                              cp_to_string(SIZE(atom_bond_list(iatom)%array1))// &
    1091           0 :                              " NOT FOUND! Check the connectivity of your system.")
    1092             :             END IF
    1093             : 
    1094     1276498 :             DO k = 1, SIZE(atom_bond_list(iatom)%array1)
    1095     1022713 :                IF (ALL(atom_bond_list(first + loc_counter - 1)%array1 - first /= &
    1096      569030 :                        atom_bond_list(iatom)%array1(k) - first_loc)) THEN
    1097             :                   CALL cp_abort(__LOCATION__, &
    1098             :                                 "different sequence of bonds for same molecule kind"// &
    1099             :                                 " molecule type  = "//cp_to_string(mol_typ)// &
    1100             :                                 " molecule number= "//cp_to_string(mol_num)// &
    1101           0 :                                 " NOT FOUND! Check the connectivity of your system.")
    1102             :                END IF
    1103             :             END DO
    1104             :          ELSE
    1105      168929 :             mol_num = atom_info%map_mol_num(iatom)
    1106      168929 :             loc_counter = 1
    1107      168929 :             first_loc = iatom
    1108             :          END IF
    1109      746801 :          IF (mol_num == 1 .AND. icheck_typ) counter = counter + 1
    1110             :       END DO
    1111        8842 :       IF (iw > 0) WRITE (iw, '(A)') "End of Molecule_Check"
    1112             : 
    1113      755643 :       DO I = 1, natom
    1114      755643 :          DEALLOCATE (atom_bond_list(I)%array1)
    1115             :       END DO
    1116        8842 :       DEALLOCATE (atom_bond_list)
    1117        8842 :       CALL timestop(handle)
    1118             :       CALL cp_print_key_finished_output(iw, logger, subsys_section, &
    1119        8842 :                                         "PRINT%TOPOLOGY_INFO/UTIL_INFO")
    1120             : 
    1121       17684 :    END SUBROUTINE topology_molecules_check
    1122             : 
    1123             : ! **************************************************************************************************
    1124             : !> \brief Check and returns the ELEMENT label
    1125             : !> \param element_in ...
    1126             : !> \param atom_name_in ...
    1127             : !> \param element_out ...
    1128             : !> \param subsys_section ...
    1129             : !> \param use_mm_map_first ...
    1130             : !> \par History
    1131             : !>      12.2005 created [teo]
    1132             : !> \author Teodoro Laino
    1133             : ! **************************************************************************************************
    1134       51062 :    SUBROUTINE check_subsys_element(element_in, atom_name_in, element_out, subsys_section, use_mm_map_first)
    1135             :       CHARACTER(len=*), INTENT(IN)                       :: element_in, atom_name_in
    1136             :       CHARACTER(len=default_string_length), INTENT(OUT)  :: element_out
    1137             :       TYPE(section_vals_type), POINTER                   :: subsys_section
    1138             :       LOGICAL                                            :: use_mm_map_first
    1139             : 
    1140             :       CHARACTER(len=default_string_length)               :: atom_name, element_symbol, keyword
    1141             :       INTEGER                                            :: i, i_rep, n_rep
    1142             :       LOGICAL                                            :: defined_kind_section, found, in_ptable
    1143             :       TYPE(section_vals_type), POINTER                   :: kind_section
    1144             : 
    1145       25531 :       found = .FALSE.
    1146       25531 :       element_symbol = element_in
    1147       25531 :       atom_name = atom_name_in
    1148       25531 :       element_out = ""
    1149       25531 :       defined_kind_section = .FALSE.
    1150             : 
    1151             :       ! First check if a KIND section is overriding the element
    1152             :       ! definition
    1153       25531 :       CALL uppercase(atom_name)
    1154       25531 :       kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
    1155       25531 :       CALL section_vals_get(kind_section, n_repetition=n_rep)
    1156       78378 :       DO i_rep = 1, n_rep
    1157             :          CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
    1158       54071 :                                    c_val=keyword, i_rep_section=i_rep)
    1159       54071 :          CALL uppercase(keyword)
    1160       78378 :          IF (TRIM(keyword) == TRIM(atom_name)) THEN
    1161             :             CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
    1162       14319 :                                       keyword_name="ELEMENT", n_rep_val=i)
    1163       14319 :             IF (i > 0) THEN
    1164             :                CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
    1165        1224 :                                          keyword_name="ELEMENT", c_val=element_symbol)
    1166             :                defined_kind_section = .TRUE.
    1167             :                EXIT
    1168             :             ELSE
    1169       13095 :                element_symbol = element_in
    1170             :                defined_kind_section = .TRUE.
    1171             :             END IF
    1172             :          END IF
    1173             :       END DO
    1174             : 
    1175             :       ! Let's check the validity of the element so far stored..
    1176             :       ! if we are not having a connectivity file, we must first match against the ptable.
    1177             :       ! this helps to resolve Ca/CA (calcium and Calpha) or Cn/CN7 (Coppernicum (112) CN) conflicts
    1178             :       ! so, in the presence of connectivity CA will be 'C', while in the absence of connectivity CA will be 'Ca'
    1179       24307 :       IF (defined_kind_section .OR. .NOT. use_mm_map_first) THEN
    1180             :          ! lengths larger than 2 should not match, because 'trailing' characters are ignored.
    1181       21055 :          in_ptable = .FALSE.
    1182       21055 :          IF (LEN_TRIM(element_symbol) <= 2) CALL get_ptable_info(element_symbol, found=in_ptable)
    1183       21055 :          IF (in_ptable) THEN
    1184       20961 :             element_out = TRIM(element_symbol)
    1185       20961 :             found = .TRUE.
    1186             :          END IF
    1187             :       END IF
    1188             : 
    1189             :       ! This is clearly a user error
    1190       25531 :       IF (.NOT. found .AND. defined_kind_section) &
    1191             :          CALL cp_abort(__LOCATION__, "Element <"//TRIM(element_symbol)// &
    1192             :                        "> provided for KIND <"//TRIM(atom_name_in)//"> "// &
    1193             :                        "which cannot be mapped with any standard element label. Please correct your "// &
    1194           0 :                        "input file!")
    1195             : 
    1196             :       ! Last chance.. are these atom_kinds of AMBER or CHARMM or GROMOS FF ?
    1197       25531 :       CALL uppercase(element_symbol)
    1198       25531 :       IF ((.NOT. found) .AND. (ASSOCIATED(amber_map))) THEN
    1199             :          ! First we go through the AMBER library
    1200      145670 :          DO i = 1, SIZE(amber_map%kind)
    1201      145670 :             IF (element_symbol == amber_map%kind(i)) THEN
    1202             :                found = .TRUE.
    1203             :                EXIT
    1204             :             END IF
    1205             :          END DO
    1206        4570 :          IF (found) THEN
    1207        4134 :             element_out = amber_map%element(i)
    1208             :          END IF
    1209             :       END IF
    1210       25531 :       IF ((.NOT. found) .AND. (ASSOCIATED(charmm_map))) THEN
    1211             :          ! Then we go through the CHARMM library
    1212       39242 :          DO i = 1, SIZE(charmm_map%kind)
    1213       39242 :             IF (element_symbol == charmm_map%kind(i)) THEN
    1214             :                found = .TRUE.
    1215             :                EXIT
    1216             :             END IF
    1217             :          END DO
    1218         436 :          IF (found) THEN
    1219         196 :             element_out = charmm_map%element(i)
    1220             :          END IF
    1221             :       END IF
    1222       25531 :       IF ((.NOT. found) .AND. (ASSOCIATED(gromos_map))) THEN
    1223             :          ! Then we go through the GROMOS library
    1224        5520 :          DO i = 1, SIZE(gromos_map%kind)
    1225        5520 :             IF (element_symbol == gromos_map%kind(i)) THEN
    1226             :                found = .TRUE.
    1227             :                EXIT
    1228             :             END IF
    1229             :          END DO
    1230         240 :          IF (found) THEN
    1231           0 :             element_out = gromos_map%element(i)
    1232             :          END IF
    1233             :       END IF
    1234             : 
    1235             :       ! final check. We has a connectivity, so we first tried to match against the ff_maps, but the element was not there.
    1236             :       ! Try again the periodic table.
    1237       25531 :       IF (.NOT. found) THEN
    1238         240 :          in_ptable = .FALSE.
    1239         240 :          IF (LEN_TRIM(element_symbol) <= 2) CALL get_ptable_info(element_symbol, found=in_ptable)
    1240         240 :          IF (in_ptable) THEN
    1241         240 :             element_out = TRIM(element_symbol)
    1242             :             found = .TRUE.
    1243             :          END IF
    1244             :       END IF
    1245             : 
    1246             :       ! If no element is found the job stops here.
    1247       25291 :       IF (.NOT. found) THEN
    1248             :          CALL cp_abort(__LOCATION__, &
    1249             :                        " Unknown element for KIND <"//TRIM(atom_name_in)//">."// &
    1250             :                        " This problem can be fixed specifying properly elements in PDB"// &
    1251             :                        " or specifying a KIND section or getting in touch with one of"// &
    1252           0 :                        " the developers!")
    1253             :       END IF
    1254             : 
    1255       25531 :    END SUBROUTINE check_subsys_element
    1256             : 
    1257           0 : END MODULE topology_util

Generated by: LCOV version 1.15