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