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