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