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