Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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_dftd3, vdw_pairpot_dftd3bj, &
52 : 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, &
106 : list_local_a_index, &
107 : list_local_b_index, &
108 : list_1d, &
109 : list_a_mol, &
110 : list_b_mol
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 36365 : 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 36365 : CALL timeset(routineN, handle)
138 104502 : DO ikind = 1, SIZE(atom2d)
139 68137 : NULLIFY (atom2d(ikind)%list)
140 68137 : IF (ASSOCIATED(atom2d(ikind)%list_local_a_index)) THEN
141 49181 : DEALLOCATE (atom2d(ikind)%list_local_a_index)
142 : END IF
143 68137 : IF (ASSOCIATED(atom2d(ikind)%list_local_b_index)) THEN
144 68079 : DEALLOCATE (atom2d(ikind)%list_local_b_index)
145 : END IF
146 68137 : IF (ASSOCIATED(atom2d(ikind)%list_a_mol)) THEN
147 49181 : DEALLOCATE (atom2d(ikind)%list_a_mol)
148 : END IF
149 68137 : IF (ASSOCIATED(atom2d(ikind)%list_b_mol)) THEN
150 68079 : DEALLOCATE (atom2d(ikind)%list_b_mol)
151 : END IF
152 104502 : IF (ASSOCIATED(atom2d(ikind)%list_1d)) THEN
153 68137 : DEALLOCATE (atom2d(ikind)%list_1d)
154 : END IF
155 : END DO
156 36365 : CALL timestop(handle)
157 :
158 36365 : 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 36365 : 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 36365 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom2mol, atom_of_kind, listindex, &
188 36365 : listsort
189 36365 : INTEGER, DIMENSION(:), POINTER :: local_cols_array, local_rows_array
190 :
191 36365 : CALL timeset(routineN, handle)
192 :
193 36365 : nkind = SIZE(atomic_kind_set)
194 36365 : natom = SIZE(particle_set)
195 36365 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
196 :
197 36365 : 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 104502 : DO ikind = 1, nkind
207 68137 : NULLIFY (atom2d(ikind)%list)
208 68137 : NULLIFY (atom2d(ikind)%list_local_a_index)
209 68137 : NULLIFY (atom2d(ikind)%list_local_b_index)
210 68137 : NULLIFY (atom2d(ikind)%list_1d)
211 68137 : NULLIFY (atom2d(ikind)%list_a_mol)
212 68137 : NULLIFY (atom2d(ikind)%list_b_mol)
213 :
214 68137 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
215 :
216 68137 : natom_a = SIZE(atom2d(ikind)%list)
217 :
218 68137 : natom_local_a = distribution_2d%n_local_rows(ikind)
219 68137 : natom_local_b = distribution_2d%n_local_cols(ikind)
220 68137 : local_rows_array => distribution_2d%local_rows(ikind)%array
221 68137 : local_cols_array => distribution_2d%local_cols(ikind)%array
222 :
223 68137 : nel = distribution_1d%n_el(ikind)
224 184454 : ALLOCATE (atom2d(ikind)%list_1d(nel))
225 157835 : DO iat = 1, nel
226 89698 : ia = distribution_1d%list(ikind)%array(iat)
227 157835 : atom2d(ikind)%list_1d(iat) = atom_of_kind(ia)
228 : END DO
229 :
230 340685 : ALLOCATE (listsort(natom_a), listindex(natom_a))
231 243709 : listsort(1:natom_a) = atom2d(ikind)%list(1:natom_a)
232 68137 : CALL sort(listsort, natom_a, listindex)
233 : ! Block rows
234 68137 : IF (natom_local_a > 0) THEN
235 147543 : ALLOCATE (atom2d(ikind)%list_local_a_index(natom_local_a))
236 147543 : ALLOCATE (atom2d(ikind)%list_a_mol(natom_local_a))
237 146669 : atom2d(ikind)%list_a_mol(:) = 0
238 :
239 : ! Build index vector for mapping
240 146669 : DO iatom_local = 1, natom_local_a
241 97488 : atom_a = local_rows_array(iatom_local)
242 97488 : iatom = locate(listsort, atom_a)
243 97488 : atom2d(ikind)%list_local_a_index(iatom_local) = listindex(iatom)
244 146669 : 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 68137 : IF (natom_local_b > 0) THEN
251 :
252 204237 : ALLOCATE (atom2d(ikind)%list_local_b_index(natom_local_b))
253 204237 : ALLOCATE (atom2d(ikind)%list_b_mol(natom_local_b))
254 243515 : atom2d(ikind)%list_b_mol(:) = 0
255 :
256 : ! Build index vector for mapping
257 243515 : DO iatom_local = 1, natom_local_b
258 175436 : atom_a = local_cols_array(iatom_local)
259 175436 : iatom = locate(listsort, atom_a)
260 175436 : atom2d(ikind)%list_local_b_index(iatom_local) = listindex(iatom)
261 243515 : IF (molecule_only) atom2d(ikind)%list_b_mol(iatom_local) = atom2mol(atom_a)
262 : END DO
263 :
264 : END IF
265 :
266 104502 : DEALLOCATE (listsort, listindex)
267 :
268 : END DO
269 :
270 36365 : CALL timestop(handle)
271 :
272 72730 : 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 22539 : 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 22539 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: all_present, aux_fit_present, aux_present, &
302 22539 : core_present, default_present, nonbond1_atom, nonbond2_atom, oce_present, orb_present, &
303 22539 : 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 22539 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: all_pot_rad, aux_fit_radius, c_radius, calpha, &
307 22539 : core_radius, oce_radius, orb_radius, ppl_radius, ppnl_radius, ri_radius, zeff
308 22539 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius, pair_radius_lb
309 : TYPE(all_potential_type), POINTER :: all_potential
310 22539 : 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 22539 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
322 22539 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
323 22539 : TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: saa_list, sab_all, sab_almo, &
324 22539 : sab_cn, sab_core, sab_gcp, sab_kp, sab_kp_nosym, sab_lrc, sab_orb, sab_scp, sab_se, &
325 22539 : sab_tbe, sab_vdw, sab_xb, sab_xtb_nonbond, sab_xtbe, sac_ae, sac_lri, sac_ppl, sap_oce, &
326 22539 : sap_ppnl, soa_list, soo_list
327 22539 : 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 22539 : 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 22539 : CALL timeset(routineN, handle)
339 22539 : NULLIFY (logger)
340 22539 : logger => cp_get_default_logger()
341 :
342 22539 : NULLIFY (atomic_kind_set, qs_kind_set, cell, neighbor_list_section, &
343 22539 : distribution_1d, distribution_2d, gth_potential, sgp_potential, orb_basis_set, &
344 22539 : particle_set, molecule_set, dft_control, ks_env)
345 :
346 22539 : NULLIFY (sab_orb)
347 22539 : NULLIFY (sac_ae)
348 22539 : NULLIFY (sac_ppl)
349 22539 : NULLIFY (sac_lri)
350 22539 : NULLIFY (sap_ppnl)
351 22539 : NULLIFY (sap_oce)
352 22539 : NULLIFY (sab_se)
353 22539 : NULLIFY (sab_lrc)
354 22539 : NULLIFY (sab_tbe)
355 22539 : NULLIFY (sab_xtbe)
356 22539 : NULLIFY (sab_core)
357 22539 : NULLIFY (sab_xb)
358 22539 : NULLIFY (sab_xtb_nonbond)
359 22539 : NULLIFY (sab_all)
360 22539 : NULLIFY (sab_vdw)
361 22539 : NULLIFY (sab_cn)
362 22539 : NULLIFY (soo_list)
363 22539 : NULLIFY (sab_scp)
364 22539 : NULLIFY (sab_almo)
365 22539 : NULLIFY (sab_kp)
366 22539 : 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 22539 : dft_control=dft_control)
379 :
380 22539 : 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 22539 : last_qs_neighbor_list_id_nr = last_qs_neighbor_list_id_nr + 1
385 22539 : 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_nonbond=sab_xtb_nonbond, &
402 : sab_scp=sab_scp, &
403 : sab_all=sab_all, &
404 : sab_almo=sab_almo, &
405 : sab_kp=sab_kp, &
406 22539 : sab_kp_nosym=sab_kp_nosym)
407 :
408 22539 : dokp = (kpoints%nkp > 0)
409 22539 : nddo = dft_control%qs_control%semi_empirical
410 22539 : dftb = dft_control%qs_control%dftb
411 22539 : xtb = dft_control%qs_control%xtb
412 22539 : almo = dft_control%qs_control%do_almo_scf
413 22539 : lrigpw = (dft_control%qs_control%method_id == do_method_lrigpw)
414 22539 : rigpw = (dft_control%qs_control%method_id == do_method_rigpw)
415 22539 : lri_optbas = dft_control%qs_control%lri_optbas
416 :
417 : ! molecular lists
418 22539 : molecule_only = .FALSE.
419 22539 : IF (PRESENT(molecular)) molecule_only = molecular
420 : ! minimum image convention (MIC)
421 22539 : mic = molecule_only
422 22539 : IF (dokp) THEN
423 : ! no MIC for kpoints
424 764 : mic = .FALSE.
425 21775 : ELSEIF (nddo) THEN
426 : ! enforce MIC for interaction lists in SE
427 5854 : mic = .TRUE.
428 : END IF
429 22539 : pdist = dft_control%qs_control%pairlist_radius
430 :
431 22539 : hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
432 22539 : CALL section_vals_get(hfx_sections, explicit=do_hfx)
433 :
434 22539 : CALL get_atomic_kind_set(atomic_kind_set, maxatom=maxatom)
435 : CALL get_qs_kind_set(qs_kind_set, paw_atom_present=paw_atom_present, &
436 : gth_potential_present=gth_potential_present, &
437 : sgp_potential_present=sgp_potential_present, &
438 22539 : all_potential_present=all_potential_present)
439 :
440 22539 : CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
441 :
442 : ! Allocate work storage
443 22539 : nkind = SIZE(atomic_kind_set)
444 : ALLOCATE (orb_present(nkind), aux_fit_present(nkind), aux_present(nkind), &
445 247929 : default_present(nkind), core_present(nkind))
446 : ALLOCATE (orb_radius(nkind), aux_fit_radius(nkind), c_radius(nkind), &
447 293007 : core_radius(nkind), calpha(nkind), zeff(nkind))
448 68095 : orb_radius(:) = 0.0_dp
449 68095 : aux_fit_radius(:) = 0.0_dp
450 68095 : c_radius(:) = 0.0_dp
451 68095 : core_radius(:) = 0.0_dp
452 68095 : calpha(:) = 0.0_dp
453 68095 : zeff(:) = 0.0_dp
454 :
455 90156 : ALLOCATE (pair_radius(nkind, nkind))
456 22539 : IF (gth_potential_present .OR. sgp_potential_present) THEN
457 47935 : ALLOCATE (ppl_present(nkind), ppl_radius(nkind))
458 26693 : ppl_radius = 0.0_dp
459 47935 : ALLOCATE (ppnl_present(nkind), ppnl_radius(nkind))
460 26693 : ppnl_radius = 0.0_dp
461 : END IF
462 22539 : IF (paw_atom_present) THEN
463 7800 : ALLOCATE (oce_present(nkind), oce_radius(nkind))
464 4638 : oce_radius = 0.0_dp
465 : END IF
466 22539 : IF (all_potential_present .OR. sgp_potential_present) THEN
467 65000 : ALLOCATE (all_present(nkind), all_pot_rad(nkind))
468 41550 : all_pot_rad = 0.0_dp
469 : END IF
470 :
471 : ! Initialize the local data structures
472 67617 : ALLOCATE (atom2d(nkind))
473 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
474 22539 : molecule_set, molecule_only, particle_set=particle_set)
475 :
476 68095 : DO ikind = 1, nkind
477 :
478 45556 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
479 :
480 45556 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
481 45556 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX")
482 45556 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis_set, basis_type="AUX_FIT")
483 :
484 : CALL get_qs_kind(qs_kind_set(ikind), &
485 : paw_proj_set=paw_proj, &
486 : paw_atom=paw_atom, &
487 : all_potential=all_potential, &
488 : gth_potential=gth_potential, &
489 45556 : sgp_potential=sgp_potential)
490 :
491 45556 : IF (dftb) THEN
492 : ! Set the interaction radius for the neighbor lists (DFTB case)
493 : ! This includes all interactions (orbitals and short range pair potential) except vdW
494 5838 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom)
495 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom, &
496 : cutoff=orb_radius(ikind), &
497 5838 : defined=orb_present(ikind))
498 : ELSE
499 39718 : IF (ASSOCIATED(orb_basis_set)) THEN
500 39716 : orb_present(ikind) = .TRUE.
501 39716 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, kind_radius=orb_radius(ikind))
502 : ELSE
503 2 : orb_present(ikind) = .FALSE.
504 : END IF
505 : END IF
506 :
507 45556 : IF (ASSOCIATED(aux_basis_set)) THEN
508 0 : aux_present(ikind) = .TRUE.
509 : ELSE
510 45556 : aux_present(ikind) = .FALSE.
511 : END IF
512 :
513 45556 : IF (ASSOCIATED(aux_fit_basis_set)) THEN
514 1482 : aux_fit_present(ikind) = .TRUE.
515 1482 : CALL get_gto_basis_set(gto_basis_set=aux_fit_basis_set, kind_radius=aux_fit_radius(ikind))
516 : ELSE
517 44074 : aux_fit_present(ikind) = .FALSE.
518 : END IF
519 :
520 : ! core overlap
521 : CALL get_qs_kind(qs_kind_set(ikind), &
522 : alpha_core_charge=calpha(ikind), &
523 : core_charge_radius=core_radius(ikind), &
524 45556 : zeff=zeff(ikind))
525 45556 : IF (zeff(ikind) /= 0._dp .AND. calpha(ikind) /= 0._dp) THEN
526 45386 : core_present(ikind) = .TRUE.
527 : ELSE
528 170 : core_present(ikind) = .FALSE.
529 : END IF
530 :
531 : ! Pseudopotentials
532 45556 : IF (gth_potential_present .OR. sgp_potential_present) THEN
533 17106 : IF (ASSOCIATED(gth_potential)) THEN
534 : CALL get_potential(potential=gth_potential, &
535 : ppl_present=ppl_present(ikind), &
536 : ppl_radius=ppl_radius(ikind), &
537 : ppnl_present=ppnl_present(ikind), &
538 16870 : ppnl_radius=ppnl_radius(ikind))
539 236 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
540 : CALL get_potential(potential=sgp_potential, &
541 : ppl_present=ppl_present(ikind), &
542 : ppl_radius=ppl_radius(ikind), &
543 : ppnl_present=ppnl_present(ikind), &
544 38 : ppnl_radius=ppnl_radius(ikind))
545 : ELSE
546 198 : ppl_present(ikind) = .FALSE.
547 198 : ppnl_present(ikind) = .FALSE.
548 : END IF
549 : END IF
550 :
551 : ! GAPW
552 45556 : IF (paw_atom_present) THEN
553 3078 : IF (paw_atom) THEN
554 2914 : oce_present(ikind) = .TRUE.
555 2914 : CALL get_paw_proj_set(paw_proj_set=paw_proj, rcprj=oce_radius(ikind))
556 : ELSE
557 164 : oce_present(ikind) = .FALSE.
558 : END IF
559 : END IF
560 :
561 : ! Check the presence of an all electron potential or ERFC potential
562 113651 : IF (all_potential_present .OR. sgp_potential_present) THEN
563 28550 : all_present(ikind) = .FALSE.
564 28550 : all_pot_rad(ikind) = 0.0_dp
565 28550 : IF (ASSOCIATED(all_potential)) THEN
566 28478 : all_present(ikind) = .TRUE.
567 28478 : CALL get_potential(potential=all_potential, core_charge_radius=all_pot_rad(ikind))
568 72 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
569 38 : IF (sgp_potential%ecp_local) THEN
570 26 : all_present(ikind) = .TRUE.
571 26 : CALL get_potential(potential=sgp_potential, core_charge_radius=all_pot_rad(ikind))
572 : END IF
573 : END IF
574 : END IF
575 :
576 : END DO
577 :
578 : ! Build the orbital-orbital overlap neighbor lists
579 22539 : IF (pdist < 0.0_dp) THEN
580 : pdist = MAX(plane_distance(1, 0, 0, cell), &
581 : plane_distance(0, 1, 0, cell), &
582 4 : plane_distance(0, 0, 1, cell))
583 : END IF
584 22539 : CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius, pdist)
585 : CALL build_neighbor_lists(sab_orb, particle_set, atom2d, cell, pair_radius, &
586 22539 : mic=mic, subcells=subcells, molecular=molecule_only, nlname="sab_orb")
587 22539 : CALL set_ks_env(ks_env=ks_env, sab_orb=sab_orb)
588 : CALL write_neighbor_lists(sab_orb, particle_set, cell, para_env, neighbor_list_section, &
589 22539 : "/SAB_ORB", "sab_orb", "ORBITAL ORBITAL")
590 :
591 : ! Build orbital-orbital list containing all the pairs, to be used with
592 : ! non-symmetric operators. Beware: the cutoff of the orbital-orbital overlap
593 : ! might not be optimal. It should be verified for each operator.
594 22539 : IF (.NOT. (nddo .OR. dftb .OR. xtb)) THEN
595 : CALL build_neighbor_lists(sab_all, particle_set, atom2d, cell, pair_radius, &
596 10259 : mic=mic, symmetric=.FALSE., subcells=subcells, molecular=molecule_only, nlname="sab_all")
597 10259 : CALL set_ks_env(ks_env=ks_env, sab_all=sab_all)
598 : END IF
599 :
600 : ! Build the core-core overlap neighbor lists
601 22539 : IF (.NOT. (nddo .OR. dftb .OR. xtb)) THEN
602 10259 : CALL pair_radius_setup(core_present, core_present, core_radius, core_radius, pair_radius)
603 : CALL build_neighbor_lists(sab_core, particle_set, atom2d, cell, pair_radius, subcells=subcells, &
604 10259 : operator_type="PP", nlname="sab_core")
605 10259 : CALL set_ks_env(ks_env=ks_env, sab_core=sab_core)
606 : CALL write_neighbor_lists(sab_core, particle_set, cell, para_env, neighbor_list_section, &
607 10259 : "/SAB_CORE", "sab_core", "CORE CORE")
608 : END IF
609 :
610 22539 : IF (dokp) THEN
611 : ! We try to guess an integration radius for K-points
612 : ! For non-HFX calculations we use the overlap list
613 : ! For HFX we use the interaction radius of kinds (ORB or ADMM basis)
614 : ! plus a range for the operator
615 764 : IF (do_hfx) THEN
616 :
617 : !case study on the HFX potential: TC, SR or Overlap?
618 70 : CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%POTENTIAL_TYPE", i_val=hfx_pot)
619 :
620 26 : SELECT CASE (hfx_pot)
621 : CASE (do_potential_id)
622 26 : roperator = 0.0_dp
623 : CASE (do_potential_truncated)
624 44 : CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", r_val=roperator)
625 : CASE (do_potential_short)
626 0 : CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%OMEGA", r_val=omega)
627 0 : CALL section_vals_val_get(hfx_sections, "SCREENING%EPS_SCHWARZ", r_val=eps_schwarz)
628 0 : CALL erfc_cutoff(eps_schwarz, omega, roperator)
629 : CASE DEFAULT
630 70 : CPABORT("HFX potential not available for K-points (NYI)")
631 : END SELECT
632 :
633 70 : IF (dft_control%do_admm) THEN
634 : CALL pair_radius_setup(aux_fit_present, aux_fit_present, aux_fit_radius, aux_fit_radius, &
635 38 : pair_radius)
636 :
637 : !We cannot accept a pair radius smaller than the ORB overlap, for sanity reasons
638 152 : ALLOCATE (pair_radius_lb(nkind, nkind))
639 38 : CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius_lb)
640 98 : DO jkind = 1, nkind
641 202 : DO ikind = 1, nkind
642 104 : IF (pair_radius(ikind, jkind) + cutoff_screen_factor*roperator .LE. pair_radius_lb(ikind, jkind)) &
643 128 : pair_radius(ikind, jkind) = pair_radius_lb(ikind, jkind) - roperator
644 : END DO
645 : END DO
646 : ELSE
647 32 : CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
648 : END IF
649 362 : pair_radius = pair_radius + cutoff_screen_factor*roperator
650 : ELSE
651 694 : CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
652 : END IF
653 : CALL build_neighbor_lists(sab_kp, particle_set, atom2d, cell, pair_radius, &
654 764 : subcells=subcells, nlname="sab_kp")
655 764 : CALL set_ks_env(ks_env=ks_env, sab_kp=sab_kp)
656 :
657 764 : IF (do_hfx) THEN
658 : CALL build_neighbor_lists(sab_kp_nosym, particle_set, atom2d, cell, pair_radius, &
659 70 : subcells=subcells, nlname="sab_kp_nosym", symmetric=.FALSE.)
660 70 : CALL set_ks_env(ks_env=ks_env, sab_kp_nosym=sab_kp_nosym)
661 : END IF
662 : END IF
663 :
664 : ! Build orbital GTH-PPL operator overlap list
665 22539 : IF (gth_potential_present .OR. sgp_potential_present) THEN
666 9677 : IF (ANY(ppl_present)) THEN
667 9585 : CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius)
668 : CALL build_neighbor_lists(sac_ppl, particle_set, atom2d, cell, pair_radius, &
669 9585 : subcells=subcells, operator_type="ABC", nlname="sac_ppl")
670 9585 : CALL set_ks_env(ks_env=ks_env, sac_ppl=sac_ppl)
671 : CALL write_neighbor_lists(sac_ppl, particle_set, cell, para_env, neighbor_list_section, &
672 9585 : "/SAC_PPL", "sac_ppl", "ORBITAL GTH-PPL")
673 9585 : IF (lrigpw) THEN
674 58 : IF (qs_env%lri_env%ppl_ri) THEN
675 : CALL build_neighbor_lists(sac_lri, particle_set, atom2d, cell, pair_radius, &
676 2 : subcells=subcells, symmetric=.FALSE., operator_type="PP", nlname="sac_lri")
677 2 : CALL set_ks_env(ks_env=ks_env, sac_lri=sac_lri)
678 : END IF
679 : END IF
680 : END IF
681 :
682 12769 : IF (ANY(ppnl_present)) THEN
683 7673 : CALL pair_radius_setup(orb_present, ppnl_present, orb_radius, ppnl_radius, pair_radius)
684 : CALL build_neighbor_lists(sap_ppnl, particle_set, atom2d, cell, pair_radius, &
685 7673 : subcells=subcells, operator_type="ABBA", nlname="sap_ppnl")
686 7673 : CALL set_ks_env(ks_env=ks_env, sap_ppnl=sap_ppnl)
687 : CALL write_neighbor_lists(sap_ppnl, particle_set, cell, para_env, neighbor_list_section, &
688 7673 : "/SAP_PPNL", "sap_ppnl", "ORBITAL GTH-PPNL")
689 : END IF
690 : END IF
691 :
692 22539 : IF (paw_atom_present) THEN
693 : ! Build orbital-GAPW projector overlap list
694 1628 : IF (ANY(oce_present)) THEN
695 1560 : CALL pair_radius_setup(orb_present, oce_present, orb_radius, oce_radius, pair_radius)
696 : CALL build_neighbor_lists(sap_oce, particle_set, atom2d, cell, pair_radius, &
697 1560 : subcells=subcells, operator_type="ABBA", nlname="sap_oce")
698 1560 : CALL set_ks_env(ks_env=ks_env, sap_oce=sap_oce)
699 : CALL write_neighbor_lists(sap_oce, particle_set, cell, para_env, neighbor_list_section, &
700 1560 : "/SAP_OCE", "sap_oce", "ORBITAL(A) PAW-PRJ")
701 : END IF
702 : END IF
703 :
704 : ! Build orbital-ERFC potential list
705 22539 : IF (.NOT. (nddo .OR. dftb .OR. xtb)) THEN
706 10259 : IF (all_potential_present .OR. sgp_potential_present) THEN
707 720 : CALL pair_radius_setup(orb_present, all_present, orb_radius, all_pot_rad, pair_radius)
708 : CALL build_neighbor_lists(sac_ae, particle_set, atom2d, cell, pair_radius, &
709 720 : subcells=subcells, operator_type="ABC", nlname="sac_ae")
710 720 : CALL set_ks_env(ks_env=ks_env, sac_ae=sac_ae)
711 : CALL write_neighbor_lists(sac_ae, particle_set, cell, para_env, neighbor_list_section, &
712 720 : "/SAC_AE", "sac_ae", "ORBITAL ERFC POTENTIAL")
713 : END IF
714 : END IF
715 :
716 22539 : IF (nddo) THEN
717 : ! Semi-empirical neighbor lists
718 18736 : default_present = .TRUE.
719 18736 : c_radius = dft_control%qs_control%se_control%cutoff_cou
720 : ! Build the neighbor lists for the Hartree terms
721 5854 : CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
722 5854 : IF (dft_control%qs_control%se_control%do_ewald_gks) THEN
723 : ! Use MIC for the periodic code of GKS
724 : CALL build_neighbor_lists(sab_se, particle_set, atom2d, cell, pair_radius, mic=mic, &
725 2 : subcells=subcells, nlname="sab_se")
726 : ELSE
727 : CALL build_neighbor_lists(sab_se, particle_set, atom2d, cell, pair_radius, &
728 5852 : subcells=subcells, nlname="sab_se")
729 : END IF
730 5854 : CALL set_ks_env(ks_env=ks_env, sab_se=sab_se)
731 : CALL write_neighbor_lists(sab_se, particle_set, cell, para_env, neighbor_list_section, &
732 5854 : "/SAB_SE", "sab_se", "HARTREE INTERACTIONS")
733 :
734 : ! If requested build the SE long-range correction neighbor list
735 5854 : IF ((dft_control%qs_control%se_control%do_ewald) .AND. &
736 : (dft_control%qs_control%se_control%integral_screening /= do_se_IS_slater)) THEN
737 328 : c_radius = dft_control%qs_control%se_control%cutoff_lrc
738 140 : CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
739 : CALL build_neighbor_lists(sab_lrc, particle_set, atom2d, cell, pair_radius, &
740 140 : subcells=subcells, nlname="sab_lrc")
741 140 : CALL set_ks_env(ks_env=ks_env, sab_lrc=sab_lrc)
742 : CALL write_neighbor_lists(sab_lrc, particle_set, cell, para_env, neighbor_list_section, &
743 140 : "/SAB_LRC", "sab_lrc", "SE LONG-RANGE CORRECTION")
744 : END IF
745 : END IF
746 :
747 22539 : IF (dftb) THEN
748 : ! Build the neighbor lists for the DFTB Ewald methods
749 2884 : IF (dft_control%qs_control%dftb_control%do_ewald) THEN
750 1074 : CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env)
751 1074 : CALL ewald_env_get(ewald_env, rcut=rcut)
752 3150 : c_radius = rcut
753 1074 : CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius)
754 : CALL build_neighbor_lists(sab_tbe, particle_set, atom2d, cell, pair_radius, mic=mic, &
755 1074 : subcells=subcells, nlname="sab_tbe")
756 1074 : CALL set_ks_env(ks_env=ks_env, sab_tbe=sab_tbe)
757 : END IF
758 :
759 : ! Build the neighbor lists for the DFTB vdW pair potential
760 2884 : IF (dft_control%qs_control%dftb_control%dispersion) THEN
761 1016 : IF (dft_control%qs_control%dftb_control%dispersion_type == dispersion_uff) THEN
762 2754 : DO ikind = 1, nkind
763 1828 : CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom)
764 2754 : CALL get_dftb_atom_param(dftb_parameter=dftb_atom, rcdisp=c_radius(ikind))
765 : END DO
766 2754 : default_present = .TRUE.
767 926 : CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
768 : CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, &
769 926 : subcells=subcells, nlname="sab_vdw")
770 926 : CALL set_ks_env(ks_env=ks_env, sab_vdw=sab_vdw)
771 : END IF
772 : END IF
773 : END IF
774 :
775 22539 : IF (xtb) THEN
776 : ! Build the neighbor lists for the xTB Ewald method
777 3542 : IF (dft_control%qs_control%xtb_control%do_ewald) THEN
778 1012 : CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env)
779 1012 : CALL ewald_env_get(ewald_env, rcut=rcut)
780 3378 : c_radius = rcut
781 1012 : CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius)
782 : CALL build_neighbor_lists(sab_tbe, particle_set, atom2d, cell, pair_radius, mic=mic, &
783 1012 : subcells=subcells, nlname="sab_tbe")
784 1012 : CALL set_ks_env(ks_env=ks_env, sab_tbe=sab_tbe)
785 : END IF
786 : ! SR part of Coulomb interaction
787 12006 : DO ikind = 1, nkind
788 8464 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom)
789 12006 : CALL get_xtb_atom_param(xtb_parameter=xtb_atom, rcut=c_radius(ikind))
790 : END DO
791 12006 : default_present = .TRUE.
792 3542 : CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
793 : CALL build_neighbor_lists(sab_xtbe, particle_set, atom2d, cell, pair_radius, &
794 3542 : subcells=subcells, nlname="sab_xtbe")
795 3542 : CALL set_ks_env(ks_env=ks_env, sab_xtbe=sab_xtbe)
796 : ! XB list
797 17710 : ALLOCATE (xb1_atom(nkind), xb2_atom(nkind))
798 12006 : c_radius = 0.5_dp*dft_control%qs_control%xtb_control%xb_radius
799 12006 : DO ikind = 1, nkind
800 8464 : CALL get_atomic_kind(atomic_kind_set(ikind), z=zat)
801 8464 : IF (zat == 17 .OR. zat == 35 .OR. zat == 53 .OR. zat == 85) THEN
802 36 : xb1_atom(ikind) = .TRUE.
803 : ELSE
804 8428 : xb1_atom(ikind) = .FALSE.
805 : END IF
806 20470 : IF (zat == 7 .OR. zat == 8 .OR. zat == 15 .OR. zat == 16) THEN
807 3652 : xb2_atom(ikind) = .TRUE.
808 : ELSE
809 4812 : xb2_atom(ikind) = .FALSE.
810 : END IF
811 : END DO
812 3542 : CALL pair_radius_setup(xb1_atom, xb2_atom, c_radius, c_radius, pair_radius)
813 : CALL build_neighbor_lists(sab_xb, particle_set, atom2d, cell, pair_radius, &
814 3542 : symmetric=.FALSE., subcells=subcells, operator_type="PP", nlname="sab_xb")
815 3542 : CALL set_ks_env(ks_env=ks_env, sab_xb=sab_xb)
816 : CALL write_neighbor_lists(sab_xb, particle_set, cell, para_env, neighbor_list_section, &
817 3542 : "/SAB_XB", "sab_xb", "XB bonding")
818 :
819 : ! nonbonded interactions list
820 3542 : IF (dft_control%qs_control%xtb_control%do_nonbonded) THEN
821 24 : ngp = SIZE(dft_control%qs_control%xtb_control%nonbonded%pot)
822 120 : ALLOCATE (nonbond1_atom(nkind), nonbond2_atom(nkind))
823 120 : nonbond1_atom = .FALSE.
824 120 : nonbond2_atom = .FALSE.
825 48 : DO ingp = 1, ngp
826 120 : DO ikind = 1, nkind
827 96 : rcut = SQRT(dft_control%qs_control%xtb_control%nonbonded%pot(ingp)%pot%rcutsq)
828 480 : c_radius = rcut
829 96 : CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=element_symbol)
830 96 : CALL uppercase(element_symbol)
831 120 : IF (TRIM(dft_control%qs_control%xtb_control%nonbonded%pot(ingp)%pot%at1) == TRIM(element_symbol)) THEN
832 24 : nonbond1_atom(ikind) = .TRUE.
833 120 : DO jkind = 1, nkind
834 96 : CALL get_atomic_kind(atomic_kind_set(jkind), element_symbol=element_symbol2)
835 96 : CALL uppercase(element_symbol2)
836 120 : IF (TRIM(dft_control%qs_control%xtb_control%nonbonded%pot(ingp)%pot%at2) == TRIM(element_symbol2)) THEN
837 24 : nonbond2_atom(jkind) = .TRUE.
838 : END IF
839 : END DO
840 : END IF
841 : END DO
842 24 : CALL pair_radius_setup(nonbond1_atom, nonbond2_atom, c_radius, c_radius, pair_radius)
843 : CALL build_neighbor_lists(sab_xtb_nonbond, particle_set, atom2d, cell, pair_radius, &
844 24 : symmetric=.FALSE., subcells=subcells, operator_type="PP", nlname="sab_xtb_nonbond")
845 24 : CALL set_ks_env(ks_env=ks_env, sab_xtb_nonbond=sab_xtb_nonbond)
846 : CALL write_neighbor_lists(sab_xtb_nonbond, particle_set, cell, para_env, neighbor_list_section, &
847 48 : "/SAB_XTB_NONBOND", "sab_xtb_nonbond", "XTB NONBONDED INTERACTIONS")
848 : END DO
849 : END IF
850 : END IF
851 :
852 : ! Build the neighbor lists for the vdW pair potential
853 22539 : CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
854 22539 : sab_vdw => dispersion_env%sab_vdw
855 22539 : sab_cn => dispersion_env%sab_cn
856 22539 : IF (dispersion_env%type == xc_vdw_fun_pairpot .OR. xtb) THEN
857 12750 : c_radius(:) = dispersion_env%rc_disp
858 12750 : default_present = .TRUE. !include all atoms in vdW (even without basis)
859 3816 : CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
860 : CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, &
861 3816 : subcells=subcells, operator_type="PP", nlname="sab_vdw")
862 3816 : dispersion_env%sab_vdw => sab_vdw
863 :
864 3816 : IF (xtb .OR. dispersion_env%pp_type == vdw_pairpot_dftd3 .OR. &
865 : dispersion_env%pp_type == vdw_pairpot_dftd3bj) THEN
866 : ! Build the neighbor lists for coordination numbers as needed by the DFT-D3 method
867 : ! This is also needed for the xTB Hamiltonian
868 12618 : DO ikind = 1, nkind
869 8846 : CALL get_atomic_kind(atomic_kind_set(ikind), z=zat)
870 12618 : c_radius(ikind) = 4._dp*ptable(zat)%covalent_radius*bohr
871 : END DO
872 3772 : CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
873 : CALL build_neighbor_lists(sab_cn, particle_set, atom2d, cell, pair_radius, &
874 3772 : subcells=subcells, operator_type="PP", nlname="sab_cn")
875 3772 : dispersion_env%sab_cn => sab_cn
876 : END IF
877 : END IF
878 :
879 : ! Build the neighbor lists for the gCP pair potential
880 22539 : NULLIFY (gcp_env)
881 22539 : CALL get_qs_env(qs_env=qs_env, gcp_env=gcp_env)
882 22539 : IF (ASSOCIATED(gcp_env)) THEN
883 10259 : IF (gcp_env%do_gcp) THEN
884 0 : sab_gcp => gcp_env%sab_gcp
885 0 : DO ikind = 1, nkind
886 0 : c_radius(ikind) = gcp_env%gcp_kind(ikind)%rcsto
887 : END DO
888 0 : CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius)
889 : CALL build_neighbor_lists(sab_gcp, particle_set, atom2d, cell, pair_radius, &
890 0 : subcells=subcells, operator_type="PP", nlname="sab_gcp")
891 0 : gcp_env%sab_gcp => sab_gcp
892 : ELSE
893 10259 : NULLIFY (gcp_env%sab_gcp)
894 : END IF
895 : END IF
896 :
897 22539 : IF (lrigpw .OR. lri_optbas) THEN
898 : ! set neighborlists in lri_env environment
899 64 : CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
900 64 : soo_list => qs_env%lri_env%soo_list
901 : CALL build_neighbor_lists(soo_list, particle_set, atom2d, cell, pair_radius, &
902 64 : mic=mic, molecular=molecule_only, subcells=subcells, nlname="soo_list")
903 64 : qs_env%lri_env%soo_list => soo_list
904 : CALL write_neighbor_lists(soo_list, particle_set, cell, para_env, neighbor_list_section, &
905 64 : "/SOO_LIST", "soo_list", "ORBITAL ORBITAL (RI)")
906 22475 : ELSEIF (rigpw) THEN
907 0 : ALLOCATE (ri_present(nkind), ri_radius(nkind))
908 0 : ri_present = .FALSE.
909 0 : ri_radius = 0.0_dp
910 0 : DO ikind = 1, nkind
911 0 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="RI_HXC")
912 0 : IF (ASSOCIATED(ri_basis_set)) THEN
913 0 : ri_present(ikind) = .TRUE.
914 0 : CALL get_gto_basis_set(gto_basis_set=ri_basis_set, kind_radius=ri_radius(ikind))
915 : ELSE
916 0 : ri_present(ikind) = .FALSE.
917 : END IF
918 : END DO
919 : ! set neighborlists in lri_env environment
920 0 : CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
921 0 : soo_list => qs_env%lri_env%soo_list
922 : CALL build_neighbor_lists(soo_list, particle_set, atom2d, cell, pair_radius, &
923 0 : mic=mic, molecular=molecule_only, subcells=subcells, nlname="soo_list")
924 0 : qs_env%lri_env%soo_list => soo_list
925 : !
926 0 : CALL pair_radius_setup(ri_present, ri_present, ri_radius, ri_radius, pair_radius)
927 0 : saa_list => qs_env%lri_env%saa_list
928 : CALL build_neighbor_lists(saa_list, particle_set, atom2d, cell, pair_radius, &
929 0 : mic=mic, molecular=molecule_only, subcells=subcells, nlname="saa_list")
930 0 : qs_env%lri_env%saa_list => saa_list
931 : !
932 0 : CALL pair_radius_setup(ri_present, orb_present, ri_radius, orb_radius, pair_radius)
933 0 : soa_list => qs_env%lri_env%soa_list
934 : CALL build_neighbor_lists(soa_list, particle_set, atom2d, cell, pair_radius, &
935 : mic=mic, symmetric=.FALSE., molecular=molecule_only, &
936 0 : subcells=subcells, operator_type="ABC", nlname="saa_list")
937 0 : qs_env%lri_env%soa_list => soa_list
938 : END IF
939 :
940 : ! Build the neighbor lists for the ALMO delocalization
941 22539 : IF (almo) THEN
942 360 : DO ikind = 1, nkind
943 244 : CALL get_atomic_kind(atomic_kind_set(ikind), rcov=almo_rcov, rvdw=almo_rvdw)
944 : ! multiply the radius by some hard-coded number
945 : c_radius(ikind) = MAX(almo_rcov, almo_rvdw)*bohr* &
946 360 : almo_max_cutoff_multiplier
947 : END DO
948 360 : default_present = .TRUE. !include all atoms (even without basis)
949 116 : CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
950 : CALL build_neighbor_lists(sab_almo, particle_set, atom2d, cell, pair_radius, &
951 116 : subcells=subcells, operator_type="PP", nlname="sab_almo")
952 116 : CALL set_ks_env(ks_env=ks_env, sab_almo=sab_almo)
953 : END IF
954 :
955 : ! Print particle distribution
956 22539 : print_key_path = "PRINT%DISTRIBUTION"
957 22539 : IF (BTEST(cp_print_key_should_output(logger%iter_info, force_env_section, &
958 : print_key_path), &
959 : cp_p_file)) THEN
960 : iw = cp_print_key_unit_nr(logger=logger, &
961 : basis_section=force_env_section, &
962 : print_key_path=print_key_path, &
963 138 : extension=".out")
964 138 : CALL write_neighbor_distribution(sab_orb, qs_kind_set, iw, para_env)
965 : CALL cp_print_key_finished_output(unit_nr=iw, &
966 : logger=logger, &
967 : basis_section=force_env_section, &
968 138 : print_key_path=print_key_path)
969 : END IF
970 :
971 : ! Release work storage
972 22539 : CALL atom2d_cleanup(atom2d)
973 :
974 22539 : DEALLOCATE (atom2d)
975 22539 : DEALLOCATE (orb_present, default_present, core_present)
976 22539 : DEALLOCATE (orb_radius, aux_fit_radius, c_radius, core_radius)
977 22539 : DEALLOCATE (calpha, zeff)
978 22539 : DEALLOCATE (pair_radius)
979 22539 : IF (gth_potential_present .OR. sgp_potential_present) THEN
980 9587 : DEALLOCATE (ppl_present, ppl_radius)
981 9587 : DEALLOCATE (ppnl_present, ppnl_radius)
982 : END IF
983 22539 : IF (paw_atom_present) THEN
984 1560 : DEALLOCATE (oce_present, oce_radius)
985 : END IF
986 22539 : IF (all_potential_present .OR. sgp_potential_present) THEN
987 13000 : DEALLOCATE (all_present, all_pot_rad)
988 : END IF
989 :
990 22539 : CALL timestop(handle)
991 :
992 67617 : END SUBROUTINE build_qs_neighbor_lists
993 :
994 : ! **************************************************************************************************
995 : !> \brief Build simple pair neighbor lists.
996 : !> \param ab_list ...
997 : !> \param particle_set ...
998 : !> \param atom ...
999 : !> \param cell ...
1000 : !> \param pair_radius ...
1001 : !> \param subcells ...
1002 : !> \param mic ...
1003 : !> \param symmetric ...
1004 : !> \param molecular ...
1005 : !> \param subset_of_mol ...
1006 : !> \param current_subset ...
1007 : !> \param operator_type ...
1008 : !> \param nlname ...
1009 : !> \param atomb_to_keep the list of atom indices to keep for pairs from the atom2d%b_list
1010 : !> \date 20.03.2002
1011 : !> \par History
1012 : !> - Major refactoring (25.07.2010,jhu)
1013 : !> - Added option to filter out atoms from list_b (08.2018, A. Bussy)
1014 : !> \author MK
1015 : !> \version 2.0
1016 : ! **************************************************************************************************
1017 104021 : SUBROUTINE build_neighbor_lists(ab_list, particle_set, atom, cell, pair_radius, subcells, &
1018 : mic, symmetric, molecular, subset_of_mol, current_subset, &
1019 104021 : operator_type, nlname, atomb_to_keep)
1020 :
1021 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1022 : POINTER :: ab_list
1023 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1024 : TYPE(local_atoms_type), DIMENSION(:), INTENT(IN) :: atom
1025 : TYPE(cell_type), POINTER :: cell
1026 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: pair_radius
1027 : REAL(dp), INTENT(IN) :: subcells
1028 : LOGICAL, INTENT(IN), OPTIONAL :: mic, symmetric, molecular
1029 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: subset_of_mol
1030 : INTEGER, OPTIONAL :: current_subset
1031 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: operator_type
1032 : CHARACTER(LEN=*), INTENT(IN) :: nlname
1033 : INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: atomb_to_keep
1034 :
1035 : CHARACTER(len=*), PARAMETER :: routineN = 'build_neighbor_lists'
1036 :
1037 : TYPE local_lists
1038 : INTEGER, DIMENSION(:), POINTER :: list
1039 : END TYPE local_lists
1040 :
1041 : INTEGER :: atom_a, atom_b, handle, i, iab, iatom, iatom_local, &
1042 : iatom_subcell, icell, ikind, j, jatom, jatom_local, jcell, jkind, k, &
1043 : kcell, maxat, mol_a, mol_b, nkind, otype, natom, inode, nnode, nentry
1044 : INTEGER, DIMENSION(3) :: cell_b, ncell, nsubcell, periodic
1045 104021 : INTEGER, DIMENSION(:), POINTER :: index_list
1046 : LOGICAL :: include_ab, my_mic, &
1047 : my_molecular, my_symmetric, my_sort_atomb
1048 104021 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: pres_a, pres_b
1049 : REAL(dp) :: rab2, rab2_max, rab_max, rabm, deth, subcell_scale
1050 : REAL(dp), DIMENSION(3) :: r, rab, ra, rb, sab_max, sb, &
1051 : sb_pbc, sb_min, sb_max, rab_pbc, pd, sab_max_guard
1052 104021 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nlista, nlistb
1053 104021 : TYPE(local_lists), DIMENSION(:), POINTER :: lista, listb
1054 : TYPE(neighbor_list_p_type), &
1055 104021 : ALLOCATABLE, DIMENSION(:) :: kind_a
1056 : TYPE(neighbor_list_set_type), POINTER :: neighbor_list_set
1057 : TYPE(subcell_type), DIMENSION(:, :, :), &
1058 104021 : POINTER :: subcell
1059 104021 : REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: r_pbc
1060 : TYPE(neighbor_list_iterator_p_type), &
1061 104021 : DIMENSION(:), POINTER :: nl_iterator
1062 :
1063 104021 : CALL timeset(routineN//"_"//TRIM(nlname), handle)
1064 :
1065 : ! input options
1066 104021 : my_mic = .FALSE.
1067 104021 : IF (PRESENT(mic)) my_mic = mic
1068 104021 : my_symmetric = .TRUE.
1069 104021 : IF (PRESENT(symmetric)) my_symmetric = symmetric
1070 104021 : my_molecular = .FALSE.
1071 : ! if we have a molecular NL, MIC has to be used
1072 104021 : IF (PRESENT(molecular)) my_molecular = molecular
1073 : ! check for operator types
1074 104021 : IF (PRESENT(operator_type)) THEN
1075 : SELECT CASE (operator_type)
1076 : CASE ("AB")
1077 11067 : otype = 1 ! simple overlap
1078 : CASE ("ABC")
1079 11067 : otype = 2 ! for three center operators
1080 11067 : CPASSERT(.NOT. my_molecular)
1081 11067 : my_symmetric = .FALSE.
1082 : CASE ("ABBA")
1083 9995 : otype = 3 ! for separable nonlocal operators
1084 9995 : my_symmetric = .FALSE.
1085 : CASE ("PP")
1086 21531 : otype = 4 ! simple atomic pair potential list
1087 : CASE default
1088 42593 : CPABORT("")
1089 : END SELECT
1090 : ELSE
1091 : ! default is a simple AB neighbor list
1092 : otype = 1
1093 : END IF
1094 104021 : my_sort_atomb = .FALSE.
1095 104021 : IF (PRESENT(atomb_to_keep)) THEN
1096 306 : my_sort_atomb = .TRUE.
1097 : END IF
1098 :
1099 104021 : nkind = SIZE(atom)
1100 : ! Deallocate the old neighbor list structure
1101 104021 : CALL release_neighbor_list_sets(ab_list)
1102 : ! Allocate and initialize the new neighbor list structure
1103 312063 : ALLOCATE (ab_list(nkind*nkind))
1104 551103 : DO iab = 1, SIZE(ab_list)
1105 447082 : NULLIFY (ab_list(iab)%neighbor_list_set)
1106 447082 : ab_list(iab)%nl_size = -1
1107 447082 : ab_list(iab)%nl_start = -1
1108 447082 : ab_list(iab)%nl_end = -1
1109 551103 : NULLIFY (ab_list(iab)%nlist_task)
1110 : END DO
1111 :
1112 : ! Allocate and initialize the kind availability
1113 520105 : ALLOCATE (pres_a(nkind), pres_b(nkind))
1114 306001 : DO ikind = 1, nkind
1115 227954 : pres_a(ikind) = ANY(pair_radius(ikind, :) > 0._dp)
1116 345381 : pres_b(ikind) = ANY(pair_radius(:, ikind) > 0._dp)
1117 : END DO
1118 :
1119 : ! create a copy of the pbc'ed coordinates
1120 104021 : natom = SIZE(particle_set)
1121 312063 : ALLOCATE (r_pbc(3, natom))
1122 653048 : DO i = 1, natom
1123 653048 : r_pbc(1:3, i) = pbc(particle_set(i)%r(1:3), cell)
1124 : END DO
1125 :
1126 : ! setup the local lists of atoms
1127 104021 : maxat = 0
1128 306001 : DO ikind = 1, nkind
1129 306001 : maxat = MAX(maxat, SIZE(atom(ikind)%list))
1130 : END DO
1131 312063 : ALLOCATE (index_list(maxat))
1132 461879 : DO i = 1, maxat
1133 461879 : index_list(i) = i
1134 : END DO
1135 936189 : ALLOCATE (lista(nkind), listb(nkind), nlista(nkind), nlistb(nkind))
1136 306001 : nlista = 0
1137 306001 : nlistb = 0
1138 306001 : DO ikind = 1, nkind
1139 201980 : NULLIFY (lista(ikind)%list, listb(ikind)%list)
1140 104021 : SELECT CASE (otype)
1141 : CASE (1)
1142 118444 : IF (ASSOCIATED(atom(ikind)%list_local_a_index)) THEN
1143 83111 : lista(ikind)%list => atom(ikind)%list_local_a_index
1144 83111 : nlista(ikind) = SIZE(lista(ikind)%list)
1145 : END IF
1146 118444 : IF (ASSOCIATED(atom(ikind)%list_local_b_index)) THEN
1147 118386 : listb(ikind)%list => atom(ikind)%list_local_b_index
1148 118386 : nlistb(ikind) = SIZE(listb(ikind)%list)
1149 : END IF
1150 : CASE (2)
1151 19510 : IF (ASSOCIATED(atom(ikind)%list_local_a_index)) THEN
1152 12912 : lista(ikind)%list => atom(ikind)%list_local_a_index
1153 12912 : nlista(ikind) = SIZE(lista(ikind)%list)
1154 : END IF
1155 19510 : nlistb(ikind) = SIZE(atom(ikind)%list)
1156 19510 : listb(ikind)%list => index_list
1157 : CASE (3)
1158 19066 : CALL combine_lists(lista(ikind)%list, nlista(ikind), ikind, atom)
1159 19066 : nlistb(ikind) = SIZE(atom(ikind)%list)
1160 19066 : listb(ikind)%list => index_list
1161 : CASE (4)
1162 44960 : nlista(ikind) = SIZE(atom(ikind)%list_1d)
1163 44960 : lista(ikind)%list => atom(ikind)%list_1d
1164 44960 : nlistb(ikind) = SIZE(atom(ikind)%list)
1165 44960 : listb(ikind)%list => index_list
1166 : CASE default
1167 201980 : CPABORT("")
1168 : END SELECT
1169 : END DO
1170 :
1171 : ! Determine max. number of local atoms
1172 104021 : maxat = 0
1173 306001 : DO ikind = 1, nkind
1174 306001 : maxat = MAX(maxat, nlista(ikind), nlistb(ikind))
1175 : END DO
1176 312063 : ALLOCATE (kind_a(2*maxat))
1177 :
1178 : ! Load informations about the simulation cell
1179 104021 : CALL get_cell(cell=cell, periodic=periodic, deth=deth)
1180 :
1181 : ! Loop over all atomic kind pairs
1182 306001 : DO ikind = 1, nkind
1183 201980 : IF (.NOT. pres_a(ikind)) CYCLE
1184 :
1185 721191 : DO jkind = 1, nkind
1186 423898 : IF (.NOT. pres_b(jkind)) CYCLE
1187 :
1188 407450 : iab = ikind + nkind*(jkind - 1)
1189 :
1190 : ! Calculate the square of the maximum interaction distance
1191 407450 : IF (pair_radius(ikind, jkind) <= 0._dp) CYCLE
1192 407406 : rab_max = pair_radius(ikind, jkind)
1193 407406 : IF (otype == 3) THEN
1194 : ! Calculate the square of the maximum interaction distance
1195 : ! for sac_max / ncell this must be the maximum over all kinds
1196 : ! to be correct for three center terms involving different kinds
1197 104484 : rabm = MAXVAL(pair_radius(:, jkind))
1198 : ELSE
1199 : rabm = rab_max
1200 : END IF
1201 407406 : rab2_max = rabm*rabm
1202 :
1203 407406 : pd(1) = plane_distance(1, 0, 0, cell)
1204 407406 : pd(2) = plane_distance(0, 1, 0, cell)
1205 407406 : pd(3) = plane_distance(0, 0, 1, cell)
1206 :
1207 1629624 : sab_max = rabm/pd
1208 1629624 : sab_max_guard = 15.0_dp/pd
1209 :
1210 : ! It makes sense to have fewer subcells for larger systems
1211 407406 : subcell_scale = ((125.0_dp**3)/deth)**(1.0_dp/6.0_dp)
1212 :
1213 : ! guess the number of subcells for optimal performance,
1214 : ! guard against crazy stuff triggered by very small rabm
1215 : nsubcell(:) = INT(MAX(1.0_dp, MIN(0.5_dp*subcells*subcell_scale/sab_max(:), &
1216 1629624 : 0.5_dp*subcells*subcell_scale/sab_max_guard(:))))
1217 :
1218 : ! number of image cells to be considered
1219 1629624 : ncell(:) = (INT(sab_max(:)) + 1)*periodic(:)
1220 :
1221 : CALL allocate_neighbor_list_set(neighbor_list_set=ab_list(iab)%neighbor_list_set, &
1222 407406 : symmetric=my_symmetric)
1223 407406 : neighbor_list_set => ab_list(iab)%neighbor_list_set
1224 :
1225 1034805 : DO iatom_local = 1, nlista(ikind)
1226 627399 : iatom = lista(ikind)%list(iatom_local)
1227 627399 : atom_a = atom(ikind)%list(iatom)
1228 : CALL add_neighbor_list(neighbor_list_set=neighbor_list_set, &
1229 : atom=atom_a, &
1230 1034805 : neighbor_list=kind_a(iatom_local)%neighbor_list)
1231 : END DO
1232 :
1233 407406 : CALL allocate_subcell(subcell, nsubcell)
1234 1034805 : DO iatom_local = 1, nlista(ikind)
1235 627399 : iatom = lista(ikind)%list(iatom_local)
1236 627399 : atom_a = atom(ikind)%list(iatom)
1237 2509596 : r = r_pbc(:, atom_a)
1238 627399 : CALL give_ijk_subcell(r, i, j, k, cell, nsubcell)
1239 1034805 : subcell(i, j, k)%natom = subcell(i, j, k)%natom + 1
1240 : END DO
1241 1185182 : DO k = 1, nsubcell(3)
1242 2826260 : DO j = 1, nsubcell(2)
1243 6701170 : DO i = 1, nsubcell(1)
1244 4282316 : maxat = subcell(i, j, k)%natom + subcell(i, j, k)%natom/10
1245 8981966 : ALLOCATE (subcell(i, j, k)%atom_list(maxat))
1246 5923394 : subcell(i, j, k)%natom = 0
1247 : END DO
1248 : END DO
1249 : END DO
1250 1034805 : DO iatom_local = 1, nlista(ikind)
1251 627399 : iatom = lista(ikind)%list(iatom_local)
1252 627399 : atom_a = atom(ikind)%list(iatom)
1253 2509596 : r = r_pbc(:, atom_a)
1254 627399 : CALL give_ijk_subcell(r, i, j, k, cell, nsubcell)
1255 627399 : subcell(i, j, k)%natom = subcell(i, j, k)%natom + 1
1256 1034805 : subcell(i, j, k)%atom_list(subcell(i, j, k)%natom) = iatom_local
1257 : END DO
1258 :
1259 1567472 : DO jatom_local = 1, nlistb(jkind)
1260 1160066 : jatom = listb(jkind)%list(jatom_local)
1261 1160066 : atom_b = atom(jkind)%list(jatom)
1262 1160066 : IF (my_sort_atomb .AND. .NOT. my_symmetric) THEN
1263 6270 : IF (.NOT. ANY(atomb_to_keep == atom_b)) CYCLE
1264 : END IF
1265 1157592 : IF (my_molecular) THEN
1266 3220 : mol_b = atom(jkind)%list_b_mol(jatom_local)
1267 3220 : IF (PRESENT(subset_of_mol)) THEN
1268 1340 : IF (subset_of_mol(mol_b) .NE. current_subset) CYCLE
1269 : END IF
1270 : END IF
1271 4627584 : r = r_pbc(:, atom_b)
1272 1156896 : CALL real_to_scaled(sb_pbc(:), r(:), cell)
1273 :
1274 4280658 : loop2_kcell: DO kcell = -ncell(3), ncell(3)
1275 2955612 : sb(3) = sb_pbc(3) + REAL(kcell, dp)
1276 2955612 : sb_min(3) = sb(3) - sab_max(3)
1277 2955612 : sb_max(3) = sb(3) + sab_max(3)
1278 2955612 : IF (periodic(3) /= 0) THEN
1279 2324840 : IF (sb_min(3) >= 0.5_dp) EXIT loop2_kcell
1280 2085584 : IF (sb_max(3) < -0.5_dp) CYCLE loop2_kcell
1281 : END IF
1282 2484802 : cell_b(3) = kcell
1283 :
1284 13741768 : loop2_jcell: DO jcell = -ncell(2), ncell(2)
1285 10859590 : sb(2) = sb_pbc(2) + REAL(jcell, dp)
1286 10859590 : sb_min(2) = sb(2) - sab_max(2)
1287 10859590 : sb_max(2) = sb(2) + sab_max(2)
1288 10859590 : IF (periodic(2) /= 0) THEN
1289 10230532 : IF (sb_min(2) >= 0.5_dp) EXIT loop2_jcell
1290 9467842 : IF (sb_max(2) < -0.5_dp) CYCLE loop2_jcell
1291 : END IF
1292 9326406 : cell_b(2) = jcell
1293 :
1294 71393954 : loop2_icell: DO icell = -ncell(1), ncell(1)
1295 63103504 : sb(1) = sb_pbc(1) + REAL(icell, dp)
1296 63103504 : sb_min(1) = sb(1) - sab_max(1)
1297 63103504 : sb_max(1) = sb(1) + sab_max(1)
1298 63103504 : IF (periodic(1) /= 0) THEN
1299 62423044 : IF (sb_min(1) >= 0.5_dp) EXIT loop2_icell
1300 58670732 : IF (sb_max(1) < -0.5_dp) CYCLE loop2_icell
1301 : END IF
1302 55992256 : cell_b(1) = icell
1303 :
1304 55992256 : CALL scaled_to_real(rb, sb, cell)
1305 :
1306 136196157 : loop_k: DO k = 1, nsubcell(3)
1307 229980960 : loop_j: DO j = 1, nsubcell(2)
1308 396927568 : loop_i: DO i = 1, nsubcell(1)
1309 :
1310 : ! FIXME for non-periodic systems, the whole subcell trick is skipped
1311 : ! yielding a Natom**2 pair list build.
1312 241936461 : IF (periodic(3) /= 0) THEN
1313 234548211 : IF (sb_max(3) < subcell(i, j, k)%s_min(3)) EXIT loop_k
1314 232436322 : IF (sb_min(3) >= subcell(i, j, k)%s_max(3)) CYCLE loop_k
1315 : END IF
1316 :
1317 236942912 : IF (periodic(2) /= 0) THEN
1318 229732028 : IF (sb_max(2) < subcell(i, j, k)%s_min(2)) EXIT loop_j
1319 226439553 : IF (sb_min(2) >= subcell(i, j, k)%s_max(2)) CYCLE loop_j
1320 : END IF
1321 :
1322 228184486 : IF (periodic(1) /= 0) THEN
1323 220858926 : IF (sb_max(1) < subcell(i, j, k)%s_min(1)) EXIT loop_i
1324 212798105 : IF (sb_min(1) >= subcell(i, j, k)%s_max(1)) CYCLE loop_i
1325 : END IF
1326 :
1327 196053243 : IF (subcell(i, j, k)%natom == 0) CYCLE
1328 :
1329 385062481 : DO iatom_subcell = 1, subcell(i, j, k)%natom
1330 227623126 : iatom_local = subcell(i, j, k)%atom_list(iatom_subcell)
1331 227623126 : iatom = lista(ikind)%list(iatom_local)
1332 227623126 : atom_a = atom(ikind)%list(iatom)
1333 227623126 : IF (my_molecular) THEN
1334 475991 : mol_a = atom(ikind)%list_a_mol(iatom_local)
1335 475991 : IF (mol_a /= mol_b) CYCLE
1336 : END IF
1337 227386639 : IF (my_symmetric) THEN
1338 216383276 : IF (atom_a > atom_b) THEN
1339 101745615 : include_ab = (MODULO(atom_a + atom_b, 2) /= 0)
1340 : ELSE
1341 114637661 : include_ab = (MODULO(atom_a + atom_b, 2) == 0)
1342 : END IF
1343 216383276 : IF (my_sort_atomb) THEN
1344 665572 : IF ((.NOT. ANY(atomb_to_keep == atom_b)) .AND. &
1345 : (.NOT. ANY(atomb_to_keep == atom_a))) THEN
1346 : include_ab = .FALSE.
1347 : END IF
1348 : END IF
1349 : ELSE
1350 : include_ab = .TRUE.
1351 : END IF
1352 436397837 : IF (include_ab) THEN
1353 501412644 : ra(:) = r_pbc(:, atom_a)
1354 501412644 : rab(:) = rb(:) - ra(:)
1355 125353161 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1356 125353161 : IF (rab2 < rab2_max) THEN
1357 37085384 : include_ab = .TRUE.
1358 37085384 : IF (my_mic) THEN
1359 : ! only if rab is minimum image the pair will be included
1360 : ! ideally the range of the pair list is < L/2 so
1361 : ! that this never triggers
1362 1334048 : rab_pbc(:) = pbc(rab(:), cell)
1363 5336192 : IF (SUM((rab_pbc - rab)**2) > EPSILON(1.0_dp)) THEN
1364 : include_ab = .FALSE.
1365 : END IF
1366 : END IF
1367 : IF (include_ab) THEN
1368 : CALL add_neighbor_node( &
1369 : neighbor_list=kind_a(iatom_local)%neighbor_list, &
1370 : neighbor=atom_b, &
1371 : cell=cell_b, &
1372 : r=rab, &
1373 36063883 : nkind=nkind)
1374 : END IF
1375 : END IF
1376 : END IF
1377 : END DO
1378 :
1379 : END DO loop_i
1380 : END DO loop_j
1381 : END DO loop_k
1382 :
1383 : END DO loop2_icell
1384 : END DO loop2_jcell
1385 : END DO loop2_kcell
1386 :
1387 : END DO
1388 :
1389 625878 : CALL deallocate_subcell(subcell)
1390 :
1391 : END DO
1392 : END DO
1393 :
1394 9995 : SELECT CASE (otype)
1395 : CASE (1:2, 4)
1396 : CASE (3)
1397 29061 : DO ikind = 1, nkind
1398 29061 : DEALLOCATE (lista(ikind)%list)
1399 : END DO
1400 : CASE default
1401 104021 : CPABORT("")
1402 : END SELECT
1403 104021 : DEALLOCATE (kind_a, pres_a, pres_b, lista, listb, nlista, nlistb)
1404 104021 : DEALLOCATE (index_list)
1405 104021 : DEALLOCATE (r_pbc)
1406 :
1407 104021 : nentry = 0
1408 104021 : CALL neighbor_list_iterator_create(nl_iterator, ab_list)
1409 36167904 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1410 36063883 : CALL get_iterator_info(nl_iterator, inode=inode, nnode=nnode)
1411 36063883 : IF (inode == 1) nentry = nentry + nnode
1412 : END DO
1413 104021 : CALL neighbor_list_iterator_release(nl_iterator)
1414 : !
1415 307121 : ALLOCATE (ab_list(1)%nlist_task(nentry))
1416 104021 : ab_list(1)%nl_size = nentry
1417 447082 : DO iab = 2, SIZE(ab_list)
1418 343061 : ab_list(iab)%nl_size = nentry
1419 447082 : ab_list(iab)%nlist_task => ab_list(1)%nlist_task
1420 : END DO
1421 : !
1422 104021 : nentry = 0
1423 104021 : CALL neighbor_list_iterator_create(nl_iterator, ab_list)
1424 36167904 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1425 36063883 : nentry = nentry + 1
1426 36063883 : CALL get_iterator_task(nl_iterator, ab_list(1)%nlist_task(nentry))
1427 36063883 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, nkind=nkind)
1428 36063883 : iab = (ikind - 1)*nkind + jkind
1429 36063883 : IF (ab_list(iab)%nl_start < 0) ab_list(iab)%nl_start = nentry
1430 36167904 : IF (ab_list(iab)%nl_end < 0) THEN
1431 275464 : ab_list(iab)%nl_end = nentry
1432 : ELSE
1433 35788419 : CPASSERT(ab_list(iab)%nl_end + 1 == nentry)
1434 35788419 : ab_list(iab)%nl_end = nentry
1435 : END IF
1436 : END DO
1437 104021 : CALL neighbor_list_iterator_release(nl_iterator)
1438 :
1439 104021 : CALL timestop(handle)
1440 :
1441 208042 : END SUBROUTINE build_neighbor_lists
1442 :
1443 : ! **************************************************************************************************
1444 : !> \brief Build a neighborlist
1445 : !> \param ab_list ...
1446 : !> \param basis_set_a ...
1447 : !> \param basis_set_b ...
1448 : !> \param qs_env ...
1449 : !> \param mic ...
1450 : !> \param symmetric ...
1451 : !> \param molecular ...
1452 : !> \param operator_type ...
1453 : !> \date 14.03.2016
1454 : !> \author JGH
1455 : ! **************************************************************************************************
1456 108 : SUBROUTINE setup_neighbor_list(ab_list, basis_set_a, basis_set_b, qs_env, &
1457 : mic, symmetric, molecular, operator_type)
1458 :
1459 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1460 : POINTER :: ab_list
1461 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_a
1462 : TYPE(gto_basis_set_p_type), DIMENSION(:), &
1463 : OPTIONAL, POINTER :: basis_set_b
1464 : TYPE(qs_environment_type), POINTER :: qs_env
1465 : LOGICAL, INTENT(IN), OPTIONAL :: mic, symmetric, molecular
1466 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: operator_type
1467 :
1468 : CHARACTER(LEN=4) :: otype
1469 : INTEGER :: ikind, nkind
1470 : LOGICAL :: my_mic, my_molecular, my_symmetric
1471 108 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: a_present, b_present
1472 108 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: a_radius, b_radius
1473 108 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius
1474 108 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1475 : TYPE(cell_type), POINTER :: cell
1476 : TYPE(distribution_1d_type), POINTER :: distribution_1d
1477 : TYPE(distribution_2d_type), POINTER :: distribution_2d
1478 108 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_a, basis_b
1479 : TYPE(gto_basis_set_type), POINTER :: abas, bbas
1480 108 : TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
1481 108 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1482 108 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1483 :
1484 108 : basis_a => basis_set_a
1485 108 : IF (PRESENT(basis_set_b)) THEN
1486 54 : basis_b => basis_set_b
1487 54 : my_symmetric = .FALSE.
1488 : ELSE
1489 54 : basis_b => basis_set_a
1490 54 : my_symmetric = .TRUE.
1491 : END IF
1492 108 : IF (PRESENT(symmetric)) my_symmetric = symmetric
1493 :
1494 108 : IF (PRESENT(mic)) THEN
1495 4 : my_mic = mic
1496 : ELSE
1497 104 : my_mic = .FALSE.
1498 : END IF
1499 :
1500 108 : IF (PRESENT(molecular)) THEN
1501 8 : my_molecular = molecular
1502 : ELSE
1503 100 : my_molecular = .FALSE.
1504 : END IF
1505 :
1506 108 : IF (PRESENT(operator_type)) THEN
1507 0 : otype = operator_type
1508 : ELSE
1509 : ! default is a simple AB neighbor list
1510 108 : otype = "AB"
1511 : END IF
1512 :
1513 108 : nkind = SIZE(basis_a)
1514 540 : ALLOCATE (a_present(nkind), b_present(nkind))
1515 332 : a_present = .FALSE.
1516 332 : b_present = .FALSE.
1517 540 : ALLOCATE (a_radius(nkind), b_radius(nkind))
1518 332 : a_radius = 0.0_dp
1519 332 : b_radius = 0.0_dp
1520 332 : DO ikind = 1, nkind
1521 224 : IF (ASSOCIATED(basis_a(ikind)%gto_basis_set)) THEN
1522 224 : a_present(ikind) = .TRUE.
1523 224 : abas => basis_a(ikind)%gto_basis_set
1524 224 : CALL get_gto_basis_set(gto_basis_set=abas, kind_radius=a_radius(ikind))
1525 : END IF
1526 332 : IF (ASSOCIATED(basis_b(ikind)%gto_basis_set)) THEN
1527 224 : b_present(ikind) = .TRUE.
1528 224 : bbas => basis_b(ikind)%gto_basis_set
1529 224 : CALL get_gto_basis_set(gto_basis_set=bbas, kind_radius=b_radius(ikind))
1530 : END IF
1531 : END DO
1532 :
1533 432 : ALLOCATE (pair_radius(nkind, nkind))
1534 812 : pair_radius = 0.0_dp
1535 108 : CALL pair_radius_setup(a_present, b_present, a_radius, b_radius, pair_radius)
1536 :
1537 : CALL get_qs_env(qs_env, &
1538 : atomic_kind_set=atomic_kind_set, &
1539 : cell=cell, &
1540 : distribution_2d=distribution_2d, &
1541 : local_particles=distribution_1d, &
1542 : particle_set=particle_set, &
1543 108 : molecule_set=molecule_set)
1544 :
1545 324 : ALLOCATE (atom2d(nkind))
1546 : CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
1547 108 : molecule_set, my_molecular, particle_set=particle_set)
1548 : CALL build_neighbor_lists(ab_list, particle_set, atom2d, cell, pair_radius, &
1549 : mic=my_mic, symmetric=my_symmetric, molecular=my_molecular, &
1550 108 : subcells=2.0_dp, nlname="AUX_NL")
1551 :
1552 108 : CALL atom2d_cleanup(atom2d)
1553 :
1554 108 : DEALLOCATE (a_present, b_present, a_radius, b_radius, pair_radius, atom2d)
1555 :
1556 108 : END SUBROUTINE setup_neighbor_list
1557 :
1558 : ! **************************************************************************************************
1559 : !> \brief ...
1560 : !> \param list ...
1561 : !> \param n ...
1562 : !> \param ikind ...
1563 : !> \param atom ...
1564 : ! **************************************************************************************************
1565 19066 : SUBROUTINE combine_lists(list, n, ikind, atom)
1566 : INTEGER, DIMENSION(:), POINTER :: list
1567 : INTEGER, INTENT(OUT) :: n
1568 : INTEGER, INTENT(IN) :: ikind
1569 : TYPE(local_atoms_type), DIMENSION(:), INTENT(IN) :: atom
1570 :
1571 : INTEGER :: i, ib, na, nb
1572 19066 : INTEGER, DIMENSION(:), POINTER :: lista, listb
1573 :
1574 0 : CPASSERT(.NOT. ASSOCIATED(list))
1575 :
1576 19066 : lista => atom(ikind)%list_local_a_index
1577 19066 : listb => atom(ikind)%list_local_b_index
1578 :
1579 19066 : IF (ASSOCIATED(lista)) THEN
1580 12021 : na = SIZE(lista)
1581 : ELSE
1582 : na = 0
1583 : END IF
1584 :
1585 19066 : IF (ASSOCIATED(listb)) THEN
1586 19066 : nb = SIZE(listb)
1587 : ELSE
1588 : nb = 0
1589 : END IF
1590 :
1591 57198 : ALLOCATE (list(na + nb))
1592 :
1593 19066 : n = na
1594 69243 : IF (na .GT. 0) list(1:na) = lista(1:na)
1595 19066 : IF (nb .GT. 0) THEN
1596 56613 : loopb: DO ib = 1, nb
1597 77567 : DO i = 1, na
1598 77567 : IF (listb(ib) == list(i)) CYCLE loopb
1599 : END DO
1600 18469 : n = n + 1
1601 56613 : list(n) = listb(ib)
1602 : END DO loopb
1603 : END IF
1604 19066 : END SUBROUTINE combine_lists
1605 :
1606 : ! **************************************************************************************************
1607 :
1608 : ! **************************************************************************************************
1609 : !> \brief ...
1610 : !> \param present_a ...
1611 : !> \param present_b ...
1612 : !> \param radius_a ...
1613 : !> \param radius_b ...
1614 : !> \param pair_radius ...
1615 : !> \param prmin ...
1616 : ! **************************************************************************************************
1617 92926 : SUBROUTINE pair_radius_setup(present_a, present_b, radius_a, radius_b, pair_radius, prmin)
1618 : LOGICAL, DIMENSION(:), INTENT(IN) :: present_a, present_b
1619 : REAL(dp), DIMENSION(:), INTENT(IN) :: radius_a, radius_b
1620 : REAL(dp), DIMENSION(:, :), INTENT(OUT) :: pair_radius
1621 : REAL(dp), INTENT(IN), OPTIONAL :: prmin
1622 :
1623 : INTEGER :: i, j, nkind
1624 : REAL(dp) :: rrmin
1625 :
1626 92926 : nkind = SIZE(present_a)
1627 :
1628 681390 : pair_radius = 0._dp
1629 :
1630 92926 : rrmin = 0.0_dp
1631 92926 : IF (PRESENT(prmin)) rrmin = prmin
1632 :
1633 274998 : DO i = 1, nkind
1634 182072 : IF (.NOT. present_a(i)) CYCLE
1635 649658 : DO j = 1, nkind
1636 383332 : IF (.NOT. present_b(j)) CYCLE
1637 366760 : pair_radius(i, j) = radius_a(i) + radius_b(j)
1638 565404 : pair_radius(i, j) = MAX(pair_radius(i, j), rrmin)
1639 : END DO
1640 : END DO
1641 :
1642 92926 : END SUBROUTINE pair_radius_setup
1643 :
1644 : ! **************************************************************************************************
1645 : !> \brief Print the distribution of the simple pair neighbor list.
1646 : !> \param ab ...
1647 : !> \param qs_kind_set ...
1648 : !> \param output_unit ...
1649 : !> \param para_env ...
1650 : !> \date 19.06.2003
1651 : !> \author MK
1652 : !> \version 1.0
1653 : ! **************************************************************************************************
1654 138 : SUBROUTINE write_neighbor_distribution(ab, qs_kind_set, output_unit, para_env)
1655 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1656 : POINTER :: ab
1657 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1658 : INTEGER, INTENT(in) :: output_unit
1659 : TYPE(mp_para_env_type), POINTER :: para_env
1660 :
1661 : CHARACTER(len=*), PARAMETER :: routineN = 'write_neighbor_distribution'
1662 : LOGICAL, PARAMETER :: full_output = .FALSE.
1663 :
1664 : INTEGER :: handle, ikind, inode, ipe, jkind, n, &
1665 : nkind, nnode
1666 : INTEGER(int_8) :: nblock_max, nblock_sum, nelement_max, &
1667 : nelement_sum, tmp(2)
1668 138 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nblock, nelement, nnsgf
1669 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1670 : TYPE(neighbor_list_iterator_p_type), &
1671 138 : DIMENSION(:), POINTER :: nl_iterator
1672 :
1673 138 : CALL timeset(routineN, handle)
1674 : ASSOCIATE (mype => para_env%mepos + 1, npe => para_env%num_pe)
1675 :
1676 : ! Allocate work storage
1677 690 : ALLOCATE (nblock(npe), nelement(npe))
1678 414 : nblock(:) = 0
1679 414 : nelement(:) = 0
1680 138 : nkind = SIZE(qs_kind_set)
1681 414 : ALLOCATE (nnsgf(nkind))
1682 372 : nnsgf = 1
1683 372 : DO ikind = 1, nkind
1684 234 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set)
1685 372 : IF (ASSOCIATED(orb_basis_set)) THEN
1686 182 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, nsgf=nnsgf(ikind))
1687 : END IF
1688 : END DO
1689 :
1690 138 : CALL neighbor_list_iterator_create(nl_iterator, ab)
1691 42740 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1692 42602 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, nnode=nnode)
1693 42740 : IF (inode == 1) THEN
1694 1065 : n = nnsgf(ikind)*nnsgf(jkind)
1695 1065 : nblock(mype) = nblock(mype) + nnode
1696 1065 : nelement(mype) = nelement(mype) + n*nnode
1697 : END IF
1698 : END DO
1699 138 : CALL neighbor_list_iterator_release(nl_iterator)
1700 :
1701 : IF (full_output) THEN
1702 : ! XXXXXXXX should gather/scatter this on ionode
1703 : CALL para_env%sum(nblock)
1704 : CALL para_env%sum(nelement)
1705 :
1706 : nblock_sum = SUM(INT(nblock, KIND=int_8))
1707 : nelement_sum = SUM(INT(nelement, KIND=int_8))
1708 : ELSE
1709 138 : nblock_sum = nblock(mype)
1710 : nblock_max = nblock(mype)
1711 138 : nelement_sum = nelement(mype)
1712 : nelement_max = nelement(mype)
1713 414 : tmp = (/nblock_sum, nelement_sum/)
1714 138 : CALL para_env%sum(tmp)
1715 138 : nblock_sum = tmp(1); nelement_sum = tmp(2)
1716 414 : tmp = (/nblock_max, nelement_max/)
1717 138 : CALL para_env%max(tmp)
1718 138 : nblock_max = tmp(1); nelement_max = tmp(2)
1719 : END IF
1720 :
1721 276 : IF (output_unit > 0) THEN
1722 : IF (full_output) THEN
1723 : WRITE (UNIT=output_unit, &
1724 : FMT="(/,/,T2,A,/,/,T3,A,/,/,(T4,I6,T27,I10,T55,I10))") &
1725 : "DISTRIBUTION OF THE NEIGHBOR LISTS", &
1726 : "Process Number of particle pairs Number of matrix elements", &
1727 : (ipe - 1, nblock(ipe), nelement(ipe), ipe=1, npe)
1728 : WRITE (UNIT=output_unit, FMT="(/,T7,A3,T27,I10,T55,I10)") &
1729 : "Sum", SUM(nblock), SUM(nelement)
1730 : ELSE
1731 69 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") "DISTRIBUTION OF THE NEIGHBOR LISTS"
1732 69 : WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Total number of particle pairs:", nblock_sum
1733 69 : WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Total number of matrix elements:", nelement_sum
1734 69 : WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Average number of particle pairs:", (nblock_sum + npe - 1)/npe
1735 69 : WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Maximum number of particle pairs:", nblock_max
1736 69 : WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Average number of matrix element:", (nelement_sum + npe - 1)/npe
1737 69 : WRITE (UNIT=output_unit, FMT="(T15,A,T68,I13)") "Maximum number of matrix elements:", nelement_max
1738 : END IF
1739 : END IF
1740 : END ASSOCIATE
1741 :
1742 : ! Release work storage
1743 :
1744 138 : DEALLOCATE (nblock, nelement, nnsgf)
1745 :
1746 138 : CALL timestop(handle)
1747 :
1748 138 : END SUBROUTINE write_neighbor_distribution
1749 :
1750 : ! **************************************************************************************************
1751 : !> \brief Write a set of neighbor lists to the output unit.
1752 : !> \param ab ...
1753 : !> \param particle_set ...
1754 : !> \param cell ...
1755 : !> \param para_env ...
1756 : !> \param neighbor_list_section ...
1757 : !> \param nl_type ...
1758 : !> \param middle_name ...
1759 : !> \param nlname ...
1760 : !> \date 04.03.2002
1761 : !> \par History
1762 : !> - Adapted to the new parallelized neighbor list version
1763 : !> (26.06.2003,MK)
1764 : !> \author MK
1765 : !> \version 1.0
1766 : ! **************************************************************************************************
1767 63880 : SUBROUTINE write_neighbor_lists(ab, particle_set, cell, para_env, neighbor_list_section, &
1768 : nl_type, middle_name, nlname)
1769 :
1770 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1771 : POINTER :: ab
1772 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1773 : TYPE(cell_type), POINTER :: cell
1774 : TYPE(mp_para_env_type), POINTER :: para_env
1775 : TYPE(section_vals_type), POINTER :: neighbor_list_section
1776 : CHARACTER(LEN=*), INTENT(IN) :: nl_type, middle_name, nlname
1777 :
1778 : CHARACTER(LEN=default_string_length) :: string, unit_str
1779 : INTEGER :: iatom, inode, iw, jatom, nneighbor, nnode
1780 : INTEGER, DIMENSION(3) :: cell_b
1781 : REAL(dp) :: dab, unit_conv
1782 : REAL(dp), DIMENSION(3) :: ra, rab, rb
1783 : TYPE(cp_logger_type), POINTER :: logger
1784 : TYPE(neighbor_list_iterator_p_type), &
1785 63880 : DIMENSION(:), POINTER :: nl_iterator
1786 :
1787 63880 : NULLIFY (logger)
1788 63880 : logger => cp_get_default_logger()
1789 63880 : IF (BTEST(cp_print_key_should_output(logger%iter_info, neighbor_list_section, &
1790 : TRIM(nl_type)), &
1791 : cp_p_file)) THEN
1792 : iw = cp_print_key_unit_nr(logger=logger, &
1793 : basis_section=neighbor_list_section, &
1794 : print_key_path=TRIM(nl_type), &
1795 : extension=".out", &
1796 : middle_name=TRIM(middle_name), &
1797 : local=.TRUE., &
1798 : log_filename=.FALSE., &
1799 4 : file_position="REWIND")
1800 : ASSOCIATE (mype => para_env%mepos)
1801 4 : CALL section_vals_val_get(neighbor_list_section, "UNIT", c_val=unit_str)
1802 4 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
1803 :
1804 : ! Print headline
1805 4 : string = ""
1806 : WRITE (UNIT=string, FMT="(A,I5,A)") &
1807 4 : TRIM(nlname)//" IN "//TRIM(unit_str)//" (PROCESS", mype, ")"
1808 4 : CALL compress(string)
1809 4 : IF (iw > 0) WRITE (UNIT=iw, FMT="(/,/,T2,A)") TRIM(string)
1810 :
1811 4 : nneighbor = 0
1812 :
1813 4 : CALL neighbor_list_iterator_create(nl_iterator, ab)
1814 16 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1815 : CALL get_iterator_info(nl_iterator, inode=inode, nnode=nnode, &
1816 12 : iatom=iatom, jatom=jatom, cell=cell_b, r=rab)
1817 12 : nneighbor = nneighbor + 1
1818 12 : ra(:) = pbc(particle_set(iatom)%r, cell)
1819 48 : rb(:) = ra(:) + rab(:)
1820 12 : dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1821 16 : IF (iw > 0) THEN
1822 12 : IF (inode == 1) THEN
1823 : WRITE (UNIT=iw, FMT="(/,T2,I5,3X,I6,3X,3F12.6)") &
1824 40 : iatom, nnode, ra(1:3)*unit_conv
1825 : END IF
1826 : WRITE (UNIT=iw, FMT="(T10,I6,3X,3I4,3F12.6,2X,F12.6)") &
1827 60 : jatom, cell_b(1:3), rb(1:3)*unit_conv, dab*unit_conv
1828 : END IF
1829 : END DO
1830 4 : CALL neighbor_list_iterator_release(nl_iterator)
1831 :
1832 4 : string = ""
1833 : WRITE (UNIT=string, FMT="(A,I12,A,I12)") &
1834 4 : "Total number of neighbor interactions for process", mype, ":", &
1835 8 : nneighbor
1836 4 : CALL compress(string)
1837 4 : IF (iw > 0) WRITE (UNIT=iw, FMT="(/,T2,A)") TRIM(string)
1838 : CALL cp_print_key_finished_output(unit_nr=iw, &
1839 : logger=logger, &
1840 : basis_section=neighbor_list_section, &
1841 : print_key_path=TRIM(nl_type), &
1842 8 : local=.TRUE.)
1843 : END ASSOCIATE
1844 : END IF
1845 :
1846 63880 : END SUBROUTINE write_neighbor_lists
1847 :
1848 0 : END MODULE qs_neighbor_lists
|
