LCOV - code coverage report
Current view: top level - src - qs_neighbor_lists.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.7 % 838 810
Test Date: 2025-07-25 12:55:17 Functions: 90.0 % 10 9

            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.
      10              : !> \par History
      11              : !>      - List rebuild for sab_orb neighbor list (10.09.2002,MK)
      12              : !>      - List rebuild for all lists (25.09.2002,MK)
      13              : !>      - Row-wise parallelized version (16.06.2003,MK)
      14              : !>      - Row- and column-wise parallelized version (19.07.2003,MK)
      15              : !>      - bug fix for non-periodic case (23.02.06,MK)
      16              : !>      - major refactoring (25.07.10,jhu)
      17              : !> \author Matthias Krack (08.10.1999,26.03.2002,16.06.2003)
      18              : ! **************************************************************************************************
      19              : MODULE qs_neighbor_lists
      20              :    USE almo_scf_types,                  ONLY: almo_max_cutoff_multiplier
      21              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      22              :                                               get_atomic_kind,&
      23              :                                               get_atomic_kind_set
      24              :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      25              :                                               gto_basis_set_p_type,&
      26              :                                               gto_basis_set_type
      27              :    USE cell_types,                      ONLY: cell_type,&
      28              :                                               get_cell,&
      29              :                                               pbc,&
      30              :                                               plane_distance,&
      31              :                                               real_to_scaled,&
      32              :                                               scaled_to_real
      33              :    USE cp_control_types,                ONLY: dft_control_type
      34              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      35              :                                               cp_logger_type
      36              :    USE cp_output_handling,              ONLY: cp_p_file,&
      37              :                                               cp_print_key_finished_output,&
      38              :                                               cp_print_key_should_output,&
      39              :                                               cp_print_key_unit_nr
      40              :    USE cp_units,                        ONLY: cp_unit_from_cp2k
      41              :    USE distribution_1d_types,           ONLY: distribution_1d_type
      42              :    USE distribution_2d_types,           ONLY: distribution_2d_type
      43              :    USE ewald_environment_types,         ONLY: ewald_env_get,&
      44              :                                               ewald_environment_type
      45              :    USE external_potential_types,        ONLY: all_potential_type,&
      46              :                                               get_potential,&
      47              :                                               gth_potential_type,&
      48              :                                               sgp_potential_type
      49              :    USE input_constants,                 ONLY: &
      50              :         dispersion_uff, do_method_lrigpw, do_method_rigpw, do_potential_id, &
      51              :         do_potential_mix_cl_trunc, do_potential_short, do_potential_truncated, do_se_IS_slater, &
      52              :         vdw_pairpot_dftd4, xc_vdw_fun_pairpot
      53              :    USE input_section_types,             ONLY: section_vals_get,&
      54              :                                               section_vals_get_subs_vals,&
      55              :                                               section_vals_type,&
      56              :                                               section_vals_val_get
      57              :    USE kinds,                           ONLY: default_string_length,&
      58              :                                               dp,&
      59              :                                               int_8
      60              :    USE kpoint_types,                    ONLY: kpoint_type
      61              :    USE libint_2c_3c,                    ONLY: cutoff_screen_factor
      62              :    USE mathlib,                         ONLY: erfc_cutoff
      63              :    USE message_passing,                 ONLY: mp_para_env_type
      64              :    USE molecule_types,                  ONLY: molecule_type
      65              :    USE particle_types,                  ONLY: particle_type
      66              :    USE paw_proj_set_types,              ONLY: get_paw_proj_set,&
      67              :                                               paw_proj_set_type
      68              :    USE periodic_table,                  ONLY: ptable
      69              :    USE physcon,                         ONLY: bohr
      70              :    USE qs_dftb_types,                   ONLY: qs_dftb_atom_type
      71              :    USE qs_dftb_utils,                   ONLY: get_dftb_atom_param
      72              :    USE qs_dispersion_types,             ONLY: qs_dispersion_type
      73              :    USE qs_environment_types,            ONLY: get_qs_env,&
      74              :                                               qs_environment_type
      75              :    USE qs_gcp_types,                    ONLY: qs_gcp_type
      76              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      77              :                                               get_qs_kind_set,&
      78              :                                               qs_kind_type
      79              :    USE qs_ks_types,                     ONLY: get_ks_env,&
      80              :                                               qs_ks_env_type,&
      81              :                                               set_ks_env
      82              :    USE qs_neighbor_list_types,          ONLY: &
      83              :         add_neighbor_list, add_neighbor_node, allocate_neighbor_list_set, get_iterator_info, &
      84              :         get_iterator_task, neighbor_list_iterate, neighbor_list_iterator_create, &
      85              :         neighbor_list_iterator_p_type, neighbor_list_iterator_release, neighbor_list_p_type, &
      86              :         neighbor_list_set_p_type, neighbor_list_set_type, release_neighbor_list_sets
      87              :    USE string_utilities,                ONLY: compress,&
      88              :                                               uppercase
      89              :    USE subcell_types,                   ONLY: allocate_subcell,&
      90              :                                               deallocate_subcell,&
      91              :                                               give_ijk_subcell,&
      92              :                                               subcell_type
      93              :    USE util,                            ONLY: locate,&
      94              :                                               sort
      95              :    USE xtb_types,                       ONLY: get_xtb_atom_param,&
      96              :                                               xtb_atom_type
      97              : #include "./base/base_uses.f90"
      98              : 
      99              :    IMPLICIT NONE
     100              : 
     101              :    PRIVATE
     102              : 
     103              : ! **************************************************************************************************
     104              :    TYPE local_atoms_type
     105              :       INTEGER, DIMENSION(:), POINTER                   :: list => NULL(), &
     106              :                                                           list_local_a_index => NULL(), &
     107              :                                                           list_local_b_index => NULL(), &
     108              :                                                           list_1d => NULL(), &
     109              :                                                           list_a_mol => NULL(), &
     110              :                                                           list_b_mol => NULL()
     111              :    END TYPE local_atoms_type
     112              : ! **************************************************************************************************
     113              : 
     114              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_neighbor_lists'
     115              : 
     116              :    ! private counter, used to version qs neighbor lists
     117              :    INTEGER, SAVE, PRIVATE :: last_qs_neighbor_list_id_nr = 0
     118              : 
     119              :    ! Public subroutines
     120              :    PUBLIC :: build_qs_neighbor_lists, local_atoms_type, atom2d_cleanup, &
     121              :              atom2d_build, build_neighbor_lists, pair_radius_setup, &
     122              :              setup_neighbor_list, write_neighbor_lists
     123              : CONTAINS
     124              : 
     125              : ! **************************************************************************************************
     126              : !> \brief   free the internals of atom2d
     127              : !> \param atom2d ...
     128              : !> \param
     129              : ! **************************************************************************************************
     130        41123 :    SUBROUTINE atom2d_cleanup(atom2d)
     131              :       TYPE(local_atoms_type), DIMENSION(:)               :: atom2d
     132              : 
     133              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'atom2d_cleanup'
     134              : 
     135              :       INTEGER                                            :: handle, ikind
     136              : 
     137        41123 :       CALL timeset(routineN, handle)
     138       118267 :       DO ikind = 1, SIZE(atom2d)
     139        77144 :          NULLIFY (atom2d(ikind)%list)
     140        77144 :          IF (ASSOCIATED(atom2d(ikind)%list_local_a_index)) THEN
     141        55557 :             DEALLOCATE (atom2d(ikind)%list_local_a_index)
     142              :          END IF
     143        77144 :          IF (ASSOCIATED(atom2d(ikind)%list_local_b_index)) THEN
     144        77086 :             DEALLOCATE (atom2d(ikind)%list_local_b_index)
     145              :          END IF
     146        77144 :          IF (ASSOCIATED(atom2d(ikind)%list_a_mol)) THEN
     147        55557 :             DEALLOCATE (atom2d(ikind)%list_a_mol)
     148              :          END IF
     149        77144 :          IF (ASSOCIATED(atom2d(ikind)%list_b_mol)) THEN
     150        77086 :             DEALLOCATE (atom2d(ikind)%list_b_mol)
     151              :          END IF
     152       118267 :          IF (ASSOCIATED(atom2d(ikind)%list_1d)) THEN
     153        77144 :             DEALLOCATE (atom2d(ikind)%list_1d)
     154              :          END IF
     155              :       END DO
     156        41123 :       CALL timestop(handle)
     157              : 
     158        41123 :    END SUBROUTINE atom2d_cleanup
     159              : 
     160              : ! **************************************************************************************************
     161              : !> \brief   Build some distribution structure of atoms, refactored from build_qs_neighbor_lists
     162              : !> \param atom2d output
     163              : !> \param distribution_1d ...
     164              : !> \param distribution_2d ...
     165              : !> \param atomic_kind_set ...
     166              : !> \param molecule_set ...
     167              : !> \param molecule_only ...
     168              : !> \param particle_set ...
     169              : !> \author  JH
     170              : ! **************************************************************************************************
     171        41123 :    SUBROUTINE atom2d_build(atom2d, distribution_1d, distribution_2d, &
     172              :                            atomic_kind_set, molecule_set, molecule_only, particle_set)
     173              :       TYPE(local_atoms_type), DIMENSION(:)               :: atom2d
     174              :       TYPE(distribution_1d_type), POINTER                :: distribution_1d
     175              :       TYPE(distribution_2d_type), POINTER                :: distribution_2d
     176              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     177              :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     178              :       LOGICAL                                            :: molecule_only
     179              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     180              : 
     181              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'atom2d_build'
     182              : 
     183              :       INTEGER                                            :: atom_a, handle, ia, iat, iatom, &
     184              :                                                             iatom_local, ikind, imol, natom, &
     185              :                                                             natom_a, natom_local_a, natom_local_b, &
     186              :                                                             nel, nkind
     187        41123 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom2mol, atom_of_kind, listindex, &
     188        41123 :                                                             listsort
     189        41123 :       INTEGER, DIMENSION(:), POINTER                     :: local_cols_array, local_rows_array
     190              : 
     191        41123 :       CALL timeset(routineN, handle)
     192              : 
     193        41123 :       nkind = SIZE(atomic_kind_set)
     194        41123 :       natom = SIZE(particle_set)
     195        41123 :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
     196              : 
     197        41123 :       IF (molecule_only) THEN
     198          876 :          ALLOCATE (atom2mol(natom))
     199          974 :          DO imol = 1, SIZE(molecule_set)
     200         2512 :             DO iat = molecule_set(imol)%first_atom, molecule_set(imol)%last_atom
     201         2220 :                atom2mol(iat) = imol
     202              :             END DO
     203              :          END DO
     204              :       END IF
     205              : 
     206       118267 :       DO ikind = 1, nkind
     207        77144 :          NULLIFY (atom2d(ikind)%list)
     208        77144 :          NULLIFY (atom2d(ikind)%list_local_a_index)
     209        77144 :          NULLIFY (atom2d(ikind)%list_local_b_index)
     210        77144 :          NULLIFY (atom2d(ikind)%list_1d)
     211        77144 :          NULLIFY (atom2d(ikind)%list_a_mol)
     212        77144 :          NULLIFY (atom2d(ikind)%list_b_mol)
     213              : 
     214        77144 :          CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
     215              : 
     216        77144 :          natom_a = SIZE(atom2d(ikind)%list)
     217              : 
     218        77144 :          natom_local_a = distribution_2d%n_local_rows(ikind)
     219        77144 :          natom_local_b = distribution_2d%n_local_cols(ikind)
     220        77144 :          local_rows_array => distribution_2d%local_rows(ikind)%array
     221        77144 :          local_cols_array => distribution_2d%local_cols(ikind)%array
     222              : 
     223        77144 :          nel = distribution_1d%n_el(ikind)
     224       208927 :          ALLOCATE (atom2d(ikind)%list_1d(nel))
     225       175070 :          DO iat = 1, nel
     226        97926 :             ia = distribution_1d%list(ikind)%array(iat)
     227       175070 :             atom2d(ikind)%list_1d(iat) = atom_of_kind(ia)
     228              :          END DO
     229              : 
     230       308576 :          ALLOCATE (listsort(natom_a), listindex(natom_a))
     231       269006 :          listsort(1:natom_a) = atom2d(ikind)%list(1:natom_a)
     232        77144 :          CALL sort(listsort, natom_a, listindex)
     233              :          ! Block rows
     234        77144 :          IF (natom_local_a > 0) THEN
     235       166671 :             ALLOCATE (atom2d(ikind)%list_local_a_index(natom_local_a))
     236       111114 :             ALLOCATE (atom2d(ikind)%list_a_mol(natom_local_a))
     237       162749 :             atom2d(ikind)%list_a_mol(:) = 0
     238              : 
     239              :             ! Build index vector for mapping
     240       162749 :             DO iatom_local = 1, natom_local_a
     241       107192 :                atom_a = local_rows_array(iatom_local)
     242       107192 :                iatom = locate(listsort, atom_a)
     243       107192 :                atom2d(ikind)%list_local_a_index(iatom_local) = listindex(iatom)
     244       162749 :                IF (molecule_only) atom2d(ikind)%list_a_mol(iatom_local) = atom2mol(atom_a)
     245              :             END DO
     246              : 
     247              :          END IF
     248              : 
     249              :          ! Block columns
     250        77144 :          IF (natom_local_b > 0) THEN
     251              : 
     252       231258 :             ALLOCATE (atom2d(ikind)%list_local_b_index(natom_local_b))
     253       154172 :             ALLOCATE (atom2d(ikind)%list_b_mol(natom_local_b))
     254       268800 :             atom2d(ikind)%list_b_mol(:) = 0
     255              : 
     256              :             ! Build index vector for mapping
     257       268800 :             DO iatom_local = 1, natom_local_b
     258       191714 :                atom_a = local_cols_array(iatom_local)
     259       191714 :                iatom = locate(listsort, atom_a)
     260       191714 :                atom2d(ikind)%list_local_b_index(iatom_local) = listindex(iatom)
     261       268800 :                IF (molecule_only) atom2d(ikind)%list_b_mol(iatom_local) = atom2mol(atom_a)
     262              :             END DO
     263              : 
     264              :          END IF
     265              : 
     266       118267 :          DEALLOCATE (listsort, listindex)
     267              : 
     268              :       END DO
     269              : 
     270        41123 :       CALL timestop(handle)
     271              : 
     272        82246 :    END SUBROUTINE atom2d_build
     273              : 
     274              : ! **************************************************************************************************
     275              : !> \brief   Build all the required neighbor lists for Quickstep.
     276              : !> \param qs_env ...
     277              : !> \param para_env ...
     278              : !> \param molecular ...
     279              : !> \param force_env_section ...
     280              : !> \date    28.08.2000
     281              : !> \par History
     282              : !>          - Major refactoring (25.07.2010,jhu)
     283              : !> \author  MK
     284              : !> \version 1.0
     285              : ! **************************************************************************************************
     286        25401 :    SUBROUTINE build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_section)
     287              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     288              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     289              :       LOGICAL, OPTIONAL                                  :: molecular
     290              :       TYPE(section_vals_type), POINTER                   :: force_env_section
     291              : 
     292              :       CHARACTER(len=*), PARAMETER :: routineN = 'build_qs_neighbor_lists'
     293              : 
     294              :       CHARACTER(LEN=2)                                   :: element_symbol, element_symbol2
     295              :       CHARACTER(LEN=default_string_length)               :: print_key_path
     296              :       INTEGER                                            :: handle, hfx_pot, ikind, ingp, iw, jkind, &
     297              :                                                             maxatom, ngp, nkind, zat
     298              :       LOGICAL :: all_potential_present, almo, dftb, do_hfx, dokp, gth_potential_present, &
     299              :          lri_optbas, lrigpw, mic, molecule_only, nddo, paw_atom, paw_atom_present, rigpw, &
     300              :          sgp_potential_present, xtb
     301        25401 :       LOGICAL, ALLOCATABLE, DIMENSION(:) :: all_present, aux_fit_present, aux_present, &
     302        25401 :          core_present, default_present, nonbond1_atom, nonbond2_atom, oce_present, orb_present, &
     303        25401 :          ppl_present, ppnl_present, ri_present, xb1_atom, xb2_atom
     304              :       REAL(dp)                                           :: almo_rcov, almo_rvdw, eps_schwarz, &
     305              :                                                             omega, pdist, rcut, roperator, subcells
     306        25401 :       REAL(dp), ALLOCATABLE, DIMENSION(:) :: all_pot_rad, aux_fit_radius, c_radius, calpha, &
     307        25401 :          core_radius, oce_radius, orb_radius, ppl_radius, ppnl_radius, ri_radius, zeff
     308        25401 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: pair_radius, pair_radius_lb
     309              :       TYPE(all_potential_type), POINTER                  :: all_potential
     310        25401 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     311              :       TYPE(cell_type), POINTER                           :: cell
     312              :       TYPE(cp_logger_type), POINTER                      :: logger
     313              :       TYPE(dft_control_type), POINTER                    :: dft_control
     314              :       TYPE(distribution_1d_type), POINTER                :: distribution_1d
     315              :       TYPE(distribution_2d_type), POINTER                :: distribution_2d
     316              :       TYPE(ewald_environment_type), POINTER              :: ewald_env
     317              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
     318              :       TYPE(gto_basis_set_type), POINTER                  :: aux_basis_set, aux_fit_basis_set, &
     319              :                                                             orb_basis_set, ri_basis_set
     320              :       TYPE(kpoint_type), POINTER                         :: kpoints
     321        25401 :       TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:)  :: atom2d
     322        25401 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     323        25401 :       TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: saa_list, sab_all, sab_almo, &
     324        25401 :          sab_cn, sab_core, sab_gcp, sab_kp, sab_kp_nosym, sab_lrc, sab_orb, sab_scp, sab_se, &
     325        25401 :          sab_tbe, sab_vdw, sab_xb, sab_xtb_nonbond, sab_xtb_pp, sab_xtbe, sac_ae, sac_lri, &
     326        25401 :          sac_ppl, sap_oce, sap_ppnl, soa_list, soo_list
     327        25401 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     328              :       TYPE(paw_proj_set_type), POINTER                   :: paw_proj
     329              :       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_atom
     330              :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
     331              :       TYPE(qs_gcp_type), POINTER                         :: gcp_env
     332        25401 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     333              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     334              :       TYPE(section_vals_type), POINTER                   :: hfx_sections, neighbor_list_section
     335              :       TYPE(sgp_potential_type), POINTER                  :: sgp_potential
     336              :       TYPE(xtb_atom_type), POINTER                       :: xtb_atom
     337              : 
     338        25401 :       CALL timeset(routineN, handle)
     339        25401 :       NULLIFY (logger)
     340        25401 :       logger => cp_get_default_logger()
     341              : 
     342        25401 :       NULLIFY (atomic_kind_set, qs_kind_set, cell, neighbor_list_section, &
     343        25401 :                distribution_1d, distribution_2d, gth_potential, sgp_potential, orb_basis_set, &
     344        25401 :                particle_set, molecule_set, dft_control, ks_env)
     345              : 
     346        25401 :       NULLIFY (sab_orb)
     347        25401 :       NULLIFY (sac_ae)
     348        25401 :       NULLIFY (sac_ppl)
     349        25401 :       NULLIFY (sac_lri)
     350        25401 :       NULLIFY (sap_ppnl)
     351        25401 :       NULLIFY (sap_oce)
     352        25401 :       NULLIFY (sab_se)
     353        25401 :       NULLIFY (sab_lrc)
     354        25401 :       NULLIFY (sab_tbe)
     355        25401 :       NULLIFY (sab_xtbe)
     356        25401 :       NULLIFY (sab_core)
     357        25401 :       NULLIFY (sab_xb)
     358        25401 :       NULLIFY (sab_xtb_pp)
     359        25401 :       NULLIFY (sab_xtb_nonbond)
     360        25401 :       NULLIFY (sab_all)
     361        25401 :       NULLIFY (sab_vdw)
     362        25401 :       NULLIFY (sab_cn)
     363        25401 :       NULLIFY (soo_list)
     364        25401 :       NULLIFY (sab_scp)
     365        25401 :       NULLIFY (sab_almo)
     366        25401 :       NULLIFY (sab_kp)
     367        25401 :       NULLIFY (sab_kp_nosym)
     368              : 
     369              :       CALL get_qs_env(qs_env, &
     370              :                       ks_env=ks_env, &
     371              :                       atomic_kind_set=atomic_kind_set, &
     372              :                       qs_kind_set=qs_kind_set, &
     373              :                       cell=cell, &
     374              :                       kpoints=kpoints, &
     375              :                       distribution_2d=distribution_2d, &
     376              :                       local_particles=distribution_1d, &
     377              :                       particle_set=particle_set, &
     378              :                       molecule_set=molecule_set, &
     379        25401 :                       dft_control=dft_control)
     380              : 
     381        25401 :       neighbor_list_section => section_vals_get_subs_vals(force_env_section, "DFT%PRINT%NEIGHBOR_LISTS")
     382              : 
     383              :       ! This sets the id number of the qs neighbor lists, new lists, means new version
     384              :       ! new version implies new sparsity of the matrices
     385        25401 :       last_qs_neighbor_list_id_nr = last_qs_neighbor_list_id_nr + 1
     386        25401 :       CALL set_ks_env(ks_env=ks_env, neighbor_list_id=last_qs_neighbor_list_id_nr)
     387              : 
     388              :       CALL get_ks_env(ks_env=ks_env, &
     389              :                       sab_orb=sab_orb, &
     390              :                       sac_ae=sac_ae, &
     391              :                       sac_ppl=sac_ppl, &
     392              :                       sac_lri=sac_lri, &
     393              :                       sab_vdw=sab_vdw, &
     394              :                       sap_ppnl=sap_ppnl, &
     395              :                       sap_oce=sap_oce, &
     396              :                       sab_se=sab_se, &
     397              :                       sab_lrc=sab_lrc, &
     398              :                       sab_tbe=sab_tbe, &
     399              :                       sab_xtbe=sab_xtbe, &
     400              :                       sab_core=sab_core, &
     401              :                       sab_xb=sab_xb, &
     402              :                       sab_xtb_pp=sab_xtb_pp, &
     403              :                       sab_xtb_nonbond=sab_xtb_nonbond, &
     404              :                       sab_scp=sab_scp, &
     405              :                       sab_all=sab_all, &
     406              :                       sab_almo=sab_almo, &
     407              :                       sab_kp=sab_kp, &
     408        25401 :                       sab_kp_nosym=sab_kp_nosym)
     409              : 
     410        25401 :       dokp = (kpoints%nkp > 0)
     411        25401 :       nddo = dft_control%qs_control%semi_empirical
     412        25401 :       dftb = dft_control%qs_control%dftb
     413        25401 :       xtb = dft_control%qs_control%xtb
     414        25401 :       almo = dft_control%qs_control%do_almo_scf
     415        25401 :       lrigpw = (dft_control%qs_control%method_id == do_method_lrigpw)
     416        25401 :       rigpw = (dft_control%qs_control%method_id == do_method_rigpw)
     417        25401 :       lri_optbas = dft_control%qs_control%lri_optbas
     418              : 
     419              :       ! molecular lists
     420        25401 :       molecule_only = .FALSE.
     421        25401 :       IF (PRESENT(molecular)) molecule_only = molecular
     422              :       ! minimum image convention (MIC)
     423        25401 :       mic = molecule_only
     424        25401 :       IF (dokp) THEN
     425              :          ! no MIC for kpoints
     426          926 :          mic = .FALSE.
     427        24475 :       ELSEIF (nddo) THEN
     428              :          ! enforce MIC for interaction lists in SE
     429         5782 :          mic = .TRUE.
     430              :       END IF
     431        25401 :       pdist = dft_control%qs_control%pairlist_radius
     432              : 
     433        25401 :       hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
     434        25401 :       CALL section_vals_get(hfx_sections, explicit=do_hfx)
     435              : 
     436        25401 :       CALL get_atomic_kind_set(atomic_kind_set, maxatom=maxatom)
     437              :       CALL get_qs_kind_set(qs_kind_set, paw_atom_present=paw_atom_present, &
     438              :                            gth_potential_present=gth_potential_present, &
     439              :                            sgp_potential_present=sgp_potential_present, &
     440        25401 :                            all_potential_present=all_potential_present)
     441              : 
     442        25401 :       CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
     443              : 
     444              :       ! Allocate work storage
     445        25401 :       nkind = SIZE(atomic_kind_set)
     446              :       ALLOCATE (orb_present(nkind), aux_fit_present(nkind), aux_present(nkind), &
     447       177807 :                 default_present(nkind), core_present(nkind))
     448              :       ALLOCATE (orb_radius(nkind), aux_fit_radius(nkind), c_radius(nkind), &
     449       203208 :                 core_radius(nkind), calpha(nkind), zeff(nkind))
     450        78143 :       orb_radius(:) = 0.0_dp
     451        78143 :       aux_fit_radius(:) = 0.0_dp
     452        78143 :       c_radius(:) = 0.0_dp
     453        78143 :       core_radius(:) = 0.0_dp
     454        78143 :       calpha(:) = 0.0_dp
     455        78143 :       zeff(:) = 0.0_dp
     456              : 
     457       101604 :       ALLOCATE (pair_radius(nkind, nkind))
     458        25401 :       IF (gth_potential_present .OR. sgp_potential_present) THEN
     459        29739 :          ALLOCATE (ppl_present(nkind), ppl_radius(nkind))
     460        27681 :          ppl_radius = 0.0_dp
     461        29739 :          ALLOCATE (ppnl_present(nkind), ppnl_radius(nkind))
     462        43169 :          ppnl_radius = 0.0_dp
     463              :       END IF
     464        25401 :       IF (paw_atom_present) THEN
     465         4776 :          ALLOCATE (oce_present(nkind), oce_radius(nkind))
     466         4716 :          oce_radius = 0.0_dp
     467              :       END IF
     468        25401 :       IF (all_potential_present .OR. sgp_potential_present) THEN
     469        46650 :          ALLOCATE (all_present(nkind), all_pot_rad(nkind))
     470        60499 :          all_pot_rad = 0.0_dp
     471              :       END IF
     472              : 
     473              :       ! Initialize the local data structures
     474       128945 :       ALLOCATE (atom2d(nkind))
     475              :       CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
     476        25401 :                         molecule_set, molecule_only, particle_set=particle_set)
     477              : 
     478        78143 :       DO ikind = 1, nkind
     479              : 
     480        52742 :          CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
     481              : 
     482        52742 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
     483        52742 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX")
     484        52742 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis_set, basis_type="AUX_FIT")
     485              : 
     486              :          CALL get_qs_kind(qs_kind_set(ikind), &
     487              :                           paw_proj_set=paw_proj, &
     488              :                           paw_atom=paw_atom, &
     489              :                           all_potential=all_potential, &
     490              :                           gth_potential=gth_potential, &
     491        52742 :                           sgp_potential=sgp_potential)
     492              : 
     493        52742 :          IF (dftb) THEN
     494              :             ! Set the interaction radius for the neighbor lists (DFTB case)
     495              :             ! This includes all interactions (orbitals and short range pair potential) except vdW
     496         5838 :             CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom)
     497              :             CALL get_dftb_atom_param(dftb_parameter=dftb_atom, &
     498              :                                      cutoff=orb_radius(ikind), &
     499         5838 :                                      defined=orb_present(ikind))
     500              :          ELSE
     501        46904 :             IF (ASSOCIATED(orb_basis_set)) THEN
     502        46902 :                orb_present(ikind) = .TRUE.
     503        46902 :                CALL get_gto_basis_set(gto_basis_set=orb_basis_set, kind_radius=orb_radius(ikind))
     504              :             ELSE
     505            2 :                orb_present(ikind) = .FALSE.
     506              :             END IF
     507              :          END IF
     508              : 
     509        52742 :          IF (ASSOCIATED(aux_basis_set)) THEN
     510            0 :             aux_present(ikind) = .TRUE.
     511              :          ELSE
     512        52742 :             aux_present(ikind) = .FALSE.
     513              :          END IF
     514              : 
     515        52742 :          IF (ASSOCIATED(aux_fit_basis_set)) THEN
     516         1500 :             aux_fit_present(ikind) = .TRUE.
     517         1500 :             CALL get_gto_basis_set(gto_basis_set=aux_fit_basis_set, kind_radius=aux_fit_radius(ikind))
     518              :          ELSE
     519        51242 :             aux_fit_present(ikind) = .FALSE.
     520              :          END IF
     521              : 
     522              :          ! core overlap
     523              :          CALL get_qs_kind(qs_kind_set(ikind), &
     524              :                           alpha_core_charge=calpha(ikind), &
     525              :                           core_charge_radius=core_radius(ikind), &
     526        52742 :                           zeff=zeff(ikind))
     527        52742 :          IF (zeff(ikind) /= 0._dp .AND. calpha(ikind) /= 0._dp) THEN
     528        52564 :             core_present(ikind) = .TRUE.
     529              :          ELSE
     530          178 :             core_present(ikind) = .FALSE.
     531              :          END IF
     532              : 
     533              :          ! Pseudopotentials
     534        52742 :          IF (gth_potential_present .OR. sgp_potential_present) THEN
     535        17768 :             IF (ASSOCIATED(gth_potential)) THEN
     536              :                CALL get_potential(potential=gth_potential, &
     537              :                                   ppl_present=ppl_present(ikind), &
     538              :                                   ppl_radius=ppl_radius(ikind), &
     539              :                                   ppnl_present=ppnl_present(ikind), &
     540        17502 :                                   ppnl_radius=ppnl_radius(ikind))
     541          266 :             ELSE IF (ASSOCIATED(sgp_potential)) THEN
     542              :                CALL get_potential(potential=sgp_potential, &
     543              :                                   ppl_present=ppl_present(ikind), &
     544              :                                   ppl_radius=ppl_radius(ikind), &
     545              :                                   ppnl_present=ppnl_present(ikind), &
     546           58 :                                   ppnl_radius=ppnl_radius(ikind))
     547              :             ELSE
     548          208 :                ppl_present(ikind) = .FALSE.
     549          208 :                ppnl_present(ikind) = .FALSE.
     550              :             END IF
     551              :          END IF
     552              : 
     553              :          ! GAPW
     554        52742 :          IF (paw_atom_present) THEN
     555         3124 :             IF (paw_atom) THEN
     556         2956 :                oce_present(ikind) = .TRUE.
     557         2956 :                CALL get_paw_proj_set(paw_proj_set=paw_proj, rcprj=oce_radius(ikind))
     558              :             ELSE
     559          168 :                oce_present(ikind) = .FALSE.
     560              :             END IF
     561              :          END IF
     562              : 
     563              :          ! Check the presence of an all electron potential or ERFC potential
     564       130885 :          IF (all_potential_present .OR. sgp_potential_present) THEN
     565        35098 :             all_present(ikind) = .FALSE.
     566        35098 :             all_pot_rad(ikind) = 0.0_dp
     567        35098 :             IF (ASSOCIATED(all_potential)) THEN
     568        35004 :                all_present(ikind) = .TRUE.
     569        35004 :                CALL get_potential(potential=all_potential, core_charge_radius=all_pot_rad(ikind))
     570           94 :             ELSE IF (ASSOCIATED(sgp_potential)) THEN
     571           58 :                IF (sgp_potential%ecp_local) THEN
     572           46 :                   all_present(ikind) = .TRUE.
     573           46 :                   CALL get_potential(potential=sgp_potential, core_charge_radius=all_pot_rad(ikind))
     574              :                END IF
     575              :             END IF
     576              :          END IF
     577              : 
     578              :       END DO
     579              : 
     580              :       ! Build the orbital-orbital overlap neighbor lists
     581        25401 :       IF (pdist < 0.0_dp) THEN
     582              :          pdist = MAX(plane_distance(1, 0, 0, cell), &
     583              :                      plane_distance(0, 1, 0, cell), &
     584            4 :                      plane_distance(0, 0, 1, cell))
     585              :       END IF
     586        25401 :       CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius, pdist)
     587              :       CALL build_neighbor_lists(sab_orb, particle_set, atom2d, cell, pair_radius, &
     588        25401 :                                 mic=mic, subcells=subcells, molecular=molecule_only, nlname="sab_orb")
     589        25401 :       CALL set_ks_env(ks_env=ks_env, sab_orb=sab_orb)
     590              :       CALL write_neighbor_lists(sab_orb, particle_set, cell, para_env, neighbor_list_section, &
     591        25401 :                                 "/SAB_ORB", "sab_orb", "ORBITAL ORBITAL")
     592              : 
     593              :       ! Build orbital-orbital list containing all the pairs, to be used with
     594              :       ! non-symmetric operators. Beware: the cutoff of the orbital-orbital overlap
     595              :       ! might not be optimal. It should be verified for each operator.
     596        25401 :       IF (.NOT. (nddo .OR. dftb .OR. xtb)) THEN
     597              :          CALL build_neighbor_lists(sab_all, particle_set, atom2d, cell, pair_radius, &
     598        10609 :                                    mic=mic, symmetric=.FALSE., subcells=subcells, molecular=molecule_only, nlname="sab_all")
     599        10609 :          CALL set_ks_env(ks_env=ks_env, sab_all=sab_all)
     600              :       END IF
     601              : 
     602              :       ! Build the core-core overlap neighbor lists
     603        25401 :       IF (.NOT. (nddo .OR. dftb .OR. xtb)) THEN
     604        10609 :          CALL pair_radius_setup(core_present, core_present, core_radius, core_radius, pair_radius)
     605              :          CALL build_neighbor_lists(sab_core, particle_set, atom2d, cell, pair_radius, subcells=subcells, &
     606        10609 :                                    operator_type="PP", nlname="sab_core")
     607        10609 :          CALL set_ks_env(ks_env=ks_env, sab_core=sab_core)
     608              :          CALL write_neighbor_lists(sab_core, particle_set, cell, para_env, neighbor_list_section, &
     609        10609 :                                    "/SAB_CORE", "sab_core", "CORE CORE")
     610              :       END IF
     611              : 
     612        25401 :       IF (dokp) THEN
     613              :          ! We try to guess an integration radius for K-points
     614              :          ! For non-HFX calculations we use the overlap list
     615              :          ! For HFX we use the interaction radius of kinds (ORB or ADMM basis)
     616              :          ! plus a range for the operator
     617          926 :          IF (do_hfx) THEN
     618              : 
     619              :             !case study on the HFX potential: TC, SR or Overlap?
     620           80 :             CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%POTENTIAL_TYPE", i_val=hfx_pot)
     621              : 
     622           26 :             SELECT CASE (hfx_pot)
     623              :             CASE (do_potential_id)
     624           26 :                roperator = 0.0_dp
     625              :             CASE (do_potential_truncated)
     626           54 :                CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", r_val=roperator)
     627              :             CASE (do_potential_mix_cl_trunc)
     628            8 :                CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", r_val=roperator)
     629              :             CASE (do_potential_short)
     630            0 :                CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%OMEGA", r_val=omega)
     631            0 :                CALL section_vals_val_get(hfx_sections, "SCREENING%EPS_SCHWARZ", r_val=eps_schwarz)
     632            0 :                CALL erfc_cutoff(eps_schwarz, omega, roperator)
     633              :             CASE DEFAULT
     634           80 :                CPABORT("HFX potential not available for K-points (NYI)")
     635              :             END SELECT
     636              : 
     637           80 :             IF (dft_control%do_admm) THEN
     638              :                CALL pair_radius_setup(aux_fit_present, aux_fit_present, aux_fit_radius, aux_fit_radius, &
     639           44 :                                       pair_radius)
     640              : 
     641              :                !We cannot accept a pair radius smaller than the ORB overlap, for sanity reasons
     642          132 :                ALLOCATE (pair_radius_lb(nkind, nkind))
     643           44 :                CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius_lb)
     644          110 :                DO jkind = 1, nkind
     645          220 :                   DO ikind = 1, nkind
     646          110 :                      IF (pair_radius(ikind, jkind) + cutoff_screen_factor*roperator .LE. pair_radius_lb(ikind, jkind)) &
     647          134 :                         pair_radius(ikind, jkind) = pair_radius_lb(ikind, jkind) - roperator
     648              :                   END DO
     649              :                END DO
     650              :             ELSE
     651           36 :                CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
     652              :             END IF
     653          400 :             pair_radius = pair_radius + cutoff_screen_factor*roperator
     654              :          ELSE
     655          846 :             CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
     656              :          END IF
     657              :          CALL build_neighbor_lists(sab_kp, particle_set, atom2d, cell, pair_radius, &
     658          926 :                                    subcells=subcells, nlname="sab_kp")
     659          926 :          CALL set_ks_env(ks_env=ks_env, sab_kp=sab_kp)
     660              : 
     661          926 :          IF (do_hfx) THEN
     662              :             CALL build_neighbor_lists(sab_kp_nosym, particle_set, atom2d, cell, pair_radius, &
     663           80 :                                       subcells=subcells, nlname="sab_kp_nosym", symmetric=.FALSE.)
     664           80 :             CALL set_ks_env(ks_env=ks_env, sab_kp_nosym=sab_kp_nosym)
     665              :          END IF
     666              :       END IF
     667              : 
     668              :       ! Build orbital GTH-PPL operator overlap list
     669        25401 :       IF (gth_potential_present .OR. sgp_potential_present) THEN
     670        10009 :          IF (ANY(ppl_present)) THEN
     671         9911 :             CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius)
     672              :             CALL build_neighbor_lists(sac_ppl, particle_set, atom2d, cell, pair_radius, &
     673         9911 :                                       subcells=subcells, operator_type="ABC", nlname="sac_ppl")
     674         9911 :             CALL set_ks_env(ks_env=ks_env, sac_ppl=sac_ppl)
     675              :             CALL write_neighbor_lists(sac_ppl, particle_set, cell, para_env, neighbor_list_section, &
     676         9911 :                                       "/SAC_PPL", "sac_ppl", "ORBITAL GTH-PPL")
     677         9911 :             IF (lrigpw) THEN
     678           58 :                IF (qs_env%lri_env%ppl_ri) THEN
     679              :                   CALL build_neighbor_lists(sac_lri, particle_set, atom2d, cell, pair_radius, &
     680            2 :                                             subcells=subcells, symmetric=.FALSE., operator_type="PP", nlname="sac_lri")
     681            2 :                   CALL set_ks_env(ks_env=ks_env, sac_lri=sac_lri)
     682              :                END IF
     683              :             END IF
     684              :          END IF
     685              : 
     686        13119 :          IF (ANY(ppnl_present)) THEN
     687         7989 :             CALL pair_radius_setup(orb_present, ppnl_present, orb_radius, ppnl_radius, pair_radius)
     688              :             CALL build_neighbor_lists(sap_ppnl, particle_set, atom2d, cell, pair_radius, &
     689         7989 :                                       subcells=subcells, operator_type="ABBA", nlname="sap_ppnl")
     690         7989 :             CALL set_ks_env(ks_env=ks_env, sap_ppnl=sap_ppnl)
     691              :             CALL write_neighbor_lists(sap_ppnl, particle_set, cell, para_env, neighbor_list_section, &
     692         7989 :                                       "/SAP_PPNL", "sap_ppnl", "ORBITAL GTH-PPNL")
     693              :          END IF
     694              :       END IF
     695              : 
     696        25401 :       IF (paw_atom_present) THEN
     697              :          ! Build orbital-GAPW projector overlap list
     698         1660 :          IF (ANY(oce_present)) THEN
     699         1592 :             CALL pair_radius_setup(orb_present, oce_present, orb_radius, oce_radius, pair_radius)
     700              :             CALL build_neighbor_lists(sap_oce, particle_set, atom2d, cell, pair_radius, &
     701         1592 :                                       subcells=subcells, operator_type="ABBA", nlname="sap_oce")
     702         1592 :             CALL set_ks_env(ks_env=ks_env, sap_oce=sap_oce)
     703              :             CALL write_neighbor_lists(sap_oce, particle_set, cell, para_env, neighbor_list_section, &
     704         1592 :                                       "/SAP_OCE", "sap_oce", "ORBITAL(A) PAW-PRJ")
     705              :          END IF
     706              :       END IF
     707              : 
     708              :       ! Build orbital-ERFC potential list
     709        25401 :       IF (.NOT. (nddo .OR. dftb .OR. xtb)) THEN
     710        10609 :          IF (all_potential_present .OR. sgp_potential_present) THEN
     711          758 :             CALL pair_radius_setup(orb_present, all_present, orb_radius, all_pot_rad, pair_radius)
     712              :             CALL build_neighbor_lists(sac_ae, particle_set, atom2d, cell, pair_radius, &
     713          758 :                                       subcells=subcells, operator_type="ABC", nlname="sac_ae")
     714          758 :             CALL set_ks_env(ks_env=ks_env, sac_ae=sac_ae)
     715              :             CALL write_neighbor_lists(sac_ae, particle_set, cell, para_env, neighbor_list_section, &
     716          758 :                                       "/SAC_AE", "sac_ae", "ORBITAL ERFC POTENTIAL")
     717              :          END IF
     718              :       END IF
     719              : 
     720        25401 :       IF (nddo) THEN
     721              :          ! Semi-empirical neighbor lists
     722        18448 :          default_present = .TRUE.
     723        18448 :          c_radius = dft_control%qs_control%se_control%cutoff_cou
     724              :          ! Build the neighbor lists for the Hartree terms
     725         5782 :          CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
     726         5782 :          IF (dft_control%qs_control%se_control%do_ewald_gks) THEN
     727              :             ! Use MIC for the periodic code of GKS
     728              :             CALL build_neighbor_lists(sab_se, particle_set, atom2d, cell, pair_radius, mic=mic, &
     729            2 :                                       subcells=subcells, nlname="sab_se")
     730              :          ELSE
     731              :             CALL build_neighbor_lists(sab_se, particle_set, atom2d, cell, pair_radius, &
     732         5780 :                                       subcells=subcells, nlname="sab_se")
     733              :          END IF
     734         5782 :          CALL set_ks_env(ks_env=ks_env, sab_se=sab_se)
     735              :          CALL write_neighbor_lists(sab_se, particle_set, cell, para_env, neighbor_list_section, &
     736         5782 :                                    "/SAB_SE", "sab_se", "HARTREE INTERACTIONS")
     737              : 
     738              :          ! If requested build the SE long-range correction neighbor list
     739         5782 :          IF ((dft_control%qs_control%se_control%do_ewald) .AND. &
     740              :              (dft_control%qs_control%se_control%integral_screening /= do_se_IS_slater)) THEN
     741          328 :             c_radius = dft_control%qs_control%se_control%cutoff_lrc
     742          140 :             CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
     743              :             CALL build_neighbor_lists(sab_lrc, particle_set, atom2d, cell, pair_radius, &
     744          140 :                                       subcells=subcells, nlname="sab_lrc")
     745          140 :             CALL set_ks_env(ks_env=ks_env, sab_lrc=sab_lrc)
     746              :             CALL write_neighbor_lists(sab_lrc, particle_set, cell, para_env, neighbor_list_section, &
     747          140 :                                       "/SAB_LRC", "sab_lrc", "SE LONG-RANGE CORRECTION")
     748              :          END IF
     749              :       END IF
     750              : 
     751        25401 :       IF (dftb) THEN
     752              :          ! Build the neighbor lists for the DFTB Ewald methods
     753         2884 :          IF (dft_control%qs_control%dftb_control%do_ewald) THEN
     754         1074 :             CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env)
     755         1074 :             CALL ewald_env_get(ewald_env, rcut=rcut)
     756         3150 :             c_radius = rcut
     757         1074 :             CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius)
     758              :             CALL build_neighbor_lists(sab_tbe, particle_set, atom2d, cell, pair_radius, mic=mic, &
     759         1074 :                                       subcells=subcells, nlname="sab_tbe")
     760         1074 :             CALL set_ks_env(ks_env=ks_env, sab_tbe=sab_tbe)
     761              :          END IF
     762              : 
     763              :          ! Build the neighbor lists for the DFTB vdW pair potential
     764         2884 :          IF (dft_control%qs_control%dftb_control%dispersion) THEN
     765         1016 :             IF (dft_control%qs_control%dftb_control%dispersion_type == dispersion_uff) THEN
     766         2754 :                DO ikind = 1, nkind
     767         1828 :                   CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom)
     768         2754 :                   CALL get_dftb_atom_param(dftb_parameter=dftb_atom, rcdisp=c_radius(ikind))
     769              :                END DO
     770         2754 :                default_present = .TRUE.
     771          926 :                CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
     772              :                CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, &
     773          926 :                                          subcells=subcells, nlname="sab_vdw")
     774          926 :                CALL set_ks_env(ks_env=ks_env, sab_vdw=sab_vdw)
     775              :             END IF
     776              :          END IF
     777              :       END IF
     778              : 
     779        25401 :       IF (xtb .AND. (.NOT. dft_control%qs_control%xtb_control%do_tblite)) THEN
     780              :          ! Build the neighbor lists for the xTB Ewald method
     781         5470 :          IF (dft_control%qs_control%xtb_control%do_ewald) THEN
     782         2008 :             CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env)
     783         2008 :             CALL ewald_env_get(ewald_env, rcut=rcut)
     784         6978 :             c_radius = rcut
     785         2008 :             CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius)
     786              :             CALL build_neighbor_lists(sab_tbe, particle_set, atom2d, cell, pair_radius, mic=mic, &
     787         2008 :                                       subcells=subcells, nlname="sab_tbe")
     788         2008 :             CALL set_ks_env(ks_env=ks_env, sab_tbe=sab_tbe)
     789              :          END IF
     790              :          ! Repulsive Potential
     791        53850 :          pair_radius(1:nkind, 1:nkind) = dft_control%qs_control%xtb_control%rcpair(1:nkind, 1:nkind)
     792        18768 :          default_present = .TRUE.
     793              :          CALL build_neighbor_lists(sab_xtb_pp, particle_set, atom2d, cell, pair_radius, &
     794         5470 :                                    subcells=subcells, nlname="sab_xtb_pp")
     795         5470 :          CALL set_ks_env(ks_env=ks_env, sab_xtb_pp=sab_xtb_pp)
     796              :          ! SR part of Coulomb interaction
     797        18768 :          DO ikind = 1, nkind
     798        13298 :             CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom)
     799        18768 :             CALL get_xtb_atom_param(xtb_parameter=xtb_atom, rcut=c_radius(ikind))
     800              :          END DO
     801        18768 :          default_present = .TRUE.
     802         5470 :          CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
     803              :          CALL build_neighbor_lists(sab_xtbe, particle_set, atom2d, cell, pair_radius, &
     804         5470 :                                    subcells=subcells, nlname="sab_xtbe")
     805         5470 :          CALL set_ks_env(ks_env=ks_env, sab_xtbe=sab_xtbe)
     806              :          ! XB list
     807        16410 :          ALLOCATE (xb1_atom(nkind), xb2_atom(nkind))
     808        18768 :          c_radius = 0.5_dp*dft_control%qs_control%xtb_control%xb_radius
     809        18768 :          DO ikind = 1, nkind
     810        13298 :             CALL get_atomic_kind(atomic_kind_set(ikind), z=zat)
     811        13298 :             IF (zat == 17 .OR. zat == 35 .OR. zat == 53 .OR. zat == 85) THEN
     812          130 :                xb1_atom(ikind) = .TRUE.
     813              :             ELSE
     814        13168 :                xb1_atom(ikind) = .FALSE.
     815              :             END IF
     816        32066 :             IF (zat == 7 .OR. zat == 8 .OR. zat == 15 .OR. zat == 16) THEN
     817         4960 :                xb2_atom(ikind) = .TRUE.
     818              :             ELSE
     819         8338 :                xb2_atom(ikind) = .FALSE.
     820              :             END IF
     821              :          END DO
     822         5470 :          CALL pair_radius_setup(xb1_atom, xb2_atom, c_radius, c_radius, pair_radius)
     823              :          CALL build_neighbor_lists(sab_xb, particle_set, atom2d, cell, pair_radius, &
     824         5470 :                                    symmetric=.FALSE., subcells=subcells, operator_type="PP", nlname="sab_xb")
     825         5470 :          CALL set_ks_env(ks_env=ks_env, sab_xb=sab_xb)
     826              :          CALL write_neighbor_lists(sab_xb, particle_set, cell, para_env, neighbor_list_section, &
     827         5470 :                                    "/SAB_XB", "sab_xb", "XB bonding")
     828              : 
     829              :          ! nonbonded interactions list
     830              :          IF (dft_control%qs_control%xtb_control%do_nonbonded &
     831         5470 :              .AND. (.NOT. dft_control%qs_control%xtb_control%do_tblite)) THEN
     832           24 :             ngp = SIZE(dft_control%qs_control%xtb_control%nonbonded%pot)
     833           72 :             ALLOCATE (nonbond1_atom(nkind), nonbond2_atom(nkind))
     834          120 :             nonbond1_atom = .FALSE.
     835          120 :             nonbond2_atom = .FALSE.
     836           48 :             DO ingp = 1, ngp
     837          120 :                DO ikind = 1, nkind
     838           96 :                   rcut = SQRT(dft_control%qs_control%xtb_control%nonbonded%pot(ingp)%pot%rcutsq)
     839          480 :                   c_radius = rcut
     840           96 :                   CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=element_symbol)
     841           96 :                   CALL uppercase(element_symbol)
     842          120 :                   IF (TRIM(dft_control%qs_control%xtb_control%nonbonded%pot(ingp)%pot%at1) == TRIM(element_symbol)) THEN
     843           24 :                      nonbond1_atom(ikind) = .TRUE.
     844          120 :                      DO jkind = 1, nkind
     845           96 :                         CALL get_atomic_kind(atomic_kind_set(jkind), element_symbol=element_symbol2)
     846           96 :                         CALL uppercase(element_symbol2)
     847          120 :                         IF (TRIM(dft_control%qs_control%xtb_control%nonbonded%pot(ingp)%pot%at2) == TRIM(element_symbol2)) THEN
     848           24 :                            nonbond2_atom(jkind) = .TRUE.
     849              :                         END IF
     850              :                      END DO
     851              :                   END IF
     852              :                END DO
     853           24 :                CALL pair_radius_setup(nonbond1_atom, nonbond2_atom, c_radius, c_radius, pair_radius)
     854              :                CALL build_neighbor_lists(sab_xtb_nonbond, particle_set, atom2d, cell, pair_radius, &
     855           24 :                                          symmetric=.FALSE., subcells=subcells, operator_type="PP", nlname="sab_xtb_nonbond")
     856           24 :                CALL set_ks_env(ks_env=ks_env, sab_xtb_nonbond=sab_xtb_nonbond)
     857              :                CALL write_neighbor_lists(sab_xtb_nonbond, particle_set, cell, para_env, neighbor_list_section, &
     858           48 :                                          "/SAB_XTB_NONBOND", "sab_xtb_nonbond", "XTB NONBONDED INTERACTIONS")
     859              :             END DO
     860              :          END IF
     861              :       END IF
     862              : 
     863              :       ! Build the neighbor lists for the vdW pair potential
     864        25401 :       IF (.NOT. dft_control%qs_control%xtb_control%do_tblite) THEN
     865        24745 :          CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
     866        24745 :          sab_vdw => dispersion_env%sab_vdw
     867        24745 :          sab_cn => dispersion_env%sab_cn
     868        24745 :          IF (dispersion_env%type == xc_vdw_fun_pairpot .OR. xtb) THEN
     869         5836 :             IF (dispersion_env%pp_type == vdw_pairpot_dftd4) THEN
     870         2928 :                c_radius(:) = dispersion_env%rc_d4
     871              :             ELSE
     872        16870 :                c_radius(:) = dispersion_env%rc_disp
     873              :             END IF
     874        19798 :             default_present = .TRUE. !include all atoms in vdW (even without basis)
     875         5836 :             CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
     876              :             CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, &
     877         5836 :                                       subcells=subcells, operator_type="PP", nlname="sab_vdw")
     878         5836 :             dispersion_env%sab_vdw => sab_vdw
     879              : 
     880              :             ! Build the neighbor lists for coordination numbers as needed by the DFT-D3/D4 method
     881              :             ! This is also needed for the xTB Hamiltonian
     882        19798 :             DO ikind = 1, nkind
     883        13962 :                CALL get_atomic_kind(atomic_kind_set(ikind), z=zat)
     884        19798 :                c_radius(ikind) = 4._dp*ptable(zat)%covalent_radius*bohr
     885              :             END DO
     886         5836 :             CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
     887              :             CALL build_neighbor_lists(sab_cn, particle_set, atom2d, cell, pair_radius, &
     888         5836 :                                       subcells=subcells, operator_type="PP", nlname="sab_cn")
     889         5836 :             dispersion_env%sab_cn => sab_cn
     890              :          END IF
     891              :       END IF
     892              : 
     893              :       ! Build the neighbor lists for the gCP pair potential
     894        25401 :       NULLIFY (gcp_env)
     895        25401 :       CALL get_qs_env(qs_env=qs_env, gcp_env=gcp_env)
     896        25401 :       IF (ASSOCIATED(gcp_env)) THEN
     897        10609 :          IF (gcp_env%do_gcp) THEN
     898            4 :             sab_gcp => gcp_env%sab_gcp
     899           10 :             DO ikind = 1, nkind
     900           10 :                c_radius(ikind) = gcp_env%gcp_kind(ikind)%rcsto
     901              :             END DO
     902            4 :             CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius)
     903              :             CALL build_neighbor_lists(sab_gcp, particle_set, atom2d, cell, pair_radius, &
     904            4 :                                       subcells=subcells, operator_type="PP", nlname="sab_gcp")
     905            4 :             gcp_env%sab_gcp => sab_gcp
     906              :          ELSE
     907        10605 :             NULLIFY (gcp_env%sab_gcp)
     908              :          END IF
     909              :       END IF
     910              : 
     911        25401 :       IF (lrigpw .OR. lri_optbas) THEN
     912              :          ! set neighborlists in lri_env environment
     913           64 :          CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
     914           64 :          soo_list => qs_env%lri_env%soo_list
     915              :          CALL build_neighbor_lists(soo_list, particle_set, atom2d, cell, pair_radius, &
     916           64 :                                    mic=mic, molecular=molecule_only, subcells=subcells, nlname="soo_list")
     917           64 :          qs_env%lri_env%soo_list => soo_list
     918              :          CALL write_neighbor_lists(soo_list, particle_set, cell, para_env, neighbor_list_section, &
     919           64 :                                    "/SOO_LIST", "soo_list", "ORBITAL ORBITAL (RI)")
     920        25337 :       ELSEIF (rigpw) THEN
     921            0 :          ALLOCATE (ri_present(nkind), ri_radius(nkind))
     922            0 :          ri_present = .FALSE.
     923            0 :          ri_radius = 0.0_dp
     924            0 :          DO ikind = 1, nkind
     925            0 :             CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="RI_HXC")
     926            0 :             IF (ASSOCIATED(ri_basis_set)) THEN
     927            0 :                ri_present(ikind) = .TRUE.
     928            0 :                CALL get_gto_basis_set(gto_basis_set=ri_basis_set, kind_radius=ri_radius(ikind))
     929              :             ELSE
     930            0 :                ri_present(ikind) = .FALSE.
     931              :             END IF
     932              :          END DO
     933              :          ! set neighborlists in lri_env environment
     934            0 :          CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
     935            0 :          soo_list => qs_env%lri_env%soo_list
     936              :          CALL build_neighbor_lists(soo_list, particle_set, atom2d, cell, pair_radius, &
     937            0 :                                    mic=mic, molecular=molecule_only, subcells=subcells, nlname="soo_list")
     938            0 :          qs_env%lri_env%soo_list => soo_list
     939              :          !
     940            0 :          CALL pair_radius_setup(ri_present, ri_present, ri_radius, ri_radius, pair_radius)
     941            0 :          saa_list => qs_env%lri_env%saa_list
     942              :          CALL build_neighbor_lists(saa_list, particle_set, atom2d, cell, pair_radius, &
     943            0 :                                    mic=mic, molecular=molecule_only, subcells=subcells, nlname="saa_list")
     944            0 :          qs_env%lri_env%saa_list => saa_list
     945              :          !
     946            0 :          CALL pair_radius_setup(ri_present, orb_present, ri_radius, orb_radius, pair_radius)
     947            0 :          soa_list => qs_env%lri_env%soa_list
     948              :          CALL build_neighbor_lists(soa_list, particle_set, atom2d, cell, pair_radius, &
     949              :                                    mic=mic, symmetric=.FALSE., molecular=molecule_only, &
     950            0 :                                    subcells=subcells, operator_type="ABC", nlname="saa_list")
     951            0 :          qs_env%lri_env%soa_list => soa_list
     952              :       END IF
     953              : 
     954              :       ! Build the neighbor lists for the ALMO delocalization
     955        25401 :       IF (almo) THEN
     956          360 :          DO ikind = 1, nkind
     957          244 :             CALL get_atomic_kind(atomic_kind_set(ikind), rcov=almo_rcov, rvdw=almo_rvdw)
     958              :             ! multiply the radius by some hard-coded number
     959              :             c_radius(ikind) = MAX(almo_rcov, almo_rvdw)*bohr* &
     960          360 :                               almo_max_cutoff_multiplier
     961              :          END DO
     962          360 :          default_present = .TRUE. !include all atoms (even without basis)
     963          116 :          CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
     964              :          CALL build_neighbor_lists(sab_almo, particle_set, atom2d, cell, pair_radius, &
     965          116 :                                    subcells=subcells, operator_type="PP", nlname="sab_almo")
     966          116 :          CALL set_ks_env(ks_env=ks_env, sab_almo=sab_almo)
     967              :       END IF
     968              : 
     969              :       ! Print particle distribution
     970        25401 :       print_key_path = "PRINT%DISTRIBUTION"
     971        25401 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, force_env_section, &
     972              :                                            print_key_path), &
     973              :                 cp_p_file)) THEN
     974              :          iw = cp_print_key_unit_nr(logger=logger, &
     975              :                                    basis_section=force_env_section, &
     976              :                                    print_key_path=print_key_path, &
     977          158 :                                    extension=".out")
     978          158 :          CALL write_neighbor_distribution(sab_orb, qs_kind_set, iw, para_env)
     979              :          CALL cp_print_key_finished_output(unit_nr=iw, &
     980              :                                            logger=logger, &
     981              :                                            basis_section=force_env_section, &
     982          158 :                                            print_key_path=print_key_path)
     983              :       END IF
     984              : 
     985              :       ! Release work storage
     986        25401 :       CALL atom2d_cleanup(atom2d)
     987              : 
     988        25401 :       DEALLOCATE (atom2d)
     989        25401 :       DEALLOCATE (orb_present, default_present, core_present)
     990        25401 :       DEALLOCATE (orb_radius, aux_fit_radius, c_radius, core_radius)
     991        25401 :       DEALLOCATE (calpha, zeff)
     992        25401 :       DEALLOCATE (pair_radius)
     993        25401 :       IF (gth_potential_present .OR. sgp_potential_present) THEN
     994         9913 :          DEALLOCATE (ppl_present, ppl_radius)
     995         9913 :          DEALLOCATE (ppnl_present, ppnl_radius)
     996              :       END IF
     997        25401 :       IF (paw_atom_present) THEN
     998         1592 :          DEALLOCATE (oce_present, oce_radius)
     999              :       END IF
    1000        25401 :       IF (all_potential_present .OR. sgp_potential_present) THEN
    1001        15550 :          DEALLOCATE (all_present, all_pot_rad)
    1002              :       END IF
    1003              : 
    1004        25401 :       CALL timestop(handle)
    1005              : 
    1006        76203 :    END SUBROUTINE build_qs_neighbor_lists
    1007              : 
    1008              : ! **************************************************************************************************
    1009              : !> \brief   Build simple pair neighbor lists.
    1010              : !> \param ab_list ...
    1011              : !> \param particle_set ...
    1012              : !> \param atom ...
    1013              : !> \param cell ...
    1014              : !> \param pair_radius ...
    1015              : !> \param subcells ...
    1016              : !> \param mic ...
    1017              : !> \param symmetric ...
    1018              : !> \param molecular ...
    1019              : !> \param subset_of_mol ...
    1020              : !> \param current_subset ...
    1021              : !> \param operator_type ...
    1022              : !> \param nlname ...
    1023              : !> \param atomb_to_keep the list of atom indices to keep for pairs from the atom2d%b_list
    1024              : !> \date    20.03.2002
    1025              : !> \par History
    1026              : !>          - Major refactoring (25.07.2010,jhu)
    1027              : !>          - Added option to filter out atoms from list_b (08.2018, A.  Bussy)
    1028              : !> \author  MK
    1029              : !> \version 2.0
    1030              : ! **************************************************************************************************
    1031       124735 :    SUBROUTINE build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, &
    1032              :                                    mic, symmetric, molecular, subset_of_mol, current_subset, &
    1033       124735 :                                    operator_type, nlname, atomb_to_keep)
    1034              : 
    1035              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1036              :          POINTER                                         :: ab_list
    1037              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1038              :       TYPE(local_atoms_type), DIMENSION(:), INTENT(IN)   :: atom
    1039              :       TYPE(cell_type), POINTER                           :: cell
    1040              :       REAL(dp), DIMENSION(:, :), INTENT(IN)              :: pair_radius
    1041              :       REAL(dp), INTENT(IN)                               :: subcells
    1042              :       LOGICAL, INTENT(IN), OPTIONAL                      :: mic, symmetric, molecular
    1043              :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: subset_of_mol
    1044              :       INTEGER, OPTIONAL                                  :: current_subset
    1045              :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: operator_type
    1046              :       CHARACTER(LEN=*), INTENT(IN)                       :: nlname
    1047              :       INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL        :: atomb_to_keep
    1048              : 
    1049              :       CHARACTER(len=*), PARAMETER :: routineN = 'build_neighbor_lists'
    1050              : 
    1051              :       TYPE local_lists
    1052              :          INTEGER, DIMENSION(:), POINTER           :: list
    1053              :       END TYPE local_lists
    1054              : 
    1055              :       INTEGER :: atom_a, atom_b, handle, i, iab, iatom, iatom_local, &
    1056              :                  iatom_subcell, icell, ikind, j, jatom, jatom_local, jcell, jkind, k, &
    1057              :                  kcell, maxat, mol_a, mol_b, nkind, otype, natom, inode, nnode, nentry
    1058              :       INTEGER, DIMENSION(3)                              :: cell_b, ncell, nsubcell, periodic
    1059       124735 :       INTEGER, DIMENSION(:), POINTER                     :: index_list
    1060              :       LOGICAL                                            :: include_ab, my_mic, &
    1061              :                                                             my_molecular, my_symmetric, my_sort_atomb
    1062       124735 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: pres_a, pres_b
    1063              :       REAL(dp)                                           :: rab2, rab2_max, rab_max, rabm, deth, subcell_scale
    1064              :       REAL(dp), DIMENSION(3)                             :: r, rab, ra, rb, sab_max, sb, &
    1065              :                                                             sb_pbc, sb_min, sb_max, rab_pbc, pd, sab_max_guard
    1066       124735 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: nlista, nlistb
    1067       124735 :       TYPE(local_lists), DIMENSION(:), POINTER           :: lista, listb
    1068              :       TYPE(neighbor_list_p_type), &
    1069       124735 :          ALLOCATABLE, DIMENSION(:)                       :: kind_a
    1070              :       TYPE(neighbor_list_set_type), POINTER              :: neighbor_list_set
    1071              :       TYPE(subcell_type), DIMENSION(:, :, :), &
    1072       124735 :          POINTER                                         :: subcell
    1073       124735 :       REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE        :: r_pbc
    1074              :       TYPE(neighbor_list_iterator_p_type), &
    1075       124735 :          DIMENSION(:), POINTER                           :: nl_iterator
    1076              : 
    1077       124735 :       CALL timeset(routineN//"_"//TRIM(nlname), handle)
    1078              : 
    1079              :       ! input options
    1080       124735 :       my_mic = .FALSE.
    1081       124735 :       IF (PRESENT(mic)) my_mic = mic
    1082       124735 :       my_symmetric = .TRUE.
    1083       124735 :       IF (PRESENT(symmetric)) my_symmetric = symmetric
    1084       124735 :       my_molecular = .FALSE.
    1085              :       ! if we have a molecular NL, MIC has to be used
    1086       124735 :       IF (PRESENT(molecular)) my_molecular = molecular
    1087              :       ! check for operator types
    1088       124735 :       IF (PRESENT(operator_type)) THEN
    1089              :          SELECT CASE (operator_type)
    1090              :          CASE ("AB")
    1091        11487 :             otype = 1 ! simple overlap
    1092              :          CASE ("ABC")
    1093        11487 :             otype = 2 ! for three center operators
    1094        11487 :             CPASSERT(.NOT. my_molecular)
    1095        11487 :             my_symmetric = .FALSE.
    1096              :          CASE ("ABBA")
    1097        10357 :             otype = 3 ! for separable nonlocal operators
    1098        10357 :             my_symmetric = .FALSE.
    1099              :          CASE ("PP")
    1100        27957 :             otype = 4 ! simple atomic pair potential list
    1101              :          CASE default
    1102        49801 :             CPABORT("")
    1103              :          END SELECT
    1104              :       ELSE
    1105              :          ! default is a simple AB neighbor list
    1106              :          otype = 1
    1107              :       END IF
    1108       124735 :       my_sort_atomb = .FALSE.
    1109       124735 :       IF (PRESENT(atomb_to_keep)) THEN
    1110          374 :          my_sort_atomb = .TRUE.
    1111              :       END IF
    1112              : 
    1113       124735 :       nkind = SIZE(atom)
    1114              :       ! Deallocate the old neighbor list structure
    1115       124735 :       CALL release_neighbor_list_sets(ab_list)
    1116              :       ! Allocate and initialize the new neighbor list structure
    1117       942152 :       ALLOCATE (ab_list(nkind*nkind))
    1118       692682 :       DO iab = 1, SIZE(ab_list)
    1119       567947 :          NULLIFY (ab_list(iab)%neighbor_list_set)
    1120       567947 :          ab_list(iab)%nl_size = -1
    1121       567947 :          ab_list(iab)%nl_start = -1
    1122       567947 :          ab_list(iab)%nl_end = -1
    1123       692682 :          NULLIFY (ab_list(iab)%nlist_task)
    1124              :       END DO
    1125              : 
    1126              :       ! Allocate and initialize the kind availability
    1127       498940 :       ALLOCATE (pres_a(nkind), pres_b(nkind))
    1128       374398 :       DO ikind = 1, nkind
    1129       288445 :          pres_a(ikind) = ANY(pair_radius(ikind, :) > 0._dp)
    1130       427254 :          pres_b(ikind) = ANY(pair_radius(:, ikind) > 0._dp)
    1131              :       END DO
    1132              : 
    1133              :       ! create a copy of the pbc'ed coordinates
    1134       124735 :       natom = SIZE(particle_set)
    1135       374205 :       ALLOCATE (r_pbc(3, natom))
    1136       776252 :       DO i = 1, natom
    1137       776252 :          r_pbc(1:3, i) = pbc(particle_set(i)%r(1:3), cell)
    1138              :       END DO
    1139              : 
    1140              :       ! setup the local lists of atoms
    1141       124735 :       maxat = 0
    1142       374398 :       DO ikind = 1, nkind
    1143       374398 :          maxat = MAX(maxat, SIZE(atom(ikind)%list))
    1144              :       END DO
    1145       374205 :       ALLOCATE (index_list(maxat))
    1146       546456 :       DO i = 1, maxat
    1147       546456 :          index_list(i) = i
    1148              :       END DO
    1149       748410 :       ALLOCATE (lista(nkind), listb(nkind), nlista(nkind), nlistb(nkind))
    1150       374398 :       nlista = 0
    1151       374398 :       nlistb = 0
    1152       374398 :       DO ikind = 1, nkind
    1153       249663 :          NULLIFY (lista(ikind)%list, listb(ikind)%list)
    1154       124735 :          SELECT CASE (otype)
    1155              :          CASE (1)
    1156       148801 :             IF (ASSOCIATED(atom(ikind)%list_local_a_index)) THEN
    1157       104066 :                lista(ikind)%list => atom(ikind)%list_local_a_index
    1158       104066 :                nlista(ikind) = SIZE(lista(ikind)%list)
    1159              :             END IF
    1160       148801 :             IF (ASSOCIATED(atom(ikind)%list_local_b_index)) THEN
    1161       148743 :                listb(ikind)%list => atom(ikind)%list_local_b_index
    1162       148743 :                nlistb(ikind) = SIZE(listb(ikind)%list)
    1163              :             END IF
    1164              :          CASE (2)
    1165        20324 :             IF (ASSOCIATED(atom(ikind)%list_local_a_index)) THEN
    1166        13404 :                lista(ikind)%list => atom(ikind)%list_local_a_index
    1167        13404 :                nlista(ikind) = SIZE(lista(ikind)%list)
    1168              :             END IF
    1169        20324 :             nlistb(ikind) = SIZE(atom(ikind)%list)
    1170        20324 :             listb(ikind)%list => index_list
    1171              :          CASE (3)
    1172        19776 :             CALL combine_lists(lista(ikind)%list, nlista(ikind), ikind, atom)
    1173        19776 :             nlistb(ikind) = SIZE(atom(ikind)%list)
    1174        19776 :             listb(ikind)%list => index_list
    1175              :          CASE (4)
    1176        60762 :             nlista(ikind) = SIZE(atom(ikind)%list_1d)
    1177        60762 :             lista(ikind)%list => atom(ikind)%list_1d
    1178        60762 :             nlistb(ikind) = SIZE(atom(ikind)%list)
    1179        60762 :             listb(ikind)%list => index_list
    1180              :          CASE default
    1181       249663 :             CPABORT("")
    1182              :          END SELECT
    1183              :       END DO
    1184              : 
    1185              :       ! Determine max. number of local atoms
    1186       124735 :       maxat = 0
    1187       374398 :       DO ikind = 1, nkind
    1188       374398 :          maxat = MAX(maxat, nlista(ikind), nlistb(ikind))
    1189              :       END DO
    1190      1217647 :       ALLOCATE (kind_a(2*maxat))
    1191              : 
    1192              :       ! Load informations about the simulation cell
    1193       124735 :       CALL get_cell(cell=cell, periodic=periodic, deth=deth)
    1194              : 
    1195              :       ! Loop over all atomic kind pairs
    1196       374398 :       DO ikind = 1, nkind
    1197       249663 :          IF (.NOT. pres_a(ikind)) CYCLE
    1198              : 
    1199       892847 :          DO jkind = 1, nkind
    1200       532003 :             IF (.NOT. pres_b(jkind)) CYCLE
    1201              : 
    1202       514825 :             iab = ikind + nkind*(jkind - 1)
    1203              : 
    1204              :             ! Calculate the square of the maximum interaction distance
    1205       514825 :             IF (pair_radius(ikind, jkind) <= 0._dp) CYCLE
    1206       514781 :             rab_max = pair_radius(ikind, jkind)
    1207       514781 :             IF (otype == 3) THEN
    1208              :                ! Calculate the square of the maximum interaction distance
    1209              :                ! for sac_max / ncell this must be the maximum over all kinds
    1210              :                ! to be correct for three center terms involving different kinds
    1211       107888 :                rabm = MAXVAL(pair_radius(:, jkind))
    1212              :             ELSE
    1213              :                rabm = rab_max
    1214              :             END IF
    1215       514781 :             rab2_max = rabm*rabm
    1216              : 
    1217       514781 :             pd(1) = plane_distance(1, 0, 0, cell)
    1218       514781 :             pd(2) = plane_distance(0, 1, 0, cell)
    1219       514781 :             pd(3) = plane_distance(0, 0, 1, cell)
    1220              : 
    1221      2059124 :             sab_max = rabm/pd
    1222      2059124 :             sab_max_guard = 15.0_dp/pd
    1223              : 
    1224              :             ! It makes sense to have fewer subcells for larger systems
    1225       514781 :             subcell_scale = ((125.0_dp**3)/deth)**(1.0_dp/6.0_dp)
    1226              : 
    1227              :             ! guess the number of subcells for optimal performance,
    1228              :             ! guard against crazy stuff triggered by very small rabm
    1229              :             nsubcell(:) = INT(MAX(1.0_dp, MIN(0.5_dp*subcells*subcell_scale/sab_max(:), &
    1230      2059124 :                                               0.5_dp*subcells*subcell_scale/sab_max_guard(:))))
    1231              : 
    1232              :             ! number of image cells to be considered
    1233      2059124 :             ncell(:) = (INT(sab_max(:)) + 1)*periodic(:)
    1234              : 
    1235              :             CALL allocate_neighbor_list_set(neighbor_list_set=ab_list(iab)%neighbor_list_set, &
    1236       514781 :                                             symmetric=my_symmetric)
    1237       514781 :             neighbor_list_set => ab_list(iab)%neighbor_list_set
    1238              : 
    1239      1260579 :             DO iatom_local = 1, nlista(ikind)
    1240       745798 :                iatom = lista(ikind)%list(iatom_local)
    1241       745798 :                atom_a = atom(ikind)%list(iatom)
    1242              :                CALL add_neighbor_list(neighbor_list_set=neighbor_list_set, &
    1243              :                                       atom=atom_a, &
    1244      1260579 :                                       neighbor_list=kind_a(iatom_local)%neighbor_list)
    1245              :             END DO
    1246              : 
    1247       514781 :             CALL allocate_subcell(subcell, nsubcell)
    1248      1260579 :             DO iatom_local = 1, nlista(ikind)
    1249       745798 :                iatom = lista(ikind)%list(iatom_local)
    1250       745798 :                atom_a = atom(ikind)%list(iatom)
    1251      2983192 :                r = r_pbc(:, atom_a)
    1252       745798 :                CALL give_ijk_subcell(r, i, j, k, cell, nsubcell)
    1253      1260579 :                subcell(i, j, k)%natom = subcell(i, j, k)%natom + 1
    1254              :             END DO
    1255      1529768 :             DO k = 1, nsubcell(3)
    1256      3790057 :                DO j = 1, nsubcell(2)
    1257      9463517 :                   DO i = 1, nsubcell(1)
    1258      6188241 :                      maxat = subcell(i, j, k)%natom + subcell(i, j, k)%natom/10
    1259     12885746 :                      ALLOCATE (subcell(i, j, k)%atom_list(maxat))
    1260      8448530 :                      subcell(i, j, k)%natom = 0
    1261              :                   END DO
    1262              :                END DO
    1263              :             END DO
    1264      1260579 :             DO iatom_local = 1, nlista(ikind)
    1265       745798 :                iatom = lista(ikind)%list(iatom_local)
    1266       745798 :                atom_a = atom(ikind)%list(iatom)
    1267      2983192 :                r = r_pbc(:, atom_a)
    1268       745798 :                CALL give_ijk_subcell(r, i, j, k, cell, nsubcell)
    1269       745798 :                subcell(i, j, k)%natom = subcell(i, j, k)%natom + 1
    1270      1260579 :                subcell(i, j, k)%atom_list(subcell(i, j, k)%natom) = iatom_local
    1271              :             END DO
    1272              : 
    1273      1906021 :             DO jatom_local = 1, nlistb(jkind)
    1274      1391240 :                jatom = listb(jkind)%list(jatom_local)
    1275      1391240 :                atom_b = atom(jkind)%list(jatom)
    1276      1391240 :                IF (my_sort_atomb .AND. .NOT. my_symmetric) THEN
    1277         6662 :                   IF (.NOT. ANY(atomb_to_keep == atom_b)) CYCLE
    1278              :                END IF
    1279      1388630 :                IF (my_molecular) THEN
    1280         3268 :                   mol_b = atom(jkind)%list_b_mol(jatom_local)
    1281         3268 :                   IF (PRESENT(subset_of_mol)) THEN
    1282         1340 :                      IF (subset_of_mol(mol_b) .NE. current_subset) CYCLE
    1283              :                   END IF
    1284              :                END IF
    1285      5551736 :                r = r_pbc(:, atom_b)
    1286      1387934 :                CALL real_to_scaled(sb_pbc(:), r(:), cell)
    1287              : 
    1288      5249121 :                loop2_kcell: DO kcell = -ncell(3), ncell(3)
    1289      3635018 :                   sb(3) = sb_pbc(3) + REAL(kcell, dp)
    1290      3635018 :                   sb_min(3) = sb(3) - sab_max(3)
    1291      3635018 :                   sb_max(3) = sb(3) + sab_max(3)
    1292      3635018 :                   IF (periodic(3) /= 0) THEN
    1293      2889866 :                      IF (sb_min(3) >= 0.5_dp) EXIT loop2_kcell
    1294      2601254 :                      IF (sb_max(3) < -0.5_dp) CYCLE loop2_kcell
    1295              :                   END IF
    1296      3064608 :                   cell_b(3) = kcell
    1297              : 
    1298     17317989 :                   loop2_jcell: DO jcell = -ncell(2), ncell(2)
    1299     13779840 :                      sb(2) = sb_pbc(2) + REAL(jcell, dp)
    1300     13779840 :                      sb_min(2) = sb(2) - sab_max(2)
    1301     13779840 :                      sb_max(2) = sb(2) + sab_max(2)
    1302     13779840 :                      IF (periodic(2) /= 0) THEN
    1303     13033146 :                         IF (sb_min(2) >= 0.5_dp) EXIT loop2_jcell
    1304     12115447 :                         IF (sb_max(2) < -0.5_dp) CYCLE loop2_jcell
    1305              :                      END IF
    1306     11913612 :                      cell_b(2) = jcell
    1307              : 
    1308     90779542 :                      loop2_icell: DO icell = -ncell(1), ncell(1)
    1309     80371462 :                         sb(1) = sb_pbc(1) + REAL(icell, dp)
    1310     80371462 :                         sb_min(1) = sb(1) - sab_max(1)
    1311     80371462 :                         sb_max(1) = sb(1) + sab_max(1)
    1312     80371462 :                         IF (periodic(1) /= 0) THEN
    1313     79566500 :                            IF (sb_min(1) >= 0.5_dp) EXIT loop2_icell
    1314     74714562 :                            IF (sb_max(1) < -0.5_dp) CYCLE loop2_icell
    1315              :                         END IF
    1316     71503024 :                         cell_b(1) = icell
    1317              : 
    1318     71503024 :                         CALL scaled_to_real(rb, sb, cell)
    1319              : 
    1320    171923925 :                         loop_k: DO k = 1, nsubcell(3)
    1321    284140047 :                            loop_j: DO j = 1, nsubcell(2)
    1322    462028728 :                               loop_i: DO i = 1, nsubcell(1)
    1323              : 
    1324              :                                  ! FIXME for non-periodic systems, the whole subcell trick is skipped
    1325              :                                  ! yielding a Natom**2 pair list build.
    1326    270792425 :                                  IF (periodic(3) /= 0) THEN
    1327    260189905 :                                     IF (sb_max(3) < subcell(i, j, k)%s_min(3)) EXIT loop_k
    1328    257707718 :                                     IF (sb_min(3) >= subcell(i, j, k)%s_max(3)) CYCLE loop_k
    1329              :                                  END IF
    1330              : 
    1331    265030358 :                                  IF (periodic(2) /= 0) THEN
    1332    254470118 :                                     IF (sb_max(2) < subcell(i, j, k)%s_min(2)) EXIT loop_j
    1333    250681105 :                                     IF (sb_min(2) >= subcell(i, j, k)%s_max(2)) CYCLE loop_j
    1334              :                                  END IF
    1335              : 
    1336    255224918 :                                  IF (periodic(1) /= 0) THEN
    1337    244514522 :                                     IF (sb_max(1) < subcell(i, j, k)%s_min(1)) EXIT loop_i
    1338    235628916 :                                     IF (sb_min(1) >= subcell(i, j, k)%s_max(1)) CYCLE loop_i
    1339              :                                  END IF
    1340              : 
    1341    221462839 :                                  IF (subcell(i, j, k)%natom == 0) CYCLE
    1342              : 
    1343    442157089 :                                  DO iatom_subcell = 1, subcell(i, j, k)%natom
    1344    251968785 :                                     iatom_local = subcell(i, j, k)%atom_list(iatom_subcell)
    1345    251968785 :                                     iatom = lista(ikind)%list(iatom_local)
    1346    251968785 :                                     atom_a = atom(ikind)%list(iatom)
    1347    251968785 :                                     IF (my_molecular) THEN
    1348       476039 :                                        mol_a = atom(ikind)%list_a_mol(iatom_local)
    1349       476039 :                                        IF (mol_a /= mol_b) CYCLE
    1350              :                                     END IF
    1351    251732298 :                                     IF (my_symmetric) THEN
    1352    240020149 :                                        IF (atom_a > atom_b) THEN
    1353    111765981 :                                           include_ab = (MODULO(atom_a + atom_b, 2) /= 0)
    1354              :                                        ELSE
    1355    128254168 :                                           include_ab = (MODULO(atom_a + atom_b, 2) == 0)
    1356              :                                        END IF
    1357    240020149 :                                        IF (my_sort_atomb) THEN
    1358       665844 :                                           IF ((.NOT. ANY(atomb_to_keep == atom_b)) .AND. &
    1359              :                                               (.NOT. ANY(atomb_to_keep == atom_a))) THEN
    1360              :                                              include_ab = .FALSE.
    1361              :                                           END IF
    1362              :                                        END IF
    1363              :                                     ELSE
    1364              :                                        include_ab = .TRUE.
    1365              :                                     END IF
    1366    486250317 :                                     IF (include_ab) THEN
    1367    558334912 :                                        ra(:) = r_pbc(:, atom_a)
    1368    558334912 :                                        rab(:) = rb(:) - ra(:)
    1369    139583728 :                                        rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
    1370    139583728 :                                        IF (rab2 < rab2_max) THEN
    1371     41704735 :                                           include_ab = .TRUE.
    1372     41704735 :                                           IF (my_mic) THEN
    1373              :                                              ! only if rab is minimum image the pair will be included
    1374              :                                              ! ideally the range of the pair list is < L/2 so
    1375              :                                              ! that this never triggers
    1376      1332764 :                                              rab_pbc(:) = pbc(rab(:), cell)
    1377      5331056 :                                              IF (SUM((rab_pbc - rab)**2) > EPSILON(1.0_dp)) THEN
    1378              :                                                 include_ab = .FALSE.
    1379              :                                              END IF
    1380              :                                           END IF
    1381              :                                           IF (include_ab) THEN
    1382              :                                              CALL add_neighbor_node( &
    1383              :                                                 neighbor_list=kind_a(iatom_local)%neighbor_list, &
    1384              :                                                 neighbor=atom_b, &
    1385              :                                                 cell=cell_b, &
    1386              :                                                 r=rab, &
    1387     40683222 :                                                 nkind=nkind)
    1388              :                                           END IF
    1389              :                                        END IF
    1390              :                                     END IF
    1391              :                                  END DO
    1392              : 
    1393              :                               END DO loop_i
    1394              :                            END DO loop_j
    1395              :                         END DO loop_k
    1396              : 
    1397              :                      END DO loop2_icell
    1398              :                   END DO loop2_jcell
    1399              :                END DO loop2_kcell
    1400              : 
    1401              :             END DO
    1402              : 
    1403       781666 :             CALL deallocate_subcell(subcell)
    1404              : 
    1405              :          END DO
    1406              :       END DO
    1407              : 
    1408        10357 :       SELECT CASE (otype)
    1409              :       CASE (1:2, 4)
    1410              :       CASE (3)
    1411        30133 :          DO ikind = 1, nkind
    1412        30133 :             DEALLOCATE (lista(ikind)%list)
    1413              :          END DO
    1414              :       CASE default
    1415       124735 :          CPABORT("")
    1416              :       END SELECT
    1417       124735 :       DEALLOCATE (kind_a, pres_a, pres_b, lista, listb, nlista, nlistb)
    1418       124735 :       DEALLOCATE (index_list)
    1419       124735 :       DEALLOCATE (r_pbc)
    1420              : 
    1421       124735 :       nentry = 0
    1422       124735 :       CALL neighbor_list_iterator_create(nl_iterator, ab_list)
    1423     40807957 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
    1424     40683222 :          CALL get_iterator_info(nl_iterator, inode=inode, nnode=nnode)
    1425     40807957 :          IF (inode == 1) nentry = nentry + nnode
    1426              :       END DO
    1427       124735 :       CALL neighbor_list_iterator_release(nl_iterator)
    1428              :       !
    1429     41798858 :       ALLOCATE (ab_list(1)%nlist_task(nentry))
    1430       124735 :       ab_list(1)%nl_size = nentry
    1431       567947 :       DO iab = 2, SIZE(ab_list)
    1432       443212 :          ab_list(iab)%nl_size = nentry
    1433       567947 :          ab_list(iab)%nlist_task => ab_list(1)%nlist_task
    1434              :       END DO
    1435              :       !
    1436       124735 :       nentry = 0
    1437       124735 :       CALL neighbor_list_iterator_create(nl_iterator, ab_list)
    1438     40807957 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
    1439     40683222 :          nentry = nentry + 1
    1440     40683222 :          CALL get_iterator_task(nl_iterator, ab_list(1)%nlist_task(nentry))
    1441     40683222 :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, nkind=nkind)
    1442     40683222 :          iab = (ikind - 1)*nkind + jkind
    1443     40683222 :          IF (ab_list(iab)%nl_start < 0) ab_list(iab)%nl_start = nentry
    1444     40807957 :          IF (ab_list(iab)%nl_end < 0) THEN
    1445       338143 :             ab_list(iab)%nl_end = nentry
    1446              :          ELSE
    1447     40345079 :             CPASSERT(ab_list(iab)%nl_end + 1 == nentry)
    1448     40345079 :             ab_list(iab)%nl_end = nentry
    1449              :          END IF
    1450              :       END DO
    1451       124735 :       CALL neighbor_list_iterator_release(nl_iterator)
    1452              : 
    1453       124735 :       CALL timestop(handle)
    1454              : 
    1455       249470 :    END SUBROUTINE build_neighbor_lists
    1456              : 
    1457              : ! **************************************************************************************************
    1458              : !> \brief Build a neighborlist
    1459              : !> \param ab_list ...
    1460              : !> \param basis_set_a ...
    1461              : !> \param basis_set_b ...
    1462              : !> \param qs_env ...
    1463              : !> \param mic ...
    1464              : !> \param symmetric ...
    1465              : !> \param molecular ...
    1466              : !> \param operator_type ...
    1467              : !> \date    14.03.2016
    1468              : !> \author  JGH
    1469              : ! **************************************************************************************************
    1470          110 :    SUBROUTINE setup_neighbor_list(ab_list, basis_set_a, basis_set_b, qs_env, &
    1471              :                                   mic, symmetric, molecular, operator_type)
    1472              : 
    1473              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1474              :          POINTER                                         :: ab_list
    1475              :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_a
    1476              :       TYPE(gto_basis_set_p_type), DIMENSION(:), &
    1477              :          OPTIONAL, POINTER                               :: basis_set_b
    1478              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1479              :       LOGICAL, INTENT(IN), OPTIONAL                      :: mic, symmetric, molecular
    1480              :       CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: operator_type
    1481              : 
    1482              :       CHARACTER(LEN=4)                                   :: otype
    1483              :       INTEGER                                            :: ikind, nkind
    1484              :       LOGICAL                                            :: my_mic, my_molecular, my_symmetric
    1485          110 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: a_present, b_present
    1486          110 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: a_radius, b_radius
    1487          110 :       REAL(dp), ALLOCATABLE, DIMENSION(:, :)             :: pair_radius
    1488          110 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1489              :       TYPE(cell_type), POINTER                           :: cell
    1490              :       TYPE(distribution_1d_type), POINTER                :: distribution_1d
    1491              :       TYPE(distribution_2d_type), POINTER                :: distribution_2d
    1492          110 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_a, basis_b
    1493              :       TYPE(gto_basis_set_type), POINTER                  :: abas, bbas
    1494          110 :       TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:)  :: atom2d
    1495          110 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
    1496          110 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1497              : 
    1498          110 :       basis_a => basis_set_a
    1499          110 :       IF (PRESENT(basis_set_b)) THEN
    1500           54 :          basis_b => basis_set_b
    1501           54 :          my_symmetric = .FALSE.
    1502              :       ELSE
    1503           56 :          basis_b => basis_set_a
    1504           56 :          my_symmetric = .TRUE.
    1505              :       END IF
    1506          110 :       IF (PRESENT(symmetric)) my_symmetric = symmetric
    1507              : 
    1508          110 :       IF (PRESENT(mic)) THEN
    1509            6 :          my_mic = mic
    1510              :       ELSE
    1511          104 :          my_mic = .FALSE.
    1512              :       END IF
    1513              : 
    1514          110 :       IF (PRESENT(molecular)) THEN
    1515            8 :          my_molecular = molecular
    1516              :       ELSE
    1517          102 :          my_molecular = .FALSE.
    1518              :       END IF
    1519              : 
    1520          110 :       IF (PRESENT(operator_type)) THEN
    1521            0 :          otype = operator_type
    1522              :       ELSE
    1523              :          ! default is a simple AB neighbor list
    1524          110 :          otype = "AB"
    1525              :       END IF
    1526              : 
    1527          110 :       nkind = SIZE(basis_a)
    1528          440 :       ALLOCATE (a_present(nkind), b_present(nkind))
    1529          338 :       a_present = .FALSE.
    1530          338 :       b_present = .FALSE.
    1531          440 :       ALLOCATE (a_radius(nkind), b_radius(nkind))
    1532          338 :       a_radius = 0.0_dp
    1533          338 :       b_radius = 0.0_dp
    1534          338 :       DO ikind = 1, nkind
    1535          228 :          IF (ASSOCIATED(basis_a(ikind)%gto_basis_set)) THEN
    1536          228 :             a_present(ikind) = .TRUE.
    1537          228 :             abas => basis_a(ikind)%gto_basis_set
    1538          228 :             CALL get_gto_basis_set(gto_basis_set=abas, kind_radius=a_radius(ikind))
    1539              :          END IF
    1540          338 :          IF (ASSOCIATED(basis_b(ikind)%gto_basis_set)) THEN
    1541          228 :             b_present(ikind) = .TRUE.
    1542          228 :             bbas => basis_b(ikind)%gto_basis_set
    1543          228 :             CALL get_gto_basis_set(gto_basis_set=bbas, kind_radius=b_radius(ikind))
    1544              :          END IF
    1545              :       END DO
    1546              : 
    1547          440 :       ALLOCATE (pair_radius(nkind, nkind))
    1548          826 :       pair_radius = 0.0_dp
    1549          110 :       CALL pair_radius_setup(a_present, b_present, a_radius, b_radius, pair_radius)
    1550              : 
    1551              :       CALL get_qs_env(qs_env, &
    1552              :                       atomic_kind_set=atomic_kind_set, &
    1553              :                       cell=cell, &
    1554              :                       distribution_2d=distribution_2d, &
    1555              :                       local_particles=distribution_1d, &
    1556              :                       particle_set=particle_set, &
    1557          110 :                       molecule_set=molecule_set)
    1558              : 
    1559          558 :       ALLOCATE (atom2d(nkind))
    1560              :       CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
    1561          110 :                         molecule_set, my_molecular, particle_set=particle_set)
    1562              :       CALL build_neighbor_lists(ab_list, particle_set, atom2d, cell, pair_radius, &
    1563              :                                 mic=my_mic, symmetric=my_symmetric, molecular=my_molecular, &
    1564          110 :                                 subcells=2.0_dp, nlname="AUX_NL")
    1565              : 
    1566          110 :       CALL atom2d_cleanup(atom2d)
    1567              : 
    1568          110 :       DEALLOCATE (a_present, b_present, a_radius, b_radius, pair_radius, atom2d)
    1569              : 
    1570          110 :    END SUBROUTINE setup_neighbor_list
    1571              : 
    1572              : ! **************************************************************************************************
    1573              : !> \brief ...
    1574              : !> \param list ...
    1575              : !> \param n ...
    1576              : !> \param ikind ...
    1577              : !> \param atom ...
    1578              : ! **************************************************************************************************
    1579        19776 :    SUBROUTINE combine_lists(list, n, ikind, atom)
    1580              :       INTEGER, DIMENSION(:), POINTER                     :: list
    1581              :       INTEGER, INTENT(OUT)                               :: n
    1582              :       INTEGER, INTENT(IN)                                :: ikind
    1583              :       TYPE(local_atoms_type), DIMENSION(:), INTENT(IN)   :: atom
    1584              : 
    1585              :       INTEGER                                            :: i, ib, na, nb
    1586        19776 :       INTEGER, DIMENSION(:), POINTER                     :: lista, listb
    1587              : 
    1588            0 :       CPASSERT(.NOT. ASSOCIATED(list))
    1589              : 
    1590        19776 :       lista => atom(ikind)%list_local_a_index
    1591        19776 :       listb => atom(ikind)%list_local_b_index
    1592              : 
    1593        19776 :       IF (ASSOCIATED(lista)) THEN
    1594        12465 :          na = SIZE(lista)
    1595              :       ELSE
    1596              :          na = 0
    1597              :       END IF
    1598              : 
    1599        19776 :       IF (ASSOCIATED(listb)) THEN
    1600        19776 :          nb = SIZE(listb)
    1601              :       ELSE
    1602              :          nb = 0
    1603              :       END IF
    1604              : 
    1605        59328 :       ALLOCATE (list(na + nb))
    1606              : 
    1607        19776 :       n = na
    1608        71883 :       IF (na .GT. 0) list(1:na) = lista(1:na)
    1609        19776 :       IF (nb .GT. 0) THEN
    1610        58809 :          loopb: DO ib = 1, nb
    1611        81693 :             DO i = 1, na
    1612        81693 :                IF (listb(ib) == list(i)) CYCLE loopb
    1613              :             END DO
    1614        19212 :             n = n + 1
    1615        58809 :             list(n) = listb(ib)
    1616              :          END DO loopb
    1617              :       END IF
    1618        19776 :    END SUBROUTINE combine_lists
    1619              : 
    1620              : ! **************************************************************************************************
    1621              : 
    1622              : ! **************************************************************************************************
    1623              : !> \brief ...
    1624              : !> \param present_a ...
    1625              : !> \param present_b ...
    1626              : !> \param radius_a ...
    1627              : !> \param radius_b ...
    1628              : !> \param pair_radius ...
    1629              : !> \param prmin ...
    1630              : ! **************************************************************************************************
    1631       107804 :    SUBROUTINE pair_radius_setup(present_a, present_b, radius_a, radius_b, pair_radius, prmin)
    1632              :       LOGICAL, DIMENSION(:), INTENT(IN)                  :: present_a, present_b
    1633              :       REAL(dp), DIMENSION(:), INTENT(IN)                 :: radius_a, radius_b
    1634              :       REAL(dp), DIMENSION(:, :), INTENT(OUT)             :: pair_radius
    1635              :       REAL(dp), INTENT(IN), OPTIONAL                     :: prmin
    1636              : 
    1637              :       INTEGER                                            :: i, j, nkind
    1638              :       REAL(dp)                                           :: rrmin
    1639              : 
    1640       107804 :       nkind = SIZE(present_a)
    1641              : 
    1642       814208 :       pair_radius = 0._dp
    1643              : 
    1644       107804 :       rrmin = 0.0_dp
    1645       107804 :       IF (PRESENT(prmin)) rrmin = prmin
    1646              : 
    1647       323541 :       DO i = 1, nkind
    1648       215737 :          IF (.NOT. present_a(i)) CYCLE
    1649       765154 :          DO j = 1, nkind
    1650       455033 :             IF (.NOT. present_b(j)) CYCLE
    1651       437551 :             pair_radius(i, j) = radius_a(i) + radius_b(j)
    1652       670770 :             pair_radius(i, j) = MAX(pair_radius(i, j), rrmin)
    1653              :          END DO
    1654              :       END DO
    1655              : 
    1656       107804 :    END SUBROUTINE pair_radius_setup
    1657              : 
    1658              : ! **************************************************************************************************
    1659              : !> \brief   Print the distribution of the simple pair neighbor list.
    1660              : !> \param ab ...
    1661              : !> \param qs_kind_set ...
    1662              : !> \param output_unit ...
    1663              : !> \param para_env ...
    1664              : !> \date    19.06.2003
    1665              : !> \author  MK
    1666              : !> \version 1.0
    1667              : ! **************************************************************************************************
    1668          158 :    SUBROUTINE write_neighbor_distribution(ab, qs_kind_set, output_unit, para_env)
    1669              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1670              :          POINTER                                         :: ab
    1671              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1672              :       INTEGER, INTENT(in)                                :: output_unit
    1673              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1674              : 
    1675              :       CHARACTER(len=*), PARAMETER :: routineN = 'write_neighbor_distribution'
    1676              :       LOGICAL, PARAMETER                                 :: full_output = .FALSE.
    1677              : 
    1678              :       INTEGER                                            :: handle, ikind, inode, ipe, jkind, n, &
    1679              :                                                             nkind, nnode
    1680              :       INTEGER(int_8)                                     :: nblock_max, nblock_sum, nelement_max, &
    1681              :                                                             nelement_sum, tmp(2)
    1682          158 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: nblock, nelement, nnsgf
    1683              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
    1684              :       TYPE(neighbor_list_iterator_p_type), &
    1685          158 :          DIMENSION(:), POINTER                           :: nl_iterator
    1686              : 
    1687          158 :       CALL timeset(routineN, handle)
    1688              :       ASSOCIATE (mype => para_env%mepos + 1, npe => para_env%num_pe)
    1689              : 
    1690              :          ! Allocate work storage
    1691          632 :          ALLOCATE (nblock(npe), nelement(npe))
    1692          474 :          nblock(:) = 0
    1693          474 :          nelement(:) = 0
    1694          158 :          nkind = SIZE(qs_kind_set)
    1695          474 :          ALLOCATE (nnsgf(nkind))
    1696          448 :          nnsgf = 1
    1697          448 :          DO ikind = 1, nkind
    1698          290 :             CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
    1699          448 :             IF (ASSOCIATED(orb_basis_set)) THEN
    1700          238 :                CALL get_gto_basis_set(gto_basis_set=orb_basis_set, nsgf=nnsgf(ikind))
    1701              :             END IF
    1702              :          END DO
    1703              : 
    1704          158 :          CALL neighbor_list_iterator_create(nl_iterator, ab)
    1705        43688 :          DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
    1706        43530 :             CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, nnode=nnode)
    1707        43688 :             IF (inode == 1) THEN
    1708         1151 :                n = nnsgf(ikind)*nnsgf(jkind)
    1709         1151 :                nblock(mype) = nblock(mype) + nnode
    1710         1151 :                nelement(mype) = nelement(mype) + n*nnode
    1711              :             END IF
    1712              :          END DO
    1713          158 :          CALL neighbor_list_iterator_release(nl_iterator)
    1714              : 
    1715              :          IF (full_output) THEN
    1716              :             ! XXXXXXXX should gather/scatter this on ionode
    1717              :             CALL para_env%sum(nblock)
    1718              :             CALL para_env%sum(nelement)
    1719              : 
    1720              :             nblock_sum = SUM(INT(nblock, KIND=int_8))
    1721              :             nelement_sum = SUM(INT(nelement, KIND=int_8))
    1722              :          ELSE
    1723          158 :             nblock_sum = nblock(mype)
    1724              :             nblock_max = nblock(mype)
    1725          158 :             nelement_sum = nelement(mype)
    1726              :             nelement_max = nelement(mype)
    1727          474 :             tmp = (/nblock_sum, nelement_sum/)
    1728          158 :             CALL para_env%sum(tmp)
    1729          158 :             nblock_sum = tmp(1); nelement_sum = tmp(2)
    1730          474 :             tmp = (/nblock_max, nelement_max/)
    1731          158 :             CALL para_env%max(tmp)
    1732          158 :             nblock_max = tmp(1); nelement_max = tmp(2)
    1733              :          END IF
    1734              : 
    1735          316 :          IF (output_unit > 0) THEN
    1736              :             IF (full_output) THEN
    1737              :                WRITE (UNIT=output_unit, &
    1738              :                       FMT="(/,/,T2,A,/,/,T3,A,/,/,(T4,I6,T27,I10,T55,I10))") &
    1739              :                   "DISTRIBUTION OF THE NEIGHBOR LISTS", &
    1740              :                   "Process   Number of particle pairs   Number of matrix elements", &
    1741              :                   (ipe - 1, nblock(ipe), nelement(ipe), ipe=1, npe)
    1742              :                WRITE (UNIT=output_unit, FMT="(/,T7,A3,T27,I10,T55,I10)") &
    1743              :                   "Sum", SUM(nblock), SUM(nelement)
    1744              :             ELSE
    1745           79 :                WRITE (UNIT=output_unit, FMT="(/,T2,A)") "DISTRIBUTION OF THE NEIGHBOR LISTS"
    1746           79 :                WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Total number of particle pairs:", nblock_sum
    1747           79 :                WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Total number of matrix elements:", nelement_sum
    1748           79 :                WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Average number of particle pairs:", (nblock_sum + npe - 1)/npe
    1749           79 :                WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Maximum number of particle pairs:", nblock_max
    1750           79 :                WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Average number of matrix element:", (nelement_sum + npe - 1)/npe
    1751           79 :                WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Maximum number of matrix elements:", nelement_max
    1752              :             END IF
    1753              :          END IF
    1754              :       END ASSOCIATE
    1755              : 
    1756              :       ! Release work storage
    1757              : 
    1758          158 :       DEALLOCATE (nblock, nelement, nnsgf)
    1759              : 
    1760          158 :       CALL timestop(handle)
    1761              : 
    1762          158 :    END SUBROUTINE write_neighbor_distribution
    1763              : 
    1764              : ! **************************************************************************************************
    1765              : !> \brief   Write a set of neighbor lists to the output unit.
    1766              : !> \param ab ...
    1767              : !> \param particle_set ...
    1768              : !> \param cell ...
    1769              : !> \param para_env ...
    1770              : !> \param neighbor_list_section ...
    1771              : !> \param nl_type ...
    1772              : !> \param middle_name ...
    1773              : !> \param nlname ...
    1774              : !> \date    04.03.2002
    1775              : !> \par History
    1776              : !>       - Adapted to the new parallelized neighbor list version
    1777              : !>         (26.06.2003,MK)
    1778              : !> \author  MK
    1779              : !> \version 1.0
    1780              : ! **************************************************************************************************
    1781        69684 :    SUBROUTINE write_neighbor_lists(ab, particle_set, cell, para_env, neighbor_list_section, &
    1782              :                                    nl_type, middle_name, nlname)
    1783              : 
    1784              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1785              :          POINTER                                         :: ab
    1786              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1787              :       TYPE(cell_type), POINTER                           :: cell
    1788              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1789              :       TYPE(section_vals_type), POINTER                   :: neighbor_list_section
    1790              :       CHARACTER(LEN=*), INTENT(IN)                       :: nl_type, middle_name, nlname
    1791              : 
    1792              :       CHARACTER(LEN=default_string_length)               :: string, unit_str
    1793              :       INTEGER                                            :: iatom, inode, iw, jatom, nneighbor, nnode
    1794              :       INTEGER, DIMENSION(3)                              :: cell_b
    1795              :       REAL(dp)                                           :: dab, unit_conv
    1796              :       REAL(dp), DIMENSION(3)                             :: ra, rab, rb
    1797              :       TYPE(cp_logger_type), POINTER                      :: logger
    1798              :       TYPE(neighbor_list_iterator_p_type), &
    1799        69684 :          DIMENSION(:), POINTER                           :: nl_iterator
    1800              : 
    1801        69684 :       NULLIFY (logger)
    1802        69684 :       logger => cp_get_default_logger()
    1803        69684 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, neighbor_list_section, &
    1804              :                                            TRIM(nl_type)), &
    1805              :                 cp_p_file)) THEN
    1806              :          iw = cp_print_key_unit_nr(logger=logger, &
    1807              :                                    basis_section=neighbor_list_section, &
    1808              :                                    print_key_path=TRIM(nl_type), &
    1809              :                                    extension=".out", &
    1810              :                                    middle_name=TRIM(middle_name), &
    1811              :                                    local=.TRUE., &
    1812              :                                    log_filename=.FALSE., &
    1813            4 :                                    file_position="REWIND")
    1814              :          ASSOCIATE (mype => para_env%mepos)
    1815            4 :             CALL section_vals_val_get(neighbor_list_section, "UNIT", c_val=unit_str)
    1816            4 :             unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
    1817              : 
    1818              :             ! Print headline
    1819            4 :             string = ""
    1820              :             WRITE (UNIT=string, FMT="(A,I5,A)") &
    1821            4 :                TRIM(nlname)//" IN "//TRIM(unit_str)//" (PROCESS", mype, ")"
    1822            4 :             CALL compress(string)
    1823            4 :             IF (iw > 0) WRITE (UNIT=iw, FMT="(/,/,T2,A)") TRIM(string)
    1824              : 
    1825            4 :             nneighbor = 0
    1826              : 
    1827            4 :             CALL neighbor_list_iterator_create(nl_iterator, ab)
    1828           16 :             DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
    1829              :                CALL get_iterator_info(nl_iterator, inode=inode, nnode=nnode, &
    1830           12 :                                       iatom=iatom, jatom=jatom, cell=cell_b, r=rab)
    1831           12 :                nneighbor = nneighbor + 1
    1832           12 :                ra(:) = pbc(particle_set(iatom)%r, cell)
    1833           48 :                rb(:) = ra(:) + rab(:)
    1834           12 :                dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
    1835           16 :                IF (iw > 0) THEN
    1836           12 :                   IF (inode == 1) THEN
    1837              :                      WRITE (UNIT=iw, FMT="(/,T2,I5,3X,I6,3X,3F12.6)") &
    1838           40 :                         iatom, nnode, ra(1:3)*unit_conv
    1839              :                   END IF
    1840              :                   WRITE (UNIT=iw, FMT="(T10,I6,3X,3I4,3F12.6,2X,F12.6)") &
    1841           60 :                      jatom, cell_b(1:3), rb(1:3)*unit_conv, dab*unit_conv
    1842              :                END IF
    1843              :             END DO
    1844            4 :             CALL neighbor_list_iterator_release(nl_iterator)
    1845              : 
    1846            4 :             string = ""
    1847              :             WRITE (UNIT=string, FMT="(A,I12,A,I12)") &
    1848            4 :                "Total number of neighbor interactions for process", mype, ":", &
    1849            8 :                nneighbor
    1850            4 :             CALL compress(string)
    1851            4 :             IF (iw > 0) WRITE (UNIT=iw, FMT="(/,T2,A)") TRIM(string)
    1852              :             CALL cp_print_key_finished_output(unit_nr=iw, &
    1853              :                                               logger=logger, &
    1854              :                                               basis_section=neighbor_list_section, &
    1855              :                                               print_key_path=TRIM(nl_type), &
    1856            8 :                                               local=.TRUE.)
    1857              :          END ASSOCIATE
    1858              :       END IF
    1859              : 
    1860        69684 :    END SUBROUTINE write_neighbor_lists
    1861              : 
    1862            0 : END MODULE qs_neighbor_lists
        

Generated by: LCOV version 2.0-1