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