-
Notifications
You must be signed in to change notification settings - Fork 1
/
qs_neighbor_lists.F
1856 lines (1665 loc) · 86.5 KB
/
qs_neighbor_lists.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Generate the atomic neighbor lists.
!> \par History
!> - List rebuild for sab_orb neighbor list (10.09.2002,MK)
!> - List rebuild for all lists (25.09.2002,MK)
!> - Row-wise parallelized version (16.06.2003,MK)
!> - Row- and column-wise parallelized version (19.07.2003,MK)
!> - bug fix for non-periodic case (23.02.06,MK)
!> - major refactoring (25.07.10,jhu)
!> \author Matthias Krack (08.10.1999,26.03.2002,16.06.2003)
! **************************************************************************************************
MODULE qs_neighbor_lists
USE almo_scf_types, ONLY: almo_max_cutoff_multiplier
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind,&
get_atomic_kind_set
USE basis_set_types, ONLY: get_gto_basis_set,&
gto_basis_set_p_type,&
gto_basis_set_type
USE cell_types, ONLY: cell_type,&
get_cell,&
pbc,&
plane_distance,&
real_to_scaled,&
scaled_to_real
USE cp_control_types, ONLY: dft_control_type
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_type
USE cp_output_handling, ONLY: cp_p_file,&
cp_print_key_finished_output,&
cp_print_key_should_output,&
cp_print_key_unit_nr
USE cp_units, ONLY: cp_unit_from_cp2k
USE distribution_1d_types, ONLY: distribution_1d_type
USE distribution_2d_types, ONLY: distribution_2d_type
USE ewald_environment_types, ONLY: ewald_env_get,&
ewald_environment_type
USE external_potential_types, ONLY: all_potential_type,&
get_potential,&
gth_potential_type,&
sgp_potential_type
USE input_constants, ONLY: &
dispersion_uff, do_method_lrigpw, do_method_rigpw, do_potential_id, do_potential_short, &
do_potential_truncated, do_se_IS_slater, vdw_pairpot_dftd4, xc_vdw_fun_pairpot
USE input_section_types, ONLY: section_vals_get,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE kinds, ONLY: default_string_length,&
dp,&
int_8
USE kpoint_types, ONLY: kpoint_type
USE libint_2c_3c, ONLY: cutoff_screen_factor
USE mathlib, ONLY: erfc_cutoff
USE message_passing, ONLY: mp_para_env_type
USE molecule_types, ONLY: molecule_type
USE particle_types, ONLY: particle_type
USE paw_proj_set_types, ONLY: get_paw_proj_set,&
paw_proj_set_type
USE periodic_table, ONLY: ptable
USE physcon, ONLY: bohr
USE qs_dftb_types, ONLY: qs_dftb_atom_type
USE qs_dftb_utils, ONLY: get_dftb_atom_param
USE qs_dispersion_types, ONLY: qs_dispersion_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_gcp_types, ONLY: qs_gcp_type
USE qs_kind_types, ONLY: get_qs_kind,&
get_qs_kind_set,&
qs_kind_type
USE qs_ks_types, ONLY: get_ks_env,&
qs_ks_env_type,&
set_ks_env
USE qs_neighbor_list_types, ONLY: &
add_neighbor_list, add_neighbor_node, allocate_neighbor_list_set, get_iterator_info, &
get_iterator_task, neighbor_list_iterate, neighbor_list_iterator_create, &
neighbor_list_iterator_p_type, neighbor_list_iterator_release, neighbor_list_p_type, &
neighbor_list_set_p_type, neighbor_list_set_type, release_neighbor_list_sets
USE string_utilities, ONLY: compress,&
uppercase
USE subcell_types, ONLY: allocate_subcell,&
deallocate_subcell,&
give_ijk_subcell,&
subcell_type
USE util, ONLY: locate,&
sort
USE xtb_types, ONLY: get_xtb_atom_param,&
xtb_atom_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
! **************************************************************************************************
TYPE local_atoms_type
INTEGER, DIMENSION(:), POINTER :: list => NULL(), &
list_local_a_index => NULL(), &
list_local_b_index => NULL(), &
list_1d => NULL(), &
list_a_mol => NULL(), &
list_b_mol => NULL()
END TYPE local_atoms_type
! **************************************************************************************************
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_neighbor_lists'
! private counter, used to version qs neighbor lists
INTEGER, SAVE, PRIVATE :: last_qs_neighbor_list_id_nr = 0
! Public subroutines
PUBLIC :: build_qs_neighbor_lists, local_atoms_type, atom2d_cleanup, &
atom2d_build, build_neighbor_lists, pair_radius_setup, &
setup_neighbor_list, write_neighbor_lists
CONTAINS
! **************************************************************************************************
!> \brief free the internals of atom2d
!> \param atom2d ...
!> \param
! **************************************************************************************************
SUBROUTINE atom2d_cleanup(atom2d)
TYPE(local_atoms_type), DIMENSION(:) :: atom2d
CHARACTER(len=*), PARAMETER :: routineN = 'atom2d_cleanup'
INTEGER :: handle, ikind
CALL timeset(routineN, handle)
DO ikind = 1, SIZE(atom2d)
NULLIFY (atom2d(ikind)%list)
IF (ASSOCIATED(atom2d(ikind)%list_local_a_index)) THEN
DEALLOCATE (atom2d(ikind)%list_local_a_index)
END IF
IF (ASSOCIATED(atom2d(ikind)%list_local_b_index)) THEN
DEALLOCATE (atom2d(ikind)%list_local_b_index)
END IF
IF (ASSOCIATED(atom2d(ikind)%list_a_mol)) THEN
DEALLOCATE (atom2d(ikind)%list_a_mol)
END IF
IF (ASSOCIATED(atom2d(ikind)%list_b_mol)) THEN
DEALLOCATE (atom2d(ikind)%list_b_mol)
END IF
IF (ASSOCIATED(atom2d(ikind)%list_1d)) THEN
DEALLOCATE (atom2d(ikind)%list_1d)
END IF
END DO
CALL timestop(handle)
END SUBROUTINE atom2d_cleanup
! **************************************************************************************************
!> \brief Build some distribution structure of atoms, refactored from build_qs_neighbor_lists
!> \param atom2d output
!> \param distribution_1d ...
!> \param distribution_2d ...
!> \param atomic_kind_set ...
!> \param molecule_set ...
!> \param molecule_only ...
!> \param particle_set ...
!> \author JH
! **************************************************************************************************
SUBROUTINE atom2d_build(atom2d, distribution_1d, distribution_2d, &
atomic_kind_set, molecule_set, molecule_only, particle_set)
TYPE(local_atoms_type), DIMENSION(:) :: atom2d
TYPE(distribution_1d_type), POINTER :: distribution_1d
TYPE(distribution_2d_type), POINTER :: distribution_2d
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
LOGICAL :: molecule_only
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
CHARACTER(len=*), PARAMETER :: routineN = 'atom2d_build'
INTEGER :: atom_a, handle, ia, iat, iatom, &
iatom_local, ikind, imol, natom, &
natom_a, natom_local_a, natom_local_b, &
nel, nkind
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom2mol, atom_of_kind, listindex, &
listsort
INTEGER, DIMENSION(:), POINTER :: local_cols_array, local_rows_array
CALL timeset(routineN, handle)
nkind = SIZE(atomic_kind_set)
natom = SIZE(particle_set)
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
IF (molecule_only) THEN
ALLOCATE (atom2mol(natom))
DO imol = 1, SIZE(molecule_set)
DO iat = molecule_set(imol)%first_atom, molecule_set(imol)%last_atom
atom2mol(iat) = imol
END DO
END DO
END IF
DO ikind = 1, nkind
NULLIFY (atom2d(ikind)%list)
NULLIFY (atom2d(ikind)%list_local_a_index)
NULLIFY (atom2d(ikind)%list_local_b_index)
NULLIFY (atom2d(ikind)%list_1d)
NULLIFY (atom2d(ikind)%list_a_mol)
NULLIFY (atom2d(ikind)%list_b_mol)
CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
natom_a = SIZE(atom2d(ikind)%list)
natom_local_a = distribution_2d%n_local_rows(ikind)
natom_local_b = distribution_2d%n_local_cols(ikind)
local_rows_array => distribution_2d%local_rows(ikind)%array
local_cols_array => distribution_2d%local_cols(ikind)%array
nel = distribution_1d%n_el(ikind)
ALLOCATE (atom2d(ikind)%list_1d(nel))
DO iat = 1, nel
ia = distribution_1d%list(ikind)%array(iat)
atom2d(ikind)%list_1d(iat) = atom_of_kind(ia)
END DO
ALLOCATE (listsort(natom_a), listindex(natom_a))
listsort(1:natom_a) = atom2d(ikind)%list(1:natom_a)
CALL sort(listsort, natom_a, listindex)
! Block rows
IF (natom_local_a > 0) THEN
ALLOCATE (atom2d(ikind)%list_local_a_index(natom_local_a))
ALLOCATE (atom2d(ikind)%list_a_mol(natom_local_a))
atom2d(ikind)%list_a_mol(:) = 0
! Build index vector for mapping
DO iatom_local = 1, natom_local_a
atom_a = local_rows_array(iatom_local)
iatom = locate(listsort, atom_a)
atom2d(ikind)%list_local_a_index(iatom_local) = listindex(iatom)
IF (molecule_only) atom2d(ikind)%list_a_mol(iatom_local) = atom2mol(atom_a)
END DO
END IF
! Block columns
IF (natom_local_b > 0) THEN
ALLOCATE (atom2d(ikind)%list_local_b_index(natom_local_b))
ALLOCATE (atom2d(ikind)%list_b_mol(natom_local_b))
atom2d(ikind)%list_b_mol(:) = 0
! Build index vector for mapping
DO iatom_local = 1, natom_local_b
atom_a = local_cols_array(iatom_local)
iatom = locate(listsort, atom_a)
atom2d(ikind)%list_local_b_index(iatom_local) = listindex(iatom)
IF (molecule_only) atom2d(ikind)%list_b_mol(iatom_local) = atom2mol(atom_a)
END DO
END IF
DEALLOCATE (listsort, listindex)
END DO
CALL timestop(handle)
END SUBROUTINE atom2d_build
! **************************************************************************************************
!> \brief Build all the required neighbor lists for Quickstep.
!> \param qs_env ...
!> \param para_env ...
!> \param molecular ...
!> \param force_env_section ...
!> \date 28.08.2000
!> \par History
!> - Major refactoring (25.07.2010,jhu)
!> \author MK
!> \version 1.0
! **************************************************************************************************
SUBROUTINE build_qs_neighbor_lists(qs_env, para_env, molecular, force_env_section)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(mp_para_env_type), POINTER :: para_env
LOGICAL, OPTIONAL :: molecular
TYPE(section_vals_type), POINTER :: force_env_section
CHARACTER(len=*), PARAMETER :: routineN = 'build_qs_neighbor_lists'
CHARACTER(LEN=2) :: element_symbol, element_symbol2
CHARACTER(LEN=default_string_length) :: print_key_path
INTEGER :: handle, hfx_pot, ikind, ingp, iw, jkind, &
maxatom, ngp, nkind, zat
LOGICAL :: all_potential_present, almo, dftb, do_hfx, dokp, gth_potential_present, &
lri_optbas, lrigpw, mic, molecule_only, nddo, paw_atom, paw_atom_present, rigpw, &
sgp_potential_present, xtb
LOGICAL, ALLOCATABLE, DIMENSION(:) :: all_present, aux_fit_present, aux_present, &
core_present, default_present, nonbond1_atom, nonbond2_atom, oce_present, orb_present, &
ppl_present, ppnl_present, ri_present, xb1_atom, xb2_atom
REAL(dp) :: almo_rcov, almo_rvdw, eps_schwarz, &
omega, pdist, rcut, roperator, subcells
REAL(dp), ALLOCATABLE, DIMENSION(:) :: all_pot_rad, aux_fit_radius, c_radius, calpha, &
core_radius, oce_radius, orb_radius, ppl_radius, ppnl_radius, ri_radius, zeff
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius, pair_radius_lb
TYPE(all_potential_type), POINTER :: all_potential
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_logger_type), POINTER :: logger
TYPE(dft_control_type), POINTER :: dft_control
TYPE(distribution_1d_type), POINTER :: distribution_1d
TYPE(distribution_2d_type), POINTER :: distribution_2d
TYPE(ewald_environment_type), POINTER :: ewald_env
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gto_basis_set_type), POINTER :: aux_basis_set, aux_fit_basis_set, &
orb_basis_set, ri_basis_set
TYPE(kpoint_type), POINTER :: kpoints
TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: saa_list, sab_all, sab_almo, &
sab_cn, sab_core, sab_gcp, sab_kp, sab_kp_nosym, sab_lrc, sab_orb, sab_scp, sab_se, &
sab_tbe, sab_vdw, sab_xb, sab_xtb_nonbond, sab_xtb_pp, sab_xtbe, sac_ae, sac_lri, &
sac_ppl, sap_oce, sap_ppnl, soa_list, soo_list
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(paw_proj_set_type), POINTER :: paw_proj
TYPE(qs_dftb_atom_type), POINTER :: dftb_atom
TYPE(qs_dispersion_type), POINTER :: dispersion_env
TYPE(qs_gcp_type), POINTER :: gcp_env
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(section_vals_type), POINTER :: hfx_sections, neighbor_list_section
TYPE(sgp_potential_type), POINTER :: sgp_potential
TYPE(xtb_atom_type), POINTER :: xtb_atom
CALL timeset(routineN, handle)
NULLIFY (logger)
logger => cp_get_default_logger()
NULLIFY (atomic_kind_set, qs_kind_set, cell, neighbor_list_section, &
distribution_1d, distribution_2d, gth_potential, sgp_potential, orb_basis_set, &
particle_set, molecule_set, dft_control, ks_env)
NULLIFY (sab_orb)
NULLIFY (sac_ae)
NULLIFY (sac_ppl)
NULLIFY (sac_lri)
NULLIFY (sap_ppnl)
NULLIFY (sap_oce)
NULLIFY (sab_se)
NULLIFY (sab_lrc)
NULLIFY (sab_tbe)
NULLIFY (sab_xtbe)
NULLIFY (sab_core)
NULLIFY (sab_xb)
NULLIFY (sab_xtb_pp)
NULLIFY (sab_xtb_nonbond)
NULLIFY (sab_all)
NULLIFY (sab_vdw)
NULLIFY (sab_cn)
NULLIFY (soo_list)
NULLIFY (sab_scp)
NULLIFY (sab_almo)
NULLIFY (sab_kp)
NULLIFY (sab_kp_nosym)
CALL get_qs_env(qs_env, &
ks_env=ks_env, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
cell=cell, &
kpoints=kpoints, &
distribution_2d=distribution_2d, &
local_particles=distribution_1d, &
particle_set=particle_set, &
molecule_set=molecule_set, &
dft_control=dft_control)
neighbor_list_section => section_vals_get_subs_vals(force_env_section, "DFT%PRINT%NEIGHBOR_LISTS")
! This sets the id number of the qs neighbor lists, new lists, means new version
! new version implies new sparsity of the matrices
last_qs_neighbor_list_id_nr = last_qs_neighbor_list_id_nr + 1
CALL set_ks_env(ks_env=ks_env, neighbor_list_id=last_qs_neighbor_list_id_nr)
CALL get_ks_env(ks_env=ks_env, &
sab_orb=sab_orb, &
sac_ae=sac_ae, &
sac_ppl=sac_ppl, &
sac_lri=sac_lri, &
sab_vdw=sab_vdw, &
sap_ppnl=sap_ppnl, &
sap_oce=sap_oce, &
sab_se=sab_se, &
sab_lrc=sab_lrc, &
sab_tbe=sab_tbe, &
sab_xtbe=sab_xtbe, &
sab_core=sab_core, &
sab_xb=sab_xb, &
sab_xtb_pp=sab_xtb_pp, &
sab_xtb_nonbond=sab_xtb_nonbond, &
sab_scp=sab_scp, &
sab_all=sab_all, &
sab_almo=sab_almo, &
sab_kp=sab_kp, &
sab_kp_nosym=sab_kp_nosym)
dokp = (kpoints%nkp > 0)
nddo = dft_control%qs_control%semi_empirical
dftb = dft_control%qs_control%dftb
xtb = dft_control%qs_control%xtb
almo = dft_control%qs_control%do_almo_scf
lrigpw = (dft_control%qs_control%method_id == do_method_lrigpw)
rigpw = (dft_control%qs_control%method_id == do_method_rigpw)
lri_optbas = dft_control%qs_control%lri_optbas
! molecular lists
molecule_only = .FALSE.
IF (PRESENT(molecular)) molecule_only = molecular
! minimum image convention (MIC)
mic = molecule_only
IF (dokp) THEN
! no MIC for kpoints
mic = .FALSE.
ELSEIF (nddo) THEN
! enforce MIC for interaction lists in SE
mic = .TRUE.
END IF
pdist = dft_control%qs_control%pairlist_radius
hfx_sections => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
CALL section_vals_get(hfx_sections, explicit=do_hfx)
CALL get_atomic_kind_set(atomic_kind_set, maxatom=maxatom)
CALL get_qs_kind_set(qs_kind_set, paw_atom_present=paw_atom_present, &
gth_potential_present=gth_potential_present, &
sgp_potential_present=sgp_potential_present, &
all_potential_present=all_potential_present)
CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
! Allocate work storage
nkind = SIZE(atomic_kind_set)
ALLOCATE (orb_present(nkind), aux_fit_present(nkind), aux_present(nkind), &
default_present(nkind), core_present(nkind))
ALLOCATE (orb_radius(nkind), aux_fit_radius(nkind), c_radius(nkind), &
core_radius(nkind), calpha(nkind), zeff(nkind))
orb_radius(:) = 0.0_dp
aux_fit_radius(:) = 0.0_dp
c_radius(:) = 0.0_dp
core_radius(:) = 0.0_dp
calpha(:) = 0.0_dp
zeff(:) = 0.0_dp
ALLOCATE (pair_radius(nkind, nkind))
IF (gth_potential_present .OR. sgp_potential_present) THEN
ALLOCATE (ppl_present(nkind), ppl_radius(nkind))
ppl_radius = 0.0_dp
ALLOCATE (ppnl_present(nkind), ppnl_radius(nkind))
ppnl_radius = 0.0_dp
END IF
IF (paw_atom_present) THEN
ALLOCATE (oce_present(nkind), oce_radius(nkind))
oce_radius = 0.0_dp
END IF
IF (all_potential_present .OR. sgp_potential_present) THEN
ALLOCATE (all_present(nkind), all_pot_rad(nkind))
all_pot_rad = 0.0_dp
END IF
! Initialize the local data structures
ALLOCATE (atom2d(nkind))
CALL atom2d_build(atom2d, distribution_1d, distribution_2d, atomic_kind_set, &
molecule_set, molecule_only, particle_set=particle_set)
DO ikind = 1, nkind
CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom2d(ikind)%list)
CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX")
CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis_set, basis_type="AUX_FIT")
CALL get_qs_kind(qs_kind_set(ikind), &
paw_proj_set=paw_proj, &
paw_atom=paw_atom, &
all_potential=all_potential, &
gth_potential=gth_potential, &
sgp_potential=sgp_potential)
IF (dftb) THEN
! Set the interaction radius for the neighbor lists (DFTB case)
! This includes all interactions (orbitals and short range pair potential) except vdW
CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom)
CALL get_dftb_atom_param(dftb_parameter=dftb_atom, &
cutoff=orb_radius(ikind), &
defined=orb_present(ikind))
ELSE
IF (ASSOCIATED(orb_basis_set)) THEN
orb_present(ikind) = .TRUE.
CALL get_gto_basis_set(gto_basis_set=orb_basis_set, kind_radius=orb_radius(ikind))
ELSE
orb_present(ikind) = .FALSE.
END IF
END IF
IF (ASSOCIATED(aux_basis_set)) THEN
aux_present(ikind) = .TRUE.
ELSE
aux_present(ikind) = .FALSE.
END IF
IF (ASSOCIATED(aux_fit_basis_set)) THEN
aux_fit_present(ikind) = .TRUE.
CALL get_gto_basis_set(gto_basis_set=aux_fit_basis_set, kind_radius=aux_fit_radius(ikind))
ELSE
aux_fit_present(ikind) = .FALSE.
END IF
! core overlap
CALL get_qs_kind(qs_kind_set(ikind), &
alpha_core_charge=calpha(ikind), &
core_charge_radius=core_radius(ikind), &
zeff=zeff(ikind))
IF (zeff(ikind) /= 0._dp .AND. calpha(ikind) /= 0._dp) THEN
core_present(ikind) = .TRUE.
ELSE
core_present(ikind) = .FALSE.
END IF
! Pseudopotentials
IF (gth_potential_present .OR. sgp_potential_present) THEN
IF (ASSOCIATED(gth_potential)) THEN
CALL get_potential(potential=gth_potential, &
ppl_present=ppl_present(ikind), &
ppl_radius=ppl_radius(ikind), &
ppnl_present=ppnl_present(ikind), &
ppnl_radius=ppnl_radius(ikind))
ELSE IF (ASSOCIATED(sgp_potential)) THEN
CALL get_potential(potential=sgp_potential, &
ppl_present=ppl_present(ikind), &
ppl_radius=ppl_radius(ikind), &
ppnl_present=ppnl_present(ikind), &
ppnl_radius=ppnl_radius(ikind))
ELSE
ppl_present(ikind) = .FALSE.
ppnl_present(ikind) = .FALSE.
END IF
END IF
! GAPW
IF (paw_atom_present) THEN
IF (paw_atom) THEN
oce_present(ikind) = .TRUE.
CALL get_paw_proj_set(paw_proj_set=paw_proj, rcprj=oce_radius(ikind))
ELSE
oce_present(ikind) = .FALSE.
END IF
END IF
! Check the presence of an all electron potential or ERFC potential
IF (all_potential_present .OR. sgp_potential_present) THEN
all_present(ikind) = .FALSE.
all_pot_rad(ikind) = 0.0_dp
IF (ASSOCIATED(all_potential)) THEN
all_present(ikind) = .TRUE.
CALL get_potential(potential=all_potential, core_charge_radius=all_pot_rad(ikind))
ELSE IF (ASSOCIATED(sgp_potential)) THEN
IF (sgp_potential%ecp_local) THEN
all_present(ikind) = .TRUE.
CALL get_potential(potential=sgp_potential, core_charge_radius=all_pot_rad(ikind))
END IF
END IF
END IF
END DO
! Build the orbital-orbital overlap neighbor lists
IF (pdist < 0.0_dp) THEN
pdist = MAX(plane_distance(1, 0, 0, cell), &
plane_distance(0, 1, 0, cell), &
plane_distance(0, 0, 1, cell))
END IF
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius, pdist)
CALL build_neighbor_lists(sab_orb, particle_set, atom2d, cell, pair_radius, &
mic=mic, subcells=subcells, molecular=molecule_only, nlname="sab_orb")
CALL set_ks_env(ks_env=ks_env, sab_orb=sab_orb)
CALL write_neighbor_lists(sab_orb, particle_set, cell, para_env, neighbor_list_section, &
"/SAB_ORB", "sab_orb", "ORBITAL ORBITAL")
! Build orbital-orbital list containing all the pairs, to be used with
! non-symmetric operators. Beware: the cutoff of the orbital-orbital overlap
! might not be optimal. It should be verified for each operator.
IF (.NOT. (nddo .OR. dftb .OR. xtb)) THEN
CALL build_neighbor_lists(sab_all, particle_set, atom2d, cell, pair_radius, &
mic=mic, symmetric=.FALSE., subcells=subcells, molecular=molecule_only, nlname="sab_all")
CALL set_ks_env(ks_env=ks_env, sab_all=sab_all)
END IF
! Build the core-core overlap neighbor lists
IF (.NOT. (nddo .OR. dftb .OR. xtb)) THEN
CALL pair_radius_setup(core_present, core_present, core_radius, core_radius, pair_radius)
CALL build_neighbor_lists(sab_core, particle_set, atom2d, cell, pair_radius, subcells=subcells, &
operator_type="PP", nlname="sab_core")
CALL set_ks_env(ks_env=ks_env, sab_core=sab_core)
CALL write_neighbor_lists(sab_core, particle_set, cell, para_env, neighbor_list_section, &
"/SAB_CORE", "sab_core", "CORE CORE")
END IF
IF (dokp) THEN
! We try to guess an integration radius for K-points
! For non-HFX calculations we use the overlap list
! For HFX we use the interaction radius of kinds (ORB or ADMM basis)
! plus a range for the operator
IF (do_hfx) THEN
!case study on the HFX potential: TC, SR or Overlap?
CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%POTENTIAL_TYPE", i_val=hfx_pot)
SELECT CASE (hfx_pot)
CASE (do_potential_id)
roperator = 0.0_dp
CASE (do_potential_truncated)
CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%CUTOFF_RADIUS", r_val=roperator)
CASE (do_potential_short)
CALL section_vals_val_get(hfx_sections, "INTERACTION_POTENTIAL%OMEGA", r_val=omega)
CALL section_vals_val_get(hfx_sections, "SCREENING%EPS_SCHWARZ", r_val=eps_schwarz)
CALL erfc_cutoff(eps_schwarz, omega, roperator)
CASE DEFAULT
CPABORT("HFX potential not available for K-points (NYI)")
END SELECT
IF (dft_control%do_admm) THEN
CALL pair_radius_setup(aux_fit_present, aux_fit_present, aux_fit_radius, aux_fit_radius, &
pair_radius)
!We cannot accept a pair radius smaller than the ORB overlap, for sanity reasons
ALLOCATE (pair_radius_lb(nkind, nkind))
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius_lb)
DO jkind = 1, nkind
DO ikind = 1, nkind
IF (pair_radius(ikind, jkind) + cutoff_screen_factor*roperator .LE. pair_radius_lb(ikind, jkind)) &
pair_radius(ikind, jkind) = pair_radius_lb(ikind, jkind) - roperator
END DO
END DO
ELSE
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
END IF
pair_radius = pair_radius + cutoff_screen_factor*roperator
ELSE
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
END IF
CALL build_neighbor_lists(sab_kp, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, nlname="sab_kp")
CALL set_ks_env(ks_env=ks_env, sab_kp=sab_kp)
IF (do_hfx) THEN
CALL build_neighbor_lists(sab_kp_nosym, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, nlname="sab_kp_nosym", symmetric=.FALSE.)
CALL set_ks_env(ks_env=ks_env, sab_kp_nosym=sab_kp_nosym)
END IF
END IF
! Build orbital GTH-PPL operator overlap list
IF (gth_potential_present .OR. sgp_potential_present) THEN
IF (ANY(ppl_present)) THEN
CALL pair_radius_setup(orb_present, ppl_present, orb_radius, ppl_radius, pair_radius)
CALL build_neighbor_lists(sac_ppl, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABC", nlname="sac_ppl")
CALL set_ks_env(ks_env=ks_env, sac_ppl=sac_ppl)
CALL write_neighbor_lists(sac_ppl, particle_set, cell, para_env, neighbor_list_section, &
"/SAC_PPL", "sac_ppl", "ORBITAL GTH-PPL")
IF (lrigpw) THEN
IF (qs_env%lri_env%ppl_ri) THEN
CALL build_neighbor_lists(sac_lri, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, symmetric=.FALSE., operator_type="PP", nlname="sac_lri")
CALL set_ks_env(ks_env=ks_env, sac_lri=sac_lri)
END IF
END IF
END IF
IF (ANY(ppnl_present)) THEN
CALL pair_radius_setup(orb_present, ppnl_present, orb_radius, ppnl_radius, pair_radius)
CALL build_neighbor_lists(sap_ppnl, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABBA", nlname="sap_ppnl")
CALL set_ks_env(ks_env=ks_env, sap_ppnl=sap_ppnl)
CALL write_neighbor_lists(sap_ppnl, particle_set, cell, para_env, neighbor_list_section, &
"/SAP_PPNL", "sap_ppnl", "ORBITAL GTH-PPNL")
END IF
END IF
IF (paw_atom_present) THEN
! Build orbital-GAPW projector overlap list
IF (ANY(oce_present)) THEN
CALL pair_radius_setup(orb_present, oce_present, orb_radius, oce_radius, pair_radius)
CALL build_neighbor_lists(sap_oce, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABBA", nlname="sap_oce")
CALL set_ks_env(ks_env=ks_env, sap_oce=sap_oce)
CALL write_neighbor_lists(sap_oce, particle_set, cell, para_env, neighbor_list_section, &
"/SAP_OCE", "sap_oce", "ORBITAL(A) PAW-PRJ")
END IF
END IF
! Build orbital-ERFC potential list
IF (.NOT. (nddo .OR. dftb .OR. xtb)) THEN
IF (all_potential_present .OR. sgp_potential_present) THEN
CALL pair_radius_setup(orb_present, all_present, orb_radius, all_pot_rad, pair_radius)
CALL build_neighbor_lists(sac_ae, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="ABC", nlname="sac_ae")
CALL set_ks_env(ks_env=ks_env, sac_ae=sac_ae)
CALL write_neighbor_lists(sac_ae, particle_set, cell, para_env, neighbor_list_section, &
"/SAC_AE", "sac_ae", "ORBITAL ERFC POTENTIAL")
END IF
END IF
IF (nddo) THEN
! Semi-empirical neighbor lists
default_present = .TRUE.
c_radius = dft_control%qs_control%se_control%cutoff_cou
! Build the neighbor lists for the Hartree terms
CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
IF (dft_control%qs_control%se_control%do_ewald_gks) THEN
! Use MIC for the periodic code of GKS
CALL build_neighbor_lists(sab_se, particle_set, atom2d, cell, pair_radius, mic=mic, &
subcells=subcells, nlname="sab_se")
ELSE
CALL build_neighbor_lists(sab_se, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, nlname="sab_se")
END IF
CALL set_ks_env(ks_env=ks_env, sab_se=sab_se)
CALL write_neighbor_lists(sab_se, particle_set, cell, para_env, neighbor_list_section, &
"/SAB_SE", "sab_se", "HARTREE INTERACTIONS")
! If requested build the SE long-range correction neighbor list
IF ((dft_control%qs_control%se_control%do_ewald) .AND. &
(dft_control%qs_control%se_control%integral_screening /= do_se_IS_slater)) THEN
c_radius = dft_control%qs_control%se_control%cutoff_lrc
CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_lrc, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, nlname="sab_lrc")
CALL set_ks_env(ks_env=ks_env, sab_lrc=sab_lrc)
CALL write_neighbor_lists(sab_lrc, particle_set, cell, para_env, neighbor_list_section, &
"/SAB_LRC", "sab_lrc", "SE LONG-RANGE CORRECTION")
END IF
END IF
IF (dftb) THEN
! Build the neighbor lists for the DFTB Ewald methods
IF (dft_control%qs_control%dftb_control%do_ewald) THEN
CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env)
CALL ewald_env_get(ewald_env, rcut=rcut)
c_radius = rcut
CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_tbe, particle_set, atom2d, cell, pair_radius, mic=mic, &
subcells=subcells, nlname="sab_tbe")
CALL set_ks_env(ks_env=ks_env, sab_tbe=sab_tbe)
END IF
! Build the neighbor lists for the DFTB vdW pair potential
IF (dft_control%qs_control%dftb_control%dispersion) THEN
IF (dft_control%qs_control%dftb_control%dispersion_type == dispersion_uff) THEN
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), dftb_parameter=dftb_atom)
CALL get_dftb_atom_param(dftb_parameter=dftb_atom, rcdisp=c_radius(ikind))
END DO
default_present = .TRUE.
CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, nlname="sab_vdw")
CALL set_ks_env(ks_env=ks_env, sab_vdw=sab_vdw)
END IF
END IF
END IF
IF (xtb) THEN
! Build the neighbor lists for the xTB Ewald method
IF (dft_control%qs_control%xtb_control%do_ewald) THEN
CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env)
CALL ewald_env_get(ewald_env, rcut=rcut)
c_radius = rcut
CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_tbe, particle_set, atom2d, cell, pair_radius, mic=mic, &
subcells=subcells, nlname="sab_tbe")
CALL set_ks_env(ks_env=ks_env, sab_tbe=sab_tbe)
END IF
! Repulsive Potential
pair_radius(1:nkind, 1:nkind) = dft_control%qs_control%xtb_control%rcpair(1:nkind, 1:nkind)
default_present = .TRUE.
CALL build_neighbor_lists(sab_xtb_pp, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, nlname="sab_xtb_pp")
CALL set_ks_env(ks_env=ks_env, sab_xtb_pp=sab_xtb_pp)
! SR part of Coulomb interaction
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom)
CALL get_xtb_atom_param(xtb_parameter=xtb_atom, rcut=c_radius(ikind))
END DO
default_present = .TRUE.
CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_xtbe, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, nlname="sab_xtbe")
CALL set_ks_env(ks_env=ks_env, sab_xtbe=sab_xtbe)
! XB list
ALLOCATE (xb1_atom(nkind), xb2_atom(nkind))
c_radius = 0.5_dp*dft_control%qs_control%xtb_control%xb_radius
DO ikind = 1, nkind
CALL get_atomic_kind(atomic_kind_set(ikind), z=zat)
IF (zat == 17 .OR. zat == 35 .OR. zat == 53 .OR. zat == 85) THEN
xb1_atom(ikind) = .TRUE.
ELSE
xb1_atom(ikind) = .FALSE.
END IF
IF (zat == 7 .OR. zat == 8 .OR. zat == 15 .OR. zat == 16) THEN
xb2_atom(ikind) = .TRUE.
ELSE
xb2_atom(ikind) = .FALSE.
END IF
END DO
CALL pair_radius_setup(xb1_atom, xb2_atom, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_xb, particle_set, atom2d, cell, pair_radius, &
symmetric=.FALSE., subcells=subcells, operator_type="PP", nlname="sab_xb")
CALL set_ks_env(ks_env=ks_env, sab_xb=sab_xb)
CALL write_neighbor_lists(sab_xb, particle_set, cell, para_env, neighbor_list_section, &
"/SAB_XB", "sab_xb", "XB bonding")
! nonbonded interactions list
IF (dft_control%qs_control%xtb_control%do_nonbonded) THEN
ngp = SIZE(dft_control%qs_control%xtb_control%nonbonded%pot)
ALLOCATE (nonbond1_atom(nkind), nonbond2_atom(nkind))
nonbond1_atom = .FALSE.
nonbond2_atom = .FALSE.
DO ingp = 1, ngp
DO ikind = 1, nkind
rcut = SQRT(dft_control%qs_control%xtb_control%nonbonded%pot(ingp)%pot%rcutsq)
c_radius = rcut
CALL get_atomic_kind(atomic_kind_set(ikind), element_symbol=element_symbol)
CALL uppercase(element_symbol)
IF (TRIM(dft_control%qs_control%xtb_control%nonbonded%pot(ingp)%pot%at1) == TRIM(element_symbol)) THEN
nonbond1_atom(ikind) = .TRUE.
DO jkind = 1, nkind
CALL get_atomic_kind(atomic_kind_set(jkind), element_symbol=element_symbol2)
CALL uppercase(element_symbol2)
IF (TRIM(dft_control%qs_control%xtb_control%nonbonded%pot(ingp)%pot%at2) == TRIM(element_symbol2)) THEN
nonbond2_atom(jkind) = .TRUE.
END IF
END DO
END IF
END DO
CALL pair_radius_setup(nonbond1_atom, nonbond2_atom, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_xtb_nonbond, particle_set, atom2d, cell, pair_radius, &
symmetric=.FALSE., subcells=subcells, operator_type="PP", nlname="sab_xtb_nonbond")
CALL set_ks_env(ks_env=ks_env, sab_xtb_nonbond=sab_xtb_nonbond)
CALL write_neighbor_lists(sab_xtb_nonbond, particle_set, cell, para_env, neighbor_list_section, &
"/SAB_XTB_NONBOND", "sab_xtb_nonbond", "XTB NONBONDED INTERACTIONS")
END DO
END IF
END IF
! Build the neighbor lists for the vdW pair potential
CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
sab_vdw => dispersion_env%sab_vdw
sab_cn => dispersion_env%sab_cn
IF (dispersion_env%type == xc_vdw_fun_pairpot .OR. xtb) THEN
IF (dispersion_env%pp_type == vdw_pairpot_dftd4) THEN
c_radius(:) = dispersion_env%rc_d4
ELSE
c_radius(:) = dispersion_env%rc_disp
END IF
default_present = .TRUE. !include all atoms in vdW (even without basis)
CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_vdw, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="PP", nlname="sab_vdw")
dispersion_env%sab_vdw => sab_vdw
! Build the neighbor lists for coordination numbers as needed by the DFT-D3/D4 method
! This is also needed for the xTB Hamiltonian
DO ikind = 1, nkind
CALL get_atomic_kind(atomic_kind_set(ikind), z=zat)
c_radius(ikind) = 4._dp*ptable(zat)%covalent_radius*bohr
END DO
CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_cn, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="PP", nlname="sab_cn")
dispersion_env%sab_cn => sab_cn
END IF
! Build the neighbor lists for the gCP pair potential
NULLIFY (gcp_env)
CALL get_qs_env(qs_env=qs_env, gcp_env=gcp_env)
IF (ASSOCIATED(gcp_env)) THEN
IF (gcp_env%do_gcp) THEN
sab_gcp => gcp_env%sab_gcp
DO ikind = 1, nkind
c_radius(ikind) = gcp_env%gcp_kind(ikind)%rcsto
END DO
CALL pair_radius_setup(orb_present, orb_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_gcp, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="PP", nlname="sab_gcp")
gcp_env%sab_gcp => sab_gcp
ELSE
NULLIFY (gcp_env%sab_gcp)
END IF
END IF
IF (lrigpw .OR. lri_optbas) THEN
! set neighborlists in lri_env environment
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
soo_list => qs_env%lri_env%soo_list
CALL build_neighbor_lists(soo_list, particle_set, atom2d, cell, pair_radius, &
mic=mic, molecular=molecule_only, subcells=subcells, nlname="soo_list")
qs_env%lri_env%soo_list => soo_list
CALL write_neighbor_lists(soo_list, particle_set, cell, para_env, neighbor_list_section, &
"/SOO_LIST", "soo_list", "ORBITAL ORBITAL (RI)")
ELSEIF (rigpw) THEN
ALLOCATE (ri_present(nkind), ri_radius(nkind))
ri_present = .FALSE.
ri_radius = 0.0_dp
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), basis_set=ri_basis_set, basis_type="RI_HXC")
IF (ASSOCIATED(ri_basis_set)) THEN
ri_present(ikind) = .TRUE.
CALL get_gto_basis_set(gto_basis_set=ri_basis_set, kind_radius=ri_radius(ikind))
ELSE
ri_present(ikind) = .FALSE.
END IF
END DO
! set neighborlists in lri_env environment
CALL pair_radius_setup(orb_present, orb_present, orb_radius, orb_radius, pair_radius)
soo_list => qs_env%lri_env%soo_list
CALL build_neighbor_lists(soo_list, particle_set, atom2d, cell, pair_radius, &
mic=mic, molecular=molecule_only, subcells=subcells, nlname="soo_list")
qs_env%lri_env%soo_list => soo_list
!
CALL pair_radius_setup(ri_present, ri_present, ri_radius, ri_radius, pair_radius)
saa_list => qs_env%lri_env%saa_list
CALL build_neighbor_lists(saa_list, particle_set, atom2d, cell, pair_radius, &
mic=mic, molecular=molecule_only, subcells=subcells, nlname="saa_list")
qs_env%lri_env%saa_list => saa_list
!
CALL pair_radius_setup(ri_present, orb_present, ri_radius, orb_radius, pair_radius)
soa_list => qs_env%lri_env%soa_list
CALL build_neighbor_lists(soa_list, particle_set, atom2d, cell, pair_radius, &
mic=mic, symmetric=.FALSE., molecular=molecule_only, &
subcells=subcells, operator_type="ABC", nlname="saa_list")
qs_env%lri_env%soa_list => soa_list
END IF
! Build the neighbor lists for the ALMO delocalization
IF (almo) THEN
DO ikind = 1, nkind
CALL get_atomic_kind(atomic_kind_set(ikind), rcov=almo_rcov, rvdw=almo_rvdw)
! multiply the radius by some hard-coded number
c_radius(ikind) = MAX(almo_rcov, almo_rvdw)*bohr* &
almo_max_cutoff_multiplier
END DO
default_present = .TRUE. !include all atoms (even without basis)
CALL pair_radius_setup(default_present, default_present, c_radius, c_radius, pair_radius)
CALL build_neighbor_lists(sab_almo, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, operator_type="PP", nlname="sab_almo")
CALL set_ks_env(ks_env=ks_env, sab_almo=sab_almo)
END IF
! Print particle distribution
print_key_path = "PRINT%DISTRIBUTION"
IF (BTEST(cp_print_key_should_output(logger%iter_info, force_env_section, &
print_key_path), &
cp_p_file)) THEN
iw = cp_print_key_unit_nr(logger=logger, &
basis_section=force_env_section, &
print_key_path=print_key_path, &
extension=".out")
CALL write_neighbor_distribution(sab_orb, qs_kind_set, iw, para_env)
CALL cp_print_key_finished_output(unit_nr=iw, &
logger=logger, &
basis_section=force_env_section, &
print_key_path=print_key_path)
END IF
! Release work storage
CALL atom2d_cleanup(atom2d)
DEALLOCATE (atom2d)
DEALLOCATE (orb_present, default_present, core_present)
DEALLOCATE (orb_radius, aux_fit_radius, c_radius, core_radius)
DEALLOCATE (calpha, zeff)
DEALLOCATE (pair_radius)
IF (gth_potential_present .OR. sgp_potential_present) THEN
DEALLOCATE (ppl_present, ppl_radius)
DEALLOCATE (ppnl_present, ppnl_radius)
END IF
IF (paw_atom_present) THEN
DEALLOCATE (oce_present, oce_radius)
END IF
IF (all_potential_present .OR. sgp_potential_present) THEN
DEALLOCATE (all_present, all_pot_rad)
END IF
CALL timestop(handle)
END SUBROUTINE build_qs_neighbor_lists