-
Notifications
You must be signed in to change notification settings - Fork 1
/
xas_tdp_methods.F
3533 lines (2928 loc) · 163 KB
/
xas_tdp_methods.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 Methods for X-Ray absorption spectroscopy (XAS) using TDDFPT
!> \author AB (11.2017)
! **************************************************************************************************
MODULE xas_tdp_methods
USE admm_types, ONLY: admm_type
USE admm_utils, ONLY: admm_correct_for_eigenvalues,&
admm_uncorrect_for_eigenvalues
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind
USE basis_set_types, ONLY: &
allocate_sto_basis_set, create_gto_from_sto_basis, deallocate_gto_basis_set, &
deallocate_sto_basis_set, get_gto_basis_set, gto_basis_set_type, init_orb_basis_set, &
set_sto_basis_set, srules, sto_basis_set_type
USE bibliography, ONLY: Bussy2021a,&
cite_reference
USE cell_types, ONLY: cell_type,&
pbc
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: &
dbcsr_add_on_diag, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_filter, &
dbcsr_finalize, dbcsr_get_info, dbcsr_get_occupation, dbcsr_multiply, dbcsr_p_type, &
dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
dbcsr_type_symmetric
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
copy_fm_to_dbcsr,&
cp_dbcsr_sm_fm_multiply
USE cp_files, ONLY: close_file,&
open_file
USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale
USE cp_fm_diag, ONLY: cp_fm_geeig,&
cp_fm_power,&
cp_fm_syevd
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
cp_fm_struct_release,&
cp_fm_struct_type
USE cp_fm_types, ONLY: &
cp_fm_copy_general, cp_fm_create, cp_fm_get_diag, cp_fm_get_info, cp_fm_get_submatrix, &
cp_fm_read_unformatted, cp_fm_release, cp_fm_set_all, cp_fm_to_fm, cp_fm_to_fm_submat, &
cp_fm_type, cp_fm_write_unformatted
USE cp_log_handling, ONLY: cp_get_default_logger,&
cp_logger_get_default_io_unit,&
cp_logger_type,&
cp_to_string
USE cp_output_handling, ONLY: cp_p_file,&
cp_print_key_finished_output,&
cp_print_key_generate_filename,&
cp_print_key_should_output,&
cp_print_key_unit_nr,&
debug_print_level
USE input_constants, ONLY: &
do_admm_purify_cauchy_subspace, do_admm_purify_mo_diag, do_admm_purify_none, do_loc_none, &
do_potential_coulomb, do_potential_id, do_potential_short, do_potential_truncated, &
op_loc_berry, state_loc_list, tddfpt_singlet, tddfpt_spin_cons, tddfpt_spin_flip, &
tddfpt_triplet, xas_1s_type, xas_2p_type, xas_2s_type, xas_dip_len, xas_dip_vel, &
xas_not_excited, xas_tdp_by_index, xas_tdp_by_kind
USE input_cp2k_loc, ONLY: create_localize_section
USE input_section_types, ONLY: section_release,&
section_type,&
section_vals_create,&
section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get,&
section_vals_val_set
USE kinds, ONLY: default_path_length,&
default_string_length,&
dp
USE libint_wrapper, ONLY: cp_libint_static_init
USE machine, ONLY: m_flush
USE mathlib, ONLY: get_diag
USE memory_utilities, ONLY: reallocate
USE message_passing, ONLY: mp_comm_type,&
mp_para_env_type
USE parallel_gemm_api, ONLY: parallel_gemm
USE parallel_rng_types, ONLY: UNIFORM,&
rng_stream_type
USE particle_methods, ONLY: get_particle_set
USE particle_types, ONLY: particle_type
USE periodic_table, ONLY: ptable
USE physcon, ONLY: a_fine,&
angstrom,&
evolt
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_interactions, ONLY: init_interaction_radii_orb_basis
USE qs_kind_types, ONLY: get_qs_kind,&
qs_kind_type
USE qs_loc_main, ONLY: qs_loc_driver
USE qs_loc_methods, ONLY: centers_spreads_berry,&
qs_print_cubes
USE qs_loc_types, ONLY: get_qs_loc_env,&
localized_wfn_control_create,&
localized_wfn_control_type,&
qs_loc_env_create,&
qs_loc_env_release,&
qs_loc_env_type
USE qs_loc_utils, ONLY: qs_loc_control_init,&
qs_loc_env_init,&
set_loc_centers
USE qs_mo_io, ONLY: write_mo_set_low
USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
USE qs_mo_types, ONLY: allocate_mo_set,&
deallocate_mo_set,&
duplicate_mo_set,&
get_mo_set,&
init_mo_set,&
mo_set_type
USE qs_operators_ao, ONLY: p_xyz_ao,&
rRc_xyz_ao
USE qs_pdos, ONLY: calculate_projected_dos
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE qs_scf_types, ONLY: ot_method_nr
USE util, ONLY: get_limit,&
locate,&
sort_unique
USE xas_methods, ONLY: calc_stogto_overlap
USE xas_tdp_atom, ONLY: init_xas_atom_env,&
integrate_fxc_atoms,&
integrate_soc_atoms
USE xas_tdp_correction, ONLY: GW2X_shift,&
get_soc_splitting
USE xas_tdp_integrals, ONLY: compute_ri_3c_coulomb,&
compute_ri_3c_exchange,&
compute_ri_coulomb2_int,&
compute_ri_exchange2_int
USE xas_tdp_types, ONLY: &
donor_state_create, donor_state_type, free_ds_memory, free_exat_memory, &
read_xas_tdp_control, set_donor_state, set_xas_tdp_env, xas_atom_env_create, &
xas_atom_env_release, xas_atom_env_type, xas_tdp_control_create, xas_tdp_control_release, &
xas_tdp_control_type, xas_tdp_env_create, xas_tdp_env_release, xas_tdp_env_type
USE xas_tdp_utils, ONLY: include_os_soc,&
include_rcs_soc,&
setup_xas_tdp_prob,&
solve_xas_tdp_prob
USE xc_write_output, ONLY: xc_write
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xas_tdp_methods'
PUBLIC :: xas_tdp, xas_tdp_init
CONTAINS
! **************************************************************************************************
!> \brief Driver for XAS TDDFT calculations.
!> \param qs_env the inherited qs_environment
!> \author AB
!> \note Empty for now...
! **************************************************************************************************
SUBROUTINE xas_tdp(qs_env)
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'xas_tdp'
CHARACTER(default_string_length) :: rst_filename
INTEGER :: handle, n_rep, output_unit
LOGICAL :: do_restart
TYPE(section_vals_type), POINTER :: xas_tdp_section
CALL timeset(routineN, handle)
! Logger initialization and XAS TDP banner printing
NULLIFY (xas_tdp_section)
xas_tdp_section => section_vals_get_subs_vals(qs_env%input, "DFT%XAS_TDP")
output_unit = cp_logger_get_default_io_unit()
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,/,T3,A,/,T3,A,/)") &
"!===========================================================================!", &
"! XAS TDP !", &
"! Starting TDDFPT driven X-rays absorption spectroscopy calculations !", &
"!===========================================================================!"
END IF
CALL cite_reference(Bussy2021a)
! Check whether this is a restart calculation, i.e. is a restart file is provided
CALL section_vals_val_get(xas_tdp_section, "RESTART_FROM_FILE", n_rep_val=n_rep)
IF (n_rep < 1) THEN
do_restart = .FALSE.
ELSE
CALL section_vals_val_get(xas_tdp_section, "RESTART_FROM_FILE", c_val=rst_filename)
do_restart = .TRUE.
END IF
! Restart the calculation if needed
IF (do_restart) THEN
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,T3,A)") &
"# This is a RESTART calculation for PDOS and/or CUBE printing"
END IF
CALL restart_calculation(rst_filename, xas_tdp_section, qs_env)
! or run the core XAS_TDP routine if not
ELSE
CALL xas_tdp_core(xas_tdp_section, qs_env)
END IF
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,/,T3,A,/)") &
"!===========================================================================!", &
"! End of TDDFPT driven X-rays absorption spectroscopy calculations !", &
"!===========================================================================!"
END IF
CALL timestop(handle)
END SUBROUTINE xas_tdp
! **************************************************************************************************
!> \brief The core workflow of the XAS_TDP method
!> \param xas_tdp_section the input values for XAS_TDP
!> \param qs_env ...
! **************************************************************************************************
SUBROUTINE xas_tdp_core(xas_tdp_section, qs_env)
TYPE(section_vals_type), POINTER :: xas_tdp_section
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(LEN=default_string_length) :: kind_name
INTEGER :: batch_size, bo(2), current_state_index, iat, iatom, ibatch, ikind, ispin, istate, &
nbatch, nex_atom, output_unit, tmp_index
INTEGER, ALLOCATABLE, DIMENSION(:) :: batch_atoms, ex_atoms_of_kind
INTEGER, DIMENSION(:), POINTER :: atoms_of_kind
LOGICAL :: do_os, end_of_batch, unique
TYPE(admm_type), POINTER :: admm_env
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
TYPE(dft_control_type), POINTER :: dft_control
TYPE(donor_state_type), POINTER :: current_state
TYPE(gto_basis_set_type), POINTER :: tmp_basis
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(xas_atom_env_type), POINTER :: xas_atom_env
TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
NULLIFY (xas_tdp_env, xas_tdp_control, atomic_kind_set, atoms_of_kind, current_state)
NULLIFY (xas_atom_env, dft_control, matrix_ks, admm_env, qs_kind_set, tmp_basis)
! Initialization
output_unit = cp_logger_get_default_io_unit()
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,T3,A)") &
"# Create and initialize the XAS_TDP environment"
END IF
CALL get_qs_env(qs_env, dft_control=dft_control)
CALL xas_tdp_init(xas_tdp_env, xas_tdp_control, qs_env)
CALL print_info(output_unit, xas_tdp_control, qs_env)
IF (output_unit > 0) THEN
IF (xas_tdp_control%check_only) THEN
CPWARN("This is a CHECK_ONLY run for donor MOs verification")
END IF
END IF
! Localization of the core orbitals if requested (used for better identification of donor states)
IF (xas_tdp_control%do_loc) THEN
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,T3,A,/)") &
"# Localizing core orbitals for better identification"
END IF
! closed shell or ROKS => myspin=1
IF (xas_tdp_control%do_uks) THEN
DO ispin = 1, dft_control%nspins
CALL qs_loc_driver(qs_env, xas_tdp_env%qs_loc_env, &
xas_tdp_control%print_loc_subsection, myspin=ispin)
END DO
ELSE
CALL qs_loc_driver(qs_env, xas_tdp_env%qs_loc_env, &
xas_tdp_control%print_loc_subsection, myspin=1)
END IF
END IF
! Find the MO centers
CALL find_mo_centers(xas_tdp_env, xas_tdp_control, qs_env)
! Assign lowest energy orbitals to excited atoms
CALL assign_mos_to_ex_atoms(xas_tdp_env, xas_tdp_control, qs_env)
! Once assigned, diagonalize the MOs wrt the KS matrix in the subspace associated to each atom
IF (xas_tdp_control%do_loc) THEN
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T5,A)") &
"# Diagonalize localized MOs wrt the KS matrix in the subspace of each excited", &
"atom for better donor state identification."
END IF
CALL diagonalize_assigned_mo_subset(xas_tdp_env, xas_tdp_control, qs_env)
! update MO centers
CALL find_mo_centers(xas_tdp_env, xas_tdp_control, qs_env)
END IF
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,T3,A,I4,A,/)") &
"# Assign the relevant subset of the ", xas_tdp_control%n_search, &
" lowest energy MOs to excited atoms"
END IF
CALL write_mos_to_ex_atoms_association(xas_tdp_env, xas_tdp_control, qs_env)
! If CHECK_ONLY run, check the donor MOs
IF (xas_tdp_control%check_only) CALL print_checks(xas_tdp_env, xas_tdp_control, qs_env)
! If not simply exact exchange, setup a xas_atom_env and compute the xc integrals on the atomic grids
! Also needed if SOC is included or XPS GW2X(). Done before looping on atoms as it's all done at once
IF ((xas_tdp_control%do_xc .OR. xas_tdp_control%do_soc .OR. xas_tdp_control%do_gw2x) &
.AND. .NOT. xas_tdp_control%check_only) THEN
IF (output_unit > 0 .AND. xas_tdp_control%do_xc) THEN
WRITE (UNIT=output_unit, FMT="(/,T3,A,I4,A)") &
"# Integrating the xc kernel on the atomic grids ..."
CALL m_flush(output_unit)
END IF
CALL xas_atom_env_create(xas_atom_env)
CALL init_xas_atom_env(xas_atom_env, xas_tdp_env, xas_tdp_control, qs_env)
do_os = xas_tdp_control%do_uks .OR. xas_tdp_control%do_roks
IF (xas_tdp_control%do_xc .AND. (.NOT. xas_tdp_control%xps_only)) THEN
CALL integrate_fxc_atoms(xas_tdp_env%ri_fxc, xas_atom_env, xas_tdp_control, qs_env)
END IF
IF (xas_tdp_control%do_soc .OR. xas_tdp_control%do_gw2x) THEN
CALL integrate_soc_atoms(xas_tdp_env%orb_soc, xas_atom_env=xas_atom_env, qs_env=qs_env)
END IF
CALL xas_atom_env_release(xas_atom_env)
END IF
! Compute the 3-center Coulomb integrals for the whole system
IF ((.NOT. (xas_tdp_control%check_only .OR. xas_tdp_control%xps_only)) .AND. &
(xas_tdp_control%do_xc .OR. xas_tdp_control%do_coulomb)) THEN
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,T3,A,I4,A)") &
"# Computing the RI 3-center Coulomb integrals ..."
CALL m_flush(output_unit)
END IF
CALL compute_ri_3c_coulomb(xas_tdp_env, qs_env)
END IF
! Loop over donor states for calculation
CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
current_state_index = 1
! Loop over atomic kinds
DO ikind = 1, SIZE(atomic_kind_set)
IF (xas_tdp_control%check_only) EXIT
IF (.NOT. ANY(xas_tdp_env%ex_kind_indices == ikind)) CYCLE
CALL get_atomic_kind(atomic_kind=atomic_kind_set(ikind), name=kind_name, &
atom_list=atoms_of_kind)
! compute the RI coulomb2 inverse for this kind, and RI exchange2 if needed
CALL compute_ri_coulomb2_int(ikind, xas_tdp_env, xas_tdp_control, qs_env)
IF (xas_tdp_control%do_hfx) THEN
CALL compute_ri_exchange2_int(ikind, xas_tdp_env, xas_tdp_control, qs_env)
END IF
!Randomly distribute excited atoms of current kinds into batches for optimal load balance
!of RI 3c exchange integrals. Take batch sizes of 2 to avoid taxing memory too much, while
!greatly improving load balance
batch_size = 2
CALL get_ri_3c_batches(ex_atoms_of_kind, nbatch, batch_size, atoms_of_kind, xas_tdp_env)
nex_atom = SIZE(ex_atoms_of_kind)
!Loop over batches
DO ibatch = 1, nbatch
bo = get_limit(nex_atom, nbatch, ibatch - 1)
batch_size = bo(2) - bo(1) + 1
ALLOCATE (batch_atoms(batch_size))
iatom = 0
DO iat = bo(1), bo(2)
iatom = iatom + 1
batch_atoms(iatom) = ex_atoms_of_kind(iat)
END DO
CALL sort_unique(batch_atoms, unique)
!compute RI 3c exchange integrals on batch, if so required
IF (xas_tdp_control%do_hfx) THEN
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,T3,A,I4,A,I4,A,I1,A,A)") &
"# Computing the RI 3-center Exchange integrals for batch ", ibatch, "(/", nbatch, ") of ", &
batch_size, " atoms of kind: ", TRIM(kind_name)
CALL m_flush(output_unit)
END IF
CALL compute_ri_3c_exchange(batch_atoms, xas_tdp_env, xas_tdp_control, qs_env)
END IF
! Loop over atoms of batch
DO iat = 1, batch_size
iatom = batch_atoms(iat)
tmp_index = locate(xas_tdp_env%ex_atom_indices, iatom)
!if dipole in length rep, compute the dipole in the AO basis for this atom
!if quadrupole is required, compute it there too (in length rep)
IF (xas_tdp_control%dipole_form == xas_dip_len .OR. xas_tdp_control%do_quad) THEN
CALL compute_lenrep_multipole(iatom, xas_tdp_env, xas_tdp_control, qs_env)
END IF
! Loop over states of excited atom of kind
DO istate = 1, SIZE(xas_tdp_env%state_types, 1)
IF (xas_tdp_env%state_types(istate, tmp_index) == xas_not_excited) CYCLE
current_state => xas_tdp_env%donor_states(current_state_index)
CALL set_donor_state(current_state, at_index=iatom, &
at_symbol=kind_name, kind_index=ikind, &
state_type=xas_tdp_env%state_types(istate, tmp_index))
! Initial write for the donor state of interest
IF (output_unit > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,T3,A,A2,A,I4,A,A,/)") &
"# Start of calculations for donor state of type ", &
xas_tdp_env%state_type_char(current_state%state_type), " for atom", &
current_state%at_index, " of kind ", TRIM(current_state%at_symbol)
CALL m_flush(output_unit)
END IF
! Assign best fitting MO(s) to current core donnor state
CALL assign_mos_to_donor_state(current_state, xas_tdp_env, xas_tdp_control, qs_env)
! Perform MO restricted Mulliken pop analysis for verification
CALL perform_mulliken_on_donor_state(current_state, qs_env)
! GW2X correction
IF (xas_tdp_control%do_gw2x) THEN
CALL GW2X_shift(current_state, xas_tdp_env, xas_tdp_control, qs_env)
END IF
! Do main XAS calculations here
IF (.NOT. xas_tdp_control%xps_only) THEN
CALL setup_xas_tdp_prob(current_state, qs_env, xas_tdp_env, xas_tdp_control)
IF (xas_tdp_control%do_spin_cons) THEN
CALL solve_xas_tdp_prob(current_state, xas_tdp_control, xas_tdp_env, qs_env, &
ex_type=tddfpt_spin_cons)
CALL compute_dipole_fosc(current_state, xas_tdp_control, xas_tdp_env)
IF (xas_tdp_control%do_quad) CALL compute_quadrupole_fosc(current_state, &
xas_tdp_control, xas_tdp_env)
CALL xas_tdp_post(tddfpt_spin_cons, current_state, xas_tdp_env, &
xas_tdp_section, qs_env)
CALL write_donor_state_restart(tddfpt_spin_cons, current_state, xas_tdp_section, qs_env)
END IF
IF (xas_tdp_control%do_spin_flip) THEN
CALL solve_xas_tdp_prob(current_state, xas_tdp_control, xas_tdp_env, qs_env, &
ex_type=tddfpt_spin_flip)
!no dipole in spin-flip (spin-forbidden)
CALL xas_tdp_post(tddfpt_spin_flip, current_state, xas_tdp_env, &
xas_tdp_section, qs_env)
CALL write_donor_state_restart(tddfpt_spin_flip, current_state, xas_tdp_section, qs_env)
END IF
IF (xas_tdp_control%do_singlet) THEN
CALL solve_xas_tdp_prob(current_state, xas_tdp_control, xas_tdp_env, qs_env, &
ex_type=tddfpt_singlet)
CALL compute_dipole_fosc(current_state, xas_tdp_control, xas_tdp_env)
IF (xas_tdp_control%do_quad) CALL compute_quadrupole_fosc(current_state, &
xas_tdp_control, xas_tdp_env)
CALL xas_tdp_post(tddfpt_singlet, current_state, xas_tdp_env, &
xas_tdp_section, qs_env)
CALL write_donor_state_restart(tddfpt_singlet, current_state, xas_tdp_section, qs_env)
END IF
IF (xas_tdp_control%do_triplet) THEN
CALL solve_xas_tdp_prob(current_state, xas_tdp_control, xas_tdp_env, qs_env, &
ex_type=tddfpt_triplet)
!no dipole for triplets by construction
CALL xas_tdp_post(tddfpt_triplet, current_state, xas_tdp_env, &
xas_tdp_section, qs_env)
CALL write_donor_state_restart(tddfpt_triplet, current_state, xas_tdp_section, qs_env)
END IF
! Include the SOC if required, only for 2p donor stataes
IF (xas_tdp_control%do_soc .AND. current_state%state_type == xas_2p_type) THEN
IF (xas_tdp_control%do_singlet .AND. xas_tdp_control%do_triplet) THEN
CALL include_rcs_soc(current_state, xas_tdp_env, xas_tdp_control, qs_env)
END IF
IF (xas_tdp_control%do_spin_cons .AND. xas_tdp_control%do_spin_flip) THEN
CALL include_os_soc(current_state, xas_tdp_env, xas_tdp_control, qs_env)
END IF
END IF
! Print the requested properties
CALL print_xas_tdp_to_file(current_state, xas_tdp_env, xas_tdp_control, xas_tdp_section)
END IF !xps_only
IF (xas_tdp_control%do_gw2x) CALL print_xps(current_state, xas_tdp_env, xas_tdp_control, qs_env)
! Free some unneeded attributes of current_state
CALL free_ds_memory(current_state)
current_state_index = current_state_index + 1
NULLIFY (current_state)
END DO ! state type
end_of_batch = .FALSE.
IF (iat == batch_size) end_of_batch = .TRUE.
CALL free_exat_memory(xas_tdp_env, iatom, end_of_batch)
END DO ! atom of batch
DEALLOCATE (batch_atoms)
END DO !ibatch
DEALLOCATE (ex_atoms_of_kind)
END DO ! kind
! Return to ususal KS matrix
IF (dft_control%do_admm) THEN
CALL get_qs_env(qs_env, matrix_ks=matrix_ks, admm_env=admm_env)
DO ispin = 1, dft_control%nspins
CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, matrix_ks(ispin)%matrix)
END DO
END IF
! Return to initial basis set radii
IF (xas_tdp_control%eps_pgf > 0.0_dp) THEN
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
DO ikind = 1, SIZE(atomic_kind_set)
CALL get_qs_kind(qs_kind_set(ikind), basis_set=tmp_basis, basis_type="ORB")
CALL init_interaction_radii_orb_basis(tmp_basis, eps_pgf_orb=dft_control%qs_control%eps_pgf_orb)
END DO
END IF
! Clean-up
CALL xas_tdp_env_release(xas_tdp_env)
CALL xas_tdp_control_release(xas_tdp_control)
END SUBROUTINE xas_tdp_core
! **************************************************************************************************
!> \brief Overall control and environment types initialization
!> \param xas_tdp_env the environment type to initialize
!> \param xas_tdp_control the control type to initialize
!> \param qs_env the inherited qs environement type
! **************************************************************************************************
SUBROUTINE xas_tdp_init(xas_tdp_env, xas_tdp_control, qs_env)
TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(LEN=default_string_length) :: kind_name
INTEGER :: at_ind, i, ispin, j, k, kind_ind, &
n_donor_states, n_kinds, nao, &
nat_of_kind, natom, nex_atoms, &
nex_kinds, nmatch, nspins
INTEGER, DIMENSION(2) :: homo, n_mo, n_moloc
INTEGER, DIMENSION(:), POINTER :: ind_of_kind
LOGICAL :: do_os, do_uks, unique
REAL(dp) :: fact
REAL(dp), DIMENSION(:), POINTER :: mo_evals
TYPE(admm_type), POINTER :: admm_env
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: at_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cp_fm_type), POINTER :: mo_coeff
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
TYPE(dbcsr_type) :: matrix_tmp
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_type), POINTER :: tmp_basis
TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_loc_env_type), POINTER :: qs_loc_env
TYPE(qs_rho_type), POINTER :: rho
TYPE(section_type), POINTER :: dummy_section
TYPE(section_vals_type), POINTER :: loc_section, xas_tdp_section
NULLIFY (xas_tdp_section, at_kind_set, ind_of_kind, dft_control, qs_kind_set, tmp_basis)
NULLIFY (qs_loc_env, loc_section, mos, particle_set, rho, rho_ao, mo_evals, cell)
NULLIFY (mo_coeff, matrix_ks, admm_env, dummy_section)
! XAS TDP control type initialization
xas_tdp_section => section_vals_get_subs_vals(qs_env%input, "DFT%XAS_TDP")
CALL get_qs_env(qs_env, dft_control=dft_control)
CALL xas_tdp_control_create(xas_tdp_control)
CALL read_xas_tdp_control(xas_tdp_control, xas_tdp_section)
! Check the qs_env for a LSD/ROKS calculation
IF (dft_control%uks) xas_tdp_control%do_uks = .TRUE.
IF (dft_control%roks) xas_tdp_control%do_roks = .TRUE.
do_uks = xas_tdp_control%do_uks
do_os = do_uks .OR. xas_tdp_control%do_roks
! XAS TDP environment type initialization
CALL xas_tdp_env_create(xas_tdp_env)
! Retrieving the excited atoms indices and correspondig state types
IF (xas_tdp_control%define_excited == xas_tdp_by_index) THEN
! simply copy indices from xas_tdp_control
nex_atoms = SIZE(xas_tdp_control%list_ex_atoms)
CALL set_xas_tdp_env(xas_tdp_env, nex_atoms=nex_atoms)
ALLOCATE (xas_tdp_env%ex_atom_indices(nex_atoms))
ALLOCATE (xas_tdp_env%state_types(SIZE(xas_tdp_control%state_types, 1), nex_atoms))
xas_tdp_env%ex_atom_indices = xas_tdp_control%list_ex_atoms
xas_tdp_env%state_types = xas_tdp_control%state_types
! Test that these indices are within the range of available atoms
CALL get_qs_env(qs_env=qs_env, natom=natom)
IF (ANY(xas_tdp_env%ex_atom_indices > natom)) THEN
CPABORT("Invalid index for the ATOM_LIST keyword.")
END IF
! Check atom kinds and fill corresponding array
ALLOCATE (xas_tdp_env%ex_kind_indices(nex_atoms))
xas_tdp_env%ex_kind_indices = 0
k = 0
CALL get_qs_env(qs_env, particle_set=particle_set)
DO i = 1, nex_atoms
at_ind = xas_tdp_env%ex_atom_indices(i)
CALL get_atomic_kind(particle_set(at_ind)%atomic_kind, kind_number=j)
IF (ALL(ABS(xas_tdp_env%ex_kind_indices - j) .NE. 0)) THEN
k = k + 1
xas_tdp_env%ex_kind_indices(k) = j
END IF
END DO
nex_kinds = k
CALL set_xas_tdp_env(xas_tdp_env, nex_kinds=nex_kinds)
CALL reallocate(xas_tdp_env%ex_kind_indices, 1, nex_kinds)
ELSE IF (xas_tdp_control%define_excited == xas_tdp_by_kind) THEN
! need to find out which atom of which kind is excited
CALL get_qs_env(qs_env=qs_env, atomic_kind_set=at_kind_set)
n_kinds = SIZE(at_kind_set)
nex_atoms = 0
nex_kinds = SIZE(xas_tdp_control%list_ex_kinds)
ALLOCATE (xas_tdp_env%ex_kind_indices(nex_kinds))
k = 0
DO i = 1, n_kinds
CALL get_atomic_kind(atomic_kind=at_kind_set(i), name=kind_name, &
natom=nat_of_kind, kind_number=kind_ind)
IF (ANY(xas_tdp_control%list_ex_kinds == kind_name)) THEN
nex_atoms = nex_atoms + nat_of_kind
k = k + 1
xas_tdp_env%ex_kind_indices(k) = kind_ind
END IF
END DO
ALLOCATE (xas_tdp_env%ex_atom_indices(nex_atoms))
ALLOCATE (xas_tdp_env%state_types(SIZE(xas_tdp_control%state_types, 1), nex_atoms))
nex_atoms = 0
nmatch = 0
DO i = 1, n_kinds
CALL get_atomic_kind(atomic_kind=at_kind_set(i), name=kind_name, &
natom=nat_of_kind, atom_list=ind_of_kind)
DO j = 1, nex_kinds
IF (xas_tdp_control%list_ex_kinds(j) == kind_name) THEN
xas_tdp_env%ex_atom_indices(nex_atoms + 1:nex_atoms + nat_of_kind) = ind_of_kind
DO k = 1, SIZE(xas_tdp_control%state_types, 1)
xas_tdp_env%state_types(k, nex_atoms + 1:nex_atoms + nat_of_kind) = &
xas_tdp_control%state_types(k, j)
END DO
nex_atoms = nex_atoms + nat_of_kind
nmatch = nmatch + 1
END IF
END DO
END DO
CALL set_xas_tdp_env(xas_tdp_env, nex_atoms=nex_atoms, nex_kinds=nex_kinds)
! Verifying that the input was valid
IF (nmatch .NE. SIZE(xas_tdp_control%list_ex_kinds)) THEN
CPABORT("Invalid kind(s) for the KIND_LIST keyword.")
END IF
END IF
! Sort the excited atoms indices (for convinience and use of locate function)
CALL sort_unique(xas_tdp_env%ex_atom_indices, unique)
IF (.NOT. unique) THEN
CPABORT("Excited atoms not uniquely defined.")
END IF
! Check for periodicity
CALL get_qs_env(qs_env, cell=cell)
IF (ALL(cell%perd == 0)) THEN
xas_tdp_control%is_periodic = .FALSE.
ELSE IF (ALL(cell%perd == 1)) THEN
xas_tdp_control%is_periodic = .TRUE.
ELSE
CPABORT("XAS TDP only implemented for full PBCs or non-PBCs")
END IF
! Allocating memory for the array of donor states
n_donor_states = COUNT(xas_tdp_env%state_types /= xas_not_excited)
ALLOCATE (xas_tdp_env%donor_states(n_donor_states))
DO i = 1, n_donor_states
CALL donor_state_create(xas_tdp_env%donor_states(i))
END DO
! In case of ADMM, for the whole duration of the XAS_TDP, we need the total KS matrix
IF (dft_control%do_admm) THEN
CALL get_qs_env(qs_env, admm_env=admm_env, matrix_ks=matrix_ks)
DO ispin = 1, dft_control%nspins
CALL admm_correct_for_eigenvalues(ispin, admm_env, matrix_ks(ispin)%matrix)
END DO
END IF
! In case of externally imposed EPS_PGF_XAS, need to update ORB and RI_XAS interaction radii
IF (xas_tdp_control%eps_pgf > 0.0_dp) THEN
CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
DO i = 1, SIZE(qs_kind_set)
CALL get_qs_kind(qs_kind_set(i), basis_set=tmp_basis, basis_type="ORB")
CALL init_interaction_radii_orb_basis(tmp_basis, eps_pgf_orb=xas_tdp_control%eps_pgf)
CALL get_qs_kind(qs_kind_set(i), basis_set=tmp_basis, basis_type="RI_XAS")
CALL init_interaction_radii_orb_basis(tmp_basis, eps_pgf_orb=xas_tdp_control%eps_pgf)
END DO
END IF
! In case of ground state OT optimization, compute the MO eigenvalues and get canonical MOs
IF (qs_env%scf_env%method == ot_method_nr) THEN
CALL get_qs_env(qs_env, mos=mos, matrix_ks=matrix_ks)
nspins = 1; IF (do_uks) nspins = 2
DO ispin = 1, nspins
CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, eigenvalues=mo_evals)
CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, evals_arg=mo_evals)
END DO
END IF
! Initializing the qs_loc_env from the LOCALIZE subsection of XAS_TDP (largely inpired by MI's XAS)
! We create the LOCALIZE subsection here, since it is completely overwritten anyways
CALL create_localize_section(dummy_section)
CALL section_vals_create(xas_tdp_control%loc_subsection, dummy_section)
CALL section_vals_val_set(xas_tdp_control%loc_subsection, "_SECTION_PARAMETERS_", &
l_val=xas_tdp_control%do_loc)
CALL section_release(dummy_section)
xas_tdp_control%print_loc_subsection => section_vals_get_subs_vals( &
xas_tdp_control%loc_subsection, "PRINT")
ALLOCATE (xas_tdp_env%qs_loc_env)
CALL qs_loc_env_create(xas_tdp_env%qs_loc_env)
qs_loc_env => xas_tdp_env%qs_loc_env
loc_section => xas_tdp_control%loc_subsection
! getting the number of MOs
CALL get_qs_env(qs_env, mos=mos)
CALL get_mo_set(mos(1), nmo=n_mo(1), homo=homo(1), nao=nao)
n_mo(2) = n_mo(1)
homo(2) = homo(1)
nspins = 1
IF (do_os) CALL get_mo_set(mos(2), nmo=n_mo(2), homo=homo(2))
IF (do_uks) nspins = 2 !in roks, same MOs for both spins
! by default, all (doubly occupied) homo are localized
IF (xas_tdp_control%n_search < 0 .OR. xas_tdp_control%n_search > MINVAL(homo)) &
xas_tdp_control%n_search = MINVAL(homo)
CALL qs_loc_control_init(qs_loc_env, loc_section, do_homo=.TRUE., do_xas=.TRUE., &
nloc_xas=xas_tdp_control%n_search, spin_xas=1)
! do_xas argument above only prepares spin-alpha localization
IF (do_uks) THEN
qs_loc_env%localized_wfn_control%nloc_states(2) = xas_tdp_control%n_search
qs_loc_env%localized_wfn_control%lu_bound_states(1, 2) = 1
qs_loc_env%localized_wfn_control%lu_bound_states(2, 2) = xas_tdp_control%n_search
END IF
! final qs_loc_env initialization. Impose Berry operator
qs_loc_env%localized_wfn_control%operator_type = op_loc_berry
qs_loc_env%localized_wfn_control%max_iter = 25000
IF (.NOT. xas_tdp_control%do_loc) THEN
qs_loc_env%localized_wfn_control%localization_method = do_loc_none
ELSE
n_moloc = qs_loc_env%localized_wfn_control%nloc_states
CALL set_loc_centers(qs_loc_env%localized_wfn_control, n_moloc, nspins)
IF (do_uks) THEN
CALL qs_loc_env_init(qs_loc_env, qs_loc_env%localized_wfn_control, &
qs_env, do_localize=.TRUE.)
ELSE
CALL qs_loc_env_init(qs_loc_env, qs_loc_env%localized_wfn_control, &
qs_env, do_localize=.TRUE., myspin=1)
END IF
END IF
! Allocating memory for the array of excited atoms MOs. Worst case senario, all searched MOs are
! associated to the same atom
ALLOCATE (xas_tdp_env%mos_of_ex_atoms(xas_tdp_control%n_search, nex_atoms, nspins))
! Compute the projector on the unoccupied, unperturbed ground state: Q = 1 - SP, sor each spin
IF (do_os) nspins = 2
CALL get_qs_env(qs_env, rho=rho, matrix_s=matrix_s)
CALL qs_rho_get(rho, rho_ao=rho_ao)
ALLOCATE (xas_tdp_env%q_projector(nspins))
ALLOCATE (xas_tdp_env%q_projector(1)%matrix)
CALL dbcsr_create(xas_tdp_env%q_projector(1)%matrix, name="Q PROJECTOR ALPHA", &
template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
IF (do_os) THEN
ALLOCATE (xas_tdp_env%q_projector(2)%matrix)
CALL dbcsr_create(xas_tdp_env%q_projector(2)%matrix, name="Q PROJECTOR BETA", &
template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
END IF
! In the case of spin-restricted calculations, rho_ao includes double occupency => 0.5 prefactor
fact = -0.5_dp; IF (do_os) fact = -1.0_dp
CALL dbcsr_multiply('N', 'N', fact, matrix_s(1)%matrix, rho_ao(1)%matrix, 0.0_dp, &
xas_tdp_env%q_projector(1)%matrix, filter_eps=xas_tdp_control%eps_filter)
CALL dbcsr_add_on_diag(xas_tdp_env%q_projector(1)%matrix, 1.0_dp)
CALL dbcsr_finalize(xas_tdp_env%q_projector(1)%matrix)
IF (do_os) THEN
CALL dbcsr_multiply('N', 'N', fact, matrix_s(1)%matrix, rho_ao(2)%matrix, 0.0_dp, &
xas_tdp_env%q_projector(2)%matrix, filter_eps=xas_tdp_control%eps_filter)
CALL dbcsr_add_on_diag(xas_tdp_env%q_projector(2)%matrix, 1.0_dp)
CALL dbcsr_finalize(xas_tdp_env%q_projector(2)%matrix)
END IF
! Create the structure for the dipole in the AO basis
ALLOCATE (xas_tdp_env%dipmat(3))
DO i = 1, 3
ALLOCATE (xas_tdp_env%dipmat(i)%matrix)
CALL dbcsr_copy(matrix_tmp, matrix_s(1)%matrix, name="XAS TDP dipole matrix")
IF (xas_tdp_control%dipole_form == xas_dip_vel) THEN
CALL dbcsr_create(xas_tdp_env%dipmat(i)%matrix, template=matrix_s(1)%matrix, &
matrix_type=dbcsr_type_antisymmetric)
CALL dbcsr_complete_redistribute(matrix_tmp, xas_tdp_env%dipmat(i)%matrix)
ELSE
CALL dbcsr_create(xas_tdp_env%dipmat(i)%matrix, template=matrix_s(1)%matrix, &
matrix_type=dbcsr_type_symmetric)
CALL dbcsr_copy(xas_tdp_env%dipmat(i)%matrix, matrix_tmp)
END IF
CALL dbcsr_set(xas_tdp_env%dipmat(i)%matrix, 0.0_dp)
CALL dbcsr_release(matrix_tmp)
END DO
! Create the structure for the electric quadrupole in the AO basis, if required
IF (xas_tdp_control%do_quad) THEN
ALLOCATE (xas_tdp_env%quadmat(6))
DO i = 1, 6
ALLOCATE (xas_tdp_env%quadmat(i)%matrix)
CALL dbcsr_copy(xas_tdp_env%quadmat(i)%matrix, matrix_s(1)%matrix, name="XAS TDP quadrupole matrix")
CALL dbcsr_set(xas_tdp_env%quadmat(i)%matrix, 0.0_dp)
END DO
END IF
! Precompute it in the velocity representation, if so chosen
IF (xas_tdp_control%dipole_form == xas_dip_vel) THEN
!enforce minimum image to avoid any PBCs related issues. Ok because very localized densities
CALL p_xyz_ao(xas_tdp_env%dipmat, qs_env, minimum_image=.TRUE.)
END IF
! Allocate SOC in AO basis matrices
IF (xas_tdp_control%do_soc .OR. xas_tdp_control%do_gw2x) THEN
ALLOCATE (xas_tdp_env%orb_soc(3))
DO i = 1, 3
ALLOCATE (xas_tdp_env%orb_soc(i)%matrix)
END DO
END IF
! Check that everything is allowed
CALL safety_check(xas_tdp_control, qs_env)
! Initialize libint for the 3-center integrals
CALL cp_libint_static_init()
! Compute LUMOs as guess for OT solver and/or for GW2X correction
IF (xas_tdp_control%do_ot .OR. xas_tdp_control%do_gw2x) THEN
CALL make_lumo_guess(xas_tdp_env, xas_tdp_control, qs_env)
END IF
END SUBROUTINE xas_tdp_init
! **************************************************************************************************
!> \brief splits the excited atoms of a kind into batches for RI 3c integrals load balance
!> \param ex_atoms_of_kind the excited atoms for the current kind, randomly shuffled
!> \param nbatch number of batches to loop over
!> \param batch_size standard size of a batch
!> \param atoms_of_kind number of atoms for the current kind (excited or not)
!> \param xas_tdp_env ...
! **************************************************************************************************
SUBROUTINE get_ri_3c_batches(ex_atoms_of_kind, nbatch, batch_size, atoms_of_kind, xas_tdp_env)
INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: ex_atoms_of_kind
INTEGER, INTENT(OUT) :: nbatch
INTEGER, INTENT(IN) :: batch_size
INTEGER, DIMENSION(:), INTENT(IN) :: atoms_of_kind
TYPE(xas_tdp_env_type), POINTER :: xas_tdp_env
INTEGER :: iat, iatom, nex_atom
TYPE(rng_stream_type), ALLOCATABLE :: rng_stream
!Get the atoms from atoms_of_kind that are excited
nex_atom = 0
DO iat = 1, SIZE(atoms_of_kind)
iatom = atoms_of_kind(iat)
IF (.NOT. ANY(xas_tdp_env%ex_atom_indices == iatom)) CYCLE
nex_atom = nex_atom + 1
END DO
ALLOCATE (ex_atoms_of_kind(nex_atom))
nex_atom = 0
DO iat = 1, SIZE(atoms_of_kind)
iatom = atoms_of_kind(iat)
IF (.NOT. ANY(xas_tdp_env%ex_atom_indices == iatom)) CYCLE
nex_atom = nex_atom + 1
ex_atoms_of_kind(nex_atom) = iatom
END DO
!We shuffle those atoms to spread them
rng_stream = rng_stream_type(name="uniform_rng", distribution_type=UNIFORM)
CALL rng_stream%shuffle(ex_atoms_of_kind(1:nex_atom))
nbatch = nex_atom/batch_size
IF (nbatch*batch_size .NE. nex_atom) nbatch = nbatch + 1
END SUBROUTINE get_ri_3c_batches
! **************************************************************************************************
!> \brief Checks for forbidden keywords combinations
!> \param xas_tdp_control ...
!> \param qs_env ...
! **************************************************************************************************
SUBROUTINE safety_check(xas_tdp_control, qs_env)
TYPE(xas_tdp_control_type), POINTER :: xas_tdp_control
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dft_control_type), POINTER :: dft_control
!PB only available without exact exchange
IF (xas_tdp_control%is_periodic .AND. xas_tdp_control%do_hfx &
.AND. xas_tdp_control%x_potential%potential_type == do_potential_coulomb) THEN
CPABORT("XAS TDP with Coulomb operator for exact exchange only supports non-periodic BCs")
END IF
!open-shell/closed-shell tests
IF (xas_tdp_control%do_roks .OR. xas_tdp_control%do_uks) THEN
IF (.NOT. (xas_tdp_control%do_spin_cons .OR. xas_tdp_control%do_spin_flip)) THEN
CPABORT("Need spin-conserving and/or spin-flip excitations for open-shell systems")
END IF
IF (xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet) THEN
CPABORT("Singlet/triplet excitations only for restricted closed-shell systems")
END IF
IF (xas_tdp_control%do_soc .AND. .NOT. &
(xas_tdp_control%do_spin_flip .AND. xas_tdp_control%do_spin_cons)) THEN
CPABORT("Both spin-conserving and spin-flip excitations are required for SOC")
END IF
ELSE
IF (.NOT. (xas_tdp_control%do_singlet .OR. xas_tdp_control%do_triplet)) THEN
CPABORT("Need singlet and/or triplet excitations for closed-shell systems")
END IF
IF (xas_tdp_control%do_spin_cons .OR. xas_tdp_control%do_spin_flip) THEN
CPABORT("Spin-conserving/spin-flip excitations only for open-shell systems")
END IF
IF (xas_tdp_control%do_soc .AND. .NOT. &
(xas_tdp_control%do_singlet .AND. xas_tdp_control%do_triplet)) THEN
CPABORT("Both singlet and triplet excitations are needed for SOC")
END IF
END IF
!Warn against using E_RANGE with SOC
IF (xas_tdp_control%do_soc .AND. xas_tdp_control%e_range > 0.0_dp) THEN
CPWARN("Using E_RANGE and SOC together may lead to crashes, use N_EXCITED for safety.")
END IF
!TDA, full-TDDFT and diagonalization
IF (.NOT. xas_tdp_control%tamm_dancoff) THEN
IF (xas_tdp_control%do_spin_flip) THEN
CPABORT("Spin-flip kernel only implemented for Tamm-Dancoff approximation")
END IF
IF (xas_tdp_control%do_ot) THEN
CPABORT("OT diagonalization only available within the Tamm-Dancoff approximation")
END IF
END IF