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