LCOV - code coverage report
Current view: top level - src - fist_neighbor_lists.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 99.0 % 392 388
Test Date: 2025-07-25 12:55:17 Functions: 80.0 % 5 4

            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 Generate the atomic neighbor lists for FIST.
      10              : !> \par History
      11              : !>      - build and update merged (11.02.2005,MK)
      12              : !>      - bug fix for PERIODIC NONE (24.02.06,MK)
      13              : !>      - Major rewriting (light memory neighbor lists): teo and joost (05.2006)
      14              : !>      - Completely new algorithm for the neighbor lists
      15              : !>        (faster and memory lighter) (Teo 08.2006)
      16              : !> \author MK (19.11.2002,24.07.2003)
      17              : !>      Teodoro Laino (08.2006) - MAJOR REWRITING
      18              : ! **************************************************************************************************
      19              : MODULE fist_neighbor_lists
      20              : 
      21              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      22              :                                               get_atomic_kind,&
      23              :                                               get_atomic_kind_set
      24              :    USE cell_types,                      ONLY: cell_type,&
      25              :                                               get_cell,&
      26              :                                               pbc,&
      27              :                                               plane_distance,&
      28              :                                               scaled_to_real
      29              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      30              :                                               cp_logger_get_default_unit_nr,&
      31              :                                               cp_logger_type
      32              :    USE cp_output_handling,              ONLY: cp_p_file,&
      33              :                                               cp_print_key_finished_output,&
      34              :                                               cp_print_key_should_output,&
      35              :                                               cp_print_key_unit_nr
      36              :    USE cp_units,                        ONLY: cp_unit_from_cp2k
      37              :    USE distribution_1d_types,           ONLY: distribution_1d_type
      38              :    USE exclusion_types,                 ONLY: exclusion_type
      39              :    USE fist_neighbor_list_types,        ONLY: fist_neighbor_add,&
      40              :                                               fist_neighbor_deallocate,&
      41              :                                               fist_neighbor_init,&
      42              :                                               fist_neighbor_type,&
      43              :                                               neighbor_kind_pairs_type
      44              :    USE input_section_types,             ONLY: section_vals_type,&
      45              :                                               section_vals_val_get
      46              :    USE kinds,                           ONLY: default_string_length,&
      47              :                                               dp
      48              :    USE memory_utilities,                ONLY: reallocate
      49              :    USE message_passing,                 ONLY: mp_para_env_type
      50              :    USE particle_types,                  ONLY: particle_type
      51              :    USE qmmm_ff_fist,                    ONLY: qmmm_ff_precond_only_qm
      52              :    USE string_utilities,                ONLY: compress
      53              :    USE subcell_types,                   ONLY: allocate_subcell,&
      54              :                                               deallocate_subcell,&
      55              :                                               give_ijk_subcell,&
      56              :                                               reorder_atoms_subcell,&
      57              :                                               subcell_type
      58              :    USE util,                            ONLY: sort
      59              : #include "./base/base_uses.f90"
      60              : 
      61              :    IMPLICIT NONE
      62              : 
      63              :    PRIVATE
      64              : 
      65              :    ! Global parameters (in this module)
      66              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_neighbor_lists'
      67              : 
      68              :    TYPE local_atoms_type
      69              :       INTEGER, DIMENSION(:), POINTER                   :: list => NULL(), &
      70              :                                                           list_local_a_index => NULL()
      71              :    END TYPE local_atoms_type
      72              : 
      73              :    ! Public subroutines
      74              :    PUBLIC :: build_fist_neighbor_lists
      75              : 
      76              : CONTAINS
      77              : 
      78              : ! **************************************************************************************************
      79              : !> \brief ...
      80              : !> \param atomic_kind_set ...
      81              : !> \param particle_set ...
      82              : !> \param local_particles ...
      83              : !> \param cell ...
      84              : !> \param r_max ...
      85              : !> \param r_minsq ...
      86              : !> \param ei_scale14 ...
      87              : !> \param vdw_scale14 ...
      88              : !> \param nonbonded ...
      89              : !> \param para_env ...
      90              : !> \param build_from_scratch ...
      91              : !> \param geo_check ...
      92              : !> \param mm_section ...
      93              : !> \param full_nl ...
      94              : !> \param exclusions ...
      95              : !> \par History
      96              : !>      08.2006 created [tlaino]
      97              : !> \author Teodoro Laino
      98              : ! **************************************************************************************************
      99        19326 :    SUBROUTINE build_fist_neighbor_lists(atomic_kind_set, particle_set, &
     100        19326 :                                         local_particles, cell, r_max, r_minsq, ei_scale14, vdw_scale14, &
     101              :                                         nonbonded, para_env, build_from_scratch, geo_check, mm_section, &
     102        19326 :                                         full_nl, exclusions)
     103              : 
     104              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     105              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     106              :       TYPE(distribution_1d_type), OPTIONAL, POINTER      :: local_particles
     107              :       TYPE(cell_type), POINTER                           :: cell
     108              :       REAL(dp), DIMENSION(:, :), INTENT(IN)              :: r_max, r_minsq
     109              :       REAL(KIND=DP), INTENT(IN)                          :: ei_scale14, vdw_scale14
     110              :       TYPE(fist_neighbor_type), POINTER                  :: nonbonded
     111              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     112              :       LOGICAL, INTENT(IN)                                :: build_from_scratch, geo_check
     113              :       TYPE(section_vals_type), POINTER                   :: mm_section
     114              :       LOGICAL, DIMENSION(:, :), OPTIONAL, POINTER        :: full_nl
     115              :       TYPE(exclusion_type), DIMENSION(:), OPTIONAL       :: exclusions
     116              : 
     117              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_fist_neighbor_lists'
     118              : 
     119              :       CHARACTER(LEN=default_string_length)               :: kind_name, print_key_path, unit_str
     120              :       INTEGER                                            :: atom_a, handle, iatom_local, ikind, iw, &
     121              :                                                             maxatom, natom_local_a, nkind, &
     122              :                                                             output_unit
     123              :       LOGICAL                                            :: present_local_particles, &
     124              :                                                             print_subcell_grid
     125        19326 :       LOGICAL, DIMENSION(:), POINTER                     :: skip_kind
     126        19326 :       LOGICAL, DIMENSION(:, :), POINTER                  :: my_full_nl
     127              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     128              :       TYPE(cp_logger_type), POINTER                      :: logger
     129              :       TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:)  :: atom
     130              : 
     131        19326 :       CALL timeset(routineN, handle)
     132        19326 :       NULLIFY (logger)
     133        19326 :       logger => cp_get_default_logger()
     134              : 
     135        19326 :       print_subcell_grid = .FALSE.
     136              :       output_unit = cp_print_key_unit_nr(logger, mm_section, "PRINT%SUBCELL", &
     137        19326 :                                          extension=".Log")
     138        19326 :       IF (output_unit > 0) print_subcell_grid = .TRUE.
     139              : 
     140              :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
     141        19326 :                                maxatom=maxatom)
     142              : 
     143        19326 :       present_local_particles = PRESENT(local_particles)
     144              : 
     145              :       ! if exclusions matters local particles are present. Seems like only the exclusions
     146              :       ! for the local particles are needed, which would imply a huge memory savings for fist
     147        19326 :       IF (PRESENT(exclusions)) THEN
     148        10839 :          CPASSERT(present_local_particles)
     149              :       END IF
     150              : 
     151              :       ! Allocate work storage
     152        19326 :       nkind = SIZE(atomic_kind_set)
     153       118257 :       ALLOCATE (atom(nkind))
     154        57978 :       ALLOCATE (skip_kind(nkind))
     155              :       ! full_nl
     156        19326 :       IF (PRESENT(full_nl)) THEN
     157        11171 :          my_full_nl => full_nl
     158              :       ELSE
     159        32620 :          ALLOCATE (my_full_nl(nkind, nkind))
     160       496909 :          my_full_nl = .FALSE.
     161              :       END IF
     162              :       ! Initialize the local data structures
     163        79605 :       DO ikind = 1, nkind
     164        60279 :          atomic_kind => atomic_kind_set(ikind)
     165        60279 :          NULLIFY (atom(ikind)%list)
     166        60279 :          NULLIFY (atom(ikind)%list_local_a_index)
     167              : 
     168              :          CALL get_atomic_kind(atomic_kind=atomic_kind, &
     169        60279 :                               atom_list=atom(ikind)%list, name=kind_name)
     170        60279 :          skip_kind(ikind) = qmmm_ff_precond_only_qm(kind_name)
     171        60279 :          IF (present_local_particles) THEN
     172        50224 :             natom_local_a = local_particles%n_el(ikind)
     173              :          ELSE
     174        10055 :             natom_local_a = SIZE(atom(ikind)%list)
     175              :          END IF
     176        79605 :          IF (natom_local_a > 0) THEN
     177       167247 :             ALLOCATE (atom(ikind)%list_local_a_index(natom_local_a))
     178              :             ! Build index vector for mapping
     179      1075400 :             DO iatom_local = 1, natom_local_a
     180      1019651 :                IF (present_local_particles) THEN
     181       703427 :                   atom_a = local_particles%list(ikind)%array(iatom_local)
     182              :                ELSE
     183       316224 :                   atom_a = atom(ikind)%list(iatom_local)
     184              :                END IF
     185      1075400 :                atom(ikind)%list_local_a_index(iatom_local) = atom_a
     186              :             END DO
     187              : 
     188              :          END IF
     189              :       END DO
     190              : 
     191        19326 :       IF (build_from_scratch) THEN
     192        10660 :          IF (ASSOCIATED(nonbonded)) THEN
     193            0 :             CALL fist_neighbor_deallocate(nonbonded)
     194              :          END IF
     195              :       END IF
     196              : 
     197              :       ! Build the nonbonded neighbor lists
     198              :       CALL build_neighbor_lists(nonbonded, particle_set, atom, cell, &
     199              :                                 print_subcell_grid, output_unit, r_max, r_minsq, &
     200              :                                 ei_scale14, vdw_scale14, geo_check, "NONBONDED", skip_kind, &
     201        27813 :                                 my_full_nl, exclusions)
     202              : 
     203              :       ! Sort the list according kinds for each cell
     204        19326 :       CALL sort_neighbor_lists(nonbonded, nkind)
     205              : 
     206        19326 :       print_key_path = "PRINT%NEIGHBOR_LISTS"
     207              : 
     208        19326 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, mm_section, print_key_path), &
     209              :                 cp_p_file)) THEN
     210              :          iw = cp_print_key_unit_nr(logger=logger, &
     211              :                                    basis_section=mm_section, &
     212              :                                    print_key_path=print_key_path, &
     213              :                                    extension=".out", &
     214              :                                    middle_name="nonbonded_nl", &
     215              :                                    local=.TRUE., &
     216              :                                    log_filename=.FALSE., &
     217          460 :                                    file_position="REWIND")
     218          460 :          CALL section_vals_val_get(mm_section, TRIM(print_key_path)//"%UNIT", c_val=unit_str)
     219              :          CALL write_neighbor_lists(nonbonded, particle_set, cell, para_env, iw, &
     220          460 :                                    "NONBONDED NEIGHBOR LISTS", unit_str)
     221              :          CALL cp_print_key_finished_output(unit_nr=iw, &
     222              :                                            logger=logger, &
     223              :                                            basis_section=mm_section, &
     224              :                                            print_key_path=print_key_path, &
     225          460 :                                            local=.TRUE.)
     226              :       END IF
     227              : 
     228              :       ! Release work storage
     229        79605 :       DO ikind = 1, nkind
     230        60279 :          NULLIFY (atom(ikind)%list)
     231        79605 :          IF (ASSOCIATED(atom(ikind)%list_local_a_index)) THEN
     232        55749 :             DEALLOCATE (atom(ikind)%list_local_a_index)
     233              :          END IF
     234              :       END DO
     235        19326 :       IF (PRESENT(full_nl)) THEN
     236              :          NULLIFY (my_full_nl)
     237              :       ELSE
     238         8155 :          DEALLOCATE (my_full_nl)
     239              :       END IF
     240        19326 :       DEALLOCATE (atom)
     241              : 
     242        19326 :       DEALLOCATE (skip_kind)
     243              : 
     244              :       CALL cp_print_key_finished_output(unit_nr=output_unit, &
     245              :                                         logger=logger, &
     246              :                                         basis_section=mm_section, &
     247        19326 :                                         print_key_path="PRINT%SUBCELL")
     248        19326 :       CALL timestop(handle)
     249              : 
     250        38652 :    END SUBROUTINE build_fist_neighbor_lists
     251              : 
     252              : ! **************************************************************************************************
     253              : !> \brief ...
     254              : !> \param nonbonded ...
     255              : !> \param particle_set ...
     256              : !> \param atom ...
     257              : !> \param cell ...
     258              : !> \param print_subcell_grid ...
     259              : !> \param output_unit ...
     260              : !> \param r_max ...
     261              : !> \param r_minsq ...
     262              : !> \param ei_scale14 ...
     263              : !> \param vdw_scale14 ...
     264              : !> \param geo_check ...
     265              : !> \param name ...
     266              : !> \param skip_kind ...
     267              : !> \param full_nl ...
     268              : !> \param exclusions ...
     269              : !> \par History
     270              : !>      08.2006 created [tlaino]
     271              : !> \author Teodoro Laino
     272              : ! **************************************************************************************************
     273        19326 :    SUBROUTINE build_neighbor_lists(nonbonded, particle_set, atom, cell, &
     274        38652 :                                    print_subcell_grid, output_unit, r_max, r_minsq, &
     275        19326 :                                    ei_scale14, vdw_scale14, geo_check, name, skip_kind, full_nl, exclusions)
     276              : 
     277              :       TYPE(fist_neighbor_type), POINTER                  :: nonbonded
     278              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     279              :       TYPE(local_atoms_type), DIMENSION(:), INTENT(IN)   :: atom
     280              :       TYPE(cell_type), POINTER                           :: cell
     281              :       LOGICAL, INTENT(IN)                                :: print_subcell_grid
     282              :       INTEGER, INTENT(IN)                                :: output_unit
     283              :       REAL(dp), DIMENSION(:, :), INTENT(IN)              :: r_max, r_minsq
     284              :       REAL(KIND=dp), INTENT(IN)                          :: ei_scale14, vdw_scale14
     285              :       LOGICAL, INTENT(IN)                                :: geo_check
     286              :       CHARACTER(LEN=*), INTENT(IN)                       :: name
     287              :       LOGICAL, DIMENSION(:), POINTER                     :: skip_kind
     288              :       LOGICAL, DIMENSION(:, :), POINTER                  :: full_nl
     289              :       TYPE(exclusion_type), DIMENSION(:), OPTIONAL       :: exclusions
     290              : 
     291              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_neighbor_lists'
     292              : 
     293              :       INTEGER :: a_i, a_j, a_k, atom_a, atom_b, b_i, b_j, b_k, b_pi, b_pj, b_pk, bg_i, bg_j, bg_k, &
     294              :          handle, i, i1, iatom_local, icell, icellmap, id_kind, ii, ii_start, ij, ij_start, ik, &
     295              :          ik_start, ikind, imap, imax_cell, invcellmap, iw, ix, j, j1, jatom_local, jcell, jkind, &
     296              :          jx, k, kcell, kx, natom_local_a, ncellmax, nkind, nkind00, tmpdim, xdim, ydim, zdim
     297        19326 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of, work
     298        19326 :       INTEGER, ALLOCATABLE, DIMENSION(:, :, :)           :: cellmap
     299              :       INTEGER, DIMENSION(3)                              :: isubcell, ncell, nsubcell, periodic
     300              :       LOGICAL                                            :: any_full, atom_order, check_spline, &
     301              :                                                             is_full, subcell000
     302        19326 :       LOGICAL, ALLOCATABLE, DIMENSION(:, :, :)           :: sphcub
     303              :       REAL(dp)                                           :: rab2, rab2_max, rab2_min, rab_max
     304              :       REAL(dp), DIMENSION(3)                             :: abc, cell_v, cv_b, rab, rb, sab_max
     305              :       REAL(KIND=dp)                                      :: ic(3), icx(3), radius, vv
     306        19326 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: coord
     307              :       TYPE(neighbor_kind_pairs_type), POINTER            :: inv_neighbor_kind_pair, &
     308              :                                                             neighbor_kind_pair
     309        19326 :       TYPE(subcell_type), DIMENSION(:, :, :), POINTER    :: subcell_a, subcell_b
     310              : 
     311        19326 :       CALL timeset(routineN, handle)
     312        19326 :       nkind = SIZE(atom)
     313        77304 :       nsubcell = 1
     314        19326 :       isubcell = 0
     315        19326 :       ncell = 0
     316      1865161 :       any_full = ANY(full_nl)
     317              :       CALL get_cell(cell=cell, &
     318              :                     abc=abc, &
     319        19326 :                     periodic=periodic)
     320              :       ! Determines the number of subcells
     321        79605 :       DO ikind = 1, nkind
     322      1004169 :          DO jkind = ikind, nkind
     323              :             ! Calculate the square of the maximum interaction distance
     324       924564 :             rab_max = r_max(ikind, jkind)
     325       924564 :             IF (skip_kind(ikind) .AND. skip_kind(jkind)) CYCLE
     326       912106 :             nsubcell(1) = MAX(nsubcell(1), CEILING(plane_distance(1, 0, 0, cell)/rab_max))
     327       912106 :             nsubcell(2) = MAX(nsubcell(2), CEILING(plane_distance(0, 1, 0, cell)/rab_max))
     328       984843 :             nsubcell(3) = MAX(nsubcell(3), CEILING(plane_distance(0, 0, 1, cell)/rab_max))
     329              :          END DO
     330              :       END DO
     331              :       ! Determines the number of periodic images and the number of interacting subcells
     332        79605 :       DO ikind = 1, nkind
     333      1004169 :          DO jkind = ikind, nkind
     334       924564 :             IF (skip_kind(ikind) .AND. skip_kind(jkind)) CYCLE
     335              :             ! Calculate the square of the maximum interaction distance
     336       912106 :             rab_max = r_max(ikind, jkind)
     337       912106 :             sab_max(1) = rab_max/plane_distance(1, 0, 0, cell)
     338       912106 :             sab_max(2) = rab_max/plane_distance(0, 1, 0, cell)
     339       912106 :             sab_max(3) = rab_max/plane_distance(0, 0, 1, cell)
     340      3648424 :             ncell = MAX(ncell(:), CEILING(sab_max(:)*periodic(:)))
     341      3721161 :             isubcell = MAX(isubcell(:), CEILING(sab_max(:)*REAL(nsubcell(:), KIND=dp)))
     342              :          END DO
     343              :       END DO
     344        19326 :       CALL fist_neighbor_init(nonbonded, ncell)
     345              :       ! Print headline
     346        19326 :       IF (print_subcell_grid) THEN
     347              :          WRITE (UNIT=output_unit, FMT="(/,/,T2,A,/)") &
     348           86 :             "SUBCELL GRID  INFO FOR THE "//TRIM(name)//" NEIGHBOR LISTS"
     349           86 :          WRITE (UNIT=output_unit, FMT="(T4,A,10X,3I10)") " NUMBER OF SUBCELLS             ::", nsubcell
     350           86 :          WRITE (UNIT=output_unit, FMT="(T4,A,10X,3I10)") " NUMBER OF PERIODIC      IMAGES ::", ncell
     351           86 :          WRITE (UNIT=output_unit, FMT="(T4,A,10X,3I10)") " NUMBER OF INTERACTING SUBCELLS ::", isubcell
     352              :       END IF
     353              :       ! Allocate subcells
     354        19326 :       CALL allocate_subcell(subcell_a, nsubcell, cell=cell)
     355        19326 :       CALL allocate_subcell(subcell_b, nsubcell, cell=cell)
     356              :       ! Let's map the sequence of the periodic images
     357        77304 :       ncellmax = MAXVAL(ncell)
     358        96630 :       ALLOCATE (cellmap(-ncellmax:ncellmax, -ncellmax:ncellmax, -ncellmax:ncellmax))
     359      1115660 :       cellmap = -1
     360        19326 :       imap = 0
     361        19326 :       nkind00 = nkind*(nkind + 1)/2
     362        56197 :       DO imax_cell = 0, ncellmax
     363       134076 :          DO kcell = -imax_cell, imax_cell
     364       367661 :             DO jcell = -imax_cell, imax_cell
     365      1867845 :                DO icell = -imax_cell, imax_cell
     366      1789966 :                   IF (cellmap(icell, jcell, kcell) == -1) THEN
     367       858560 :                      imap = imap + 1
     368       858560 :                      cellmap(icell, jcell, kcell) = imap
     369       858560 :                      CPASSERT(imap <= nonbonded%nlists)
     370       858560 :                      neighbor_kind_pair => nonbonded%neighbor_kind_pairs(imap)
     371              : 
     372       858560 :                      neighbor_kind_pair%cell_vector(1) = icell
     373       858560 :                      neighbor_kind_pair%cell_vector(2) = jcell
     374       858560 :                      neighbor_kind_pair%cell_vector(3) = kcell
     375              :                   END IF
     376              :                END DO
     377              :             END DO
     378              :          END DO
     379              :       END DO
     380              :       ! Mapping the spherical interaction between subcells
     381              :       ALLOCATE (sphcub(-isubcell(1):isubcell(1), &
     382              :                        -isubcell(2):isubcell(2), &
     383        96630 :                        -isubcell(3):isubcell(3)))
     384      3113004 :       sphcub = .FALSE.
     385        77298 :       IF (ALL(isubcell /= 0)) THEN
     386              :          radius = REAL(isubcell(1), KIND=dp)**2 + REAL(isubcell(2), KIND=dp)**2 + &
     387        19324 :                   REAL(isubcell(3), KIND=dp)**2
     388       113340 :          loop1: DO k = -isubcell(3), isubcell(3)
     389       584628 :             loop2: DO j = -isubcell(2), isubcell(2)
     390      3093672 :                loop3: DO i = -isubcell(1), isubcell(1)
     391     10113472 :                   ic = REAL((/i, j, k/), KIND=dp)
     392              :                   ! subcell cube vertex
     393      3000034 :                   DO kx = -1, 1
     394      2528746 :                      icx(3) = ic(3) + SIGN(0.5_dp, REAL(kx, KIND=dp))
     395      2622052 :                      DO jx = -1, 1
     396      2621674 :                         icx(2) = ic(2) + SIGN(0.5_dp, REAL(jx, KIND=dp))
     397      3068784 :                         DO ix = -1, 1
     398      2975100 :                            icx(1) = ic(1) + SIGN(0.5_dp, REAL(ix, KIND=dp))
     399      2975100 :                            vv = icx(1)*icx(1) + icx(2)*icx(2) + icx(3)*icx(3)
     400      2975100 :                            vv = vv/radius
     401      3068406 :                            IF (vv <= 1.0_dp) THEN
     402      2528368 :                               sphcub(i, j, k) = .TRUE.
     403      2528368 :                               CYCLE loop3
     404              :                            END IF
     405              :                         END DO
     406              :                      END DO
     407              :                   END DO
     408              :                END DO loop3
     409              :             END DO loop2
     410              :          END DO loop1
     411              :       END IF
     412              :       ! Mapping locally all atoms in the zeroth cell
     413        57978 :       ALLOCATE (coord(3, SIZE(particle_set)))
     414      1501066 :       DO atom_a = 1, SIZE(particle_set)
     415      1501066 :          coord(:, atom_a) = pbc(particle_set(atom_a)%r, cell)
     416              :       END DO
     417              :       ! Associate particles to subcells (local particles)
     418        79605 :       DO ikind = 1, nkind
     419        60279 :          IF (.NOT. ASSOCIATED(atom(ikind)%list_local_a_index)) CYCLE
     420        55749 :          natom_local_a = SIZE(atom(ikind)%list_local_a_index)
     421      1094726 :          DO iatom_local = 1, natom_local_a
     422      1019651 :             atom_a = atom(ikind)%list_local_a_index(iatom_local)
     423      1019651 :             CALL give_ijk_subcell(coord(:, atom_a), i, j, k, cell, nsubcell)
     424      1079930 :             subcell_a(i, j, k)%natom = subcell_a(i, j, k)%natom + 1
     425              :          END DO
     426              :       END DO
     427        78480 :       DO k = 1, nsubcell(3)
     428       336705 :          DO j = 1, nsubcell(2)
     429      2272858 :             DO i = 1, nsubcell(1)
     430      4099529 :                ALLOCATE (subcell_a(i, j, k)%atom_list(subcell_a(i, j, k)%natom))
     431      2213704 :                subcell_a(i, j, k)%natom = 0
     432              :             END DO
     433              :          END DO
     434              :       END DO
     435        79605 :       DO ikind = 1, nkind
     436        60279 :          IF (.NOT. ASSOCIATED(atom(ikind)%list_local_a_index)) CYCLE
     437        55749 :          natom_local_a = SIZE(atom(ikind)%list_local_a_index)
     438      1094726 :          DO iatom_local = 1, natom_local_a
     439      1019651 :             atom_a = atom(ikind)%list_local_a_index(iatom_local)
     440      1019651 :             CALL give_ijk_subcell(coord(:, atom_a), i, j, k, cell, nsubcell)
     441      1019651 :             subcell_a(i, j, k)%natom = subcell_a(i, j, k)%natom + 1
     442      1079930 :             subcell_a(i, j, k)%atom_list(subcell_a(i, j, k)%natom) = atom_a
     443              :          END DO
     444              :       END DO
     445              :       ! Associate particles to subcells (distributed particles)
     446      1501066 :       DO atom_b = 1, SIZE(particle_set)
     447      1481740 :          CALL give_ijk_subcell(coord(:, atom_b), i, j, k, cell, nsubcell)
     448      1501066 :          subcell_b(i, j, k)%natom = subcell_b(i, j, k)%natom + 1
     449              :       END DO
     450        78480 :       DO k = 1, nsubcell(3)
     451       336705 :          DO j = 1, nsubcell(2)
     452      2272858 :             DO i = 1, nsubcell(1)
     453      4105981 :                ALLOCATE (subcell_b(i, j, k)%atom_list(subcell_b(i, j, k)%natom))
     454      2213704 :                subcell_b(i, j, k)%natom = 0
     455              :             END DO
     456              :          END DO
     457              :       END DO
     458      1501066 :       DO atom_b = 1, SIZE(particle_set)
     459      1481740 :          CALL give_ijk_subcell(coord(:, atom_b), i, j, k, cell, nsubcell)
     460      1481740 :          subcell_b(i, j, k)%natom = subcell_b(i, j, k)%natom + 1
     461      1501066 :          subcell_b(i, j, k)%atom_list(subcell_b(i, j, k)%natom) = atom_b
     462              :       END DO
     463              :       ! Reorder atoms associated to subcells
     464      2292184 :       tmpdim = MAXVAL(subcell_a(:, :, :)%natom)
     465      2292184 :       tmpdim = MAX(tmpdim, MAXVAL(subcell_b(:, :, :)%natom))
     466        57978 :       ALLOCATE (work(3*tmpdim))
     467        57978 :       ALLOCATE (kind_of(SIZE(particle_set)))
     468      1501066 :       DO i = 1, SIZE(particle_set)
     469      1501066 :          kind_of(i) = particle_set(i)%atomic_kind%kind_number
     470              :       END DO
     471        78480 :       DO k = 1, nsubcell(3)
     472       336705 :          DO j = 1, nsubcell(2)
     473      2272858 :             DO i = 1, nsubcell(1)
     474      1955479 :                CALL reorder_atoms_subcell(subcell_a(i, j, k)%atom_list, kind_of, work)
     475      2213704 :                CALL reorder_atoms_subcell(subcell_b(i, j, k)%atom_list, kind_of, work)
     476              :             END DO
     477              :          END DO
     478              :       END DO
     479        19326 :       DEALLOCATE (work, kind_of)
     480        19326 :       zdim = nsubcell(3)
     481        19326 :       ydim = nsubcell(2)
     482        19326 :       xdim = nsubcell(1)
     483        19326 :       is_full = .FALSE.
     484              :       ! We can skip until ik>=0.. this prescreens the order of the subcells
     485        19326 :       ik_start = -isubcell(3)
     486        19326 :       IF (.NOT. any_full) ik_start = 0
     487              :       ! Loop over first subcell
     488        78480 :       loop_a_k: DO a_k = 1, nsubcell(3)
     489       336705 :       loop_a_j: DO a_j = 1, nsubcell(2)
     490      2272858 :       loop_a_i: DO a_i = 1, nsubcell(1)
     491      1955479 :          IF (subcell_a(a_i, a_j, a_k)%natom == 0) CYCLE
     492              :          ! Loop over second subcell
     493      1008890 :          loop_b_k: DO ik = ik_start, isubcell(3)
     494       562094 :             bg_k = a_k + ik
     495       562094 :             b_k = MOD(bg_k, zdim)
     496       562094 :             b_pk = bg_k/zdim
     497       562094 :             IF (b_k <= 0) THEN
     498       123788 :                b_k = zdim + b_k
     499       123788 :                b_pk = b_pk - 1
     500              :             END IF
     501       562094 :             IF ((periodic(3) == 0) .AND. (ABS(b_pk) > 0)) CYCLE
     502              :             ! Setup the starting point.. this prescreens the order of the subcells
     503       554031 :             ij_start = -isubcell(2)
     504       554031 :             IF ((ik == 0) .AND. (ik_start == 0)) ij_start = 0
     505      4893012 :             loop_b_j: DO ij = ij_start, isubcell(2)
     506      2383502 :                bg_j = a_j + ij
     507      2383502 :                b_j = MOD(bg_j, ydim)
     508      2383502 :                b_pj = bg_j/ydim
     509      2383502 :                IF (b_j <= 0) THEN
     510       600156 :                   b_j = ydim + b_j
     511       600156 :                   b_pj = b_pj - 1
     512              :                END IF
     513      2383502 :                IF ((periodic(2) == 0) .AND. (ABS(b_pj) > 0)) CYCLE
     514              :                ! Setup the starting point.. this prescreens the order of the subcells
     515      2367348 :                ii_start = -isubcell(1)
     516      2367348 :                IF ((ij == 0) .AND. (ij_start == 0)) ii_start = 0
     517     14438039 :                loop_b_i: DO ii = ii_start, isubcell(1)
     518              :                   ! Ellipsoidal screening of subcells
     519     11508597 :                   IF (.NOT. sphcub(ii, ij, ik)) CYCLE
     520     11508596 :                   bg_i = a_i + ii
     521     11508596 :                   b_i = MOD(bg_i, xdim)
     522     11508596 :                   b_pi = bg_i/xdim
     523     11508596 :                   IF (b_i <= 0) THEN
     524      3077687 :                      b_i = xdim + b_i
     525      3077687 :                      b_pi = b_pi - 1
     526              :                   END IF
     527     11508596 :                   IF ((periodic(1) == 0) .AND. (ABS(b_pi) > 0)) CYCLE
     528     11470469 :                   IF (subcell_b(b_i, b_j, b_k)%natom == 0) CYCLE
     529              :                   ! Find the proper neighbor kind pair
     530      8106774 :                   icellmap = cellmap(b_pi, b_pj, b_pk)
     531      8106774 :                   neighbor_kind_pair => nonbonded%neighbor_kind_pairs(icellmap)
     532              :                   ! Find the replica vector
     533      8106774 :                   cell_v = 0.0_dp
     534      8106774 :                   IF ((b_pi /= 0) .OR. (b_pj /= 0) .OR. (b_pk /= 0)) THEN
     535      3854334 :                      cv_b(1) = b_pi; cv_b(2) = b_pj; cv_b(3) = b_pk
     536      3854334 :                      CALL scaled_to_real(cell_v, cv_b, cell)
     537              :                   END IF
     538      8106774 :                   subcell000 = (a_k == bg_k) .AND. (a_j == bg_j) .AND. (a_i == bg_i)
     539              :                   ! Loop over particles inside subcell_a and subcell_b
     540     93317802 :                   DO jatom_local = 1, subcell_b(b_i, b_j, b_k)%natom
     541     82827526 :                      atom_b = subcell_b(b_i, b_j, b_k)%atom_list(jatom_local)
     542     82827526 :                      jkind = particle_set(atom_b)%atomic_kind%kind_number
     543     82827526 :                      rb(1) = coord(1, atom_b) + cell_v(1)
     544     82827526 :                      rb(2) = coord(2, atom_b) + cell_v(2)
     545     82827526 :                      rb(3) = coord(3, atom_b) + cell_v(3)
     546   2488256709 :                      DO iatom_local = 1, subcell_a(a_i, a_j, a_k)%natom
     547   2393920586 :                         atom_a = subcell_a(a_i, a_j, a_k)%atom_list(iatom_local)
     548   2393920586 :                         ikind = particle_set(atom_a)%atomic_kind%kind_number
     549              :                         ! Screen interaction to avoid double counting
     550   2393920586 :                         atom_order = (atom_a <= atom_b)
     551              :                         ! Special case for kind combination requiring the full NL
     552   2393920586 :                         IF (any_full) THEN
     553     19336474 :                            is_full = full_nl(ikind, jkind)
     554     19336474 :                            IF (is_full) THEN
     555      9066349 :                               atom_order = (atom_a == atom_b)
     556              :                            ELSE
     557     10270125 :                               IF (ik < 0) CYCLE
     558      6900294 :                               IF (ik == 0 .AND. ij < 0) CYCLE
     559      5491910 :                               IF (ij == 0 .AND. ii < 0) CYCLE
     560              :                            END IF
     561              :                         END IF
     562   2388586438 :                         IF (subcell000 .AND. atom_order) CYCLE
     563   2368390332 :                         rab(1) = rb(1) - coord(1, atom_a)
     564   2368390332 :                         rab(2) = rb(2) - coord(2, atom_a)
     565   2368390332 :                         rab(3) = rb(3) - coord(3, atom_a)
     566   2368390332 :                         rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
     567   2368390332 :                         rab_max = r_max(ikind, jkind)
     568   2368390332 :                         rab2_max = rab_max*rab_max
     569   2451217858 :                         IF (rab2 < rab2_max) THEN
     570              :                            ! Diagonal storage
     571    151495343 :                            j1 = MIN(ikind, jkind)
     572    151495343 :                            i1 = MAX(ikind, jkind) - j1 + 1
     573    151495343 :                            j1 = nkind - j1 + 1
     574    151495343 :                            id_kind = nkind00 - (j1*(j1 + 1)/2) + i1
     575              :                            ! Store the pair
     576              :                            CALL fist_neighbor_add(neighbor_kind_pair, atom_a, atom_b, &
     577              :                                                   rab=rab, &
     578              :                                                   check_spline=check_spline, id_kind=id_kind, &
     579              :                                                   skip=(skip_kind(ikind) .AND. skip_kind(jkind)), &
     580              :                                                   cell=cell, ei_scale14=ei_scale14, &
     581    310480404 :                                                   vdw_scale14=vdw_scale14, exclusions=exclusions)
     582              :                            ! This is to handle properly when interaction radius is larger than cell size
     583    151495343 :                            IF ((atom_a == atom_b) .AND. (ik_start == 0)) THEN
     584       352913 :                               invcellmap = cellmap(-b_pi, -b_pj, -b_pk)
     585       352913 :                               inv_neighbor_kind_pair => nonbonded%neighbor_kind_pairs(invcellmap)
     586      1411652 :                               rab = rab - 2.0_dp*cell_v
     587              :                               CALL fist_neighbor_add(inv_neighbor_kind_pair, atom_a, atom_b, &
     588              :                                                      rab=rab, &
     589              :                                                      check_spline=check_spline, id_kind=id_kind, &
     590              :                                                      skip=(skip_kind(ikind) .AND. skip_kind(jkind)), &
     591              :                                                      cell=cell, ei_scale14=ei_scale14, &
     592       725693 :                                                      vdw_scale14=vdw_scale14, exclusions=exclusions)
     593              :                            END IF
     594              :                            ! Check for too close hits
     595    151495343 :                            IF (check_spline) THEN
     596    150834642 :                               rab2_min = r_minsq(ikind, jkind)
     597    150834642 :                               IF (rab2 < rab2_min) THEN
     598            1 :                                  iw = cp_logger_get_default_unit_nr()
     599            1 :                                  WRITE (iw, '(T2,A,2I7,2(A,F15.8),A)') "WARNING| Particles: ", &
     600            1 :                                     atom_a, atom_b, &
     601            1 :                                     " at distance [au]:", SQRT(rab2), " less than: ", &
     602            1 :                                     SQRT(rab2_min), &
     603            2 :                                     "; increase EMAX_SPLINE."
     604            1 :                                  IF (rab2 < rab2_min/(1.06_dp)**2) THEN
     605            0 :                                     IF (geo_check) THEN
     606            0 :                                        CPABORT("GEOMETRY wrong or EMAX_SPLINE too small!")
     607              :                                     END IF
     608              :                                  END IF
     609              :                               END IF
     610              :                            END IF
     611              :                         END IF
     612              :                      END DO
     613              :                   END DO
     614              :                END DO loop_b_i
     615              :             END DO loop_b_j
     616              :          END DO loop_b_k
     617              :       END DO loop_a_i
     618              :       END DO loop_a_j
     619              :       END DO loop_a_k
     620        19326 :       DEALLOCATE (coord)
     621        19326 :       DEALLOCATE (cellmap)
     622        19326 :       DEALLOCATE (sphcub)
     623        19326 :       CALL deallocate_subcell(subcell_a)
     624        19326 :       CALL deallocate_subcell(subcell_b)
     625              : 
     626        19326 :       CALL timestop(handle)
     627        19326 :    END SUBROUTINE build_neighbor_lists
     628              : 
     629              : ! **************************************************************************************************
     630              : !> \brief Write a set of neighbor lists to the output unit.
     631              : !> \param nonbonded ...
     632              : !> \param particle_set ...
     633              : !> \param cell ...
     634              : !> \param para_env ...
     635              : !> \param output_unit ...
     636              : !> \param name ...
     637              : !> \param unit_str ...
     638              : !> \par History
     639              : !>      08.2006 created [tlaino]
     640              : !> \author Teodoro Laino
     641              : ! **************************************************************************************************
     642          460 :    SUBROUTINE write_neighbor_lists(nonbonded, particle_set, cell, para_env, output_unit, &
     643              :                                    name, unit_str)
     644              : 
     645              :       TYPE(fist_neighbor_type), POINTER                  :: nonbonded
     646              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     647              :       TYPE(cell_type), POINTER                           :: cell
     648              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     649              :       INTEGER, INTENT(IN)                                :: output_unit
     650              :       CHARACTER(LEN=*), INTENT(IN)                       :: name, unit_str
     651              : 
     652              :       CHARACTER(LEN=default_string_length)               :: string
     653              :       INTEGER                                            :: atom_a, atom_b, iab, ilist, nneighbor
     654              :       LOGICAL                                            :: print_headline
     655              :       REAL(dp)                                           :: conv, dab
     656              :       REAL(dp), DIMENSION(3)                             :: cell_v, ra, rab, rb
     657              :       TYPE(neighbor_kind_pairs_type), POINTER            :: neighbor_kind_pair
     658              : 
     659              :       ! Print headline
     660          460 :       string = ""
     661              :       WRITE (UNIT=string, FMT="(A,I5,A)") &
     662          460 :          TRIM(name)//" IN "//TRIM(unit_str)//" (PROCESS", para_env%mepos, ")"
     663          460 :       CALL compress(string)
     664          460 :       IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") TRIM(string)
     665              : 
     666          460 :       print_headline = .TRUE.
     667          460 :       nneighbor = 0
     668          460 :       conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
     669        22116 :       DO iab = 1, SIZE(nonbonded%neighbor_kind_pairs)
     670        21656 :          neighbor_kind_pair => nonbonded%neighbor_kind_pairs(iab)
     671       346496 :          cell_v = MATMUL(cell%hmat, REAL(neighbor_kind_pair%cell_vector, KIND=dp))
     672       274699 :          DO ilist = 1, neighbor_kind_pair%npairs
     673       252583 :             nneighbor = nneighbor + 1
     674       274239 :             IF (output_unit > 0) THEN
     675              :                ! Print second part of the headline
     676       252583 :                atom_a = neighbor_kind_pair%list(1, ilist)
     677       252583 :                atom_b = neighbor_kind_pair%list(2, ilist)
     678       252583 :                IF (print_headline) THEN
     679              :                   WRITE (UNIT=output_unit, FMT="(T3,2(A6,3(5X,A,5X)),1X,A11,10X,A8,A5,A10,A9)") &
     680          294 :                      "Atom-A", "X", "Y", "Z", "Atom-B", "X", "Y", "Z", "Cell(i,j,k)", &
     681          588 :                      "Distance", "ONFO", "VDW-scale", "EI-scale"
     682          294 :                   print_headline = .FALSE.
     683              :                END IF
     684              : 
     685       252583 :                ra(:) = pbc(particle_set(atom_a)%r, cell)
     686       252583 :                rb(:) = pbc(particle_set(atom_b)%r, cell)
     687      1010332 :                rab = rb(:) - ra(:) + cell_v
     688      1010332 :                dab = SQRT(DOT_PRODUCT(rab, rab))
     689       252583 :                IF (ilist <= neighbor_kind_pair%nscale) THEN
     690              :                   WRITE (UNIT=output_unit, FMT="(T3,2(I6,3(1X,F10.6)),3(1X,I3),10X,F8.4,L4,F11.5,F9.5)") &
     691        27248 :                      atom_a, ra(1:3)*conv, &
     692        27248 :                      atom_b, rb(1:3)*conv, &
     693         6812 :                      neighbor_kind_pair%cell_vector, &
     694         6812 :                      dab*conv, &
     695         6812 :                      neighbor_kind_pair%is_onfo(ilist), &
     696         6812 :                      neighbor_kind_pair%vdw_scale(ilist), &
     697        13624 :                      neighbor_kind_pair%ei_scale(ilist)
     698              :                ELSE
     699              :                   WRITE (UNIT=output_unit, FMT="(T3,2(I6,3(1X,F10.6)),3(1X,I3),10X,F8.4)") &
     700       983084 :                      atom_a, ra(1:3)*conv, &
     701       983084 :                      atom_b, rb(1:3)*conv, &
     702       245771 :                      neighbor_kind_pair%cell_vector, &
     703       491542 :                      dab*conv
     704              :                END IF
     705              :             END IF
     706              :          END DO ! ilist
     707              :       END DO ! iab
     708              : 
     709          460 :       string = ""
     710              :       WRITE (UNIT=string, FMT="(A,I12,A,I12)") &
     711          460 :          "Total number of neighbor interactions for process", para_env%mepos, ":", &
     712          920 :          nneighbor
     713          460 :       CALL compress(string)
     714          460 :       IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="(/,T2,A)") TRIM(string)
     715              : 
     716          460 :    END SUBROUTINE write_neighbor_lists
     717              : 
     718              : ! **************************************************************************************************
     719              : !> \brief Sort the generated neighbor list according the kind
     720              : !> \param nonbonded ...
     721              : !> \param nkinds ...
     722              : !> \par History
     723              : !>      09.2007 created [tlaino] University of Zurich - Reducing memory usage
     724              : !>              for the FIST neighbor lists
     725              : !> \author Teodoro Laino - University of Zurich
     726              : ! **************************************************************************************************
     727        19326 :    SUBROUTINE sort_neighbor_lists(nonbonded, nkinds)
     728              : 
     729              :       TYPE(fist_neighbor_type), POINTER                  :: nonbonded
     730              :       INTEGER, INTENT(IN)                                :: nkinds
     731              : 
     732              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'sort_neighbor_lists'
     733              : 
     734              :       INTEGER                                            :: handle, iab, id_kind, ikind, ipair, &
     735              :                                                             jkind, max_alloc_size, npairs, nscale, &
     736              :                                                             tmp
     737        19326 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: indj
     738        19326 :       INTEGER, DIMENSION(:), POINTER                     :: work
     739        19326 :       INTEGER, DIMENSION(:, :), POINTER                  :: list_copy
     740              :       TYPE(neighbor_kind_pairs_type), POINTER            :: neighbor_kind_pair
     741              : 
     742        19326 :       NULLIFY (neighbor_kind_pair)
     743        19326 :       CALL timeset(routineN, handle)
     744              :       ! define a lookup table to get jkind for a given id_kind
     745        57978 :       ALLOCATE (indj(nkinds*(nkinds + 1)/2))
     746        19326 :       id_kind = 0
     747        79605 :       DO jkind = 1, nkinds
     748      1004169 :          DO ikind = jkind, nkinds
     749       924564 :             id_kind = id_kind + 1
     750       984843 :             indj(id_kind) = jkind
     751              :          END DO
     752              :       END DO
     753              :       ! loop over all nlists and sort the pairs within each list.
     754       877886 :       DO iab = 1, nonbonded%nlists
     755       858560 :          neighbor_kind_pair => nonbonded%neighbor_kind_pairs(iab)
     756       858560 :          npairs = neighbor_kind_pair%npairs
     757       858560 :          nscale = neighbor_kind_pair%nscale
     758       858560 :          IF (npairs /= 0) THEN
     759       193417 :             IF (npairs > nscale) THEN
     760              :                ! 1) Sort the atom pairs by id_kind. Pairs whose interactions are
     761              :                ! scaled (possibly to zero for exclusion) are not touched. They
     762              :                ! stay packed in the beginning. Sorting is skipped altogether when
     763              :                ! all pairs have scaled interactions.
     764       574512 :                ALLOCATE (work(1:npairs - nscale))
     765       574512 :                ALLOCATE (list_copy(2, 1:npairs - nscale))
     766              :                ! Copy of the pair list is required to perform the permutation below
     767              :                ! correctly.
     768    906571516 :                list_copy = neighbor_kind_pair%list(:, nscale + 1:npairs)
     769       191504 :                CALL sort(neighbor_kind_pair%id_kind(nscale + 1:npairs), npairs - nscale, work)
     770              :                ! Reorder atoms using the same permutation that was used to sort
     771              :                ! the array id_kind.
     772    151222922 :                DO ipair = nscale + 1, npairs
     773    151031418 :                   tmp = work(ipair - nscale)
     774    151031418 :                   neighbor_kind_pair%list(1, ipair) = list_copy(1, tmp)
     775    151222922 :                   neighbor_kind_pair%list(2, ipair) = list_copy(2, tmp)
     776              :                END DO
     777       191504 :                DEALLOCATE (work)
     778       191504 :                DEALLOCATE (list_copy)
     779              :             END IF
     780              :             ! 2) determine the intervals (groups) in the pair list that correspond
     781              :             !    to a certain id_kind. also store the corresponding ikind and
     782              :             !    jkind. Note that this part does not assume ikind to be sorted,
     783              :             !    but it only makes sense when contiguous blobs of the same ikind
     784              :             !    are present.
     785              :             ! Allocate sufficient memory in case all pairs of atom kinds are
     786              :             ! present, and also provide storage for the pairs with exclusion
     787              :             ! flags, which are unsorted.
     788       193417 :             max_alloc_size = nkinds*(nkinds + 1)/2 + nscale
     789       193417 :             IF (ASSOCIATED(neighbor_kind_pair%grp_kind_start)) THEN
     790        59769 :                DEALLOCATE (neighbor_kind_pair%grp_kind_start)
     791              :             END IF
     792       580251 :             ALLOCATE (neighbor_kind_pair%grp_kind_start(max_alloc_size))
     793       193417 :             IF (ASSOCIATED(neighbor_kind_pair%grp_kind_end)) THEN
     794        59769 :                DEALLOCATE (neighbor_kind_pair%grp_kind_end)
     795              :             END IF
     796       386834 :             ALLOCATE (neighbor_kind_pair%grp_kind_end(max_alloc_size))
     797       193417 :             IF (ASSOCIATED(neighbor_kind_pair%ij_kind)) THEN
     798        59769 :                DEALLOCATE (neighbor_kind_pair%ij_kind)
     799              :             END IF
     800       580251 :             ALLOCATE (neighbor_kind_pair%ij_kind(2, max_alloc_size))
     801              :             ! Start the first interval.
     802       193417 :             ipair = 1
     803       193417 :             neighbor_kind_pair%ngrp_kind = 1
     804       193417 :             neighbor_kind_pair%grp_kind_start(neighbor_kind_pair%ngrp_kind) = ipair
     805              :             ! Get ikind and jkind corresponding to id_kind.
     806       193417 :             id_kind = neighbor_kind_pair%id_kind(ipair)
     807       193417 :             jkind = indj(id_kind)
     808       193417 :             tmp = nkinds - jkind
     809       193417 :             ikind = nkinds + id_kind - nkinds*(nkinds + 1)/2 + (tmp*(tmp + 1)/2)
     810       193417 :             neighbor_kind_pair%ij_kind(1, neighbor_kind_pair%ngrp_kind) = ikind
     811       193417 :             neighbor_kind_pair%ij_kind(2, neighbor_kind_pair%ngrp_kind) = jkind
     812              :             ! Define the remaining intervals.
     813    151821095 :             DO ipair = 2, npairs
     814    151821095 :                IF (neighbor_kind_pair%id_kind(ipair) /= neighbor_kind_pair%id_kind(ipair - 1)) THEN
     815      1116612 :                   neighbor_kind_pair%grp_kind_end(neighbor_kind_pair%ngrp_kind) = ipair - 1
     816      1116612 :                   neighbor_kind_pair%ngrp_kind = neighbor_kind_pair%ngrp_kind + 1
     817      1116612 :                   neighbor_kind_pair%grp_kind_start(neighbor_kind_pair%ngrp_kind) = ipair
     818              :                   ! Get ikind and jkind corresponding to id_kind.
     819      1116612 :                   id_kind = neighbor_kind_pair%id_kind(ipair)
     820      1116612 :                   jkind = indj(id_kind)
     821      1116612 :                   tmp = nkinds - jkind
     822      1116612 :                   ikind = nkinds + id_kind - nkinds*(nkinds + 1)/2 + (tmp*(tmp + 1)/2)
     823      1116612 :                   neighbor_kind_pair%ij_kind(1, neighbor_kind_pair%ngrp_kind) = ikind
     824      1116612 :                   neighbor_kind_pair%ij_kind(2, neighbor_kind_pair%ngrp_kind) = jkind
     825              :                END IF
     826              :             END DO
     827              :             ! Finish the last interval.
     828       193417 :             neighbor_kind_pair%grp_kind_end(neighbor_kind_pair%ngrp_kind) = npairs
     829              :             ! Reduce the grp arrays to the actual size because not all pairs of
     830              :             ! atom types have to be present in this pair list.
     831       193417 :             CALL reallocate(neighbor_kind_pair%grp_kind_start, 1, neighbor_kind_pair%ngrp_kind)
     832       193417 :             CALL reallocate(neighbor_kind_pair%grp_kind_end, 1, neighbor_kind_pair%ngrp_kind)
     833       193417 :             CALL reallocate(neighbor_kind_pair%ij_kind, 1, 2, 1, neighbor_kind_pair%ngrp_kind)
     834              :          END IF
     835              :          ! Clean the memory..
     836       877886 :          DEALLOCATE (neighbor_kind_pair%id_kind)
     837              :       END DO
     838        19326 :       DEALLOCATE (indj)
     839        19326 :       CALL timestop(handle)
     840        19326 :    END SUBROUTINE sort_neighbor_lists
     841              : 
     842            0 : END MODULE fist_neighbor_lists
        

Generated by: LCOV version 2.0-1