(git:374b731)
Loading...
Searching...
No Matches
hfx_ri.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief RI-methods for HFX
10! **************************************************************************************************
11
12MODULE hfx_ri
13
14 USE omp_lib, ONLY: omp_get_num_threads,&
15 omp_get_thread_num
16 USE arnoldi_api, ONLY: arnoldi_extremal
17 USE atomic_kind_types, ONLY: atomic_kind_type,&
18 get_atomic_kind_set
19 USE basis_set_types, ONLY: gto_basis_set_p_type,&
20 gto_basis_set_type
21 USE cell_types, ONLY: cell_type
22 USE cp_blacs_env, ONLY: cp_blacs_env_type
23 USE cp_control_types, ONLY: dft_control_type
24 USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
25 cp_dbcsr_cholesky_invert
26 USE cp_dbcsr_diag, ONLY: cp_dbcsr_power
27 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 copy_fm_to_dbcsr,&
29 cp_dbcsr_dist2d_to_dist,&
30 dbcsr_deallocate_matrix_set
31 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
32 cp_fm_struct_release,&
33 cp_fm_struct_type
34 USE cp_fm_types, ONLY: cp_fm_create,&
35 cp_fm_release,&
36 cp_fm_type,&
37 cp_fm_write_unformatted
38 USE cp_log_handling, ONLY: cp_get_default_logger,&
39 cp_logger_type
40 USE cp_output_handling, ONLY: cp_p_file,&
41 cp_print_key_finished_output,&
42 cp_print_key_should_output,&
43 cp_print_key_unit_nr
44 USE dbcsr_api, ONLY: &
45 dbcsr_add, dbcsr_add_on_diag, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, &
46 dbcsr_desymmetrize, dbcsr_distribution_get, dbcsr_distribution_release, &
47 dbcsr_distribution_type, dbcsr_dot, dbcsr_filter, dbcsr_frobenius_norm, dbcsr_get_info, &
48 dbcsr_get_num_blocks, dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, &
49 dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
50 USE dbt_api, ONLY: &
51 dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_clear, dbt_contract, &
52 dbt_copy, dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, &
53 dbt_default_distvec, dbt_destroy, dbt_distribution_destroy, dbt_distribution_new, &
54 dbt_distribution_type, dbt_filter, dbt_get_block, dbt_get_info, dbt_get_num_blocks_total, &
55 dbt_iterator_blocks_left, dbt_iterator_next_block, dbt_iterator_start, dbt_iterator_stop, &
56 dbt_iterator_type, dbt_mp_environ_pgrid, dbt_nd_mp_comm, dbt_pgrid_create, &
57 dbt_pgrid_destroy, dbt_pgrid_type, dbt_put_block, dbt_reserve_blocks, dbt_scale, dbt_type
58 USE distribution_2d_types, ONLY: distribution_2d_type
59 USE hfx_types, ONLY: alloc_containers,&
60 block_ind_type,&
61 dealloc_containers,&
62 hfx_compression_type,&
63 hfx_ri_init,&
64 hfx_ri_release,&
65 hfx_ri_type
66 USE input_constants, ONLY: hfx_ri_do_2c_cholesky,&
67 hfx_ri_do_2c_diag,&
68 hfx_ri_do_2c_iter
69 USE input_cp2k_hfx, ONLY: ri_mo,&
70 ri_pmat
71 USE input_section_types, ONLY: section_vals_get_subs_vals,&
72 section_vals_type,&
73 section_vals_val_get
74 USE iterate_matrix, ONLY: invert_hotelling,&
75 matrix_sqrt_newton_schulz
76 USE kinds, ONLY: default_string_length,&
77 dp,&
78 int_8
79 USE machine, ONLY: m_walltime
80 USE message_passing, ONLY: mp_cart_type,&
81 mp_comm_type,&
82 mp_para_env_type
83 USE particle_methods, ONLY: get_particle_set
84 USE particle_types, ONLY: particle_type
85 USE qs_environment_types, ONLY: get_qs_env,&
86 qs_environment_type
87 USE qs_force_types, ONLY: qs_force_type
88 USE qs_integral_utils, ONLY: basis_set_list_setup
89 USE qs_interactions, ONLY: init_interaction_radii_orb_basis
90 USE qs_kind_types, ONLY: qs_kind_type
91 USE qs_ks_types, ONLY: qs_ks_env_type
92 USE qs_mo_types, ONLY: get_mo_set,&
93 mo_set_type
94 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
95 release_neighbor_list_sets
96 USE qs_rho_types, ONLY: qs_rho_get,&
97 qs_rho_type
98 USE qs_tensors, ONLY: &
99 build_2c_derivatives, build_2c_integrals, build_2c_neighbor_lists, build_3c_derivatives, &
100 build_3c_integrals, build_3c_neighbor_lists, calc_2c_virial, calc_3c_virial, &
101 compress_tensor, decompress_tensor, get_tensor_occupancy, neighbor_list_3c_destroy
102 USE qs_tensors_types, ONLY: create_2c_tensor,&
103 create_3c_tensor,&
104 create_tensor_batches,&
105 distribution_3d_create,&
106 distribution_3d_type,&
107 neighbor_list_3c_type,&
108 split_block_sizes
109 USE util, ONLY: sort
110 USE virial_types, ONLY: virial_type
111#include "./base/base_uses.f90"
112
113!$ USE OMP_LIB, ONLY: omp_get_num_threads
114
115 IMPLICIT NONE
116 PRIVATE
117
118 PUBLIC :: hfx_ri_update_ks, hfx_ri_update_forces, get_force_from_3c_trace, get_2c_der_force, &
119 get_idx_to_atom, print_ri_hfx, hfx_ri_pre_scf_calc_tensors, check_inverse
120
121 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_ri'
122CONTAINS
123
124! **************************************************************************************************
125!> \brief Switches the RI_FLAVOR from MO to RHO or vice-versa
126!> \param ri_data ...
127!> \param qs_env ...
128!> \note As a side product, the ri_data is mostly reinitialized and the integrals recomputed
129! **************************************************************************************************
130 SUBROUTINE switch_ri_flavor(ri_data, qs_env)
131 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
132 TYPE(qs_environment_type), POINTER :: qs_env
133
134 CHARACTER(LEN=*), PARAMETER :: routineN = 'switch_ri_flavor'
135
136 INTEGER :: handle, n_mem, new_flavor
137 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
138 TYPE(dft_control_type), POINTER :: dft_control
139 TYPE(mp_para_env_type), POINTER :: para_env
140 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
141 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
142
143 NULLIFY (qs_kind_set, particle_set, atomic_kind_set, para_env, dft_control)
144
145 CALL timeset(routinen, handle)
146
147 CALL get_qs_env(qs_env, para_env=para_env, dft_control=dft_control, atomic_kind_set=atomic_kind_set, &
148 particle_set=particle_set, qs_kind_set=qs_kind_set)
149
150 CALL hfx_ri_release(ri_data, write_stats=.false.)
151
152 IF (ri_data%flavor == ri_pmat) THEN
153 new_flavor = ri_mo
154 ELSE
155 new_flavor = ri_pmat
156 END IF
157 ri_data%flavor = new_flavor
158
159 n_mem = ri_data%n_mem_input
160 ri_data%n_mem_input = ri_data%n_mem_flavor_switch
161 ri_data%n_mem_flavor_switch = n_mem
162
163 CALL hfx_ri_init(ri_data, qs_kind_set, particle_set, atomic_kind_set, para_env)
164
165 !Need to recalculate the integrals in the new flavor
166 !TODO: should we backup the integrals and symmetrize/desymmetrize them instead of recomputing ?!?
167 ! that only makes sense if actual integral calculation is not negligible
168 IF (ri_data%flavor == ri_pmat) THEN
169 CALL hfx_ri_pre_scf_pmat(qs_env, ri_data)
170 ELSE
171 CALL hfx_ri_pre_scf_mo(qs_env, ri_data, dft_control%nspins)
172 END IF
173
174 IF (ri_data%unit_nr > 0) THEN
175 IF (ri_data%flavor == ri_pmat) THEN
176 WRITE (ri_data%unit_nr, '(T2,A)') "HFX_RI_INFO| temporarily switched to RI_FLAVOR RHO"
177 ELSE
178 WRITE (ri_data%unit_nr, '(T2,A)') "HFX_RI_INFO| temporarily switched to RI_FLAVOR MO"
179 END IF
180 END IF
181
182 CALL timestop(handle)
183
184 END SUBROUTINE switch_ri_flavor
185
186! **************************************************************************************************
187!> \brief Pre-SCF steps in MO flavor of RI HFX
188!>
189!> Calculate 2-center & 3-center integrals (see hfx_ri_pre_scf_calc_tensors) and contract
190!> K(P, S) = sum_R K_2(P, R)^{-1} K_1(R, S)^{1/2}
191!> B(mu, lambda, R) = sum_P int_3c(mu, lambda, P) K(P, R)
192!> \param qs_env ...
193!> \param ri_data ...
194!> \param nspins ...
195! **************************************************************************************************
196 SUBROUTINE hfx_ri_pre_scf_mo(qs_env, ri_data, nspins)
197 TYPE(qs_environment_type), POINTER :: qs_env
198 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
199 INTEGER, INTENT(IN) :: nspins
200
201 CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_ri_pre_scf_mo'
202
203 INTEGER :: handle, handle2, ispin, n_dependent, &
204 unit_nr, unit_nr_dbcsr
205 REAL(KIND=dp) :: threshold
206 TYPE(cp_blacs_env_type), POINTER :: blacs_env
207 TYPE(dbcsr_type), DIMENSION(1) :: dbcsr_work_1, dbcsr_work_2, t_2c_int_mat, t_2c_op_pot, &
208 t_2c_op_pot_sqrt, t_2c_op_pot_sqrt_inv, t_2c_op_RI, t_2c_op_RI_inv
209 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_2c_int, t_2c_work
210 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_int
211 TYPE(mp_para_env_type), POINTER :: para_env
212
213 CALL timeset(routinen, handle)
214
215 unit_nr_dbcsr = ri_data%unit_nr_dbcsr
216 unit_nr = ri_data%unit_nr
217
218 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
219
220 CALL timeset(routinen//"_int", handle2)
221
222 ALLOCATE (t_2c_int(1), t_2c_work(1), t_3c_int(1, 1))
223 CALL hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_op_ri, t_2c_op_pot, t_3c_int)
224
225 CALL timestop(handle2)
226
227 CALL timeset(routinen//"_2c", handle2)
228 IF (.NOT. ri_data%same_op) THEN
229 SELECT CASE (ri_data%t2c_method)
230 CASE (hfx_ri_do_2c_iter)
231 CALL dbcsr_create(t_2c_op_ri_inv(1), template=t_2c_op_ri(1), matrix_type=dbcsr_type_no_symmetry)
232 threshold = max(ri_data%filter_eps, 1.0e-12_dp)
233 CALL invert_hotelling(t_2c_op_ri_inv(1), t_2c_op_ri(1), threshold=threshold, silent=.false.)
234 CASE (hfx_ri_do_2c_cholesky)
235 CALL dbcsr_copy(t_2c_op_ri_inv(1), t_2c_op_ri(1))
236 CALL cp_dbcsr_cholesky_decompose(t_2c_op_ri_inv(1), para_env=para_env, blacs_env=blacs_env)
237 CALL cp_dbcsr_cholesky_invert(t_2c_op_ri_inv(1), para_env=para_env, blacs_env=blacs_env, upper_to_full=.true.)
238 CASE (hfx_ri_do_2c_diag)
239 CALL dbcsr_copy(t_2c_op_ri_inv(1), t_2c_op_ri(1))
240 CALL cp_dbcsr_power(t_2c_op_ri_inv(1), -1.0_dp, ri_data%eps_eigval, n_dependent, &
241 para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
242 END SELECT
243
244 IF (ri_data%check_2c_inv) THEN
245 CALL check_inverse(t_2c_op_ri_inv(1), t_2c_op_ri(1), unit_nr=unit_nr)
246 END IF
247
248 CALL dbcsr_release(t_2c_op_ri(1))
249
250 SELECT CASE (ri_data%t2c_method)
251 CASE (hfx_ri_do_2c_iter)
252 CALL dbcsr_create(t_2c_op_pot_sqrt(1), template=t_2c_op_pot(1), matrix_type=dbcsr_type_symmetric)
253 CALL dbcsr_create(t_2c_op_pot_sqrt_inv(1), template=t_2c_op_pot(1), matrix_type=dbcsr_type_symmetric)
254 CALL matrix_sqrt_newton_schulz(t_2c_op_pot_sqrt(1), t_2c_op_pot_sqrt_inv(1), t_2c_op_pot(1), &
255 ri_data%filter_eps, ri_data%t2c_sqrt_order, ri_data%eps_lanczos, &
256 ri_data%max_iter_lanczos)
257
258 CALL dbcsr_release(t_2c_op_pot_sqrt_inv(1))
259 CASE (hfx_ri_do_2c_diag, hfx_ri_do_2c_cholesky)
260 CALL dbcsr_copy(t_2c_op_pot_sqrt(1), t_2c_op_pot(1))
261 CALL cp_dbcsr_power(t_2c_op_pot_sqrt(1), 0.5_dp, ri_data%eps_eigval, n_dependent, &
262 para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
263 END SELECT
264
265 !We need S^-1 and (P|Q) for the forces.
266 CALL dbt_create(t_2c_op_ri_inv(1), t_2c_work(1))
267 CALL dbt_copy_matrix_to_tensor(t_2c_op_ri_inv(1), t_2c_work(1))
268 CALL dbt_copy(t_2c_work(1), ri_data%t_2c_inv(1, 1), move_data=.true.)
269 CALL dbt_destroy(t_2c_work(1))
270 CALL dbt_filter(ri_data%t_2c_inv(1, 1), ri_data%filter_eps)
271
272 CALL dbt_create(t_2c_op_pot(1), t_2c_work(1))
273 CALL dbt_copy_matrix_to_tensor(t_2c_op_pot(1), t_2c_work(1))
274 CALL dbt_copy(t_2c_work(1), ri_data%t_2c_pot(1, 1), move_data=.true.)
275 CALL dbt_destroy(t_2c_work(1))
276 CALL dbt_filter(ri_data%t_2c_pot(1, 1), ri_data%filter_eps)
277
278 IF (ri_data%check_2c_inv) THEN
279 CALL check_sqrt(t_2c_op_pot(1), matrix_sqrt=t_2c_op_pot_sqrt(1), unit_nr=unit_nr)
280 END IF
281 CALL dbcsr_create(t_2c_int_mat(1), template=t_2c_op_pot(1), matrix_type=dbcsr_type_no_symmetry)
282 CALL dbcsr_multiply("N", "N", 1.0_dp, t_2c_op_ri_inv(1), t_2c_op_pot_sqrt(1), &
283 0.0_dp, t_2c_int_mat(1), filter_eps=ri_data%filter_eps)
284 CALL dbcsr_release(t_2c_op_ri_inv(1))
285 CALL dbcsr_release(t_2c_op_pot_sqrt(1))
286 ELSE
287 SELECT CASE (ri_data%t2c_method)
288 CASE (hfx_ri_do_2c_iter)
289 CALL dbcsr_create(t_2c_int_mat(1), template=t_2c_op_pot(1), matrix_type=dbcsr_type_symmetric)
290 CALL dbcsr_create(t_2c_op_pot_sqrt(1), template=t_2c_op_pot(1), matrix_type=dbcsr_type_symmetric)
291 CALL matrix_sqrt_newton_schulz(t_2c_op_pot_sqrt(1), t_2c_int_mat(1), t_2c_op_pot(1), &
292 ri_data%filter_eps, ri_data%t2c_sqrt_order, ri_data%eps_lanczos, &
293 ri_data%max_iter_lanczos)
294 CALL dbcsr_release(t_2c_op_pot_sqrt(1))
295 CASE (hfx_ri_do_2c_diag, hfx_ri_do_2c_cholesky)
296 CALL dbcsr_copy(t_2c_int_mat(1), t_2c_op_pot(1))
297 CALL cp_dbcsr_power(t_2c_int_mat(1), -0.5_dp, ri_data%eps_eigval, n_dependent, &
298 para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
299 END SELECT
300 IF (ri_data%check_2c_inv) THEN
301 CALL check_sqrt(t_2c_op_pot(1), matrix_sqrt_inv=t_2c_int_mat(1), unit_nr=unit_nr)
302 END IF
303
304 !We need (P|Q)^-1 for the forces
305 CALL dbcsr_copy(dbcsr_work_1(1), t_2c_int_mat(1))
306 CALL dbcsr_create(dbcsr_work_2(1), template=t_2c_int_mat(1))
307 CALL dbcsr_multiply("N", "N", 1.0_dp, dbcsr_work_1(1), t_2c_int_mat(1), 0.0_dp, dbcsr_work_2(1))
308 CALL dbcsr_release(dbcsr_work_1(1))
309 CALL dbt_create(dbcsr_work_2(1), t_2c_work(1))
310 CALL dbt_copy_matrix_to_tensor(dbcsr_work_2(1), t_2c_work(1))
311 CALL dbcsr_release(dbcsr_work_2(1))
312 CALL dbt_copy(t_2c_work(1), ri_data%t_2c_inv(1, 1), move_data=.true.)
313 CALL dbt_destroy(t_2c_work(1))
314 CALL dbt_filter(ri_data%t_2c_inv(1, 1), ri_data%filter_eps)
315 END IF
316
317 CALL dbcsr_release(t_2c_op_pot(1))
318
319 CALL dbt_create(t_2c_int_mat(1), t_2c_int(1), name="(RI|RI)")
320 CALL dbt_copy_matrix_to_tensor(t_2c_int_mat(1), t_2c_int(1))
321 CALL dbcsr_release(t_2c_int_mat(1))
322 DO ispin = 1, nspins
323 CALL dbt_copy(t_2c_int(1), ri_data%t_2c_int(ispin, 1))
324 END DO
325 CALL dbt_destroy(t_2c_int(1))
326 CALL timestop(handle2)
327
328 CALL timeset(routinen//"_3c", handle2)
329 CALL dbt_copy(t_3c_int(1, 1), ri_data%t_3c_int_ctr_1(1, 1), order=[2, 1, 3], move_data=.true.)
330 CALL dbt_filter(ri_data%t_3c_int_ctr_1(1, 1), ri_data%filter_eps)
331 CALL dbt_copy(ri_data%t_3c_int_ctr_1(1, 1), ri_data%t_3c_int_ctr_2(1, 1))
332 CALL dbt_destroy(t_3c_int(1, 1))
333 CALL timestop(handle2)
334
335 CALL timestop(handle)
336
337 END SUBROUTINE
338
339! **************************************************************************************************
340!> \brief ...
341!> \param matrix_1 ...
342!> \param matrix_2 ...
343!> \param name ...
344!> \param unit_nr ...
345! **************************************************************************************************
346 SUBROUTINE check_inverse(matrix_1, matrix_2, name, unit_nr)
347 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_1, matrix_2
348 CHARACTER(len=*), INTENT(IN), OPTIONAL :: name
349 INTEGER, INTENT(IN) :: unit_nr
350
351 CHARACTER(len=default_string_length) :: name_prv
352 REAL(kind=dp) :: error, frob_matrix, frob_matrix_base
353 TYPE(dbcsr_type) :: matrix_tmp
354
355 IF (PRESENT(name)) THEN
356 name_prv = name
357 ELSE
358 CALL dbcsr_get_info(matrix_1, name=name_prv)
359 END IF
360
361 CALL dbcsr_create(matrix_tmp, template=matrix_1)
362 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_1, matrix_2, &
363 0.0_dp, matrix_tmp)
364 frob_matrix_base = dbcsr_frobenius_norm(matrix_tmp)
365 CALL dbcsr_add_on_diag(matrix_tmp, -1.0_dp)
366 frob_matrix = dbcsr_frobenius_norm(matrix_tmp)
367 error = frob_matrix/frob_matrix_base
368 IF (unit_nr > 0) THEN
369 WRITE (unit=unit_nr, fmt="(T3,A,A,A,T73,ES8.1)") &
370 "HFX_RI_INFO| Error for INV(", trim(name_prv), "):", error
371 END IF
372
373 CALL dbcsr_release(matrix_tmp)
374 END SUBROUTINE
375
376! **************************************************************************************************
377!> \brief ...
378!> \param matrix ...
379!> \param matrix_sqrt ...
380!> \param matrix_sqrt_inv ...
381!> \param name ...
382!> \param unit_nr ...
383! **************************************************************************************************
384 SUBROUTINE check_sqrt(matrix, matrix_sqrt, matrix_sqrt_inv, name, unit_nr)
385 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
386 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix_sqrt, matrix_sqrt_inv
387 CHARACTER(len=*), INTENT(IN), OPTIONAL :: name
388 INTEGER, INTENT(IN) :: unit_nr
389
390 CHARACTER(len=default_string_length) :: name_prv
391 REAL(kind=dp) :: frob_matrix
392 TYPE(dbcsr_type) :: matrix_copy, matrix_tmp
393
394 IF (PRESENT(name)) THEN
395 name_prv = name
396 ELSE
397 CALL dbcsr_get_info(matrix, name=name_prv)
398 END IF
399 IF (PRESENT(matrix_sqrt)) THEN
400 CALL dbcsr_create(matrix_tmp, template=matrix)
401 CALL dbcsr_copy(matrix_copy, matrix_sqrt)
402 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sqrt, matrix_copy, &
403 0.0_dp, matrix_tmp)
404 CALL dbcsr_add(matrix_tmp, matrix, 1.0_dp, -1.0_dp)
405 frob_matrix = dbcsr_frobenius_norm(matrix_tmp)
406 IF (unit_nr > 0) THEN
407 WRITE (unit=unit_nr, fmt="(T3,A,A,A,T73,ES8.1)") &
408 "HFX_RI_INFO| Error for SQRT(", trim(name_prv), "):", frob_matrix
409 END IF
410 CALL dbcsr_release(matrix_tmp)
411 CALL dbcsr_release(matrix_copy)
412 END IF
413
414 IF (PRESENT(matrix_sqrt_inv)) THEN
415 CALL dbcsr_create(matrix_tmp, template=matrix)
416 CALL dbcsr_copy(matrix_copy, matrix_sqrt_inv)
417 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_sqrt_inv, matrix_copy, &
418 0.0_dp, matrix_tmp)
419 CALL check_inverse(matrix_tmp, matrix, name="SQRT("//trim(name_prv)//")", unit_nr=unit_nr)
420 CALL dbcsr_release(matrix_tmp)
421 CALL dbcsr_release(matrix_copy)
422 END IF
423
424 END SUBROUTINE
425
426! **************************************************************************************************
427!> \brief Calculate 2-center and 3-center integrals
428!>
429!> 2c: K_1(P, R) = (P|v1|R) and K_2(P, R) = (P|v2|R)
430!> 3c: int_3c(mu, lambda, P) = (mu lambda |v2| P)
431!> v_1 is HF operator, v_2 is RI metric
432!> \param qs_env ...
433!> \param ri_data ...
434!> \param t_2c_int_RI K_2(P, R) note: even with k-point, only need on central cell
435!> \param t_2c_int_pot K_1(P, R)
436!> \param t_3c_int int_3c(mu, lambda, P)
437!> \param do_kpoints ...
438!> \notes the integral tensor arrays are already allocated on entry
439! **************************************************************************************************
440 SUBROUTINE hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_int_RI, t_2c_int_pot, t_3c_int, do_kpoints)
441 TYPE(qs_environment_type), POINTER :: qs_env
442 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
443 TYPE(dbcsr_type), DIMENSION(:), INTENT(OUT) :: t_2c_int_ri, t_2c_int_pot
444 TYPE(dbt_type), DIMENSION(:, :) :: t_3c_int
445 LOGICAL, INTENT(IN), OPTIONAL :: do_kpoints
446
447 CHARACTER(LEN=*), PARAMETER :: routinen = 'hfx_ri_pre_scf_calc_tensors'
448
449 CHARACTER :: symm
450 INTEGER :: handle, i_img, i_mem, ibasis, j_img, &
451 n_mem, natom, nblks, nimg, nkind, &
452 nthreads
453 INTEGER(int_8) :: nze
454 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist_ao_1, dist_ao_2, dist_ri, dist_ri_ext, &
455 ends_array_mc_block_int, ends_array_mc_int, sizes_ao, sizes_ri, sizes_ri_ext, &
456 sizes_ri_ext_split, starts_array_mc_block_int, starts_array_mc_int
457 INTEGER, DIMENSION(3) :: pcoord, pdims
458 INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
459 LOGICAL :: converged, do_kpoints_prv
460 REAL(dp) :: max_ev, min_ev, occ, ri_range
461 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
462 TYPE(dbcsr_distribution_type) :: dbcsr_dist
463 TYPE(dbt_distribution_type) :: t_dist
464 TYPE(dbt_pgrid_type) :: pgrid
465 TYPE(dbt_type) :: t_3c_tmp
466 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_int_batched
467 TYPE(dft_control_type), POINTER :: dft_control
468 TYPE(distribution_2d_type), POINTER :: dist_2d
469 TYPE(distribution_3d_type) :: dist_3d
470 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
471 DIMENSION(:), TARGET :: basis_set_ao, basis_set_ri
472 TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
473 TYPE(mp_cart_type) :: mp_comm_t3c
474 TYPE(mp_para_env_type), POINTER :: para_env
475 TYPE(neighbor_list_3c_type) :: nl_3c
476 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
477 POINTER :: nl_2c_pot, nl_2c_ri
478 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
479 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
480 TYPE(qs_ks_env_type), POINTER :: ks_env
481
482 CALL timeset(routinen, handle)
483 NULLIFY (col_bsize, row_bsize, dist_2d, nl_2c_pot, nl_2c_ri, &
484 particle_set, qs_kind_set, ks_env, para_env)
485
486 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, particle_set=particle_set, &
487 distribution_2d=dist_2d, ks_env=ks_env, dft_control=dft_control, para_env=para_env)
488
489 ri_range = 0.0_dp
490 do_kpoints_prv = .false.
491 IF (PRESENT(do_kpoints)) do_kpoints_prv = do_kpoints
492 nimg = 1
493 IF (do_kpoints_prv) THEN
494 nimg = ri_data%nimg
495 ri_range = ri_data%kp_RI_range
496 END IF
497
498 ALLOCATE (sizes_ri(natom), sizes_ao(natom))
499 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
500 CALL basis_set_list_setup(basis_set_ri, ri_data%ri_basis_type, qs_kind_set)
501 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ri, basis=basis_set_ri)
502 CALL basis_set_list_setup(basis_set_ao, ri_data%orb_basis_type, qs_kind_set)
503 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ao, basis=basis_set_ao)
504
505 DO ibasis = 1, SIZE(basis_set_ao)
506 orb_basis => basis_set_ao(ibasis)%gto_basis_set
507 ri_basis => basis_set_ri(ibasis)%gto_basis_set
508 ! interaction radii should be based on eps_pgf_orb controlled in RI section
509 ! (since hartree-fock needs very tight eps_pgf_orb for Kohn-Sham/Fock matrix but eps_pgf_orb
510 ! can be much looser in RI HFX since no systematic error is introduced with tensor sparsity)
511 CALL init_interaction_radii_orb_basis(orb_basis, ri_data%eps_pgf_orb)
512 CALL init_interaction_radii_orb_basis(ri_basis, ri_data%eps_pgf_orb)
513 END DO
514
515 n_mem = ri_data%n_mem
516 CALL create_tensor_batches(sizes_ri, n_mem, starts_array_mc_int, ends_array_mc_int, &
517 starts_array_mc_block_int, ends_array_mc_block_int)
518
519 DEALLOCATE (starts_array_mc_int, ends_array_mc_int)
520
521 !We separate the K-points and standard 3c integrals, because handle quite differently
522 IF (.NOT. do_kpoints_prv) THEN
523
524 nthreads = 1
525!$ nthreads = omp_get_num_threads()
526 pdims = 0
527 CALL dbt_pgrid_create(para_env, pdims, pgrid, tensor_dims=[max(1, natom/(n_mem*nthreads)), natom, natom])
528
529 ALLOCATE (t_3c_int_batched(1, 1))
530 CALL create_3c_tensor(t_3c_int_batched(1, 1), dist_ri, dist_ao_1, dist_ao_2, pgrid, &
531 sizes_ri, sizes_ao, sizes_ao, map1=[1], map2=[2, 3], &
532 name="(RI | AO AO)")
533
534 CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set, atomic_kind_set=atomic_kind_set)
535 CALL dbt_mp_environ_pgrid(pgrid, pdims, pcoord)
536 CALL mp_comm_t3c%create(pgrid%mp_comm_2d, 3, pdims)
537 CALL distribution_3d_create(dist_3d, dist_ri, dist_ao_1, dist_ao_2, &
538 nkind, particle_set, mp_comm_t3c, own_comm=.true.)
539 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
540
541 CALL create_3c_tensor(t_3c_int(1, 1), dist_ri, dist_ao_1, dist_ao_2, ri_data%pgrid, &
542 ri_data%bsizes_RI_split, ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
543 map1=[1], map2=[2, 3], &
544 name="O (RI AO | AO)")
545
546 CALL build_3c_neighbor_lists(nl_3c, basis_set_ri, basis_set_ao, basis_set_ao, dist_3d, ri_data%ri_metric, &
547 "HFX_3c_nl", qs_env, op_pos=1, sym_jk=.true., own_dist=.true.)
548
549 DO i_mem = 1, n_mem
550 CALL build_3c_integrals(t_3c_int_batched, ri_data%filter_eps/2, qs_env, nl_3c, &
551 basis_set_ri, basis_set_ao, basis_set_ao, &
552 ri_data%ri_metric, int_eps=ri_data%eps_schwarz, op_pos=1, &
553 desymmetrize=.false., &
554 bounds_i=[starts_array_mc_block_int(i_mem), ends_array_mc_block_int(i_mem)])
555 CALL dbt_copy(t_3c_int_batched(1, 1), t_3c_int(1, 1), summation=.true., move_data=.true.)
556 CALL dbt_filter(t_3c_int(1, 1), ri_data%filter_eps/2)
557 END DO
558
559 CALL dbt_destroy(t_3c_int_batched(1, 1))
560
561 CALL neighbor_list_3c_destroy(nl_3c)
562
563 CALL dbt_create(t_3c_int(1, 1), t_3c_tmp)
564
565 IF (ri_data%flavor == ri_pmat) THEN ! desymmetrize
566 ! desymmetrize
567 CALL dbt_copy(t_3c_int(1, 1), t_3c_tmp)
568 CALL dbt_copy(t_3c_tmp, t_3c_int(1, 1), order=[1, 3, 2], summation=.true., move_data=.true.)
569
570 ! For RI-RHO filter_eps_storage is reserved for screening tensor contracted with RI-metric
571 ! with RI metric but not to bare integral tensor
572 CALL dbt_filter(t_3c_int(1, 1), ri_data%filter_eps)
573 ELSE
574 CALL dbt_filter(t_3c_int(1, 1), ri_data%filter_eps_storage/2)
575 END IF
576
577 CALL dbt_destroy(t_3c_tmp)
578
579 ELSE !do_kpoints
580
581 nthreads = 1
582!$ nthreads = omp_get_num_threads()
583 pdims = 0
584 CALL dbt_pgrid_create(para_env, pdims, pgrid, tensor_dims=[natom, natom, max(1, natom/(n_mem*nthreads))])
585
586 !In k-points HFX, we stack all images along the RI direction in the same tensors, in order
587 !to avoid storing nimg x ncell_RI different tensors (very memory intensive)
588 nblks = SIZE(ri_data%bsizes_RI_split)
589 ALLOCATE (sizes_ri_ext(natom*ri_data%ncell_RI), sizes_ri_ext_split(nblks*ri_data%ncell_RI))
590 DO i_img = 1, ri_data%ncell_RI
591 sizes_ri_ext((i_img - 1)*natom + 1:i_img*natom) = sizes_ri(:)
592 sizes_ri_ext_split((i_img - 1)*nblks + 1:i_img*nblks) = ri_data%bsizes_RI_split(:)
593 END DO
594
595 CALL create_3c_tensor(t_3c_tmp, dist_ao_1, dist_ao_2, dist_ri, &
596 pgrid, sizes_ao, sizes_ao, sizes_ri, map1=[1, 2], map2=[3], &
597 name="(AO AO | RI)")
598 CALL dbt_destroy(t_3c_tmp)
599
600 !For the integrals to work, the distribution along the RI direction must be repeated
601 ALLOCATE (dist_ri_ext(natom*ri_data%ncell_RI))
602 DO i_img = 1, ri_data%ncell_RI
603 dist_ri_ext((i_img - 1)*natom + 1:i_img*natom) = dist_ri(:)
604 END DO
605
606 ALLOCATE (t_3c_int_batched(nimg, 1))
607 CALL dbt_distribution_new(t_dist, pgrid, dist_ao_1, dist_ao_2, dist_ri_ext)
608 CALL dbt_create(t_3c_int_batched(1, 1), "KP_3c_ints", t_dist, [1, 2], [3], &
609 sizes_ao, sizes_ao, sizes_ri_ext)
610 DO i_img = 2, nimg
611 CALL dbt_create(t_3c_int_batched(1, 1), t_3c_int_batched(i_img, 1))
612 END DO
613 CALL dbt_distribution_destroy(t_dist)
614
615 CALL get_qs_env(qs_env, nkind=nkind, particle_set=particle_set, atomic_kind_set=atomic_kind_set)
616 CALL dbt_mp_environ_pgrid(pgrid, pdims, pcoord)
617 CALL mp_comm_t3c%create(pgrid%mp_comm_2d, 3, pdims)
618 CALL distribution_3d_create(dist_3d, dist_ao_1, dist_ao_2, dist_ri, &
619 nkind, particle_set, mp_comm_t3c, own_comm=.true.)
620 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
621
622 ! create 3c tensor for storage of ints
623 CALL build_3c_neighbor_lists(nl_3c, basis_set_ao, basis_set_ao, basis_set_ri, dist_3d, ri_data%ri_metric, &
624 "HFX_3c_nl", qs_env, op_pos=2, sym_ij=.false., own_dist=.true.)
625
626 CALL create_3c_tensor(t_3c_int(1, 1), dist_ri, dist_ao_1, dist_ao_2, ri_data%pgrid, &
627 sizes_ri_ext_split, ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
628 map1=[1], map2=[2, 3], &
629 name="O (RI AO | AO)")
630 DO j_img = 2, nimg
631 CALL dbt_create(t_3c_int(1, 1), t_3c_int(1, j_img))
632 END DO
633
634 CALL dbt_create(t_3c_int(1, 1), t_3c_tmp)
635 DO i_mem = 1, n_mem
636 CALL build_3c_integrals(t_3c_int_batched, ri_data%filter_eps, qs_env, nl_3c, &
637 basis_set_ao, basis_set_ao, basis_set_ri, &
638 ri_data%ri_metric, int_eps=ri_data%eps_schwarz, op_pos=2, &
639 desymmetrize=.false., do_kpoints=.true., do_hfx_kpoints=.true., &
640 bounds_k=[starts_array_mc_block_int(i_mem), ends_array_mc_block_int(i_mem)], &
641 ri_range=ri_range, img_to_ri_cell=ri_data%img_to_RI_cell)
642
643 DO i_img = 1, nimg
644 !we start with (mu^0 sigma^i | P^j) and finish with (P^i | mu^0 sigma^j)
645 CALL get_tensor_occupancy(t_3c_int_batched(i_img, 1), nze, occ)
646 IF (nze > 0) THEN
647 CALL dbt_copy(t_3c_int_batched(i_img, 1), t_3c_tmp, order=[3, 2, 1], move_data=.true.)
648 CALL dbt_filter(t_3c_tmp, ri_data%filter_eps)
649 CALL dbt_copy(t_3c_tmp, t_3c_int(1, i_img), order=[1, 3, 2], &
650 summation=.true., move_data=.true.)
651 END IF
652 END DO
653 END DO
654
655 DO i_img = 1, nimg
656 CALL dbt_destroy(t_3c_int_batched(i_img, 1))
657 END DO
658 DEALLOCATE (t_3c_int_batched)
659 CALL neighbor_list_3c_destroy(nl_3c)
660 CALL dbt_destroy(t_3c_tmp)
661 END IF
662 CALL dbt_pgrid_destroy(pgrid)
663
664 CALL build_2c_neighbor_lists(nl_2c_pot, basis_set_ri, basis_set_ri, ri_data%hfx_pot, &
665 "HFX_2c_nl_pot", qs_env, sym_ij=.NOT. do_kpoints_prv, &
666 dist_2d=dist_2d)
667
668 CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
669 ALLOCATE (row_bsize(SIZE(sizes_ri)))
670 ALLOCATE (col_bsize(SIZE(sizes_ri)))
671 row_bsize(:) = sizes_ri
672 col_bsize(:) = sizes_ri
673
674 !Use non-symmetric nl and matrices for k-points
675 symm = dbcsr_type_symmetric
676 IF (do_kpoints_prv) symm = dbcsr_type_no_symmetry
677
678 CALL dbcsr_create(t_2c_int_pot(1), "(R|P) HFX", dbcsr_dist, symm, row_bsize, col_bsize)
679 DO i_img = 2, nimg
680 CALL dbcsr_create(t_2c_int_pot(i_img), template=t_2c_int_pot(1))
681 END DO
682
683 IF (.NOT. ri_data%same_op) THEN
684 CALL dbcsr_create(t_2c_int_ri(1), "(R|P) HFX", dbcsr_dist, symm, row_bsize, col_bsize)
685 DO i_img = 2, nimg
686 CALL dbcsr_create(t_2c_int_ri(i_img), template=t_2c_int_ri(1))
687 END DO
688 END IF
689 DEALLOCATE (col_bsize, row_bsize)
690
691 CALL dbcsr_distribution_release(dbcsr_dist)
692
693 CALL build_2c_integrals(t_2c_int_pot, ri_data%filter_eps_2c, qs_env, nl_2c_pot, basis_set_ri, basis_set_ri, &
694 ri_data%hfx_pot, do_kpoints=do_kpoints_prv, do_hfx_kpoints=do_kpoints_prv)
695 CALL release_neighbor_list_sets(nl_2c_pot)
696
697 IF (.NOT. ri_data%same_op) THEN
698 CALL build_2c_neighbor_lists(nl_2c_ri, basis_set_ri, basis_set_ri, ri_data%ri_metric, &
699 "HFX_2c_nl_RI", qs_env, sym_ij=.NOT. do_kpoints_prv, &
700 dist_2d=dist_2d)
701
702 CALL build_2c_integrals(t_2c_int_ri, ri_data%filter_eps_2c, qs_env, nl_2c_ri, basis_set_ri, basis_set_ri, &
703 ri_data%ri_metric, do_kpoints=do_kpoints_prv, do_hfx_kpoints=do_kpoints_prv)
704
705 CALL release_neighbor_list_sets(nl_2c_ri)
706 END IF
707
708 DO ibasis = 1, SIZE(basis_set_ao)
709 orb_basis => basis_set_ao(ibasis)%gto_basis_set
710 ri_basis => basis_set_ri(ibasis)%gto_basis_set
711 ! reset interaction radii of orb basis
712 CALL init_interaction_radii_orb_basis(orb_basis, dft_control%qs_control%eps_pgf_orb)
713 CALL init_interaction_radii_orb_basis(ri_basis, dft_control%qs_control%eps_pgf_orb)
714 END DO
715
716 IF (ri_data%calc_condnum) THEN
717 CALL arnoldi_extremal(t_2c_int_pot(1), max_ev, min_ev, threshold=ri_data%eps_lanczos, &
718 max_iter=ri_data%max_iter_lanczos, converged=converged)
719
720 IF (.NOT. converged) THEN
721 cpwarn("Condition number estimate of (P|Q) (HFX potential) is not reliable (not converged).")
722 END IF
723
724 IF (ri_data%unit_nr > 0) THEN
725 WRITE (ri_data%unit_nr, '(T2,A)') "2-Norm Condition Number of (P|Q) integrals (HFX potential)"
726 IF (min_ev > 0) THEN
727 WRITE (ri_data%unit_nr, '(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
728 "CN : max/min ev: ", max_ev, " / ", min_ev, "=", max_ev/min_ev, "Log(2-CN):", log10(max_ev/min_ev)
729 ELSE
730 WRITE (ri_data%unit_nr, '(T4,A,ES11.3E3,T32,A,ES11.3E3,T63,A)') &
731 "CN : max/min ev: ", max_ev, " / ", min_ev, "Log(CN): infinity"
732 END IF
733 END IF
734
735 IF (.NOT. ri_data%same_op) THEN
736 CALL arnoldi_extremal(t_2c_int_ri(1), max_ev, min_ev, threshold=ri_data%eps_lanczos, &
737 max_iter=ri_data%max_iter_lanczos, converged=converged)
738
739 IF (.NOT. converged) THEN
740 cpwarn("Condition number estimate of (P|Q) matrix (RI metric) is not reliable (not converged).")
741 END IF
742
743 IF (ri_data%unit_nr > 0) THEN
744 WRITE (ri_data%unit_nr, '(T2,A)') "2-Norm Condition Number of (P|Q) integrals (RI metric)"
745 IF (min_ev > 0) THEN
746 WRITE (ri_data%unit_nr, '(T4,A,ES11.3E3,T32,A,ES11.3E3,A4,ES11.3E3,T63,A,F8.4)') &
747 "CN : max/min ev: ", max_ev, " / ", min_ev, "=", max_ev/min_ev, "Log(2-CN):", log10(max_ev/min_ev)
748 ELSE
749 WRITE (ri_data%unit_nr, '(T4,A,ES11.3E3,T32,A,ES11.3E3,T63,A)') &
750 "CN : max/min ev: ", max_ev, " / ", min_ev, "Log(CN): infinity"
751 END IF
752 END IF
753 END IF
754 END IF
755
756 CALL timestop(handle)
757 END SUBROUTINE
758
759! **************************************************************************************************
760!> \brief Pre-SCF steps in rho flavor of RI HFX
761!>
762!> K(P, S) = sum_{R,Q} K_2(P, R)^{-1} K_1(R, Q) K_2(Q, S)^{-1}
763!> Calculate B(mu, lambda, R) = sum_P int_3c(mu, lambda, P) K(P, R)
764!> \param qs_env ...
765!> \param ri_data ...
766! **************************************************************************************************
767 SUBROUTINE hfx_ri_pre_scf_pmat(qs_env, ri_data)
768 TYPE(qs_environment_type), POINTER :: qs_env
769 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
770
771 CHARACTER(LEN=*), PARAMETER :: routinen = 'hfx_ri_pre_scf_Pmat'
772
773 INTEGER :: handle, handle2, i_mem, j_mem, &
774 n_dependent, unit_nr, unit_nr_dbcsr
775 INTEGER(int_8) :: nflop, nze, nze_o
776 INTEGER, ALLOCATABLE, DIMENSION(:) :: batch_ranges_ao, batch_ranges_ri
777 INTEGER, DIMENSION(2, 1) :: bounds_i
778 INTEGER, DIMENSION(2, 2) :: bounds_j
779 INTEGER, DIMENSION(3) :: dims_3c
780 REAL(kind=dp) :: compression_factor, memory_3c, occ, &
781 threshold
782 TYPE(cp_blacs_env_type), POINTER :: blacs_env
783 TYPE(dbcsr_type), DIMENSION(1) :: t_2c_int_mat, t_2c_op_pot, t_2c_op_ri, &
784 t_2c_tmp, t_2c_tmp_2
785 TYPE(dbt_type) :: t_3c_2
786 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_2c_int, t_2c_work
787 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_int_1
788 TYPE(mp_para_env_type), POINTER :: para_env
789
790 CALL timeset(routinen, handle)
791
792 unit_nr_dbcsr = ri_data%unit_nr_dbcsr
793 unit_nr = ri_data%unit_nr
794
795 CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
796
797 CALL timeset(routinen//"_int", handle2)
798
799 ALLOCATE (t_2c_int(1), t_2c_work(1), t_3c_int_1(1, 1))
800 CALL hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_op_ri, t_2c_op_pot, t_3c_int_1)
801
802 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.)
803
804 CALL dbt_destroy(t_3c_int_1(1, 1))
805
806 CALL timestop(handle2)
807
808 CALL timeset(routinen//"_2c", handle2)
809
810 IF (ri_data%same_op) t_2c_op_ri(1) = t_2c_op_pot(1)
811 CALL dbcsr_create(t_2c_int_mat(1), template=t_2c_op_ri(1), matrix_type=dbcsr_type_no_symmetry)
812 threshold = max(ri_data%filter_eps, 1.0e-12_dp)
813
814 SELECT CASE (ri_data%t2c_method)
815 CASE (hfx_ri_do_2c_iter)
816 CALL invert_hotelling(t_2c_int_mat(1), t_2c_op_ri(1), &
817 threshold=threshold, silent=.false.)
818 CASE (hfx_ri_do_2c_cholesky)
819 CALL dbcsr_copy(t_2c_int_mat(1), t_2c_op_ri(1))
820 CALL cp_dbcsr_cholesky_decompose(t_2c_int_mat(1), para_env=para_env, blacs_env=blacs_env)
821 CALL cp_dbcsr_cholesky_invert(t_2c_int_mat(1), para_env=para_env, blacs_env=blacs_env, upper_to_full=.true.)
822 CASE (hfx_ri_do_2c_diag)
823 CALL dbcsr_copy(t_2c_int_mat(1), t_2c_op_ri(1))
824 CALL cp_dbcsr_power(t_2c_int_mat(1), -1.0_dp, ri_data%eps_eigval, n_dependent, &
825 para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
826 END SELECT
827
828 IF (ri_data%check_2c_inv) THEN
829 CALL check_inverse(t_2c_int_mat(1), t_2c_op_ri(1), unit_nr=unit_nr)
830 END IF
831
832 !Need to save the (P|Q)^-1 tensor for forces (inverse metric if not same_op)
833 CALL dbt_create(t_2c_int_mat(1), t_2c_work(1))
834 CALL dbt_copy_matrix_to_tensor(t_2c_int_mat(1), t_2c_work(1))
835 CALL dbt_copy(t_2c_work(1), ri_data%t_2c_inv(1, 1), move_data=.true.)
836 CALL dbt_destroy(t_2c_work(1))
837 CALL dbt_filter(ri_data%t_2c_inv(1, 1), ri_data%filter_eps)
838 IF (.NOT. ri_data%same_op) THEN
839 !Also save the RI (P|Q) integral
840 CALL dbt_create(t_2c_op_pot(1), t_2c_work(1))
841 CALL dbt_copy_matrix_to_tensor(t_2c_op_pot(1), t_2c_work(1))
842 CALL dbt_copy(t_2c_work(1), ri_data%t_2c_pot(1, 1), move_data=.true.)
843 CALL dbt_destroy(t_2c_work(1))
844 CALL dbt_filter(ri_data%t_2c_pot(1, 1), ri_data%filter_eps)
845 END IF
846
847 IF (ri_data%same_op) THEN
848 CALL dbcsr_release(t_2c_op_pot(1))
849 ELSE
850 CALL dbcsr_create(t_2c_tmp(1), template=t_2c_op_ri(1), matrix_type=dbcsr_type_no_symmetry)
851 CALL dbcsr_create(t_2c_tmp_2(1), template=t_2c_op_ri(1), matrix_type=dbcsr_type_no_symmetry)
852 CALL dbcsr_release(t_2c_op_ri(1))
853 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), &
854 filter_eps=ri_data%filter_eps)
855
856 CALL dbcsr_release(t_2c_op_pot(1))
857 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), &
858 filter_eps=ri_data%filter_eps)
859 CALL dbcsr_release(t_2c_tmp(1))
860 CALL dbcsr_release(t_2c_int_mat(1))
861 t_2c_int_mat(1) = t_2c_tmp_2(1)
862 END IF
863
864 CALL dbt_create(t_2c_int_mat(1), t_2c_int(1), name="(RI|RI)")
865 CALL dbt_copy_matrix_to_tensor(t_2c_int_mat(1), t_2c_int(1))
866 CALL dbcsr_release(t_2c_int_mat(1))
867 CALL dbt_copy(t_2c_int(1), ri_data%t_2c_int(1, 1), move_data=.true.)
868 CALL dbt_destroy(t_2c_int(1))
869 CALL dbt_filter(ri_data%t_2c_int(1, 1), ri_data%filter_eps)
870
871 CALL timestop(handle2)
872
873 CALL dbt_create(ri_data%t_3c_int_ctr_3(1, 1), t_3c_2)
874
875 CALL dbt_get_info(ri_data%t_3c_int_ctr_3(1, 1), nfull_total=dims_3c)
876
877 memory_3c = 0.0_dp
878 nze_o = 0
879
880 ALLOCATE (batch_ranges_ri(ri_data%n_mem_RI + 1))
881 ALLOCATE (batch_ranges_ao(ri_data%n_mem + 1))
882 batch_ranges_ri(:ri_data%n_mem_RI) = ri_data%starts_array_RI_mem_block(:)
883 batch_ranges_ri(ri_data%n_mem_RI + 1) = ri_data%ends_array_RI_mem_block(ri_data%n_mem_RI) + 1
884 batch_ranges_ao(:ri_data%n_mem) = ri_data%starts_array_mem_block(:)
885 batch_ranges_ao(ri_data%n_mem + 1) = ri_data%ends_array_mem_block(ri_data%n_mem) + 1
886
887 CALL dbt_batched_contract_init(ri_data%t_3c_int_ctr_3(1, 1), batch_range_1=batch_ranges_ri, &
888 batch_range_2=batch_ranges_ao)
889 CALL dbt_batched_contract_init(t_3c_2, batch_range_1=batch_ranges_ri, &
890 batch_range_2=batch_ranges_ao)
891
892 DO i_mem = 1, ri_data%n_mem_RI
893 bounds_i(:, 1) = [ri_data%starts_array_RI_mem(i_mem), ri_data%ends_array_RI_mem(i_mem)]
894
895 CALL dbt_batched_contract_init(ri_data%t_2c_int(1, 1))
896 DO j_mem = 1, ri_data%n_mem
897 bounds_j(:, 1) = [ri_data%starts_array_mem(j_mem), ri_data%ends_array_mem(j_mem)]
898 bounds_j(:, 2) = [1, dims_3c(3)]
899 CALL timeset(routinen//"_RIx3C", handle2)
900 CALL dbt_contract(1.0_dp, ri_data%t_2c_int(1, 1), ri_data%t_3c_int_ctr_3(1, 1), &
901 0.0_dp, t_3c_2, &
902 contract_1=[2], notcontract_1=[1], &
903 contract_2=[1], notcontract_2=[2, 3], &
904 map_1=[1], map_2=[2, 3], filter_eps=ri_data%filter_eps_storage, &
905 bounds_2=bounds_i, &
906 bounds_3=bounds_j, &
907 unit_nr=unit_nr_dbcsr, &
908 flop=nflop)
909
910 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
911 CALL timestop(handle2)
912
913 CALL timeset(routinen//"_copy_2", handle2)
914 CALL dbt_copy(t_3c_2, ri_data%t_3c_int_ctr_1(1, 1), order=[2, 1, 3], move_data=.true.)
915
916 CALL get_tensor_occupancy(ri_data%t_3c_int_ctr_1(1, 1), nze, occ)
917 nze_o = nze_o + nze
918
919 CALL compress_tensor(ri_data%t_3c_int_ctr_1(1, 1), ri_data%blk_indices(j_mem, i_mem)%ind, &
920 ri_data%store_3c(j_mem, i_mem), ri_data%filter_eps_storage, memory_3c)
921
922 CALL timestop(handle2)
923 END DO
924 CALL dbt_batched_contract_finalize(ri_data%t_2c_int(1, 1))
925 END DO
926 CALL dbt_batched_contract_finalize(t_3c_2)
927 CALL dbt_batched_contract_finalize(ri_data%t_3c_int_ctr_3(1, 1))
928
929 CALL para_env%sum(memory_3c)
930 compression_factor = real(nze_o, dp)*1.0e-06*8.0_dp/memory_3c
931
932 IF (unit_nr > 0) THEN
933 WRITE (unit=unit_nr, fmt="((T3,A,T66,F11.2,A4))") &
934 "MEMORY_INFO| Memory for 3-center HFX integrals (compressed):", memory_3c, ' MiB'
935
936 WRITE (unit=unit_nr, fmt="((T3,A,T60,F21.2))") &
937 "MEMORY_INFO| Compression factor: ", compression_factor
938 END IF
939
940 CALL dbt_clear(ri_data%t_2c_int(1, 1))
941 CALL dbt_destroy(t_3c_2)
942
943 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.)
944
945 CALL timestop(handle)
946 END SUBROUTINE
947
948! **************************************************************************************************
949!> \brief Sorts 2d indices w.r.t. rows and columns
950!> \param blk_ind ...
951! **************************************************************************************************
952 SUBROUTINE sort_unique_blkind_2d(blk_ind)
953 INTEGER, ALLOCATABLE, DIMENSION(:, :), &
954 INTENT(INOUT) :: blk_ind
955
956 INTEGER :: end_ind, iblk, iblk_all, irow, nblk, &
957 ncols, start_ind
958 INTEGER, ALLOCATABLE, DIMENSION(:) :: ind_1, ind_2, sort_1, sort_2
959 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: blk_ind_tmp
960
961 nblk = SIZE(blk_ind, 1)
962
963 ALLOCATE (sort_1(nblk))
964 ALLOCATE (ind_1(nblk))
965
966 sort_1(:) = blk_ind(:, 1)
967 CALL sort(sort_1, nblk, ind_1)
968
969 blk_ind(:, :) = blk_ind(ind_1, :)
970
971 start_ind = 1
972
973 DO WHILE (start_ind <= nblk)
974 irow = blk_ind(start_ind, 1)
975 end_ind = start_ind
976
977 IF (end_ind + 1 <= nblk) THEN
978 DO WHILE (blk_ind(end_ind + 1, 1) == irow)
979 end_ind = end_ind + 1
980 IF (end_ind + 1 > nblk) EXIT
981 END DO
982 END IF
983
984 ncols = end_ind - start_ind + 1
985 ALLOCATE (sort_2(ncols))
986 ALLOCATE (ind_2(ncols))
987 sort_2(:) = blk_ind(start_ind:end_ind, 2)
988 CALL sort(sort_2, ncols, ind_2)
989 ind_2 = ind_2 + start_ind - 1
990
991 blk_ind(start_ind:end_ind, :) = blk_ind(ind_2, :)
992 start_ind = end_ind + 1
993
994 DEALLOCATE (sort_2, ind_2)
995 END DO
996
997 ALLOCATE (blk_ind_tmp(nblk, 2))
998 blk_ind_tmp = 0
999
1000 iblk = 0
1001 DO iblk_all = 1, nblk
1002 IF (iblk >= 1) THEN
1003 IF (all(blk_ind_tmp(iblk, :) == blk_ind(iblk_all, :))) THEN
1004 cycle
1005 END IF
1006 END IF
1007 iblk = iblk + 1
1008 blk_ind_tmp(iblk, :) = blk_ind(iblk_all, :)
1009 END DO
1010 nblk = iblk
1011
1012 DEALLOCATE (blk_ind)
1013 ALLOCATE (blk_ind(nblk, 2))
1014
1015 blk_ind(:, :) = blk_ind_tmp(:nblk, :)
1016
1017 END SUBROUTINE
1018
1019! **************************************************************************************************
1020!> \brief ...
1021!> \param qs_env ...
1022!> \param ri_data ...
1023!> \param ks_matrix ...
1024!> \param ehfx ...
1025!> \param mos ...
1026!> \param rho_ao ...
1027!> \param geometry_did_change ...
1028!> \param nspins ...
1029!> \param hf_fraction ...
1030! **************************************************************************************************
1031 SUBROUTINE hfx_ri_update_ks(qs_env, ri_data, ks_matrix, ehfx, mos, rho_ao, &
1032 geometry_did_change, nspins, hf_fraction)
1033 TYPE(qs_environment_type), POINTER :: qs_env
1034 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
1035 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
1036 REAL(kind=dp), INTENT(OUT) :: ehfx
1037 TYPE(mo_set_type), DIMENSION(:), INTENT(IN), &
1038 OPTIONAL :: mos
1039 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
1040 LOGICAL, INTENT(IN) :: geometry_did_change
1041 INTEGER, INTENT(IN) :: nspins
1042 REAL(kind=dp), INTENT(IN) :: hf_fraction
1043
1044 CHARACTER(LEN=*), PARAMETER :: routinen = 'hfx_ri_update_ks'
1045
1046 CHARACTER(1) :: mtype
1047 INTEGER :: handle, handle2, i, ispin, j
1048 INTEGER(int_8) :: nblks
1049 INTEGER, DIMENSION(2) :: homo
1050 LOGICAL :: is_antisymmetric
1051 REAL(dp) :: etmp, fac
1052 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
1053 TYPE(cp_fm_type), POINTER :: mo_coeff
1054 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: my_ks_matrix, my_rho_ao
1055 TYPE(dbcsr_type), DIMENSION(2) :: mo_coeff_b
1056 TYPE(dbcsr_type), POINTER :: mo_coeff_b_tmp
1057 TYPE(mp_para_env_type), POINTER :: para_env
1058
1059 CALL timeset(routinen, handle)
1060
1061 IF (nspins == 1) THEN
1062 fac = 0.5_dp*hf_fraction
1063 ELSE
1064 fac = 1.0_dp*hf_fraction
1065 END IF
1066
1067 !If incoming assymetric matrices, need to convert to normal
1068 NULLIFY (my_ks_matrix, my_rho_ao)
1069 CALL dbcsr_get_info(ks_matrix(1, 1)%matrix, matrix_type=mtype)
1070 is_antisymmetric = mtype == dbcsr_type_antisymmetric
1071 IF (is_antisymmetric) THEN
1072 ALLOCATE (my_ks_matrix(SIZE(ks_matrix, 1), SIZE(ks_matrix, 2)))
1073 ALLOCATE (my_rho_ao(SIZE(rho_ao, 1), SIZE(rho_ao, 2)))
1074
1075 DO i = 1, SIZE(ks_matrix, 1)
1076 DO j = 1, SIZE(ks_matrix, 2)
1077 ALLOCATE (my_ks_matrix(i, j)%matrix, my_rho_ao(i, j)%matrix)
1078 CALL dbcsr_create(my_ks_matrix(i, j)%matrix, template=ks_matrix(i, j)%matrix, &
1079 matrix_type=dbcsr_type_no_symmetry)
1080 CALL dbcsr_desymmetrize(ks_matrix(i, j)%matrix, my_ks_matrix(i, j)%matrix)
1081 CALL dbcsr_create(my_rho_ao(i, j)%matrix, template=rho_ao(i, j)%matrix, &
1082 matrix_type=dbcsr_type_no_symmetry)
1083 CALL dbcsr_desymmetrize(rho_ao(i, j)%matrix, my_rho_ao(i, j)%matrix)
1084 END DO
1085 END DO
1086 ELSE
1087 my_ks_matrix => ks_matrix
1088 my_rho_ao => rho_ao
1089 END IF
1090
1091 !Case analysis on RI_FLAVOR: we switch if the input flavor is MO, there is no provided MO, and
1092 ! the current flavor is not yet RHO. We switch back to MO if there are
1093 ! MOs available and the current flavor is actually RHO
1094 IF (ri_data%input_flavor == ri_mo .AND. (.NOT. PRESENT(mos)) .AND. ri_data%flavor == ri_mo) THEN
1095 CALL switch_ri_flavor(ri_data, qs_env)
1096 ELSE IF (ri_data%input_flavor == ri_mo .AND. PRESENT(mos) .AND. ri_data%flavor == ri_pmat) THEN
1097 CALL switch_ri_flavor(ri_data, qs_env)
1098 END IF
1099
1100 SELECT CASE (ri_data%flavor)
1101 CASE (ri_mo)
1102 cpassert(PRESENT(mos))
1103 CALL timeset(routinen//"_MO", handle2)
1104
1105 DO ispin = 1, nspins
1106 NULLIFY (mo_coeff_b_tmp)
1107 cpassert(mos(ispin)%uniform_occupation)
1108 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, eigenvalues=mo_eigenvalues, mo_coeff_b=mo_coeff_b_tmp)
1109
1110 IF (.NOT. mos(ispin)%use_mo_coeff_b) CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b_tmp)
1111 CALL dbcsr_copy(mo_coeff_b(ispin), mo_coeff_b_tmp)
1112 END DO
1113
1114 DO ispin = 1, nspins
1115 CALL dbcsr_scale(mo_coeff_b(ispin), sqrt(mos(ispin)%maxocc))
1116 homo(ispin) = mos(ispin)%homo
1117 END DO
1118
1119 CALL timestop(handle2)
1120
1121 CALL hfx_ri_update_ks_mo(qs_env, ri_data, my_ks_matrix, mo_coeff_b, homo, &
1122 geometry_did_change, nspins, fac)
1123 CASE (ri_pmat)
1124
1125 NULLIFY (para_env)
1126 CALL get_qs_env(qs_env, para_env=para_env)
1127 DO ispin = 1, SIZE(my_rho_ao, 1)
1128 nblks = dbcsr_get_num_blocks(my_rho_ao(ispin, 1)%matrix)
1129 CALL para_env%sum(nblks)
1130 IF (nblks == 0) THEN
1131 cpabort("received empty density matrix")
1132 END IF
1133 END DO
1134
1135 CALL hfx_ri_update_ks_pmat(qs_env, ri_data, my_ks_matrix, my_rho_ao, &
1136 geometry_did_change, nspins, fac)
1137
1138 END SELECT
1139
1140 DO ispin = 1, nspins
1141 CALL dbcsr_release(mo_coeff_b(ispin))
1142 END DO
1143
1144 DO ispin = 1, nspins
1145 CALL dbcsr_filter(my_ks_matrix(ispin, 1)%matrix, ri_data%filter_eps)
1146 END DO
1147
1148 CALL timeset(routinen//"_energy", handle2)
1149 ! Calculate the exchange energy
1150 ehfx = 0.0_dp
1151 DO ispin = 1, nspins
1152 CALL dbcsr_dot(my_ks_matrix(ispin, 1)%matrix, my_rho_ao(ispin, 1)%matrix, &
1153 etmp)
1154 ehfx = ehfx + 0.5_dp*etmp
1155
1156 END DO
1157 CALL timestop(handle2)
1158
1159 !Anti-symmetric case
1160 IF (is_antisymmetric) THEN
1161 DO i = 1, SIZE(ks_matrix, 1)
1162 DO j = 1, SIZE(ks_matrix, 2)
1163 CALL dbcsr_complete_redistribute(my_ks_matrix(i, j)%matrix, ks_matrix(i, j)%matrix)
1164 CALL dbcsr_complete_redistribute(my_rho_ao(i, j)%matrix, rho_ao(i, j)%matrix)
1165 END DO
1166 END DO
1167 CALL dbcsr_deallocate_matrix_set(my_ks_matrix)
1168 CALL dbcsr_deallocate_matrix_set(my_rho_ao)
1169 END IF
1170
1171 CALL timestop(handle)
1172 END SUBROUTINE
1173
1174! **************************************************************************************************
1175!> \brief Calculate Fock (AKA Kohn-Sham) matrix in MO flavor
1176!>
1177!> C(mu, i) (MO coefficients)
1178!> M(mu, i, R) = sum_nu B(mu, nu, R) C(nu, i)
1179!> KS(mu, lambda) = sum_{i,R} M(mu, i, R) M(lambda, i, R)
1180!> \param qs_env ...
1181!> \param ri_data ...
1182!> \param ks_matrix ...
1183!> \param mo_coeff C(mu, i)
1184!> \param homo ...
1185!> \param geometry_did_change ...
1186!> \param nspins ...
1187!> \param fac ...
1188! **************************************************************************************************
1189 SUBROUTINE hfx_ri_update_ks_mo(qs_env, ri_data, ks_matrix, mo_coeff, &
1190 homo, geometry_did_change, nspins, fac)
1191 TYPE(qs_environment_type), POINTER :: qs_env
1192 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
1193 TYPE(dbcsr_p_type), DIMENSION(:, :) :: ks_matrix
1194 TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: mo_coeff
1195 INTEGER, DIMENSION(:) :: homo
1196 LOGICAL, INTENT(IN) :: geometry_did_change
1197 INTEGER, INTENT(IN) :: nspins
1198 REAL(dp), INTENT(IN) :: fac
1199
1200 CHARACTER(LEN=*), PARAMETER :: routinen = 'hfx_ri_update_ks_mo'
1201
1202 INTEGER :: bsize, bsum, comm_2d_handle, handle, &
1203 handle2, i_mem, iblock, iproc, ispin, &
1204 n_mem, n_mos, nblock, unit_nr_dbcsr
1205 INTEGER(int_8) :: nblks, nflop
1206 INTEGER, ALLOCATABLE, DIMENSION(:) :: batch_ranges_1, batch_ranges_2, dist1, dist2, dist3, &
1207 mem_end, mem_end_block_1, mem_end_block_2, mem_size, mem_start, mem_start_block_1, &
1208 mem_start_block_2, mo_bsizes_1, mo_bsizes_2
1209 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: bounds
1210 INTEGER, DIMENSION(2) :: pdims_2d
1211 INTEGER, DIMENSION(3) :: pdims
1212 LOGICAL :: do_initialize
1213 REAL(dp) :: t1, t2
1214 TYPE(dbcsr_distribution_type) :: ks_dist
1215 TYPE(dbt_pgrid_type) :: pgrid, pgrid_2d
1216 TYPE(dbt_type) :: ks_t, ks_t_mat, mo_coeff_t, &
1217 mo_coeff_t_split
1218 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_int_mo_1, t_3c_int_mo_2
1219 TYPE(mp_comm_type) :: comm_2d
1220 TYPE(mp_para_env_type), POINTER :: para_env
1221
1222 CALL timeset(routinen, handle)
1223
1224 cpassert(SIZE(ks_matrix, 2) == 1)
1225
1226 unit_nr_dbcsr = ri_data%unit_nr_dbcsr
1227
1228 IF (geometry_did_change) THEN
1229 CALL hfx_ri_pre_scf_mo(qs_env, ri_data, nspins)
1230 END IF
1231
1232 nblks = dbt_get_num_blocks_total(ri_data%t_3c_int_ctr_1(1, 1))
1233 IF (nblks == 0) THEN
1234 cpabort("3-center integrals are not available (first call requires geometry_did_change=.TRUE.)")
1235 END IF
1236
1237 DO ispin = 1, nspins
1238 nblks = dbt_get_num_blocks_total(ri_data%t_2c_int(ispin, 1))
1239 IF (nblks == 0) THEN
1240 cpabort("2-center integrals are not available (first call requires geometry_did_change=.TRUE.)")
1241 END IF
1242 END DO
1243
1244 IF (.NOT. ALLOCATED(ri_data%t_3c_int_mo)) THEN
1245 do_initialize = .true.
1246 cpassert(.NOT. ALLOCATED(ri_data%t_3c_ctr_RI))
1247 cpassert(.NOT. ALLOCATED(ri_data%t_3c_ctr_KS))
1248 cpassert(.NOT. ALLOCATED(ri_data%t_3c_ctr_KS_copy))
1249 ALLOCATE (ri_data%t_3c_int_mo(nspins, 1, 1))
1250 ALLOCATE (ri_data%t_3c_ctr_RI(nspins, 1, 1))
1251 ALLOCATE (ri_data%t_3c_ctr_KS(nspins, 1, 1))
1252 ALLOCATE (ri_data%t_3c_ctr_KS_copy(nspins, 1, 1))
1253 ELSE
1254 do_initialize = .false.
1255 END IF
1256
1257 CALL get_qs_env(qs_env, para_env=para_env)
1258
1259 ALLOCATE (bounds(2, 1))
1260
1261 CALL dbcsr_get_info(ks_matrix(1, 1)%matrix, distribution=ks_dist)
1262 CALL dbcsr_distribution_get(ks_dist, group=comm_2d_handle, nprows=pdims_2d(1), npcols=pdims_2d(2))
1263
1264 CALL comm_2d%set_handle(comm_2d_handle)
1265 pgrid_2d = dbt_nd_mp_comm(comm_2d, [1], [2], pdims_2d=pdims_2d)
1266
1267 CALL create_2c_tensor(ks_t, dist1, dist2, pgrid_2d, ri_data%bsizes_AO_fit, ri_data%bsizes_AO_fit, &
1268 name="(AO | AO)")
1269
1270 DEALLOCATE (dist1, dist2)
1271
1272 CALL para_env%sync()
1273 t1 = m_walltime()
1274
1275 ALLOCATE (t_3c_int_mo_1(1, 1), t_3c_int_mo_2(1, 1))
1276 DO ispin = 1, nspins
1277
1278 CALL dbcsr_get_info(mo_coeff(ispin), nfullcols_total=n_mos)
1279 ALLOCATE (mo_bsizes_2(n_mos))
1280 mo_bsizes_2 = 1
1281
1282 CALL create_tensor_batches(mo_bsizes_2, ri_data%n_mem, mem_start, mem_end, &
1283 mem_start_block_2, mem_end_block_2)
1284 n_mem = ri_data%n_mem
1285 ALLOCATE (mem_size(n_mem))
1286
1287 DO i_mem = 1, n_mem
1288 bsize = sum(mo_bsizes_2(mem_start_block_2(i_mem):mem_end_block_2(i_mem)))
1289 mem_size(i_mem) = bsize
1290 END DO
1291
1292 CALL split_block_sizes(mem_size, mo_bsizes_1, ri_data%max_bsize_MO)
1293 ALLOCATE (mem_start_block_1(n_mem))
1294 ALLOCATE (mem_end_block_1(n_mem))
1295 nblock = SIZE(mo_bsizes_1)
1296 iblock = 0
1297 DO i_mem = 1, n_mem
1298 bsum = 0
1299 DO
1300 iblock = iblock + 1
1301 cpassert(iblock <= nblock)
1302 bsum = bsum + mo_bsizes_1(iblock)
1303 IF (bsum == mem_size(i_mem)) THEN
1304 IF (i_mem == 1) THEN
1305 mem_start_block_1(i_mem) = 1
1306 ELSE
1307 mem_start_block_1(i_mem) = mem_end_block_1(i_mem - 1) + 1
1308 END IF
1309 mem_end_block_1(i_mem) = iblock
1310 EXIT
1311 END IF
1312 END DO
1313 END DO
1314
1315 ALLOCATE (batch_ranges_1(ri_data%n_mem + 1))
1316 batch_ranges_1(:ri_data%n_mem) = mem_start_block_1(:)
1317 batch_ranges_1(ri_data%n_mem + 1) = mem_end_block_1(ri_data%n_mem) + 1
1318
1319 ALLOCATE (batch_ranges_2(ri_data%n_mem + 1))
1320 batch_ranges_2(:ri_data%n_mem) = mem_start_block_2(:)
1321 batch_ranges_2(ri_data%n_mem + 1) = mem_end_block_2(ri_data%n_mem) + 1
1322
1323 iproc = para_env%mepos
1324
1325 CALL create_3c_tensor(t_3c_int_mo_1(1, 1), dist1, dist2, dist3, ri_data%pgrid_1, &
1326 ri_data%bsizes_AO_split, ri_data%bsizes_RI_split, mo_bsizes_1, &
1327 [1, 2], [3], &
1328 name="(AO RI | MO)")
1329
1330 DEALLOCATE (dist1, dist2, dist3)
1331
1332 CALL create_3c_tensor(t_3c_int_mo_2(1, 1), dist1, dist2, dist3, ri_data%pgrid_2, &
1333 mo_bsizes_1, ri_data%bsizes_RI_split, ri_data%bsizes_AO_split, &
1334 [1], [2, 3], &
1335 name="(MO | RI AO)")
1336
1337 DEALLOCATE (dist1, dist2, dist3)
1338
1339 CALL create_2c_tensor(mo_coeff_t_split, dist1, dist2, pgrid_2d, ri_data%bsizes_AO_split, mo_bsizes_1, &
1340 name="(AO | MO)")
1341
1342 DEALLOCATE (dist1, dist2)
1343
1344 cpassert(homo(ispin)/ri_data%n_mem > 0)
1345
1346 IF (do_initialize) THEN
1347 pdims(:) = 0
1348
1349 CALL dbt_pgrid_create(para_env, pdims, pgrid, &
1350 tensor_dims=[SIZE(ri_data%bsizes_RI_fit), &
1351 (homo(ispin) - 1)/ri_data%n_mem + 1, &
1352 SIZE(ri_data%bsizes_AO_fit)])
1353 CALL create_3c_tensor(ri_data%t_3c_int_mo(ispin, 1, 1), dist1, dist2, dist3, pgrid, &
1354 ri_data%bsizes_RI_fit, mo_bsizes_2, ri_data%bsizes_AO_fit, &
1355 [1], [2, 3], &
1356 name="(RI | MO AO)")
1357
1358 DEALLOCATE (dist1, dist2, dist3)
1359
1360 CALL create_3c_tensor(ri_data%t_3c_ctr_KS(ispin, 1, 1), dist1, dist2, dist3, pgrid, &
1361 ri_data%bsizes_RI_fit, mo_bsizes_2, ri_data%bsizes_AO_fit, &
1362 [1, 2], [3], &
1363 name="(RI MO | AO)")
1364 DEALLOCATE (dist1, dist2, dist3)
1365 CALL dbt_pgrid_destroy(pgrid)
1366
1367 CALL dbt_create(ri_data%t_3c_int_mo(ispin, 1, 1), ri_data%t_3c_ctr_RI(ispin, 1, 1), name="(RI | MO AO)")
1368 CALL dbt_create(ri_data%t_3c_ctr_KS(ispin, 1, 1), ri_data%t_3c_ctr_KS_copy(ispin, 1, 1))
1369 END IF
1370
1371 CALL dbt_create(mo_coeff(ispin), mo_coeff_t, name="MO coeffs")
1372 CALL dbt_copy_matrix_to_tensor(mo_coeff(ispin), mo_coeff_t)
1373 CALL dbt_copy(mo_coeff_t, mo_coeff_t_split, move_data=.true.)
1374 CALL dbt_filter(mo_coeff_t_split, ri_data%filter_eps_mo)
1375 CALL dbt_destroy(mo_coeff_t)
1376
1377 CALL dbt_batched_contract_init(ks_t)
1378 CALL dbt_batched_contract_init(ri_data%t_3c_ctr_KS(ispin, 1, 1), batch_range_2=batch_ranges_2)
1379 CALL dbt_batched_contract_init(ri_data%t_3c_ctr_KS_copy(ispin, 1, 1), batch_range_2=batch_ranges_2)
1380
1381 CALL dbt_batched_contract_init(ri_data%t_2c_int(ispin, 1))
1382 CALL dbt_batched_contract_init(ri_data%t_3c_int_mo(ispin, 1, 1), batch_range_2=batch_ranges_2)
1383 CALL dbt_batched_contract_init(ri_data%t_3c_ctr_RI(ispin, 1, 1), batch_range_2=batch_ranges_2)
1384
1385 DO i_mem = 1, n_mem
1386
1387 bounds(:, 1) = [mem_start(i_mem), mem_end(i_mem)]
1388
1389 CALL dbt_batched_contract_init(mo_coeff_t_split)
1390 CALL dbt_batched_contract_init(ri_data%t_3c_int_ctr_1(1, 1))
1391 CALL dbt_batched_contract_init(t_3c_int_mo_1(1, 1), &
1392 batch_range_3=batch_ranges_1)
1393 CALL timeset(routinen//"_MOx3C_R", handle2)
1394 CALL dbt_contract(1.0_dp, mo_coeff_t_split, ri_data%t_3c_int_ctr_1(1, 1), &
1395 0.0_dp, t_3c_int_mo_1(1, 1), &
1396 contract_1=[1], notcontract_1=[2], &
1397 contract_2=[3], notcontract_2=[1, 2], &
1398 map_1=[3], map_2=[1, 2], &
1399 bounds_2=bounds, &
1400 filter_eps=ri_data%filter_eps_mo/2, &
1401 unit_nr=unit_nr_dbcsr, &
1402 move_data=.false., &
1403 flop=nflop)
1404
1405 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1406
1407 CALL timestop(handle2)
1408 CALL dbt_batched_contract_finalize(mo_coeff_t_split)
1409 CALL dbt_batched_contract_finalize(ri_data%t_3c_int_ctr_1(1, 1))
1410 CALL dbt_batched_contract_finalize(t_3c_int_mo_1(1, 1))
1411
1412 CALL timeset(routinen//"_copy_1", handle2)
1413 CALL dbt_copy(t_3c_int_mo_1(1, 1), ri_data%t_3c_int_mo(ispin, 1, 1), order=[3, 1, 2], move_data=.true.)
1414 CALL timestop(handle2)
1415
1416 CALL dbt_batched_contract_init(mo_coeff_t_split)
1417 CALL dbt_batched_contract_init(ri_data%t_3c_int_ctr_2(1, 1))
1418 CALL dbt_batched_contract_init(t_3c_int_mo_2(1, 1), &
1419 batch_range_1=batch_ranges_1)
1420
1421 CALL timeset(routinen//"_MOx3C_L", handle2)
1422 CALL dbt_contract(1.0_dp, mo_coeff_t_split, ri_data%t_3c_int_ctr_2(1, 1), &
1423 0.0_dp, t_3c_int_mo_2(1, 1), &
1424 contract_1=[1], notcontract_1=[2], &
1425 contract_2=[1], notcontract_2=[2, 3], &
1426 map_1=[1], map_2=[2, 3], &
1427 bounds_2=bounds, &
1428 filter_eps=ri_data%filter_eps_mo/2, &
1429 unit_nr=unit_nr_dbcsr, &
1430 move_data=.false., &
1431 flop=nflop)
1432
1433 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1434
1435 CALL timestop(handle2)
1436
1437 CALL dbt_batched_contract_finalize(mo_coeff_t_split)
1438 CALL dbt_batched_contract_finalize(ri_data%t_3c_int_ctr_2(1, 1))
1439 CALL dbt_batched_contract_finalize(t_3c_int_mo_2(1, 1))
1440
1441 CALL timeset(routinen//"_copy_1", handle2)
1442 CALL dbt_copy(t_3c_int_mo_2(1, 1), ri_data%t_3c_int_mo(ispin, 1, 1), order=[2, 1, 3], &
1443 summation=.true., move_data=.true.)
1444
1445 CALL dbt_filter(ri_data%t_3c_int_mo(ispin, 1, 1), ri_data%filter_eps_mo)
1446 CALL timestop(handle2)
1447
1448 CALL timeset(routinen//"_RIx3C", handle2)
1449
1450 CALL dbt_contract(1.0_dp, ri_data%t_2c_int(ispin, 1), ri_data%t_3c_int_mo(ispin, 1, 1), &
1451 0.0_dp, ri_data%t_3c_ctr_RI(ispin, 1, 1), &
1452 contract_1=[1], notcontract_1=[2], &
1453 contract_2=[1], notcontract_2=[2, 3], &
1454 map_1=[1], map_2=[2, 3], filter_eps=ri_data%filter_eps, &
1455 unit_nr=unit_nr_dbcsr, &
1456 flop=nflop)
1457
1458 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1459
1460 CALL timestop(handle2)
1461
1462 CALL timeset(routinen//"_copy_2", handle2)
1463
1464 ! note: this copy should not involve communication (same block sizes, same 3d distribution on same process grid)
1465 CALL dbt_copy(ri_data%t_3c_ctr_RI(ispin, 1, 1), ri_data%t_3c_ctr_KS(ispin, 1, 1), move_data=.true.)
1466 CALL dbt_copy(ri_data%t_3c_ctr_KS(ispin, 1, 1), ri_data%t_3c_ctr_KS_copy(ispin, 1, 1))
1467 CALL timestop(handle2)
1468
1469 CALL timeset(routinen//"_3Cx3C", handle2)
1470 CALL dbt_contract(-fac, ri_data%t_3c_ctr_KS(ispin, 1, 1), ri_data%t_3c_ctr_KS_copy(ispin, 1, 1), &
1471 1.0_dp, ks_t, &
1472 contract_1=[1, 2], notcontract_1=[3], &
1473 contract_2=[1, 2], notcontract_2=[3], &
1474 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps/n_mem, &
1475 unit_nr=unit_nr_dbcsr, move_data=.true., &
1476 flop=nflop)
1477
1478 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1479
1480 CALL timestop(handle2)
1481 END DO
1482
1483 CALL dbt_batched_contract_finalize(ks_t)
1484 CALL dbt_batched_contract_finalize(ri_data%t_3c_ctr_KS(ispin, 1, 1))
1485 CALL dbt_batched_contract_finalize(ri_data%t_3c_ctr_KS_copy(ispin, 1, 1))
1486
1487 CALL dbt_batched_contract_finalize(ri_data%t_2c_int(ispin, 1))
1488 CALL dbt_batched_contract_finalize(ri_data%t_3c_int_mo(ispin, 1, 1))
1489 CALL dbt_batched_contract_finalize(ri_data%t_3c_ctr_RI(ispin, 1, 1))
1490
1491 CALL dbt_destroy(t_3c_int_mo_1(1, 1))
1492 CALL dbt_destroy(t_3c_int_mo_2(1, 1))
1493 CALL dbt_clear(ri_data%t_3c_int_mo(ispin, 1, 1))
1494
1495 CALL dbt_destroy(mo_coeff_t_split)
1496
1497 CALL dbt_filter(ks_t, ri_data%filter_eps)
1498
1499 CALL dbt_create(ks_matrix(ispin, 1)%matrix, ks_t_mat)
1500 CALL dbt_copy(ks_t, ks_t_mat, move_data=.true.)
1501 CALL dbt_copy_tensor_to_matrix(ks_t_mat, ks_matrix(ispin, 1)%matrix, summation=.true.)
1502 CALL dbt_destroy(ks_t_mat)
1503
1504 DEALLOCATE (mem_end, mem_start, mo_bsizes_2, mem_size, mem_start_block_1, mem_end_block_1, &
1505 mem_start_block_2, mem_end_block_2, batch_ranges_1, batch_ranges_2)
1506
1507 END DO
1508
1509 CALL dbt_pgrid_destroy(pgrid_2d)
1510 CALL dbt_destroy(ks_t)
1511
1512 CALL para_env%sync()
1513 t2 = m_walltime()
1514
1515 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
1516
1517 CALL timestop(handle)
1518
1519 END SUBROUTINE
1520
1521! **************************************************************************************************
1522!> \brief Calculate Fock (AKA Kohn-Sham) matrix in rho flavor
1523!>
1524!> M(mu, lambda, R) = sum_{nu} int_3c(mu, nu, R) P(nu, lambda)
1525!> KS(mu, lambda) = sum_{nu,R} B(mu, nu, R) M(lambda, nu, R)
1526!> \param qs_env ...
1527!> \param ri_data ...
1528!> \param ks_matrix ...
1529!> \param rho_ao ...
1530!> \param geometry_did_change ...
1531!> \param nspins ...
1532!> \param fac ...
1533! **************************************************************************************************
1534 SUBROUTINE hfx_ri_update_ks_pmat(qs_env, ri_data, ks_matrix, rho_ao, &
1535 geometry_did_change, nspins, fac)
1536 TYPE(qs_environment_type), POINTER :: qs_env
1537 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
1538 TYPE(dbcsr_p_type), DIMENSION(:, :) :: ks_matrix, rho_ao
1539 LOGICAL, INTENT(IN) :: geometry_did_change
1540 INTEGER, INTENT(IN) :: nspins
1541 REAL(dp), INTENT(IN) :: fac
1542
1543 CHARACTER(LEN=*), PARAMETER :: routinen = 'hfx_ri_update_ks_Pmat'
1544
1545 INTEGER :: handle, handle2, i_mem, ispin, j_mem, &
1546 n_mem, n_mem_ri, unit_nr, unit_nr_dbcsr
1547 INTEGER(int_8) :: flops_ks_max, flops_p_max, nblks, nflop, &
1548 nze, nze_3c, nze_3c_1, nze_3c_2, &
1549 nze_ks, nze_rho
1550 INTEGER, ALLOCATABLE, DIMENSION(:) :: batch_ranges_ao, batch_ranges_ri, dist1, &
1551 dist2
1552 INTEGER, DIMENSION(2, 1) :: bounds_i
1553 INTEGER, DIMENSION(2, 2) :: bounds_ij, bounds_j
1554 INTEGER, DIMENSION(3) :: dims_3c
1555 REAL(dp) :: memory_3c, occ, occ_3c, occ_3c_1, &
1556 occ_3c_2, occ_ks, occ_rho, t1, t2, &
1557 unused
1558 TYPE(dbt_type) :: ks_t, ks_tmp, rho_ao_tmp, t_3c_1, &
1559 t_3c_3, tensor_old
1560 TYPE(mp_para_env_type), POINTER :: para_env
1561
1562 IF (.NOT. fac > epsilon(0.0_dp)) RETURN
1563
1564 CALL timeset(routinen, handle)
1565
1566 NULLIFY (para_env)
1567
1568 ! get a useful output_unit
1569 unit_nr_dbcsr = ri_data%unit_nr_dbcsr
1570 unit_nr = ri_data%unit_nr
1571
1572 CALL get_qs_env(qs_env, para_env=para_env)
1573
1574 cpassert(SIZE(ks_matrix, 2) == 1)
1575
1576 IF (geometry_did_change) THEN
1577 CALL hfx_ri_pre_scf_pmat(qs_env, ri_data)
1578 DO ispin = 1, nspins
1579 CALL dbt_scale(ri_data%rho_ao_t(ispin, 1), 0.0_dp)
1580 CALL dbt_scale(ri_data%ks_t(ispin, 1), 0.0_dp)
1581 END DO
1582 END IF
1583
1584 nblks = dbt_get_num_blocks_total(ri_data%t_3c_int_ctr_2(1, 1))
1585 IF (nblks == 0) THEN
1586 cpabort("3-center integrals are not available (first call requires geometry_did_change=.TRUE.)")
1587 END IF
1588
1589 n_mem = ri_data%n_mem
1590 n_mem_ri = ri_data%n_mem_RI
1591
1592 CALL dbt_create(ks_matrix(1, 1)%matrix, ks_tmp)
1593 CALL dbt_create(rho_ao(1, 1)%matrix, rho_ao_tmp)
1594
1595 CALL create_2c_tensor(ks_t, dist1, dist2, ri_data%pgrid_2d, &
1596 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
1597 name="(AO | AO)")
1598 DEALLOCATE (dist1, dist2)
1599
1600 CALL dbt_create(ri_data%t_3c_int_ctr_2(1, 1), t_3c_1)
1601 CALL dbt_create(ri_data%t_3c_int_ctr_1(1, 1), t_3c_3)
1602
1603 CALL para_env%sync()
1604 t1 = m_walltime()
1605
1606 flops_ks_max = 0; flops_p_max = 0
1607
1608 ALLOCATE (batch_ranges_ri(ri_data%n_mem_RI + 1))
1609 ALLOCATE (batch_ranges_ao(ri_data%n_mem + 1))
1610 batch_ranges_ri(:ri_data%n_mem_RI) = ri_data%starts_array_RI_mem_block(:)
1611 batch_ranges_ri(ri_data%n_mem_RI + 1) = ri_data%ends_array_RI_mem_block(ri_data%n_mem_RI) + 1
1612 batch_ranges_ao(:ri_data%n_mem) = ri_data%starts_array_mem_block(:)
1613 batch_ranges_ao(ri_data%n_mem + 1) = ri_data%ends_array_mem_block(ri_data%n_mem) + 1
1614
1615 memory_3c = 0.0_dp
1616 DO ispin = 1, nspins
1617
1618 CALL get_tensor_occupancy(ri_data%t_3c_int_ctr_2(1, 1), nze_3c, occ_3c)
1619
1620 nze_rho = 0
1621 occ_rho = 0.0_dp
1622 nze_3c_1 = 0
1623 occ_3c_1 = 0.0_dp
1624 nze_3c_2 = 0
1625 occ_3c_2 = 0.0_dp
1626
1627 CALL dbt_copy_matrix_to_tensor(rho_ao(ispin, 1)%matrix, rho_ao_tmp)
1628
1629 !We work with Delta P: the diff between previous SCF step and this one, for increased sparsity
1630 CALL dbt_scale(ri_data%rho_ao_t(ispin, 1), -1.0_dp)
1631 CALL dbt_copy(rho_ao_tmp, ri_data%rho_ao_t(ispin, 1), summation=.true., move_data=.true.)
1632
1633 CALL get_tensor_occupancy(ri_data%rho_ao_t(ispin, 1), nze_rho, occ_rho)
1634
1635 CALL dbt_batched_contract_init(ri_data%t_3c_int_ctr_1(1, 1), batch_range_1=batch_ranges_ao, &
1636 batch_range_2=batch_ranges_ri)
1637 CALL dbt_batched_contract_init(t_3c_3, batch_range_1=batch_ranges_ao, batch_range_2=batch_ranges_ri)
1638
1639 CALL dbt_create(ri_data%t_3c_int_ctr_1(1, 1), tensor_old)
1640
1641 DO i_mem = 1, n_mem
1642
1643 CALL dbt_batched_contract_init(ri_data%rho_ao_t(ispin, 1))
1644 CALL dbt_batched_contract_init(ri_data%t_3c_int_ctr_2(1, 1), batch_range_2=batch_ranges_ri, &
1645 batch_range_3=batch_ranges_ao)
1646 CALL dbt_batched_contract_init(t_3c_1, batch_range_2=batch_ranges_ri, batch_range_3=batch_ranges_ao)
1647 DO j_mem = 1, n_mem_ri
1648
1649 CALL timeset(routinen//"_Px3C", handle2)
1650
1651 CALL dbt_get_info(t_3c_1, nfull_total=dims_3c)
1652 bounds_i(:, 1) = [ri_data%starts_array_mem(i_mem), ri_data%ends_array_mem(i_mem)]
1653 bounds_j(:, 1) = [1, dims_3c(1)]
1654 bounds_j(:, 2) = [ri_data%starts_array_RI_mem(j_mem), ri_data%ends_array_RI_mem(j_mem)]
1655
1656 CALL dbt_contract(1.0_dp, ri_data%rho_ao_t(ispin, 1), ri_data%t_3c_int_ctr_2(1, 1), &
1657 0.0_dp, t_3c_1, &
1658 contract_1=[2], notcontract_1=[1], &
1659 contract_2=[3], notcontract_2=[1, 2], &
1660 map_1=[3], map_2=[1, 2], filter_eps=ri_data%filter_eps, &
1661 bounds_2=bounds_i, &
1662 bounds_3=bounds_j, &
1663 unit_nr=unit_nr_dbcsr, &
1664 flop=nflop)
1665
1666 CALL timestop(handle2)
1667
1668 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1669
1670 CALL get_tensor_occupancy(t_3c_1, nze, occ)
1671 nze_3c_1 = nze_3c_1 + nze
1672 occ_3c_1 = occ_3c_1 + occ
1673
1674 CALL timeset(routinen//"_copy_2", handle2)
1675 CALL dbt_copy(t_3c_1, t_3c_3, order=[3, 2, 1], move_data=.true.)
1676 CALL timestop(handle2)
1677
1678 bounds_ij(:, 1) = [ri_data%starts_array_mem(i_mem), ri_data%ends_array_mem(i_mem)]
1679 bounds_ij(:, 2) = [ri_data%starts_array_RI_mem(j_mem), ri_data%ends_array_RI_mem(j_mem)]
1680
1681 CALL decompress_tensor(tensor_old, ri_data%blk_indices(i_mem, j_mem)%ind, &
1682 ri_data%store_3c(i_mem, j_mem), ri_data%filter_eps_storage)
1683
1684 CALL dbt_copy(tensor_old, ri_data%t_3c_int_ctr_1(1, 1), move_data=.true.)
1685
1686 CALL get_tensor_occupancy(ri_data%t_3c_int_ctr_1(1, 1), nze, occ)
1687 nze_3c_2 = nze_3c_2 + nze
1688 occ_3c_2 = occ_3c_2 + occ
1689 CALL timeset(routinen//"_KS", handle2)
1690 CALL dbt_batched_contract_init(ks_t)
1691 CALL dbt_contract(-fac, ri_data%t_3c_int_ctr_1(1, 1), t_3c_3, &
1692 1.0_dp, ks_t, &
1693 contract_1=[1, 2], notcontract_1=[3], &
1694 contract_2=[1, 2], notcontract_2=[3], &
1695 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps/n_mem, &
1696 bounds_1=bounds_ij, &
1697 unit_nr=unit_nr_dbcsr, &
1698 flop=nflop, move_data=.true.)
1699
1700 CALL dbt_batched_contract_finalize(ks_t, unit_nr=unit_nr_dbcsr)
1701 CALL timestop(handle2)
1702
1703 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1704
1705 END DO
1706 CALL dbt_batched_contract_finalize(ri_data%rho_ao_t(ispin, 1), unit_nr=unit_nr_dbcsr)
1707 CALL dbt_batched_contract_finalize(ri_data%t_3c_int_ctr_2(1, 1))
1708 CALL dbt_batched_contract_finalize(t_3c_1)
1709 END DO
1710 CALL dbt_batched_contract_finalize(ri_data%t_3c_int_ctr_1(1, 1))
1711 CALL dbt_batched_contract_finalize(t_3c_3)
1712
1713 DO i_mem = 1, n_mem
1714 DO j_mem = 1, n_mem_ri
1715 associate(blk_indices => ri_data%blk_indices(i_mem, j_mem), t_3c => ri_data%t_3c_int_ctr_1(1, 1))
1716 CALL decompress_tensor(tensor_old, blk_indices%ind, &
1717 ri_data%store_3c(i_mem, j_mem), ri_data%filter_eps_storage)
1718 CALL dbt_copy(tensor_old, t_3c, move_data=.true.)
1719
1720 unused = 0
1721 CALL compress_tensor(t_3c, blk_indices%ind, ri_data%store_3c(i_mem, j_mem), &
1722 ri_data%filter_eps_storage, unused)
1723 END associate
1724 END DO
1725 END DO
1726
1727 CALL dbt_destroy(tensor_old)
1728
1729 CALL get_tensor_occupancy(ks_t, nze_ks, occ_ks)
1730
1731 !rho_ao_t holds the density difference, and ks_t is built upon it => need the full picture
1732 CALL dbt_copy_matrix_to_tensor(rho_ao(ispin, 1)%matrix, rho_ao_tmp)
1733 CALL dbt_copy(rho_ao_tmp, ri_data%rho_ao_t(ispin, 1), move_data=.true.)
1734 CALL dbt_copy(ks_t, ri_data%ks_t(ispin, 1), summation=.true., move_data=.true.)
1735
1736 CALL dbt_copy(ri_data%ks_t(ispin, 1), ks_tmp)
1737 CALL dbt_copy_tensor_to_matrix(ks_tmp, ks_matrix(ispin, 1)%matrix, summation=.true.)
1738 CALL dbt_clear(ks_tmp)
1739
1740 IF (unit_nr > 0 .AND. geometry_did_change) THEN
1741 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1742 'Occupancy of density matrix P:', real(nze_rho, dp), '/', occ_rho*100, '%'
1743 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1744 'Occupancy of 3c ints:', real(nze_3c, dp), '/', occ_3c*100, '%'
1745 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1746 'Occupancy after contraction with K:', real(nze_3c_2, dp), '/', occ_3c_2*100, '%'
1747 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1748 'Occupancy after contraction with P:', real(nze_3c_1, dp), '/', occ_3c_1*100, '%'
1749 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1750 'Occupancy of Kohn-Sham matrix:', real(nze_ks, dp), '/', occ_ks*100, '%'
1751 END IF
1752
1753 END DO
1754
1755 CALL para_env%sync()
1756 t2 = m_walltime()
1757
1758 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
1759
1760 CALL dbt_destroy(t_3c_1)
1761 CALL dbt_destroy(t_3c_3)
1762
1763 CALL dbt_destroy(rho_ao_tmp)
1764 CALL dbt_destroy(ks_t)
1765 CALL dbt_destroy(ks_tmp)
1766
1767 CALL timestop(handle)
1768
1769 END SUBROUTINE
1770
1771! **************************************************************************************************
1772!> \brief Implementation based on the MO flavor
1773!> \param qs_env ...
1774!> \param ri_data ...
1775!> \param nspins ...
1776!> \param hf_fraction ...
1777!> \param mo_coeff ...
1778!> \param use_virial ...
1779!> \note There is no response code for forces with the MO flavor
1780! **************************************************************************************************
1781 SUBROUTINE hfx_ri_forces_mo(qs_env, ri_data, nspins, hf_fraction, mo_coeff, use_virial)
1782
1783 TYPE(qs_environment_type), POINTER :: qs_env
1784 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
1785 INTEGER, INTENT(IN) :: nspins
1786 REAL(dp), INTENT(IN) :: hf_fraction
1787 TYPE(dbcsr_type), DIMENSION(:), INTENT(IN) :: mo_coeff
1788 LOGICAL, INTENT(IN), OPTIONAL :: use_virial
1789
1790 CHARACTER(LEN=*), PARAMETER :: routinen = 'hfx_ri_forces_mo'
1791
1792 INTEGER :: dummy_int, handle, i_mem, i_xyz, ibasis, ispin, j_mem, k_mem, n_mem, n_mem_input, &
1793 n_mem_input_ri, n_mem_ri, n_mem_ri_fit, n_mos, natom, nkind, unit_nr_dbcsr
1794 INTEGER(int_8) :: nflop
1795 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, batch_blk_end, batch_blk_start, &
1796 batch_end, batch_end_ri, batch_end_ri_fit, batch_ranges, batch_ranges_ri, &
1797 batch_ranges_ri_fit, batch_start, batch_start_ri, batch_start_ri_fit, bsizes_mo, dist1, &
1798 dist2, dist3, idx_to_at_ao, idx_to_at_ri, kind_of
1799 INTEGER, DIMENSION(2, 1) :: bounds_ctr_1d
1800 INTEGER, DIMENSION(2, 2) :: bounds_ctr_2d
1801 INTEGER, DIMENSION(3) :: pdims
1802 LOGICAL :: use_virial_prv
1803 REAL(dp) :: pref, spin_fac, t1, t2
1804 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1805 TYPE(block_ind_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_der_ao_ind, t_3c_der_ri_ind
1806 TYPE(cell_type), POINTER :: cell
1807 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1808 TYPE(dbt_pgrid_type) :: pgrid_1, pgrid_2
1809 TYPE(dbt_type) :: t_2c_ri, t_2c_ri_inv, t_2c_ri_met, t_2c_ri_pq, t_2c_tmp, t_3c_0, t_3c_1, &
1810 t_3c_2, t_3c_3, t_3c_4, t_3c_5, t_3c_6, t_3c_ao_ri_ao, t_3c_ao_ri_mo, t_3c_desymm, &
1811 t_3c_mo_ri_ao, t_3c_mo_ri_mo, t_3c_ri_ao_ao, t_3c_ri_ctr, t_3c_ri_mo_mo, &
1812 t_3c_ri_mo_mo_fit, t_3c_work, t_mo_coeff, t_mo_cpy
1813 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_2c_der_metric, t_2c_der_ri, t_2c_mo_ao, &
1814 t_2c_mo_ao_ctr, t_3c_der_ao, t_3c_der_ao_ctr_1, t_3c_der_ri, t_3c_der_ri_ctr_1, &
1815 t_3c_der_ri_ctr_2
1816 TYPE(dft_control_type), POINTER :: dft_control
1817 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1818 DIMENSION(:), TARGET :: basis_set_ao, basis_set_ri
1819 TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
1820 TYPE(hfx_compression_type), ALLOCATABLE, &
1821 DIMENSION(:, :) :: t_3c_der_ao_comp, t_3c_der_ri_comp
1822 TYPE(mp_para_env_type), POINTER :: para_env
1823 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1824 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1825 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1826
1827 ! 1) Precompute the derivatives that are needed (3c, 3c RI and metric)
1828 ! 2) Go over batched of occupied MOs so as to save memory and optimize contractions
1829 ! 3) Contract all 3c integrals and derivatives with MO coeffs
1830 ! 4) Contract relevant quantities with the inverse 2c RI (metric or pot)
1831 ! 5) First force contribution with the 2c RI derivative d/dx (Q|R)
1832 ! 6) If metric, do the additional contraction with S_pq^-1 (Q|R)
1833 ! 7) Do the force contribution due to 3c integrals (a'b|P) and (ab|P')
1834 ! 8) If metric, do the last force contribution due to d/dx S^-1 (First contract (ab|P), then S^-1)
1835
1836 use_virial_prv = .false.
1837 IF (PRESENT(use_virial)) use_virial_prv = use_virial
1838 IF (use_virial_prv) THEN
1839 cpabort("Stress tensor with RI-HFX MO flavor not implemented.")
1840 END IF
1841
1842 unit_nr_dbcsr = ri_data%unit_nr_dbcsr
1843
1844 CALL get_qs_env(qs_env, natom=natom, particle_set=particle_set, nkind=nkind, &
1845 atomic_kind_set=atomic_kind_set, cell=cell, force=force, &
1846 matrix_s=matrix_s, para_env=para_env, dft_control=dft_control, &
1847 qs_kind_set=qs_kind_set)
1848
1849 pdims(:) = 0
1850 CALL dbt_pgrid_create(para_env, pdims, pgrid_1, tensor_dims=[SIZE(ri_data%bsizes_AO_split), &
1851 SIZE(ri_data%bsizes_RI_split), &
1852 SIZE(ri_data%bsizes_AO_split)])
1853 pdims(:) = 0
1854 CALL dbt_pgrid_create(para_env, pdims, pgrid_2, tensor_dims=[SIZE(ri_data%bsizes_RI_split), &
1855 SIZE(ri_data%bsizes_AO_split), &
1856 SIZE(ri_data%bsizes_AO_split)])
1857
1858 CALL create_3c_tensor(t_3c_ao_ri_ao, dist1, dist2, dist3, pgrid_1, &
1859 ri_data%bsizes_AO_split, ri_data%bsizes_RI_split, ri_data%bsizes_AO_split, &
1860 [1, 2], [3], name="(AO RI | AO)")
1861 DEALLOCATE (dist1, dist2, dist3)
1862 CALL create_3c_tensor(t_3c_ri_ao_ao, dist1, dist2, dist3, pgrid_2, &
1863 ri_data%bsizes_RI_split, ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
1864 [1], [2, 3], name="(RI | AO AO)")
1865 DEALLOCATE (dist1, dist2, dist3)
1866
1867 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
1868 CALL basis_set_list_setup(basis_set_ri, ri_data%ri_basis_type, qs_kind_set)
1869 CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_ri)
1870 CALL basis_set_list_setup(basis_set_ao, ri_data%orb_basis_type, qs_kind_set)
1871 CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_ao)
1872
1873 DO ibasis = 1, SIZE(basis_set_ao)
1874 orb_basis => basis_set_ao(ibasis)%gto_basis_set
1875 CALL init_interaction_radii_orb_basis(orb_basis, ri_data%eps_pgf_orb)
1876 ri_basis => basis_set_ri(ibasis)%gto_basis_set
1877 CALL init_interaction_radii_orb_basis(ri_basis, ri_data%eps_pgf_orb)
1878 END DO
1879
1880 ALLOCATE (t_2c_der_metric(3), t_2c_der_ri(3), t_2c_mo_ao(3), t_2c_mo_ao_ctr(3), t_3c_der_ao(3), &
1881 t_3c_der_ao_ctr_1(3), t_3c_der_ri(3), t_3c_der_ri_ctr_1(3), t_3c_der_ri_ctr_2(3))
1882
1883 ! 1) Precompute the derivatives
1884 CALL precalc_derivatives(t_3c_der_ri_comp, t_3c_der_ao_comp, t_3c_der_ri_ind, t_3c_der_ao_ind, &
1885 t_2c_der_ri, t_2c_der_metric, t_3c_ri_ao_ao, &
1886 basis_set_ao, basis_set_ri, ri_data, qs_env)
1887
1888 DO ibasis = 1, SIZE(basis_set_ao)
1889 orb_basis => basis_set_ao(ibasis)%gto_basis_set
1890 ri_basis => basis_set_ri(ibasis)%gto_basis_set
1891 CALL init_interaction_radii_orb_basis(orb_basis, dft_control%qs_control%eps_pgf_orb)
1892 CALL init_interaction_radii_orb_basis(ri_basis, dft_control%qs_control%eps_pgf_orb)
1893 END DO
1894
1895 n_mem = SIZE(t_3c_der_ri_comp, 1)
1896 DO i_xyz = 1, 3
1897 CALL dbt_create(t_3c_ao_ri_ao, t_3c_der_ri(i_xyz))
1898 CALL dbt_create(t_3c_ao_ri_ao, t_3c_der_ao(i_xyz))
1899
1900 DO i_mem = 1, n_mem
1901 CALL decompress_tensor(t_3c_ri_ao_ao, t_3c_der_ri_ind(i_mem, i_xyz)%ind, &
1902 t_3c_der_ri_comp(i_mem, i_xyz), ri_data%filter_eps_storage)
1903 CALL dbt_copy(t_3c_ri_ao_ao, t_3c_der_ri(i_xyz), order=[2, 1, 3], move_data=.true., summation=.true.)
1904
1905 CALL decompress_tensor(t_3c_ri_ao_ao, t_3c_der_ao_ind(i_mem, i_xyz)%ind, &
1906 t_3c_der_ao_comp(i_mem, i_xyz), ri_data%filter_eps_storage)
1907 CALL dbt_copy(t_3c_ri_ao_ao, t_3c_der_ao(i_xyz), order=[2, 1, 3], move_data=.true., summation=.true.)
1908 END DO
1909 END DO
1910
1911 DO i_xyz = 1, 3
1912 DO i_mem = 1, n_mem
1913 CALL dealloc_containers(t_3c_der_ao_comp(i_mem, i_xyz), dummy_int)
1914 CALL dealloc_containers(t_3c_der_ri_comp(i_mem, i_xyz), dummy_int)
1915 END DO
1916 END DO
1917 DEALLOCATE (t_3c_der_ao_ind, t_3c_der_ri_ind)
1918
1919 ! Get the 3c integrals (desymmetrized)
1920 CALL dbt_create(t_3c_ao_ri_ao, t_3c_desymm)
1921 CALL dbt_copy(ri_data%t_3c_int_ctr_1(1, 1), t_3c_desymm)
1922 CALL dbt_copy(ri_data%t_3c_int_ctr_1(1, 1), t_3c_desymm, order=[3, 2, 1], &
1923 summation=.true., move_data=.true.)
1924
1925 CALL dbt_destroy(t_3c_ao_ri_ao)
1926 CALL dbt_destroy(t_3c_ri_ao_ao)
1927
1928 ! Some utilities
1929 spin_fac = 0.5_dp
1930 IF (nspins == 2) spin_fac = 1.0_dp
1931
1932 ALLOCATE (idx_to_at_ri(SIZE(ri_data%bsizes_RI_split)))
1933 CALL get_idx_to_atom(idx_to_at_ri, ri_data%bsizes_RI_split, ri_data%bsizes_RI)
1934
1935 ALLOCATE (idx_to_at_ao(SIZE(ri_data%bsizes_AO_split)))
1936 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
1937
1938 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
1939
1940 ! 2-center RI tensors
1941 CALL create_2c_tensor(t_2c_ri, dist1, dist2, ri_data%pgrid_2d, &
1942 ri_data%bsizes_RI_split, ri_data%bsizes_RI_split, name="(RI | RI)")
1943 DEALLOCATE (dist1, dist2)
1944
1945 CALL create_2c_tensor(t_2c_ri_pq, dist1, dist2, ri_data%pgrid_2d, &
1946 ri_data%bsizes_RI_fit, ri_data%bsizes_RI_fit, name="(RI | RI)")
1947 DEALLOCATE (dist1, dist2)
1948
1949 IF (.NOT. ri_data%same_op) THEN
1950 !precompute the (P|Q)*S^-1 product
1951 CALL dbt_create(t_2c_ri_pq, t_2c_ri_inv)
1952 CALL dbt_create(t_2c_ri_pq, t_2c_ri_met)
1953 CALL dbt_create(ri_data%t_2c_inv(1, 1), t_2c_tmp)
1954
1955 CALL dbt_contract(1.0_dp, ri_data%t_2c_inv(1, 1), ri_data%t_2c_pot(1, 1), &
1956 0.0_dp, t_2c_tmp, &
1957 contract_1=[2], notcontract_1=[1], &
1958 contract_2=[1], notcontract_2=[2], &
1959 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps, &
1960 unit_nr=unit_nr_dbcsr, flop=nflop)
1961 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
1962
1963 CALL dbt_copy(t_2c_tmp, t_2c_ri_inv, move_data=.true.)
1964 CALL dbt_destroy(t_2c_tmp)
1965 END IF
1966
1967 !3 loops in MO force evaluations. To be consistent with input MEMORY_CUT, need to take the cubic root
1968 !No need to cut memory further because SCF tensors alrady dense
1969 n_mem_input = floor((ri_data%n_mem_input - 0.1_dp)**(1._dp/3._dp)) + 1
1970 n_mem_input_ri = floor((ri_data%n_mem_input - 0.1_dp)/n_mem_input**2) + 1
1971
1972 !batches on RI_split and RI_fit blocks
1973 n_mem_ri = n_mem_input_ri
1974 CALL create_tensor_batches(ri_data%bsizes_RI_split, n_mem_ri, batch_start_ri, batch_end_ri, &
1975 batch_blk_start, batch_blk_end)
1976 ALLOCATE (batch_ranges_ri(n_mem_ri + 1))
1977 batch_ranges_ri(1:n_mem_ri) = batch_blk_start(1:n_mem_ri)
1978 batch_ranges_ri(n_mem_ri + 1) = batch_blk_end(n_mem_ri) + 1
1979 DEALLOCATE (batch_blk_start, batch_blk_end)
1980
1981 n_mem_ri_fit = n_mem_input_ri
1982 CALL create_tensor_batches(ri_data%bsizes_RI_fit, n_mem_ri_fit, batch_start_ri_fit, batch_end_ri_fit, &
1983 batch_blk_start, batch_blk_end)
1984 ALLOCATE (batch_ranges_ri_fit(n_mem_ri_fit + 1))
1985 batch_ranges_ri_fit(1:n_mem_ri_fit) = batch_blk_start(1:n_mem_ri_fit)
1986 batch_ranges_ri_fit(n_mem_ri_fit + 1) = batch_blk_end(n_mem_ri_fit) + 1
1987 DEALLOCATE (batch_blk_start, batch_blk_end)
1988
1989 DO ispin = 1, nspins
1990
1991 ! 2 )Prepare the batches for this spin
1992 CALL dbcsr_get_info(mo_coeff(ispin), nfullcols_total=n_mos)
1993 !note: optimized GPU block size for SCF is 64x1x64. Here we do 8x8x64
1994 CALL split_block_sizes([n_mos], bsizes_mo, max_size=floor(sqrt(ri_data%max_bsize_MO - 0.1)) + 1)
1995
1996 !batching on MO blocks
1997 n_mem = n_mem_input
1998 CALL create_tensor_batches(bsizes_mo, n_mem, batch_start, batch_end, &
1999 batch_blk_start, batch_blk_end)
2000 ALLOCATE (batch_ranges(n_mem + 1))
2001 batch_ranges(1:n_mem) = batch_blk_start(1:n_mem)
2002 batch_ranges(n_mem + 1) = batch_blk_end(n_mem) + 1
2003 DEALLOCATE (batch_blk_start, batch_blk_end)
2004
2005 ! Initialize the different tensors needed (Note: keep MO coeffs as (MO | AO) for less transpose)
2006 CALL create_2c_tensor(t_mo_coeff, dist1, dist2, ri_data%pgrid_2d, bsizes_mo, &
2007 ri_data%bsizes_AO_split, name="MO coeffs")
2008 DEALLOCATE (dist1, dist2)
2009 CALL dbt_create(mo_coeff(ispin), t_2c_tmp, name="MO coeffs")
2010 CALL dbt_copy_matrix_to_tensor(mo_coeff(ispin), t_2c_tmp)
2011 CALL dbt_copy(t_2c_tmp, t_mo_coeff, order=[2, 1], move_data=.true.)
2012 CALL dbt_destroy(t_2c_tmp)
2013
2014 CALL dbt_create(t_mo_coeff, t_mo_cpy)
2015 CALL dbt_copy(t_mo_coeff, t_mo_cpy)
2016 DO i_xyz = 1, 3
2017 CALL dbt_create(t_mo_coeff, t_2c_mo_ao_ctr(i_xyz))
2018 CALL dbt_create(t_mo_coeff, t_2c_mo_ao(i_xyz))
2019 END DO
2020
2021 CALL create_3c_tensor(t_3c_ao_ri_mo, dist1, dist2, dist3, pgrid_1, ri_data%bsizes_AO_split, &
2022 ri_data%bsizes_RI_split, bsizes_mo, [1, 2], [3], name="(AO RI| MO)")
2023 DEALLOCATE (dist1, dist2, dist3)
2024
2025 CALL dbt_create(t_3c_ao_ri_mo, t_3c_0)
2026 CALL dbt_destroy(t_3c_ao_ri_mo)
2027
2028 CALL create_3c_tensor(t_3c_mo_ri_ao, dist1, dist2, dist3, pgrid_1, bsizes_mo, ri_data%bsizes_RI_split, &
2029 ri_data%bsizes_AO_split, [1, 2], [3], name="(MO RI | AO)")
2030 DEALLOCATE (dist1, dist2, dist3)
2031 CALL dbt_create(t_3c_mo_ri_ao, t_3c_1)
2032
2033 DO i_xyz = 1, 3
2034 CALL dbt_create(t_3c_mo_ri_ao, t_3c_der_ri_ctr_1(i_xyz))
2035 CALL dbt_create(t_3c_mo_ri_ao, t_3c_der_ao_ctr_1(i_xyz))
2036 END DO
2037
2038 CALL create_3c_tensor(t_3c_mo_ri_mo, dist1, dist2, dist3, pgrid_1, bsizes_mo, &
2039 ri_data%bsizes_RI_split, bsizes_mo, [1, 2], [3], name="(MO RI | MO)")
2040 DEALLOCATE (dist1, dist2, dist3)
2041 CALL dbt_create(t_3c_mo_ri_mo, t_3c_work)
2042
2043 CALL create_3c_tensor(t_3c_ri_mo_mo, dist1, dist2, dist3, pgrid_2, ri_data%bsizes_RI_split, &
2044 bsizes_mo, bsizes_mo, [1], [2, 3], name="(RI| MO MO)")
2045 DEALLOCATE (dist1, dist2, dist3)
2046
2047 CALL dbt_create(t_3c_ri_mo_mo, t_3c_2)
2048 CALL dbt_create(t_3c_ri_mo_mo, t_3c_3)
2049 CALL dbt_create(t_3c_ri_mo_mo, t_3c_ri_ctr)
2050 DO i_xyz = 1, 3
2051 CALL dbt_create(t_3c_ri_mo_mo, t_3c_der_ri_ctr_2(i_xyz))
2052 END DO
2053
2054 !Very large RI_fit blocks => new pgrid to make sure distribution is ideal
2055 pdims(:) = 0
2056 CALL create_3c_tensor(t_3c_ri_mo_mo_fit, dist1, dist2, dist3, pgrid_2, ri_data%bsizes_RI_fit, &
2057 bsizes_mo, bsizes_mo, [1], [2, 3], name="(RI| MO MO)")
2058 DEALLOCATE (dist1, dist2, dist3)
2059
2060 CALL dbt_create(t_3c_ri_mo_mo_fit, t_3c_4)
2061 CALL dbt_create(t_3c_ri_mo_mo_fit, t_3c_5)
2062 CALL dbt_create(t_3c_ri_mo_mo_fit, t_3c_6)
2063
2064 CALL dbt_batched_contract_init(t_3c_desymm, batch_range_2=batch_ranges_ri)
2065 CALL dbt_batched_contract_init(t_3c_0, batch_range_2=batch_ranges_ri, batch_range_3=batch_ranges)
2066
2067 DO i_xyz = 1, 3
2068 CALL dbt_batched_contract_init(t_3c_der_ao(i_xyz), batch_range_2=batch_ranges_ri)
2069 CALL dbt_batched_contract_init(t_3c_der_ri(i_xyz), batch_range_2=batch_ranges_ri)
2070 END DO
2071
2072 CALL para_env%sync()
2073 t1 = m_walltime()
2074
2075 ! 2) Loop over batches
2076 DO i_mem = 1, n_mem
2077
2078 bounds_ctr_1d(1, 1) = batch_start(i_mem)
2079 bounds_ctr_1d(2, 1) = batch_end(i_mem)
2080
2081 bounds_ctr_2d(1, 1) = 1
2082 bounds_ctr_2d(2, 1) = sum(ri_data%bsizes_AO)
2083
2084 ! 3) Do the first AO to MO contraction here
2085 CALL timeset(routinen//"_AO2MO_1", handle)
2086 CALL dbt_batched_contract_init(t_mo_coeff)
2087 DO k_mem = 1, n_mem_ri
2088 bounds_ctr_2d(1, 2) = batch_start_ri(k_mem)
2089 bounds_ctr_2d(2, 2) = batch_end_ri(k_mem)
2090
2091 CALL dbt_contract(1.0_dp, t_mo_coeff, t_3c_desymm, &
2092 1.0_dp, t_3c_0, &
2093 contract_1=[2], notcontract_1=[1], &
2094 contract_2=[3], notcontract_2=[1, 2], &
2095 map_1=[3], map_2=[1, 2], filter_eps=ri_data%filter_eps, &
2096 bounds_2=bounds_ctr_1d, &
2097 bounds_3=bounds_ctr_2d, &
2098 unit_nr=unit_nr_dbcsr, flop=nflop)
2099 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2100 END DO
2101 CALL dbt_copy(t_3c_0, t_3c_1, order=[3, 2, 1], move_data=.true.)
2102
2103 DO i_xyz = 1, 3
2104 DO k_mem = 1, n_mem_ri
2105 bounds_ctr_2d(1, 2) = batch_start_ri(k_mem)
2106 bounds_ctr_2d(2, 2) = batch_end_ri(k_mem)
2107
2108 CALL dbt_contract(1.0_dp, t_mo_coeff, t_3c_der_ao(i_xyz), &
2109 1.0_dp, t_3c_0, &
2110 contract_1=[2], notcontract_1=[1], &
2111 contract_2=[3], notcontract_2=[1, 2], &
2112 map_1=[3], map_2=[1, 2], filter_eps=ri_data%filter_eps, &
2113 bounds_2=bounds_ctr_1d, &
2114 bounds_3=bounds_ctr_2d, &
2115 unit_nr=unit_nr_dbcsr, flop=nflop)
2116 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2117 END DO
2118 CALL dbt_copy(t_3c_0, t_3c_der_ao_ctr_1(i_xyz), order=[3, 2, 1], move_data=.true.)
2119
2120 DO k_mem = 1, n_mem_ri
2121 bounds_ctr_2d(1, 2) = batch_start_ri(k_mem)
2122 bounds_ctr_2d(2, 2) = batch_end_ri(k_mem)
2123
2124 CALL dbt_contract(1.0_dp, t_mo_coeff, t_3c_der_ri(i_xyz), &
2125 1.0_dp, t_3c_0, &
2126 contract_1=[2], notcontract_1=[1], &
2127 contract_2=[3], notcontract_2=[1, 2], &
2128 map_1=[3], map_2=[1, 2], filter_eps=ri_data%filter_eps, &
2129 bounds_2=bounds_ctr_1d, &
2130 bounds_3=bounds_ctr_2d, &
2131 unit_nr=unit_nr_dbcsr, flop=nflop)
2132 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2133 END DO
2134 CALL dbt_copy(t_3c_0, t_3c_der_ri_ctr_1(i_xyz), order=[3, 2, 1], move_data=.true.)
2135 END DO
2136 CALL dbt_batched_contract_finalize(t_mo_coeff)
2137 CALL timestop(handle)
2138
2139 CALL dbt_batched_contract_init(t_3c_1, batch_range_1=batch_ranges, batch_range_2=batch_ranges_ri)
2140 CALL dbt_batched_contract_init(t_3c_work, batch_range_1=batch_ranges, batch_range_2=batch_ranges_ri, &
2141 batch_range_3=batch_ranges)
2142 CALL dbt_batched_contract_init(t_3c_2, batch_range_2=batch_ranges, batch_range_3=batch_ranges)
2143 CALL dbt_batched_contract_init(t_3c_3, batch_range_1=batch_ranges_ri, &
2144 batch_range_2=batch_ranges, batch_range_3=batch_ranges)
2145
2146 CALL dbt_batched_contract_init(t_3c_4, batch_range_1=batch_ranges_ri_fit, &
2147 batch_range_2=batch_ranges, batch_range_3=batch_ranges)
2148 CALL dbt_batched_contract_init(t_3c_5, batch_range_2=batch_ranges, batch_range_3=batch_ranges)
2149
2150 DO i_xyz = 1, 3
2151 CALL dbt_batched_contract_init(t_3c_der_ri_ctr_1(i_xyz), batch_range_1=batch_ranges, &
2152 batch_range_2=batch_ranges_ri)
2153 CALL dbt_batched_contract_init(t_3c_der_ao_ctr_1(i_xyz), batch_range_1=batch_ranges, &
2154 batch_range_2=batch_ranges_ri)
2155
2156 END DO
2157
2158 IF (.NOT. ri_data%same_op) THEN
2159 CALL dbt_batched_contract_init(t_3c_6, batch_range_2=batch_ranges, batch_range_3=batch_ranges)
2160 END IF
2161
2162 DO j_mem = 1, n_mem
2163
2164 bounds_ctr_1d(1, 1) = batch_start(j_mem)
2165 bounds_ctr_1d(2, 1) = batch_end(j_mem)
2166
2167 bounds_ctr_2d(1, 1) = batch_start(i_mem)
2168 bounds_ctr_2d(2, 1) = batch_end(i_mem)
2169
2170 ! 3) Do the second AO to MO contraction here, followed by the S^-1 contraction
2171 CALL timeset(routinen//"_AO2MO_2", handle)
2172 CALL dbt_batched_contract_init(t_mo_coeff)
2173 DO k_mem = 1, n_mem_ri
2174 bounds_ctr_2d(1, 2) = batch_start_ri(k_mem)
2175 bounds_ctr_2d(2, 2) = batch_end_ri(k_mem)
2176
2177 CALL dbt_contract(1.0_dp, t_mo_coeff, t_3c_1, &
2178 1.0_dp, t_3c_work, &
2179 contract_1=[2], notcontract_1=[1], &
2180 contract_2=[3], notcontract_2=[1, 2], &
2181 map_1=[3], map_2=[1, 2], filter_eps=ri_data%filter_eps, &
2182 bounds_2=bounds_ctr_1d, &
2183 bounds_3=bounds_ctr_2d, &
2184 unit_nr=unit_nr_dbcsr, flop=nflop)
2185 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2186 END DO
2187 CALL dbt_batched_contract_finalize(t_mo_coeff)
2188 CALL timestop(handle)
2189
2190 bounds_ctr_2d(1, 1) = batch_start(i_mem)
2191 bounds_ctr_2d(2, 1) = batch_end(i_mem)
2192 bounds_ctr_2d(1, 2) = batch_start(j_mem)
2193 bounds_ctr_2d(2, 2) = batch_end(j_mem)
2194
2195 ! 4) Contract 3c MO integrals with S^-1 as well
2196 CALL timeset(routinen//"_2c_inv", handle)
2197 CALL dbt_copy(t_3c_work, t_3c_3, order=[2, 1, 3], move_data=.true.)
2198 DO k_mem = 1, n_mem_ri
2199 bounds_ctr_1d(1, 1) = batch_start_ri(k_mem)
2200 bounds_ctr_1d(2, 1) = batch_end_ri(k_mem)
2201
2202 CALL dbt_batched_contract_init(ri_data%t_2c_inv(1, 1))
2203 CALL dbt_contract(1.0_dp, ri_data%t_2c_inv(1, 1), t_3c_3, &
2204 1.0_dp, t_3c_2, &
2205 contract_1=[2], notcontract_1=[1], &
2206 contract_2=[1], notcontract_2=[2, 3], &
2207 map_1=[1], map_2=[2, 3], filter_eps=ri_data%filter_eps, &
2208 bounds_1=bounds_ctr_1d, &
2209 bounds_3=bounds_ctr_2d, &
2210 unit_nr=unit_nr_dbcsr, flop=nflop)
2211 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2212 CALL dbt_batched_contract_finalize(ri_data%t_2c_inv(1, 1))
2213 END DO
2214 CALL dbt_copy(t_3c_ri_mo_mo, t_3c_3)
2215 CALL timestop(handle)
2216
2217 !Only contract (ab|P') with MO coeffs since need AO rep for the force of (a'b|P)
2218 bounds_ctr_1d(1, 1) = batch_start(j_mem)
2219 bounds_ctr_1d(2, 1) = batch_end(j_mem)
2220
2221 bounds_ctr_2d(1, 1) = batch_start(i_mem)
2222 bounds_ctr_2d(2, 1) = batch_end(i_mem)
2223
2224 CALL timeset(routinen//"_AO2MO_2", handle)
2225 CALL dbt_batched_contract_init(t_mo_coeff)
2226 DO i_xyz = 1, 3
2227 DO k_mem = 1, n_mem_ri
2228 bounds_ctr_2d(1, 2) = batch_start_ri(k_mem)
2229 bounds_ctr_2d(2, 2) = batch_end_ri(k_mem)
2230
2231 CALL dbt_contract(1.0_dp, t_mo_coeff, t_3c_der_ri_ctr_1(i_xyz), &
2232 1.0_dp, t_3c_work, &
2233 contract_1=[2], notcontract_1=[1], &
2234 contract_2=[3], notcontract_2=[1, 2], &
2235 map_1=[3], map_2=[1, 2], filter_eps=ri_data%filter_eps, &
2236 bounds_2=bounds_ctr_1d, &
2237 bounds_3=bounds_ctr_2d, &
2238 unit_nr=unit_nr_dbcsr, flop=nflop)
2239 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2240 END DO
2241 CALL dbt_copy(t_3c_work, t_3c_der_ri_ctr_2(i_xyz), order=[2, 1, 3], move_data=.true.)
2242 END DO
2243 CALL dbt_batched_contract_finalize(t_mo_coeff)
2244 CALL timestop(handle)
2245
2246 bounds_ctr_2d(1, 1) = batch_start(i_mem)
2247 bounds_ctr_2d(2, 1) = batch_end(i_mem)
2248 bounds_ctr_2d(1, 2) = batch_start(j_mem)
2249 bounds_ctr_2d(2, 2) = batch_end(j_mem)
2250
2251 ! 5) Force due to d/dx (P|Q)
2252 CALL timeset(routinen//"_PQ_der", handle)
2253 CALL dbt_copy(t_3c_2, t_3c_4, move_data=.true.)
2254 CALL dbt_copy(t_3c_4, t_3c_5)
2255 DO k_mem = 1, n_mem_ri_fit
2256 bounds_ctr_1d(1, 1) = batch_start_ri_fit(k_mem)
2257 bounds_ctr_1d(2, 1) = batch_end_ri_fit(k_mem)
2258
2259 CALL dbt_batched_contract_init(t_2c_ri_pq)
2260 CALL dbt_contract(1.0_dp, t_3c_4, t_3c_5, &
2261 1.0_dp, t_2c_ri_pq, &
2262 contract_1=[2, 3], notcontract_1=[1], &
2263 contract_2=[2, 3], notcontract_2=[1], &
2264 bounds_1=bounds_ctr_2d, &
2265 bounds_2=bounds_ctr_1d, &
2266 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps, &
2267 unit_nr=unit_nr_dbcsr, flop=nflop)
2268 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2269 CALL dbt_batched_contract_finalize(t_2c_ri_pq)
2270 END DO
2271 CALL timestop(handle)
2272
2273 ! 6) If metric, do the additional contraction with S_pq^-1 (Q|R) (not on the derivatives)
2274 IF (.NOT. ri_data%same_op) THEN
2275 CALL timeset(routinen//"_metric", handle)
2276 DO k_mem = 1, n_mem_ri_fit
2277 bounds_ctr_1d(1, 1) = batch_start_ri_fit(k_mem)
2278 bounds_ctr_1d(2, 1) = batch_end_ri_fit(k_mem)
2279
2280 CALL dbt_batched_contract_init(t_2c_ri_inv)
2281 CALL dbt_contract(1.0_dp, t_2c_ri_inv, t_3c_4, &
2282 1.0_dp, t_3c_6, &
2283 contract_1=[2], notcontract_1=[1], &
2284 contract_2=[1], notcontract_2=[2, 3], &
2285 bounds_1=bounds_ctr_1d, &
2286 bounds_3=bounds_ctr_2d, &
2287 map_1=[1], map_2=[2, 3], filter_eps=ri_data%filter_eps, &
2288 unit_nr=unit_nr_dbcsr, flop=nflop)
2289 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2290 CALL dbt_batched_contract_finalize(t_2c_ri_inv)
2291 END DO
2292 CALL dbt_copy(t_3c_6, t_3c_4, move_data=.true.)
2293
2294 ! 8) and get the force due to d/dx S^-1
2295 DO k_mem = 1, n_mem_ri_fit
2296 bounds_ctr_1d(1, 1) = batch_start_ri_fit(k_mem)
2297 bounds_ctr_1d(2, 1) = batch_end_ri_fit(k_mem)
2298
2299 CALL dbt_batched_contract_init(t_2c_ri_met)
2300 CALL dbt_contract(1.0_dp, t_3c_4, t_3c_5, &
2301 1.0_dp, t_2c_ri_met, &
2302 contract_1=[2, 3], notcontract_1=[1], &
2303 contract_2=[2, 3], notcontract_2=[1], &
2304 bounds_1=bounds_ctr_2d, &
2305 bounds_2=bounds_ctr_1d, &
2306 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps, &
2307 unit_nr=unit_nr_dbcsr, flop=nflop)
2308 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2309 CALL dbt_batched_contract_finalize(t_2c_ri_met)
2310 END DO
2311 CALL timestop(handle)
2312 END IF
2313 CALL dbt_copy(t_3c_ri_mo_mo_fit, t_3c_5)
2314
2315 ! 7) Do the force contribution due to 3c integrals (a'b|P) and (ab|P')
2316
2317 ! (ab|P')
2318 CALL timeset(routinen//"_3c_RI", handle)
2319 pref = -0.5_dp*2.0_dp*hf_fraction*spin_fac
2320 CALL dbt_copy(t_3c_4, t_3c_ri_ctr, move_data=.true.)
2321 CALL get_force_from_3c_trace(force, t_3c_ri_ctr, t_3c_der_ri_ctr_2, atom_of_kind, kind_of, &
2322 idx_to_at_ri, pref)
2323 CALL timestop(handle)
2324
2325 ! (a'b|P) Note that derivative remains in AO rep until the actual force evaluation,
2326 ! which also prevents doing a direct 3-center trace
2327 bounds_ctr_2d(1, 1) = batch_start(i_mem)
2328 bounds_ctr_2d(2, 1) = batch_end(i_mem)
2329
2330 bounds_ctr_1d(1, 1) = batch_start(j_mem)
2331 bounds_ctr_1d(2, 1) = batch_end(j_mem)
2332
2333 CALL timeset(routinen//"_3c_AO", handle)
2334 CALL dbt_copy(t_3c_ri_ctr, t_3c_work, order=[2, 1, 3], move_data=.true.)
2335 DO i_xyz = 1, 3
2336
2337 CALL dbt_batched_contract_init(t_2c_mo_ao_ctr(i_xyz))
2338 DO k_mem = 1, n_mem_ri
2339 bounds_ctr_2d(1, 2) = batch_start_ri(k_mem)
2340 bounds_ctr_2d(2, 2) = batch_end_ri(k_mem)
2341
2342 CALL dbt_contract(1.0_dp, t_3c_work, t_3c_der_ao_ctr_1(i_xyz), &
2343 1.0_dp, t_2c_mo_ao_ctr(i_xyz), &
2344 contract_1=[1, 2], notcontract_1=[3], &
2345 contract_2=[1, 2], notcontract_2=[3], &
2346 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps, &
2347 bounds_1=bounds_ctr_2d, &
2348 bounds_2=bounds_ctr_1d, &
2349 unit_nr=unit_nr_dbcsr, flop=nflop)
2350 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2351 END DO
2352 CALL dbt_batched_contract_finalize(t_2c_mo_ao_ctr(i_xyz))
2353 END DO
2354 CALL timestop(handle)
2355
2356 END DO !j_mem
2357 CALL dbt_batched_contract_finalize(t_3c_1)
2358 CALL dbt_batched_contract_finalize(t_3c_work)
2359 CALL dbt_batched_contract_finalize(t_3c_2)
2360 CALL dbt_batched_contract_finalize(t_3c_3)
2361 CALL dbt_batched_contract_finalize(t_3c_4)
2362 CALL dbt_batched_contract_finalize(t_3c_5)
2363
2364 DO i_xyz = 1, 3
2365 CALL dbt_batched_contract_finalize(t_3c_der_ri_ctr_1(i_xyz))
2366 CALL dbt_batched_contract_finalize(t_3c_der_ao_ctr_1(i_xyz))
2367 END DO
2368
2369 IF (.NOT. ri_data%same_op) THEN
2370 CALL dbt_batched_contract_finalize(t_3c_6)
2371 END IF
2372
2373 END DO !i_mem
2374 CALL dbt_batched_contract_finalize(t_3c_desymm)
2375 CALL dbt_batched_contract_finalize(t_3c_0)
2376
2377 DO i_xyz = 1, 3
2378 CALL dbt_batched_contract_finalize(t_3c_der_ao(i_xyz))
2379 CALL dbt_batched_contract_finalize(t_3c_der_ri(i_xyz))
2380 END DO
2381
2382 !Force contribution due to 3-center AO derivatives (a'b|P)
2383 pref = -0.5_dp*4.0_dp*hf_fraction*spin_fac
2384 DO i_xyz = 1, 3
2385 CALL dbt_copy(t_2c_mo_ao_ctr(i_xyz), t_2c_mo_ao(i_xyz), move_data=.true.) !ensures matching distributions
2386 CALL get_mo_ao_force(force, t_mo_cpy, t_2c_mo_ao(i_xyz), atom_of_kind, kind_of, idx_to_at_ao, pref, i_xyz)
2387 CALL dbt_clear(t_2c_mo_ao(i_xyz))
2388 END DO
2389
2390 !Force contribution of d/dx (P|Q)
2391 pref = 0.5_dp*hf_fraction*spin_fac
2392 IF (.NOT. ri_data%same_op) pref = -pref
2393
2394 !Making sure dists of the t_2c_RI tensors match
2395 CALL dbt_copy(t_2c_ri_pq, t_2c_ri, move_data=.true.)
2396 CALL get_2c_der_force(force, t_2c_ri, t_2c_der_ri, atom_of_kind, &
2397 kind_of, idx_to_at_ri, pref)
2398 CALL dbt_clear(t_2c_ri)
2399
2400 !Force contribution due to the inverse metric
2401 IF (.NOT. ri_data%same_op) THEN
2402 pref = 0.5_dp*2.0_dp*hf_fraction*spin_fac
2403
2404 CALL dbt_copy(t_2c_ri_met, t_2c_ri, move_data=.true.)
2405 CALL get_2c_der_force(force, t_2c_ri, t_2c_der_metric, atom_of_kind, &
2406 kind_of, idx_to_at_ri, pref)
2407 CALL dbt_clear(t_2c_ri)
2408 END IF
2409
2410 CALL dbt_destroy(t_3c_0)
2411 CALL dbt_destroy(t_3c_1)
2412 CALL dbt_destroy(t_3c_2)
2413 CALL dbt_destroy(t_3c_3)
2414 CALL dbt_destroy(t_3c_4)
2415 CALL dbt_destroy(t_3c_5)
2416 CALL dbt_destroy(t_3c_6)
2417 CALL dbt_destroy(t_3c_work)
2418 CALL dbt_destroy(t_3c_ri_ctr)
2419 CALL dbt_destroy(t_3c_mo_ri_ao)
2420 CALL dbt_destroy(t_3c_mo_ri_mo)
2421 CALL dbt_destroy(t_3c_ri_mo_mo)
2422 CALL dbt_destroy(t_3c_ri_mo_mo_fit)
2423 CALL dbt_destroy(t_mo_coeff)
2424 CALL dbt_destroy(t_mo_cpy)
2425 DO i_xyz = 1, 3
2426 CALL dbt_destroy(t_2c_mo_ao(i_xyz))
2427 CALL dbt_destroy(t_2c_mo_ao_ctr(i_xyz))
2428 CALL dbt_destroy(t_3c_der_ri_ctr_1(i_xyz))
2429 CALL dbt_destroy(t_3c_der_ao_ctr_1(i_xyz))
2430 CALL dbt_destroy(t_3c_der_ri_ctr_2(i_xyz))
2431 END DO
2432 DEALLOCATE (batch_ranges, batch_start, batch_end)
2433 END DO !ispin
2434
2435 ! Clean-up
2436 CALL dbt_pgrid_destroy(pgrid_1)
2437 CALL dbt_pgrid_destroy(pgrid_2)
2438 CALL dbt_destroy(t_3c_desymm)
2439 CALL dbt_destroy(t_2c_ri)
2440 CALL dbt_destroy(t_2c_ri_pq)
2441 IF (.NOT. ri_data%same_op) THEN
2442 CALL dbt_destroy(t_2c_ri_met)
2443 CALL dbt_destroy(t_2c_ri_inv)
2444 END IF
2445 DO i_xyz = 1, 3
2446 CALL dbt_destroy(t_3c_der_ao(i_xyz))
2447 CALL dbt_destroy(t_3c_der_ri(i_xyz))
2448 CALL dbt_destroy(t_2c_der_ri(i_xyz))
2449 IF (.NOT. ri_data%same_op) CALL dbt_destroy(t_2c_der_metric(i_xyz))
2450 END DO
2451 CALL dbt_copy(ri_data%t_3c_int_ctr_2(1, 1), ri_data%t_3c_int_ctr_1(1, 1))
2452
2453 CALL para_env%sync()
2454 t2 = m_walltime()
2455
2456 ri_data%dbcsr_time = ri_data%dbcsr_time + t2 - t1
2457
2458 END SUBROUTINE hfx_ri_forces_mo
2459
2460! **************************************************************************************************
2461!> \brief New sparser implementation
2462!> \param qs_env ...
2463!> \param ri_data ...
2464!> \param nspins ...
2465!> \param hf_fraction ...
2466!> \param rho_ao ...
2467!> \param rho_ao_resp ...
2468!> \param use_virial ...
2469!> \param resp_only ...
2470!> \param rescale_factor ...
2471! **************************************************************************************************
2472 SUBROUTINE hfx_ri_forces_pmat(qs_env, ri_data, nspins, hf_fraction, rho_ao, rho_ao_resp, &
2473 use_virial, resp_only, rescale_factor)
2474
2475 TYPE(qs_environment_type), POINTER :: qs_env
2476 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
2477 INTEGER, INTENT(IN) :: nspins
2478 REAL(dp), INTENT(IN) :: hf_fraction
2479 TYPE(dbcsr_p_type), DIMENSION(:, :) :: rho_ao
2480 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL :: rho_ao_resp
2481 LOGICAL, INTENT(IN), OPTIONAL :: use_virial, resp_only
2482 REAL(dp), INTENT(IN), OPTIONAL :: rescale_factor
2483
2484 CHARACTER(LEN=*), PARAMETER :: routinen = 'hfx_ri_forces_Pmat'
2485
2486 INTEGER :: dummy_int, handle, i_mem, i_spin, i_xyz, &
2487 ibasis, j_mem, j_xyz, k_mem, k_xyz, &
2488 n_mem, n_mem_ri, natom, nkind, &
2489 unit_nr_dbcsr
2490 INTEGER(int_8) :: nflop
2491 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, batch_end, batch_end_ri, batch_ranges, &
2492 batch_ranges_ri, batch_start, batch_start_ri, dist1, dist2, dist3, idx_to_at_ao, &
2493 idx_to_at_ri, kind_of
2494 INTEGER, DIMENSION(2, 1) :: ibounds, jbounds, kbounds
2495 INTEGER, DIMENSION(2, 2) :: ijbounds
2496 INTEGER, DIMENSION(2, 3) :: bounds_cpy
2497 INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
2498 LOGICAL :: do_resp, resp_only_prv, use_virial_prv
2499 REAL(dp) :: pref, spin_fac, t1, t2
2500 REAL(dp), DIMENSION(3, 3) :: work_virial
2501 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2502 TYPE(block_ind_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_der_ao_ind, t_3c_der_ri_ind
2503 TYPE(cell_type), POINTER :: cell
2504 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
2505 TYPE(dbcsr_type) :: dbcsr_tmp, virial_trace
2506 TYPE(dbt_type) :: rho_ao_1, rho_ao_2, t_2c_ri, t_2c_ri_tmp, t_2c_tmp, t_2c_virial, t_3c_1, &
2507 t_3c_2, t_3c_3, t_3c_4, t_3c_5, t_3c_ao_ri_ao, t_3c_help_1, t_3c_help_2, t_3c_int, &
2508 t_3c_int_2, t_3c_ri_ao_ao, t_3c_sparse, t_3c_virial, t_r, t_svs
2509 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_2c_der_metric, t_2c_der_ri, &
2510 t_3c_der_ao, t_3c_der_ri
2511 TYPE(dft_control_type), POINTER :: dft_control
2512 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
2513 DIMENSION(:), TARGET :: basis_set_ao, basis_set_ri
2514 TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
2515 TYPE(hfx_compression_type), ALLOCATABLE, &
2516 DIMENSION(:, :) :: t_3c_der_ao_comp, t_3c_der_ri_comp
2517 TYPE(mp_para_env_type), POINTER :: para_env
2518 TYPE(neighbor_list_3c_type) :: nl_3c
2519 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2520 POINTER :: nl_2c_met, nl_2c_pot
2521 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2522 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2523 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2524 TYPE(virial_type), POINTER :: virial
2525
2526 !The idea is the following: we need to compute the gradients
2527 ! d/dx [P_ab P_cd (acP) S^-1_PQ (Q|R) S^-1_RS (Sbd)]
2528 ! Which we do in a few steps:
2529 ! 1) Contract the density matrices with the 3c integrals: M_acS = P_ab P_cd (Sbd)
2530 ! 2) Calculate the 3c contributions: d/dx (acP) [S^-1_PQ (Q|R) S^-1_RS M_acS]
2531 ! For maximum perf, we first multiply all 2c matrices together, than contract with retain_sparsity
2532 ! 3) Contract the 3c integrals and the M tensor together in order to only work with 2c quantities:
2533 ! R_PS = (acP) M_acS
2534 ! 4) From there, we can easily calculate the 2c contributions to the force:
2535 ! Potential: [S^-1*R*S^-1]_QR d/dx (Q|R)
2536 ! Metric: [S^-1*R*S^-1*(Q|R)*S^-1]_UV d/dx S_UV
2537
2538 NULLIFY (particle_set, virial, cell, force, atomic_kind_set, nl_2c_pot, nl_2c_met)
2539 NULLIFY (orb_basis, ri_basis, qs_kind_set, particle_set, dft_control, dbcsr_dist)
2540
2541 use_virial_prv = .false.
2542 IF (PRESENT(use_virial)) use_virial_prv = use_virial
2543
2544 do_resp = .false.
2545 IF (PRESENT(rho_ao_resp)) THEN
2546 IF (ASSOCIATED(rho_ao_resp(1)%matrix)) do_resp = .true.
2547 END IF
2548
2549 resp_only_prv = .false.
2550 IF (PRESENT(resp_only)) resp_only_prv = resp_only
2551
2552 unit_nr_dbcsr = ri_data%unit_nr_dbcsr
2553
2554 CALL get_qs_env(qs_env, natom=natom, particle_set=particle_set, nkind=nkind, &
2555 atomic_kind_set=atomic_kind_set, virial=virial, &
2556 cell=cell, force=force, para_env=para_env, dft_control=dft_control, &
2557 qs_kind_set=qs_kind_set, dbcsr_dist=dbcsr_dist)
2558
2559 CALL create_3c_tensor(t_3c_ao_ri_ao, dist1, dist2, dist3, ri_data%pgrid_1, &
2560 ri_data%bsizes_AO_split, ri_data%bsizes_RI_split, ri_data%bsizes_AO_split, &
2561 [1, 2], [3], name="(AO RI | AO)")
2562 DEALLOCATE (dist1, dist2, dist3)
2563
2564 CALL create_3c_tensor(t_3c_ri_ao_ao, dist1, dist2, dist3, ri_data%pgrid_2, &
2565 ri_data%bsizes_RI_split, ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
2566 [1], [2, 3], name="(RI | AO AO)")
2567 DEALLOCATE (dist1, dist2, dist3)
2568
2569 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
2570 CALL basis_set_list_setup(basis_set_ri, ri_data%ri_basis_type, qs_kind_set)
2571 CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_ri)
2572 CALL basis_set_list_setup(basis_set_ao, ri_data%orb_basis_type, qs_kind_set)
2573 CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_ao)
2574
2575 DO ibasis = 1, SIZE(basis_set_ao)
2576 orb_basis => basis_set_ao(ibasis)%gto_basis_set
2577 CALL init_interaction_radii_orb_basis(orb_basis, ri_data%eps_pgf_orb)
2578 ri_basis => basis_set_ri(ibasis)%gto_basis_set
2579 CALL init_interaction_radii_orb_basis(ri_basis, ri_data%eps_pgf_orb)
2580 END DO
2581
2582 ! Precompute the derivatives
2583 ALLOCATE (t_2c_der_metric(3), t_2c_der_ri(3), t_3c_der_ao(3), t_3c_der_ri(3))
2584 IF (use_virial) THEN
2585 CALL precalc_derivatives(t_3c_der_ri_comp, t_3c_der_ao_comp, t_3c_der_ri_ind, t_3c_der_ao_ind, &
2586 t_2c_der_ri, t_2c_der_metric, t_3c_ri_ao_ao, &
2587 basis_set_ao, basis_set_ri, ri_data, qs_env, &
2588 nl_2c_pot=nl_2c_pot, nl_2c_met=nl_2c_met, &
2589 nl_3c_out=nl_3c, t_3c_virial=t_3c_virial)
2590
2591 ALLOCATE (col_bsize(natom), row_bsize(natom))
2592 col_bsize(:) = ri_data%bsizes_RI
2593 row_bsize(:) = ri_data%bsizes_RI
2594 CALL dbcsr_create(virial_trace, "virial_trace", dbcsr_dist, dbcsr_type_no_symmetry, row_bsize, col_bsize)
2595 CALL dbt_create(virial_trace, t_2c_virial)
2596 DEALLOCATE (col_bsize, row_bsize)
2597 ELSE
2598 CALL precalc_derivatives(t_3c_der_ri_comp, t_3c_der_ao_comp, t_3c_der_ri_ind, t_3c_der_ao_ind, &
2599 t_2c_der_ri, t_2c_der_metric, t_3c_ri_ao_ao, &
2600 basis_set_ao, basis_set_ri, ri_data, qs_env)
2601 END IF
2602
2603 ! Keep track of derivative sparsity to be able to use retain_sparsity in contraction
2604 CALL dbt_create(t_3c_ri_ao_ao, t_3c_sparse)
2605 DO i_xyz = 1, 3
2606 DO i_mem = 1, SIZE(t_3c_der_ri_comp, 1)
2607 CALL decompress_tensor(t_3c_ri_ao_ao, t_3c_der_ri_ind(i_mem, i_xyz)%ind, &
2608 t_3c_der_ri_comp(i_mem, i_xyz), ri_data%filter_eps_storage)
2609 CALL dbt_copy(t_3c_ri_ao_ao, t_3c_sparse, summation=.true., move_data=.true.)
2610
2611 CALL decompress_tensor(t_3c_ri_ao_ao, t_3c_der_ao_ind(i_mem, i_xyz)%ind, &
2612 t_3c_der_ao_comp(i_mem, i_xyz), ri_data%filter_eps_storage)
2613 CALL dbt_copy(t_3c_ri_ao_ao, t_3c_sparse, summation=.true.)
2614 CALL dbt_copy(t_3c_ri_ao_ao, t_3c_sparse, order=[1, 3, 2], summation=.true., move_data=.true.)
2615 END DO
2616 END DO
2617
2618 DO i_xyz = 1, 3
2619 CALL dbt_create(t_3c_ri_ao_ao, t_3c_der_ri(i_xyz))
2620 CALL dbt_create(t_3c_ri_ao_ao, t_3c_der_ao(i_xyz))
2621 END DO
2622
2623 ! Some utilities
2624 spin_fac = 0.5_dp
2625 IF (nspins == 2) spin_fac = 1.0_dp
2626 IF (PRESENT(rescale_factor)) spin_fac = spin_fac*rescale_factor
2627
2628 ALLOCATE (idx_to_at_ri(SIZE(ri_data%bsizes_RI_split)))
2629 CALL get_idx_to_atom(idx_to_at_ri, ri_data%bsizes_RI_split, ri_data%bsizes_RI)
2630
2631 ALLOCATE (idx_to_at_ao(SIZE(ri_data%bsizes_AO_split)))
2632 CALL get_idx_to_atom(idx_to_at_ao, ri_data%bsizes_AO_split, ri_data%bsizes_AO)
2633
2634 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
2635
2636 ! Go over batches of the 2 AO indices to save memory
2637 n_mem = ri_data%n_mem
2638 ALLOCATE (batch_start(n_mem), batch_end(n_mem))
2639 batch_start(:) = ri_data%starts_array_mem(:)
2640 batch_end(:) = ri_data%ends_array_mem(:)
2641
2642 ALLOCATE (batch_ranges(n_mem + 1))
2643 batch_ranges(:n_mem) = ri_data%starts_array_mem_block(:)
2644 batch_ranges(n_mem + 1) = ri_data%ends_array_mem_block(n_mem) + 1
2645
2646 n_mem_ri = ri_data%n_mem_RI
2647 ALLOCATE (batch_start_ri(n_mem_ri), batch_end_ri(n_mem_ri))
2648 batch_start_ri(:) = ri_data%starts_array_RI_mem(:)
2649 batch_end_ri(:) = ri_data%ends_array_RI_mem(:)
2650
2651 ALLOCATE (batch_ranges_ri(n_mem_ri + 1))
2652 batch_ranges_ri(:n_mem_ri) = ri_data%starts_array_RI_mem_block(:)
2653 batch_ranges_ri(n_mem_ri + 1) = ri_data%ends_array_RI_mem_block(n_mem_ri) + 1
2654
2655 ! Pre-create all the needed tensors
2656 CALL create_2c_tensor(rho_ao_1, dist1, dist2, ri_data%pgrid_2d, &
2657 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
2658 name="(AO | AO)")
2659 DEALLOCATE (dist1, dist2)
2660 CALL dbt_create(rho_ao_1, rho_ao_2)
2661
2662 CALL create_2c_tensor(t_2c_ri, dist1, dist2, ri_data%pgrid_2d, &
2663 ri_data%bsizes_RI_split, ri_data%bsizes_RI_split, name="(RI | RI)")
2664 DEALLOCATE (dist1, dist2)
2665 CALL dbt_create(t_2c_ri, t_svs)
2666 CALL dbt_create(t_2c_ri, t_r)
2667 CALL dbt_create(t_2c_ri, t_2c_ri_tmp)
2668
2669 CALL dbt_create(t_3c_ao_ri_ao, t_3c_1)
2670 CALL dbt_create(t_3c_ao_ri_ao, t_3c_2)
2671 CALL dbt_create(t_3c_ri_ao_ao, t_3c_3)
2672 CALL dbt_create(t_3c_ri_ao_ao, t_3c_4)
2673 CALL dbt_create(t_3c_ri_ao_ao, t_3c_5)
2674 CALL dbt_create(t_3c_ri_ao_ao, t_3c_help_1)
2675 CALL dbt_create(t_3c_ri_ao_ao, t_3c_help_2)
2676
2677 CALL dbt_create(t_3c_ao_ri_ao, t_3c_int)
2678 CALL dbt_copy(ri_data%t_3c_int_ctr_2(1, 1), t_3c_int)
2679
2680 CALL dbt_create(t_3c_ri_ao_ao, t_3c_int_2)
2681
2682 CALL para_env%sync()
2683 t1 = m_walltime()
2684
2685 !Pre-calculate the necessary 2-center quantities
2686 IF (.NOT. ri_data%same_op) THEN
2687 !S^-1 * V * S^-1
2688 CALL dbt_contract(1.0_dp, ri_data%t_2c_inv(1, 1), ri_data%t_2c_pot(1, 1), 0.0_dp, t_2c_ri, &
2689 contract_1=[2], notcontract_1=[1], &
2690 contract_2=[1], notcontract_2=[2], &
2691 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps, &
2692 unit_nr=unit_nr_dbcsr, flop=nflop)
2693 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2694
2695 CALL dbt_contract(1.0_dp, t_2c_ri, ri_data%t_2c_inv(1, 1), 0.0_dp, t_svs, &
2696 contract_1=[2], notcontract_1=[1], &
2697 contract_2=[1], notcontract_2=[2], &
2698 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps, &
2699 unit_nr=unit_nr_dbcsr, flop=nflop)
2700 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2701 ELSE
2702 ! Simply V^-1
2703 CALL dbt_copy(ri_data%t_2c_inv(1, 1), t_svs)
2704 END IF
2705
2706 CALL dbt_batched_contract_init(t_3c_int, batch_range_1=batch_ranges, batch_range_3=batch_ranges)
2707 CALL dbt_batched_contract_init(t_3c_int_2, batch_range_1=batch_ranges_ri, &
2708 batch_range_2=batch_ranges, batch_range_3=batch_ranges)
2709 CALL dbt_batched_contract_init(t_3c_1, batch_range_1=batch_ranges, batch_range_3=batch_ranges)
2710 CALL dbt_batched_contract_init(t_3c_2, batch_range_1=batch_ranges, batch_range_3=batch_ranges)
2711 CALL dbt_batched_contract_init(t_3c_3, batch_range_1=batch_ranges_ri, &
2712 batch_range_2=batch_ranges, batch_range_3=batch_ranges)
2713 CALL dbt_batched_contract_init(t_3c_4, batch_range_1=batch_ranges_ri, &
2714 batch_range_2=batch_ranges, batch_range_3=batch_ranges)
2715 CALL dbt_batched_contract_init(t_3c_5, batch_range_1=batch_ranges_ri, &
2716 batch_range_2=batch_ranges, batch_range_3=batch_ranges)
2717 CALL dbt_batched_contract_init(t_3c_sparse, batch_range_1=batch_ranges_ri, &
2718 batch_range_2=batch_ranges, batch_range_3=batch_ranges)
2719
2720 DO i_spin = 1, nspins
2721
2722 !Prepare Pmat in tensor format
2723 CALL dbt_create(rho_ao(i_spin, 1)%matrix, t_2c_tmp)
2724 CALL dbt_copy_matrix_to_tensor(rho_ao(i_spin, 1)%matrix, t_2c_tmp)
2725 CALL dbt_copy(t_2c_tmp, rho_ao_1, move_data=.true.)
2726 CALL dbt_destroy(t_2c_tmp)
2727
2728 IF (.NOT. do_resp) THEN
2729 CALL dbt_copy(rho_ao_1, rho_ao_2)
2730 ELSE IF (do_resp .AND. resp_only_prv) THEN
2731
2732 CALL dbt_create(rho_ao_resp(i_spin)%matrix, t_2c_tmp)
2733 CALL dbt_copy_matrix_to_tensor(rho_ao_resp(i_spin)%matrix, t_2c_tmp)
2734 CALL dbt_copy(t_2c_tmp, rho_ao_2)
2735 !symmetry allows to take 2*P_resp rasther than explicitely take all cross products
2736 CALL dbt_copy(t_2c_tmp, rho_ao_2, summation=.true., move_data=.true.)
2737 CALL dbt_destroy(t_2c_tmp)
2738 ELSE
2739
2740 !if not resp_only, need P-P_resp and P+P_resp
2741 CALL dbt_copy(rho_ao_1, rho_ao_2)
2742 CALL dbcsr_create(dbcsr_tmp, template=rho_ao_resp(i_spin)%matrix)
2743 CALL dbcsr_add(dbcsr_tmp, rho_ao_resp(i_spin)%matrix, 0.0_dp, -1.0_dp)
2744 CALL dbt_create(dbcsr_tmp, t_2c_tmp)
2745 CALL dbt_copy_matrix_to_tensor(dbcsr_tmp, t_2c_tmp)
2746 CALL dbt_copy(t_2c_tmp, rho_ao_1, summation=.true., move_data=.true.)
2747 CALL dbcsr_release(dbcsr_tmp)
2748
2749 CALL dbt_copy_matrix_to_tensor(rho_ao_resp(i_spin)%matrix, t_2c_tmp)
2750 CALL dbt_copy(t_2c_tmp, rho_ao_2, summation=.true., move_data=.true.)
2751 CALL dbt_destroy(t_2c_tmp)
2752
2753 END IF
2754 work_virial = 0.0_dp
2755
2756 CALL timeset(routinen//"_3c", handle)
2757 !Start looping of the batches
2758 DO i_mem = 1, n_mem
2759 ibounds(:, 1) = [batch_start(i_mem), batch_end(i_mem)]
2760
2761 CALL dbt_batched_contract_init(rho_ao_1)
2762 CALL dbt_contract(1.0_dp, rho_ao_1, t_3c_int, 0.0_dp, t_3c_1, &
2763 contract_1=[1], notcontract_1=[2], &
2764 contract_2=[3], notcontract_2=[1, 2], &
2765 map_1=[3], map_2=[1, 2], filter_eps=ri_data%filter_eps, &
2766 bounds_2=ibounds, unit_nr=unit_nr_dbcsr, flop=nflop)
2767 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2768 CALL dbt_batched_contract_finalize(rho_ao_1)
2769
2770 CALL dbt_copy(t_3c_1, t_3c_2, order=[3, 2, 1], move_data=.true.)
2771
2772 DO j_mem = 1, n_mem
2773 jbounds(:, 1) = [batch_start(j_mem), batch_end(j_mem)]
2774
2775 CALL dbt_batched_contract_init(rho_ao_2)
2776 CALL dbt_contract(1.0_dp, rho_ao_2, t_3c_2, 0.0_dp, t_3c_1, &
2777 contract_1=[1], notcontract_1=[2], &
2778 contract_2=[3], notcontract_2=[1, 2], &
2779 map_1=[3], map_2=[1, 2], filter_eps=ri_data%filter_eps, &
2780 bounds_2=jbounds, unit_nr=unit_nr_dbcsr, flop=nflop)
2781 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2782 CALL dbt_batched_contract_finalize(rho_ao_2)
2783
2784 bounds_cpy(:, 1) = [batch_start(i_mem), batch_end(i_mem)]
2785 bounds_cpy(:, 2) = [1, sum(ri_data%bsizes_RI)]
2786 bounds_cpy(:, 3) = [batch_start(j_mem), batch_end(j_mem)]
2787 CALL dbt_copy(t_3c_int, t_3c_int_2, order=[2, 1, 3], bounds=bounds_cpy)
2788 CALL dbt_copy(t_3c_1, t_3c_3, order=[2, 1, 3], move_data=.true.)
2789
2790 DO k_mem = 1, n_mem_ri
2791 kbounds(:, 1) = [batch_start_ri(k_mem), batch_end_ri(k_mem)]
2792
2793 bounds_cpy(:, 1) = [batch_start_ri(k_mem), batch_end_ri(k_mem)]
2794 bounds_cpy(:, 2) = [batch_start(i_mem), batch_end(i_mem)]
2795 bounds_cpy(:, 3) = [batch_start(j_mem), batch_end(j_mem)]
2796 CALL dbt_copy(t_3c_sparse, t_3c_4, bounds=bounds_cpy)
2797
2798 !Contract with the 2-center product S^-1 * V * S^-1 while keeping sparsity of derivatives
2799 CALL dbt_batched_contract_init(t_svs)
2800 CALL dbt_contract(1.0_dp, t_svs, t_3c_3, 0.0_dp, t_3c_4, &
2801 contract_1=[2], notcontract_1=[1], &
2802 contract_2=[1], notcontract_2=[2, 3], &
2803 map_1=[1], map_2=[2, 3], filter_eps=ri_data%filter_eps, &
2804 retain_sparsity=.true., unit_nr=unit_nr_dbcsr, flop=nflop)
2805 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2806 CALL dbt_batched_contract_finalize(t_svs)
2807
2808 CALL dbt_copy(t_3c_4, t_3c_5, summation=.true., move_data=.true.)
2809
2810 ijbounds(:, 1) = ibounds(:, 1)
2811 ijbounds(:, 2) = jbounds(:, 1)
2812
2813 !Contract R_PS = (acP) M_acS
2814 CALL dbt_batched_contract_init(t_r)
2815 CALL dbt_contract(1.0_dp, t_3c_int_2, t_3c_3, 1.0_dp, t_r, &
2816 contract_1=[2, 3], notcontract_1=[1], &
2817 contract_2=[2, 3], notcontract_2=[1], &
2818 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps, &
2819 bounds_1=ijbounds, bounds_3=kbounds, &
2820 unit_nr=unit_nr_dbcsr, flop=nflop)
2821 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2822 CALL dbt_batched_contract_finalize(t_r)
2823
2824 END DO !k_mem
2825 END DO !j_mem
2826
2827 CALL dbt_copy(t_3c_5, t_3c_help_1, move_data=.true.)
2828
2829 !The force from the 3c derivatives
2830 pref = -0.5_dp*2.0_dp*hf_fraction*spin_fac
2831
2832 DO k_mem = 1, SIZE(t_3c_der_ri_comp, 1)
2833 DO i_xyz = 1, 3
2834 CALL dbt_clear(t_3c_der_ri(i_xyz))
2835 CALL decompress_tensor(t_3c_der_ri(i_xyz), t_3c_der_ri_ind(k_mem, i_xyz)%ind, &
2836 t_3c_der_ri_comp(k_mem, i_xyz), ri_data%filter_eps_storage)
2837 END DO
2838 CALL get_force_from_3c_trace(force, t_3c_help_1, t_3c_der_ri, atom_of_kind, kind_of, &
2839 idx_to_at_ri, pref)
2840 END DO
2841
2842 pref = -0.5_dp*4.0_dp*hf_fraction*spin_fac
2843 IF (do_resp) THEN
2844 pref = 0.5_dp*pref
2845 CALL dbt_copy(t_3c_help_1, t_3c_help_2, order=[1, 3, 2])
2846 END IF
2847
2848 DO k_mem = 1, SIZE(t_3c_der_ao_comp, 1)
2849 DO i_xyz = 1, 3
2850 CALL dbt_clear(t_3c_der_ao(i_xyz))
2851 CALL decompress_tensor(t_3c_der_ao(i_xyz), t_3c_der_ao_ind(k_mem, i_xyz)%ind, &
2852 t_3c_der_ao_comp(k_mem, i_xyz), ri_data%filter_eps_storage)
2853 END DO
2854 CALL get_force_from_3c_trace(force, t_3c_help_1, t_3c_der_ao, atom_of_kind, kind_of, &
2855 idx_to_at_ao, pref, deriv_dim=2)
2856
2857 IF (do_resp) THEN
2858 CALL get_force_from_3c_trace(force, t_3c_help_2, t_3c_der_ao, atom_of_kind, kind_of, &
2859 idx_to_at_ao, pref, deriv_dim=2)
2860 END IF
2861 END DO
2862
2863 !The 3c virial contribution. Note: only fraction of integrals correspondig to i_mem calculated
2864 IF (use_virial) THEN
2865 pref = -0.5_dp*2.0_dp*hf_fraction*spin_fac
2866 CALL dbt_copy(t_3c_help_1, t_3c_virial, move_data=.true.)
2867 CALL calc_3c_virial(work_virial, t_3c_virial, pref, qs_env, nl_3c, basis_set_ri, &
2868 basis_set_ao, basis_set_ao, ri_data%ri_metric, &
2869 der_eps=ri_data%eps_schwarz_forces, op_pos=1)
2870
2871 CALL dbt_clear(t_3c_virial)
2872 END IF
2873
2874 CALL dbt_clear(t_3c_help_1)
2875 CALL dbt_clear(t_3c_help_2)
2876 END DO !i_mem
2877 CALL timestop(handle)
2878
2879 CALL timeset(routinen//"_2c", handle)
2880 !Now we deal with all the 2-center quantities
2881 !Calculate S^-1 * R * S^-1
2882 CALL dbt_contract(1.0_dp, ri_data%t_2c_inv(1, 1), t_r, 0.0_dp, t_2c_ri_tmp, &
2883 contract_1=[2], notcontract_1=[1], &
2884 contract_2=[1], notcontract_2=[2], &
2885 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps, &
2886 unit_nr=unit_nr_dbcsr, flop=nflop)
2887 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2888
2889 CALL dbt_contract(1.0_dp, t_2c_ri_tmp, ri_data%t_2c_inv(1, 1), 0.0_dp, t_2c_ri, &
2890 contract_1=[2], notcontract_1=[1], &
2891 contract_2=[1], notcontract_2=[2], &
2892 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps, &
2893 unit_nr=unit_nr_dbcsr, flop=nflop)
2894 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2895
2896 !Calculate the potential contribution to the force: [S^-1*R*S^-1]_QR d/dx (Q|R)
2897 pref = 0.5_dp*hf_fraction*spin_fac
2898 IF (.NOT. ri_data%same_op) pref = -pref
2899 CALL get_2c_der_force(force, t_2c_ri, t_2c_der_ri, atom_of_kind, kind_of, idx_to_at_ri, pref)
2900
2901 !Calculate the contribution to the virial on the fly
2902 IF (use_virial_prv) THEN
2903 CALL dbt_copy(t_2c_ri, t_2c_virial)
2904 CALL dbt_copy_tensor_to_matrix(t_2c_virial, virial_trace)
2905 CALL calc_2c_virial(work_virial, virial_trace, pref, qs_env, nl_2c_pot, &
2906 basis_set_ri, basis_set_ri, ri_data%hfx_pot)
2907 END IF
2908
2909 !And that from the metric: [S^-1*R*S^-1*(Q|R)*S^-1]_UV d/dx S_UV
2910 IF (.NOT. ri_data%same_op) THEN
2911 CALL dbt_contract(1.0_dp, t_2c_ri, ri_data%t_2c_pot(1, 1), 0.0_dp, t_2c_ri_tmp, &
2912 contract_1=[2], notcontract_1=[1], &
2913 contract_2=[1], notcontract_2=[2], &
2914 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps, &
2915 unit_nr=unit_nr_dbcsr, flop=nflop)
2916 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2917
2918 CALL dbt_contract(1.0_dp, t_2c_ri_tmp, ri_data%t_2c_inv(1, 1), 0.0_dp, t_2c_ri, &
2919 contract_1=[2], notcontract_1=[1], &
2920 contract_2=[1], notcontract_2=[2], &
2921 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps, &
2922 unit_nr=unit_nr_dbcsr, flop=nflop)
2923 ri_data%dbcsr_nflop = ri_data%dbcsr_nflop + nflop
2924
2925 pref = 0.5_dp*2.0_dp*hf_fraction*spin_fac
2926 CALL get_2c_der_force(force, t_2c_ri, t_2c_der_metric, atom_of_kind, kind_of, idx_to_at_ri, pref)
2927
2928 IF (use_virial_prv) THEN
2929 CALL dbt_copy(t_2c_ri, t_2c_virial)
2930 CALL dbt_copy_tensor_to_matrix(t_2c_virial, virial_trace)
2931 CALL calc_2c_virial(work_virial, virial_trace, pref, qs_env, nl_2c_met, &
2932 basis_set_ri, basis_set_ri, ri_data%ri_metric)
2933 END IF
2934 END IF
2935 CALL dbt_clear(t_2c_ri)
2936 CALL dbt_clear(t_2c_ri_tmp)
2937 CALL dbt_clear(t_r)
2938 CALL dbt_clear(t_3c_help_1)
2939 CALL dbt_clear(t_3c_help_2)
2940 CALL timestop(handle)
2941
2942 IF (use_virial_prv) THEN
2943 DO k_xyz = 1, 3
2944 DO j_xyz = 1, 3
2945 DO i_xyz = 1, 3
2946 virial%pv_fock_4c(i_xyz, j_xyz) = virial%pv_fock_4c(i_xyz, j_xyz) &
2947 + work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
2948 END DO
2949 END DO
2950 END DO
2951 END IF
2952 END DO !i_spin
2953
2954 CALL dbt_batched_contract_finalize(t_3c_int)
2955 CALL dbt_batched_contract_finalize(t_3c_int_2)
2956 CALL dbt_batched_contract_finalize(t_3c_1)
2957 CALL dbt_batched_contract_finalize(t_3c_2)
2958 CALL dbt_batched_contract_finalize(t_3c_3)
2959 CALL dbt_batched_contract_finalize(t_3c_4)
2960 CALL dbt_batched_contract_finalize(t_3c_5)
2961 CALL dbt_batched_contract_finalize(t_3c_sparse)
2962
2963 CALL para_env%sync()
2964 t2 = m_walltime()
2965
2966 CALL dbt_copy(t_3c_int, ri_data%t_3c_int_ctr_2(1, 1), move_data=.true.)
2967
2968 !clean-up
2969 CALL dbt_destroy(rho_ao_1)
2970 CALL dbt_destroy(rho_ao_2)
2971 CALL dbt_destroy(t_3c_ao_ri_ao)
2972 CALL dbt_destroy(t_3c_ri_ao_ao)
2973 CALL dbt_destroy(t_3c_int)
2974 CALL dbt_destroy(t_3c_int_2)
2975 CALL dbt_destroy(t_3c_1)
2976 CALL dbt_destroy(t_3c_2)
2977 CALL dbt_destroy(t_3c_3)
2978 CALL dbt_destroy(t_3c_4)
2979 CALL dbt_destroy(t_3c_5)
2980 CALL dbt_destroy(t_3c_help_1)
2981 CALL dbt_destroy(t_3c_help_2)
2982 CALL dbt_destroy(t_3c_sparse)
2983 CALL dbt_destroy(t_svs)
2984 CALL dbt_destroy(t_r)
2985 CALL dbt_destroy(t_2c_ri)
2986 CALL dbt_destroy(t_2c_ri_tmp)
2987
2988 DO i_xyz = 1, 3
2989 CALL dbt_destroy(t_3c_der_ri(i_xyz))
2990 CALL dbt_destroy(t_3c_der_ao(i_xyz))
2991 CALL dbt_destroy(t_2c_der_ri(i_xyz))
2992 IF (.NOT. ri_data%same_op) CALL dbt_destroy(t_2c_der_metric(i_xyz))
2993 END DO
2994
2995 DO i_xyz = 1, 3
2996 DO i_mem = 1, SIZE(t_3c_der_ao_comp, 1)
2997 CALL dealloc_containers(t_3c_der_ao_comp(i_mem, i_xyz), dummy_int)
2998 CALL dealloc_containers(t_3c_der_ri_comp(i_mem, i_xyz), dummy_int)
2999 END DO
3000 END DO
3001 DEALLOCATE (t_3c_der_ao_ind, t_3c_der_ri_ind)
3002
3003 DO ibasis = 1, SIZE(basis_set_ao)
3004 orb_basis => basis_set_ao(ibasis)%gto_basis_set
3005 ri_basis => basis_set_ri(ibasis)%gto_basis_set
3006 CALL init_interaction_radii_orb_basis(orb_basis, dft_control%qs_control%eps_pgf_orb)
3007 CALL init_interaction_radii_orb_basis(ri_basis, dft_control%qs_control%eps_pgf_orb)
3008 END DO
3009
3010 IF (use_virial) THEN
3011 CALL release_neighbor_list_sets(nl_2c_met)
3012 CALL release_neighbor_list_sets(nl_2c_pot)
3013 CALL neighbor_list_3c_destroy(nl_3c)
3014 CALL dbcsr_release(virial_trace)
3015 CALL dbt_destroy(t_2c_virial)
3016 CALL dbt_destroy(t_3c_virial)
3017 END IF
3018
3019 END SUBROUTINE hfx_ri_forces_pmat
3020
3021! **************************************************************************************************
3022!> \brief the general routine that calls the relevant force code
3023!> \param qs_env ...
3024!> \param ri_data ...
3025!> \param nspins ...
3026!> \param hf_fraction ...
3027!> \param rho_ao ...
3028!> \param rho_ao_resp ...
3029!> \param mos ...
3030!> \param use_virial ...
3031!> \param resp_only ...
3032!> \param rescale_factor ...
3033! **************************************************************************************************
3034 SUBROUTINE hfx_ri_update_forces(qs_env, ri_data, nspins, hf_fraction, rho_ao, rho_ao_resp, &
3035 mos, use_virial, resp_only, rescale_factor)
3036
3037 TYPE(qs_environment_type), POINTER :: qs_env
3038 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
3039 INTEGER, INTENT(IN) :: nspins
3040 REAL(kind=dp), INTENT(IN) :: hf_fraction
3041 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL :: rho_ao
3042 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL :: rho_ao_resp
3043 TYPE(mo_set_type), DIMENSION(:), INTENT(IN), &
3044 OPTIONAL :: mos
3045 LOGICAL, INTENT(IN), OPTIONAL :: use_virial, resp_only
3046 REAL(dp), INTENT(IN), OPTIONAL :: rescale_factor
3047
3048 CHARACTER(LEN=*), PARAMETER :: routinen = 'hfx_ri_update_forces'
3049
3050 INTEGER :: handle, ispin
3051 INTEGER, DIMENSION(2) :: homo
3052 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
3053 TYPE(cp_fm_type), POINTER :: mo_coeff
3054 TYPE(dbcsr_type), DIMENSION(2) :: mo_coeff_b
3055 TYPE(dbcsr_type), POINTER :: mo_coeff_b_tmp
3056
3057 CALL timeset(routinen, handle)
3058
3059 SELECT CASE (ri_data%flavor)
3060 CASE (ri_mo)
3061
3062 DO ispin = 1, nspins
3063 NULLIFY (mo_coeff_b_tmp)
3064 cpassert(mos(ispin)%uniform_occupation)
3065 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, eigenvalues=mo_eigenvalues, mo_coeff_b=mo_coeff_b_tmp)
3066
3067 IF (.NOT. mos(ispin)%use_mo_coeff_b) CALL copy_fm_to_dbcsr(mo_coeff, mo_coeff_b_tmp)
3068 CALL dbcsr_copy(mo_coeff_b(ispin), mo_coeff_b_tmp)
3069 END DO
3070
3071 DO ispin = 1, nspins
3072 CALL dbcsr_scale(mo_coeff_b(ispin), sqrt(mos(ispin)%maxocc))
3073 homo(ispin) = mos(ispin)%homo
3074 END DO
3075
3076 CALL hfx_ri_forces_mo(qs_env, ri_data, nspins, hf_fraction, mo_coeff_b, use_virial)
3077
3078 CASE (ri_pmat)
3079
3080 CALL hfx_ri_forces_pmat(qs_env, ri_data, nspins, hf_fraction, rho_ao, rho_ao_resp, use_virial, &
3081 resp_only, rescale_factor)
3082 END SELECT
3083
3084 DO ispin = 1, nspins
3085 CALL dbcsr_release(mo_coeff_b(ispin))
3086 END DO
3087
3088 CALL timestop(handle)
3089
3090 END SUBROUTINE hfx_ri_update_forces
3091
3092! **************************************************************************************************
3093!> \brief Calculate the derivatives tensors for the force, in a format fit for contractions
3094!> \param t_3c_der_RI_comp compressed RI derivatives
3095!> \param t_3c_der_AO_comp compressed AO derivatives
3096!> \param t_3c_der_RI_ind ...
3097!> \param t_3c_der_AO_ind ...
3098!> \param t_2c_der_RI format based on standard atomic block sizes
3099!> \param t_2c_der_metric format based on standard atomic block sizes
3100!> \param ri_ao_ao_template ...
3101!> \param basis_set_AO ...
3102!> \param basis_set_RI ...
3103!> \param ri_data ...
3104!> \param qs_env ...
3105!> \param nl_2c_pot ...
3106!> \param nl_2c_met ...
3107!> \param nl_3c_out ...
3108!> \param t_3c_virial ...
3109! **************************************************************************************************
3110 SUBROUTINE precalc_derivatives(t_3c_der_RI_comp, t_3c_der_AO_comp, t_3c_der_RI_ind, t_3c_der_AO_ind, &
3111 t_2c_der_RI, t_2c_der_metric, ri_ao_ao_template, &
3112 basis_set_AO, basis_set_RI, ri_data, qs_env, &
3113 nl_2c_pot, nl_2c_met, nl_3c_out, t_3c_virial)
3114
3115 TYPE(hfx_compression_type), ALLOCATABLE, &
3116 DIMENSION(:, :), INTENT(INOUT) :: t_3c_der_ri_comp, t_3c_der_ao_comp
3117 TYPE(block_ind_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_der_ri_ind, t_3c_der_ao_ind
3118 TYPE(dbt_type), DIMENSION(3), INTENT(OUT) :: t_2c_der_ri, t_2c_der_metric
3119 TYPE(dbt_type), INTENT(INOUT) :: ri_ao_ao_template
3120 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
3121 DIMENSION(:), TARGET :: basis_set_ao, basis_set_ri
3122 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
3123 TYPE(qs_environment_type), POINTER :: qs_env
3124 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3125 OPTIONAL, POINTER :: nl_2c_pot, nl_2c_met
3126 TYPE(neighbor_list_3c_type), OPTIONAL :: nl_3c_out
3127 TYPE(dbt_type), INTENT(INOUT), OPTIONAL :: t_3c_virial
3128
3129 CHARACTER(LEN=*), PARAMETER :: routinen = 'precalc_derivatives'
3130
3131 INTEGER :: handle, i_mem, i_xyz, n_mem, natom, &
3132 nkind, nthreads
3133 INTEGER(int_8) :: nze, nze_tot
3134 INTEGER, ALLOCATABLE, DIMENSION(:) :: dist1, dist2, dist_ao_1, dist_ao_2, &
3135 dist_ri, dummy_end, dummy_start, &
3136 end_blocks, start_blocks
3137 INTEGER, DIMENSION(3) :: pcoord, pdims
3138 INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
3139 REAL(dp) :: compression_factor, memory, occ
3140 TYPE(dbcsr_distribution_type) :: dbcsr_dist
3141 TYPE(dbcsr_type), DIMENSION(1, 3) :: t_2c_der_metric_prv, t_2c_der_ri_prv
3142 TYPE(dbt_pgrid_type) :: pgrid
3143 TYPE(dbt_type) :: t_2c_template, t_2c_tmp, t_3c_template
3144 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :, :) :: t_3c_der_ao_prv, t_3c_der_ri_prv
3145 TYPE(dft_control_type), POINTER :: dft_control
3146 TYPE(distribution_2d_type), POINTER :: dist_2d
3147 TYPE(distribution_3d_type) :: dist_3d, dist_3d_out
3148 TYPE(mp_cart_type) :: mp_comm_t3c, mp_comm_t3c_out
3149 TYPE(mp_para_env_type), POINTER :: para_env
3150 TYPE(neighbor_list_3c_type) :: nl_3c
3151 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3152 POINTER :: nl_2c
3153 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3154 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3155
3156 NULLIFY (qs_kind_set, dist_2d, nl_2c, particle_set, dft_control, para_env)
3157
3158 CALL timeset(routinen, handle)
3159
3160 CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, distribution_2d=dist_2d, natom=natom, &
3161 particle_set=particle_set, dft_control=dft_control, para_env=para_env)
3162
3163 !TODO: is such a pgrid correct?
3164 !Dealing with the 3c derivatives
3165 nthreads = 1
3166!$ nthreads = omp_get_num_threads()
3167 pdims = 0
3168 CALL dbt_pgrid_create(para_env, pdims, pgrid, tensor_dims=[max(1, natom/(ri_data%n_mem*nthreads)), natom, natom])
3169
3170 CALL create_3c_tensor(t_3c_template, dist_ri, dist_ao_1, dist_ao_2, pgrid, &
3171 ri_data%bsizes_RI, ri_data%bsizes_AO, ri_data%bsizes_AO, &
3172 map1=[1], map2=[2, 3], &
3173 name="der (RI AO | AO)")
3174
3175 ALLOCATE (t_3c_der_ao_prv(1, 1, 3), t_3c_der_ri_prv(1, 1, 3))
3176 DO i_xyz = 1, 3
3177 CALL dbt_create(t_3c_template, t_3c_der_ri_prv(1, 1, i_xyz))
3178 CALL dbt_create(t_3c_template, t_3c_der_ao_prv(1, 1, i_xyz))
3179 END DO
3180 IF (PRESENT(t_3c_virial)) THEN
3181 CALL dbt_create(t_3c_template, t_3c_virial)
3182 END IF
3183 CALL dbt_destroy(t_3c_template)
3184
3185 CALL dbt_mp_environ_pgrid(pgrid, pdims, pcoord)
3186 CALL mp_comm_t3c%create(pgrid%mp_comm_2d, 3, pdims)
3187 CALL distribution_3d_create(dist_3d, dist_ri, dist_ao_1, dist_ao_2, &
3188 nkind, particle_set, mp_comm_t3c, own_comm=.true.)
3189
3190 CALL build_3c_neighbor_lists(nl_3c, basis_set_ri, basis_set_ao, basis_set_ao, dist_3d, ri_data%ri_metric, &
3191 "HFX_3c_nl", qs_env, op_pos=1, sym_jk=.true., own_dist=.true.)
3192
3193 IF (PRESENT(nl_3c_out)) THEN
3194 CALL dbt_mp_environ_pgrid(pgrid, pdims, pcoord)
3195 CALL mp_comm_t3c_out%create(pgrid%mp_comm_2d, 3, pdims)
3196 CALL distribution_3d_create(dist_3d_out, dist_ri, dist_ao_1, dist_ao_2, &
3197 nkind, particle_set, mp_comm_t3c_out, own_comm=.true.)
3198 CALL build_3c_neighbor_lists(nl_3c_out, basis_set_ri, basis_set_ao, basis_set_ao, dist_3d_out, &
3199 ri_data%ri_metric, "HFX_3c_nl", qs_env, op_pos=1, sym_jk=.false., &
3200 own_dist=.true.)
3201 END IF
3202 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
3203 CALL dbt_pgrid_destroy(pgrid)
3204
3205 n_mem = ri_data%n_mem
3206 CALL create_tensor_batches(ri_data%bsizes_RI, n_mem, dummy_start, dummy_end, &
3207 start_blocks, end_blocks)
3208 DEALLOCATE (dummy_start, dummy_end)
3209
3210 ALLOCATE (t_3c_der_ao_comp(n_mem, 3), t_3c_der_ri_comp(n_mem, 3))
3211 ALLOCATE (t_3c_der_ao_ind(n_mem, 3), t_3c_der_ri_ind(n_mem, 3))
3212
3213 memory = 0.0_dp
3214 nze_tot = 0
3215 DO i_mem = 1, n_mem
3216 CALL build_3c_derivatives(t_3c_der_ri_prv, t_3c_der_ao_prv, ri_data%filter_eps, qs_env, &
3217 nl_3c, basis_set_ri, basis_set_ao, basis_set_ao, &
3218 ri_data%ri_metric, der_eps=ri_data%eps_schwarz_forces, op_pos=1, &
3219 bounds_i=[start_blocks(i_mem), end_blocks(i_mem)])
3220
3221 DO i_xyz = 1, 3
3222 CALL dbt_copy(t_3c_der_ri_prv(1, 1, i_xyz), ri_ao_ao_template, move_data=.true.)
3223 CALL dbt_filter(ri_ao_ao_template, ri_data%filter_eps)
3224 CALL get_tensor_occupancy(ri_ao_ao_template, nze, occ)
3225 nze_tot = nze_tot + nze
3226
3227 CALL alloc_containers(t_3c_der_ri_comp(i_mem, i_xyz), 1)
3228 CALL compress_tensor(ri_ao_ao_template, t_3c_der_ri_ind(i_mem, i_xyz)%ind, &
3229 t_3c_der_ri_comp(i_mem, i_xyz), ri_data%filter_eps_storage, memory)
3230 CALL dbt_clear(ri_ao_ao_template)
3231
3232 !put AO derivative as middle index
3233 CALL dbt_copy(t_3c_der_ao_prv(1, 1, i_xyz), ri_ao_ao_template, order=[1, 3, 2], move_data=.true.)
3234 CALL dbt_filter(ri_ao_ao_template, ri_data%filter_eps)
3235 CALL get_tensor_occupancy(ri_ao_ao_template, nze, occ)
3236 nze_tot = nze_tot + nze
3237
3238 CALL alloc_containers(t_3c_der_ao_comp(i_mem, i_xyz), 1)
3239 CALL compress_tensor(ri_ao_ao_template, t_3c_der_ao_ind(i_mem, i_xyz)%ind, &
3240 t_3c_der_ao_comp(i_mem, i_xyz), ri_data%filter_eps_storage, memory)
3241 CALL dbt_clear(ri_ao_ao_template)
3242 END DO
3243 END DO
3244
3245 CALL neighbor_list_3c_destroy(nl_3c)
3246 DO i_xyz = 1, 3
3247 CALL dbt_destroy(t_3c_der_ri_prv(1, 1, i_xyz))
3248 CALL dbt_destroy(t_3c_der_ao_prv(1, 1, i_xyz))
3249 END DO
3250
3251 CALL para_env%sum(memory)
3252 compression_factor = real(nze_tot, dp)*1.0e-06*8.0_dp/memory
3253 IF (ri_data%unit_nr > 0) THEN
3254 WRITE (unit=ri_data%unit_nr, fmt="((T3,A,T66,F11.2,A4))") &
3255 "MEMORY_INFO| Memory for 3-center HFX derivatives (compressed):", memory, ' MiB'
3256
3257 WRITE (unit=ri_data%unit_nr, fmt="((T3,A,T60,F21.2))") &
3258 "MEMORY_INFO| Compression factor: ", compression_factor
3259 END IF
3260
3261 !Deal with the 2-center derivatives
3262 CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
3263 ALLOCATE (row_bsize(SIZE(ri_data%bsizes_RI)))
3264 ALLOCATE (col_bsize(SIZE(ri_data%bsizes_RI)))
3265 row_bsize(:) = ri_data%bsizes_RI
3266 col_bsize(:) = ri_data%bsizes_RI
3267
3268 CALL build_2c_neighbor_lists(nl_2c, basis_set_ri, basis_set_ri, ri_data%hfx_pot, &
3269 "HFX_2c_nl_pot", qs_env, sym_ij=.true., dist_2d=dist_2d)
3270
3271 DO i_xyz = 1, 3
3272 CALL dbcsr_create(t_2c_der_ri_prv(1, i_xyz), "(R|P) HFX der", dbcsr_dist, &
3273 dbcsr_type_antisymmetric, row_bsize, col_bsize)
3274 END DO
3275
3276 CALL build_2c_derivatives(t_2c_der_ri_prv, ri_data%filter_eps_2c, qs_env, nl_2c, basis_set_ri, &
3277 basis_set_ri, ri_data%hfx_pot)
3278 CALL release_neighbor_list_sets(nl_2c)
3279
3280 IF (PRESENT(nl_2c_pot)) THEN
3281 NULLIFY (nl_2c_pot)
3282 CALL build_2c_neighbor_lists(nl_2c_pot, basis_set_ri, basis_set_ri, ri_data%hfx_pot, &
3283 "HFX_2c_nl_pot", qs_env, sym_ij=.false., dist_2d=dist_2d)
3284 END IF
3285
3286 !copy 2c derivative tensor into the standard format
3287 CALL create_2c_tensor(t_2c_template, dist1, dist2, ri_data%pgrid_2d, ri_data%bsizes_RI_split, &
3288 ri_data%bsizes_RI_split, name='(RI| RI)')
3289 DEALLOCATE (dist1, dist2)
3290
3291 DO i_xyz = 1, 3
3292 CALL dbt_create(t_2c_der_ri_prv(1, i_xyz), t_2c_tmp)
3293 CALL dbt_copy_matrix_to_tensor(t_2c_der_ri_prv(1, i_xyz), t_2c_tmp)
3294
3295 CALL dbt_create(t_2c_template, t_2c_der_ri(i_xyz))
3296 CALL dbt_copy(t_2c_tmp, t_2c_der_ri(i_xyz), move_data=.true.)
3297
3298 CALL dbt_destroy(t_2c_tmp)
3299 CALL dbcsr_release(t_2c_der_ri_prv(1, i_xyz))
3300 END DO
3301
3302 !Repeat with the metric, if required
3303 IF (.NOT. ri_data%same_op) THEN
3304
3305 CALL build_2c_neighbor_lists(nl_2c, basis_set_ri, basis_set_ri, ri_data%ri_metric, &
3306 "HFX_2c_nl_RI", qs_env, sym_ij=.true., dist_2d=dist_2d)
3307
3308 DO i_xyz = 1, 3
3309 CALL dbcsr_create(t_2c_der_metric_prv(1, i_xyz), "(R|P) HFX der", dbcsr_dist, &
3310 dbcsr_type_antisymmetric, row_bsize, col_bsize)
3311 END DO
3312
3313 CALL build_2c_derivatives(t_2c_der_metric_prv, ri_data%filter_eps_2c, qs_env, nl_2c, &
3314 basis_set_ri, basis_set_ri, ri_data%ri_metric)
3315 CALL release_neighbor_list_sets(nl_2c)
3316
3317 IF (PRESENT(nl_2c_met)) THEN
3318 NULLIFY (nl_2c_met)
3319 CALL build_2c_neighbor_lists(nl_2c_met, basis_set_ri, basis_set_ri, ri_data%ri_metric, &
3320 "HFX_2c_nl_RI", qs_env, sym_ij=.false., dist_2d=dist_2d)
3321 END IF
3322
3323 DO i_xyz = 1, 3
3324 CALL dbt_create(t_2c_der_metric_prv(1, i_xyz), t_2c_tmp)
3325 CALL dbt_copy_matrix_to_tensor(t_2c_der_metric_prv(1, i_xyz), t_2c_tmp)
3326
3327 CALL dbt_create(t_2c_template, t_2c_der_metric(i_xyz))
3328 CALL dbt_copy(t_2c_tmp, t_2c_der_metric(i_xyz), move_data=.true.)
3329
3330 CALL dbt_destroy(t_2c_tmp)
3331 CALL dbcsr_release(t_2c_der_metric_prv(1, i_xyz))
3332 END DO
3333
3334 END IF
3335
3336 CALL dbt_destroy(t_2c_template)
3337 CALL dbcsr_distribution_release(dbcsr_dist)
3338 DEALLOCATE (row_bsize, col_bsize)
3339
3340 CALL timestop(handle)
3341
3342 END SUBROUTINE precalc_derivatives
3343
3344! **************************************************************************************************
3345!> \brief This routines calculates the force contribution from a trace over 3D tensors, i.e.
3346!> force = sum_ijk A_ijk B_ijk. An iteration over the blocks is made, which index determin
3347!> the atom on which the force acts
3348!> \param force ...
3349!> \param t_3c_contr ...
3350!> \param t_3c_der ...
3351!> \param atom_of_kind ...
3352!> \param kind_of ...
3353!> \param idx_to_at ...
3354!> \param pref ...
3355!> \param do_mp2 ...
3356!> \param deriv_dim the dimension of the tensor corresponding to the derivative (by default 1)
3357! **************************************************************************************************
3358 SUBROUTINE get_force_from_3c_trace(force, t_3c_contr, t_3c_der, atom_of_kind, kind_of, idx_to_at, &
3359 pref, do_mp2, deriv_dim)
3360
3361 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3362 TYPE(dbt_type), INTENT(INOUT) :: t_3c_contr
3363 TYPE(dbt_type), DIMENSION(3), INTENT(INOUT) :: t_3c_der
3364 INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind, kind_of, idx_to_at
3365 REAL(dp), INTENT(IN) :: pref
3366 LOGICAL, INTENT(IN), OPTIONAL :: do_mp2
3367 INTEGER, INTENT(IN), OPTIONAL :: deriv_dim
3368
3369 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_force_from_3c_trace'
3370
3371 INTEGER :: deriv_dim_prv, handle, i_xyz, iat, &
3372 iat_of_kind, ikind, j_xyz
3373 INTEGER, DIMENSION(3) :: ind
3374 LOGICAL :: do_mp2_prv, found
3375 REAL(dp) :: new_force
3376 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :), TARGET :: contr_blk, der_blk
3377 REAL(dp), DIMENSION(3) :: scoord
3378 TYPE(dbt_iterator_type) :: iter
3379
3380 CALL timeset(routinen, handle)
3381
3382 do_mp2_prv = .false.
3383 IF (PRESENT(do_mp2)) do_mp2_prv = do_mp2
3384
3385 deriv_dim_prv = 1
3386 IF (PRESENT(deriv_dim)) deriv_dim_prv = deriv_dim
3387
3388!$OMP PARALLEL DEFAULT(NONE) &
3389!$OMP SHARED(t_3c_der,t_3c_contr,force,do_mp2_prv,deriv_dim_prv,pref,idx_to_at,atom_of_kind,kind_of) &
3390!$OMP PRIVATE(i_xyz,j_xyz,iter,ind,der_blk,contr_blk,found,new_force,iat,iat_of_kind,ikind,scoord)
3391 DO i_xyz = 1, 3
3392 CALL dbt_iterator_start(iter, t_3c_der(i_xyz))
3393 DO WHILE (dbt_iterator_blocks_left(iter))
3394 CALL dbt_iterator_next_block(iter, ind)
3395
3396 CALL dbt_get_block(t_3c_der(i_xyz), ind, der_blk, found)
3397 cpassert(found)
3398 CALL dbt_get_block(t_3c_contr, ind, contr_blk, found)
3399
3400 IF (found) THEN
3401
3402 !take the trace of the blocks
3403 new_force = pref*sum(der_blk(:, :, :)*contr_blk(:, :, :))
3404
3405 !the first index of the derivative tensor defines the atom
3406 iat = idx_to_at(ind(deriv_dim_prv))
3407 iat_of_kind = atom_of_kind(iat)
3408 ikind = kind_of(iat)
3409
3410 IF (.NOT. do_mp2_prv) THEN
3411!$OMP ATOMIC
3412 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) &
3413 + new_force
3414 ELSE
3415!$OMP ATOMIC
3416 force(ikind)%mp2_non_sep(i_xyz, iat_of_kind) = force(ikind)%mp2_non_sep(i_xyz, iat_of_kind) &
3417 + new_force
3418 END IF
3419
3420 DEALLOCATE (contr_blk)
3421 END IF
3422 DEALLOCATE (der_blk)
3423 END DO !iter
3424 CALL dbt_iterator_stop(iter)
3425 END DO
3426!$OMP END PARALLEL
3427 CALL timestop(handle)
3428
3429 END SUBROUTINE get_force_from_3c_trace
3430
3431! **************************************************************************************************
3432!> \brief Update the forces due to the derivative of the a 2-center product d/dR (Q|R)
3433!> \param force ...
3434!> \param t_2c_contr A precontracted tensor containing sum_abcdPS (ab|P)(P|Q)^-1 (R|S)^-1 (S|cd) P_ac P_bd
3435!> \param t_2c_der the d/dR (Q|R) tensor, in all 3 cartesian directions
3436!> \param atom_of_kind ...
3437!> \param kind_of ...
3438!> \param idx_to_at ...
3439!> \param pref ...
3440!> \param do_mp2 ...
3441!> \param do_ovlp ...
3442!> \note IMPORTANT: t_tc_contr and t_2c_der need to have the same distribution
3443! **************************************************************************************************
3444 SUBROUTINE get_2c_der_force(force, t_2c_contr, t_2c_der, atom_of_kind, kind_of, idx_to_at, &
3445 pref, do_mp2, do_ovlp)
3446
3447 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3448 TYPE(dbt_type), INTENT(INOUT) :: t_2c_contr
3449 TYPE(dbt_type), DIMENSION(3), INTENT(INOUT) :: t_2c_der
3450 INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind, kind_of, idx_to_at
3451 REAL(dp), INTENT(IN) :: pref
3452 LOGICAL, INTENT(IN), OPTIONAL :: do_mp2, do_ovlp
3453
3454 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_2c_der_force'
3455
3456 INTEGER :: handle, i_xyz, iat, iat_of_kind, ikind, &
3457 j_xyz, jat, jat_of_kind, jkind
3458 INTEGER, DIMENSION(2) :: ind
3459 LOGICAL :: do_mp2_prv, do_ovlp_prv, found
3460 REAL(dp) :: new_force
3461 REAL(dp), ALLOCATABLE, DIMENSION(:, :), TARGET :: contr_blk, der_blk
3462 REAL(dp), DIMENSION(3) :: scoord
3463 TYPE(dbt_iterator_type) :: iter
3464
3465 !Loop over the blocks of d/dR (Q|R), contract with the corresponding block of t_2c_contr and
3466 !update the relevant force
3467
3468 CALL timeset(routinen, handle)
3469
3470 do_mp2_prv = .false.
3471 IF (PRESENT(do_mp2)) do_mp2_prv = do_mp2
3472
3473 do_ovlp_prv = .false.
3474 IF (PRESENT(do_ovlp)) do_ovlp_prv = do_ovlp
3475
3476!$OMP PARALLEL DEFAULT(NONE) &
3477!$OMP SHARED(t_2c_der,t_2c_contr,force,do_mp2_prv,do_ovlp_prv,pref,idx_to_at,atom_of_kind,kind_of) &
3478!$OMP PRIVATE(i_xyz,j_xyz,iter,ind,der_blk,contr_blk,found,new_force) &
3479!$OMP PRIVATE(iat,jat,iat_of_kind,jat_of_kind,ikind,jkind,scoord)
3480 DO i_xyz = 1, 3
3481 CALL dbt_iterator_start(iter, t_2c_der(i_xyz))
3482 DO WHILE (dbt_iterator_blocks_left(iter))
3483 CALL dbt_iterator_next_block(iter, ind)
3484
3485 IF (ind(1) == ind(2)) cycle
3486
3487 CALL dbt_get_block(t_2c_der(i_xyz), ind, der_blk, found)
3488 cpassert(found)
3489 CALL dbt_get_block(t_2c_contr, ind, contr_blk, found)
3490
3491 IF (found) THEN
3492
3493 !an element of d/dR (Q|R) corresponds to 2 things because of translational invariance
3494 !(Q'| R) = - (Q| R'), once wrt the center on Q, and once on R
3495 new_force = pref*sum(der_blk(:, :)*contr_blk(:, :))
3496
3497 iat = idx_to_at(ind(1))
3498 iat_of_kind = atom_of_kind(iat)
3499 ikind = kind_of(iat)
3500
3501 IF (do_mp2_prv) THEN
3502!$OMP ATOMIC
3503 force(ikind)%mp2_non_sep(i_xyz, iat_of_kind) = force(ikind)%mp2_non_sep(i_xyz, iat_of_kind) &
3504 + new_force
3505 ELSE IF (do_ovlp_prv) THEN
3506!$OMP ATOMIC
3507 force(ikind)%overlap(i_xyz, iat_of_kind) = force(ikind)%overlap(i_xyz, iat_of_kind) &
3508 + new_force
3509 ELSE
3510!$OMP ATOMIC
3511 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) &
3512 + new_force
3513 END IF
3514
3515 jat = idx_to_at(ind(2))
3516 jat_of_kind = atom_of_kind(jat)
3517 jkind = kind_of(jat)
3518
3519 IF (do_mp2_prv) THEN
3520!$OMP ATOMIC
3521 force(jkind)%mp2_non_sep(i_xyz, jat_of_kind) = force(jkind)%mp2_non_sep(i_xyz, jat_of_kind) &
3522 - new_force
3523 ELSE IF (do_ovlp_prv) THEN
3524!$OMP ATOMIC
3525 force(jkind)%overlap(i_xyz, jat_of_kind) = force(jkind)%overlap(i_xyz, jat_of_kind) &
3526 - new_force
3527 ELSE
3528!$OMP ATOMIC
3529 force(jkind)%fock_4c(i_xyz, jat_of_kind) = force(jkind)%fock_4c(i_xyz, jat_of_kind) &
3530 - new_force
3531 END IF
3532
3533 DEALLOCATE (contr_blk)
3534 END IF
3535
3536 DEALLOCATE (der_blk)
3537 END DO !iter
3538 CALL dbt_iterator_stop(iter)
3539
3540 END DO !i_xyz
3541!$OMP END PARALLEL
3542 CALL timestop(handle)
3543
3544 END SUBROUTINE get_2c_der_force
3545
3546! **************************************************************************************************
3547!> \brief Get the force from a contraction of type SUM_a,beta (a|beta') C_a,beta, where beta is an AO
3548!> and a is a MO
3549!> \param force ...
3550!> \param t_mo_coeff ...
3551!> \param t_2c_MO_AO ...
3552!> \param atom_of_kind ...
3553!> \param kind_of ...
3554!> \param idx_to_at ...
3555!> \param pref ...
3556!> \param i_xyz ...
3557! **************************************************************************************************
3558 SUBROUTINE get_mo_ao_force(force, t_mo_coeff, t_2c_MO_AO, atom_of_kind, kind_of, idx_to_at, pref, i_xyz)
3559
3560 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3561 TYPE(dbt_type), INTENT(INOUT) :: t_mo_coeff, t_2c_mo_ao
3562 INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind, kind_of, idx_to_at
3563 REAL(dp), INTENT(IN) :: pref
3564 INTEGER, INTENT(IN) :: i_xyz
3565
3566 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_MO_AO_force'
3567
3568 INTEGER :: handle, iat, iat_of_kind, ikind, j_xyz
3569 INTEGER, DIMENSION(2) :: ind
3570 LOGICAL :: found
3571 REAL(dp) :: new_force
3572 REAL(dp), ALLOCATABLE, DIMENSION(:, :), TARGET :: mo_ao_blk, mo_coeff_blk
3573 REAL(dp), DIMENSION(3) :: scoord
3574 TYPE(dbt_iterator_type) :: iter
3575
3576 CALL timeset(routinen, handle)
3577
3578!$OMP PARALLEL DEFAULT(NONE) &
3579!$OMP SHARED(t_2c_MO_AO,t_mo_coeff,pref,force,idx_to_at,atom_of_kind,kind_of,i_xyz) &
3580!$OMP PRIVATE(iter,ind,mo_ao_blk,mo_coeff_blk,found,new_force,iat,iat_of_kind,ikind,scoord,j_xyz)
3581 CALL dbt_iterator_start(iter, t_2c_mo_ao)
3582 DO WHILE (dbt_iterator_blocks_left(iter))
3583 CALL dbt_iterator_next_block(iter, ind)
3584
3585 CALL dbt_get_block(t_2c_mo_ao, ind, mo_ao_blk, found)
3586 cpassert(found)
3587 CALL dbt_get_block(t_mo_coeff, ind, mo_coeff_blk, found)
3588
3589 IF (found) THEN
3590
3591 new_force = pref*sum(mo_ao_blk(:, :)*mo_coeff_blk(:, :))
3592
3593 iat = idx_to_at(ind(2)) !AO index is column index
3594 iat_of_kind = atom_of_kind(iat)
3595 ikind = kind_of(iat)
3596
3597!$OMP ATOMIC
3598 force(ikind)%fock_4c(i_xyz, iat_of_kind) = force(ikind)%fock_4c(i_xyz, iat_of_kind) &
3599 + new_force
3600
3601 DEALLOCATE (mo_coeff_blk)
3602 END IF
3603
3604 DEALLOCATE (mo_ao_blk)
3605 END DO !iter
3606 CALL dbt_iterator_stop(iter)
3607!$OMP END PARALLEL
3608
3609 CALL timestop(handle)
3610
3611 END SUBROUTINE get_mo_ao_force
3612
3613! **************************************************************************************************
3614!> \brief Print RI-HFX quantities, as required by the PRINT subsection
3615!> \param ri_data ...
3616!> \param qs_env ...
3617! **************************************************************************************************
3618 SUBROUTINE print_ri_hfx(ri_data, qs_env)
3619
3620 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
3621 TYPE(qs_environment_type), POINTER :: qs_env
3622
3623 INTEGER :: i_ri, ibasis, nkind, nspins, unit_nr
3624 INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
3625 LOGICAL :: mult_by_s, print_density, print_ri_metric
3626 REAL(dp), ALLOCATABLE, DIMENSION(:) :: density_coeffs, density_coeffs_2
3627 TYPE(cp_blacs_env_type), POINTER :: blacs_env
3628 TYPE(cp_fm_struct_type), POINTER :: fm_struct
3629 TYPE(cp_fm_type) :: matrix_s_fm
3630 TYPE(cp_logger_type), POINTER :: logger
3631 TYPE(dbcsr_distribution_type) :: dbcsr_dist
3632 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
3633 TYPE(dbcsr_type), DIMENSION(1) :: matrix_s
3634 TYPE(dft_control_type), POINTER :: dft_control
3635 TYPE(distribution_2d_type), POINTER :: dist_2d
3636 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
3637 DIMENSION(:), TARGET :: basis_set_ao, basis_set_ri
3638 TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
3639 TYPE(mp_para_env_type), POINTER :: para_env
3640 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3641 POINTER :: nl_2c
3642 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3643 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3644 TYPE(qs_rho_type), POINTER :: rho
3645 TYPE(section_vals_type), POINTER :: input, print_section
3646
3647 NULLIFY (rho_ao, input, print_section, logger, rho, particle_set, qs_kind_set, ri_basis, nl_2c, &
3648 dist_2d, col_bsize, row_bsize, para_env, blacs_env, fm_struct, orb_basis, dft_control)
3649
3650 CALL get_qs_env(qs_env, input=input, dft_control=dft_control)
3651 logger => cp_get_default_logger()
3652 print_density = .false.
3653 print_ri_metric = .false.
3654
3655 !Do we print the RI density coeffs and/or RI_metric 2c integrals?
3656 print_section => section_vals_get_subs_vals(input, "DFT%XC%HF%RI%PRINT")
3657 IF (btest(cp_print_key_should_output(logger%iter_info, print_section, "RI_DENSITY_COEFFS"), &
3658 cp_p_file)) print_density = .true.
3659 IF (btest(cp_print_key_should_output(logger%iter_info, print_section, "RI_METRIC_2C_INTS"), &
3660 cp_p_file)) print_ri_metric = .true.
3661
3662 !common stuff
3663 IF (print_density .OR. print_ri_metric) THEN
3664
3665 !Re-calculate the 2-center RI_metric integrals (because not stored and cheap)
3666 !Recalculated the RI_metric 2c-integrals, as it is cheap, and not stored
3667 CALL get_qs_env(qs_env, nkind=nkind, qs_kind_set=qs_kind_set, particle_set=particle_set, &
3668 distribution_2d=dist_2d, para_env=para_env, blacs_env=blacs_env)
3669 ALLOCATE (basis_set_ri(nkind), basis_set_ao(nkind))
3670 CALL basis_set_list_setup(basis_set_ri, ri_data%ri_basis_type, qs_kind_set)
3671 CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_ri)
3672 CALL basis_set_list_setup(basis_set_ao, ri_data%orb_basis_type, qs_kind_set)
3673 CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_ao)
3674
3675 DO ibasis = 1, nkind
3676 ri_basis => basis_set_ri(ibasis)%gto_basis_set
3677 CALL init_interaction_radii_orb_basis(ri_basis, ri_data%eps_pgf_orb)
3678 orb_basis => basis_set_ao(ibasis)%gto_basis_set
3679 CALL init_interaction_radii_orb_basis(orb_basis, ri_data%eps_pgf_orb)
3680 END DO
3681
3682 CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
3683 ALLOCATE (row_bsize(SIZE(ri_data%bsizes_RI)))
3684 ALLOCATE (col_bsize(SIZE(ri_data%bsizes_RI)))
3685 row_bsize(:) = ri_data%bsizes_RI
3686 col_bsize(:) = ri_data%bsizes_RI
3687
3688 CALL dbcsr_create(matrix_s(1), "RI metric", dbcsr_dist, dbcsr_type_symmetric, row_bsize, col_bsize)
3689
3690 CALL build_2c_neighbor_lists(nl_2c, basis_set_ri, basis_set_ri, ri_data%ri_metric, &
3691 "HFX_2c_nl_pot", qs_env, sym_ij=.true., dist_2d=dist_2d)
3692 CALL build_2c_integrals(matrix_s, ri_data%filter_eps_2c, qs_env, nl_2c, basis_set_ri, &
3693 basis_set_ri, ri_data%ri_metric)
3694
3695 CALL release_neighbor_list_sets(nl_2c)
3696 CALL dbcsr_distribution_release(dbcsr_dist)
3697 END IF
3698
3699 IF (print_density) THEN
3700 CALL get_qs_env(qs_env, rho=rho)
3701 CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
3702 nspins = SIZE(rho_ao, 1)
3703
3704 CALL section_vals_val_get(print_section, "RI_DENSITY_COEFFS%MULTIPLY_BY_RI_2C_INTEGRALS", l_val=mult_by_s)
3705
3706 CALL get_ri_density_coeffs(density_coeffs, matrix_s(1), rho_ao, 1, basis_set_ao, basis_set_ri, &
3707 mult_by_s, ri_data, qs_env)
3708 IF (nspins == 2) &
3709 CALL get_ri_density_coeffs(density_coeffs_2, matrix_s(1), rho_ao, 2, basis_set_ao, &
3710 basis_set_ri, mult_by_s, ri_data, qs_env)
3711
3712 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%XC%HF%RI%PRINT%RI_DENSITY_COEFFS", &
3713 extension=".dat", file_status="REPLACE", &
3714 file_action="WRITE", file_form="FORMATTED")
3715
3716 IF (unit_nr > 0) THEN
3717 IF (nspins == 1) THEN
3718 WRITE (unit_nr, fmt="(A,A,A)") &
3719 "# Coefficients of the electronic density projected on the RI_HFX basis for ", &
3720 trim(logger%iter_info%project_name), " project"
3721 DO i_ri = 1, SIZE(density_coeffs)
3722 WRITE (unit_nr, fmt="(F20.12)") density_coeffs(i_ri)
3723 END DO
3724 ELSE
3725 WRITE (unit_nr, fmt="(A,A,A)") &
3726 "# Coefficients of the electronic density projected on the RI_HFX basis for ", &
3727 trim(logger%iter_info%project_name), " project. Spin up, spin down"
3728 DO i_ri = 1, SIZE(density_coeffs)
3729 WRITE (unit_nr, fmt="(F20.12,F20.12)") density_coeffs(i_ri), density_coeffs_2(i_ri)
3730 END DO
3731 END IF
3732 END IF
3733
3734 CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%XC%HF%RI%PRINT%RI_DENSITY_COEFFS")
3735 END IF
3736
3737 IF (print_ri_metric) THEN
3738
3739 !convert 2c integrals to fm for dumping
3740 CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env, &
3741 nrow_global=sum(row_bsize), ncol_global=sum(col_bsize))
3742 CALL cp_fm_create(matrix_s_fm, fm_struct)
3743
3744 CALL copy_dbcsr_to_fm(matrix_s(1), matrix_s_fm)
3745
3746 unit_nr = cp_print_key_unit_nr(logger, input, "DFT%XC%HF%RI%PRINT%RI_METRIC_2C_INTS", &
3747 extension=".fm", file_status="REPLACE", &
3748 file_action="WRITE", file_form="UNFORMATTED")
3749 CALL cp_fm_write_unformatted(matrix_s_fm, unit_nr)
3750
3751 CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%XC%HF%RI%PRINT%RI_METRIC_2C_INTS")
3752
3753 CALL cp_fm_struct_release(fm_struct)
3754 CALL cp_fm_release(matrix_s_fm)
3755 END IF
3756
3757 !clean-up
3758 IF (print_density .OR. print_ri_metric) THEN
3759 DO ibasis = 1, nkind
3760 ri_basis => basis_set_ri(ibasis)%gto_basis_set
3761 CALL init_interaction_radii_orb_basis(ri_basis, dft_control%qs_control%eps_pgf_orb)
3762 orb_basis => basis_set_ao(ibasis)%gto_basis_set
3763 CALL init_interaction_radii_orb_basis(orb_basis, dft_control%qs_control%eps_pgf_orb)
3764 END DO
3765
3766 CALL dbcsr_release(matrix_s(1))
3767 DEALLOCATE (row_bsize, col_bsize)
3768 END IF
3769
3770 END SUBROUTINE print_ri_hfx
3771
3772! **************************************************************************************************
3773!> \brief Projects the density on the RI basis and return the array of the RI coefficients
3774!> \param density_coeffs ...
3775!> \param ri_2c_ints ...
3776!> \param rho_ao ...
3777!> \param ispin ...
3778!> \param basis_set_AO ...
3779!> \param basis_set_RI ...
3780!> \param multiply_by_S ...
3781!> \param ri_data ...
3782!> \param qs_env ...
3783! **************************************************************************************************
3784 SUBROUTINE get_ri_density_coeffs(density_coeffs, ri_2c_ints, rho_ao, ispin, basis_set_AO, &
3785 basis_set_RI, multiply_by_S, ri_data, qs_env)
3786
3787 REAL(dp), ALLOCATABLE, DIMENSION(:) :: density_coeffs
3788 TYPE(dbcsr_type), INTENT(INOUT) :: ri_2c_ints
3789 TYPE(dbcsr_p_type), DIMENSION(:, :) :: rho_ao
3790 INTEGER, INTENT(IN) :: ispin
3791 TYPE(gto_basis_set_p_type), DIMENSION(:) :: basis_set_ao, basis_set_ri
3792 LOGICAL, INTENT(IN) :: multiply_by_s
3793 TYPE(hfx_ri_type), INTENT(INOUT) :: ri_data
3794 TYPE(qs_environment_type), POINTER :: qs_env
3795
3796 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_RI_density_coeffs'
3797
3798 INTEGER :: a, b, handle, i_mem, idx, n_dependent, &
3799 n_mem, n_mem_ri, natom, &
3800 nblk_per_thread, nblks, nkind
3801 INTEGER(int_8) :: nze
3802 INTEGER, ALLOCATABLE, DIMENSION(:) :: batch_block_end, batch_block_start, &
3803 dist1, dist2, dist3, dummy1, dummy2, &
3804 idx1, idx2, idx3
3805 INTEGER, DIMENSION(2) :: ind, pdims_2d
3806 INTEGER, DIMENSION(2, 3) :: bounds_cpy
3807 INTEGER, DIMENSION(3) :: dims_3c, pcoord_3d, pdims_3d
3808 LOGICAL :: calc_ints, found
3809 REAL(dp) :: occ, threshold
3810 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: blk
3811 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: blk_3d
3812 TYPE(cp_blacs_env_type), POINTER :: blacs_env
3813 TYPE(dbcsr_type) :: ri_2c_inv
3814 TYPE(dbt_distribution_type) :: dist_2d, dist_3d
3815 TYPE(dbt_iterator_type) :: iter
3816 TYPE(dbt_pgrid_type) :: pgrid_2d, pgrid_3d
3817 TYPE(dbt_type) :: density_coeffs_t, density_tmp, rho_ao_t, &
3818 rho_ao_t_3d, rho_ao_tmp, t2c_ri_ints, &
3819 t2c_ri_inv, t2c_ri_tmp
3820 TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_int_batched
3821 TYPE(dft_control_type), POINTER :: dft_control
3822 TYPE(distribution_3d_type) :: dist_nl3c
3823 TYPE(mp_cart_type) :: mp_comm_t3c
3824 TYPE(mp_para_env_type), POINTER :: para_env
3825 TYPE(neighbor_list_3c_type) :: nl_3c
3826 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3827 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3828
3829 NULLIFY (dft_control, para_env, blacs_env, particle_set, qs_kind_set)
3830
3831 CALL timeset(routinen, handle)
3832
3833 ! Projection of the density on the RI basis: n(r) = sum_pq sum_munu P_pq (pq|mu) (mu|nu)^-1 nu(r)
3834 ! = sum_nu d_nu nu(r)
3835 ! the (pq|mu) (mu|nu)^-1 contraction is already stored in compressed format
3836
3837 IF (.NOT. ri_data%flavor == ri_pmat) THEN
3838 cpabort("Can only calculate and print the RI density coefficients within the RHO flavor of RI-HFX")
3839 END IF
3840
3841 CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env, blacs_env=blacs_env, nkind=nkind, &
3842 particle_set=particle_set, qs_kind_set=qs_kind_set, natom=natom)
3843 n_mem = ri_data%n_mem
3844 n_mem_ri = ri_data%n_mem_RI
3845
3846 ! The RI 2c int tensor and its inverse
3847 CALL dbcsr_create(ri_2c_inv, template=ri_2c_ints, matrix_type=dbcsr_type_no_symmetry)
3848
3849 SELECT CASE (ri_data%t2c_method)
3850 CASE (hfx_ri_do_2c_iter)
3851 threshold = max(ri_data%filter_eps, 1.0e-12_dp)
3852 CALL invert_hotelling(ri_2c_inv, ri_2c_ints, threshold=threshold, silent=.false.)
3853 CASE (hfx_ri_do_2c_cholesky)
3854 CALL dbcsr_copy(ri_2c_inv, ri_2c_ints)
3855 CALL cp_dbcsr_cholesky_decompose(ri_2c_inv, para_env=para_env, blacs_env=blacs_env)
3856 CALL cp_dbcsr_cholesky_invert(ri_2c_inv, para_env=para_env, blacs_env=blacs_env, upper_to_full=.true.)
3857 CASE (hfx_ri_do_2c_diag)
3858 CALL dbcsr_copy(ri_2c_inv, ri_2c_ints)
3859 CALL cp_dbcsr_power(ri_2c_inv, -1.0_dp, ri_data%eps_eigval, n_dependent, &
3860 para_env, blacs_env, verbose=ri_data%unit_nr_dbcsr > 0)
3861 END SELECT
3862
3863 CALL dbt_create(ri_2c_ints, t2c_ri_tmp)
3864 CALL create_2c_tensor(t2c_ri_ints, dist1, dist2, ri_data%pgrid_2d, &
3865 ri_data%bsizes_RI_split, ri_data%bsizes_RI_split, &
3866 name="(RI | RI)")
3867 CALL dbt_create(t2c_ri_ints, t2c_ri_inv)
3868
3869 CALL dbt_copy_matrix_to_tensor(ri_2c_ints, t2c_ri_tmp)
3870 CALL dbt_copy(t2c_ri_tmp, t2c_ri_ints, move_data=.true.)
3871 CALL dbt_filter(t2c_ri_ints, ri_data%filter_eps)
3872
3873 CALL dbt_copy_matrix_to_tensor(ri_2c_inv, t2c_ri_tmp)
3874 CALL dbt_copy(t2c_ri_tmp, t2c_ri_inv, move_data=.true.)
3875 CALL dbt_filter(t2c_ri_inv, ri_data%filter_eps)
3876
3877 CALL dbcsr_release(ri_2c_inv)
3878 CALL dbt_destroy(t2c_ri_tmp)
3879 DEALLOCATE (dist1, dist2)
3880
3881 ! The AO density tensor
3882 CALL dbt_create(rho_ao(ispin, 1)%matrix, rho_ao_tmp)
3883 CALL create_2c_tensor(rho_ao_t, dist1, dist2, ri_data%pgrid_2d, &
3884 ri_data%bsizes_AO_split, ri_data%bsizes_AO_split, &
3885 name="(AO | AO)")
3886 DEALLOCATE (dist1, dist2)
3887
3888 CALL dbt_copy_matrix_to_tensor(rho_ao(ispin, 1)%matrix, rho_ao_tmp)
3889 CALL dbt_copy(rho_ao_tmp, rho_ao_t, move_data=.true.)
3890 CALL dbt_filter(rho_ao_t, ri_data%filter_eps)
3891 CALL dbt_destroy(rho_ao_tmp)
3892
3893 ! Put in in 3D
3894 ALLOCATE (dist1(SIZE(ri_data%bsizes_AO_split)), dist2(SIZE(ri_data%bsizes_AO_split)), dist3(1))
3895 dist3(1) = 0
3896 CALL dbt_get_info(rho_ao_t, pdims=pdims_2d, proc_dist_1=dist1, proc_dist_2=dist2)
3897 CALL dbt_default_distvec(1, 1, [1], dist3)
3898
3899 pdims_3d(1) = pdims_2d(1)
3900 pdims_3d(2) = pdims_2d(2)
3901 pdims_3d(3) = 1
3902
3903 CALL dbt_pgrid_create(para_env, pdims_3d, pgrid_3d)
3904 CALL dbt_distribution_new(dist_3d, pgrid_3d, dist1, dist2, dist3)
3905 CALL dbt_create(rho_ao_t_3d, "rho_ao_3d", dist_3d, [1, 2], [3], ri_data%bsizes_AO_split, &
3906 ri_data%bsizes_AO_split, [1])
3907 DEALLOCATE (dist1, dist2, dist3)
3908 CALL dbt_pgrid_destroy(pgrid_3d)
3909 CALL dbt_distribution_destroy(dist_3d)
3910
3911 ! copy density
3912 nblks = 0
3913!$OMP PARALLEL DEFAULT(NONE) &
3914!$OMP SHARED(rho_ao_t,nblks) &
3915!$OMP PRIVATE(iter,ind,blk,found)
3916 CALL dbt_iterator_start(iter, rho_ao_t)
3917 DO WHILE (dbt_iterator_blocks_left(iter))
3918 CALL dbt_iterator_next_block(iter, ind)
3919 CALL dbt_get_block(rho_ao_t, ind, blk, found)
3920 IF (found) THEN
3921!$OMP ATOMIC
3922 nblks = nblks + 1
3923 DEALLOCATE (blk)
3924 END IF
3925 END DO
3926 CALL dbt_iterator_stop(iter)
3927!$OMP END PARALLEL
3928
3929 ALLOCATE (idx1(nblks), idx2(nblks), idx3(nblks))
3930 idx3 = 1
3931 nblks = 0
3932!$OMP PARALLEL DEFAULT(NONE) &
3933!$OMP SHARED(rho_ao_t,nblks,idx1,idx2) &
3934!$OMP PRIVATE(iter,ind,blk,found)
3935 CALL dbt_iterator_start(iter, rho_ao_t)
3936 DO WHILE (dbt_iterator_blocks_left(iter))
3937 CALL dbt_iterator_next_block(iter, ind)
3938 CALL dbt_get_block(rho_ao_t, ind, blk, found)
3939 IF (found) THEN
3940!$OMP CRITICAL
3941 nblks = nblks + 1
3942 idx1(nblks) = ind(1)
3943 idx2(nblks) = ind(2)
3944!$OMP END CRITICAL
3945 DEALLOCATE (blk)
3946 END IF
3947 END DO
3948 CALL dbt_iterator_stop(iter)
3949!$OMP END PARALLEL
3950
3951!$OMP PARALLEL DEFAULT(NONE) SHARED(rho_ao_t_3d,nblks,idx1,idx2,idx3) PRIVATE(nblk_per_thread,A,b)
3952 nblk_per_thread = nblks/omp_get_num_threads() + 1
3953 a = omp_get_thread_num()*nblk_per_thread + 1
3954 b = min(a + nblk_per_thread, nblks)
3955 CALL dbt_reserve_blocks(rho_ao_t_3d, idx1(a:b), idx2(a:b), idx3(a:b))
3956!$OMP END PARALLEL
3957
3958!$OMP PARALLEL DEFAULT(NONE) &
3959!$OMP SHARED(rho_ao_t,rho_ao_t_3d) &
3960!$OMP PRIVATE(iter,ind,blk,found,blk_3d)
3961 CALL dbt_iterator_start(iter, rho_ao_t)
3962 DO WHILE (dbt_iterator_blocks_left(iter))
3963 CALL dbt_iterator_next_block(iter, ind)
3964 CALL dbt_get_block(rho_ao_t, ind, blk, found)
3965 IF (found) THEN
3966 ALLOCATE (blk_3d(SIZE(blk, 1), SIZE(blk, 2), 1))
3967 blk_3d(:, :, 1) = blk(:, :)
3968!$OMP CRITICAL
3969 CALL dbt_put_block(rho_ao_t_3d, [ind(1), ind(2), 1], [SIZE(blk, 1), SIZE(blk, 2), 1], blk_3d)
3970!$OMP END CRITICAL
3971 DEALLOCATE (blk, blk_3d)
3972 END IF
3973 END DO
3974 CALL dbt_iterator_stop(iter)
3975!$OMP END PARALLEL
3976
3977 ! The 1D tensor with the density coeffs
3978 pdims_2d(1) = para_env%num_pe
3979 pdims_2d(2) = 1
3980
3981 ALLOCATE (dist1(SIZE(ri_data%bsizes_RI_split)), dist2(1))
3982 CALL dbt_default_distvec(SIZE(ri_data%bsizes_RI_split), pdims_2d(1), ri_data%bsizes_RI_split, dist1)
3983 CALL dbt_default_distvec(1, pdims_2d(2), [1], dist2)
3984
3985 CALL dbt_pgrid_create(para_env, pdims_2d, pgrid_2d)
3986 CALL dbt_distribution_new(dist_2d, pgrid_2d, dist1, dist2)
3987 CALL dbt_create(density_coeffs_t, "density_coeffs", dist_2d, [1], [2], ri_data%bsizes_RI_split, [1])
3988 CALL dbt_create(density_coeffs_t, density_tmp)
3989 DEALLOCATE (dist1, dist2)
3990 CALL dbt_pgrid_destroy(pgrid_2d)
3991 CALL dbt_distribution_destroy(dist_2d)
3992
3993 CALL dbt_get_info(ri_data%t_3c_int_ctr_3(1, 1), nfull_total=dims_3c)
3994
3995 ! The 3c integrals tensor, in case we compute them here
3996 pdims_3d = 0
3997 CALL dbt_pgrid_create(para_env, pdims_3d, pgrid_3d, tensor_dims=[max(1, natom/n_mem), natom, natom])
3998 ALLOCATE (t_3c_int_batched(1, 1))
3999 CALL create_3c_tensor(t_3c_int_batched(1, 1), dist1, dist2, dist3, pgrid_3d, &
4000 ri_data%bsizes_RI, ri_data%bsizes_AO, ri_data%bsizes_AO, map1=[1], map2=[2, 3], &
4001 name="(RI | AO AO)")
4002
4003 CALL dbt_mp_environ_pgrid(pgrid_3d, pdims_3d, pcoord_3d)
4004 CALL mp_comm_t3c%create(pgrid_3d%mp_comm_2d, 3, pdims_3d)
4005 CALL distribution_3d_create(dist_nl3c, dist1, dist2, dist3, nkind, particle_set, &
4006 mp_comm_t3c, own_comm=.true.)
4007 DEALLOCATE (dist1, dist2, dist3)
4008 CALL dbt_pgrid_destroy(pgrid_3d)
4009
4010 CALL build_3c_neighbor_lists(nl_3c, basis_set_ri, basis_set_ao, basis_set_ao, dist_nl3c, ri_data%ri_metric, &
4011 "HFX_3c_nl", qs_env, op_pos=1, sym_jk=.true., own_dist=.true.)
4012
4013 n_mem = ri_data%n_mem
4014 CALL create_tensor_batches(ri_data%bsizes_RI, n_mem, dummy1, dummy2, batch_block_start, batch_block_end)
4015
4016 calc_ints = .false.
4017 CALL get_tensor_occupancy(ri_data%t_3c_int_ctr_2(1, 1), nze, occ)
4018 IF (nze == 0) calc_ints = .true.
4019
4020 DO i_mem = 1, n_mem
4021 IF (calc_ints) THEN
4022 CALL build_3c_integrals(t_3c_int_batched, ri_data%filter_eps, qs_env, nl_3c, &
4023 basis_set_ri, basis_set_ao, basis_set_ao, &
4024 ri_data%ri_metric, int_eps=ri_data%eps_schwarz, op_pos=1, &
4025 desymmetrize=.false., &
4026 bounds_i=[batch_block_start(i_mem), batch_block_end(i_mem)])
4027 CALL dbt_copy(t_3c_int_batched(1, 1), ri_data%t_3c_int_ctr_3(1, 1), order=[1, 3, 2])
4028 CALL dbt_copy(t_3c_int_batched(1, 1), ri_data%t_3c_int_ctr_3(1, 1), move_data=.true., summation=.true.)
4029 CALL dbt_filter(ri_data%t_3c_int_ctr_3(1, 1), ri_data%filter_eps)
4030 ELSE
4031 bounds_cpy(:, 2) = [sum(ri_data%bsizes_RI(1:batch_block_start(i_mem) - 1)) + 1, &
4032 sum(ri_data%bsizes_RI(1:batch_block_end(i_mem)))]
4033 bounds_cpy(:, 1) = [1, sum(ri_data%bsizes_AO)]
4034 bounds_cpy(:, 3) = [1, sum(ri_data%bsizes_AO)]
4035 CALL dbt_copy(ri_data%t_3c_int_ctr_2(1, 1), ri_data%t_3c_int_ctr_3(1, 1), &
4036 order=[2, 1, 3], bounds=bounds_cpy)
4037 END IF
4038
4039 !contract the integrals with the density P_pq (pq|R)
4040 CALL dbt_contract(1.0_dp, ri_data%t_3c_int_ctr_3(1, 1), rho_ao_t_3d, 0.0_dp, density_tmp, &
4041 contract_1=[2, 3], notcontract_1=[1], &
4042 contract_2=[1, 2], notcontract_2=[3], &
4043 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps)
4044 CALL dbt_clear(ri_data%t_3c_int_ctr_3(1, 1))
4045
4046 !contract the above vector with the inverse metric
4047 CALL dbt_contract(1.0_dp, t2c_ri_inv, density_tmp, 1.0_dp, density_coeffs_t, &
4048 contract_1=[2], notcontract_1=[1], &
4049 contract_2=[1], notcontract_2=[2], &
4050 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps)
4051
4052 END DO
4053 CALL neighbor_list_3c_destroy(nl_3c)
4054
4055 IF (multiply_by_s) THEN
4056 CALL dbt_contract(1.0_dp, t2c_ri_ints, density_coeffs_t, 0.0_dp, density_tmp, &
4057 contract_1=[2], notcontract_1=[1], &
4058 contract_2=[1], notcontract_2=[2], &
4059 map_1=[1], map_2=[2], filter_eps=ri_data%filter_eps)
4060 CALL dbt_copy(density_tmp, density_coeffs_t, move_data=.true.)
4061 END IF
4062
4063 ALLOCATE (density_coeffs(sum(ri_data%bsizes_RI)))
4064 density_coeffs = 0.0
4065
4066!$OMP PARALLEL DEFAULT(NONE) &
4067!$OMP SHARED(density_coeffs_t,ri_data,density_coeffs) &
4068!$OMP PRIVATE(iter,ind,blk,found,idx)
4069 CALL dbt_iterator_start(iter, density_coeffs_t)
4070 DO WHILE (dbt_iterator_blocks_left(iter))
4071 CALL dbt_iterator_next_block(iter, ind)
4072 CALL dbt_get_block(density_coeffs_t, ind, blk, found)
4073 IF (found) THEN
4074
4075 idx = sum(ri_data%bsizes_RI_split(1:ind(1) - 1))
4076!$OMP CRITICAL
4077 density_coeffs(idx + 1:idx + ri_data%bsizes_RI_split(ind(1))) = blk(:, 1)
4078!$OMP END CRITICAL
4079 DEALLOCATE (blk)
4080 END IF
4081 END DO
4082 CALL dbt_iterator_stop(iter)
4083!$OMP END PARALLEL
4084 CALL para_env%sum(density_coeffs)
4085
4086 CALL dbt_destroy(t2c_ri_ints)
4087 CALL dbt_destroy(t2c_ri_inv)
4088 CALL dbt_destroy(density_tmp)
4089 CALL dbt_destroy(rho_ao_t)
4090 CALL dbt_destroy(rho_ao_t_3d)
4091 CALL dbt_destroy(density_coeffs_t)
4092 CALL dbt_destroy(t_3c_int_batched(1, 1))
4093
4094 CALL timestop(handle)
4095
4096 END SUBROUTINE get_ri_density_coeffs
4097
4098! **************************************************************************************************
4099!> \brief a small utility function that returns the atom corresponding to a block of a split tensor
4100!> \param idx_to_at ...
4101!> \param bsizes_split ...
4102!> \param bsizes_orig ...
4103!> \return ...
4104! **************************************************************************************************
4105 SUBROUTINE get_idx_to_atom(idx_to_at, bsizes_split, bsizes_orig)
4106 INTEGER, DIMENSION(:), INTENT(INOUT) :: idx_to_at
4107 INTEGER, DIMENSION(:), INTENT(IN) :: bsizes_split, bsizes_orig
4108
4109 INTEGER :: full_sum, iat, iblk, split_sum
4110
4111 iat = 1
4112 full_sum = bsizes_orig(iat)
4113 split_sum = 0
4114 DO iblk = 1, SIZE(bsizes_split)
4115 split_sum = split_sum + bsizes_split(iblk)
4116
4117 IF (split_sum .GT. full_sum) THEN
4118 iat = iat + 1
4119 full_sum = full_sum + bsizes_orig(iat)
4120 END IF
4121
4122 idx_to_at(iblk) = iat
4123 END DO
4124
4125 END SUBROUTINE get_idx_to_atom
4126
4127! **************************************************************************************************
4128!> \brief Function for calculating sqrt of a matrix
4129!> \param values ...
4130!> \return ...
4131! **************************************************************************************************
4132 FUNCTION my_sqrt(values)
4133 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: values
4134 REAL(kind=dp), DIMENSION(SIZE(values)) :: my_sqrt
4135
4136 my_sqrt = sqrt(values)
4137 END FUNCTION
4138
4139! **************************************************************************************************
4140!> \brief Function for calculation inverse sqrt of a matrix
4141!> \param values ...
4142!> \return ...
4143! **************************************************************************************************
4144 FUNCTION my_invsqrt(values)
4145 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: values
4146 REAL(kind=dp), DIMENSION(SIZE(values)) :: my_invsqrt
4147
4148 my_invsqrt = sqrt(1.0_dp/values)
4149 END FUNCTION
4150END MODULE
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:48
idx
static GRID_HOST_DEVICE int idx(const orbital a)
Return coset index of given orbital angular momentum.
Definition grid_common.h:153
idx2
#define idx2(a, i, j)
Definition grid_dgemm_tensor_local.h:258
idx3
#define idx3(a, i, j, k)
Definition grid_dgemm_tensor_local.h:257
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:89
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:15
arnoldi_api::arnoldi_extremal
subroutine, public arnoldi_extremal(matrix_a, max_ev, min_ev, converged, threshold, max_iter)
simple wrapper to estimate extremal eigenvalues with arnoldi, using the old lanczos interface this hi...
Definition arnoldi_api.F:254
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
basis_set_types
Definition basis_set_types.F:15
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_cholesky
Interface to (sca)lapack for the Cholesky based procedures.
Definition cp_dbcsr_cholesky.F:17
cp_dbcsr_cholesky::cp_dbcsr_cholesky_decompose
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
Definition cp_dbcsr_cholesky.F:61
cp_dbcsr_cholesky::cp_dbcsr_cholesky_invert
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, upper_to_full)
used to replace the cholesky decomposition by the inverse
Definition cp_dbcsr_cholesky.F:114
cp_dbcsr_diag
Interface to (sca)lapack for the Cholesky based procedures.
Definition cp_dbcsr_diag.F:17
cp_dbcsr_diag::cp_dbcsr_power
subroutine, public cp_dbcsr_power(matrix, exponent, threshold, n_dependent, para_env, blacs_env, verbose, eigenvectors, eigenvalues)
...
Definition cp_dbcsr_diag.F:233
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_dist2d_to_dist
subroutine, public cp_dbcsr_dist2d_to_dist(dist2d, dist)
Creates a DBCSR distribution from a distribution_2d.
Definition cp_dbcsr_operations.F:611
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:208
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:118
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_write_unformatted
subroutine, public cp_fm_write_unformatted(fm, unit)
...
Definition cp_fm_types.F:2147
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
distribution_2d_types
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
Definition distribution_2d_types.F:15
hfx_ri
RI-methods for HFX.
Definition hfx_ri.F:12
hfx_ri::hfx_ri_update_ks
subroutine, public hfx_ri_update_ks(qs_env, ri_data, ks_matrix, ehfx, mos, rho_ao, geometry_did_change, nspins, hf_fraction)
...
Definition hfx_ri.F:1033
hfx_ri::check_inverse
subroutine, public check_inverse(matrix_1, matrix_2, name, unit_nr)
...
Definition hfx_ri.F:347
hfx_ri::get_force_from_3c_trace
subroutine, public get_force_from_3c_trace(force, t_3c_contr, t_3c_der, atom_of_kind, kind_of, idx_to_at, pref, do_mp2, deriv_dim)
This routines calculates the force contribution from a trace over 3D tensors, i.e....
Definition hfx_ri.F:3360
hfx_ri::get_idx_to_atom
subroutine, public get_idx_to_atom(idx_to_at, bsizes_split, bsizes_orig)
a small utility function that returns the atom corresponding to a block of a split tensor
Definition hfx_ri.F:4106
hfx_ri::hfx_ri_pre_scf_calc_tensors
subroutine, public hfx_ri_pre_scf_calc_tensors(qs_env, ri_data, t_2c_int_ri, t_2c_int_pot, t_3c_int, do_kpoints)
Calculate 2-center and 3-center integrals.
Definition hfx_ri.F:441
hfx_ri::hfx_ri_update_forces
subroutine, public hfx_ri_update_forces(qs_env, ri_data, nspins, hf_fraction, rho_ao, rho_ao_resp, mos, use_virial, resp_only, rescale_factor)
the general routine that calls the relevant force code
Definition hfx_ri.F:3036
hfx_ri::print_ri_hfx
subroutine, public print_ri_hfx(ri_data, qs_env)
Print RI-HFX quantities, as required by the PRINT subsection.
Definition hfx_ri.F:3619
hfx_ri::get_2c_der_force
subroutine, public get_2c_der_force(force, t_2c_contr, t_2c_der, atom_of_kind, kind_of, idx_to_at, pref, do_mp2, do_ovlp)
Update the forces due to the derivative of the a 2-center product d/dR (Q|R)
Definition hfx_ri.F:3446
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:15
hfx_types::alloc_containers
subroutine, public alloc_containers(data, bin_size)
...
Definition hfx_types.F:2906
hfx_types::dealloc_containers
subroutine, public dealloc_containers(data, memory_usage)
...
Definition hfx_types.F:2874
hfx_types::hfx_ri_init
subroutine, public hfx_ri_init(ri_data, qs_kind_set, particle_set, atomic_kind_set, para_env)
...
Definition hfx_types.F:1199
hfx_types::hfx_ri_release
subroutine, public hfx_ri_release(ri_data, write_stats)
...
Definition hfx_types.F:1464
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::hfx_ri_do_2c_cholesky
integer, parameter, public hfx_ri_do_2c_cholesky
Definition input_constants.F:827
input_constants::hfx_ri_do_2c_diag
integer, parameter, public hfx_ri_do_2c_diag
Definition input_constants.F:827
input_constants::hfx_ri_do_2c_iter
integer, parameter, public hfx_ri_do_2c_iter
Definition input_constants.F:827
input_cp2k_hfx
function that builds the hartree fock exchange section of the input
Definition input_cp2k_hfx.F:14
input_cp2k_hfx::ri_pmat
integer, parameter, public ri_pmat
Definition input_cp2k_hfx.F:44
input_cp2k_hfx::ri_mo
integer, parameter, public ri_mo
Definition input_cp2k_hfx.F:44
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
iterate_matrix
Routines useful for iterative matrix calculations.
Definition iterate_matrix.F:13
iterate_matrix::matrix_sqrt_newton_schulz
subroutine, public matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged)
compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations the...
Definition iterate_matrix.F:1624
iterate_matrix::invert_hotelling
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
Definition iterate_matrix.F:472
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:123
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_force_types
Definition qs_force_types.F:13
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition qs_integral_utils.F:149
qs_interactions
Calculate the interaction radii for the operator matrix calculation.
Definition qs_interactions.F:17
qs_interactions::init_interaction_radii_orb_basis
subroutine, public init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_short)
...
Definition qs_interactions.F:420
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_ks_types
Definition qs_ks_types.F:15
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:397
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::release_neighbor_list_sets
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Definition qs_neighbor_list_types.F:1097
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_tensors_types
Utility methods to build 3-center integral tensors of various types.
Definition qs_tensors_types.F:12
qs_tensors_types::distribution_3d_create
subroutine, public distribution_3d_create(dist_3d, dist1, dist2, dist3, nkind, particle_set, mp_comm_3d, own_comm)
Create a 3d distribution.
Definition qs_tensors_types.F:91
qs_tensors_types::create_2c_tensor
subroutine, public create_2c_tensor(t2c, dist_1, dist_2, pgrid, sizes_1, sizes_2, order, name)
...
Definition qs_tensors_types.F:382
qs_tensors_types::split_block_sizes
subroutine, public split_block_sizes(blk_sizes, blk_sizes_split, max_size)
...
Definition qs_tensors_types.F:428
qs_tensors_types::create_tensor_batches
subroutine, public create_tensor_batches(sizes, nbatches, starts_array, ends_array, starts_array_block, ends_array_block)
...
Definition qs_tensors_types.F:522
qs_tensors_types::create_3c_tensor
subroutine, public create_3c_tensor(t3c, dist_1, dist_2, dist_3, pgrid, sizes_1, sizes_2, sizes_3, map1, map2, name)
...
Definition qs_tensors_types.F:335
qs_tensors
Utility methods to build 3-center integral tensors of various types.
Definition qs_tensors.F:11
qs_tensors::build_2c_integrals
subroutine, public build_2c_integrals(t2c, filter_eps, qs_env, nl_2c, basis_i, basis_j, potential_parameter, do_kpoints, do_hfx_kpoints, ext_kpoints, regularization_ri)
...
Definition qs_tensors.F:3406
qs_tensors::calc_3c_virial
subroutine, public calc_3c_virial(work_virial, t3c_trace, pref, qs_env, nl_3c, basis_i, basis_j, basis_k, potential_parameter, der_eps, op_pos)
Calculates the 3c virial contributions on the fly.
Definition qs_tensors.F:1868
qs_tensors::build_3c_derivatives
subroutine, public build_3c_derivatives(t3c_der_i, t3c_der_k, filter_eps, qs_env, nl_3c, basis_i, basis_j, basis_k, potential_parameter, der_eps, op_pos, do_kpoints, do_hfx_kpoints, bounds_i, bounds_j, bounds_k, ri_range, img_to_ri_cell)
Build 3-center derivative tensors.
Definition qs_tensors.F:1147
qs_tensors::build_2c_neighbor_lists
subroutine, public build_2c_neighbor_lists(ij_list, basis_i, basis_j, potential_parameter, name, qs_env, sym_ij, molecular, dist_2d, pot_to_rad)
Build 2-center neighborlists adapted to different operators This mainly wraps build_neighbor_lists fo...
Definition qs_tensors.F:143
qs_tensors::compress_tensor
subroutine, public compress_tensor(tensor, blk_indices, compressed, eps, memory)
...
Definition qs_tensors.F:3660
qs_tensors::neighbor_list_3c_destroy
subroutine, public neighbor_list_3c_destroy(ijk_list)
Destroy 3c neighborlist.
Definition qs_tensors.F:383
qs_tensors::calc_2c_virial
subroutine, public calc_2c_virial(work_virial, t2c_trace, pref, qs_env, nl_2c, basis_i, basis_j, potential_parameter)
Calculates the virial coming from 2c derivatives on the fly.
Definition qs_tensors.F:3226
qs_tensors::decompress_tensor
subroutine, public decompress_tensor(tensor, blk_indices, compressed, eps)
...
Definition qs_tensors.F:3788
qs_tensors::build_2c_derivatives
subroutine, public build_2c_derivatives(t2c_der, filter_eps, qs_env, nl_2c, basis_i, basis_j, potential_parameter, do_kpoints)
Calculates the derivatives of 2-center integrals, wrt to the first center.
Definition qs_tensors.F:2983
qs_tensors::get_tensor_occupancy
subroutine, public get_tensor_occupancy(tensor, nze, occ)
...
Definition qs_tensors.F:3879
qs_tensors::build_3c_neighbor_lists
subroutine, public build_3c_neighbor_lists(ijk_list, basis_i, basis_j, basis_k, dist_3d, potential_parameter, name, qs_env, sym_ij, sym_jk, sym_ik, molecular, op_pos, own_dist)
Build a 3-center neighbor list.
Definition qs_tensors.F:282
qs_tensors::build_3c_integrals
subroutine, public build_3c_integrals(t3c, filter_eps, qs_env, nl_3c, basis_i, basis_j, basis_k, potential_parameter, int_eps, op_pos, do_kpoints, do_hfx_kpoints, desymmetrize, bounds_i, bounds_j, bounds_k, ri_range, img_to_ri_cell)
Build 3-center integral tensor.
Definition qs_tensors.F:2376
util
All kind of helpful little routines.
Definition util.F:14
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:93
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
distribution_2d_types::distribution_2d_type
distributes pairs on a 2d grid of processors
Definition distribution_2d_types.F:62
hfx_types::block_ind_type
Definition hfx_types.F:367
hfx_types::hfx_compression_type
Definition hfx_types.F:351
hfx_types::hfx_ri_type
Definition hfx_types.F:371
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_cart_type
Definition message_passing.F:740
message_passing::mp_comm_type
Definition message_passing.F:189
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:135
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_tensors_types::distribution_3d_type
Definition qs_tensors_types.F:57
qs_tensors_types::neighbor_list_3c_type
Definition qs_tensors_types.F:63
virial_types::virial_type
Definition virial_types.F:26