LCOV - code coverage report
Current view: top level - src - topology_coordinate_util.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.6 % 422 412
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 2 2

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 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_coordinate_util
      14              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      15              :                                               get_atomic_kind,&
      16              :                                               set_atomic_kind
      17              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      18              :                                               cp_logger_type
      19              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      20              :                                               cp_print_key_unit_nr
      21              :    USE exclusion_types,                 ONLY: exclusion_type
      22              :    USE external_potential_types,        ONLY: allocate_potential,&
      23              :                                               fist_potential_type,&
      24              :                                               get_potential,&
      25              :                                               set_potential
      26              :    USE input_constants,                 ONLY: do_fist,&
      27              :                                               do_skip_12,&
      28              :                                               do_skip_13,&
      29              :                                               do_skip_14
      30              :    USE input_section_types,             ONLY: section_vals_get,&
      31              :                                               section_vals_get_subs_vals,&
      32              :                                               section_vals_type,&
      33              :                                               section_vals_val_get
      34              :    USE kinds,                           ONLY: default_string_length,&
      35              :                                               dp
      36              :    USE memory_utilities,                ONLY: reallocate
      37              :    USE molecule_kind_types,             ONLY: atom_type,&
      38              :                                               get_molecule_kind,&
      39              :                                               molecule_kind_type,&
      40              :                                               set_molecule_kind
      41              :    USE molecule_types,                  ONLY: get_molecule,&
      42              :                                               molecule_type
      43              :    USE particle_types,                  ONLY: allocate_particle_set,&
      44              :                                               particle_type
      45              :    USE physcon,                         ONLY: massunit
      46              :    USE qmmm_types_low,                  ONLY: qmmm_env_mm_type
      47              :    USE string_table,                    ONLY: id2str,&
      48              :                                               s2s,&
      49              :                                               str2id
      50              :    USE topology_types,                  ONLY: atom_info_type,&
      51              :                                               connectivity_info_type,&
      52              :                                               topology_parameters_type
      53              :    USE topology_util,                   ONLY: array1_list_type,&
      54              :                                               reorder_structure
      55              : #include "./base/base_uses.f90"
      56              : 
      57              :    IMPLICIT NONE
      58              : 
      59              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_coordinate_util'
      60              : 
      61              :    PRIVATE
      62              :    PUBLIC :: topology_coordinate_pack
      63              : 
      64              : CONTAINS
      65              : 
      66              : ! **************************************************************************************************
      67              : !> \brief Take info readin from different file format and stuff it into
      68              : !>      compatible data structure in cp2k
      69              : !> \param particle_set ...
      70              : !> \param atomic_kind_set ...
      71              : !> \param molecule_kind_set ...
      72              : !> \param molecule_set ...
      73              : !> \param topology ...
      74              : !> \param qmmm ...
      75              : !> \param qmmm_env ...
      76              : !> \param subsys_section ...
      77              : !> \param force_env_section ...
      78              : !> \param exclusions ...
      79              : !> \param ignore_outside_box ...
      80              : !> \par History
      81              : !>      Teodoro Laino - modified in order to optimize the list of molecules
      82              : !>                      to build the exclusion lists
      83              : ! **************************************************************************************************
      84        10230 :    SUBROUTINE topology_coordinate_pack(particle_set, atomic_kind_set, &
      85              :                                        molecule_kind_set, molecule_set, topology, qmmm, qmmm_env, &
      86              :                                        subsys_section, force_env_section, exclusions, ignore_outside_box)
      87              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      88              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      89              :       TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
      90              :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
      91              :       TYPE(topology_parameters_type), INTENT(INOUT)      :: topology
      92              :       LOGICAL, INTENT(IN), OPTIONAL                      :: qmmm
      93              :       TYPE(qmmm_env_mm_type), OPTIONAL, POINTER          :: qmmm_env
      94              :       TYPE(section_vals_type), POINTER                   :: subsys_section, force_env_section
      95              :       TYPE(exclusion_type), DIMENSION(:), OPTIONAL, &
      96              :          POINTER                                         :: exclusions
      97              :       LOGICAL, INTENT(IN), OPTIONAL                      :: ignore_outside_box
      98              : 
      99              :       CHARACTER(len=*), PARAMETER :: routineN = 'topology_coordinate_pack'
     100              : 
     101              :       CHARACTER(LEN=default_string_length)               :: atmname
     102              :       INTEGER                                            :: atom_i, atom_j, counter, dim0, dim1, &
     103              :                                                             dim2, dim3, first, handle, handle2, i, &
     104              :                                                             iatom, ikind, iw, j, k, last, &
     105              :                                                             method_name_id, n, natom
     106        10230 :       INTEGER, DIMENSION(:), POINTER                     :: iatomlist, id_element, id_work, kind_of, &
     107        10230 :                                                             list, list2, molecule_list, &
     108        10230 :                                                             natom_of_kind, wlist
     109        10230 :       INTEGER, DIMENSION(:, :), POINTER                  :: pairs
     110              :       LOGICAL :: autogen, check, disable_exclusion_lists, do_center, explicit, found, &
     111              :          my_ignore_outside_box, my_qmmm, present_12_excl_ei_list, present_12_excl_vdw_list
     112              :       REAL(KIND=dp)                                      :: bounds(2, 3), cdims(3), dims(3), qeff, &
     113              :                                                             vec(3)
     114        10230 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: charge, cpoint, mass
     115        10230 :       TYPE(array1_list_type), DIMENSION(:), POINTER      :: ex_bend_list, ex_bond_list, &
     116        10230 :                                                             ex_bond_list_ei, ex_bond_list_vdw, &
     117        10230 :                                                             ex_onfo_list
     118              :       TYPE(atom_info_type), POINTER                      :: atom_info
     119        10230 :       TYPE(atom_type), DIMENSION(:), POINTER             :: atom_list
     120              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     121              :       TYPE(connectivity_info_type), POINTER              :: conn_info
     122              :       TYPE(cp_logger_type), POINTER                      :: logger
     123              :       TYPE(fist_potential_type), POINTER                 :: fist_potential
     124              :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     125              :       TYPE(molecule_type), POINTER                       :: molecule
     126              :       TYPE(section_vals_type), POINTER                   :: exclude_section, topology_section
     127              : 
     128        10230 :       NULLIFY (logger)
     129        20460 :       logger => cp_get_default_logger()
     130              :       iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
     131        10230 :                                 extension=".subsysLog")
     132        10230 :       topology_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY")
     133        10230 :       CALL timeset(routineN, handle)
     134              : 
     135        10230 :       my_qmmm = .FALSE.
     136        10230 :       IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
     137        10230 :       atom_info => topology%atom_info
     138        10230 :       conn_info => topology%conn_info
     139              :       !-----------------------------------------------------------------------------
     140              :       !-----------------------------------------------------------------------------
     141              :       ! 1. Determine topology%[natom_type,atom_names] and save mass(natom_type)
     142              :       !    and element(natom_type)
     143              :       !-----------------------------------------------------------------------------
     144        10230 :       CALL timeset(routineN//'_1', handle2)
     145        10230 :       counter = 0
     146        10230 :       NULLIFY (id_work, mass, id_element, charge)
     147        30690 :       ALLOCATE (id_work(topology%natoms))
     148        30690 :       ALLOCATE (mass(topology%natoms))
     149        20460 :       ALLOCATE (id_element(topology%natoms))
     150        20460 :       ALLOCATE (charge(topology%natoms))
     151       764757 :       id_work = str2id(s2s(""))
     152        10230 :       IF (iw > 0) WRITE (iw, *) "molecule_kind_set ::", SIZE(molecule_kind_set)
     153       146661 :       DO i = 1, SIZE(molecule_kind_set)
     154       136431 :          j = molecule_kind_set(i)%molecule_list(1)
     155       136431 :          molecule => molecule_set(j)
     156       136431 :          molecule_kind => molecule_set(j)%molecule_kind
     157       136431 :          IF (iw > 0) WRITE (iw, *) "molecule number ::", j, " has molecule kind number ::", i
     158              :          CALL get_molecule_kind(molecule_kind=molecule_kind, &
     159       136431 :                                 natom=natom, atom_list=atom_list)
     160              :          CALL get_molecule(molecule=molecule, &
     161       136431 :                            first_atom=first, last_atom=last)
     162       136431 :          IF (iw > 0) WRITE (iw, *) "boundaries of molecules (first, last) ::", first, last
     163       544358 :          DO j = 1, natom
     164     18758351 :             IF (.NOT. ANY(id_work(1:counter) .EQ. atom_list(j)%id_name)) THEN
     165        25649 :                counter = counter + 1
     166        25649 :                id_work(counter) = atom_list(j)%id_name
     167        25649 :                mass(counter) = atom_info%atm_mass(first + j - 1)
     168        25649 :                id_element(counter) = atom_info%id_element(first + j - 1)
     169        25649 :                charge(counter) = atom_info%atm_charge(first + j - 1)
     170        25649 :                IF (iw > 0) WRITE (iw, '(7X,A,1X,A5,F10.5,5X,A2,5X,F10.5)') &
     171           83 :                   "NEW ATOMIC KIND", id2str(id_work(counter)), mass(counter), id2str(id_element(counter)), charge(counter)
     172              :             ELSE
     173     18448285 :                found = .FALSE.
     174     18448285 :                DO k = 1, counter
     175     18448285 :                   IF ((id_work(k) == atom_list(j)%id_name) .AND. (charge(k) == atom_info%atm_charge(first + j - 1))) THEN
     176              :                      found = .TRUE.
     177              :                      EXIT
     178              :                   END IF
     179              :                END DO
     180       235617 :                IF (.NOT. found) THEN
     181          524 :                   counter = counter + 1
     182          524 :                   id_work(counter) = atom_list(j)%id_name
     183          524 :                   mass(counter) = atom_info%atm_mass(first + j - 1)
     184          524 :                   id_element(counter) = atom_info%id_element(first + j - 1)
     185          524 :                   charge(counter) = atom_info%atm_charge(first + j - 1)
     186          524 :                   IF (iw > 0) WRITE (iw, '(7X,A,1X,A5,F10.5,5X,A2,5X,F10.5)') &
     187            0 :                      "NEW ATOMIC KIND", id2str(id_work(counter)), mass(counter), id2str(id_element(counter)), charge(counter)
     188              :                END IF
     189              :             END IF
     190              :          END DO
     191              :       END DO
     192        10230 :       topology%natom_type = counter
     193        30690 :       ALLOCATE (atom_info%id_atom_names(topology%natom_type))
     194        36403 :       DO k = 1, counter
     195        36403 :          atom_info%id_atom_names(k) = id_work(k)
     196              :       END DO
     197        10230 :       DEALLOCATE (id_work)
     198        10230 :       CALL reallocate(mass, 1, counter)
     199        10230 :       CALL reallocate(id_element, 1, counter)
     200        10230 :       CALL reallocate(charge, 1, counter)
     201        10230 :       IF (iw > 0) &
     202           26 :          WRITE (iw, '(5X,A,I3)') "Total Number of Atomic Kinds = ", topology%natom_type
     203        10230 :       CALL timestop(handle2)
     204              : 
     205              :       !-----------------------------------------------------------------------------
     206              :       !-----------------------------------------------------------------------------
     207              :       ! 2. Allocate the data structure for the atomic kind information
     208              :       !-----------------------------------------------------------------------------
     209        10230 :       CALL timeset(routineN//'_2', handle2)
     210        10230 :       NULLIFY (atomic_kind_set)
     211        56863 :       ALLOCATE (atomic_kind_set(topology%natom_type))
     212        10230 :       CALL timestop(handle2)
     213              : 
     214              :       !-----------------------------------------------------------------------------
     215              :       !-----------------------------------------------------------------------------
     216              :       ! 3.  Allocate the data structure for the atomic information
     217              :       !-----------------------------------------------------------------------------
     218        10230 :       CALL timeset(routineN//'_3', handle2)
     219        10230 :       NULLIFY (particle_set)
     220        10230 :       CALL allocate_particle_set(particle_set, topology%natoms)
     221        10230 :       CALL timestop(handle2)
     222              : 
     223              :       !-----------------------------------------------------------------------------
     224              :       !-----------------------------------------------------------------------------
     225              :       ! 4. Set the atomic_kind_set(ikind)%[name,kind_number,mass]
     226              :       !-----------------------------------------------------------------------------
     227        10230 :       CALL timeset(routineN//'_4', handle2)
     228        36403 :       DO i = 1, topology%natom_type
     229        26173 :          atomic_kind => atomic_kind_set(i)
     230        26173 :          mass(i) = mass(i)*massunit
     231              :          CALL set_atomic_kind(atomic_kind=atomic_kind, &
     232              :                               kind_number=i, &
     233              :                               name=id2str(atom_info%id_atom_names(i)), &
     234              :                               element_symbol=id2str(id_element(i)), &
     235        26173 :                               mass=mass(i))
     236        36403 :          IF (iw > 0) THEN
     237           83 :             WRITE (iw, '(A,I5,A,I5,4A)') "Atomic Kind n.:", i, " out of:", topology%natom_type, &
     238           83 :                " name:   ", TRIM(id2str(atom_info%id_atom_names(i))), "   element:   ", &
     239          166 :                TRIM(id2str(id_element(i)))
     240              :          END IF
     241              :       END DO
     242        10230 :       DEALLOCATE (mass)
     243        10230 :       DEALLOCATE (id_element)
     244        10230 :       CALL timestop(handle2)
     245              : 
     246              :       !-----------------------------------------------------------------------------
     247              :       !-----------------------------------------------------------------------------
     248              :       ! 5. Determine number of atom of each kind (ie natom_of_kind and kind_of)
     249              :       !-----------------------------------------------------------------------------
     250        10230 :       CALL timeset(routineN//'_5', handle2)
     251        30690 :       ALLOCATE (kind_of(topology%natoms))
     252        30690 :       ALLOCATE (natom_of_kind(topology%natom_type))
     253       764757 :       kind_of(:) = 0
     254        36403 :       natom_of_kind(:) = 0
     255        36403 :       DO i = 1, topology%natom_type
     256     37089735 :          DO j = 1, topology%natoms
     257     37079505 :             IF ((atom_info%id_atom_names(i) == atom_info%id_atmname(j)) .AND. (charge(i) == atom_info%atm_charge(j))) THEN
     258       754527 :                natom_of_kind(i) = natom_of_kind(i) + 1
     259       754527 :                IF (kind_of(j) == 0) kind_of(j) = i
     260              :             END IF
     261              :          END DO
     262              :       END DO
     263       764757 :       IF (ANY(kind_of == 0)) THEN
     264            0 :          DO i = 1, topology%natoms
     265            0 :             IF (kind_of(i) == 0) THEN
     266            0 :                WRITE (*, *) i, kind_of(i)
     267            0 :                WRITE (*, *) "Two molecules have been defined as identical molecules but atoms mismatch charges!"
     268              :             END IF
     269              :          END DO
     270            0 :          CPABORT("")
     271              :       END IF
     272        10230 :       CALL timestop(handle2)
     273              : 
     274              :       !-----------------------------------------------------------------------------
     275              :       !-----------------------------------------------------------------------------
     276              :       ! 6. Set the atom_kind_set(ikind)%[natom,atom_list]
     277              :       !-----------------------------------------------------------------------------
     278        10230 :       CALL timeset(routineN//'_6', handle2)
     279        36403 :       DO i = 1, topology%natom_type
     280        26173 :          atomic_kind => atomic_kind_set(i)
     281              :          NULLIFY (iatomlist)
     282        78519 :          ALLOCATE (iatomlist(natom_of_kind(i)))
     283        26173 :          counter = 0
     284     37079505 :          DO j = 1, topology%natoms
     285     37079505 :             IF (kind_of(j) == i) THEN
     286       754527 :                counter = counter + 1
     287       754527 :                iatomlist(counter) = j
     288              :             END IF
     289              :          END DO
     290        26173 :          IF (iw > 0) THEN
     291           83 :             WRITE (iw, '(A,I6,A)') "      Atomic kind ", i, " contains particles"
     292         1614 :             DO J = 1, SIZE(iatomlist)
     293         1614 :                IF (MOD(J, 5) .EQ. 0) THEN ! split long lines
     294          271 :                   WRITE (iw, '(I12)') iatomlist(J)
     295              :                ELSE
     296         1260 :                   WRITE (iw, '(I12)', ADVANCE="NO") iatomlist(J)
     297              :                END IF
     298              :             END DO
     299           83 :             WRITE (iw, *)
     300              :          END IF
     301              :          CALL set_atomic_kind(atomic_kind=atomic_kind, &
     302              :                               natom=natom_of_kind(i), &
     303        26173 :                               atom_list=iatomlist)
     304        36403 :          DEALLOCATE (iatomlist)
     305              :       END DO
     306        10230 :       DEALLOCATE (natom_of_kind)
     307        10230 :       CALL timestop(handle2)
     308              : 
     309              :       !-----------------------------------------------------------------------------
     310              :       !-----------------------------------------------------------------------------
     311              :       ! 7. Possibly center the coordinates and fill in coordinates in particle_set
     312              :       !-----------------------------------------------------------------------------
     313              :       CALL section_vals_val_get(subsys_section, &
     314        10230 :                                 "TOPOLOGY%CENTER_COORDINATES%_SECTION_PARAMETERS_", l_val=do_center)
     315        10230 :       CALL timeset(routineN//'_7a', handle2)
     316       774987 :       bounds(1, 1) = MINVAL(atom_info%r(1, :))
     317       774987 :       bounds(2, 1) = MAXVAL(atom_info%r(1, :))
     318              : 
     319       774987 :       bounds(1, 2) = MINVAL(atom_info%r(2, :))
     320       774987 :       bounds(2, 2) = MAXVAL(atom_info%r(2, :))
     321              : 
     322       774987 :       bounds(1, 3) = MINVAL(atom_info%r(3, :))
     323       774987 :       bounds(2, 3) = MAXVAL(atom_info%r(3, :))
     324              : 
     325        40920 :       dims = bounds(2, :) - bounds(1, :)
     326        10230 :       cdims(1) = topology%cell%hmat(1, 1)
     327        10230 :       cdims(2) = topology%cell%hmat(2, 2)
     328        10230 :       cdims(3) = topology%cell%hmat(3, 3)
     329        10230 :       IF (iw > 0) THEN
     330           26 :          WRITE (iw, '(A,3F12.6)') "System sizes: ", dims, "Cell sizes (diagonal): ", cdims
     331              :       END IF
     332        10230 :       check = .TRUE.
     333        40920 :       DO i = 1, 3
     334        40920 :          IF (topology%cell%perd(i) == 0) THEN
     335        10880 :             check = check .AND. (dims(i) < cdims(i))
     336              :          END IF
     337              :       END DO
     338        10230 :       my_ignore_outside_box = .FALSE.
     339        10230 :       IF (PRESENT(ignore_outside_box)) my_ignore_outside_box = ignore_outside_box
     340        10230 :       IF (.NOT. my_ignore_outside_box .AND. .NOT. check) &
     341              :          CALL cp_abort(__LOCATION__, &
     342              :                        "A non-periodic calculation has been requested but the system size "// &
     343            0 :                        "exceeds the cell size in at least one of the non-periodic directions!")
     344        10230 :       IF (do_center) THEN
     345              :          CALL section_vals_val_get(subsys_section, &
     346         2206 :                                    "TOPOLOGY%CENTER_COORDINATES%CENTER_POINT", explicit=explicit)
     347         2206 :          IF (explicit) THEN
     348              :             CALL section_vals_val_get(subsys_section, &
     349            0 :                                       "TOPOLOGY%CENTER_COORDINATES%CENTER_POINT", r_vals=cpoint)
     350            0 :             vec = cpoint
     351              :          ELSE
     352         8824 :             vec = cdims/2.0_dp
     353              :          END IF
     354         8824 :          dims = (bounds(2, :) + bounds(1, :))/2.0_dp - vec
     355              :       ELSE
     356         8024 :          dims = 0.0_dp
     357              :       END IF
     358        10230 :       CALL timestop(handle2)
     359        10230 :       CALL timeset(routineN//'_7b', handle2)
     360       764757 :       DO i = 1, topology%natoms
     361       754527 :          ikind = kind_of(i)
     362       754527 :          IF (iw > 0) THEN
     363         1531 :             WRITE (iw, *) "atom number :: ", i, "kind number ::", ikind
     364              :          END IF
     365       754527 :          particle_set(i)%atomic_kind => atomic_kind_set(ikind)
     366      5281689 :          particle_set(i)%r(:) = atom_info%r(:, i) - dims
     367       764757 :          particle_set(i)%atom_index = i
     368              :       END DO
     369        10230 :       CALL timestop(handle2)
     370        10230 :       DEALLOCATE (kind_of)
     371              : 
     372              :       !-----------------------------------------------------------------------------
     373              :       !-----------------------------------------------------------------------------
     374              :       ! 8. Fill in the exclusions%list_exclude_vdw
     375              :       ! 9. Fill in the exclusions%list_exclude_ei
     376              :       ! 10. Fill in the exclusions%list_onfo
     377              :       !-----------------------------------------------------------------------------
     378        10230 :       CALL timeset(routineN//'_89', handle2)
     379        10230 :       CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id)
     380              :       CALL section_vals_val_get(subsys_section, "TOPOLOGY%DISABLE_EXCLUSION_LISTS", &
     381        10230 :                                 l_val=disable_exclusion_lists)
     382        10230 :       IF ((method_name_id == do_fist) .AND. (.NOT. disable_exclusion_lists)) THEN
     383         2498 :          CPASSERT(PRESENT(exclusions))
     384         2498 :          natom = topology%natoms
     385              :          ! allocate exclusions. Most likely they would only be needed for the local_particles
     386       640830 :          ALLOCATE (exclusions(natom))
     387       635834 :          DO I = 1, natom
     388       633336 :             NULLIFY (exclusions(i)%list_exclude_vdw)
     389       633336 :             NULLIFY (exclusions(i)%list_exclude_ei)
     390       635834 :             NULLIFY (exclusions(i)%list_onfo)
     391              :          END DO
     392              :          ! Reorder bonds
     393       640830 :          ALLOCATE (ex_bond_list(natom))
     394       635834 :          DO I = 1, natom
     395       635834 :             ALLOCATE (ex_bond_list(I)%array1(0))
     396              :          END DO
     397         2498 :          N = 0
     398         2498 :          IF (ASSOCIATED(conn_info%bond_a)) THEN
     399         2498 :             N = SIZE(conn_info%bond_a)
     400         2498 :             CALL reorder_structure(ex_bond_list, conn_info%bond_a, conn_info%bond_b, N)
     401              :          END IF
     402              : 
     403              :          ! Check if a list of 1-2 exclusion bonds is defined.. if not use all bonds
     404         2498 :          NULLIFY (ex_bond_list_vdw, ex_bond_list_ei)
     405              :          ! VdW
     406         2498 :          exclude_section => section_vals_get_subs_vals(topology_section, "EXCLUDE_VDW_LIST")
     407         2498 :          CALL section_vals_get(exclude_section, explicit=explicit)
     408         2498 :          present_12_excl_vdw_list = .FALSE.
     409         2498 :          IF (explicit) present_12_excl_vdw_list = .TRUE.
     410              :          IF (present_12_excl_vdw_list) THEN
     411           40 :             ALLOCATE (ex_bond_list_vdw(natom))
     412           32 :             DO I = 1, natom
     413           32 :                ALLOCATE (ex_bond_list_vdw(I)%array1(0))
     414              :             END DO
     415              :             CALL setup_exclusion_list(exclude_section, "BOND", ex_bond_list, ex_bond_list_vdw, &
     416            8 :                                       particle_set)
     417              :          ELSE
     418         2490 :             ex_bond_list_vdw => ex_bond_list
     419              :          END IF
     420              :          ! EI
     421         2498 :          exclude_section => section_vals_get_subs_vals(topology_section, "EXCLUDE_EI_LIST")
     422         2498 :          CALL section_vals_get(exclude_section, explicit=explicit)
     423         2498 :          present_12_excl_ei_list = .FALSE.
     424         2498 :          IF (explicit) present_12_excl_ei_list = .TRUE.
     425              :          IF (present_12_excl_ei_list) THEN
     426           50 :             ALLOCATE (ex_bond_list_ei(natom))
     427           40 :             DO I = 1, natom
     428           40 :                ALLOCATE (ex_bond_list_ei(I)%array1(0))
     429              :             END DO
     430              :             CALL setup_exclusion_list(exclude_section, "BOND", ex_bond_list, ex_bond_list_ei, &
     431           10 :                                       particle_set)
     432              :          ELSE
     433         2488 :             ex_bond_list_ei => ex_bond_list
     434              :          END IF
     435              : 
     436              :          CALL section_vals_val_get(topology_section, "AUTOGEN_EXCLUDE_LISTS", &
     437         2498 :                                    l_val=autogen)
     438              :          ! Reorder bends
     439       638332 :          ALLOCATE (ex_bend_list(natom))
     440       635834 :          DO I = 1, natom
     441       635834 :             ALLOCATE (ex_bend_list(I)%array1(0))
     442              :          END DO
     443         2498 :          IF (autogen) THEN
     444              :             ! Construct autogenerated 1-3 pairs, i.e. all possible 1-3 pairs instead
     445              :             ! of only the bends that are present in the topology.
     446            4 :             ALLOCATE (pairs(0, 2))
     447            4 :             N = 0
     448           26 :             DO iatom = 1, natom
     449           62 :                DO i = 1, SIZE(ex_bond_list(iatom)%array1)
     450              :                   ! a neighboring atom of iatom:
     451           36 :                   atom_i = ex_bond_list(iatom)%array1(i)
     452           92 :                   DO j = 1, i - 1
     453              :                      ! another neighboring atom of iatom
     454           34 :                      atom_j = ex_bond_list(iatom)%array1(j)
     455              :                      ! It is only a true bend if there is no shorter path.
     456              :                      ! No need to check if i and j correspond to the same atom.
     457              :                      ! Check if i and j are not involved in a bond:
     458           34 :                      check = .FALSE.
     459           70 :                      DO counter = 1, SIZE(ex_bond_list(atom_i)%array1)
     460           70 :                         IF (ex_bond_list(atom_i)%array1(counter) == atom_j) THEN
     461              :                            check = .TRUE.
     462              :                            EXIT
     463              :                         END IF
     464              :                      END DO
     465           34 :                      IF (check) CYCLE
     466              :                      ! Add the genuine 1-3 pair
     467           34 :                      N = N + 1
     468           34 :                      IF (SIZE(pairs, dim=1) <= N) THEN
     469            8 :                         CALL reallocate(pairs, 1, N + 5, 1, 2)
     470              :                      END IF
     471           34 :                      pairs(N, 1) = atom_i
     472           70 :                      pairs(N, 2) = atom_j
     473              :                   END DO
     474              :                END DO
     475              :             END DO
     476            4 :             CALL reorder_structure(ex_bend_list, pairs(:, 1), pairs(:, 2), N)
     477            4 :             DEALLOCATE (pairs)
     478              :          ELSE
     479         2494 :             IF (ASSOCIATED(conn_info%theta_a)) THEN
     480         2494 :                N = SIZE(conn_info%theta_a)
     481         2494 :                CALL reorder_structure(ex_bend_list, conn_info%theta_a, conn_info%theta_c, N)
     482              :             END IF
     483              :          END IF
     484              : 
     485              :          ! Reorder onfo
     486       638332 :          ALLOCATE (ex_onfo_list(natom))
     487       635834 :          DO I = 1, natom
     488       635834 :             ALLOCATE (ex_onfo_list(I)%array1(0))
     489              :          END DO
     490         2498 :          IF (autogen) THEN
     491              :             ! Construct autogenerated 1-4 pairs, i.e. all possible 1-4 pairs instead
     492              :             ! of only the onfo's that are present in the topology.
     493            4 :             ALLOCATE (pairs(0, 2))
     494            4 :             N = 0
     495           26 :             DO iatom = 1, natom
     496           62 :                DO i = 1, SIZE(ex_bond_list(iatom)%array1)
     497              :                   ! a neighboring atom of iatom:
     498           36 :                   atom_i = ex_bond_list(iatom)%array1(i)
     499          130 :                   DO j = 1, SIZE(ex_bend_list(iatom)%array1)
     500              :                      ! a next neighboring atom of iatom:
     501           72 :                      atom_j = ex_bend_list(iatom)%array1(j)
     502              :                      ! It is only a true onfo if there is no shorter path.
     503              :                      ! check if i and j are not the same atom
     504           72 :                      IF (atom_i == atom_j) CYCLE
     505              :                      ! check if i and j are not involved in a bond
     506           72 :                      check = .FALSE.
     507          230 :                      DO counter = 1, SIZE(ex_bond_list(atom_i)%array1)
     508          230 :                         IF (ex_bond_list(atom_i)%array1(counter) == atom_j) THEN
     509              :                            check = .TRUE.
     510              :                            EXIT
     511              :                         END IF
     512              :                      END DO
     513           72 :                      IF (check) CYCLE
     514              :                      ! check if i and j are not involved in a bend
     515            4 :                      check = .FALSE.
     516            8 :                      DO counter = 1, SIZE(ex_bend_list(atom_i)%array1)
     517            8 :                         IF (ex_bend_list(atom_i)%array1(counter) == atom_j) THEN
     518              :                            check = .TRUE.
     519              :                            EXIT
     520              :                         END IF
     521              :                      END DO
     522            4 :                      IF (check) CYCLE
     523              :                      ! Add the true onfo.
     524            4 :                      N = N + 1
     525            4 :                      IF (SIZE(pairs, dim=1) <= N) THEN
     526            2 :                         CALL reallocate(pairs, 1, N + 5, 1, 2)
     527              :                      END IF
     528            4 :                      pairs(N, 1) = atom_i
     529          108 :                      pairs(N, 2) = atom_j
     530              :                   END DO
     531              :                END DO
     532              :             END DO
     533            4 :             CALL reorder_structure(ex_onfo_list, pairs(:, 1), pairs(:, 2), N)
     534            4 :             DEALLOCATE (pairs)
     535              :          ELSE
     536         2494 :             IF (ASSOCIATED(conn_info%onfo_a)) THEN
     537         2488 :                N = SIZE(conn_info%onfo_a)
     538         2488 :                CALL reorder_structure(ex_onfo_list, conn_info%onfo_a, conn_info%onfo_b, N)
     539              :             END IF
     540              :          END IF
     541              : 
     542              :          ! Build the exclusion (and onfo) list per atom.
     543       635834 :          DO iatom = 1, SIZE(particle_set)
     544              :             ! Setup exclusion list for VDW: always exclude itself
     545       633336 :             dim0 = 1
     546              :             ! exclude bond-neighbors (only if do_skip_12 .OR. do_skip_13 .OR. do_skip_14)
     547       633336 :             dim1 = 0
     548              :             IF (topology%exclude_vdw == do_skip_12 .OR. &
     549       633336 :                 topology%exclude_vdw == do_skip_13 .OR. &
     550       632742 :                 topology%exclude_vdw == do_skip_14) dim1 = SIZE(ex_bond_list_vdw(iatom)%array1)
     551       633336 :             dim1 = dim1 + dim0
     552       633336 :             dim2 = 0
     553       633336 :             IF (topology%exclude_vdw == do_skip_13 .OR. &
     554       632136 :                 topology%exclude_vdw == do_skip_14) dim2 = SIZE(ex_bend_list(iatom)%array1)
     555       633336 :             dim2 = dim1 + dim2
     556       633336 :             dim3 = 0
     557       633336 :             IF (topology%exclude_vdw == do_skip_14) dim3 = SIZE(ex_onfo_list(iatom)%array1)
     558       633336 :             dim3 = dim2 + dim3
     559       633336 :             IF (dim3 /= 0) THEN
     560       633336 :                NULLIFY (list, wlist)
     561      1900008 :                ALLOCATE (wlist(dim3))
     562      1266672 :                wlist(dim0:dim0) = iatom
     563      1603470 :                IF (dim1 > dim0) wlist(dim0 + 1:dim1) = ex_bond_list_vdw(iatom)%array1
     564      1153846 :                IF (dim2 > dim1) wlist(dim1 + 1:dim2) = ex_bend_list(iatom)%array1
     565       634884 :                IF (dim3 > dim2) wlist(dim2 + 1:dim3) = ex_onfo_list(iatom)%array1
     566              :                ! Get a unique list
     567      2125528 :                DO i = 1, SIZE(wlist) - 1
     568      1492192 :                   IF (wlist(i) == 0) CYCLE
     569      5406430 :                   DO j = i + 1, SIZE(wlist)
     570      4773164 :                      IF (wlist(j) == wlist(i)) wlist(j) = 0
     571              :                   END DO
     572              :                END DO
     573      2758864 :                dim3 = SIZE(wlist) - COUNT(wlist == 0)
     574      1900008 :                ALLOCATE (list(dim3))
     575       633336 :                j = 0
     576      2758864 :                DO i = 1, SIZE(wlist)
     577      2125528 :                   IF (wlist(i) == 0) CYCLE
     578      2039644 :                   j = j + 1
     579      2758864 :                   list(j) = wlist(i)
     580              :                END DO
     581       633336 :                DEALLOCATE (wlist)
     582              :                ! Unique list completed
     583       633336 :                NULLIFY (list2)
     584              :                IF ((topology%exclude_vdw == topology%exclude_ei) .AND. &
     585       633336 :                    (.NOT. present_12_excl_ei_list) .AND. (.NOT. present_12_excl_vdw_list)) THEN
     586              :                   list2 => list
     587              :                ELSE
     588              :                   ! Setup exclusion list for EI : always exclude itself
     589         1770 :                   dim0 = 1
     590              :                   ! exclude bond-neighbors (only if do_skip_12 .OR. do_skip_13 .OR. do_skip_14)
     591         1770 :                   dim1 = 0
     592              :                   IF (topology%exclude_ei == do_skip_12 .OR. &
     593         1770 :                       topology%exclude_ei == do_skip_13 .OR. &
     594         1326 :                       topology%exclude_ei == do_skip_14) dim1 = SIZE(ex_bond_list_ei(iatom)%array1)
     595         1770 :                   dim1 = dim1 + dim0
     596         1770 :                   dim2 = 0
     597         1770 :                   IF (topology%exclude_ei == do_skip_13 .OR. &
     598          864 :                       topology%exclude_ei == do_skip_14) dim2 = SIZE(ex_bend_list(iatom)%array1)
     599         1770 :                   dim2 = dim1 + dim2
     600         1770 :                   dim3 = 0
     601         1770 :                   IF (topology%exclude_ei == do_skip_14) dim3 = SIZE(ex_onfo_list(iatom)%array1)
     602         1770 :                   dim3 = dim2 + dim3
     603              : 
     604         1770 :                   IF (dim3 /= 0) THEN
     605         5310 :                      ALLOCATE (wlist(dim3))
     606         3540 :                      wlist(dim0:dim0) = iatom
     607         4098 :                      IF (dim1 > dim0) wlist(dim0 + 1:dim1) = ex_bond_list_ei(iatom)%array1
     608         4266 :                      IF (dim2 > dim1) wlist(dim1 + 1:dim2) = ex_bend_list(iatom)%array1
     609         2922 :                      IF (dim3 > dim2) wlist(dim2 + 1:dim3) = ex_onfo_list(iatom)%array1
     610              :                      ! Get a unique list
     611         7746 :                      DO i = 1, SIZE(wlist) - 1
     612         5976 :                         IF (wlist(i) == 0) CYCLE
     613        28896 :                         DO j = i + 1, SIZE(wlist)
     614        27126 :                            IF (wlist(j) == wlist(i)) wlist(j) = 0
     615              :                         END DO
     616              :                      END DO
     617         9516 :                      dim3 = SIZE(wlist) - COUNT(wlist == 0)
     618         5310 :                      ALLOCATE (list2(dim3))
     619         1770 :                      j = 0
     620         9516 :                      DO i = 1, SIZE(wlist)
     621         7746 :                         IF (wlist(i) == 0) CYCLE
     622         7746 :                         j = j + 1
     623         9516 :                         list2(j) = wlist(i)
     624              :                      END DO
     625         1770 :                      DEALLOCATE (wlist)
     626              :                      ! Unique list completed
     627              :                   END IF
     628              :                END IF
     629              :             END IF
     630       633336 :             exclusions(iatom)%list_exclude_vdw => list
     631       633336 :             exclusions(iatom)%list_exclude_ei => list2
     632              :             ! Keep a list of onfo atoms for proper selection of specialized 1-4
     633              :             ! potentials instead of conventional nonbonding potentials.
     634      1341211 :             ALLOCATE (exclusions(iatom)%list_onfo(SIZE(ex_onfo_list(iatom)%array1)))
     635              :             ! copy of data, not copy of pointer
     636      1245156 :             exclusions(iatom)%list_onfo = ex_onfo_list(iatom)%array1
     637       635834 :             IF (iw > 0) THEN
     638         1243 :                IF (ASSOCIATED(list)) &
     639         1243 :                   WRITE (iw, *) "exclusion list_vdw :: ", &
     640         1243 :                   "atom num :", iatom, "exclusion list ::", &
     641         6003 :                   list
     642         1243 :                IF (topology%exclude_vdw /= topology%exclude_ei) THEN
     643            9 :                   IF (ASSOCIATED(list2)) &
     644            9 :                      WRITE (iw, *) "exclusion list_ei :: ", &
     645            9 :                      "atom num :", iatom, "exclusion list ::", &
     646           35 :                      list2
     647              :                END IF
     648         1243 :                IF (ASSOCIATED(exclusions(iatom)%list_onfo)) &
     649         1243 :                   WRITE (iw, *) "onfo list :: ", &
     650         1243 :                   "atom num :", iatom, "onfo list ::", &
     651         3320 :                   exclusions(iatom)%list_onfo
     652              :             END IF
     653              :          END DO
     654              :          ! deallocate onfo
     655       635834 :          DO I = 1, natom
     656       635834 :             DEALLOCATE (ex_onfo_list(I)%array1)
     657              :          END DO
     658         2498 :          DEALLOCATE (ex_onfo_list)
     659              :          ! deallocate bends
     660       635834 :          DO I = 1, natom
     661       635834 :             DEALLOCATE (ex_bend_list(I)%array1)
     662              :          END DO
     663         2498 :          DEALLOCATE (ex_bend_list)
     664              :          ! deallocate bonds
     665         2498 :          IF (present_12_excl_ei_list) THEN
     666           40 :             DO I = 1, natom
     667           40 :                DEALLOCATE (ex_bond_list_ei(I)%array1)
     668              :             END DO
     669           10 :             DEALLOCATE (ex_bond_list_ei)
     670              :          ELSE
     671         2488 :             NULLIFY (ex_bond_list_ei)
     672              :          END IF
     673         2498 :          IF (present_12_excl_vdw_list) THEN
     674           32 :             DO I = 1, natom
     675           32 :                DEALLOCATE (ex_bond_list_vdw(I)%array1)
     676              :             END DO
     677            8 :             DEALLOCATE (ex_bond_list_vdw)
     678              :          ELSE
     679         2490 :             NULLIFY (ex_bond_list_vdw)
     680              :          END IF
     681       635834 :          DO I = 1, natom
     682       635834 :             DEALLOCATE (ex_bond_list(I)%array1)
     683              :          END DO
     684         9992 :          DEALLOCATE (ex_bond_list)
     685              :       END IF
     686        10230 :       CALL timestop(handle2)
     687              :       !-----------------------------------------------------------------------------
     688              :       !-----------------------------------------------------------------------------
     689              :       ! 11. Set the atomic_kind_set()%fist_potential%[qeff] (PART 1)
     690              :       !-----------------------------------------------------------------------------
     691        10230 :       CALL timeset(routineN//'_10', handle2)
     692        10230 :       CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id)
     693        10230 :       IF (method_name_id == do_fist) THEN
     694        13928 :          DO i = 1, SIZE(atomic_kind_set)
     695        11282 :             atomic_kind => atomic_kind_set(i)
     696        11282 :             CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
     697        11282 :             qeff = charge(i)
     698        11282 :             NULLIFY (fist_potential)
     699        11282 :             CALL allocate_potential(fist_potential)
     700        11282 :             CALL set_potential(potential=fist_potential, qeff=qeff)
     701        13928 :             CALL set_atomic_kind(atomic_kind=atomic_kind, fist_potential=fist_potential)
     702              :          END DO
     703              :       END IF
     704        10230 :       DEALLOCATE (charge)
     705        10230 :       CALL timestop(handle2)
     706              : 
     707              :       !-----------------------------------------------------------------------------
     708              :       !-----------------------------------------------------------------------------
     709              :       ! 12. Set the atom_list for molecule_kind in molecule_kind_set (PART 2)
     710              :       !-----------------------------------------------------------------------------
     711        10230 :       CALL timeset(routineN//'_11', handle2)
     712       146661 :       DO i = 1, SIZE(molecule_kind_set)
     713       136431 :          molecule_kind => molecule_kind_set(i)
     714              :          CALL get_molecule_kind(molecule_kind=molecule_kind, &
     715              :                                 natom=natom, molecule_list=molecule_list, &
     716       136431 :                                 atom_list=atom_list)
     717       136431 :          molecule => molecule_set(molecule_list(1))
     718              :          CALL get_molecule(molecule=molecule, &
     719       136431 :                            first_atom=first, last_atom=last)
     720       397697 :          DO j = 1, natom
     721     18850527 :             DO k = 1, SIZE(atomic_kind_set)
     722     18714096 :                atomic_kind => atomic_kind_set(k)
     723     18714096 :                CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
     724     18714096 :                IF (method_name_id == do_fist) THEN
     725     18387061 :                   CALL get_atomic_kind(atomic_kind=atomic_kind, fist_potential=fist_potential)
     726     18387061 :                   CALL get_potential(potential=fist_potential, qeff=qeff)
     727     18387061 :                   IF ((id2str(atom_list(j)%id_name) == atmname) .AND. (qeff == atom_info%atm_charge(first + j - 1))) THEN
     728       195761 :                      atom_list(j)%atomic_kind => atomic_kind_set(k)
     729     18582822 :                      EXIT
     730              :                   END IF
     731              :                ELSE
     732       327035 :                   IF (id2str(atom_list(j)%id_name) == atmname) THEN
     733        65505 :                      atom_list(j)%atomic_kind => atomic_kind_set(k)
     734       392540 :                      EXIT
     735              :                   END IF
     736              :                END IF
     737              :             END DO
     738              :          END DO
     739       283092 :          CALL set_molecule_kind(molecule_kind=molecule_kind, atom_list=atom_list)
     740              :       END DO
     741        10230 :       CALL timestop(handle2)
     742              : 
     743        10230 :       CALL timestop(handle)
     744              :       CALL cp_print_key_finished_output(iw, logger, subsys_section, &
     745        10230 :                                         "PRINT%TOPOLOGY_INFO/UTIL_INFO")
     746       122760 :    END SUBROUTINE topology_coordinate_pack
     747              : 
     748              : ! **************************************************************************************************
     749              : !> \brief Builds the exclusion list for VDW and EI if an explicit list of terms
     750              : !>        is provided by the user. Otherwise all possibilities are excluded
     751              : !> \param exclude_section ...
     752              : !> \param keyword ...
     753              : !> \param ex_bond_list ...
     754              : !> \param ex_bond_list_w ...
     755              : !> \param particle_set ...
     756              : !> \par History
     757              : !>      Teodoro Laino [tlaino] - 12.2009
     758              : ! **************************************************************************************************
     759           18 :    SUBROUTINE setup_exclusion_list(exclude_section, keyword, ex_bond_list, &
     760              :                                    ex_bond_list_w, particle_set)
     761              :       TYPE(section_vals_type), POINTER                   :: exclude_section
     762              :       CHARACTER(LEN=*), INTENT(IN)                       :: keyword
     763              :       TYPE(array1_list_type), DIMENSION(:), POINTER      :: ex_bond_list, ex_bond_list_w
     764              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     765              : 
     766              :       CHARACTER(LEN=default_string_length)               :: flag1, flag2
     767              :       CHARACTER(LEN=default_string_length), &
     768           18 :          DIMENSION(:), POINTER                           :: names
     769              :       INTEGER                                            :: i, ind, j, k, l, m, n_rep
     770              : 
     771            0 :       CPASSERT(ASSOCIATED(ex_bond_list))
     772           18 :       CPASSERT(ASSOCIATED(ex_bond_list_w))
     773           18 :       SELECT CASE (keyword)
     774              :       CASE ("BOND")
     775           18 :          CALL section_vals_val_get(exclude_section, keyword, n_rep_val=n_rep)
     776           72 :          DO j = 1, SIZE(ex_bond_list)
     777           54 :             CPASSERT(ASSOCIATED(ex_bond_list(j)%array1))
     778           54 :             CPASSERT(ASSOCIATED(ex_bond_list_w(j)%array1))
     779              : 
     780           54 :             flag1 = particle_set(j)%atomic_kind%name
     781           54 :             m = SIZE(ex_bond_list(j)%array1)
     782           54 :             CALL reallocate(ex_bond_list_w(j)%array1, 1, m)
     783              : 
     784           54 :             l = 0
     785          126 :             DO k = 1, m
     786           72 :                ind = ex_bond_list(j)%array1(k)
     787           72 :                flag2 = particle_set(ind)%atomic_kind%name
     788          150 :                DO i = 1, n_rep
     789              :                   CALL section_vals_val_get(exclude_section, keyword, i_rep_val=i, &
     790           24 :                                             c_vals=names)
     791           24 :                   IF (((TRIM(names(1)) == TRIM(flag1)) .AND. (TRIM(names(2)) == TRIM(flag2))) .OR. &
     792           72 :                       ((TRIM(names(1)) == TRIM(flag2)) .AND. (TRIM(names(2)) == TRIM(flag1)))) THEN
     793           24 :                      l = l + 1
     794           24 :                      ex_bond_list_w(j)%array1(l) = ind
     795              :                   END IF
     796              :                END DO
     797              :             END DO
     798           72 :             CALL reallocate(ex_bond_list_w(j)%array1, 1, l)
     799              :          END DO
     800              :       CASE DEFAULT
     801           18 :          CPABORT("")
     802              :       END SELECT
     803              : 
     804           18 :    END SUBROUTINE setup_exclusion_list
     805              : 
     806              : END MODULE topology_coordinate_util
        

Generated by: LCOV version 2.0-1