-
Notifications
You must be signed in to change notification settings - Fork 1
/
hfx_ri.F
4150 lines (3453 loc) · 198 KB
/
hfx_ri.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 RI-methods for HFX
! **************************************************************************************************
MODULE hfx_ri
USE OMP_LIB, ONLY: omp_get_num_threads,&
omp_get_thread_num
USE arnoldi_api, ONLY: arnoldi_extremal
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind_set
USE basis_set_types, ONLY: gto_basis_set_p_type,&
gto_basis_set_type
USE cell_types, ONLY: cell_type
USE cp_blacs_env, ONLY: cp_blacs_env_type
USE cp_control_types, ONLY: dft_control_type
USE cp_dbcsr_api, ONLY: &
dbcsr_add, dbcsr_add_on_diag, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, &
dbcsr_desymmetrize, dbcsr_distribution_get, dbcsr_distribution_release, &
dbcsr_distribution_type, dbcsr_dot, dbcsr_filter, dbcsr_frobenius_norm, dbcsr_get_info, &
dbcsr_get_num_blocks, dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, &
dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
cp_dbcsr_cholesky_invert
USE cp_dbcsr_diag, ONLY: cp_dbcsr_power
USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
copy_fm_to_dbcsr,&
cp_dbcsr_dist2d_to_dist,&
dbcsr_deallocate_matrix_set
USE cp_fm_struct, ONLY: cp_fm_struct_create,&
cp_fm_struct_release,&
cp_fm_struct_type
USE cp_fm_types, ONLY: cp_fm_create,&
cp_fm_release,&
cp_fm_type,&
cp_fm_write_unformatted
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 dbt_api, ONLY: &
dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_clear, dbt_contract, &
dbt_copy, dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, &
dbt_default_distvec, dbt_destroy, dbt_distribution_destroy, dbt_distribution_new, &
dbt_distribution_type, dbt_filter, dbt_get_block, dbt_get_info, dbt_get_num_blocks_total, &
dbt_iterator_blocks_left, dbt_iterator_next_block, dbt_iterator_start, dbt_iterator_stop, &
dbt_iterator_type, dbt_mp_environ_pgrid, dbt_nd_mp_comm, dbt_pgrid_create, &
dbt_pgrid_destroy, dbt_pgrid_type, dbt_put_block, dbt_reserve_blocks, dbt_scale, dbt_type
USE distribution_2d_types, ONLY: distribution_2d_type
USE hfx_types, ONLY: alloc_containers,&
block_ind_type,&
dealloc_containers,&
hfx_compression_type,&
hfx_ri_init,&
hfx_ri_release,&
hfx_ri_type
USE input_constants, ONLY: hfx_ri_do_2c_cholesky,&
hfx_ri_do_2c_diag,&
hfx_ri_do_2c_iter
USE input_cp2k_hfx, ONLY: ri_mo,&
ri_pmat
USE input_section_types, ONLY: section_vals_get_subs_vals,&
section_vals_type,&
section_vals_val_get
USE iterate_matrix, ONLY: invert_hotelling,&
matrix_sqrt_newton_schulz
USE kinds, ONLY: default_string_length,&
dp,&
int_8
USE machine, ONLY: m_walltime
USE message_passing, ONLY: mp_cart_type,&
mp_comm_type,&
mp_para_env_type
USE particle_methods, ONLY: get_particle_set
USE particle_types, ONLY: particle_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_force_types, ONLY: qs_force_type
USE qs_integral_utils, ONLY: basis_set_list_setup
USE qs_interactions, ONLY: init_interaction_radii_orb_basis
USE qs_kind_types, ONLY: qs_kind_type
USE qs_ks_types, ONLY: qs_ks_env_type
USE qs_mo_types, ONLY: get_mo_set,&
mo_set_type
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
release_neighbor_list_sets
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE qs_tensors, ONLY: &
build_2c_derivatives, build_2c_integrals, build_2c_neighbor_lists, build_3c_derivatives, &
build_3c_integrals, build_3c_neighbor_lists, calc_2c_virial, calc_3c_virial, &
compress_tensor, decompress_tensor, get_tensor_occupancy, neighbor_list_3c_destroy
USE qs_tensors_types, ONLY: create_2c_tensor,&
create_3c_tensor,&
create_tensor_batches,&
distribution_3d_create,&
distribution_3d_type,&
neighbor_list_3c_type,&
split_block_sizes
USE util, ONLY: sort
USE virial_types, ONLY: virial_type
#include "./base/base_uses.f90"
!$ USE OMP_LIB, ONLY: omp_get_num_threads
IMPLICIT NONE
PRIVATE
PUBLIC :: hfx_ri_update_ks, hfx_ri_update_forces, get_force_from_3c_trace, get_2c_der_force, &
get_idx_to_atom, print_ri_hfx, hfx_ri_pre_scf_calc_tensors, check_inverse
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_ri'
CONTAINS
! **************************************************************************************************
!> \brief Switches the RI_FLAVOR from MO to RHO or vice-versa
!> \param ri_data ...
!> \param qs_env ...
!> \note As a side product, the ri_data is mostly reinitialized and the integrals recomputed
! **************************************************************************************************
SUBROUTINE switch_ri_flavor(ri_data, qs_env)
TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'switch_ri_flavor'
INTEGER :: handle, n_mem, new_flavor
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
NULLIFY (qs_kind_set, particle_set, atomic_kind_set, para_env, dft_control)
CALL timeset(routineN, handle)
CALL get_qs_env(qs_env, para_env=para_env, dft_control=dft_control, atomic_kind_set=atomic_kind_set, &
particle_set=particle_set, qs_kind_set=qs_kind_set)
CALL hfx_ri_release(ri_data, write_stats=.FALSE.)
IF (ri_data%flavor == ri_pmat) THEN
new_flavor = ri_mo
ELSE
new_flavor = ri_pmat
END IF
ri_data%flavor = new_flavor
n_mem = ri_data%n_mem_input
ri_data%n_mem_input = ri_data%n_mem_flavor_switch
ri_data%n_mem_flavor_switch = n_mem
CALL hfx_ri_init(ri_data, qs_kind_set, particle_set, atomic_kind_set, para_env)
!Need to recalculate the integrals in the new flavor
!TODO: should we backup the integrals and symmetrize/desymmetrize them instead of recomputing ?!?
! that only makes sense if actual integral calculation is not negligible
IF (ri_data%flavor == ri_pmat) THEN
CALL hfx_ri_pre_scf_Pmat(qs_env, ri_data)
ELSE
CALL hfx_ri_pre_scf_mo(qs_env, ri_data, dft_control%nspins)
END IF
IF (ri_data%unit_nr > 0) THEN
IF (ri_data%flavor == ri_pmat) THEN
WRITE (ri_data%unit_nr, '(T2,A)') "HFX_RI_INFO| temporarily switched to RI_FLAVOR RHO"
ELSE
WRITE (ri_data%unit_nr, '(T2,A)') "HFX_RI_INFO| temporarily switched to RI_FLAVOR MO"
END IF
END IF
CALL timestop(handle)
END SUBROUTINE switch_ri_flavor
! **************************************************************************************************
!> \brief Pre-SCF steps in MO flavor of RI HFX
!>
!> Calculate 2-center & 3-center integrals (see hfx_ri_pre_scf_calc_tensors) and contract
!> K(P, S) = sum_R K_2(P, R)^{-1} K_1(R, S)^{1/2}
!> B(mu, lambda, R) = sum_P int_3c(mu, lambda, P) K(P, R)
!> \param qs_env ...
!> \param ri_data ...
!> \param nspins ...
! **************************************************************************************************
SUBROUTINE hfx_ri_pre_scf_mo(qs_env, ri_data, nspins)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
INTEGER, INTENT(IN) :: nspins
CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ri_pre_scf_mo'
INTEGER :: handle, handle2, ispin, n_dependent, &
unit_nr, unit_nr_dbcsr
REAL(KIND=dp) :: threshold
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(dbcsr_type), DIMENSION(1) :: dbcsr_work_1, dbcsr_work_2, t_2c_int_mat, t_2c_op_pot, &
t_2c_op_pot_sqrt, t_2c_op_pot_sqrt_inv, t_2c_op_RI, t_2c_op_RI_inv
TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_2c_int, t_2c_work
TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_int
TYPE(mp_para_env_type), POINTER :: para_env
CALL timeset(routineN, handle)
unit_nr_dbcsr = ri_data%unit_nr_dbcsr
unit_nr = ri_data%unit_nr
CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
CALL timeset(routineN//"_int", handle2)
ALLOCATE (t_2c_int(1), t_2c_work(1), t_3c_int(1, 1))
CALL hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_op_RI, t_2c_op_pot, t_3c_int)
CALL timestop(handle2)
CALL timeset(routineN//"_2c", handle2)
IF (.NOT. ri_data%same_op) THEN
SELECT CASE (ri_data%t2c_method)
CASE (hfx_ri_do_2c_iter)
CALL dbcsr_create(t_2c_op_RI_inv(1), template=t_2c_op_RI(1), matrix_type=dbcsr_type_no_symmetry)
threshold = MAX(ri_data%filter_eps, 1.0e-12_dp)
CALL invert_hotelling(t_2c_op_RI_inv(1), t_2c_op_RI(1), threshold=threshold, silent=.FALSE.)
CASE (hfx_ri_do_2c_cholesky)
CALL dbcsr_copy(t_2c_op_RI_inv(1), t_2c_op_RI(1))
CALL cp_dbcsr_cholesky_decompose(t_2c_op_RI_inv(1), para_env=para_env, blacs_env=blacs_env)
CALL cp_dbcsr_cholesky_invert(t_2c_op_RI_inv(1), para_env=para_env, blacs_env=blacs_env, upper_to_full=.TRUE.)
CASE (hfx_ri_do_2c_diag)
CALL dbcsr_copy(t_2c_op_RI_inv(1), t_2c_op_RI(1))
CALL cp_dbcsr_power(t_2c_op_RI_inv(1), -1.0_dp, ri_data%eps_eigval, n_dependent, &
para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
END SELECT
IF (ri_data%check_2c_inv) THEN
CALL check_inverse(t_2c_op_RI_inv(1), t_2c_op_RI(1), unit_nr=unit_nr)
END IF
CALL dbcsr_release(t_2c_op_RI(1))
SELECT CASE (ri_data%t2c_method)
CASE (hfx_ri_do_2c_iter)
CALL dbcsr_create(t_2c_op_pot_sqrt(1), template=t_2c_op_pot(1), matrix_type=dbcsr_type_symmetric)
CALL dbcsr_create(t_2c_op_pot_sqrt_inv(1), template=t_2c_op_pot(1), matrix_type=dbcsr_type_symmetric)
CALL matrix_sqrt_newton_schulz(t_2c_op_pot_sqrt(1), t_2c_op_pot_sqrt_inv(1), t_2c_op_pot(1), &
ri_data%filter_eps, ri_data%t2c_sqrt_order, ri_data%eps_lanczos, &
ri_data%max_iter_lanczos)
CALL dbcsr_release(t_2c_op_pot_sqrt_inv(1))
CASE (hfx_ri_do_2c_diag, hfx_ri_do_2c_cholesky)
CALL dbcsr_copy(t_2c_op_pot_sqrt(1), t_2c_op_pot(1))
CALL cp_dbcsr_power(t_2c_op_pot_sqrt(1), 0.5_dp, ri_data%eps_eigval, n_dependent, &
para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
END SELECT
!We need S^-1 and (P|Q) for the forces.
CALL dbt_create(t_2c_op_RI_inv(1), t_2c_work(1))
CALL dbt_copy_matrix_to_tensor(t_2c_op_RI_inv(1), t_2c_work(1))
CALL dbt_copy(t_2c_work(1), ri_data%t_2c_inv(1, 1), move_data=.TRUE.)
CALL dbt_destroy(t_2c_work(1))
CALL dbt_filter(ri_data%t_2c_inv(1, 1), ri_data%filter_eps)
CALL dbt_create(t_2c_op_pot(1), t_2c_work(1))
CALL dbt_copy_matrix_to_tensor(t_2c_op_pot(1), t_2c_work(1))
CALL dbt_copy(t_2c_work(1), ri_data%t_2c_pot(1, 1), move_data=.TRUE.)
CALL dbt_destroy(t_2c_work(1))
CALL dbt_filter(ri_data%t_2c_pot(1, 1), ri_data%filter_eps)
IF (ri_data%check_2c_inv) THEN
CALL check_sqrt(t_2c_op_pot(1), matrix_sqrt=t_2c_op_pot_sqrt(1), unit_nr=unit_nr)
END IF
CALL dbcsr_create(t_2c_int_mat(1), template=t_2c_op_pot(1), matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_multiply("N", "N", 1.0_dp, t_2c_op_RI_inv(1), t_2c_op_pot_sqrt(1), &
0.0_dp, t_2c_int_mat(1), filter_eps=ri_data%filter_eps)
CALL dbcsr_release(t_2c_op_RI_inv(1))
CALL dbcsr_release(t_2c_op_pot_sqrt(1))
ELSE
SELECT CASE (ri_data%t2c_method)
CASE (hfx_ri_do_2c_iter)
CALL dbcsr_create(t_2c_int_mat(1), template=t_2c_op_pot(1), matrix_type=dbcsr_type_symmetric)
CALL dbcsr_create(t_2c_op_pot_sqrt(1), template=t_2c_op_pot(1), matrix_type=dbcsr_type_symmetric)
CALL matrix_sqrt_newton_schulz(t_2c_op_pot_sqrt(1), t_2c_int_mat(1), t_2c_op_pot(1), &
ri_data%filter_eps, ri_data%t2c_sqrt_order, ri_data%eps_lanczos, &
ri_data%max_iter_lanczos)
CALL dbcsr_release(t_2c_op_pot_sqrt(1))
CASE (hfx_ri_do_2c_diag, hfx_ri_do_2c_cholesky)
CALL dbcsr_copy(t_2c_int_mat(1), t_2c_op_pot(1))
CALL cp_dbcsr_power(t_2c_int_mat(1), -0.5_dp, ri_data%eps_eigval, n_dependent, &
para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
END SELECT
IF (ri_data%check_2c_inv) THEN
CALL check_sqrt(t_2c_op_pot(1), matrix_sqrt_inv=t_2c_int_mat(1), unit_nr=unit_nr)
END IF
!We need (P|Q)^-1 for the forces
CALL dbcsr_copy(dbcsr_work_1(1), t_2c_int_mat(1))
CALL dbcsr_create(dbcsr_work_2(1), template=t_2c_int_mat(1))
CALL dbcsr_multiply("N", "N", 1.0_dp, dbcsr_work_1(1), t_2c_int_mat(1), 0.0_dp, dbcsr_work_2(1))
CALL dbcsr_release(dbcsr_work_1(1))
CALL dbt_create(dbcsr_work_2(1), t_2c_work(1))
CALL dbt_copy_matrix_to_tensor(dbcsr_work_2(1), t_2c_work(1))
CALL dbcsr_release(dbcsr_work_2(1))
CALL dbt_copy(t_2c_work(1), ri_data%t_2c_inv(1, 1), move_data=.TRUE.)
CALL dbt_destroy(t_2c_work(1))
CALL dbt_filter(ri_data%t_2c_inv(1, 1), ri_data%filter_eps)
END IF
CALL dbcsr_release(t_2c_op_pot(1))
CALL dbt_create(t_2c_int_mat(1), t_2c_int(1), name="(RI|RI)")
CALL dbt_copy_matrix_to_tensor(t_2c_int_mat(1), t_2c_int(1))
CALL dbcsr_release(t_2c_int_mat(1))
DO ispin = 1, nspins
CALL dbt_copy(t_2c_int(1), ri_data%t_2c_int(ispin, 1))
END DO
CALL dbt_destroy(t_2c_int(1))
CALL timestop(handle2)
CALL timeset(routineN//"_3c", handle2)
CALL dbt_copy(t_3c_int(1, 1), ri_data%t_3c_int_ctr_1(1, 1), order=[2, 1, 3], move_data=.TRUE.)
CALL dbt_filter(ri_data%t_3c_int_ctr_1(1, 1), ri_data%filter_eps)
CALL dbt_copy(ri_data%t_3c_int_ctr_1(1, 1), ri_data%t_3c_int_ctr_2(1, 1))
CALL dbt_destroy(t_3c_int(1, 1))
CALL timestop(handle2)
CALL timestop(handle)
END SUBROUTINE
! **************************************************************************************************
!> \brief ...
!> \param matrix_1 ...
!> \param matrix_2 ...
!> \param name ...
!> \param unit_nr ...
! **************************************************************************************************
SUBROUTINE check_inverse(matrix_1, matrix_2, name, unit_nr)
TYPE(dbcsr_type), INTENT(INOUT) :: matrix_1, matrix_2
CHARACTER(len=*), INTENT(IN), OPTIONAL :: name
INTEGER, INTENT(IN) :: unit_nr
CHARACTER(len=default_string_length) :: name_prv
REAL(KIND=dp) :: error, frob_matrix, frob_matrix_base
TYPE(dbcsr_type) :: matrix_tmp
IF (PRESENT(name)) THEN
name_prv = name
ELSE
CALL dbcsr_get_info(matrix_1, name=name_prv)
END IF
CALL dbcsr_create(matrix_tmp, template=matrix_1)
CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_1, matrix_2, &
0.0_dp, matrix_tmp)
frob_matrix_base = dbcsr_frobenius_norm(matrix_tmp)
CALL dbcsr_add_on_diag(matrix_tmp, -1.0_dp)
frob_matrix = dbcsr_frobenius_norm(matrix_tmp)
error = frob_matrix/frob_matrix_base
IF (unit_nr > 0) THEN
WRITE (UNIT=unit_nr, FMT="(T3,A,A,A,T73,ES8.1)") &
"HFX_RI_INFO| Error for INV(", TRIM(name_prv), "):", error
END IF
CALL dbcsr_release(matrix_tmp)
END SUBROUTINE
! **************************************************************************************************
!> \brief ...
!> \param matrix ...
!> \param matrix_sqrt ...
!> \param matrix_sqrt_inv ...
!> \param name ...
!> \param unit_nr ...
! **************************************************************************************************
SUBROUTINE check_sqrt(matrix, matrix_sqrt, matrix_sqrt_inv, name, unit_nr)
TYPE(dbcsr_type), INTENT(INOUT) :: matrix
TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix_sqrt, matrix_sqrt_inv
CHARACTER(len=*), INTENT(IN), OPTIONAL :: name
INTEGER, INTENT(IN) :: unit_nr
CHARACTER(len=default_string_length) :: name_prv
REAL(KIND=dp) :: frob_matrix
TYPE(dbcsr_type) :: matrix_copy, matrix_tmp
IF (PRESENT(name)) THEN
name_prv = name
ELSE
CALL dbcsr_get_info(matrix, name=name_prv)
END IF
IF (PRESENT(matrix_sqrt)) THEN
CALL dbcsr_create(matrix_tmp, template=matrix)
CALL dbcsr_copy(matrix_copy, matrix_sqrt)
CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sqrt, matrix_copy, &
0.0_dp, matrix_tmp)
CALL dbcsr_add(matrix_tmp, matrix, 1.0_dp, -1.0_dp)
frob_matrix = dbcsr_frobenius_norm(matrix_tmp)
IF (unit_nr > 0) THEN
WRITE (UNIT=unit_nr, FMT="(T3,A,A,A,T73,ES8.1)") &
"HFX_RI_INFO| Error for SQRT(", TRIM(name_prv), "):", frob_matrix
END IF
CALL dbcsr_release(matrix_tmp)
CALL dbcsr_release(matrix_copy)
END IF
IF (PRESENT(matrix_sqrt_inv)) THEN
CALL dbcsr_create(matrix_tmp, template=matrix)
CALL dbcsr_copy(matrix_copy, matrix_sqrt_inv)
CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sqrt_inv, matrix_copy, &
0.0_dp, matrix_tmp)
CALL check_inverse(matrix_tmp, matrix, name="SQRT("//TRIM(name_prv)//")", unit_nr=unit_nr)
CALL dbcsr_release(matrix_tmp)
CALL dbcsr_release(matrix_copy)
END IF
END SUBROUTINE
! **************************************************************************************************
!> \brief Calculate 2-center and 3-center integrals
!>
!> 2c: K_1(P, R) = (P|v1|R) and K_2(P, R) = (P|v2|R)
!> 3c: int_3c(mu, lambda, P) = (mu lambda |v2| P)
!> v_1 is HF operator, v_2 is RI metric
!> \param qs_env ...
!> \param ri_data ...
!> \param t_2c_int_RI K_2(P, R) note: even with k-point, only need on central cell
!> \param t_2c_int_pot K_1(P, R)
!> \param t_3c_int int_3c(mu, lambda, P)
!> \param do_kpoints ...
!> \notes the integral tensor arrays are already allocated on entry
! **************************************************************************************************
SUBROUTINE hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_int_RI, t_2c_int_pot, t_3c_int, do_kpoints)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
TYPE(dbcsr_type), DIMENSION(:), INTENT(OUT) :: t_2c_int_RI, t_2c_int_pot
TYPE(dbt_type), DIMENSION(:, :) :: t_3c_int
LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints
CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ri_pre_scf_calc_tensors'
CHARACTER :: symm
INTEGER :: handle, i_img, i_mem, ibasis, j_img, &
n_mem, natom, nblks, nimg, nkind, &
nthreads
INTEGER(int_8) :: nze
INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_AO_1, dist_AO_2, dist_RI, dist_RI_ext, &
ends_array_mc_block_int, ends_array_mc_int, sizes_AO, sizes_RI, sizes_RI_ext, &
sizes_RI_ext_split, starts_array_mc_block_int, starts_array_mc_int
INTEGER, DIMENSION(3) :: pcoord, pdims
INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
LOGICAL :: converged, do_kpoints_prv
REAL(dp) :: max_ev, min_ev, occ, RI_range
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dbcsr_distribution_type) :: dbcsr_dist
TYPE(dbt_distribution_type) :: t_dist
TYPE(dbt_pgrid_type) :: pgrid
TYPE(dbt_type) :: t_3c_tmp
TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_int_batched
TYPE(dft_control_type), POINTER :: dft_control
TYPE(distribution_2d_type), POINTER :: dist_2d
TYPE(distribution_3d_type) :: dist_3d
TYPE(gto_basis_set_p_type), ALLOCATABLE, &
DIMENSION(:), TARGET :: basis_set_AO, basis_set_RI
TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
TYPE(mp_cart_type) :: mp_comm_t3c
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(neighbor_list_3c_type) :: nl_3c
TYPE(neighbor_list_set_p_type), DIMENSION(:), &
POINTER :: nl_2c_pot, nl_2c_RI
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_ks_env_type), POINTER :: ks_env
CALL timeset(routineN, handle)
NULLIFY (col_bsize, row_bsize, dist_2d, nl_2c_pot, nl_2c_RI, &
particle_set, qs_kind_set, ks_env, para_env)
CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, particle_set=particle_set, &
distribution_2d=dist_2d, ks_env=ks_env, dft_control=dft_control, para_env=para_env)
RI_range = 0.0_dp
do_kpoints_prv = .FALSE.
IF (PRESENT(do_kpoints)) do_kpoints_prv = do_kpoints
nimg = 1
IF (do_kpoints_prv) THEN
nimg = ri_data%nimg
RI_range = ri_data%kp_RI_range
END IF
ALLOCATE (sizes_RI(natom), sizes_AO(natom))
ALLOCATE (basis_set_RI(nkind), basis_set_AO(nkind))
CALL basis_set_list_setup(basis_set_RI, ri_data%ri_basis_type, qs_kind_set)
CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_RI, basis=basis_set_RI)
CALL basis_set_list_setup(basis_set_AO, ri_data%orb_basis_type, qs_kind_set)
CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_AO, basis=basis_set_AO)
DO ibasis = 1, SIZE(basis_set_AO)
orb_basis => basis_set_AO(ibasis)%gto_basis_set
ri_basis => basis_set_RI(ibasis)%gto_basis_set
! interaction radii should be based on eps_pgf_orb controlled in RI section
! (since hartree-fock needs very tight eps_pgf_orb for Kohn-Sham/Fock matrix but eps_pgf_orb
! can be much looser in RI HFX since no systematic error is introduced with tensor sparsity)
CALL init_interaction_radii_orb_basis(orb_basis, ri_data%eps_pgf_orb)
CALL init_interaction_radii_orb_basis(ri_basis, ri_data%eps_pgf_orb)
END DO
n_mem = ri_data%n_mem
CALL create_tensor_batches(sizes_RI, n_mem, starts_array_mc_int, ends_array_mc_int, &
starts_array_mc_block_int, ends_array_mc_block_int)
DEALLOCATE (starts_array_mc_int, ends_array_mc_int)
!We separate the K-points and standard 3c integrals, because handle quite differently
IF (.NOT. do_kpoints_prv) THEN
nthreads = 1
!$ nthreads = omp_get_num_threads()
pdims = 0
CALL dbt_pgrid_create(para_env, pdims, pgrid, tensor_dims=[MAX(1, natom/(n_mem*nthreads)), natom, natom])
ALLOCATE (t_3c_int_batched(1, 1))
CALL create_3c_tensor(t_3c_int_batched(1, 1), dist_RI, dist_AO_1, dist_AO_2, pgrid, &
sizes_RI, sizes_AO, sizes_AO, map1=[1], map2=[2, 3], &
name="(RI | AO AO)")
CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set, atomic_kind_set=atomic_kind_set)
CALL dbt_mp_environ_pgrid(pgrid, pdims, pcoord)
CALL mp_comm_t3c%create(pgrid%mp_comm_2d, 3, pdims)
CALL distribution_3d_create(dist_3d, dist_RI, dist_AO_1, dist_AO_2, &
nkind, particle_set, mp_comm_t3c, own_comm=.TRUE.)
DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2)
CALL create_3c_tensor(t_3c_int(1, 1), dist_RI, dist_AO_1, dist_AO_2, ri_data%pgrid, &
ri_data%bsizes_RI_split, ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
map1=[1], map2=[2, 3], &
name="O (RI AO | AO)")
CALL build_3c_neighbor_lists(nl_3c, basis_set_RI, basis_set_AO, basis_set_AO, dist_3d, ri_data%ri_metric, &
"HFX_3c_nl", qs_env, op_pos=1, sym_jk=.TRUE., own_dist=.TRUE.)
DO i_mem = 1, n_mem
CALL build_3c_integrals(t_3c_int_batched, ri_data%filter_eps/2, qs_env, nl_3c, &
basis_set_RI, basis_set_AO, basis_set_AO, &
ri_data%ri_metric, int_eps=ri_data%eps_schwarz, op_pos=1, &
desymmetrize=.FALSE., &
bounds_i=[starts_array_mc_block_int(i_mem), ends_array_mc_block_int(i_mem)])
CALL dbt_copy(t_3c_int_batched(1, 1), t_3c_int(1, 1), summation=.TRUE., move_data=.TRUE.)
CALL dbt_filter(t_3c_int(1, 1), ri_data%filter_eps/2)
END DO
CALL dbt_destroy(t_3c_int_batched(1, 1))
CALL neighbor_list_3c_destroy(nl_3c)
CALL dbt_create(t_3c_int(1, 1), t_3c_tmp)
IF (ri_data%flavor == ri_pmat) THEN ! desymmetrize
! desymmetrize
CALL dbt_copy(t_3c_int(1, 1), t_3c_tmp)
CALL dbt_copy(t_3c_tmp, t_3c_int(1, 1), order=[1, 3, 2], summation=.TRUE., move_data=.TRUE.)
! For RI-RHO filter_eps_storage is reserved for screening tensor contracted with RI-metric
! with RI metric but not to bare integral tensor
CALL dbt_filter(t_3c_int(1, 1), ri_data%filter_eps)
ELSE
CALL dbt_filter(t_3c_int(1, 1), ri_data%filter_eps_storage/2)
END IF
CALL dbt_destroy(t_3c_tmp)
ELSE !do_kpoints
nthreads = 1
!$ nthreads = omp_get_num_threads()
pdims = 0
CALL dbt_pgrid_create(para_env, pdims, pgrid, tensor_dims=[natom, natom, MAX(1, natom/(n_mem*nthreads))])
!In k-points HFX, we stack all images along the RI direction in the same tensors, in order
!to avoid storing nimg x ncell_RI different tensors (very memory intensive)
nblks = SIZE(ri_data%bsizes_RI_split)
ALLOCATE (sizes_RI_ext(natom*ri_data%ncell_RI), sizes_RI_ext_split(nblks*ri_data%ncell_RI))
DO i_img = 1, ri_data%ncell_RI
sizes_RI_ext((i_img - 1)*natom + 1:i_img*natom) = sizes_RI(:)
sizes_RI_ext_split((i_img - 1)*nblks + 1:i_img*nblks) = ri_data%bsizes_RI_split(:)
END DO
CALL create_3c_tensor(t_3c_tmp, dist_AO_1, dist_AO_2, dist_RI, &
pgrid, sizes_AO, sizes_AO, sizes_RI, map1=[1, 2], map2=[3], &
name="(AO AO | RI)")
CALL dbt_destroy(t_3c_tmp)
!For the integrals to work, the distribution along the RI direction must be repeated
ALLOCATE (dist_RI_ext(natom*ri_data%ncell_RI))
DO i_img = 1, ri_data%ncell_RI
dist_RI_ext((i_img - 1)*natom + 1:i_img*natom) = dist_RI(:)
END DO
ALLOCATE (t_3c_int_batched(nimg, 1))
CALL dbt_distribution_new(t_dist, pgrid, dist_AO_1, dist_AO_2, dist_RI_ext)
CALL dbt_create(t_3c_int_batched(1, 1), "KP_3c_ints", t_dist, [1, 2], [3], &
sizes_AO, sizes_AO, sizes_RI_ext)
DO i_img = 2, nimg
CALL dbt_create(t_3c_int_batched(1, 1), t_3c_int_batched(i_img, 1))
END DO
CALL dbt_distribution_destroy(t_dist)
CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set, atomic_kind_set=atomic_kind_set)
CALL dbt_mp_environ_pgrid(pgrid, pdims, pcoord)
CALL mp_comm_t3c%create(pgrid%mp_comm_2d, 3, pdims)
CALL distribution_3d_create(dist_3d, dist_AO_1, dist_AO_2, dist_RI, &
nkind, particle_set, mp_comm_t3c, own_comm=.TRUE.)
DEALLOCATE (dist_RI, dist_AO_1, dist_AO_2)
! create 3c tensor for storage of ints
CALL build_3c_neighbor_lists(nl_3c, basis_set_AO, basis_set_AO, basis_set_RI, dist_3d, ri_data%ri_metric, &
"HFX_3c_nl", qs_env, op_pos=2, sym_ij=.FALSE., own_dist=.TRUE.)
CALL create_3c_tensor(t_3c_int(1, 1), dist_RI, dist_AO_1, dist_AO_2, ri_data%pgrid, &
sizes_RI_ext_split, ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
map1=[1], map2=[2, 3], &
name="O (RI AO | AO)")
DO j_img = 2, nimg
CALL dbt_create(t_3c_int(1, 1), t_3c_int(1, j_img))
END DO
CALL dbt_create(t_3c_int(1, 1), t_3c_tmp)
DO i_mem = 1, n_mem
CALL build_3c_integrals(t_3c_int_batched, ri_data%filter_eps, qs_env, nl_3c, &
basis_set_AO, basis_set_AO, basis_set_RI, &
ri_data%ri_metric, int_eps=ri_data%eps_schwarz, op_pos=2, &
desymmetrize=.FALSE., do_kpoints=.TRUE., do_hfx_kpoints=.TRUE., &
bounds_k=[starts_array_mc_block_int(i_mem), ends_array_mc_block_int(i_mem)], &
RI_range=RI_range, img_to_RI_cell=ri_data%img_to_RI_cell)
DO i_img = 1, nimg
!we start with (mu^0 sigma^i | P^j) and finish with (P^i | mu^0 sigma^j)
CALL get_tensor_occupancy(t_3c_int_batched(i_img, 1), nze, occ)
IF (nze > 0) THEN
CALL dbt_copy(t_3c_int_batched(i_img, 1), t_3c_tmp, order=[3, 2, 1], move_data=.TRUE.)
CALL dbt_filter(t_3c_tmp, ri_data%filter_eps)
CALL dbt_copy(t_3c_tmp, t_3c_int(1, i_img), order=[1, 3, 2], &
summation=.TRUE., move_data=.TRUE.)
END IF
END DO
END DO
DO i_img = 1, nimg
CALL dbt_destroy(t_3c_int_batched(i_img, 1))
END DO
DEALLOCATE (t_3c_int_batched)
CALL neighbor_list_3c_destroy(nl_3c)
CALL dbt_destroy(t_3c_tmp)
END IF
CALL dbt_pgrid_destroy(pgrid)
CALL build_2c_neighbor_lists(nl_2c_pot, basis_set_RI, basis_set_RI, ri_data%hfx_pot, &
"HFX_2c_nl_pot", qs_env, sym_ij=.NOT. do_kpoints_prv, &
dist_2d=dist_2d)
CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
ALLOCATE (row_bsize(SIZE(sizes_RI)))
ALLOCATE (col_bsize(SIZE(sizes_RI)))
row_bsize(:) = sizes_RI
col_bsize(:) = sizes_RI
!Use non-symmetric nl and matrices for k-points
symm = dbcsr_type_symmetric
IF (do_kpoints_prv) symm = dbcsr_type_no_symmetry
CALL dbcsr_create(t_2c_int_pot(1), "(R|P) HFX", dbcsr_dist, symm, row_bsize, col_bsize)
DO i_img = 2, nimg
CALL dbcsr_create(t_2c_int_pot(i_img), template=t_2c_int_pot(1))
END DO
IF (.NOT. ri_data%same_op) THEN
CALL dbcsr_create(t_2c_int_RI(1), "(R|P) HFX", dbcsr_dist, symm, row_bsize, col_bsize)
DO i_img = 2, nimg
CALL dbcsr_create(t_2c_int_RI(i_img), template=t_2c_int_RI(1))
END DO
END IF
DEALLOCATE (col_bsize, row_bsize)
CALL dbcsr_distribution_release(dbcsr_dist)
CALL build_2c_integrals(t_2c_int_pot, ri_data%filter_eps_2c, qs_env, nl_2c_pot, basis_set_RI, basis_set_RI, &
ri_data%hfx_pot, do_kpoints=do_kpoints_prv, do_hfx_kpoints=do_kpoints_prv)
CALL release_neighbor_list_sets(nl_2c_pot)
IF (.NOT. ri_data%same_op) THEN
CALL build_2c_neighbor_lists(nl_2c_RI, basis_set_RI, basis_set_RI, ri_data%ri_metric, &
"HFX_2c_nl_RI", qs_env, sym_ij=.NOT. do_kpoints_prv, &
dist_2d=dist_2d)
CALL build_2c_integrals(t_2c_int_RI, ri_data%filter_eps_2c, qs_env, nl_2c_RI, basis_set_RI, basis_set_RI, &
ri_data%ri_metric, do_kpoints=do_kpoints_prv, do_hfx_kpoints=do_kpoints_prv)
CALL release_neighbor_list_sets(nl_2c_RI)
END IF
DO ibasis = 1, SIZE(basis_set_AO)
orb_basis => basis_set_AO(ibasis)%gto_basis_set
ri_basis => basis_set_RI(ibasis)%gto_basis_set
! reset interaction radii of orb basis
CALL init_interaction_radii_orb_basis(orb_basis, dft_control%qs_control%eps_pgf_orb)
CALL init_interaction_radii_orb_basis(ri_basis, dft_control%qs_control%eps_pgf_orb)
END DO
IF (ri_data%calc_condnum) THEN
CALL arnoldi_extremal(t_2c_int_pot(1), max_ev, min_ev, threshold=ri_data%eps_lanczos, &
max_iter=ri_data%max_iter_lanczos, converged=converged)
IF (.NOT. converged) THEN
CPWARN("Not converged: unreliable condition number estimate of (P|Q) matrix (HFX potential).")
END IF
IF (ri_data%unit_nr > 0) THEN
WRITE (ri_data%unit_nr, '(T2,A)') "2-Norm Condition Number of (P|Q) integrals (HFX potential)"
IF (min_ev > 0) THEN
WRITE (ri_data%unit_nr, '(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
"CN : max/min ev: ", max_ev, " / ", min_ev, "=", max_ev/min_ev, "Log(2-CN):", LOG10(max_ev/min_ev)
ELSE
WRITE (ri_data%unit_nr, '(T4,A,ES11.3E3,T32,A,ES11.3E3,T63,A)') &
"CN : max/min ev: ", max_ev, " / ", min_ev, "Log(CN): infinity"
END IF
END IF
IF (.NOT. ri_data%same_op) THEN
CALL arnoldi_extremal(t_2c_int_RI(1), max_ev, min_ev, threshold=ri_data%eps_lanczos, &
max_iter=ri_data%max_iter_lanczos, converged=converged)
IF (.NOT. converged) THEN
CPWARN("Not converged: unreliable condition number estimate of (P|Q) matrix (RI metric).")
END IF
IF (ri_data%unit_nr > 0) THEN
WRITE (ri_data%unit_nr, '(T2,A)') "2-Norm Condition Number of (P|Q) integrals (RI metric)"
IF (min_ev > 0) THEN
WRITE (ri_data%unit_nr, '(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
"CN : max/min ev: ", max_ev, " / ", min_ev, "=", max_ev/min_ev, "Log(2-CN):", LOG10(max_ev/min_ev)
ELSE
WRITE (ri_data%unit_nr, '(T4,A,ES11.3E3,T32,A,ES11.3E3,T63,A)') &
"CN : max/min ev: ", max_ev, " / ", min_ev, "Log(CN): infinity"
END IF
END IF
END IF
END IF
CALL timestop(handle)
END SUBROUTINE
! **************************************************************************************************
!> \brief Pre-SCF steps in rho flavor of RI HFX
!>
!> K(P, S) = sum_{R,Q} K_2(P, R)^{-1} K_1(R, Q) K_2(Q, S)^{-1}
!> Calculate B(mu, lambda, R) = sum_P int_3c(mu, lambda, P) K(P, R)
!> \param qs_env ...
!> \param ri_data ...
! **************************************************************************************************
SUBROUTINE hfx_ri_pre_scf_Pmat(qs_env, ri_data)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ri_pre_scf_Pmat'
INTEGER :: handle, handle2, i_mem, j_mem, &
n_dependent, unit_nr, unit_nr_dbcsr
INTEGER(int_8) :: nflop, nze, nze_O
INTEGER, ALLOCATABLE, DIMENSION(:) :: batch_ranges_AO, batch_ranges_RI
INTEGER, DIMENSION(2, 1) :: bounds_i
INTEGER, DIMENSION(2, 2) :: bounds_j
INTEGER, DIMENSION(3) :: dims_3c
REAL(KIND=dp) :: compression_factor, memory_3c, occ, &
threshold
TYPE(cp_blacs_env_type), POINTER :: blacs_env
TYPE(dbcsr_type), DIMENSION(1) :: t_2c_int_mat, t_2c_op_pot, t_2c_op_RI, &
t_2c_tmp, t_2c_tmp_2
TYPE(dbt_type) :: t_3c_2
TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_2c_int, t_2c_work
TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_int_1
TYPE(mp_para_env_type), POINTER :: para_env
CALL timeset(routineN, handle)
unit_nr_dbcsr = ri_data%unit_nr_dbcsr
unit_nr = ri_data%unit_nr
CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
CALL timeset(routineN//"_int", handle2)
ALLOCATE (t_2c_int(1), t_2c_work(1), t_3c_int_1(1, 1))
CALL hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_op_RI, t_2c_op_pot, t_3c_int_1)
CALL dbt_copy(t_3c_int_1(1, 1), ri_data%t_3c_int_ctr_3(1, 1), order=[1, 2, 3], move_data=.TRUE.)
CALL dbt_destroy(t_3c_int_1(1, 1))
CALL timestop(handle2)
CALL timeset(routineN//"_2c", handle2)
IF (ri_data%same_op) t_2c_op_RI(1) = t_2c_op_pot(1)
CALL dbcsr_create(t_2c_int_mat(1), template=t_2c_op_RI(1), matrix_type=dbcsr_type_no_symmetry)
threshold = MAX(ri_data%filter_eps, 1.0e-12_dp)
SELECT CASE (ri_data%t2c_method)
CASE (hfx_ri_do_2c_iter)
CALL invert_hotelling(t_2c_int_mat(1), t_2c_op_RI(1), &
threshold=threshold, silent=.FALSE.)
CASE (hfx_ri_do_2c_cholesky)
CALL dbcsr_copy(t_2c_int_mat(1), t_2c_op_RI(1))
CALL cp_dbcsr_cholesky_decompose(t_2c_int_mat(1), para_env=para_env, blacs_env=blacs_env)
CALL cp_dbcsr_cholesky_invert(t_2c_int_mat(1), para_env=para_env, blacs_env=blacs_env, upper_to_full=.TRUE.)
CASE (hfx_ri_do_2c_diag)
CALL dbcsr_copy(t_2c_int_mat(1), t_2c_op_RI(1))
CALL cp_dbcsr_power(t_2c_int_mat(1), -1.0_dp, ri_data%eps_eigval, n_dependent, &
para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
END SELECT
IF (ri_data%check_2c_inv) THEN
CALL check_inverse(t_2c_int_mat(1), t_2c_op_RI(1), unit_nr=unit_nr)
END IF
!Need to save the (P|Q)^-1 tensor for forces (inverse metric if not same_op)
CALL dbt_create(t_2c_int_mat(1), t_2c_work(1))
CALL dbt_copy_matrix_to_tensor(t_2c_int_mat(1), t_2c_work(1))
CALL dbt_copy(t_2c_work(1), ri_data%t_2c_inv(1, 1), move_data=.TRUE.)
CALL dbt_destroy(t_2c_work(1))
CALL dbt_filter(ri_data%t_2c_inv(1, 1), ri_data%filter_eps)
IF (.NOT. ri_data%same_op) THEN
!Also save the RI (P|Q) integral
CALL dbt_create(t_2c_op_pot(1), t_2c_work(1))
CALL dbt_copy_matrix_to_tensor(t_2c_op_pot(1), t_2c_work(1))
CALL dbt_copy(t_2c_work(1), ri_data%t_2c_pot(1, 1), move_data=.TRUE.)
CALL dbt_destroy(t_2c_work(1))
CALL dbt_filter(ri_data%t_2c_pot(1, 1), ri_data%filter_eps)
END IF
IF (ri_data%same_op) THEN
CALL dbcsr_release(t_2c_op_pot(1))
ELSE
CALL dbcsr_create(t_2c_tmp(1), template=t_2c_op_RI(1), matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_create(t_2c_tmp_2(1), template=t_2c_op_RI(1), matrix_type=dbcsr_type_no_symmetry)
CALL dbcsr_release(t_2c_op_RI(1))
CALL dbcsr_multiply('N', 'N', 1.0_dp, t_2c_int_mat(1), t_2c_op_pot(1), 0.0_dp, t_2c_tmp(1), &
filter_eps=ri_data%filter_eps)
CALL dbcsr_release(t_2c_op_pot(1))
CALL dbcsr_multiply('N', 'N', 1.0_dp, t_2c_tmp(1), t_2c_int_mat(1), 0.0_dp, t_2c_tmp_2(1), &
filter_eps=ri_data%filter_eps)
CALL dbcsr_release(t_2c_tmp(1))
CALL dbcsr_release(t_2c_int_mat(1))
t_2c_int_mat(1) = t_2c_tmp_2(1)
END IF
CALL dbt_create(t_2c_int_mat(1), t_2c_int(1), name="(RI|RI)")
CALL dbt_copy_matrix_to_tensor(t_2c_int_mat(1), t_2c_int(1))
CALL dbcsr_release(t_2c_int_mat(1))
CALL dbt_copy(t_2c_int(1), ri_data%t_2c_int(1, 1), move_data=.TRUE.)
CALL dbt_destroy(t_2c_int(1))
CALL dbt_filter(ri_data%t_2c_int(1, 1), ri_data%filter_eps)
CALL timestop(handle2)
CALL dbt_create(ri_data%t_3c_int_ctr_3(1, 1), t_3c_2)
CALL dbt_get_info(ri_data%t_3c_int_ctr_3(1, 1), nfull_total=dims_3c)
memory_3c = 0.0_dp
nze_O = 0
ALLOCATE (batch_ranges_RI(ri_data%n_mem_RI + 1))
ALLOCATE (batch_ranges_AO(ri_data%n_mem + 1))
batch_ranges_RI(:ri_data%n_mem_RI) = ri_data%starts_array_RI_mem_block(:)
batch_ranges_RI(ri_data%n_mem_RI + 1) = ri_data%ends_array_RI_mem_block(ri_data%n_mem_RI) + 1
batch_ranges_AO(:ri_data%n_mem) = ri_data%starts_array_mem_block(:)
batch_ranges_AO(ri_data%n_mem + 1) = ri_data%ends_array_mem_block(ri_data%n_mem) + 1
CALL dbt_batched_contract_init(ri_data%t_3c_int_ctr_3(1, 1), batch_range_1=batch_ranges_RI, &
batch_range_2=batch_ranges_AO)
CALL dbt_batched_contract_init(t_3c_2, batch_range_1=batch_ranges_RI, &
batch_range_2=batch_ranges_AO)
DO i_mem = 1, ri_data%n_mem_RI
bounds_i(:, 1) = [ri_data%starts_array_RI_mem(i_mem), ri_data%ends_array_RI_mem(i_mem)]
CALL dbt_batched_contract_init(ri_data%t_2c_int(1, 1))
DO j_mem = 1, ri_data%n_mem
bounds_j(:, 1) = [ri_data%starts_array_mem(j_mem), ri_data%ends_array_mem(j_mem)]
bounds_j(:, 2) = [1, dims_3c(3)]
CALL timeset(routineN//"_RIx3C", handle2)
CALL dbt_contract(1.0_dp, ri_data%t_2c_int(1, 1), ri_data%t_3c_int_ctr_3(1, 1), &
0.0_dp, t_3c_2, &
contract_1=[2], notcontract_1=[1], &
contract_2=[1], notcontract_2=[2, 3], &
map_1=[1], map_2=[2, 3], filter_eps=ri_data%filter_eps_storage, &
bounds_2=bounds_i, &
bounds_3=bounds_j, &
unit_nr=unit_nr_dbcsr, &
flop=nflop)
ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
CALL timestop(handle2)
CALL timeset(routineN//"_copy_2", handle2)
CALL dbt_copy(t_3c_2, ri_data%t_3c_int_ctr_1(1, 1), order=[2, 1, 3], move_data=.TRUE.)
CALL get_tensor_occupancy(ri_data%t_3c_int_ctr_1(1, 1), nze, occ)
nze_O = nze_O + nze
CALL compress_tensor(ri_data%t_3c_int_ctr_1(1, 1), ri_data%blk_indices(j_mem, i_mem)%ind, &
ri_data%store_3c(j_mem, i_mem), ri_data%filter_eps_storage, memory_3c)
CALL timestop(handle2)
END DO
CALL dbt_batched_contract_finalize(ri_data%t_2c_int(1, 1))
END DO
CALL dbt_batched_contract_finalize(t_3c_2)
CALL dbt_batched_contract_finalize(ri_data%t_3c_int_ctr_3(1, 1))
CALL para_env%sum(memory_3c)
compression_factor = REAL(nze_O, dp)*1.0E-06*8.0_dp/memory_3c
IF (unit_nr > 0) THEN
WRITE (UNIT=unit_nr, FMT="((T3,A,T66,F11.2,A4))") &
"MEMORY_INFO| Memory for 3-center HFX integrals (compressed):", memory_3c, ' MiB'
WRITE (UNIT=unit_nr, FMT="((T3,A,T60,F21.2))") &
"MEMORY_INFO| Compression factor: ", compression_factor
END IF
CALL dbt_clear(ri_data%t_2c_int(1, 1))
CALL dbt_destroy(t_3c_2)
CALL dbt_copy(ri_data%t_3c_int_ctr_3(1, 1), ri_data%t_3c_int_ctr_2(1, 1), order=[2, 1, 3], move_data=.TRUE.)
CALL timestop(handle)
END SUBROUTINE
! **************************************************************************************************
!> \brief Sorts 2d indices w.r.t. rows and columns
!> \param blk_ind ...
! **************************************************************************************************
SUBROUTINE sort_unique_blkind_2d(blk_ind)
INTEGER, ALLOCATABLE, DIMENSION(:, :), &
INTENT(INOUT) :: blk_ind
INTEGER :: end_ind, iblk, iblk_all, irow, nblk, &
ncols, start_ind
INTEGER, ALLOCATABLE, DIMENSION(:) :: ind_1, ind_2, sort_1, sort_2
INTEGER, ALLOCATABLE, DIMENSION(:, :) :: blk_ind_tmp
nblk = SIZE(blk_ind, 1)
ALLOCATE (sort_1(nblk))
ALLOCATE (ind_1(nblk))
sort_1(:) = blk_ind(:, 1)
CALL sort(sort_1, nblk, ind_1)
blk_ind(:, :) = blk_ind(ind_1, :)
start_ind = 1
DO WHILE (start_ind <= nblk)
irow = blk_ind(start_ind, 1)
end_ind = start_ind
IF (end_ind + 1 <= nblk) THEN
DO WHILE (blk_ind(end_ind + 1, 1) == irow)
end_ind = end_ind + 1
IF (end_ind + 1 > nblk) EXIT
END DO
END IF
ncols = end_ind - start_ind + 1
ALLOCATE (sort_2(ncols))
ALLOCATE (ind_2(ncols))
sort_2(:) = blk_ind(start_ind:end_ind, 2)
CALL sort(sort_2, ncols, ind_2)
ind_2 = ind_2 + start_ind - 1
blk_ind(start_ind:end_ind, :) = blk_ind(ind_2, :)
start_ind = end_ind + 1
DEALLOCATE (sort_2, ind_2)
END DO
ALLOCATE (blk_ind_tmp(nblk, 2))
blk_ind_tmp = 0
iblk = 0