(git:0de0cc2)
rpa_im_time_force_methods.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Routines needed for cubic-scaling RPA and SOS-Laplace-MP2 forces
10 !> \author Augustin Bussy
11 ! **************************************************************************************************
14  USE admm_types, ONLY: admm_type,&
17  USE atomic_kind_types, ONLY: atomic_kind_type,&
19  USE basis_set_types, ONLY: gto_basis_set_p_type,&
20  gto_basis_set_type
21  USE bibliography, ONLY: bussy2023,&
22  cite_reference
23  USE cell_types, ONLY: cell_type,&
24  pbc
25  USE core_ae, ONLY: build_core_ae
26  USE core_ppl, ONLY: build_core_ppl
27  USE core_ppnl, ONLY: build_core_ppnl
28  USE cp_blacs_env, ONLY: cp_blacs_env_type
29  USE cp_control_types, ONLY: dft_control_type
32  USE cp_dbcsr_diag, ONLY: cp_dbcsr_power
42  cp_fm_struct_type
43  USE cp_fm_types, ONLY: cp_fm_create,&
44  cp_fm_release,&
46  cp_fm_to_fm,&
47  cp_fm_type
48  USE dbcsr_api, ONLY: &
49  dbcsr_add, dbcsr_add_on_diag, dbcsr_clear, dbcsr_complete_redistribute, dbcsr_copy, &
50  dbcsr_create, dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, &
51  dbcsr_frobenius_norm, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
52  dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
53  dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, &
54  dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
55  USE dbt_api, ONLY: &
56  dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_clear, dbt_contract, &
57  dbt_copy, dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, dbt_destroy, &
58  dbt_filter, dbt_get_info, dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, &
59  dbt_pgrid_type, dbt_scale, dbt_type
60  USE distribution_2d_types, ONLY: distribution_2d_type
61  USE ec_methods, ONLY: create_kernel
65  USE hfx_exx, ONLY: add_exx_to_rhs
66  USE hfx_ri, ONLY: get_2c_der_force,&
70  USE hfx_types, ONLY: alloc_containers,&
71  block_ind_type,&
73  hfx_compression_type,&
74  hfx_type
76  do_eri_gpw,&
77  do_eri_mme,&
82  section_vals_type,&
85  USE kinds, ONLY: dp,&
86  int_8
87  USE libint_2c_3c, ONLY: libint_potential_type
88  USE machine, ONLY: m_flush,&
90  USE mathconstants, ONLY: fourpi
91  USE message_passing, ONLY: mp_cart_type,&
93  mp_para_env_type
94  USE mp2_eri, ONLY: integrate_set_2c
95  USE mp2_eri_gpw, ONLY: calc_potential_gpw,&
96  cleanup_gpw,&
97  prepare_gpw,&
99  USE mp2_types, ONLY: mp2_type
100  USE orbital_pointers, ONLY: ncoset
101  USE parallel_gemm_api, ONLY: parallel_gemm
103  USE particle_types, ONLY: particle_type
104  USE pw_env_types, ONLY: pw_env_get,&
105  pw_env_type
106  USE pw_methods, ONLY: pw_axpy,&
107  pw_copy,&
108  pw_integral_ab,&
109  pw_scale,&
110  pw_transfer,&
111  pw_zero
112  USE pw_poisson_methods, ONLY: pw_poisson_solve
113  USE pw_poisson_types, ONLY: pw_poisson_type
114  USE pw_pool_types, ONLY: pw_pool_type
115  USE pw_types, ONLY: pw_c1d_gs_type,&
116  pw_r3d_rs_type
121  USE qs_environment_types, ONLY: get_qs_env,&
122  qs_environment_type,&
123  set_qs_env
124  USE qs_force_types, ONLY: qs_force_type
126  USE qs_integrate_potential, ONLY: integrate_pgf_product,&
127  integrate_v_core_rspace,&
128  integrate_v_rspace
130  USE qs_kind_types, ONLY: qs_kind_type
134  USE qs_ks_types, ONLY: set_ks_env
135  USE qs_linres_types, ONLY: linres_control_type
136  USE qs_mo_types, ONLY: get_mo_set,&
137  mo_set_type
138  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
141  USE qs_p_env_methods, ONLY: p_env_create,&
143  USE qs_p_env_types, ONLY: p_env_release,&
144  qs_p_env_type
145  USE qs_rho_types, ONLY: qs_rho_get,&
146  qs_rho_type
147  USE qs_tensors, ONLY: &
155  distribution_3d_type,&
156  neighbor_list_3c_type
158  realspace_grid_type
161  USE rpa_im_time_force_types, ONLY: im_time_force_type
163  USE task_list_types, ONLY: task_list_type
164  USE virial_types, ONLY: virial_type
165 #include "./base/base_uses.f90"
166 
167  IMPLICIT NONE
168 
169  PRIVATE
170 
171  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_im_time_force_methods'
172 
175 
176 CONTAINS
177 
178 ! **************************************************************************************************
179 !> \brief Initializes and pre-calculates all needed tensors for the forces
180 !> \param force_data ...
181 !> \param fm_matrix_PQ ...
182 !> \param t_3c_M the 3-center M tensor to be used as a template
183 !> \param unit_nr ...
184 !> \param mp2_env ...
185 !> \param qs_env ...
186 ! **************************************************************************************************
187  SUBROUTINE init_im_time_forces(force_data, fm_matrix_PQ, t_3c_M, unit_nr, mp2_env, qs_env)
188 
189  TYPE(im_time_force_type), INTENT(INOUT) :: force_data
190  TYPE(cp_fm_type), INTENT(IN) :: fm_matrix_pq
191  TYPE(dbt_type), INTENT(INOUT) :: t_3c_m
192  INTEGER, INTENT(IN) :: unit_nr
193  TYPE(mp2_type) :: mp2_env
194  TYPE(qs_environment_type), POINTER :: qs_env
195 
196  CHARACTER(LEN=*), PARAMETER :: routinen = 'init_im_time_forces'
197 
198  INTEGER :: handle, i_mem, i_xyz, ibasis, ispin, &
199  n_dependent, n_mem, n_rep, natom, &
200  nkind, nspins
201  INTEGER(int_8) :: nze, nze_tot
202  INTEGER, ALLOCATABLE, DIMENSION(:) :: dist1, dist2, dist_ao_1, dist_ao_2, &
203  dist_ri, dummy_end, dummy_start, &
204  end_blocks, sizes_ao, sizes_ri, &
205  start_blocks
206  INTEGER, DIMENSION(2) :: pdims_t2c
207  INTEGER, DIMENSION(3) :: nblks_total, pcoord, pdims, pdims_t3c
208  INTEGER, DIMENSION(:), POINTER :: col_bsize, row_bsize
209  LOGICAL :: do_periodic, use_virial
210  REAL(dp) :: compression_factor, eps_pgf_orb, &
211  eps_pgf_orb_old, memory, occ
212  TYPE(cell_type), POINTER :: cell
213  TYPE(cp_blacs_env_type), POINTER :: blacs_env
214  TYPE(dbcsr_distribution_type) :: dbcsr_dist
215  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
216  TYPE(dbcsr_type) :: dbcsr_work, dbcsr_work2, dbcsr_work3
217  TYPE(dbcsr_type), DIMENSION(1) :: t_2c_int_tmp
218  TYPE(dbcsr_type), DIMENSION(1, 3) :: t_2c_der_tmp
219  TYPE(dbt_pgrid_type) :: pgrid_t2c, pgrid_t3c
220  TYPE(dbt_type) :: t_2c_template, t_2c_tmp, t_3c_template
221  TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :, :) :: t_3c_der_ao_prv, t_3c_der_ri_prv
222  TYPE(dft_control_type), POINTER :: dft_control
223  TYPE(distribution_2d_type), POINTER :: dist_2d
224  TYPE(distribution_3d_type) :: dist_3d, dist_vir
225  TYPE(gto_basis_set_p_type), ALLOCATABLE, &
226  DIMENSION(:), TARGET :: basis_set_ao, basis_set_ri_aux
227  TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
228  TYPE(libint_potential_type) :: identity_pot
229  TYPE(mp_cart_type) :: mp_comm_t3c, mp_comm_vir
230  TYPE(mp_para_env_type), POINTER :: para_env
231  TYPE(neighbor_list_3c_type) :: nl_3c
232  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
233  POINTER :: nl_2c
234  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
235  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
236  TYPE(qs_rho_type), POINTER :: rho
237  TYPE(section_vals_type), POINTER :: qs_section
238  TYPE(virial_type), POINTER :: virial
239 
240  NULLIFY (dft_control, para_env, particle_set, qs_kind_set, dist_2d, nl_2c, blacs_env, matrix_s, &
241  rho, rho_ao, cell, qs_section, orb_basis, ri_basis, virial)
242 
243  CALL cite_reference(bussy2023)
244 
245  CALL timeset(routinen, handle)
246 
247  CALL get_qs_env(qs_env, natom=natom, nkind=nkind, dft_control=dft_control, para_env=para_env, &
248  particle_set=particle_set, qs_kind_set=qs_kind_set, cell=cell, virial=virial)
249  IF (dft_control%qs_control%gapw) THEN
250  cpabort("Low-scaling RPA/SOS-MP2 forces only available with GPW")
251  END IF
252 
253  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
254 
255  do_periodic = .false.
256  IF (any(cell%perd == 1)) do_periodic = .true.
257  force_data%do_periodic = do_periodic
258 
259  !Dealing with the 3-center derivatives
260  pdims_t3c = 0
261  CALL dbt_pgrid_create(para_env, pdims_t3c, pgrid_t3c)
262 
263  !Make sure we use the proper QS EPS_PGF_ORB values
264  qs_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS")
265  CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", n_rep_val=n_rep)
266  IF (n_rep /= 0) THEN
267  CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", r_val=eps_pgf_orb)
268  ELSE
269  CALL section_vals_val_get(qs_section, "EPS_DEFAULT", r_val=eps_pgf_orb)
270  eps_pgf_orb = sqrt(eps_pgf_orb)
271  END IF
272  eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
273 
274  ALLOCATE (sizes_ri(natom), sizes_ao(natom))
275  ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind))
276  CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
277  CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ri, basis=basis_set_ri_aux)
278  CALL basis_set_list_setup(basis_set_ao, "ORB", qs_kind_set)
279  CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ao, basis=basis_set_ao)
280 
281  DO ibasis = 1, SIZE(basis_set_ao)
282  orb_basis => basis_set_ao(ibasis)%gto_basis_set
283  CALL init_interaction_radii_orb_basis(orb_basis, eps_pgf_orb)
284  ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
285  CALL init_interaction_radii_orb_basis(ri_basis, eps_pgf_orb)
286  END DO
287 
288  CALL create_3c_tensor(t_3c_template, dist_ri, dist_ao_1, dist_ao_2, pgrid_t3c, &
289  sizes_ri, sizes_ao, sizes_ao, map1=[1], map2=[2, 3], name="der (RI AO | AO)")
290 
291  ALLOCATE (t_3c_der_ri_prv(1, 1, 3), t_3c_der_ao_prv(1, 1, 3))
292  DO i_xyz = 1, 3
293  CALL dbt_create(t_3c_template, t_3c_der_ri_prv(1, 1, i_xyz))
294  CALL dbt_create(t_3c_template, t_3c_der_ao_prv(1, 1, i_xyz))
295  END DO
296 
297  IF (use_virial) THEN
298  ALLOCATE (force_data%t_3c_virial, force_data%t_3c_virial_split)
299  CALL dbt_create(t_3c_template, force_data%t_3c_virial)
300  CALL dbt_create(t_3c_m, force_data%t_3c_virial_split)
301  END IF
302  CALL dbt_destroy(t_3c_template)
303 
304  CALL dbt_mp_environ_pgrid(pgrid_t3c, pdims, pcoord)
305  CALL mp_comm_t3c%create(pgrid_t3c%mp_comm_2d, 3, pdims)
306  CALL distribution_3d_create(dist_3d, dist_ri, dist_ao_1, dist_ao_2, &
307  nkind, particle_set, mp_comm_t3c, own_comm=.true.)
308 
309  !In case of virial, we need to store the 3c_nl
310  IF (use_virial) THEN
311  ALLOCATE (force_data%nl_3c)
312  CALL mp_comm_vir%create(pgrid_t3c%mp_comm_2d, 3, pdims)
313  CALL distribution_3d_create(dist_vir, dist_ri, dist_ao_1, dist_ao_2, &
314  nkind, particle_set, mp_comm_vir, own_comm=.true.)
315  CALL build_3c_neighbor_lists(force_data%nl_3c, basis_set_ri_aux, basis_set_ao, basis_set_ao, &
316  dist_vir, mp2_env%ri_metric, "RPA_3c_nl", qs_env, op_pos=1, &
317  sym_jk=.false., own_dist=.true.)
318  END IF
319 
320  CALL build_3c_neighbor_lists(nl_3c, basis_set_ri_aux, basis_set_ao, basis_set_ao, dist_3d, &
321  mp2_env%ri_metric, "RPA_3c_nl", qs_env, op_pos=1, sym_jk=.true., &
322  own_dist=.true.)
323  DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
324 
325  !Prepare the resulting 3c tensors in the format of t_3c_M for compatible traces: (RI|AO AO), split blocks
326  CALL dbt_get_info(t_3c_m, nblks_total=nblks_total)
327  ALLOCATE (force_data%bsizes_RI_split(nblks_total(1)), force_data%bsizes_AO_split(nblks_total(2)))
328  CALL dbt_get_info(t_3c_m, blk_size_1=force_data%bsizes_RI_split, blk_size_2=force_data%bsizes_AO_split)
329  DO i_xyz = 1, 3
330  CALL dbt_create(t_3c_m, force_data%t_3c_der_RI(i_xyz))
331  CALL dbt_create(t_3c_m, force_data%t_3c_der_AO(i_xyz))
332  END DO
333 
334  !Keep track of atom index corresponding to split blocks
335  ALLOCATE (force_data%idx_to_at_RI(nblks_total(1)))
336  CALL get_idx_to_atom(force_data%idx_to_at_RI, force_data%bsizes_RI_split, sizes_ri)
337 
338  ALLOCATE (force_data%idx_to_at_AO(nblks_total(2)))
339  CALL get_idx_to_atom(force_data%idx_to_at_AO, force_data%bsizes_AO_split, sizes_ao)
340 
341  n_mem = mp2_env%ri_rpa_im_time%cut_memory
342  CALL create_tensor_batches(sizes_ri, n_mem, dummy_start, dummy_end, start_blocks, end_blocks)
343  DEALLOCATE (dummy_start, dummy_end)
344 
345  ALLOCATE (force_data%t_3c_der_AO_comp(n_mem, 3), force_data%t_3c_der_RI_comp(n_mem, 3))
346  ALLOCATE (force_data%t_3c_der_AO_ind(n_mem, 3), force_data%t_3c_der_RI_ind(n_mem, 3))
347 
348  memory = 0.0_dp
349  nze_tot = 0
350  DO i_mem = 1, n_mem
351  CALL build_3c_derivatives(t_3c_der_ri_prv, t_3c_der_ao_prv, mp2_env%ri_rpa_im_time%eps_filter, &
352  qs_env, nl_3c, basis_set_ri_aux, basis_set_ao, basis_set_ao, &
353  mp2_env%ri_metric, der_eps=mp2_env%ri_rpa_im_time%eps_filter, op_pos=1, &
354  bounds_i=[start_blocks(i_mem), end_blocks(i_mem)])
355 
356  DO i_xyz = 1, 3
357  CALL dbt_copy(t_3c_der_ri_prv(1, 1, i_xyz), force_data%t_3c_der_RI(i_xyz), move_data=.true.)
358  CALL dbt_filter(force_data%t_3c_der_RI(i_xyz), mp2_env%ri_rpa_im_time%eps_filter)
359  CALL get_tensor_occupancy(force_data%t_3c_der_RI(i_xyz), nze, occ)
360  nze_tot = nze_tot + nze
361 
362  CALL alloc_containers(force_data%t_3c_der_RI_comp(i_mem, i_xyz), 1)
363  CALL compress_tensor(force_data%t_3c_der_RI(i_xyz), force_data%t_3c_der_RI_ind(i_mem, i_xyz)%ind, &
364  force_data%t_3c_der_RI_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress, memory)
365  CALL dbt_clear(force_data%t_3c_der_RI(i_xyz))
366 
367  CALL dbt_copy(t_3c_der_ao_prv(1, 1, i_xyz), force_data%t_3c_der_AO(i_xyz), move_data=.true.)
368  CALL dbt_filter(force_data%t_3c_der_AO(i_xyz), mp2_env%ri_rpa_im_time%eps_filter)
369  CALL get_tensor_occupancy(force_data%t_3c_der_AO(i_xyz), nze, occ)
370  nze_tot = nze_tot + nze
371 
372  CALL alloc_containers(force_data%t_3c_der_AO_comp(i_mem, i_xyz), 1)
373  CALL compress_tensor(force_data%t_3c_der_AO(i_xyz), force_data%t_3c_der_AO_ind(i_mem, i_xyz)%ind, &
374  force_data%t_3c_der_AO_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress, memory)
375  CALL dbt_clear(force_data%t_3c_der_AO(i_xyz))
376  END DO
377  END DO
378  CALL neighbor_list_3c_destroy(nl_3c)
379  DO i_xyz = 1, 3
380  CALL dbt_destroy(t_3c_der_ri_prv(1, 1, i_xyz))
381  CALL dbt_destroy(t_3c_der_ao_prv(1, 1, i_xyz))
382  END DO
383 
384  CALL para_env%sum(memory)
385  compression_factor = real(nze_tot, dp)*1.0e-06*8.0_dp/memory
386  IF (unit_nr > 0) THEN
387  WRITE (unit=unit_nr, fmt="((T3,A,T66,F11.2,A4))") &
388  "MEMORY_INFO| Memory for 3-center derivatives (compressed):", memory, ' MiB'
389 
390  WRITE (unit=unit_nr, fmt="((T3,A,T60,F21.2))") &
391  "MEMORY_INFO| Compression factor: ", compression_factor
392  END IF
393 
394  !Dealing with the 2-center derivatives
395  CALL get_qs_env(qs_env, distribution_2d=dist_2d, blacs_env=blacs_env, matrix_s=matrix_s)
396  CALL cp_dbcsr_dist2d_to_dist(dist_2d, dbcsr_dist)
397  ALLOCATE (row_bsize(SIZE(sizes_ri)))
398  ALLOCATE (col_bsize(SIZE(sizes_ri)))
399  row_bsize(:) = sizes_ri(:)
400  col_bsize(:) = sizes_ri(:)
401 
402  pdims_t2c = 0
403  CALL dbt_pgrid_create(para_env, pdims_t2c, pgrid_t2c)
404  CALL create_2c_tensor(t_2c_template, dist1, dist2, pgrid_t2c, force_data%bsizes_RI_split, &
405  force_data%bsizes_RI_split, name='(RI| RI)')
406  DEALLOCATE (dist1, dist2)
407 
408  CALL dbcsr_create(t_2c_int_tmp(1), "(P|Q) RPA", dbcsr_dist, dbcsr_type_symmetric, row_bsize, col_bsize)
409  DO i_xyz = 1, 3
410  CALL dbcsr_create(t_2c_der_tmp(1, i_xyz), "(P|Q) RPA der", dbcsr_dist, &
411  dbcsr_type_antisymmetric, row_bsize, col_bsize)
412  END DO
413 
414  IF (use_virial) THEN
415  ALLOCATE (force_data%RI_virial_pot, force_data%RI_virial_met)
416  CALL dbcsr_create(force_data%RI_virial_pot, "RI_virial", dbcsr_dist, &
417  dbcsr_type_no_symmetry, row_bsize, col_bsize)
418  CALL dbcsr_create(force_data%RI_virial_met, "RI_virial", dbcsr_dist, &
419  dbcsr_type_no_symmetry, row_bsize, col_bsize)
420  END IF
421 
422  ! Main (P|Q) integrals and derivatives
423  ! Integrals are passed as a full matrix => convert to DBCSR
424  CALL dbcsr_create(dbcsr_work, template=t_2c_int_tmp(1))
425  CALL copy_fm_to_dbcsr(fm_matrix_pq, dbcsr_work)
426 
427  ! We need the +/- square root of (P|Q)
428  CALL dbcsr_create(dbcsr_work2, template=t_2c_int_tmp(1))
429  CALL dbcsr_create(dbcsr_work3, template=t_2c_int_tmp(1))
430  CALL dbcsr_copy(dbcsr_work2, dbcsr_work)
431  CALL cp_dbcsr_power(dbcsr_work, -0.5_dp, 1.0e-7_dp, n_dependent, para_env, blacs_env) !1.0E-7 ev qunenching thresh
432 
433  ! Transfer to tensor format with split blocks
434  CALL dbt_create(dbcsr_work, t_2c_tmp)
435  CALL dbt_copy_matrix_to_tensor(dbcsr_work, t_2c_tmp)
436  CALL dbt_create(t_2c_template, force_data%t_2c_pot_msqrt)
437  CALL dbt_copy(t_2c_tmp, force_data%t_2c_pot_msqrt, move_data=.true.)
438  CALL dbt_filter(force_data%t_2c_pot_msqrt, mp2_env%ri_rpa_im_time%eps_filter)
439 
440  CALL dbcsr_multiply('N', 'N', 1.0_dp, dbcsr_work2, dbcsr_work, 0.0_dp, dbcsr_work3)
441  CALL dbt_copy_matrix_to_tensor(dbcsr_work3, t_2c_tmp)
442  CALL dbt_create(t_2c_template, force_data%t_2c_pot_psqrt)
443  CALL dbt_copy(t_2c_tmp, force_data%t_2c_pot_psqrt, move_data=.true.)
444  CALL dbt_filter(force_data%t_2c_pot_psqrt, mp2_env%ri_rpa_im_time%eps_filter)
445  CALL dbt_destroy(t_2c_tmp)
446  CALL dbcsr_release(dbcsr_work2)
447  CALL dbcsr_release(dbcsr_work3)
448  CALL dbcsr_clear(dbcsr_work)
449 
450  ! Deal with the 2c potential derivatives. Only precompute if not in PBCs
451  IF (.NOT. do_periodic) THEN
452  CALL build_2c_neighbor_lists(nl_2c, basis_set_ri_aux, basis_set_ri_aux, mp2_env%potential_parameter, &
453  "RPA_2c_nl_pot", qs_env, sym_ij=.true., dist_2d=dist_2d)
454  CALL build_2c_derivatives(t_2c_der_tmp, mp2_env%ri_rpa_im_time%eps_filter, qs_env, nl_2c, &
455  basis_set_ri_aux, basis_set_ri_aux, mp2_env%potential_parameter)
456  CALL release_neighbor_list_sets(nl_2c)
457 
458  DO i_xyz = 1, 3
459  CALL dbt_create(t_2c_der_tmp(1, i_xyz), t_2c_tmp)
460  CALL dbt_copy_matrix_to_tensor(t_2c_der_tmp(1, i_xyz), t_2c_tmp)
461  CALL dbt_create(t_2c_template, force_data%t_2c_der_pot(i_xyz))
462  CALL dbt_copy(t_2c_tmp, force_data%t_2c_der_pot(i_xyz), move_data=.true.)
463  CALL dbt_filter(force_data%t_2c_der_pot(i_xyz), mp2_env%ri_rpa_im_time%eps_filter)
464  CALL dbt_destroy(t_2c_tmp)
465  CALL dbcsr_clear(t_2c_der_tmp(1, i_xyz))
466  END DO
467 
468  IF (use_virial) THEN
469  CALL build_2c_neighbor_lists(force_data%nl_2c_pot, basis_set_ri_aux, basis_set_ri_aux, &
470  mp2_env%potential_parameter, "RPA_2c_nl_pot", qs_env, &
471  sym_ij=.false., dist_2d=dist_2d)
472  END IF
473  END IF
474  ! Create a G_PQ matrix to collect the terms for the force trace in the periodic case
475  CALL dbcsr_create(force_data%G_PQ, "G_PQ", dbcsr_dist, dbcsr_type_no_symmetry, row_bsize, col_bsize)
476 
477  ! we need the RI metric derivatives and the inverse of the integrals
478  CALL build_2c_neighbor_lists(nl_2c, basis_set_ri_aux, basis_set_ri_aux, mp2_env%ri_metric, &
479  "RPA_2c_nl_metric", qs_env, sym_ij=.true., dist_2d=dist_2d)
480  CALL build_2c_integrals(t_2c_int_tmp, mp2_env%ri_rpa_im_time%eps_filter, qs_env, nl_2c, &
481  basis_set_ri_aux, basis_set_ri_aux, mp2_env%ri_metric)
482  CALL build_2c_derivatives(t_2c_der_tmp, mp2_env%ri_rpa_im_time%eps_filter, qs_env, nl_2c, &
483  basis_set_ri_aux, basis_set_ri_aux, mp2_env%ri_metric)
484  CALL release_neighbor_list_sets(nl_2c)
485 
486  IF (use_virial) THEN
487  CALL build_2c_neighbor_lists(force_data%nl_2c_met, basis_set_ri_aux, basis_set_ri_aux, &
488  mp2_env%ri_metric, "RPA_2c_nl_metric", qs_env, sym_ij=.false., &
489  dist_2d=dist_2d)
490  END IF
491 
492  CALL dbcsr_copy(dbcsr_work, t_2c_int_tmp(1))
493  CALL cp_dbcsr_cholesky_decompose(dbcsr_work, para_env=para_env, blacs_env=blacs_env)
494  CALL cp_dbcsr_cholesky_invert(dbcsr_work, para_env=para_env, blacs_env=blacs_env, upper_to_full=.true.)
495 
496  CALL dbt_create(dbcsr_work, t_2c_tmp)
497  CALL dbt_copy_matrix_to_tensor(dbcsr_work, t_2c_tmp)
498  CALL dbt_create(t_2c_template, force_data%t_2c_inv_metric)
499  CALL dbt_copy(t_2c_tmp, force_data%t_2c_inv_metric, move_data=.true.)
500  CALL dbt_filter(force_data%t_2c_inv_metric, mp2_env%ri_rpa_im_time%eps_filter)
501  CALL dbt_destroy(t_2c_tmp)
502  CALL dbcsr_clear(dbcsr_work)
503  CALL dbcsr_clear(t_2c_int_tmp(1))
504 
505  DO i_xyz = 1, 3
506  CALL dbt_create(t_2c_der_tmp(1, i_xyz), t_2c_tmp)
507  CALL dbt_copy_matrix_to_tensor(t_2c_der_tmp(1, i_xyz), t_2c_tmp)
508  CALL dbt_create(t_2c_template, force_data%t_2c_der_metric(i_xyz))
509  CALL dbt_copy(t_2c_tmp, force_data%t_2c_der_metric(i_xyz), move_data=.true.)
510  CALL dbt_filter(force_data%t_2c_der_metric(i_xyz), mp2_env%ri_rpa_im_time%eps_filter)
511  CALL dbt_destroy(t_2c_tmp)
512  CALL dbcsr_clear(t_2c_der_tmp(1, i_xyz))
513  END DO
514 
515  !Pre-calculate matrix K = metric^-1 * V^0.5
516  CALL dbt_create(t_2c_template, force_data%t_2c_K)
517  CALL dbt_contract(1.0_dp, force_data%t_2c_inv_metric, force_data%t_2c_pot_psqrt, &
518  0.0_dp, force_data%t_2c_K, &
519  contract_1=[2], notcontract_1=[1], &
520  contract_2=[1], notcontract_2=[2], &
521  map_1=[1], map_2=[2], filter_eps=mp2_env%ri_rpa_im_time%eps_filter)
522 
523  ! Finally, we need the overlap matrix derivative and the inverse of the integrals
524  CALL dbt_destroy(t_2c_template)
525  CALL dbcsr_release(dbcsr_work)
526  CALL dbcsr_release(t_2c_int_tmp(1))
527  DO i_xyz = 1, 3
528  CALL dbcsr_release(t_2c_der_tmp(1, i_xyz))
529  END DO
530 
531  DEALLOCATE (row_bsize, col_bsize)
532  ALLOCATE (row_bsize(SIZE(sizes_ao)))
533  ALLOCATE (col_bsize(SIZE(sizes_ao)))
534  row_bsize(:) = sizes_ao(:)
535  col_bsize(:) = sizes_ao(:)
536 
537  CALL create_2c_tensor(t_2c_template, dist1, dist2, pgrid_t2c, force_data%bsizes_AO_split, &
538  force_data%bsizes_AO_split, name='(AO| AO)')
539  DEALLOCATE (dist1, dist2)
540 
541  DO i_xyz = 1, 3
542  CALL dbcsr_create(t_2c_der_tmp(1, i_xyz), "(P|Q) RPA der", dbcsr_dist, &
543  dbcsr_type_antisymmetric, row_bsize, col_bsize)
544  END DO
545 
546  identity_pot%potential_type = do_potential_id
547  CALL build_2c_neighbor_lists(nl_2c, basis_set_ao, basis_set_ao, identity_pot, &
548  "RPA_2c_nl_metric", qs_env, sym_ij=.true., dist_2d=dist_2d)
549  CALL build_2c_derivatives(t_2c_der_tmp, mp2_env%ri_rpa_im_time%eps_filter, qs_env, nl_2c, &
550  basis_set_ao, basis_set_ao, identity_pot)
551  CALL release_neighbor_list_sets(nl_2c)
552 
553  IF (use_virial) THEN
554  CALL build_2c_neighbor_lists(force_data%nl_2c_ovlp, basis_set_ao, basis_set_ao, identity_pot, &
555  "RPA_2c_nl_metric", qs_env, sym_ij=.false., dist_2d=dist_2d)
556  END IF
557 
558  CALL dbcsr_create(force_data%inv_ovlp, template=matrix_s(1)%matrix)
559  CALL dbcsr_copy(force_data%inv_ovlp, matrix_s(1)%matrix)
560  CALL cp_dbcsr_cholesky_decompose(force_data%inv_ovlp, para_env=para_env, blacs_env=blacs_env)
561  CALL cp_dbcsr_cholesky_invert(force_data%inv_ovlp, para_env=para_env, blacs_env=blacs_env, upper_to_full=.true.)
562 
563  DO i_xyz = 1, 3
564  CALL dbt_create(t_2c_der_tmp(1, i_xyz), t_2c_tmp)
565  CALL dbt_copy_matrix_to_tensor(t_2c_der_tmp(1, i_xyz), t_2c_tmp)
566  CALL dbt_create(t_2c_template, force_data%t_2c_der_ovlp(i_xyz))
567  CALL dbt_copy(t_2c_tmp, force_data%t_2c_der_ovlp(i_xyz), move_data=.true.)
568  CALL dbt_filter(force_data%t_2c_der_ovlp(i_xyz), mp2_env%ri_rpa_im_time%eps_filter)
569  CALL dbt_destroy(t_2c_tmp)
570  CALL dbcsr_clear(t_2c_der_tmp(1, i_xyz))
571  END DO
572 
573  !Create the rest of the 2-center AO tensors
574  nspins = dft_control%nspins
575  ALLOCATE (force_data%P_virt(nspins), force_data%P_occ(nspins))
576  ALLOCATE (force_data%sum_YP_tau(nspins), force_data%sum_O_tau(nspins))
577  DO ispin = 1, nspins
578  ALLOCATE (force_data%P_virt(ispin)%matrix, force_data%P_occ(ispin)%matrix)
579  ALLOCATE (force_data%sum_YP_tau(ispin)%matrix, force_data%sum_O_tau(ispin)%matrix)
580  CALL dbcsr_create(force_data%P_virt(ispin)%matrix, template=matrix_s(1)%matrix)
581  CALL dbcsr_create(force_data%P_occ(ispin)%matrix, template=matrix_s(1)%matrix)
582  CALL dbcsr_create(force_data%sum_O_tau(ispin)%matrix, template=matrix_s(1)%matrix)
583  CALL dbcsr_create(force_data%sum_YP_tau(ispin)%matrix, template=matrix_s(1)%matrix)
584 
585  CALL dbcsr_copy(force_data%sum_O_tau(ispin)%matrix, matrix_s(1)%matrix)
586  CALL dbcsr_copy(force_data%sum_YP_tau(ispin)%matrix, matrix_s(1)%matrix)
587 
588  CALL dbcsr_set(force_data%sum_O_tau(ispin)%matrix, 0.0_dp)
589  CALL dbcsr_set(force_data%sum_YP_tau(ispin)%matrix, 0.0_dp)
590  END DO
591 
592  !Populate the density matrices: 1 = P_virt*S +P_occ*S ==> P_virt = S^-1 - P_occ
593  CALL get_qs_env(qs_env, rho=rho)
594  CALL qs_rho_get(rho, rho_ao=rho_ao)
595  CALL dbcsr_copy(force_data%P_occ(1)%matrix, rho_ao(1)%matrix)
596  IF (nspins == 1) THEN
597  CALL dbcsr_scale(force_data%P_occ(1)%matrix, 0.5_dp) !because double occupency
598  ELSE
599  CALL dbcsr_copy(force_data%P_occ(2)%matrix, rho_ao(2)%matrix)
600  END IF
601  DO ispin = 1, nspins
602  CALL dbcsr_copy(force_data%P_virt(ispin)%matrix, force_data%inv_ovlp)
603  CALL dbcsr_add(force_data%P_virt(ispin)%matrix, force_data%P_occ(ispin)%matrix, 1.0_dp, -1.0_dp)
604  END DO
605 
606  DO ibasis = 1, SIZE(basis_set_ao)
607  orb_basis => basis_set_ao(ibasis)%gto_basis_set
608  CALL init_interaction_radii_orb_basis(orb_basis, eps_pgf_orb_old)
609  ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
610  CALL init_interaction_radii_orb_basis(ri_basis, eps_pgf_orb_old)
611  END DO
612 
613  CALL dbt_destroy(t_2c_template)
614  CALL dbcsr_release(dbcsr_work)
615  DO i_xyz = 1, 3
616  CALL dbcsr_release(t_2c_der_tmp(1, i_xyz))
617  END DO
618  DEALLOCATE (row_bsize, col_bsize)
619  CALL dbt_pgrid_destroy(pgrid_t3c)
620  CALL dbt_pgrid_destroy(pgrid_t2c)
621  CALL dbcsr_distribution_release(dbcsr_dist)
622  CALL timestop(handle)
623 
624  END SUBROUTINE init_im_time_forces
625 
626 ! **************************************************************************************************
627 !> \brief Updates the cubic-scaling SOS-Laplace-MP2 contribution to the forces at each quadrature point
628 !> \param force_data ...
629 !> \param mat_P_omega ...
630 !> \param t_3c_M ...
631 !> \param t_3c_O ...
632 !> \param t_3c_O_compressed ...
633 !> \param t_3c_O_ind ...
634 !> \param fm_scaled_dm_occ_tau ...
635 !> \param fm_scaled_dm_virt_tau ...
636 !> \param fm_mo_coeff_occ ...
637 !> \param fm_mo_coeff_virt ...
638 !> \param fm_mo_coeff_occ_scaled ...
639 !> \param fm_mo_coeff_virt_scaled ...
640 !> \param starts_array_mc ...
641 !> \param ends_array_mc ...
642 !> \param starts_array_mc_block ...
643 !> \param ends_array_mc_block ...
644 !> \param num_integ_points ...
645 !> \param nmo ...
646 !> \param Eigenval ...
647 !> \param tau_tj ...
648 !> \param tau_wj ...
649 !> \param cut_memory ...
650 !> \param Pspin ...
651 !> \param Qspin ...
652 !> \param open_shell ...
653 !> \param unit_nr ...
654 !> \param dbcsr_time ...
655 !> \param dbcsr_nflop ...
656 !> \param mp2_env ...
657 !> \param qs_env ...
658 !> \note In open-shell, we need to take Q from one spin, and everything from the other
659 ! **************************************************************************************************
660  SUBROUTINE calc_laplace_loop_forces(force_data, mat_P_omega, t_3c_M, t_3c_O, t_3c_O_compressed, &
661  t_3c_O_ind, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, &
662  fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, &
663  fm_mo_coeff_virt_scaled, starts_array_mc, ends_array_mc, &
664  starts_array_mc_block, ends_array_mc_block, num_integ_points, &
665  nmo, Eigenval, tau_tj, tau_wj, cut_memory, Pspin, Qspin, &
666  open_shell, unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
667 
668  TYPE(im_time_force_type), INTENT(INOUT) :: force_data
669  TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_p_omega
670  TYPE(dbt_type), INTENT(INOUT) :: t_3c_m, t_3c_o
671  TYPE(hfx_compression_type), DIMENSION(:) :: t_3c_o_compressed
672  TYPE(block_ind_type), DIMENSION(:), INTENT(INOUT) :: t_3c_o_ind
673  TYPE(cp_fm_type), INTENT(IN) :: fm_scaled_dm_occ_tau, &
674  fm_scaled_dm_virt_tau
675  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mo_coeff_occ, fm_mo_coeff_virt
676  TYPE(cp_fm_type), INTENT(IN) :: fm_mo_coeff_occ_scaled, &
677  fm_mo_coeff_virt_scaled
678  INTEGER, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
679  starts_array_mc_block, &
680  ends_array_mc_block
681  INTEGER, INTENT(IN) :: num_integ_points, nmo
682  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: eigenval
683  REAL(kind=dp), DIMENSION(0:num_integ_points), &
684  INTENT(IN) :: tau_tj
685  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
686  INTENT(IN) :: tau_wj
687  INTEGER, INTENT(IN) :: cut_memory, pspin, qspin
688  LOGICAL, INTENT(IN) :: open_shell
689  INTEGER, INTENT(IN) :: unit_nr
690  REAL(dp), INTENT(INOUT) :: dbcsr_time
691  INTEGER(int_8), INTENT(INOUT) :: dbcsr_nflop
692  TYPE(mp2_type) :: mp2_env
693  TYPE(qs_environment_type), POINTER :: qs_env
694 
695  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_laplace_loop_forces'
696 
697  INTEGER :: dummy_int, handle, handle2, i_mem, i_xyz, ibasis, ispin, j_xyz, jquad, k_xyz, &
698  n_mem_ri, n_rep, natom, nkind, nspins, unit_nr_dbcsr
699  INTEGER(int_8) :: flop, nze, nze_ddint, nze_der_ao, &
700  nze_der_ri, nze_kqk
701  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, batch_blk_end, &
702  batch_blk_start, batch_end_ri, &
703  batch_start_ri, kind_of, mc_ranges, &
704  mc_ranges_ri
705  INTEGER, DIMENSION(:, :), POINTER :: dummy_ptr
706  LOGICAL :: memory_info, use_virial
707  REAL(dp) :: eps_filter, eps_pgf_orb, &
708  eps_pgf_orb_old, fac, occ, occ_ddint, &
709  occ_der_ao, occ_der_ri, occ_kqk, &
710  omega, pref, t1, t2, tau
711  REAL(dp), DIMENSION(3, 3) :: work_virial, work_virial_ovlp
712  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
713  TYPE(cell_type), POINTER :: cell
714  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
715  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_occ, mat_dm_virt
716  TYPE(dbcsr_type) :: dbcsr_work1, dbcsr_work2, dbcsr_work3, &
717  exp_occ, exp_virt, r_occ, r_virt, &
718  virial_ovlp, y_1, y_2
719  TYPE(dbt_type) :: t_2c_ao, t_2c_ri, t_2c_ri_2, t_2c_tmp, t_3c_0, t_3c_1, t_3c_3, t_3c_4, &
720  t_3c_5, t_3c_6, t_3c_7, t_3c_8, t_3c_help_1, t_3c_help_2, t_3c_ints, t_3c_sparse, &
721  t_3c_work, t_dm_occ, t_dm_virt, t_kqkt, t_m_occ, t_m_virt, t_q, t_r_occ, t_r_virt
722  TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_p
723  TYPE(dft_control_type), POINTER :: dft_control
724  TYPE(gto_basis_set_p_type), ALLOCATABLE, &
725  DIMENSION(:), TARGET :: basis_set_ao, basis_set_ri_aux
726  TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
727  TYPE(libint_potential_type) :: identity_pot
728  TYPE(mp_para_env_type), POINTER :: para_env
729  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
730  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
731  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
732  TYPE(section_vals_type), POINTER :: qs_section
733  TYPE(virial_type), POINTER :: virial
734 
735  NULLIFY (matrix_s, dummy_ptr, atomic_kind_set, force, matrix_s, matrix_ks, mat_dm_occ, mat_dm_virt)
736  NULLIFY (dft_control, virial, particle_set, cell, para_env, orb_basis, ri_basis, qs_section)
737  NULLIFY (qs_kind_set)
738 
739  CALL timeset(routinen, handle)
740 
741  CALL get_qs_env(qs_env, matrix_s=matrix_s, natom=natom, atomic_kind_set=atomic_kind_set, &
742  force=force, matrix_ks=matrix_ks, dft_control=dft_control, virial=virial, &
743  particle_set=particle_set, cell=cell, para_env=para_env, nkind=nkind, &
744  qs_kind_set=qs_kind_set)
745  eps_filter = mp2_env%ri_rpa_im_time%eps_filter
746  nspins = dft_control%nspins
747 
748  memory_info = mp2_env%ri_rpa_im_time%memory_info
749  IF (memory_info) THEN
750  unit_nr_dbcsr = unit_nr
751  ELSE
752  unit_nr_dbcsr = 0
753  END IF
754 
755  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
756 
757  IF (use_virial) virial%pv_calculate = .true.
758 
759  IF (use_virial) THEN
760  qs_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS")
761  CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", n_rep_val=n_rep)
762  IF (n_rep /= 0) THEN
763  CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", r_val=eps_pgf_orb)
764  ELSE
765  CALL section_vals_val_get(qs_section, "EPS_DEFAULT", r_val=eps_pgf_orb)
766  eps_pgf_orb = sqrt(eps_pgf_orb)
767  END IF
768  eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
769 
770  ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind))
771  CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
772  CALL basis_set_list_setup(basis_set_ao, "ORB", qs_kind_set)
773 
774  DO ibasis = 1, SIZE(basis_set_ao)
775  orb_basis => basis_set_ao(ibasis)%gto_basis_set
776  CALL init_interaction_radii_orb_basis(orb_basis, eps_pgf_orb)
777  ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
778  CALL init_interaction_radii_orb_basis(ri_basis, eps_pgf_orb)
779  END DO
780  END IF
781 
782  !We follow the general logic of the compute_mat_P_omega routine
783  ALLOCATE (t_p(nspins))
784  CALL dbt_create(force_data%t_2c_K, t_2c_ri)
785  CALL dbt_create(force_data%t_2c_K, t_2c_ri_2)
786  CALL dbt_create(force_data%t_2c_der_ovlp(1), t_2c_ao)
787 
788  CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
789 
790  ! Always do the batching of the MO on mu and sigma, such that it is consistent between
791  ! the occupied and the virtual quantities
792  ALLOCATE (mc_ranges(cut_memory + 1))
793  mc_ranges(:cut_memory) = starts_array_mc_block(:)
794  mc_ranges(cut_memory + 1) = ends_array_mc_block(cut_memory) + 1
795 
796  ! Also need some batching on the RI, because it loses sparsity at some point
797  n_mem_ri = cut_memory
798  CALL create_tensor_batches(force_data%bsizes_RI_split, n_mem_ri, batch_start_ri, batch_end_ri, &
799  batch_blk_start, batch_blk_end)
800  ALLOCATE (mc_ranges_ri(n_mem_ri + 1))
801  mc_ranges_ri(1:n_mem_ri) = batch_blk_start(1:n_mem_ri)
802  mc_ranges_ri(n_mem_ri + 1) = batch_blk_end(n_mem_ri) + 1
803  DEALLOCATE (batch_blk_start, batch_blk_end)
804 
805  !Pre-allocate all required tensors and matrices
806  DO ispin = 1, nspins
807  CALL dbt_create(t_2c_ri, t_p(ispin))
808  END DO
809  CALL dbt_create(t_2c_ri, t_q)
810  CALL dbt_create(t_2c_ri, t_kqkt)
811  CALL dbt_create(t_2c_ao, t_dm_occ)
812  CALL dbt_create(t_2c_ao, t_dm_virt)
813 
814  !note: t_3c_O and t_3c_M have different mappings (map_1d, map_2d)
815  CALL dbt_create(t_3c_o, t_m_occ)
816  CALL dbt_create(t_3c_o, t_m_virt)
817  CALL dbt_create(t_3c_o, t_3c_0)
818 
819  CALL dbt_create(t_3c_o, t_3c_1)
820  CALL dbt_create(t_3c_o, t_3c_3)
821  CALL dbt_create(t_3c_o, t_3c_4)
822  CALL dbt_create(t_3c_o, t_3c_5)
823  CALL dbt_create(t_3c_m, t_3c_6)
824  CALL dbt_create(t_3c_m, t_3c_7)
825  CALL dbt_create(t_3c_m, t_3c_8)
826  CALL dbt_create(t_3c_m, t_3c_sparse)
827  CALL dbt_create(t_3c_o, t_3c_help_1)
828  CALL dbt_create(t_3c_o, t_3c_help_2)
829  CALL dbt_create(t_2c_ao, t_r_occ)
830  CALL dbt_create(t_2c_ao, t_r_virt)
831  CALL dbt_create(t_3c_m, t_3c_ints)
832  CALL dbt_create(t_3c_m, t_3c_work)
833 
834  !Pre-define the sparsity of t_3c_4 as a function of the derivatives
835  occ_der_ao = 0; nze_der_ao = 0
836  occ_der_ri = 0; nze_der_ri = 0
837  DO i_xyz = 1, 3
838  DO i_mem = 1, cut_memory
839  CALL decompress_tensor(force_data%t_3c_der_RI(i_xyz), force_data%t_3c_der_RI_ind(i_mem, i_xyz)%ind, &
840  force_data%t_3c_der_RI_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
841  CALL get_tensor_occupancy(force_data%t_3c_der_RI(i_xyz), nze, occ)
842  occ_der_ri = occ_der_ri + occ
843  nze_der_ri = nze_der_ri + nze
844  CALL dbt_copy(force_data%t_3c_der_RI(i_xyz), t_3c_sparse, summation=.true., move_data=.true.)
845 
846  CALL decompress_tensor(force_data%t_3c_der_AO(i_xyz), force_data%t_3c_der_AO_ind(i_mem, i_xyz)%ind, &
847  force_data%t_3c_der_AO_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
848  CALL get_tensor_occupancy(force_data%t_3c_der_AO(i_xyz), nze, occ)
849  occ_der_ao = occ_der_ao + occ
850  nze_der_ao = nze_der_ao + nze
851  CALL dbt_copy(force_data%t_3c_der_AO(i_xyz), t_3c_sparse, order=[1, 3, 2], summation=.true.)
852  CALL dbt_copy(force_data%t_3c_der_AO(i_xyz), t_3c_sparse, summation=.true., move_data=.true.)
853  END DO
854  END DO
855  occ_der_ri = occ_der_ri/3.0_dp
856  occ_der_ao = occ_der_ao/3.0_dp
857  nze_der_ri = nze_der_ri/3
858  nze_der_ao = nze_der_ao/3
859 
860  CALL dbcsr_create(r_occ, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
861  CALL dbcsr_create(r_virt, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
862  CALL dbcsr_create(dbcsr_work1, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
863  CALL dbcsr_create(dbcsr_work2, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
864  CALL dbcsr_create(dbcsr_work3, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
865  CALL dbcsr_create(exp_occ, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
866  CALL dbcsr_create(exp_virt, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
867  IF (use_virial) CALL dbcsr_create(virial_ovlp, template=dbcsr_work1)
868 
869  CALL dbt_batched_contract_init(t_3c_0, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
870  CALL dbt_batched_contract_init(t_3c_1, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
871  CALL dbt_batched_contract_init(t_3c_3, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
872  CALL dbt_batched_contract_init(t_m_occ, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
873  CALL dbt_batched_contract_init(t_m_virt, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
874 
875  CALL dbt_batched_contract_init(t_3c_ints, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
876  CALL dbt_batched_contract_init(t_3c_work, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
877 
878  CALL dbt_batched_contract_init(t_3c_4, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
879  batch_range_3=mc_ranges)
880  CALL dbt_batched_contract_init(t_3c_5, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
881  batch_range_3=mc_ranges)
882  CALL dbt_batched_contract_init(t_3c_6, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
883  batch_range_3=mc_ranges)
884  CALL dbt_batched_contract_init(t_3c_7, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
885  batch_range_3=mc_ranges)
886  CALL dbt_batched_contract_init(t_3c_8, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
887  batch_range_3=mc_ranges)
888  CALL dbt_batched_contract_init(t_3c_sparse, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
889  batch_range_3=mc_ranges)
890 
891  work_virial = 0.0_dp
892  work_virial_ovlp = 0.0_dp
893  DO jquad = 1, num_integ_points
894  tau = tau_tj(jquad)
895  omega = tau_wj(jquad)
896  fac = -2.0_dp*omega*mp2_env%scale_S
897  IF (open_shell) fac = 0.5_dp*fac
898  occ_ddint = 0; nze_ddint = 0
899 
900  CALL para_env%sync()
901  t1 = m_walltime()
902 
903  !Deal with the force contributions where there is no explicit 3-center quantities, i.e. the
904  !forces due to the metric and potential derivatives
905  DO ispin = 1, nspins
906  CALL dbt_create(mat_p_omega(jquad, ispin)%matrix, t_2c_tmp)
907  CALL dbt_copy_matrix_to_tensor(mat_p_omega(jquad, ispin)%matrix, t_2c_tmp)
908  CALL dbt_copy(t_2c_tmp, t_p(ispin), move_data=.true.)
909  CALL dbt_filter(t_p(ispin), eps_filter)
910  CALL dbt_destroy(t_2c_tmp)
911  END DO
912 
913  !Q = K^T*P*K, open-shell: Q is from one spin, everything else from the other
914  CALL dbt_contract(1.0_dp, t_p(qspin), force_data%t_2c_K, 0.0_dp, t_2c_ri, &
915  contract_1=[2], notcontract_1=[1], &
916  contract_2=[1], notcontract_2=[2], &
917  map_1=[1], map_2=[2], filter_eps=eps_filter, &
918  flop=flop, unit_nr=unit_nr_dbcsr)
919  dbcsr_nflop = dbcsr_nflop + flop
920  CALL dbt_contract(1.0_dp, force_data%t_2c_K, t_2c_ri, 0.0_dp, t_q, &
921  contract_1=[1], notcontract_1=[2], &
922  contract_2=[1], notcontract_2=[2], &
923  map_1=[1], map_2=[2], filter_eps=eps_filter, &
924  flop=flop, unit_nr=unit_nr_dbcsr)
925  dbcsr_nflop = dbcsr_nflop + flop
926  CALL dbt_clear(t_2c_ri)
927 
928  CALL perform_2c_ops(force, t_kqkt, force_data, fac, t_q, t_p(pspin), t_2c_ri, t_2c_ri_2, &
929  use_virial, atom_of_kind, kind_of, eps_filter, dbcsr_nflop, unit_nr_dbcsr)
930  CALL get_tensor_occupancy(t_kqkt, nze_kqk, occ_kqk)
931 
932  !Calculate the pseudo-density matrix in tensor form. There are a few useless arguments for SOS-MP2
933  CALL compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, tau_tj, num_integ_points, &
934  nmo, fm_mo_coeff_occ(pspin), fm_mo_coeff_virt(pspin), &
935  fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm_occ, mat_dm_virt, &
936  matrix_s, pspin, eigenval(:, pspin), 0.0_dp, eps_filter, &
937  mp2_env%ri_rpa_im_time%memory_info, unit_nr, &
938  jquad, .false., .false., qs_env, dummy_int, dummy_ptr, para_env)
939 
940  CALL dbt_create(mat_dm_occ(jquad, 1)%matrix, t_2c_tmp)
941  CALL dbt_copy_matrix_to_tensor(mat_dm_occ(jquad, 1)%matrix, t_2c_tmp)
942  CALL dbt_copy(t_2c_tmp, t_dm_occ, move_data=.true.)
943  CALL dbt_filter(t_dm_occ, eps_filter)
944  CALL dbt_destroy(t_2c_tmp)
945 
946  CALL dbt_create(mat_dm_virt(jquad, 1)%matrix, t_2c_tmp)
947  CALL dbt_copy_matrix_to_tensor(mat_dm_virt(jquad, 1)%matrix, t_2c_tmp)
948  CALL dbt_copy(t_2c_tmp, t_dm_virt, move_data=.true.)
949  CALL dbt_filter(t_dm_virt, eps_filter)
950  CALL dbt_destroy(t_2c_tmp)
951 
952  !Deal with the 3-center quantities.
953  CALL perform_3c_ops(force, t_r_occ, t_r_virt, force_data, fac, cut_memory, n_mem_ri, &
954  t_kqkt, t_dm_occ, t_dm_virt, t_3c_o, t_3c_m, t_m_occ, t_m_virt, t_3c_0, t_3c_1, &
955  t_3c_3, t_3c_4, t_3c_5, t_3c_6, t_3c_7, t_3c_8, t_3c_sparse, t_3c_help_1, t_3c_help_2, &
956  t_3c_ints, t_3c_work, starts_array_mc, ends_array_mc, batch_start_ri, &
957  batch_end_ri, t_3c_o_compressed, t_3c_o_ind, use_virial, &
958  atom_of_kind, kind_of, eps_filter, occ_ddint, nze_ddint, dbcsr_nflop, &
959  unit_nr_dbcsr, mp2_env)
960 
961  CALL timeset(routinen//"_dbcsr", handle2)
962  !We go back to DBCSR matrices from now on
963  !Note: R matrices are in fact symmetric, but use a normal type for convenience
964  CALL dbt_create(matrix_s(1)%matrix, t_2c_tmp)
965  CALL dbt_copy(t_r_occ, t_2c_tmp, move_data=.true.)
966  CALL dbt_copy_tensor_to_matrix(t_2c_tmp, r_occ)
967 
968  CALL dbt_copy(t_r_virt, t_2c_tmp, move_data=.true.)
969  CALL dbt_copy_tensor_to_matrix(t_2c_tmp, r_virt)
970 
971  !Iteratively calculate the Y1 and Y2 matrices
972  CALL dbcsr_multiply('N', 'N', tau, force_data%P_occ(pspin)%matrix, &
973  matrix_ks(pspin)%matrix, 0.0_dp, dbcsr_work1)
974  CALL build_y_matrix(y_1, dbcsr_work1, force_data%P_occ(pspin)%matrix, r_virt, eps_filter)
975  CALL matrix_exponential(exp_occ, dbcsr_work1, 1.0_dp, 1.0_dp, eps_filter)
976 
977  CALL dbcsr_multiply('N', 'N', -tau, force_data%P_virt(pspin)%matrix, &
978  matrix_ks(pspin)%matrix, 0.0_dp, dbcsr_work1)
979  CALL build_y_matrix(y_2, dbcsr_work1, force_data%P_virt(pspin)%matrix, r_occ, eps_filter)
980  CALL matrix_exponential(exp_virt, dbcsr_work1, 1.0_dp, 1.0_dp, eps_filter)
981 
982  !The force contribution coming from [-S^-1*(e^-tau*P_virt*F)^T*R_occ*S^-1
983  ! +tau*S^-1*Y_2^T*F*S^-1] * der_S
984  CALL dbcsr_multiply('N', 'N', 1.0_dp, r_occ, force_data%inv_ovlp, 0.0_dp, dbcsr_work1)
985  CALL dbcsr_multiply('T', 'N', 1.0_dp, exp_virt, dbcsr_work1, 0.0_dp, dbcsr_work3)
986  CALL dbcsr_multiply('N', 'N', 1.0_dp, force_data%inv_ovlp, dbcsr_work3, 0.0_dp, dbcsr_work2)
987 
988  CALL dbcsr_multiply('N', 'T', tau, force_data%inv_ovlp, y_2, 0.0_dp, dbcsr_work3)
989  CALL dbcsr_multiply('N', 'N', 1.0_dp, dbcsr_work3, matrix_ks(pspin)%matrix, 0.0_dp, dbcsr_work1)
990  CALL dbcsr_multiply('N', 'N', 1.0_dp, dbcsr_work1, force_data%inv_ovlp, 0.0_dp, dbcsr_work3)
991 
992  CALL dbcsr_add(dbcsr_work2, dbcsr_work3, 1.0_dp, -1.0_dp)
993 
994  CALL dbt_copy_matrix_to_tensor(dbcsr_work2, t_2c_tmp)
995  CALL dbt_copy(t_2c_tmp, t_2c_ao, move_data=.true.)
996 
997  pref = -1.0_dp*fac
998  CALL get_2c_der_force(force, t_2c_ao, force_data%t_2c_der_ovlp, atom_of_kind, &
999  kind_of, force_data%idx_to_at_AO, pref, do_ovlp=.true.)
1000 
1001  IF (use_virial) CALL dbcsr_add(virial_ovlp, dbcsr_work2, 1.0_dp, pref)
1002 
1003  !The final contribution from Tr[(tau*Y_1*P_occ - tau*Y_2*P_virt) * der_F]
1004  CALL dbcsr_multiply('N', 'N', tau*fac, y_1, force_data%P_occ(pspin)%matrix, 1.0_dp, &
1005  force_data%sum_YP_tau(pspin)%matrix, retain_sparsity=.true.)
1006  CALL dbcsr_multiply('N', 'N', -tau*fac, y_2, force_data%P_virt(pspin)%matrix, 1.0_dp, &
1007  force_data%sum_YP_tau(pspin)%matrix, retain_sparsity=.true.)
1008 
1009  !Build-up the RHS of the response equation.
1010  pref = -omega*mp2_env%scale_S
1011  CALL dbcsr_multiply('N', 'N', pref, r_virt, exp_occ, 1.0_dp, &
1012  force_data%sum_O_tau(pspin)%matrix, retain_sparsity=.true.)
1013  CALL dbcsr_multiply('N', 'N', -pref, r_occ, exp_virt, 1.0_dp, &
1014  force_data%sum_O_tau(pspin)%matrix, retain_sparsity=.true.)
1015  CALL dbcsr_multiply('N', 'N', pref*tau, matrix_ks(pspin)%matrix, y_1, 1.0_dp, &
1016  force_data%sum_O_tau(pspin)%matrix, retain_sparsity=.true.)
1017  CALL dbcsr_multiply('N', 'N', pref*tau, matrix_ks(pspin)%matrix, y_2, 1.0_dp, &
1018  force_data%sum_O_tau(pspin)%matrix, retain_sparsity=.true.)
1019 
1020  CALL timestop(handle2)
1021 
1022  !Print some info
1023  CALL para_env%sync()
1024  t2 = m_walltime()
1025  dbcsr_time = dbcsr_time + t2 - t1
1026 
1027  IF (unit_nr > 0) THEN
1028  WRITE (unit_nr, '(/T3,A,1X,I3,A)') &
1029  'RPA_LOW_SCALING_INFO| Info for time point', jquad, ' (gradients)'
1030  WRITE (unit_nr, '(T6,A,T56,F25.6)') &
1031  'Execution time (s):', t2 - t1
1032  WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1033  'Occupancy of 3c AO derivs:', real(nze_der_ao, dp), '/', occ_der_ao*100, '%'
1034  WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1035  'Occupancy of 3c RI derivs:', real(nze_der_ri, dp), '/', occ_der_ri*100, '%'
1036  WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1037  'Occupancy of the Docc * Dvirt * 3c-int tensor', real(nze_ddint, dp), '/', occ_ddint*100, '%'
1038  WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1039  'Occupancy of KQK^T 2c-tensor:', real(nze_kqk, dp), '/', occ_kqk*100, '%'
1040  CALL m_flush(unit_nr)
1041  END IF
1042 
1043  !intermediate clean-up
1044  CALL dbcsr_release(y_1)
1045  CALL dbcsr_release(y_2)
1046  CALL dbt_destroy(t_2c_tmp)
1047  END DO !jquad
1048 
1049  CALL dbt_batched_contract_finalize(t_3c_0)
1050  CALL dbt_batched_contract_finalize(t_3c_1)
1051  CALL dbt_batched_contract_finalize(t_3c_3)
1052  CALL dbt_batched_contract_finalize(t_m_occ)
1053  CALL dbt_batched_contract_finalize(t_m_virt)
1054 
1055  CALL dbt_batched_contract_finalize(t_3c_ints)
1056  CALL dbt_batched_contract_finalize(t_3c_work)
1057 
1058  CALL dbt_batched_contract_finalize(t_3c_4)
1059  CALL dbt_batched_contract_finalize(t_3c_5)
1060  CALL dbt_batched_contract_finalize(t_3c_6)
1061  CALL dbt_batched_contract_finalize(t_3c_7)
1062  CALL dbt_batched_contract_finalize(t_3c_8)
1063  CALL dbt_batched_contract_finalize(t_3c_sparse)
1064 
1065  !Calculate the 2c and 3c contributions to the virial
1066  IF (use_virial) THEN
1067  CALL dbt_copy(force_data%t_3c_virial_split, force_data%t_3c_virial, move_data=.true.)
1068  CALL calc_3c_virial(work_virial, force_data%t_3c_virial, 1.0_dp, qs_env, force_data%nl_3c, &
1069  basis_set_ri_aux, basis_set_ao, basis_set_ao, mp2_env%ri_metric, &
1070  der_eps=mp2_env%ri_rpa_im_time%eps_filter, op_pos=1)
1071 
1072  CALL calc_2c_virial(work_virial, force_data%RI_virial_met, 1.0_dp, qs_env, force_data%nl_2c_met, &
1073  basis_set_ri_aux, basis_set_ri_aux, mp2_env%ri_metric)
1074  CALL dbcsr_clear(force_data%RI_virial_met)
1075 
1076  IF (.NOT. force_data%do_periodic) THEN
1077  CALL calc_2c_virial(work_virial, force_data%RI_virial_pot, 1.0_dp, qs_env, force_data%nl_2c_pot, &
1078  basis_set_ri_aux, basis_set_ri_aux, mp2_env%potential_parameter)
1079  CALL dbcsr_clear(force_data%RI_virial_pot)
1080  END IF
1081 
1082  identity_pot%potential_type = do_potential_id
1083  CALL calc_2c_virial(work_virial_ovlp, virial_ovlp, 1.0_dp, qs_env, force_data%nl_2c_ovlp, &
1084  basis_set_ao, basis_set_ao, identity_pot)
1085  CALL dbcsr_release(virial_ovlp)
1086 
1087  DO k_xyz = 1, 3
1088  DO j_xyz = 1, 3
1089  DO i_xyz = 1, 3
1090  virial%pv_mp2(i_xyz, j_xyz) = virial%pv_mp2(i_xyz, j_xyz) &
1091  - work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1092  virial%pv_overlap(i_xyz, j_xyz) = virial%pv_overlap(i_xyz, j_xyz) &
1093  - work_virial_ovlp(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1094  virial%pv_virial(i_xyz, j_xyz) = virial%pv_virial(i_xyz, j_xyz) &
1095  - work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz) &
1096  - work_virial_ovlp(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1097  END DO
1098  END DO
1099  END DO
1100  END IF
1101 
1102  !Calculate the periodic contributions of (P|Q) to the force and the virial
1103  work_virial = 0.0_dp
1104  IF (force_data%do_periodic) THEN
1105  IF (mp2_env%eri_method == do_eri_gpw) THEN
1106  CALL get_2c_gpw_forces(force_data%G_PQ, force, work_virial, use_virial, mp2_env, qs_env)
1107  ELSE IF (mp2_env%eri_method == do_eri_mme) THEN
1108  CALL get_2c_mme_forces(force_data%G_PQ, force, mp2_env, qs_env)
1109  IF (use_virial) cpabort("Stress tensor not available with MME intrgrals")
1110  ELSE
1111  cpabort("Periodic case not possible with OS integrals")
1112  END IF
1113  CALL dbcsr_clear(force_data%G_PQ)
1114  END IF
1115 
1116  IF (use_virial) THEN
1117  virial%pv_mp2 = virial%pv_mp2 + work_virial
1118  virial%pv_virial = virial%pv_virial + work_virial
1119  virial%pv_calculate = .false.
1120 
1121  DO ibasis = 1, SIZE(basis_set_ao)
1122  orb_basis => basis_set_ao(ibasis)%gto_basis_set
1123  CALL init_interaction_radii_orb_basis(orb_basis, eps_pgf_orb_old)
1124  ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
1125  CALL init_interaction_radii_orb_basis(ri_basis, eps_pgf_orb_old)
1126  END DO
1127  END IF
1128 
1129  !clean-up
1130  IF (ASSOCIATED(dummy_ptr)) DEALLOCATE (dummy_ptr)
1131  DO ispin = 1, nspins
1132  CALL dbt_destroy(t_p(ispin))
1133  END DO
1134  CALL dbt_destroy(t_3c_0)
1135  CALL dbt_destroy(t_3c_1)
1136  CALL dbt_destroy(t_3c_3)
1137  CALL dbt_destroy(t_3c_4)
1138  CALL dbt_destroy(t_3c_5)
1139  CALL dbt_destroy(t_3c_6)
1140  CALL dbt_destroy(t_3c_7)
1141  CALL dbt_destroy(t_3c_8)
1142  CALL dbt_destroy(t_3c_sparse)
1143  CALL dbt_destroy(t_3c_help_1)
1144  CALL dbt_destroy(t_3c_help_2)
1145  CALL dbt_destroy(t_3c_ints)
1146  CALL dbt_destroy(t_3c_work)
1147  CALL dbt_destroy(t_r_occ)
1148  CALL dbt_destroy(t_r_virt)
1149  CALL dbt_destroy(t_dm_occ)
1150  CALL dbt_destroy(t_dm_virt)
1151  CALL dbt_destroy(t_q)
1152  CALL dbt_destroy(t_kqkt)
1153  CALL dbt_destroy(t_m_occ)
1154  CALL dbt_destroy(t_m_virt)
1155  CALL dbcsr_release(r_occ)
1156  CALL dbcsr_release(r_virt)
1157  CALL dbcsr_release(dbcsr_work1)
1158  CALL dbcsr_release(dbcsr_work2)
1159  CALL dbcsr_release(dbcsr_work3)
1160  CALL dbcsr_release(exp_occ)
1161  CALL dbcsr_release(exp_virt)
1162 
1163  CALL dbt_destroy(t_2c_ri)
1164  CALL dbt_destroy(t_2c_ri_2)
1165  CALL dbt_destroy(t_2c_ao)
1166  CALL dbcsr_deallocate_matrix_set(mat_dm_occ)
1167  CALL dbcsr_deallocate_matrix_set(mat_dm_virt)
1168 
1169  CALL timestop(handle)
1170 
1171  END SUBROUTINE calc_laplace_loop_forces
1172 
1173 ! **************************************************************************************************
1174 !> \brief Updates the cubic-scaling RPA contribution to the forces at each quadrature point. This
1175 !> routine is adapted from the corresponding Laplace SOS-MP2 loop force one.
1176 !> \param force_data ...
1177 !> \param mat_P_omega ...
1178 !> \param t_3c_M ...
1179 !> \param t_3c_O ...
1180 !> \param t_3c_O_compressed ...
1181 !> \param t_3c_O_ind ...
1182 !> \param fm_scaled_dm_occ_tau ...
1183 !> \param fm_scaled_dm_virt_tau ...
1184 !> \param fm_mo_coeff_occ ...
1185 !> \param fm_mo_coeff_virt ...
1186 !> \param fm_mo_coeff_occ_scaled ...
1187 !> \param fm_mo_coeff_virt_scaled ...
1188 !> \param starts_array_mc ...
1189 !> \param ends_array_mc ...
1190 !> \param starts_array_mc_block ...
1191 !> \param ends_array_mc_block ...
1192 !> \param num_integ_points ...
1193 !> \param nmo ...
1194 !> \param Eigenval ...
1195 !> \param e_fermi ...
1196 !> \param weights_cos_tf_t_to_w ...
1197 !> \param weights_cos_tf_w_to_t ...
1198 !> \param tj ...
1199 !> \param wj ...
1200 !> \param tau_tj ...
1201 !> \param cut_memory ...
1202 !> \param ispin ...
1203 !> \param open_shell ...
1204 !> \param unit_nr ...
1205 !> \param dbcsr_time ...
1206 !> \param dbcsr_nflop ...
1207 !> \param mp2_env ...
1208 !> \param qs_env ...
1209 ! **************************************************************************************************
1210  SUBROUTINE calc_rpa_loop_forces(force_data, mat_P_omega, t_3c_M, t_3c_O, t_3c_O_compressed, &
1211  t_3c_O_ind, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, &
1212  fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, &
1213  fm_mo_coeff_virt_scaled, starts_array_mc, ends_array_mc, &
1214  starts_array_mc_block, ends_array_mc_block, num_integ_points, &
1215  nmo, Eigenval, e_fermi, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, &
1216  tj, wj, tau_tj, cut_memory, ispin, open_shell, unit_nr, dbcsr_time, &
1217  dbcsr_nflop, mp2_env, qs_env)
1218 
1219  TYPE(im_time_force_type), INTENT(INOUT) :: force_data
1220  TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: mat_p_omega
1221  TYPE(dbt_type), INTENT(INOUT) :: t_3c_m, t_3c_o
1222  TYPE(hfx_compression_type), DIMENSION(:) :: t_3c_o_compressed
1223  TYPE(block_ind_type), DIMENSION(:), INTENT(INOUT) :: t_3c_o_ind
1224  TYPE(cp_fm_type), INTENT(IN) :: fm_scaled_dm_occ_tau, &
1225  fm_scaled_dm_virt_tau
1226  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mo_coeff_occ, fm_mo_coeff_virt
1227  TYPE(cp_fm_type), INTENT(IN) :: fm_mo_coeff_occ_scaled, &
1228  fm_mo_coeff_virt_scaled
1229  INTEGER, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
1230  starts_array_mc_block, &
1231  ends_array_mc_block
1232  INTEGER, INTENT(IN) :: num_integ_points, nmo
1233  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: eigenval
1234  REAL(kind=dp), INTENT(IN) :: e_fermi
1235  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :), &
1236  INTENT(IN) :: weights_cos_tf_t_to_w, &
1237  weights_cos_tf_w_to_t
1238  REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
1239  INTENT(IN) :: tj, wj
1240  REAL(kind=dp), DIMENSION(0:num_integ_points), &
1241  INTENT(IN) :: tau_tj
1242  INTEGER, INTENT(IN) :: cut_memory, ispin
1243  LOGICAL, INTENT(IN) :: open_shell
1244  INTEGER, INTENT(IN) :: unit_nr
1245  REAL(dp), INTENT(INOUT) :: dbcsr_time
1246  INTEGER(int_8), INTENT(INOUT) :: dbcsr_nflop
1247  TYPE(mp2_type) :: mp2_env
1248  TYPE(qs_environment_type), POINTER :: qs_env
1249 
1250  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_rpa_loop_forces'
1251 
1252  INTEGER :: dummy_int, handle, handle2, i_mem, i_xyz, ibasis, iquad, j_xyz, jquad, k_xyz, &
1253  n_mem_ri, n_rep, natom, nkind, nspins, unit_nr_dbcsr
1254  INTEGER(int_8) :: flop, nze, nze_ddint, nze_der_ao, &
1255  nze_der_ri, nze_kbk
1256  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, batch_blk_end, &
1257  batch_blk_start, batch_end_ri, &
1258  batch_start_ri, kind_of, mc_ranges, &
1259  mc_ranges_ri
1260  INTEGER, DIMENSION(:, :), POINTER :: dummy_ptr
1261  LOGICAL :: memory_info, use_virial
1262  REAL(dp) :: eps_filter, eps_pgf_orb, eps_pgf_orb_old, fac, occ, occ_ddint, occ_der_ao, &
1263  occ_der_ri, occ_kbk, omega, pref, spin_fac, t1, t2, tau, weight
1264  REAL(dp), DIMENSION(3, 3) :: work_virial, work_virial_ovlp
1265  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1266  TYPE(cell_type), POINTER :: cell
1267  TYPE(cp_blacs_env_type), POINTER :: blacs_env
1268  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_p_tau, matrix_ks, matrix_s
1269  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_dm_occ, mat_dm_virt
1270  TYPE(dbcsr_type) :: dbcsr_work1, dbcsr_work2, dbcsr_work3, &
1271  dbcsr_work_symm, exp_occ, exp_virt, &
1272  r_occ, r_virt, virial_ovlp, y_1, y_2
1273  TYPE(dbt_type) :: t_2c_ao, t_2c_ri, t_2c_ri_2, t_2c_tmp, t_3c_0, t_3c_1, t_3c_3, t_3c_4, &
1274  t_3c_5, t_3c_6, t_3c_7, t_3c_8, t_3c_help_1, t_3c_help_2, t_3c_ints, t_3c_sparse, &
1275  t_3c_work, t_dm_occ, t_dm_virt, t_kbkt, t_m_occ, t_m_virt, t_p, t_r_occ, t_r_virt
1276  TYPE(dbt_type), ALLOCATABLE, DIMENSION(:) :: t_b
1277  TYPE(dft_control_type), POINTER :: dft_control
1278  TYPE(gto_basis_set_p_type), ALLOCATABLE, &
1279  DIMENSION(:), TARGET :: basis_set_ao, basis_set_ri_aux
1280  TYPE(gto_basis_set_type), POINTER :: orb_basis, ri_basis
1281  TYPE(libint_potential_type) :: identity_pot
1282  TYPE(mp_para_env_type), POINTER :: para_env
1283  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1284  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1285  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1286  TYPE(section_vals_type), POINTER :: qs_section
1287  TYPE(virial_type), POINTER :: virial
1288 
1289  NULLIFY (matrix_s, dummy_ptr, atomic_kind_set, force, matrix_s, matrix_ks, mat_dm_occ, mat_dm_virt)
1290  NULLIFY (dft_control, virial, particle_set, cell, blacs_env, para_env, orb_basis, ri_basis)
1291  NULLIFY (qs_kind_set)
1292 
1293  CALL timeset(routinen, handle)
1294 
1295  CALL get_qs_env(qs_env, matrix_s=matrix_s, natom=natom, atomic_kind_set=atomic_kind_set, &
1296  force=force, matrix_ks=matrix_ks, dft_control=dft_control, virial=virial, &
1297  particle_set=particle_set, cell=cell, blacs_env=blacs_env, para_env=para_env, &
1298  qs_kind_set=qs_kind_set, nkind=nkind)
1299  eps_filter = mp2_env%ri_rpa_im_time%eps_filter
1300  nspins = dft_control%nspins
1301 
1302  memory_info = mp2_env%ri_rpa_im_time%memory_info
1303  IF (memory_info) THEN
1304  unit_nr_dbcsr = unit_nr
1305  ELSE
1306  unit_nr_dbcsr = 0
1307  END IF
1308 
1309  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1310 
1311  IF (use_virial) virial%pv_calculate = .true.
1312 
1313  IF (use_virial) THEN
1314  qs_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS")
1315  CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", n_rep_val=n_rep)
1316  IF (n_rep /= 0) THEN
1317  CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", r_val=eps_pgf_orb)
1318  ELSE
1319  CALL section_vals_val_get(qs_section, "EPS_DEFAULT", r_val=eps_pgf_orb)
1320  eps_pgf_orb = sqrt(eps_pgf_orb)
1321  END IF
1322  eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
1323 
1324  ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind))
1325  CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
1326  CALL basis_set_list_setup(basis_set_ao, "ORB", qs_kind_set)
1327 
1328  DO ibasis = 1, SIZE(basis_set_ao)
1329  orb_basis => basis_set_ao(ibasis)%gto_basis_set
1330  CALL init_interaction_radii_orb_basis(orb_basis, eps_pgf_orb)
1331  ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
1332  CALL init_interaction_radii_orb_basis(ri_basis, eps_pgf_orb)
1333  END DO
1334  END IF
1335 
1336  !We follow the general logic of the compute_mat_P_omega routine
1337  CALL dbt_create(force_data%t_2c_K, t_2c_ri)
1338  CALL dbt_create(force_data%t_2c_K, t_2c_ri_2)
1339  CALL dbt_create(force_data%t_2c_der_ovlp(1), t_2c_ao)
1340 
1341  CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
1342 
1343  ! Always do the batching of the MO on mu and sigma, such that it is consistent between
1344  ! the occupied and the virtual quantities
1345  ALLOCATE (mc_ranges(cut_memory + 1))
1346  mc_ranges(:cut_memory) = starts_array_mc_block(:)
1347  mc_ranges(cut_memory + 1) = ends_array_mc_block(cut_memory) + 1
1348 
1349  ! Also need some batching on the RI, because it loses sparsity at some point
1350  n_mem_ri = cut_memory
1351  CALL create_tensor_batches(force_data%bsizes_RI_split, n_mem_ri, batch_start_ri, batch_end_ri, &
1352  batch_blk_start, batch_blk_end)
1353  ALLOCATE (mc_ranges_ri(n_mem_ri + 1))
1354  mc_ranges_ri(1:n_mem_ri) = batch_blk_start(1:n_mem_ri)
1355  mc_ranges_ri(n_mem_ri + 1) = batch_blk_end(n_mem_ri) + 1
1356  DEALLOCATE (batch_blk_start, batch_blk_end)
1357 
1358  !Pre-allocate all required tensors and matrices
1359  CALL dbt_create(t_2c_ri, t_p)
1360  CALL dbt_create(t_2c_ri, t_kbkt)
1361  CALL dbt_create(t_2c_ao, t_dm_occ)
1362  CALL dbt_create(t_2c_ao, t_dm_virt)
1363 
1364  !note: t_3c_O and t_3c_M have different mappings (map_1d, map_2d)
1365  CALL dbt_create(t_3c_o, t_m_occ)
1366  CALL dbt_create(t_3c_o, t_m_virt)
1367  CALL dbt_create(t_3c_o, t_3c_0)
1368 
1369  CALL dbt_create(t_3c_o, t_3c_1)
1370  CALL dbt_create(t_3c_o, t_3c_3)
1371  CALL dbt_create(t_3c_o, t_3c_4)
1372  CALL dbt_create(t_3c_o, t_3c_5)
1373  CALL dbt_create(t_3c_m, t_3c_6)
1374  CALL dbt_create(t_3c_m, t_3c_7)
1375  CALL dbt_create(t_3c_m, t_3c_8)
1376  CALL dbt_create(t_3c_m, t_3c_sparse)
1377  CALL dbt_create(t_3c_o, t_3c_help_1)
1378  CALL dbt_create(t_3c_o, t_3c_help_2)
1379  CALL dbt_create(t_2c_ao, t_r_occ)
1380  CALL dbt_create(t_2c_ao, t_r_virt)
1381  CALL dbt_create(t_3c_m, t_3c_ints)
1382  CALL dbt_create(t_3c_m, t_3c_work)
1383 
1384  !Before entring the loop, need to compute the 2c tensors B = (1 + Q(w))^-1 - 1, for each
1385  !frequency grid point, before doing the transformation to the time grid
1386  ALLOCATE (t_b(num_integ_points))
1387  DO jquad = 1, num_integ_points
1388  CALL dbt_create(t_2c_ri, t_b(jquad))
1389  END DO
1390 
1391  ALLOCATE (mat_p_tau(num_integ_points))
1392  DO jquad = 1, num_integ_points
1393  ALLOCATE (mat_p_tau(jquad)%matrix)
1394  CALL dbcsr_create(mat_p_tau(jquad)%matrix, template=mat_p_omega(jquad, ispin)%matrix)
1395  END DO
1396 
1397  CALL dbcsr_create(dbcsr_work_symm, template=force_data%G_PQ, matrix_type=dbcsr_type_symmetric)
1398  CALL dbt_create(dbcsr_work_symm, t_2c_tmp)
1399 
1400  !loop over freqeuncies
1401  DO iquad = 1, num_integ_points
1402  omega = tj(iquad)
1403 
1404  !calculate (1 + Q(w))^-1 - 1 for the given freq.
1405  !Always take spin alpha (get 2*alpha in closed shell, and alpha+beta in open-shell)
1406  CALL dbcsr_copy(dbcsr_work_symm, mat_p_omega(iquad, 1)%matrix)
1407  CALL dbt_copy_matrix_to_tensor(dbcsr_work_symm, t_2c_tmp)
1408  CALL dbt_copy(t_2c_tmp, t_2c_ri, move_data=.true.)
1409 
1410  CALL dbt_contract(1.0_dp, t_2c_ri, force_data%t_2c_K, 0.0_dp, t_2c_ri_2, &
1411  contract_1=[2], notcontract_1=[1], &
1412  contract_2=[1], notcontract_2=[2], &
1413  map_1=[1], map_2=[2], filter_eps=eps_filter, &
1414  flop=flop, unit_nr=unit_nr_dbcsr)
1415  dbcsr_nflop = dbcsr_nflop + flop
1416  CALL dbt_contract(1.0_dp, force_data%t_2c_K, t_2c_ri_2, 0.0_dp, t_2c_ri, &
1417  contract_1=[1], notcontract_1=[2], &
1418  contract_2=[1], notcontract_2=[2], &
1419  map_1=[1], map_2=[2], filter_eps=eps_filter, &
1420  flop=flop, unit_nr=unit_nr_dbcsr)
1421  CALL dbt_copy(t_2c_ri, t_2c_tmp, move_data=.true.)
1422  CALL dbt_copy_tensor_to_matrix(t_2c_tmp, dbcsr_work_symm)
1423  CALL dbcsr_add_on_diag(dbcsr_work_symm, 1.0_dp)
1424 
1425  CALL cp_dbcsr_cholesky_decompose(dbcsr_work_symm, para_env=para_env, blacs_env=blacs_env)
1426  CALL cp_dbcsr_cholesky_invert(dbcsr_work_symm, para_env=para_env, blacs_env=blacs_env, upper_to_full=.true.)
1427 
1428  CALL dbcsr_add_on_diag(dbcsr_work_symm, -1.0_dp)
1429 
1430  DO jquad = 1, num_integ_points
1431  tau = tau_tj(jquad)
1432 
1433  !the P matrix to time.
1434  weight = weights_cos_tf_w_to_t(jquad, iquad)*cos(tau*omega)
1435  IF (open_shell) THEN
1436  IF (ispin == 1) THEN
1437  !mat_P_omega contains the sum of alpha and beta spin => we only want alpha
1438  CALL dbcsr_add(mat_p_tau(jquad)%matrix, mat_p_omega(iquad, 1)%matrix, 1.0_dp, weight)
1439  CALL dbcsr_add(mat_p_tau(jquad)%matrix, mat_p_omega(iquad, 2)%matrix, 1.0_dp, -weight)
1440  ELSE
1441  CALL dbcsr_add(mat_p_tau(jquad)%matrix, mat_p_omega(iquad, 2)%matrix, 1.0_dp, weight)
1442  END IF
1443  ELSE
1444  !factor 0.5 because originam matrix Q is scaled by 2 in RPA (spin)
1445  weight = 0.5_dp*weight
1446  CALL dbcsr_add(mat_p_tau(jquad)%matrix, mat_p_omega(iquad, 1)%matrix, 1.0_dp, weight)
1447  END IF
1448 
1449  !convert B matrix to time
1450  weight = weights_cos_tf_t_to_w(iquad, jquad)*cos(tau*omega)*wj(iquad)
1451  CALL dbt_copy_matrix_to_tensor(dbcsr_work_symm, t_2c_tmp)
1452  CALL dbt_scale(t_2c_tmp, weight)
1453  CALL dbt_copy(t_2c_tmp, t_b(jquad), summation=.true., move_data=.true.)
1454  END DO
1455  END DO
1456  CALL dbt_destroy(t_2c_tmp)
1457  CALL dbcsr_release(dbcsr_work_symm)
1458  CALL dbt_clear(t_2c_ri)
1459  CALL dbt_clear(t_2c_ri_2)
1460 
1461  !Pre-define the sparsity of t_3c_4 as a function of the derivatives
1462  occ_der_ao = 0; nze_der_ao = 0
1463  occ_der_ri = 0; nze_der_ri = 0
1464  DO i_xyz = 1, 3
1465  DO i_mem = 1, cut_memory
1466  CALL decompress_tensor(force_data%t_3c_der_RI(i_xyz), force_data%t_3c_der_RI_ind(i_mem, i_xyz)%ind, &
1467  force_data%t_3c_der_RI_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
1468  CALL get_tensor_occupancy(force_data%t_3c_der_RI(i_xyz), nze, occ)
1469  occ_der_ri = occ_der_ri + occ
1470  nze_der_ri = nze_der_ri + nze
1471  CALL dbt_copy(force_data%t_3c_der_RI(i_xyz), t_3c_sparse, summation=.true., move_data=.true.)
1472 
1473  CALL decompress_tensor(force_data%t_3c_der_AO(i_xyz), force_data%t_3c_der_AO_ind(i_mem, i_xyz)%ind, &
1474  force_data%t_3c_der_AO_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
1475  CALL get_tensor_occupancy(force_data%t_3c_der_AO(i_xyz), nze, occ)
1476  occ_der_ao = occ_der_ao + occ
1477  nze_der_ao = nze_der_ao + nze
1478  CALL dbt_copy(force_data%t_3c_der_AO(i_xyz), t_3c_sparse, order=[1, 3, 2], summation=.true.)
1479  CALL dbt_copy(force_data%t_3c_der_AO(i_xyz), t_3c_sparse, summation=.true., move_data=.true.)
1480  END DO
1481  END DO
1482  occ_der_ri = occ_der_ri/3.0_dp
1483  occ_der_ao = occ_der_ao/3.0_dp
1484  nze_der_ri = nze_der_ri/3
1485  nze_der_ao = nze_der_ao/3
1486 
1487  CALL dbcsr_create(r_occ, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1488  CALL dbcsr_create(r_virt, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1489  CALL dbcsr_create(dbcsr_work_symm, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_symmetric)
1490  CALL dbcsr_create(dbcsr_work1, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1491  CALL dbcsr_create(dbcsr_work2, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1492  CALL dbcsr_create(dbcsr_work3, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1493  CALL dbcsr_create(exp_occ, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1494  CALL dbcsr_create(exp_virt, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1495  IF (use_virial) CALL dbcsr_create(virial_ovlp, template=dbcsr_work1)
1496 
1497  CALL dbt_batched_contract_init(t_3c_0, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
1498  CALL dbt_batched_contract_init(t_3c_1, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
1499  CALL dbt_batched_contract_init(t_3c_3, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1500  CALL dbt_batched_contract_init(t_m_occ, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1501  CALL dbt_batched_contract_init(t_m_virt, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1502 
1503  CALL dbt_batched_contract_init(t_3c_ints, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1504  CALL dbt_batched_contract_init(t_3c_work, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1505 
1506  CALL dbt_batched_contract_init(t_3c_4, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1507  batch_range_3=mc_ranges)
1508  CALL dbt_batched_contract_init(t_3c_5, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1509  batch_range_3=mc_ranges)
1510  CALL dbt_batched_contract_init(t_3c_6, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1511  batch_range_3=mc_ranges)
1512  CALL dbt_batched_contract_init(t_3c_7, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1513  batch_range_3=mc_ranges)
1514  CALL dbt_batched_contract_init(t_3c_8, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1515  batch_range_3=mc_ranges)
1516  CALL dbt_batched_contract_init(t_3c_sparse, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1517  batch_range_3=mc_ranges)
1518 
1519  fac = 1.0_dp/fourpi*mp2_env%ri_rpa%scale_rpa
1520  IF (open_shell) fac = 0.5_dp*fac
1521 
1522  work_virial = 0.0_dp
1523  work_virial_ovlp = 0.0_dp
1524  DO jquad = 1, num_integ_points
1525  tau = tau_tj(jquad)
1526  occ_ddint = 0; nze_ddint = 0
1527 
1528  CALL para_env%sync()
1529  t1 = m_walltime()
1530 
1531  !Deal with the force contributions where there is no explicit 3-center quantities, i.e. the
1532  !forces due to the metric and potential derivatives
1533  CALL dbt_create(mat_p_tau(jquad)%matrix, t_2c_tmp)
1534  CALL dbt_copy_matrix_to_tensor(mat_p_tau(jquad)%matrix, t_2c_tmp)
1535  CALL dbt_copy(t_2c_tmp, t_p, move_data=.true.)
1536  CALL dbt_filter(t_p, eps_filter)
1537  CALL dbt_destroy(t_2c_tmp)
1538 
1539  CALL perform_2c_ops(force, t_kbkt, force_data, fac, t_b(jquad), t_p, t_2c_ri, t_2c_ri_2, &
1540  use_virial, atom_of_kind, kind_of, eps_filter, dbcsr_nflop, unit_nr_dbcsr)
1541  CALL get_tensor_occupancy(t_kbkt, nze_kbk, occ_kbk)
1542 
1543  !Calculate the pseudo-density matrix in tensor form. There are a few useless arguments for SOS-MP2
1544  CALL compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, tau_tj, num_integ_points, &
1545  nmo, fm_mo_coeff_occ(ispin), fm_mo_coeff_virt(ispin), &
1546  fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm_occ, mat_dm_virt, &
1547  matrix_s, ispin, eigenval(:, ispin), e_fermi, eps_filter, &
1548  mp2_env%ri_rpa_im_time%memory_info, unit_nr, &
1549  jquad, .false., .false., qs_env, dummy_int, dummy_ptr, para_env)
1550 
1551  CALL dbt_create(mat_dm_occ(jquad, 1)%matrix, t_2c_tmp)
1552  CALL dbt_copy_matrix_to_tensor(mat_dm_occ(jquad, 1)%matrix, t_2c_tmp)
1553  CALL dbt_copy(t_2c_tmp, t_dm_occ, move_data=.true.)
1554  CALL dbt_filter(t_dm_occ, eps_filter)
1555  CALL dbt_destroy(t_2c_tmp)
1556 
1557  CALL dbt_create(mat_dm_virt(jquad, 1)%matrix, t_2c_tmp)
1558  CALL dbt_copy_matrix_to_tensor(mat_dm_virt(jquad, 1)%matrix, t_2c_tmp)
1559  CALL dbt_copy(t_2c_tmp, t_dm_virt, move_data=.true.)
1560  CALL dbt_filter(t_dm_virt, eps_filter)
1561  CALL dbt_destroy(t_2c_tmp)
1562 
1563  !Deal with the 3-center quantities.
1564  CALL perform_3c_ops(force, t_r_occ, t_r_virt, force_data, fac, cut_memory, n_mem_ri, &
1565  t_kbkt, t_dm_occ, t_dm_virt, t_3c_o, t_3c_m, t_m_occ, t_m_virt, t_3c_0, t_3c_1, &
1566  t_3c_3, t_3c_4, t_3c_5, t_3c_6, t_3c_7, t_3c_8, t_3c_sparse, t_3c_help_1, t_3c_help_2, &
1567  t_3c_ints, t_3c_work, starts_array_mc, ends_array_mc, batch_start_ri, &
1568  batch_end_ri, t_3c_o_compressed, t_3c_o_ind, use_virial, &
1569  atom_of_kind, kind_of, eps_filter, occ_ddint, nze_ddint, dbcsr_nflop, &
1570  unit_nr_dbcsr, mp2_env)
1571 
1572  CALL timeset(routinen//"_dbcsr", handle2)
1573  !We go back to DBCSR matrices from now on
1574  !Note: R matrices are in fact symmetric, but use a normal type for convenience
1575  CALL dbt_create(matrix_s(1)%matrix, t_2c_tmp)
1576  CALL dbt_copy(t_r_occ, t_2c_tmp, move_data=.true.)
1577  CALL dbt_copy_tensor_to_matrix(t_2c_tmp, r_occ)
1578 
1579  CALL dbt_copy(t_r_virt, t_2c_tmp, move_data=.true.)
1580  CALL dbt_copy_tensor_to_matrix(t_2c_tmp, r_virt)
1581 
1582  !Iteratively calculate the Y1 and Y2 matrices
1583  CALL dbcsr_copy(dbcsr_work_symm, matrix_ks(ispin)%matrix)
1584  CALL dbcsr_add(dbcsr_work_symm, matrix_s(1)%matrix, 1.0_dp, -e_fermi)
1585  CALL dbcsr_multiply('N', 'N', tau, force_data%P_occ(ispin)%matrix, &
1586  dbcsr_work_symm, 0.0_dp, dbcsr_work1)
1587  CALL build_y_matrix(y_1, dbcsr_work1, force_data%P_occ(ispin)%matrix, r_virt, eps_filter)
1588  CALL matrix_exponential(exp_occ, dbcsr_work1, 1.0_dp, 1.0_dp, eps_filter)
1589 
1590  CALL dbcsr_multiply('N', 'N', -tau, force_data%P_virt(ispin)%matrix, &
1591  dbcsr_work_symm, 0.0_dp, dbcsr_work1)
1592  CALL build_y_matrix(y_2, dbcsr_work1, force_data%P_virt(ispin)%matrix, r_occ, eps_filter)
1593  CALL matrix_exponential(exp_virt, dbcsr_work1, 1.0_dp, 1.0_dp, eps_filter)
1594 
1595  !The force contribution coming from [-S^-1*(e^-tau*P_virt*F)^T*R_occ*S^-1
1596  ! +tau*S^-1*Y_2^T*F*S^-1] * der_S
1597  !as well as -tau*e_fermi*Y_1*P^occ + tau*e_fermi*Y_2*P^virt
1598  CALL dbcsr_multiply('N', 'N', 1.0_dp, r_occ, force_data%inv_ovlp, 0.0_dp, dbcsr_work1)
1599  CALL dbcsr_multiply('T', 'N', 1.0_dp, exp_virt, dbcsr_work1, 0.0_dp, dbcsr_work3)
1600  CALL dbcsr_multiply('N', 'N', 1.0_dp, force_data%inv_ovlp, dbcsr_work3, 0.0_dp, dbcsr_work2)
1601 
1602  CALL dbcsr_multiply('N', 'T', tau, force_data%inv_ovlp, y_2, 0.0_dp, dbcsr_work3)
1603  CALL dbcsr_multiply('N', 'N', 1.0_dp, dbcsr_work3, dbcsr_work_symm, 0.0_dp, dbcsr_work1)
1604  CALL dbcsr_multiply('N', 'N', -1.0_dp, dbcsr_work1, force_data%inv_ovlp, 1.0_dp, dbcsr_work2)
1605 
1606  CALL dbcsr_multiply('N', 'T', tau*e_fermi, force_data%P_occ(ispin)%matrix, y_1, 1.0_dp, dbcsr_work2)
1607  CALL dbcsr_multiply('N', 'T', -tau*e_fermi, force_data%P_virt(ispin)%matrix, y_2, 1.0_dp, dbcsr_work2)
1608 
1609  CALL dbt_copy_matrix_to_tensor(dbcsr_work2, t_2c_tmp)
1610  CALL dbt_copy(t_2c_tmp, t_2c_ao, move_data=.true.)
1611 
1612  pref = -1.0_dp*fac
1613  CALL get_2c_der_force(force, t_2c_ao, force_data%t_2c_der_ovlp, atom_of_kind, &
1614  kind_of, force_data%idx_to_at_AO, pref, do_ovlp=.true.)
1615 
1616  IF (use_virial) CALL dbcsr_add(virial_ovlp, dbcsr_work2, 1.0_dp, pref)
1617 
1618  !The final contribution from Tr[(tau*Y_1*P_occ - tau*Y_2*P_virt) * der_F]
1619  CALL dbcsr_multiply('N', 'N', fac*tau, y_1, force_data%P_occ(ispin)%matrix, 1.0_dp, &
1620  force_data%sum_YP_tau(ispin)%matrix, retain_sparsity=.true.)
1621  CALL dbcsr_multiply('N', 'N', -fac*tau, y_2, force_data%P_virt(ispin)%matrix, 1.0_dp, &
1622  force_data%sum_YP_tau(ispin)%matrix, retain_sparsity=.true.)
1623 
1624  spin_fac = 0.5_dp*fac
1625  IF (open_shell) spin_fac = 2.0_dp*spin_fac
1626  !Build-up the RHS of the response equation.
1627  CALL dbcsr_multiply('N', 'N', 1.0_dp*spin_fac, r_virt, exp_occ, 1.0_dp, &
1628  force_data%sum_O_tau(ispin)%matrix, retain_sparsity=.true.)
1629  CALL dbcsr_multiply('N', 'N', -1.0_dp*spin_fac, r_occ, exp_virt, 1.0_dp, &
1630  force_data%sum_O_tau(ispin)%matrix, retain_sparsity=.true.)
1631  CALL dbcsr_multiply('N', 'N', tau*spin_fac, dbcsr_work_symm, y_1, 1.0_dp, &
1632  force_data%sum_O_tau(ispin)%matrix, retain_sparsity=.true.)
1633  CALL dbcsr_multiply('N', 'N', tau*spin_fac, dbcsr_work_symm, y_2, 1.0_dp, &
1634  force_data%sum_O_tau(ispin)%matrix, retain_sparsity=.true.)
1635 
1636  CALL timestop(handle2)
1637 
1638  !Print some info
1639  CALL para_env%sync()
1640  t2 = m_walltime()
1641  dbcsr_time = dbcsr_time + t2 - t1
1642 
1643  IF (unit_nr > 0) THEN
1644  WRITE (unit_nr, '(/T3,A,1X,I3,A)') &
1645  'RPA_LOW_SCALING_INFO| Info for time point', jquad, ' (gradients)'
1646  WRITE (unit_nr, '(T6,A,T56,F25.6)') &
1647  'Time:', t2 - t1
1648  WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1649  'Occupancy of 3c AO derivs:', real(nze_der_ao, dp), '/', occ_der_ao*100, '%'
1650  WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1651  'Occupancy of 3c RI derivs:', real(nze_der_ri, dp), '/', occ_der_ri*100, '%'
1652  WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1653  'Occupancy of the Docc * Dvirt * 3c-int tensor', real(nze_ddint, dp), '/', occ_ddint*100, '%'
1654  WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1655  'Occupancy of KBK^T 2c-tensor:', real(nze_kbk, dp), '/', occ_kbk*100, '%'
1656  CALL m_flush(unit_nr)
1657  END IF
1658 
1659  !intermediate clean-up
1660  CALL dbcsr_release(y_1)
1661  CALL dbcsr_release(y_2)
1662  CALL dbt_destroy(t_2c_tmp)
1663 
1664  END DO !jquad
1665 
1666  CALL dbt_batched_contract_finalize(t_3c_0)
1667  CALL dbt_batched_contract_finalize(t_3c_1)
1668  CALL dbt_batched_contract_finalize(t_3c_3)
1669  CALL dbt_batched_contract_finalize(t_m_occ)
1670  CALL dbt_batched_contract_finalize(t_m_virt)
1671 
1672  CALL dbt_batched_contract_finalize(t_3c_ints)
1673  CALL dbt_batched_contract_finalize(t_3c_work)
1674 
1675  CALL dbt_batched_contract_finalize(t_3c_4)
1676  CALL dbt_batched_contract_finalize(t_3c_5)
1677  CALL dbt_batched_contract_finalize(t_3c_6)
1678  CALL dbt_batched_contract_finalize(t_3c_7)
1679  CALL dbt_batched_contract_finalize(t_3c_8)
1680  CALL dbt_batched_contract_finalize(t_3c_sparse)
1681 
1682  !Calculate the 2c and 3c contributions to the virial
1683  IF (use_virial) THEN
1684  CALL dbt_copy(force_data%t_3c_virial_split, force_data%t_3c_virial, move_data=.true.)
1685  CALL calc_3c_virial(work_virial, force_data%t_3c_virial, 1.0_dp, qs_env, force_data%nl_3c, &
1686  basis_set_ri_aux, basis_set_ao, basis_set_ao, mp2_env%ri_metric, &
1687  der_eps=mp2_env%ri_rpa_im_time%eps_filter, op_pos=1)
1688 
1689  CALL calc_2c_virial(work_virial, force_data%RI_virial_met, 1.0_dp, qs_env, force_data%nl_2c_met, &
1690  basis_set_ri_aux, basis_set_ri_aux, mp2_env%ri_metric)
1691  CALL dbcsr_clear(force_data%RI_virial_met)
1692 
1693  IF (.NOT. force_data%do_periodic) THEN
1694  CALL calc_2c_virial(work_virial, force_data%RI_virial_pot, 1.0_dp, qs_env, force_data%nl_2c_pot, &
1695  basis_set_ri_aux, basis_set_ri_aux, mp2_env%potential_parameter)
1696  CALL dbcsr_clear(force_data%RI_virial_pot)
1697  END IF
1698 
1699  identity_pot%potential_type = do_potential_id
1700  CALL calc_2c_virial(work_virial_ovlp, virial_ovlp, 1.0_dp, qs_env, force_data%nl_2c_ovlp, &
1701  basis_set_ao, basis_set_ao, identity_pot)
1702  CALL dbcsr_release(virial_ovlp)
1703 
1704  DO k_xyz = 1, 3
1705  DO j_xyz = 1, 3
1706  DO i_xyz = 1, 3
1707  virial%pv_mp2(i_xyz, j_xyz) = virial%pv_mp2(i_xyz, j_xyz) &
1708  - work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1709  virial%pv_overlap(i_xyz, j_xyz) = virial%pv_overlap(i_xyz, j_xyz) &
1710  - work_virial_ovlp(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1711  virial%pv_virial(i_xyz, j_xyz) = virial%pv_virial(i_xyz, j_xyz) &
1712  - work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz) &
1713  - work_virial_ovlp(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1714  END DO
1715  END DO
1716  END DO
1717  END IF
1718 
1719  !Calculate the periodic contributions of (P|Q) to the force and the virial
1720  work_virial = 0.0_dp
1721  IF (force_data%do_periodic) THEN
1722  IF (mp2_env%eri_method == do_eri_gpw) THEN
1723  CALL get_2c_gpw_forces(force_data%G_PQ, force, work_virial, use_virial, mp2_env, qs_env)
1724  ELSE IF (mp2_env%eri_method == do_eri_mme) THEN
1725  CALL get_2c_mme_forces(force_data%G_PQ, force, mp2_env, qs_env)
1726  IF (use_virial) cpabort("Stress tensor not available with MME intrgrals")
1727  ELSE
1728  cpabort("Periodic case not possible with OS integrals")
1729  END IF
1730  CALL dbcsr_clear(force_data%G_PQ)
1731  END IF
1732 
1733  IF (use_virial) THEN
1734  virial%pv_mp2 = virial%pv_mp2 + work_virial
1735  virial%pv_virial = virial%pv_virial + work_virial
1736  virial%pv_calculate = .false.
1737 
1738  DO ibasis = 1, SIZE(basis_set_ao)
1739  orb_basis => basis_set_ao(ibasis)%gto_basis_set
1740  CALL init_interaction_radii_orb_basis(orb_basis, eps_pgf_orb_old)
1741  ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
1742  CALL init_interaction_radii_orb_basis(ri_basis, eps_pgf_orb_old)
1743  END DO
1744  END IF
1745 
1746  !clean-up
1747  IF (ASSOCIATED(dummy_ptr)) DEALLOCATE (dummy_ptr)
1748  DO jquad = 1, num_integ_points
1749  CALL dbt_destroy(t_b(jquad))
1750  END DO
1751  CALL dbt_destroy(t_p)
1752  CALL dbt_destroy(t_3c_0)
1753  CALL dbt_destroy(t_3c_1)
1754  CALL dbt_destroy(t_3c_3)
1755  CALL dbt_destroy(t_3c_4)
1756  CALL dbt_destroy(t_3c_5)
1757  CALL dbt_destroy(t_3c_6)
1758  CALL dbt_destroy(t_3c_7)
1759  CALL dbt_destroy(t_3c_8)
1760  CALL dbt_destroy(t_3c_sparse)
1761  CALL dbt_destroy(t_3c_help_1)
1762  CALL dbt_destroy(t_3c_help_2)
1763  CALL dbt_destroy(t_3c_ints)
1764  CALL dbt_destroy(t_3c_work)
1765  CALL dbt_destroy(t_r_occ)
1766  CALL dbt_destroy(t_r_virt)
1767  CALL dbt_destroy(t_dm_occ)
1768  CALL dbt_destroy(t_dm_virt)
1769  CALL dbt_destroy(t_kbkt)
1770  CALL dbt_destroy(t_m_occ)
1771  CALL dbt_destroy(t_m_virt)
1772  CALL dbcsr_release(r_occ)
1773  CALL dbcsr_release(r_virt)
1774  CALL dbcsr_release(dbcsr_work_symm)
1775  CALL dbcsr_release(dbcsr_work1)
1776  CALL dbcsr_release(dbcsr_work2)
1777  CALL dbcsr_release(dbcsr_work3)
1778  CALL dbcsr_release(exp_occ)
1779  CALL dbcsr_release(exp_virt)
1780 
1781  CALL dbt_destroy(t_2c_ri)
1782  CALL dbt_destroy(t_2c_ri_2)
1783  CALL dbt_destroy(t_2c_ao)
1784  CALL dbcsr_deallocate_matrix_set(mat_dm_occ)
1785  CALL dbcsr_deallocate_matrix_set(mat_dm_virt)
1786  CALL dbcsr_deallocate_matrix_set(mat_p_tau)
1787 
1788  CALL timestop(handle)
1789 
1790  END SUBROUTINE calc_rpa_loop_forces
1791 
1792 ! **************************************************************************************************
1793 !> \brief This subroutines performs the 2c tensor operations that are common accros low-scaling RPA
1794 !> and SOS-MP2, including forces and virial
1795 !> \param force ...
1796 !> \param t_KBKT returns the 2c tensor product of K*B*K^T
1797 !> \param force_data ...
1798 !> \param fac ...
1799 !> \param t_B depending on RPA or SOS-MP2, t_B contains (1 + Q)^-1 - 1 or simply Q, respectively
1800 !> \param t_P ...
1801 !> \param t_2c_RI ...
1802 !> \param t_2c_RI_2 ...
1803 !> \param use_virial ...
1804 !> \param atom_of_kind ...
1805 !> \param kind_of ...
1806 !> \param eps_filter ...
1807 !> \param dbcsr_nflop ...
1808 !> \param unit_nr_dbcsr ...
1809 ! **************************************************************************************************
1810  SUBROUTINE perform_2c_ops(force, t_KBKT, force_data, fac, t_B, t_P, t_2c_RI, t_2c_RI_2, use_virial, &
1811  atom_of_kind, kind_of, eps_filter, dbcsr_nflop, unit_nr_dbcsr)
1812 
1813  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1814  TYPE(dbt_type), INTENT(INOUT) :: t_kbkt
1815  TYPE(im_time_force_type), INTENT(INOUT) :: force_data
1816  REAL(dp), INTENT(IN) :: fac
1817  TYPE(dbt_type), INTENT(INOUT) :: t_b, t_p, t_2c_ri, t_2c_ri_2
1818  LOGICAL, INTENT(IN) :: use_virial
1819  INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind, kind_of
1820  REAL(dp), INTENT(IN) :: eps_filter
1821  INTEGER(int_8), INTENT(INOUT) :: dbcsr_nflop
1822  INTEGER, INTENT(IN) :: unit_nr_dbcsr
1823 
1824  CHARACTER(LEN=*), PARAMETER :: routinen = 'perform_2c_ops'
1825 
1826  INTEGER :: handle
1827  INTEGER(int_8) :: flop
1828  REAL(dp) :: pref
1829  TYPE(dbt_type) :: t_2c_tmp, t_2c_virial
1830 
1831  CALL timeset(routinen, handle)
1832 
1833  IF (use_virial) CALL dbt_create(force_data%RI_virial_pot, t_2c_virial)
1834 
1835  !P^T*K*B + P*K*B^T (note we calculate and save K*B*K^T for later, and P=P^T)
1836  CALL dbt_contract(1.0_dp, force_data%t_2c_K, t_b, 0.0_dp, t_2c_ri, &
1837  contract_1=[2], notcontract_1=[1], &
1838  contract_2=[1], notcontract_2=[2], &
1839  map_1=[1], map_2=[2], filter_eps=eps_filter, &
1840  flop=flop, unit_nr=unit_nr_dbcsr)
1841  dbcsr_nflop = dbcsr_nflop + flop
1842 
1843  CALL dbt_contract(1.0_dp, t_2c_ri, force_data%t_2c_K, 0.0_dp, t_kbkt, &
1844  contract_1=[2], notcontract_1=[1], &
1845  contract_2=[2], notcontract_2=[1], &
1846  map_1=[1], map_2=[2], filter_eps=eps_filter, &
1847  flop=flop, unit_nr=unit_nr_dbcsr)
1848  dbcsr_nflop = dbcsr_nflop + flop
1849 
1850  CALL dbt_contract(2.0_dp, t_p, t_2c_ri, 0.0_dp, t_2c_ri_2, & !t_2c_RI_2 holds P^T*K*B
1851  contract_1=[2], notcontract_1=[1], &
1852  contract_2=[1], notcontract_2=[2], &
1853  map_1=[1], map_2=[2], filter_eps=eps_filter, &
1854  flop=flop, unit_nr=unit_nr_dbcsr)
1855  dbcsr_nflop = dbcsr_nflop + flop
1856  CALL dbt_clear(t_2c_ri)
1857  !t_2c_RI_2 currently holds 2*P^T*K*B = P^T*K*B + P*K*B^T (because of symmetry)
1858 
1859  !For the metric contribution, we need S^-1*(P^T*K*B + P*K*B^T)*K^T
1860  CALL dbt_contract(1.0_dp, force_data%t_2c_inv_metric, t_2c_ri_2, 0.0_dp, t_2c_ri, &
1861  contract_1=[2], notcontract_1=[1], &
1862  contract_2=[1], notcontract_2=[2], &
1863  map_1=[1], map_2=[2], filter_eps=eps_filter, &
1864  flop=flop, unit_nr=unit_nr_dbcsr)
1865  dbcsr_nflop = dbcsr_nflop + flop
1866 
1867  CALL dbt_contract(1.0_dp, t_2c_ri, force_data%t_2c_K, 0.0_dp, t_2c_ri_2, &
1868  contract_1=[2], notcontract_1=[1], &
1869  contract_2=[2], notcontract_2=[1], &
1870  map_1=[1], map_2=[2], filter_eps=eps_filter, &
1871  flop=flop, unit_nr=unit_nr_dbcsr)
1872  dbcsr_nflop = dbcsr_nflop + flop
1873 
1874  !Here we do the trace for the force
1875  pref = -1.0_dp*fac
1876  CALL get_2c_der_force(force, t_2c_ri_2, force_data%t_2c_der_metric, atom_of_kind, &
1877  kind_of, force_data%idx_to_at_RI, pref, do_mp2=.true.)
1878  IF (use_virial) THEN
1879  CALL dbt_copy(t_2c_ri_2, t_2c_virial)
1880  CALL dbt_scale(t_2c_virial, pref)
1881  CALL dbt_copy_tensor_to_matrix(t_2c_virial, force_data%RI_virial_met, summation=.true.)
1882  CALL dbt_clear(t_2c_virial)
1883  END IF
1884 
1885  !For the potential contribution, we need S^-1*(P^T*K*B + P*K*B^T)*V^-0.5
1886  !some of it is still in t_2c_RI: ( S^-1*(P^T*K*B + P*K*B^T) )
1887  CALL dbt_contract(1.0_dp, t_2c_ri, force_data%t_2c_pot_msqrt, 0.0_dp, t_2c_ri_2, &
1888  contract_1=[2], notcontract_1=[1], &
1889  contract_2=[1], notcontract_2=[2], &
1890  map_1=[1], map_2=[2], filter_eps=eps_filter, &
1891  flop=flop, unit_nr=unit_nr_dbcsr)
1892  dbcsr_nflop = dbcsr_nflop + flop
1893 
1894  !Here we do the trace for the force. In the periodic case, we store the matrix in G_PQ for later
1895  pref = 0.5_dp*fac
1896  IF (force_data%do_periodic) THEN
1897  CALL dbt_scale(t_2c_ri_2, pref)
1898  CALL dbt_create(force_data%G_PQ, t_2c_tmp)
1899  CALL dbt_copy(t_2c_ri_2, t_2c_tmp, move_data=.true.)
1900  CALL dbt_copy_tensor_to_matrix(t_2c_tmp, force_data%G_PQ, summation=.true.)
1901  CALL dbt_destroy(t_2c_tmp)
1902  ELSE
1903  CALL get_2c_der_force(force, t_2c_ri_2, force_data%t_2c_der_pot, atom_of_kind, &
1904  kind_of, force_data%idx_to_at_RI, pref, do_mp2=.true.)
1905 
1906  IF (use_virial) THEN
1907  CALL dbt_copy(t_2c_ri_2, t_2c_virial)
1908  CALL dbt_scale(t_2c_virial, pref)
1909  CALL dbt_copy_tensor_to_matrix(t_2c_virial, force_data%RI_virial_pot, summation=.true.)
1910  CALL dbt_clear(t_2c_virial)
1911  END IF
1912  END IF
1913 
1914  CALL dbt_clear(t_2c_ri)
1915  CALL dbt_clear(t_2c_ri_2)
1916 
1917  IF (use_virial) CALL dbt_destroy(t_2c_virial)
1918 
1919  CALL timestop(handle)
1920 
1921  END SUBROUTINE perform_2c_ops
1922 
1923 ! **************************************************************************************************
1924 !> \brief This subroutines performs the 3c tensor operations that are common accros low-scaling RPA
1925 !> and SOS-MP2, including forces and virial
1926 !> \param force ...
1927 !> \param t_R_occ ...
1928 !> \param t_R_virt ...
1929 !> \param force_data ...
1930 !> \param fac ...
1931 !> \param cut_memory ...
1932 !> \param n_mem_RI ...
1933 !> \param t_KBKT ...
1934 !> \param t_dm_occ ...
1935 !> \param t_dm_virt ...
1936 !> \param t_3c_O ...
1937 !> \param t_3c_M ...
1938 !> \param t_M_occ ...
1939 !> \param t_M_virt ...
1940 !> \param t_3c_0 ...
1941 !> \param t_3c_1 ...
1942 !> \param t_3c_3 ...
1943 !> \param t_3c_4 ...
1944 !> \param t_3c_5 ...
1945 !> \param t_3c_6 ...
1946 !> \param t_3c_7 ...
1947 !> \param t_3c_8 ...
1948 !> \param t_3c_sparse ...
1949 !> \param t_3c_help_1 ...
1950 !> \param t_3c_help_2 ...
1951 !> \param t_3c_ints ...
1952 !> \param t_3c_work ...
1953 !> \param starts_array_mc ...
1954 !> \param ends_array_mc ...
1955 !> \param batch_start_RI ...
1956 !> \param batch_end_RI ...
1957 !> \param t_3c_O_compressed ...
1958 !> \param t_3c_O_ind ...
1959 !> \param use_virial ...
1960 !> \param atom_of_kind ...
1961 !> \param kind_of ...
1962 !> \param eps_filter ...
1963 !> \param occ_ddint ...
1964 !> \param nze_ddint ...
1965 !> \param dbcsr_nflop ...
1966 !> \param unit_nr_dbcsr ...
1967 !> \param mp2_env ...
1968 ! **************************************************************************************************
1969  SUBROUTINE perform_3c_ops(force, t_R_occ, t_R_virt, force_data, fac, cut_memory, n_mem_RI, &
1970  t_KBKT, t_dm_occ, t_dm_virt, t_3c_O, t_3c_M, t_M_occ, t_M_virt, t_3c_0, t_3c_1, &
1971  t_3c_3, t_3c_4, t_3c_5, t_3c_6, t_3c_7, t_3c_8, t_3c_sparse, t_3c_help_1, t_3c_help_2, &
1972  t_3c_ints, t_3c_work, starts_array_mc, ends_array_mc, batch_start_RI, &
1973  batch_end_RI, t_3c_O_compressed, t_3c_O_ind, use_virial, &
1974  atom_of_kind, kind_of, eps_filter, occ_ddint, nze_ddint, dbcsr_nflop, &
1975  unit_nr_dbcsr, mp2_env)
1976 
1977  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1978  TYPE(dbt_type), INTENT(INOUT) :: t_r_occ, t_r_virt
1979  TYPE(im_time_force_type), INTENT(INOUT) :: force_data
1980  REAL(dp), INTENT(IN) :: fac
1981  INTEGER, INTENT(IN) :: cut_memory, n_mem_ri
1982  TYPE(dbt_type), INTENT(INOUT) :: t_kbkt, t_dm_occ, t_dm_virt, t_3c_o, t_3c_m, t_m_occ, &
1983  t_m_virt, t_3c_0, t_3c_1, t_3c_3, t_3c_4, t_3c_5, t_3c_6, t_3c_7, t_3c_8, t_3c_sparse, &
1984  t_3c_help_1, t_3c_help_2, t_3c_ints, t_3c_work
1985  INTEGER, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
1986  batch_start_ri, batch_end_ri
1987  TYPE(hfx_compression_type), DIMENSION(:) :: t_3c_o_compressed
1988  TYPE(block_ind_type), DIMENSION(:), INTENT(INOUT) :: t_3c_o_ind
1989  LOGICAL, INTENT(IN) :: use_virial
1990  INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind, kind_of
1991  REAL(dp), INTENT(IN) :: eps_filter
1992  REAL(dp), INTENT(INOUT) :: occ_ddint
1993  INTEGER(int_8), INTENT(INOUT) :: nze_ddint, dbcsr_nflop
1994  INTEGER, INTENT(IN) :: unit_nr_dbcsr
1995  TYPE(mp2_type) :: mp2_env
1996 
1997  CHARACTER(LEN=*), PARAMETER :: routinen = 'perform_3c_ops'
1998 
1999  INTEGER :: dummy_int, handle, handle2, i_mem, &
2000  i_xyz, j_mem, k_mem
2001  INTEGER(int_8) :: flop, nze
2002  INTEGER, DIMENSION(2, 1) :: ibounds, jbounds, kbounds
2003  INTEGER, DIMENSION(2, 2) :: bounds_2c
2004  INTEGER, DIMENSION(2, 3) :: bounds_cpy
2005  INTEGER, DIMENSION(3) :: bounds_3c
2006  REAL(dp) :: memory, occ, pref
2007  TYPE(block_ind_type), ALLOCATABLE, DIMENSION(:, :) :: blk_indices
2008  TYPE(hfx_compression_type), ALLOCATABLE, &
2009  DIMENSION(:, :) :: store_3c
2010 
2011  CALL timeset(routinen, handle)
2012 
2013  CALL dbt_get_info(t_3c_m, nfull_total=bounds_3c)
2014 
2015  !Pre-compute and compress KBK^T * (pq|R)
2016  ALLOCATE (store_3c(n_mem_ri, cut_memory))
2017  ALLOCATE (blk_indices(n_mem_ri, cut_memory))
2018  memory = 0.0_dp
2019  CALL timeset(routinen//"_pre_3c", handle2)
2020  !temporarily build the full int 3c tensor
2021  CALL dbt_copy(t_3c_o, t_3c_0)
2022  DO i_mem = 1, cut_memory
2023  CALL decompress_tensor(t_3c_o, t_3c_o_ind(i_mem)%ind, t_3c_o_compressed(i_mem), &
2024  mp2_env%ri_rpa_im_time%eps_compress)
2025  CALL dbt_copy(t_3c_o, t_3c_ints)
2026  CALL dbt_copy(t_3c_o, t_3c_0, move_data=.true., summation=.true.)
2027 
2028  DO k_mem = 1, n_mem_ri
2029  kbounds(:, 1) = [batch_start_ri(k_mem), batch_end_ri(k_mem)]
2030 
2031  CALL alloc_containers(store_3c(k_mem, i_mem), 1)
2032 
2033  !contract witht KBK^T over the RI index and store
2034  CALL dbt_batched_contract_init(t_kbkt)
2035  CALL dbt_contract(1.0_dp, t_kbkt, t_3c_ints, 0.0_dp, t_3c_work, &
2036  contract_1=[2], notcontract_1=[1], &
2037  contract_2=[1], notcontract_2=[2, 3], &
2038  map_1=[1], map_2=[2, 3], filter_eps=eps_filter, &
2039  bounds_2=kbounds, flop=flop, unit_nr=unit_nr_dbcsr)
2040  CALL dbt_batched_contract_finalize(t_kbkt)
2041  dbcsr_nflop = dbcsr_nflop + flop
2042 
2043  CALL dbt_copy(t_3c_work, t_3c_m, move_data=.true.)
2044  CALL compress_tensor(t_3c_m, blk_indices(k_mem, i_mem)%ind, store_3c(k_mem, i_mem), &
2045  mp2_env%ri_rpa_im_time%eps_compress, memory)
2046  END DO
2047  END DO !i_mem
2048  CALL dbt_clear(t_3c_m)
2049  CALL dbt_copy(t_3c_m, t_3c_ints)
2050  CALL timestop(handle2)
2051 
2052  CALL dbt_batched_contract_init(t_r_occ)
2053  CALL dbt_batched_contract_init(t_r_virt)
2054  DO i_mem = 1, cut_memory
2055  ibounds(:, 1) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
2056 
2057  !Compute the matrices M (integrals in t_3c_0)
2058  CALL timeset(routinen//"_3c_M", handle2)
2059  CALL dbt_batched_contract_init(t_dm_occ)
2060  CALL dbt_contract(1.0_dp, t_3c_0, t_dm_occ, 0.0_dp, t_3c_1, &
2061  contract_1=[3], notcontract_1=[1, 2], &
2062  contract_2=[1], notcontract_2=[2], &
2063  map_1=[1, 2], map_2=[3], filter_eps=eps_filter, &
2064  bounds_3=ibounds, flop=flop, unit_nr=unit_nr_dbcsr)
2065  dbcsr_nflop = dbcsr_nflop + flop
2066  CALL dbt_batched_contract_finalize(t_dm_occ)
2067  CALL dbt_copy(t_3c_1, t_m_occ, order=[1, 3, 2], move_data=.true.)
2068 
2069  CALL dbt_batched_contract_init(t_dm_virt)
2070  CALL dbt_contract(1.0_dp, t_3c_0, t_dm_virt, 0.0_dp, t_3c_1, &
2071  contract_1=[3], notcontract_1=[1, 2], &
2072  contract_2=[1], notcontract_2=[2], &
2073  map_1=[1, 2], map_2=[3], filter_eps=eps_filter, &
2074  bounds_3=ibounds, flop=flop, unit_nr=unit_nr_dbcsr)
2075  dbcsr_nflop = dbcsr_nflop + flop
2076  CALL dbt_batched_contract_finalize(t_dm_virt)
2077  CALL dbt_copy(t_3c_1, t_m_virt, order=[1, 3, 2], move_data=.true.)
2078  CALL timestop(handle2)
2079 
2080  !Compute the R matrices
2081  CALL timeset(routinen//"_3c_R", handle2)
2082  DO k_mem = 1, n_mem_ri
2083  CALL decompress_tensor(t_3c_m, blk_indices(k_mem, i_mem)%ind, store_3c(k_mem, i_mem), &
2084  mp2_env%ri_rpa_im_time%eps_compress)
2085  CALL dbt_copy(t_3c_m, t_3c_3, move_data=.true.)
2086 
2087  CALL dbt_contract(1.0_dp, t_m_occ, t_3c_3, 1.0_dp, t_r_occ, &
2088  contract_1=[1, 2], notcontract_1=[3], &
2089  contract_2=[1, 2], notcontract_2=[3], &
2090  map_1=[1], map_2=[2], filter_eps=eps_filter, &
2091  flop=flop, unit_nr=unit_nr_dbcsr)
2092  dbcsr_nflop = dbcsr_nflop + flop
2093 
2094  CALL dbt_contract(1.0_dp, t_m_virt, t_3c_3, 1.0_dp, t_r_virt, &
2095  contract_1=[1, 2], notcontract_1=[3], &
2096  contract_2=[1, 2], notcontract_2=[3], &
2097  map_1=[1], map_2=[2], filter_eps=eps_filter, &
2098  flop=flop, unit_nr=unit_nr_dbcsr)
2099  dbcsr_nflop = dbcsr_nflop + flop
2100  END DO
2101  CALL dbt_copy(t_3c_m, t_3c_3)
2102  CALL dbt_copy(t_3c_m, t_m_virt)
2103  CALL timestop(handle2)
2104 
2105  CALL dbt_copy(t_m_occ, t_3c_4, move_data=.true.)
2106 
2107  DO j_mem = 1, cut_memory
2108  jbounds(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
2109 
2110  bounds_cpy(:, 1) = [1, bounds_3c(1)]
2111  bounds_cpy(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
2112  bounds_cpy(:, 3) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
2113  CALL dbt_copy(t_3c_sparse, t_3c_7, bounds=bounds_cpy)
2114 
2115  CALL dbt_batched_contract_init(t_dm_virt)
2116  DO k_mem = 1, n_mem_ri
2117  bounds_2c(:, 1) = [batch_start_ri(k_mem), batch_end_ri(k_mem)]
2118  bounds_2c(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
2119 
2120  CALL timeset(routinen//"_3c_dm", handle2)
2121 
2122  !Calculate (mu nu| P) * D_occ * D_virt
2123  !Note: technically need M_occ*D_virt + M_virt*D_occ, but it is equivalent to 2*M_occ*D_virt
2124  CALL dbt_contract(2.0_dp, t_3c_4, t_dm_virt, 0.0_dp, t_3c_5, &
2125  contract_1=[3], notcontract_1=[1, 2], &
2126  contract_2=[1], notcontract_2=[2], &
2127  map_1=[1, 2], map_2=[3], filter_eps=eps_filter, &
2128  bounds_2=bounds_2c, bounds_3=jbounds, flop=flop, unit_nr=unit_nr_dbcsr)
2129  dbcsr_nflop = dbcsr_nflop + flop
2130 
2131  CALL get_tensor_occupancy(t_3c_5, nze, occ)
2132  nze_ddint = nze_ddint + nze
2133  occ_ddint = occ_ddint + occ
2134 
2135  CALL dbt_copy(t_3c_5, t_3c_6, move_data=.true.)
2136  CALL timestop(handle2)
2137 
2138  !Calculate the contraction of the above with K*B*K^T
2139  CALL timeset(routinen//"_3c_KBK", handle2)
2140  CALL dbt_batched_contract_init(t_kbkt)
2141  CALL dbt_contract(1.0_dp, t_kbkt, t_3c_6, 0.0_dp, t_3c_7, &
2142  contract_1=[2], notcontract_1=[1], &
2143  contract_2=[1], notcontract_2=[2, 3], &
2144  map_1=[1], map_2=[2, 3], &
2145  retain_sparsity=.true., flop=flop, unit_nr=unit_nr_dbcsr)
2146  dbcsr_nflop = dbcsr_nflop + flop
2147  CALL dbt_batched_contract_finalize(t_kbkt)
2148  CALL timestop(handle2)
2149  CALL dbt_copy(t_3c_7, t_3c_8, summation=.true.)
2150 
2151  END DO !k_mem
2152  CALL dbt_batched_contract_finalize(t_dm_virt)
2153  END DO !j_mem
2154 
2155  CALL dbt_copy(t_3c_8, t_3c_help_1, move_data=.true.)
2156 
2157  pref = 1.0_dp*fac
2158  DO k_mem = 1, cut_memory
2159  DO i_xyz = 1, 3
2160  CALL dbt_clear(force_data%t_3c_der_RI(i_xyz))
2161  CALL decompress_tensor(force_data%t_3c_der_RI(i_xyz), force_data%t_3c_der_RI_ind(k_mem, i_xyz)%ind, &
2162  force_data%t_3c_der_RI_comp(k_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
2163  END DO
2164  CALL get_force_from_3c_trace(force, t_3c_help_1, force_data%t_3c_der_RI, atom_of_kind, kind_of, &
2165  force_data%idx_to_at_RI, pref, do_mp2=.true., deriv_dim=1)
2166  END DO
2167 
2168  IF (use_virial) THEN
2169  CALL dbt_copy(t_3c_help_1, t_3c_help_2)
2170  CALL dbt_scale(t_3c_help_2, pref)
2171  CALL dbt_copy(t_3c_help_2, force_data%t_3c_virial_split, summation=.true., move_data=.true.)
2172  END IF
2173 
2174  CALL dbt_copy(t_3c_help_1, t_3c_help_2)
2175  CALL dbt_copy(t_3c_help_1, t_3c_help_2, order=[1, 3, 2], move_data=.true., summation=.true.)
2176  DO k_mem = 1, cut_memory
2177  DO i_xyz = 1, 3
2178  CALL dbt_clear(force_data%t_3c_der_AO(i_xyz))
2179  CALL decompress_tensor(force_data%t_3c_der_AO(i_xyz), force_data%t_3c_der_AO_ind(k_mem, i_xyz)%ind, &
2180  force_data%t_3c_der_AO_comp(k_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
2181  END DO
2182  CALL get_force_from_3c_trace(force, t_3c_help_2, force_data%t_3c_der_AO, atom_of_kind, kind_of, &
2183  force_data%idx_to_at_AO, pref, do_mp2=.true., deriv_dim=3)
2184  END DO
2185 
2186  CALL dbt_clear(t_3c_help_2)
2187  END DO !i_mem
2188  CALL dbt_batched_contract_finalize(t_r_occ)
2189  CALL dbt_batched_contract_finalize(t_r_virt)
2190 
2191  DO k_mem = 1, n_mem_ri
2192  DO i_mem = 1, cut_memory
2193  CALL dealloc_containers(store_3c(k_mem, i_mem), dummy_int)
2194  END DO
2195  END DO
2196  DEALLOCATE (store_3c, blk_indices)
2197 
2198  CALL timestop(handle)
2199 
2200  END SUBROUTINE perform_3c_ops
2201 
2202 ! **************************************************************************************************
2203 !> \brief All the forces that can be calculated after the loop on the Laplace quaradture, using
2204 !> terms collected during the said loop. This inludes the z-vector equation and its reponse
2205 !> forces, as well as the force coming from the trace with the derivative of the KS matrix
2206 !> \param force_data ...
2207 !> \param unit_nr ...
2208 !> \param qs_env ...
2209 ! **************************************************************************************************
2210  SUBROUTINE calc_post_loop_forces(force_data, unit_nr, qs_env)
2211 
2212  TYPE(im_time_force_type), INTENT(INOUT) :: force_data
2213  INTEGER, INTENT(IN) :: unit_nr
2214  TYPE(qs_environment_type), POINTER :: qs_env
2215 
2216  CHARACTER(len=*), PARAMETER :: routinen = 'calc_post_loop_forces'
2217 
2218  INTEGER :: handle, ispin, nao, nao_aux, nocc, nspins
2219  LOGICAL :: do_exx
2220  REAL(dp) :: focc
2221  TYPE(admm_type), POINTER :: admm_env
2222  TYPE(cp_fm_struct_type), POINTER :: fm_struct
2223  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: cpmos, mo_occ
2224  TYPE(cp_fm_type), POINTER :: mo_coeff
2225  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_p_work, matrix_p_mp2, &
2226  matrix_p_mp2_admm, matrix_s, &
2227  matrix_s_aux, work_admm, yp_admm
2228  TYPE(dft_control_type), POINTER :: dft_control
2229  TYPE(linres_control_type), POINTER :: linres_control
2230  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2231  TYPE(qs_p_env_type), POINTER :: p_env
2232  TYPE(section_vals_type), POINTER :: hfx_section, lr_section
2233 
2234  NULLIFY (linres_control, p_env, dft_control, matrix_s, mos, mo_coeff, fm_struct, lr_section, &
2235  dbcsr_p_work, yp_admm, matrix_p_mp2, admm_env, work_admm, matrix_s_aux, matrix_p_mp2_admm)
2236 
2237  CALL timeset(routinen, handle)
2238 
2239  CALL get_qs_env(qs_env, dft_control=dft_control, matrix_s=matrix_s, mos=mos)
2240  nspins = dft_control%nspins
2241 
2242  ! Setting up for the z-vector equation
2243 
2244  ! Initialize linres_control
2245  lr_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%LOW_SCALING%CPHF")
2246 
2247  ALLOCATE (linres_control)
2248  CALL section_vals_val_get(lr_section, "MAX_ITER", i_val=linres_control%max_iter)
2249  CALL section_vals_val_get(lr_section, "EPS_CONV", r_val=linres_control%eps)
2250  CALL section_vals_val_get(lr_section, "PRECONDITIONER", i_val=linres_control%preconditioner_type)
2251  CALL section_vals_val_get(lr_section, "ENERGY_GAP", r_val=linres_control%energy_gap)
2252 
2253  linres_control%do_kernel = .true.
2254  linres_control%lr_triplet = .false.
2255  linres_control%converged = .false.
2256  linres_control%eps_filter = qs_env%mp2_env%ri_rpa_im_time%eps_filter
2257 
2258  CALL set_qs_env(qs_env, linres_control=linres_control)
2259 
2260  IF (unit_nr > 0) THEN
2261  WRITE (unit_nr, *)
2262  WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations'
2263  WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Convergence threshold:', linres_control%eps
2264  WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Maximum number of iterations: ', linres_control%max_iter
2265  END IF
2266 
2267  ALLOCATE (p_env)
2268  CALL p_env_create(p_env, qs_env, orthogonal_orbitals=.true., linres_control=linres_control)
2269  CALL p_env_psi0_changed(p_env, qs_env)
2270 
2271  ! Matrix allocation
2272  CALL dbcsr_allocate_matrix_set(p_env%p1, nspins)
2273  CALL dbcsr_allocate_matrix_set(p_env%w1, nspins)
2274  CALL dbcsr_allocate_matrix_set(dbcsr_p_work, nspins)
2275  DO ispin = 1, nspins
2276  ALLOCATE (p_env%p1(ispin)%matrix, p_env%w1(ispin)%matrix, dbcsr_p_work(ispin)%matrix)
2277  CALL dbcsr_create(matrix=p_env%p1(ispin)%matrix, template=matrix_s(1)%matrix)
2278  CALL dbcsr_create(matrix=p_env%w1(ispin)%matrix, template=matrix_s(1)%matrix)
2279  CALL dbcsr_create(matrix=dbcsr_p_work(ispin)%matrix, template=matrix_s(1)%matrix)
2280  CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_s(1)%matrix)
2281  CALL dbcsr_copy(p_env%w1(ispin)%matrix, matrix_s(1)%matrix)
2282  CALL dbcsr_copy(dbcsr_p_work(ispin)%matrix, matrix_s(1)%matrix)
2283  CALL dbcsr_set(p_env%p1(ispin)%matrix, 0.0_dp)
2284  CALL dbcsr_set(p_env%w1(ispin)%matrix, 0.0_dp)
2285  CALL dbcsr_set(dbcsr_p_work(ispin)%matrix, 0.0_dp)
2286  END DO
2287 
2288  IF (dft_control%do_admm) THEN
2289  CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux)
2290  CALL dbcsr_allocate_matrix_set(p_env%p1_admm, nspins)
2291  CALL dbcsr_allocate_matrix_set(work_admm, nspins)
2292  DO ispin = 1, nspins
2293  ALLOCATE (p_env%p1_admm(ispin)%matrix, work_admm(ispin)%matrix)
2294  CALL dbcsr_create(p_env%p1_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2295  CALL dbcsr_copy(p_env%p1_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2296  CALL dbcsr_set(p_env%p1_admm(ispin)%matrix, 0.0_dp)
2297  CALL dbcsr_create(work_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2298  CALL dbcsr_copy(work_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2299  CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp)
2300  END DO
2301  END IF
2302 
2303  ! Preparing the RHS of the z-vector equation
2304  CALL prepare_for_response(force_data, qs_env)
2305  ALLOCATE (cpmos(nspins), mo_occ(nspins))
2306  DO ispin = 1, nspins
2307  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, homo=nocc)
2308  NULLIFY (fm_struct)
2309  CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, &
2310  template_fmstruct=mo_coeff%matrix_struct)
2311  CALL cp_fm_create(cpmos(ispin), fm_struct)
2312  CALL cp_fm_set_all(cpmos(ispin), 0.0_dp)
2313  CALL cp_fm_create(mo_occ(ispin), fm_struct)
2314  CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc)
2315  CALL cp_fm_struct_release(fm_struct)
2316  END DO
2317 
2318  ! in case of EXX, need to add the HF Hamiltonian to the RHS of the Z-vector equation
2319  ! Strategy: we take the ks_matrix, remove the current xc contribution, and then add the RPA HF one
2320  do_exx = .false.
2321  IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
2322  hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
2323  CALL section_vals_get(hfx_section, explicit=do_exx)
2324  END IF
2325 
2326  IF (do_exx) THEN
2327  CALL add_exx_to_rhs(rhs=force_data%sum_O_tau, &
2328  qs_env=qs_env, &
2329  ext_hfx_section=hfx_section, &
2330  x_data=qs_env%mp2_env%ri_rpa%x_data, &
2331  recalc_integrals=.false., &
2332  do_admm=qs_env%mp2_env%ri_rpa%do_admm, &
2333  do_exx=do_exx, &
2334  reuse_hfx=qs_env%mp2_env%ri_rpa%reuse_hfx)
2335  END IF
2336 
2337  focc = 2.0_dp
2338  IF (nspins == 1) focc = 4.0_dp
2339  DO ispin = 1, nspins
2340  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
2341  CALL cp_dbcsr_sm_fm_multiply(force_data%sum_O_tau(ispin)%matrix, mo_occ(ispin), &
2342  cpmos(ispin), nocc, &
2343  alpha=focc, beta=0.0_dp)
2344  END DO
2345 
2346  ! The z-vector equation and associated forces
2347  CALL response_equation_new(qs_env, p_env, cpmos, unit_nr)
2348 
2349  ! Save the mp2 density matrix
2350  CALL get_qs_env(qs_env, matrix_p_mp2=matrix_p_mp2)
2351  IF (ASSOCIATED(matrix_p_mp2)) CALL dbcsr_deallocate_matrix_set(matrix_p_mp2)
2352  DO ispin = 1, nspins
2353  CALL dbcsr_copy(dbcsr_p_work(ispin)%matrix, p_env%p1(ispin)%matrix)
2354  CALL dbcsr_add(dbcsr_p_work(ispin)%matrix, force_data%sum_YP_tau(ispin)%matrix, 1.0_dp, 1.0_dp)
2355  END DO
2356  CALL set_ks_env(qs_env%ks_env, matrix_p_mp2=dbcsr_p_work)
2357 
2358  IF (dft_control%do_admm) THEN
2359  CALL dbcsr_allocate_matrix_set(yp_admm, nspins)
2360  CALL get_qs_env(qs_env, matrix_p_mp2_admm=matrix_p_mp2_admm, admm_env=admm_env)
2361  nao = admm_env%nao_orb
2362  nao_aux = admm_env%nao_aux_fit
2363  IF (ASSOCIATED(matrix_p_mp2_admm)) CALL dbcsr_deallocate_matrix_set(matrix_p_mp2_admm)
2364  DO ispin = 1, nspins
2365 
2366  !sum_YP_tau in the auxiliary basis
2367  CALL copy_dbcsr_to_fm(force_data%sum_YP_tau(ispin)%matrix, admm_env%work_orb_orb)
2368  CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, admm_env%work_orb_orb, &
2369  0.0_dp, admm_env%work_aux_orb)
2370  CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
2371  0.0_dp, admm_env%work_aux_aux)
2372  CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_admm(ispin)%matrix, keep_sparsity=.true.)
2373 
2374  !save the admm representation od sum_YP_tau
2375  ALLOCATE (yp_admm(ispin)%matrix)
2376  CALL dbcsr_create(yp_admm(ispin)%matrix, template=work_admm(ispin)%matrix)
2377  CALL dbcsr_copy(yp_admm(ispin)%matrix, work_admm(ispin)%matrix)
2378 
2379  CALL dbcsr_add(work_admm(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
2380 
2381  END DO
2382  CALL set_ks_env(qs_env%ks_env, matrix_p_mp2_admm=work_admm)
2383  END IF
2384 
2385  !Calculate the response force and the force from the trace with F
2386  CALL update_im_time_forces(p_env, force_data%sum_O_tau, force_data%sum_YP_tau, yp_admm, qs_env)
2387 
2388  !clean-up
2389  IF (dft_control%do_admm) CALL dbcsr_deallocate_matrix_set(yp_admm)
2390 
2391  CALL cp_fm_release(cpmos)
2392  CALL cp_fm_release(mo_occ)
2393  CALL p_env_release(p_env)
2394  DEALLOCATE (p_env)
2395 
2396  CALL timestop(handle)
2397 
2398  END SUBROUTINE calc_post_loop_forces
2399 
2400 ! **************************************************************************************************
2401 !> \brief Prepares the RHS of the z-vector equation. Apply the xc and HFX kernel on the previously
2402 !> stored sum_YP_tau density, and add it to the final force_data%sum_O_tau quantity
2403 !> \param force_data ...
2404 !> \param qs_env ...
2405 ! **************************************************************************************************
2406  SUBROUTINE prepare_for_response(force_data, qs_env)
2407 
2408  TYPE(im_time_force_type), INTENT(INOUT) :: force_data
2409  TYPE(qs_environment_type), POINTER :: qs_env
2410 
2411  CHARACTER(LEN=*), PARAMETER :: routinen = 'prepare_for_response'
2412 
2413  INTEGER :: handle, ispin, nao, nao_aux, nspins
2414  LOGICAL :: do_hfx, do_tau, do_tau_admm
2415  REAL(dp) :: ehartree
2416  TYPE(admm_type), POINTER :: admm_env
2417  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_p_work, ker_tau_admm, matrix_s, &
2418  matrix_s_aux, work_admm
2419  TYPE(dbcsr_type) :: dbcsr_work
2420  TYPE(dft_control_type), POINTER :: dft_control
2421  TYPE(pw_c1d_gs_type) :: rhoz_tot_gspace, zv_hartree_gspace
2422  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rhoz_g
2423  TYPE(pw_env_type), POINTER :: pw_env
2424  TYPE(pw_poisson_type), POINTER :: poisson_env
2425  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2426  TYPE(pw_r3d_rs_type) :: zv_hartree_rspace
2427  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rhoz_r, tauz_r, v_xc, v_xc_tau
2428  TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
2429  TYPE(section_vals_type), POINTER :: hfx_section, xc_section
2430  TYPE(task_list_type), POINTER :: task_list_aux_fit
2431 
2432  NULLIFY (pw_env, rhoz_r, rhoz_g, tauz_r, v_xc, v_xc_tau, &
2433  poisson_env, auxbas_pw_pool, dft_control, admm_env, xc_section, rho, rho_aux_fit, &
2434  task_list_aux_fit, ker_tau_admm, work_admm, dbcsr_p_work, matrix_s, hfx_section)
2435 
2436  CALL timeset(routinen, handle)
2437 
2438  CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env, rho=rho, matrix_s=matrix_s)
2439  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
2440  nspins = dft_control%nspins
2441 
2442  CALL dbcsr_allocate_matrix_set(dbcsr_p_work, nspins)
2443  DO ispin = 1, nspins
2444  ALLOCATE (dbcsr_p_work(ispin)%matrix)
2445  CALL dbcsr_create(matrix=dbcsr_p_work(ispin)%matrix, template=matrix_s(1)%matrix)
2446  CALL dbcsr_copy(dbcsr_p_work(ispin)%matrix, matrix_s(1)%matrix)
2447  CALL dbcsr_set(dbcsr_p_work(ispin)%matrix, 0.0_dp)
2448  END DO
2449 
2450  !Apply the kernel on the density saved in force_data%sum_YP_tau
2451  ALLOCATE (rhoz_r(nspins), rhoz_g(nspins))
2452  DO ispin = 1, nspins
2453  CALL auxbas_pw_pool%create_pw(rhoz_r(ispin))
2454  CALL auxbas_pw_pool%create_pw(rhoz_g(ispin))
2455  END DO
2456  CALL auxbas_pw_pool%create_pw(rhoz_tot_gspace)
2457  CALL auxbas_pw_pool%create_pw(zv_hartree_rspace)
2458  CALL auxbas_pw_pool%create_pw(zv_hartree_gspace)
2459 
2460  CALL pw_zero(rhoz_tot_gspace)
2461  DO ispin = 1, nspins
2462  CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=force_data%sum_YP_tau(ispin)%matrix, &
2463  rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin))
2464  CALL pw_axpy(rhoz_g(ispin), rhoz_tot_gspace)
2465  END DO
2466 
2467  CALL pw_poisson_solve(poisson_env, rhoz_tot_gspace, ehartree, &
2468  zv_hartree_gspace)
2469 
2470  CALL pw_transfer(zv_hartree_gspace, zv_hartree_rspace)
2471  CALL pw_scale(zv_hartree_rspace, zv_hartree_rspace%pw_grid%dvol)
2472 
2473  CALL qs_rho_get(rho, tau_r_valid=do_tau)
2474  IF (do_tau) THEN
2475  block
2476  TYPE(pw_c1d_gs_type) :: tauz_g
2477  ALLOCATE (tauz_r(nspins))
2478  CALL auxbas_pw_pool%create_pw(tauz_g)
2479  DO ispin = 1, nspins
2480  CALL auxbas_pw_pool%create_pw(tauz_r(ispin))
2481 
2482  CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=force_data%sum_YP_tau(ispin)%matrix, &
2483  rho=tauz_r(ispin), rho_gspace=tauz_g, compute_tau=.true.)
2484  END DO
2485  CALL auxbas_pw_pool%give_back_pw(tauz_g)
2486  END block
2487  END IF
2488 
2489  IF (dft_control%do_admm) THEN
2490  CALL get_qs_env(qs_env, admm_env=admm_env)
2491  xc_section => admm_env%xc_section_primary
2492  ELSE
2493  xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
2494  END IF
2495 
2496  !Primary XC kernel
2497  CALL create_kernel(qs_env, v_xc, v_xc_tau, rho, rhoz_r, rhoz_g, tauz_r, xc_section)
2498 
2499  DO ispin = 1, nspins
2500  CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2501  CALL pw_axpy(zv_hartree_rspace, v_xc(ispin))
2502  CALL integrate_v_rspace(qs_env=qs_env, &
2503  v_rspace=v_xc(ispin), &
2504  hmat=dbcsr_p_work(ispin), &
2505  calculate_forces=.false.)
2506 
2507  CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2508  END DO
2509  CALL auxbas_pw_pool%give_back_pw(rhoz_tot_gspace)
2510  CALL auxbas_pw_pool%give_back_pw(zv_hartree_rspace)
2511  CALL auxbas_pw_pool%give_back_pw(zv_hartree_gspace)
2512  DEALLOCATE (v_xc)
2513 
2514  IF (do_tau) THEN
2515  DO ispin = 1, nspins
2516  CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2517  CALL integrate_v_rspace(qs_env=qs_env, &
2518  v_rspace=v_xc_tau(ispin), &
2519  hmat=dbcsr_p_work(ispin), &
2520  compute_tau=.true., &
2521  calculate_forces=.false.)
2522  CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2523  END DO
2524  DEALLOCATE (v_xc_tau)
2525  END IF
2526 
2527  !Auxiliary xc kernel (admm)
2528  IF (dft_control%do_admm) THEN
2529  CALL get_qs_env(qs_env, admm_env=admm_env)
2530  CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux, &
2531  task_list_aux_fit=task_list_aux_fit, rho_aux_fit=rho_aux_fit)
2532 
2533  CALL qs_rho_get(rho_aux_fit, tau_r_valid=do_tau_admm)
2534 
2535  CALL dbcsr_allocate_matrix_set(work_admm, nspins)
2536  CALL dbcsr_allocate_matrix_set(ker_tau_admm, nspins)
2537  DO ispin = 1, nspins
2538  ALLOCATE (work_admm(ispin)%matrix, ker_tau_admm(ispin)%matrix)
2539  CALL dbcsr_create(work_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2540  CALL dbcsr_copy(work_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2541  CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp)
2542  CALL dbcsr_create(ker_tau_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2543  CALL dbcsr_copy(ker_tau_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2544  CALL dbcsr_set(ker_tau_admm(ispin)%matrix, 0.0_dp)
2545  END DO
2546 
2547  !get the density in the auxuliary density
2548  cpassert(ASSOCIATED(admm_env%work_orb_orb))
2549  cpassert(ASSOCIATED(admm_env%work_aux_orb))
2550  cpassert(ASSOCIATED(admm_env%work_aux_aux))
2551  nao = admm_env%nao_orb
2552  nao_aux = admm_env%nao_aux_fit
2553  DO ispin = 1, nspins
2554  CALL copy_dbcsr_to_fm(force_data%sum_YP_tau(ispin)%matrix, admm_env%work_orb_orb)
2555  CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, admm_env%work_orb_orb, &
2556  0.0_dp, admm_env%work_aux_orb)
2557  CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
2558  0.0_dp, admm_env%work_aux_aux)
2559  CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, ker_tau_admm(ispin)%matrix, keep_sparsity=.true.)
2560  END DO
2561 
2562  IF (.NOT. qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
2563  DO ispin = 1, nspins
2564  CALL pw_zero(rhoz_r(ispin))
2565  CALL pw_zero(rhoz_g(ispin))
2566  CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=ker_tau_admm(ispin)%matrix, &
2567  rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin), &
2568  basis_type="AUX_FIT", task_list_external=task_list_aux_fit)
2569  END DO
2570 
2571  IF (do_tau_admm) THEN
2572  block
2573  TYPE(pw_c1d_gs_type) :: tauz_g
2574  CALL auxbas_pw_pool%create_pw(tauz_g)
2575  DO ispin = 1, nspins
2576  CALL pw_zero(tauz_r(ispin))
2577  CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=ker_tau_admm(ispin)%matrix, &
2578  rho=tauz_r(ispin), rho_gspace=tauz_g, &
2579  basis_type="AUX_FIT", task_list_external=task_list_aux_fit, &
2580  compute_tau=.true.)
2581  END DO
2582  CALL auxbas_pw_pool%give_back_pw(tauz_g)
2583  END block
2584  END IF
2585 
2586  xc_section => admm_env%xc_section_aux
2587  CALL create_kernel(qs_env, v_xc, v_xc_tau, rho_aux_fit, rhoz_r, rhoz_g, tauz_r, xc_section)
2588 
2589  DO ispin = 1, nspins
2590  CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2591  CALL integrate_v_rspace(qs_env=qs_env, &
2592  v_rspace=v_xc(ispin), &
2593  hmat=work_admm(ispin), &
2594  calculate_forces=.false., &
2595  basis_type="AUX_FIT", &
2596  task_list_external=task_list_aux_fit)
2597  CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2598  END DO
2599  DEALLOCATE (v_xc)
2600 
2601  IF (do_tau_admm) THEN
2602  DO ispin = 1, nspins
2603  CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2604  CALL integrate_v_rspace(qs_env=qs_env, &
2605  v_rspace=v_xc_tau(ispin), &
2606  hmat=work_admm(ispin), &
2607  calculate_forces=.false., &
2608  basis_type="AUX_FIT", &
2609  task_list_external=task_list_aux_fit, &
2610  compute_tau=.true.)
2611  CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2612  END DO
2613  DEALLOCATE (v_xc_tau)
2614  END IF
2615  END IF !admm
2616  END IF
2617 
2618  DO ispin = 1, nspins
2619  CALL auxbas_pw_pool%give_back_pw(rhoz_r(ispin))
2620  CALL auxbas_pw_pool%give_back_pw(rhoz_g(ispin))
2621  END DO
2622  DEALLOCATE (rhoz_r, rhoz_g)
2623 
2624  IF (do_tau) THEN
2625  DO ispin = 1, nspins
2626  CALL auxbas_pw_pool%give_back_pw(tauz_r(ispin))
2627  END DO
2628  DEALLOCATE (tauz_r)
2629  END IF
2630 
2631  !HFX kernel
2632  hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
2633  CALL section_vals_get(hfx_section, explicit=do_hfx)
2634  IF (do_hfx) THEN
2635  IF (dft_control%do_admm) THEN
2636  CALL tddft_hfx_matrix(work_admm, ker_tau_admm, qs_env, .false., .false.)
2637 
2638  !Going back to primary basis
2639  CALL dbcsr_create(dbcsr_work, template=dbcsr_p_work(1)%matrix)
2640  CALL dbcsr_copy(dbcsr_work, dbcsr_p_work(1)%matrix)
2641  CALL dbcsr_set(dbcsr_work, 0.0_dp)
2642  DO ispin = 1, nspins
2643  CALL copy_dbcsr_to_fm(work_admm(ispin)%matrix, admm_env%work_aux_aux)
2644  CALL parallel_gemm('N', 'N', nao_aux, nao, nao_aux, 1.0_dp, admm_env%work_aux_aux, admm_env%A, &
2645  0.0_dp, admm_env%work_aux_orb)
2646  CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, admm_env%work_aux_orb, &
2647  0.0_dp, admm_env%work_orb_orb)
2648  CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_work, keep_sparsity=.true.)
2649  CALL dbcsr_add(dbcsr_p_work(ispin)%matrix, dbcsr_work, 1.0_dp, 1.0_dp)
2650  END DO
2651  CALL dbcsr_release(dbcsr_work)
2652  CALL dbcsr_deallocate_matrix_set(ker_tau_admm)
2653  ELSE
2654  CALL tddft_hfx_matrix(dbcsr_p_work, force_data%sum_YP_tau, qs_env, .false., .false.)
2655  END IF
2656  END IF
2657 
2658  DO ispin = 1, nspins
2659  CALL dbcsr_add(force_data%sum_O_tau(ispin)%matrix, dbcsr_p_work(ispin)%matrix, 1.0_dp, 1.0_dp)
2660  END DO
2661 
2662  CALL dbcsr_deallocate_matrix_set(dbcsr_p_work)
2663  CALL dbcsr_deallocate_matrix_set(work_admm)
2664 
2665  CALL timestop(handle)
2666 
2667  END SUBROUTINE prepare_for_response
2668 
2669 ! **************************************************************************************************
2670 !> \brief Calculate the force and virial due to the (P|Q) GPW integral derivatives
2671 !> \param G_PQ ...
2672 !> \param force ...
2673 !> \param h_stress ...
2674 !> \param use_virial ...
2675 !> \param mp2_env ...
2676 !> \param qs_env ...
2677 ! **************************************************************************************************
2678  SUBROUTINE get_2c_gpw_forces(G_PQ, force, h_stress, use_virial, mp2_env, qs_env)
2679 
2680  TYPE(dbcsr_type), INTENT(INOUT) :: g_pq
2681  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2682  REAL(dp), DIMENSION(3, 3), INTENT(INOUT) :: h_stress
2683  LOGICAL, INTENT(IN) :: use_virial
2684  TYPE(mp2_type), INTENT(INOUT) :: mp2_env
2685  TYPE(qs_environment_type), POINTER :: qs_env
2686 
2687  CHARACTER(len=*), PARAMETER :: routinen = 'get_2c_gpw_forces'
2688 
2689  INTEGER :: atom_a, color, handle, i, i_ri, i_xyz, iatom, igrid_level, ikind, ipgf, iset, j, &
2690  j_ri, jatom, lb_ri, n_ri, natom, ncoa, ncoms, nkind, nproc, nseta, o1, offset, pdims(2), &
2691  sgfa, ub_ri
2692  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, iproc_map, kind_of, &
2693  sizes_ri
2694  INTEGER, DIMENSION(:), POINTER :: col_dist, la_max, la_min, npgfa, nsgfa, &
2695  row_dist
2696  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, pgrid
2697  LOGICAL :: found, one_proc_group
2698  REAL(dp) :: cutoff_old, radius, relative_cutoff_old
2699  REAL(dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old, wf_vector
2700  REAL(dp), DIMENSION(3) :: force_a, force_b, ra
2701  REAL(dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
2702  REAL(kind=dp), DIMENSION(:, :), POINTER :: h_tmp, i_ab, pab, pblock, sphi_a, zeta
2703  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2704  TYPE(cell_type), POINTER :: cell
2705  TYPE(dbcsr_distribution_type) :: dbcsr_dist
2706  TYPE(dbcsr_type) :: tmp_g_pq
2707  TYPE(dft_control_type), POINTER :: dft_control
2708  TYPE(gto_basis_set_p_type), ALLOCATABLE, &
2709  DIMENSION(:), TARGET :: basis_set_ri_aux
2710  TYPE(gto_basis_set_type), POINTER :: basis_set_a
2711  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2712  TYPE(mp_para_env_type), POINTER :: para_env, para_env_ext
2713  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2714  POINTER :: sab_orb
2715  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2716  TYPE(pw_c1d_gs_type) :: dvg(3), pot_g, rho_g, rho_g_copy
2717  TYPE(pw_env_type), POINTER :: pw_env_ext
2718  TYPE(pw_poisson_type), POINTER :: poisson_env
2719  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2720  TYPE(pw_r3d_rs_type) :: psi_l, rho_r
2721  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2722  TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
2723  TYPE(task_list_type), POINTER :: task_list_ext
2724 
2725  NULLIFY (sab_orb, task_list_ext, particle_set, qs_kind_set, dft_control, pw_env_ext, auxbas_pw_pool, &
2726  poisson_env, atomic_kind_set, para_env, cell, rs_v, mos, basis_set_a)
2727 
2728  CALL timeset(routinen, handle)
2729 
2730  CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env, sab_orb=sab_orb, &
2731  natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, particle_set=particle_set, &
2732  mos=mos, cell=cell, atomic_kind_set=atomic_kind_set)
2733 
2734  !The idea is to use GPW to compute the integrals and derivatives. Because the potential needs
2735  !to be calculated for each phi_j (column) of all AO pairs, and because that is expensive, we want
2736  !to minimize the amount of time we do that. Therefore, we work with a special distribution, where
2737  !each column of the resulting DBCSR matrix is mapped to a sub-communicator.
2738 
2739  !Try to get the optimal pdims (we want a grid that is flat: many cols, few rows)
2740  IF (para_env%num_pe <= natom) THEN
2741  pdims(1) = 1
2742  pdims(2) = para_env%num_pe
2743  ELSE
2744  DO i = natom, 1, -1
2745  IF (modulo(para_env%num_pe, i) == 0) THEN
2746  pdims(1) = para_env%num_pe/i
2747  pdims(2) = i
2748  EXIT
2749  END IF
2750  END DO
2751  END IF
2752 
2753  ALLOCATE (row_dist(natom), col_dist(natom))
2754  DO iatom = 1, natom
2755  row_dist(iatom) = modulo(iatom, pdims(1))
2756  END DO
2757  DO jatom = 1, natom
2758  col_dist(jatom) = modulo(jatom, pdims(2))
2759  END DO
2760 
2761  ALLOCATE (pgrid(0:pdims(1) - 1, 0:pdims(2) - 1))
2762  nproc = 0
2763  DO i = 0, pdims(1) - 1
2764  DO j = 0, pdims(2) - 1
2765  pgrid(i, j) = nproc
2766  nproc = nproc + 1
2767  END DO
2768  END DO
2769 
2770  CALL dbcsr_distribution_new(dbcsr_dist, group=para_env%get_handle(), pgrid=pgrid, row_dist=row_dist, col_dist=col_dist)
2771 
2772  !The temporary DBCSR integrals and derivatives matrices in this flat distribution
2773  CALL dbcsr_create(tmp_g_pq, template=g_pq, matrix_type=dbcsr_type_no_symmetry, dist=dbcsr_dist)
2774  CALL dbcsr_complete_redistribute(g_pq, tmp_g_pq)
2775 
2776  ALLOCATE (basis_set_ri_aux(nkind), sizes_ri(natom))
2777  CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
2778  CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ri, basis=basis_set_ri_aux)
2779  n_ri = sum(sizes_ri)
2780 
2781  one_proc_group = mp2_env%mp2_num_proc == 1
2782  ALLOCATE (para_env_ext)
2783  IF (one_proc_group) THEN
2784  !one subgroup per proc
2785  CALL para_env_ext%from_split(para_env, para_env%mepos)
2786  ELSE
2787  !Split the communicator accross the columns of the matrix
2788  ncoms = min(pdims(2), para_env%num_pe/mp2_env%mp2_num_proc)
2789  DO i = 0, pdims(1) - 1
2790  DO j = 0, pdims(2) - 1
2791  IF (pgrid(i, j) == para_env%mepos) color = modulo(j + 1, ncoms)
2792  END DO
2793  END DO
2794  CALL para_env_ext%from_split(para_env, color)
2795  END IF
2796 
2797  !sab_orb and task_list_ext are essentially dummies
2798  CALL prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_ext, pw_env_ext, &
2799  auxbas_pw_pool, poisson_env, task_list_ext, rho_r, rho_g, pot_g, psi_l, sab_orb)
2800 
2801  IF (use_virial) THEN
2802  CALL auxbas_pw_pool%create_pw(rho_g_copy)
2803  DO i_xyz = 1, 3
2804  CALL auxbas_pw_pool%create_pw(dvg(i_xyz))
2805  END DO
2806  END IF
2807 
2808  ALLOCATE (wf_vector(n_ri))
2809 
2810  CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
2811 
2812  ALLOCATE (iproc_map(natom))
2813 
2814  !Loop over the atomic blocks
2815  DO jatom = 1, natom
2816 
2817  !Only calculate if on the correct sub-communicator/proc
2818  IF (one_proc_group) THEN
2819  iproc_map = 0
2820  DO iatom = 1, natom
2821  IF (pgrid(row_dist(iatom), col_dist(jatom)) == para_env%mepos) iproc_map(iatom) = 1
2822  END DO
2823  IF (.NOT. any(iproc_map == 1)) cycle
2824  ELSE
2825  IF (.NOT. modulo(col_dist(jatom) + 1, ncoms) == color) cycle
2826  END IF
2827 
2828  lb_ri = sum(sizes_ri(1:jatom - 1))
2829  ub_ri = lb_ri + sizes_ri(jatom)
2830  DO j_ri = lb_ri + 1, ub_ri
2831 
2832  wf_vector = 0.0_dp
2833  wf_vector(j_ri) = 1.0_dp
2834 
2835  CALL calculate_wavefunction(mos(1)%mo_coeff, 1, psi_l, rho_g, atomic_kind_set, &
2836  qs_kind_set, cell, dft_control, particle_set, pw_env_ext, &
2837  basis_type="RI_AUX", external_vector=wf_vector)
2838 
2839  IF (use_virial) THEN
2840  CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, mp2_env%potential_parameter, dvg)
2841 
2842  wf_vector = 0.0_dp
2843  DO iatom = 1, natom
2844  !only compute if i,j atom pair on correct proc
2845  IF (one_proc_group) THEN
2846  IF (.NOT. iproc_map(iatom) == 1) cycle
2847  END IF
2848 
2849  CALL dbcsr_get_block_p(tmp_g_pq, iatom, jatom, pblock, found)
2850  IF (.NOT. found) cycle
2851 
2852  i_ri = sum(sizes_ri(1:iatom - 1))
2853  wf_vector(i_ri + 1:i_ri + sizes_ri(iatom)) = pblock(:, j_ri - lb_ri)
2854  END DO
2855 
2856  CALL pw_copy(rho_g, rho_g_copy)
2857  CALL calculate_wavefunction(mos(1)%mo_coeff, 1, psi_l, rho_g, atomic_kind_set, &
2858  qs_kind_set, cell, dft_control, particle_set, pw_env_ext, &
2859  basis_type="RI_AUX", external_vector=wf_vector)
2860 
2861  CALL calc_potential_gpw(psi_l, rho_g, poisson_env, pot_g, mp2_env%potential_parameter, &
2862  no_transfer=.true.)
2863  CALL virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, &
2864  mp2_env%potential_parameter, para_env_ext)
2865  ELSE
2866  CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, mp2_env%potential_parameter)
2867  END IF
2868 
2869  NULLIFY (rs_v)
2870  CALL pw_env_get(pw_env_ext, rs_grids=rs_v)
2871  CALL potential_pw2rs(rs_v, rho_r, pw_env_ext)
2872 
2873  DO iatom = 1, natom
2874 
2875  !only compute if i,j atom pair on correct proc
2876  IF (one_proc_group) THEN
2877  IF (.NOT. iproc_map(iatom) == 1) cycle
2878  END IF
2879 
2880  force_a(:) = 0.0_dp
2881  force_b(:) = 0.0_dp
2882  IF (use_virial) THEN
2883  my_virial_a = 0.0_dp
2884  my_virial_b = 0.0_dp
2885  END IF
2886 
2887  ikind = kind_of(iatom)
2888  atom_a = atom_of_kind(iatom)
2889 
2890  basis_set_a => basis_set_ri_aux(ikind)%gto_basis_set
2891  first_sgfa => basis_set_a%first_sgf
2892  la_max => basis_set_a%lmax
2893  la_min => basis_set_a%lmin
2894  nseta = basis_set_a%nset
2895  nsgfa => basis_set_a%nsgf_set
2896  sphi_a => basis_set_a%sphi
2897  zeta => basis_set_a%zet
2898  npgfa => basis_set_a%npgf
2899 
2900  ra(:) = pbc(particle_set(iatom)%r, cell)
2901 
2902  CALL dbcsr_get_block_p(tmp_g_pq, iatom, jatom, pblock, found)
2903  IF (.NOT. found) cycle
2904 
2905  offset = 0
2906  DO iset = 1, nseta
2907  ncoa = npgfa(iset)*ncoset(la_max(iset))
2908  sgfa = first_sgfa(1, iset)
2909 
2910  ALLOCATE (h_tmp(ncoa, 1)); h_tmp = 0.0_dp
2911  ALLOCATE (i_ab(nsgfa(iset), 1)); i_ab = 0.0_dp
2912  ALLOCATE (pab(ncoa, 1)); pab = 0.0_dp
2913 
2914  i_ab(1:nsgfa(iset), 1) = 2.0_dp*pblock(offset + 1:offset + nsgfa(iset), j_ri - lb_ri)
2915  CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
2916  i_ab(1, 1), nsgfa(iset), 0.0_dp, pab(1, 1), ncoa)
2917 
2918  igrid_level = gaussian_gridlevel(pw_env_ext%gridlevel_info, minval(zeta(:, iset)))
2919 
2920  ! The last three parameters are used to check whether a given function is within the own range.
2921  ! Here, it is always the case, so let's enforce it because mod(0, 1)==0
2922  IF (map_gaussian_here(rs_v(igrid_level), cell%h_inv, ra, 0, 1, 0)) THEN
2923  DO ipgf = 1, npgfa(iset)
2924  o1 = (ipgf - 1)*ncoset(la_max(iset))
2925  igrid_level = gaussian_gridlevel(pw_env_ext%gridlevel_info, zeta(ipgf, iset))
2926 
2927  radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
2928  lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
2929  zetp=zeta(ipgf, iset), &
2930  eps=dft_control%qs_control%eps_gvg_rspace, &
2931  prefactor=1.0_dp, cutoff=1.0_dp)
2932 
2933  CALL integrate_pgf_product( &
2934  la_max=la_max(iset), zeta=zeta(ipgf, iset), la_min=la_min(iset), &
2935  lb_max=0, zetb=0.0_dp, lb_min=0, &
2936  ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
2937  rsgrid=rs_v(igrid_level), &
2938  hab=h_tmp, pab=pab, &
2939  o1=o1, &
2940  o2=0, &
2941  radius=radius, &
2942  calculate_forces=.true., &
2943  force_a=force_a, force_b=force_b, &
2944  use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
2945 
2946  END DO
2947 
2948  END IF
2949 
2950  offset = offset + nsgfa(iset)
2951  DEALLOCATE (pab, h_tmp, i_ab)
2952  END DO !iset
2953 
2954  force(ikind)%mp2_non_sep(:, atom_a) = force(ikind)%mp2_non_sep(:, atom_a) + force_a + force_b
2955  IF (use_virial) h_stress = h_stress + my_virial_a + my_virial_b
2956 
2957  END DO !iatom
2958  END DO !j_RI
2959  END DO !jatom
2960 
2961  IF (use_virial) THEN
2962  CALL auxbas_pw_pool%give_back_pw(rho_g_copy)
2963  DO i_xyz = 1, 3
2964  CALL auxbas_pw_pool%give_back_pw(dvg(i_xyz))
2965  END DO
2966  END IF
2967 
2968  CALL cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_ext, pw_env_ext, &
2969  task_list_ext, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_l)
2970 
2971  CALL dbcsr_release(tmp_g_pq)
2972  CALL dbcsr_distribution_release(dbcsr_dist)
2973  DEALLOCATE (col_dist, row_dist, pgrid)
2974 
2975  CALL mp_para_env_release(para_env_ext)
2976 
2977  CALL timestop(handle)
2978 
2979  END SUBROUTINE get_2c_gpw_forces
2980 
2981 ! **************************************************************************************************
2982 !> \brief Calculate the forces due to the (P|Q) MME integral derivatives
2983 !> \param G_PQ ...
2984 !> \param force ...
2985 !> \param mp2_env ...
2986 !> \param qs_env ...
2987 ! **************************************************************************************************
2988  SUBROUTINE get_2c_mme_forces(G_PQ, force, mp2_env, qs_env)
2989 
2990  TYPE(dbcsr_type), INTENT(INOUT) :: g_pq
2991  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2992  TYPE(mp2_type), INTENT(INOUT) :: mp2_env
2993  TYPE(qs_environment_type), POINTER :: qs_env
2994 
2995  CHARACTER(len=*), PARAMETER :: routinen = 'get_2c_mme_forces'
2996 
2997  INTEGER :: atom_a, atom_b, blk, g_count, handle, i_xyz, iatom, ikind, iset, jatom, jkind, &
2998  jset, natom, nkind, nseta, nsetb, offset_hab_a, offset_hab_b, r_count, sgfa, sgfb
2999  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
3000  INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
3001  npgfb, nsgfa, nsgfb
3002  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
3003  LOGICAL :: found
3004  REAL(dp) :: new_force, pref
3005  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: hab
3006  REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: hdab
3007  REAL(dp), DIMENSION(:, :), POINTER :: pblock
3008  REAL(kind=dp), DIMENSION(3) :: ra, rb
3009  REAL(kind=dp), DIMENSION(:, :), POINTER :: sphi_a, sphi_b, zeta, zetb
3010  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3011  TYPE(cell_type), POINTER :: cell
3012  TYPE(dbcsr_iterator_type) :: iter
3013  TYPE(gto_basis_set_p_type), ALLOCATABLE, &
3014  DIMENSION(:), TARGET :: basis_set_ri_aux
3015  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
3016  TYPE(mp_para_env_type), POINTER :: para_env
3017  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3018  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3019 
3020  NULLIFY (qs_kind_set, basis_set_a, basis_set_b, pblock, particle_set, &
3021  cell, la_max, la_min, lb_min, npgfa, lb_max, npgfb, nsgfa, &
3022  nsgfb, first_sgfa, first_sgfb, sphi_a, sphi_b, zeta, zetb, &
3023  atomic_kind_set, para_env)
3024 
3025  CALL timeset(routinen, handle)
3026 
3027  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind, particle_set=particle_set, &
3028  cell=cell, atomic_kind_set=atomic_kind_set, natom=natom, para_env=para_env)
3029 
3030  CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
3031 
3032  ALLOCATE (basis_set_ri_aux(nkind))
3033  CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
3034 
3035  g_count = 0; r_count = 0
3036 
3037  CALL dbcsr_iterator_start(iter, g_pq)
3038  DO WHILE (dbcsr_iterator_blocks_left(iter))
3039 
3040  CALL dbcsr_iterator_next_block(iter, row=iatom, column=jatom, blk=blk)
3041  CALL dbcsr_get_block_p(g_pq, iatom, jatom, pblock, found)
3042  IF (.NOT. found) cycle
3043  IF (iatom > jatom) cycle
3044  pref = 2.0_dp
3045  IF (iatom == jatom) pref = 1.0_dp
3046 
3047  ikind = kind_of(iatom)
3048  jkind = kind_of(jatom)
3049 
3050  atom_a = atom_of_kind(iatom)
3051  atom_b = atom_of_kind(jatom)
3052 
3053  basis_set_a => basis_set_ri_aux(ikind)%gto_basis_set
3054  first_sgfa => basis_set_a%first_sgf
3055  la_max => basis_set_a%lmax
3056  la_min => basis_set_a%lmin
3057  nseta = basis_set_a%nset
3058  nsgfa => basis_set_a%nsgf_set
3059  sphi_a => basis_set_a%sphi
3060  zeta => basis_set_a%zet
3061  npgfa => basis_set_a%npgf
3062 
3063  basis_set_b => basis_set_ri_aux(jkind)%gto_basis_set
3064  first_sgfb => basis_set_b%first_sgf
3065  lb_max => basis_set_b%lmax
3066  lb_min => basis_set_b%lmin
3067  nsetb = basis_set_b%nset
3068  nsgfb => basis_set_b%nsgf_set
3069  sphi_b => basis_set_b%sphi
3070  zetb => basis_set_b%zet
3071  npgfb => basis_set_b%npgf
3072 
3073  ra(:) = pbc(particle_set(iatom)%r, cell)
3074  rb(:) = pbc(particle_set(jatom)%r, cell)
3075 
3076  ALLOCATE (hab(basis_set_a%nsgf, basis_set_b%nsgf))
3077  ALLOCATE (hdab(3, basis_set_a%nsgf, basis_set_b%nsgf))
3078  hab(:, :) = 0.0_dp
3079  hdab(:, :, :) = 0.0_dp
3080 
3081  offset_hab_a = 0
3082  DO iset = 1, nseta
3083  sgfa = first_sgfa(1, iset)
3084 
3085  offset_hab_b = 0
3086  DO jset = 1, nsetb
3087  sgfb = first_sgfb(1, jset)
3088 
3089  CALL integrate_set_2c(mp2_env%eri_mme_param%par, mp2_env%potential_parameter, la_min(iset), &
3090  la_max(iset), lb_min(jset), lb_max(jset), npgfa(iset), npgfb(jset), &
3091  zeta(:, iset), zetb(:, jset), ra, rb, hab, nsgfa(iset), nsgfb(jset), &
3092  offset_hab_a, offset_hab_b, 0, 0, sphi_a, sphi_b, sgfa, sgfb, &
3093  nsgfa(iset), nsgfb(jset), do_eri_mme, hdab=hdab, &
3094  g_count=g_count, r_count=r_count)
3095 
3096  offset_hab_b = offset_hab_b + nsgfb(jset)
3097  END DO
3098  offset_hab_a = offset_hab_a + nsgfa(iset)
3099  END DO
3100 
3101  DO i_xyz = 1, 3
3102  new_force = pref*sum(pblock(:, :)*hdab(i_xyz, :, :))
3103  force(ikind)%mp2_non_sep(i_xyz, atom_a) = force(ikind)%mp2_non_sep(i_xyz, atom_a) + new_force
3104  force(jkind)%mp2_non_sep(i_xyz, atom_b) = force(jkind)%mp2_non_sep(i_xyz, atom_b) - new_force
3105  END DO
3106 
3107  DEALLOCATE (hab, hdab)
3108  END DO
3109  CALL dbcsr_iterator_stop(iter)
3110 
3111  CALL cp_eri_mme_update_local_counts(mp2_env%eri_mme_param, para_env, g_count_2c=g_count, r_count_2c=r_count)
3112 
3113  CALL timestop(handle)
3114 
3115  END SUBROUTINE get_2c_mme_forces
3116 
3117 ! **************************************************************************************************
3118 !> \brief This routines gather all the force updates due to the response density and the trace with F
3119 !> Also update the forces due to the SCF density for XC and exact exchange
3120 !> \param p_env the p_env coming from the response calculation
3121 !> \param matrix_hz the matrix going into the RHS of the response equation
3122 !> \param matrix_p_F the density matrix with which we evaluate Trace[P*F]
3123 !> \param matrix_p_F_admm ...
3124 !> \param qs_env ...
3125 !> \note very much inspired from the response_force routine in response_solver.F, especially for virial
3126 ! **************************************************************************************************
3127  SUBROUTINE update_im_time_forces(p_env, matrix_hz, matrix_p_F, matrix_p_F_admm, qs_env)
3128 
3129  TYPE(qs_p_env_type), POINTER :: p_env
3130  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz, matrix_p_f, matrix_p_f_admm
3131  TYPE(qs_environment_type), POINTER :: qs_env
3132 
3133  CHARACTER(len=*), PARAMETER :: routinen = 'update_im_time_forces'
3134 
3135  INTEGER :: handle, i, idens, ispin, n_rep_hf, nao, &
3136  nao_aux, nder, nimages, nocc, nspins
3137  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3138  LOGICAL :: do_exx, do_hfx, do_tau, do_tau_admm, &
3139  use_virial
3140  REAL(dp) :: dummy_real1, dummy_real2, ehartree, &
3141  eps_ppnl, exc, focc
3142  REAL(dp), DIMENSION(3, 3) :: h_stress, pv_loc
3143  TYPE(admm_type), POINTER :: admm_env
3144  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3145  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: current_density, current_density_admm, &
3146  current_mat_h, matrix_p_mp2, matrix_p_mp2_admm, matrix_s, matrix_s_aux_fit, matrix_w, &
3147  rho_ao, rho_ao_aux, scrm, scrm_admm
3148  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: dbcsr_work_h, dbcsr_work_p
3149  TYPE(dbcsr_type) :: dbcsr_work
3150  TYPE(dft_control_type), POINTER :: dft_control
3151  TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
3152  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3153  TYPE(mp_para_env_type), POINTER :: para_env
3154  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3155  POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
3156  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3157  TYPE(pw_c1d_gs_type) :: rho_tot_gspace, rhoz_tot_gspace, &
3158  zv_hartree_gspace
3159  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rhoz_g
3160  TYPE(pw_env_type), POINTER :: pw_env
3161  TYPE(pw_poisson_type), POINTER :: poisson_env
3162  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
3163  TYPE(pw_r3d_rs_type) :: vh_rspace, vhxc_rspace, zv_hartree_rspace
3164  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rhoz_r, tauz_r, v_xc, v_xc_tau, &
3165  vadmm_rspace, vtau_rspace, vxc_rspace
3166  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3167  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3168  TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
3169  TYPE(section_vals_type), POINTER :: hfx_section, xc_section
3170  TYPE(task_list_type), POINTER :: task_list_aux_fit
3171  TYPE(virial_type), POINTER :: virial
3172 
3173  NULLIFY (scrm, rho, dft_control, matrix_p_mp2, matrix_s, matrix_p_mp2_admm, admm_env, sab_orb, &
3174  cell_to_index, dbcsr_work_p, dbcsr_work_h, sac_ae, sac_ppl, sap_ppnl, force, virial, &
3175  qs_kind_set, atomic_kind_set, particle_set, pw_env, poisson_env, auxbas_pw_pool, &
3176  task_list_aux_fit, matrix_s_aux_fit, scrm_admm, rho_aux_fit, rho_ao_aux, x_data, &
3177  hfx_section, xc_section, para_env, rhoz_g, rhoz_r, tauz_r, v_xc, v_xc_tau, &
3178  vxc_rspace, vtau_rspace, vadmm_rspace, rho_ao, matrix_w)
3179 
3180  CALL timeset(routinen, handle)
3181 
3182  CALL get_qs_env(qs_env, rho=rho, dft_control=dft_control, matrix_s=matrix_s, admm_env=admm_env, &
3183  sab_orb=sab_orb, sac_ae=sac_ae, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl, force=force, &
3184  virial=virial, particle_set=particle_set, qs_kind_set=qs_kind_set, &
3185  atomic_kind_set=atomic_kind_set, x_data=x_data, para_env=para_env)
3186  nspins = dft_control%nspins
3187 
3188  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
3189  IF (use_virial) virial%pv_calculate = .true.
3190 
3191  !Whether we replace the force/energy of SCF XC with HF in RPA
3192  do_exx = .false.
3193  IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
3194  hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
3195  CALL section_vals_get(hfx_section, explicit=do_exx)
3196  END IF
3197 
3198  !Get the mp2 density matrix which is p_env%p1 + matrix_p_F
3199  CALL get_qs_env(qs_env, matrix_p_mp2=matrix_p_mp2, matrix_p_mp2_admm=matrix_p_mp2_admm)
3200 
3201  !The kinetic term (only response density)
3202  NULLIFY (scrm)
3203  IF (nspins == 2) CALL dbcsr_add(matrix_p_mp2(1)%matrix, matrix_p_mp2(2)%matrix, 1.0_dp, 1.0_dp)
3204  CALL build_kinetic_matrix(qs_env%ks_env, matrix_t=scrm, &
3205  matrix_name="KINETIC ENERGY MATRIX", &
3206  basis_type="ORB", &
3207  sab_nl=sab_orb, calculate_forces=.true., &
3208  matrix_p=matrix_p_mp2(1)%matrix)
3209  IF (nspins == 2) CALL dbcsr_add(matrix_p_mp2(1)%matrix, matrix_p_mp2(2)%matrix, 1.0_dp, -1.0_dp)
3210  CALL dbcsr_deallocate_matrix_set(scrm)
3211 
3212  !The pseudo-potential terms (only reponse density)
3213  CALL dbcsr_allocate_matrix_set(scrm, nspins)
3214  DO ispin = 1, nspins
3215  ALLOCATE (scrm(ispin)%matrix)
3216  CALL dbcsr_create(scrm(ispin)%matrix, template=matrix_s(1)%matrix)
3217  CALL dbcsr_copy(scrm(ispin)%matrix, matrix_s(1)%matrix)
3218  CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
3219  END DO
3220 
3221  nder = 1
3222  nimages = 1
3223  NULLIFY (cell_to_index)
3224  ALLOCATE (dbcsr_work_p(nspins, 1), dbcsr_work_h(nspins, 1))
3225  DO ispin = 1, nspins
3226  dbcsr_work_p(ispin, 1)%matrix => matrix_p_mp2(ispin)%matrix
3227  dbcsr_work_h(ispin, 1)%matrix => scrm(ispin)%matrix
3228  END DO
3229 
3230  IF (ASSOCIATED(sac_ae)) THEN
3231  CALL build_core_ae(dbcsr_work_h, dbcsr_work_p, force, &
3232  virial, .true., use_virial, nder, &
3233  qs_kind_set, atomic_kind_set, particle_set, &
3234  sab_orb, sac_ae, nimages, cell_to_index)
3235  END IF
3236 
3237  IF (ASSOCIATED(sac_ppl)) THEN
3238  CALL build_core_ppl(dbcsr_work_h, dbcsr_work_p, force, &
3239  virial, .true., use_virial, nder, &
3240  qs_kind_set, atomic_kind_set, particle_set, &
3241  sab_orb, sac_ppl, nimages, cell_to_index, "ORB")
3242  END IF
3243 
3244  IF (ASSOCIATED(sap_ppnl)) THEN
3245  eps_ppnl = dft_control%qs_control%eps_ppnl
3246  CALL build_core_ppnl(dbcsr_work_h, dbcsr_work_p, force, &
3247  virial, .true., use_virial, nder, &
3248  qs_kind_set, atomic_kind_set, particle_set, &
3249  sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, "ORB")
3250  END IF
3251  DEALLOCATE (dbcsr_work_p, dbcsr_work_h)
3252 
3253  IF (use_virial) THEN
3254  h_stress = 0.0_dp
3255  virial%pv_xc = 0.0_dp
3256  NULLIFY (vxc_rspace, vtau_rspace, vadmm_rspace)
3257  CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, &
3258  dummy_real1, dummy_real2, h_stress)
3259  virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
3260  virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
3261  IF (.NOT. do_exx) THEN
3262  !if RPA EXX, then do not consider XC virial (replaced by RPA%HF virial)
3263  virial%pv_exc = virial%pv_exc - virial%pv_xc
3264  virial%pv_virial = virial%pv_virial - virial%pv_xc
3265  END IF
3266  ELSE
3267  CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, dummy_real1, dummy_real2)
3268  END IF
3269  do_tau = ASSOCIATED(vtau_rspace)
3270 
3271  !Core forces from the SCF
3272  CALL integrate_v_core_rspace(vh_rspace, qs_env)
3273 
3274  !The Hartree-xc potential term, P*dVHxc (mp2 + SCF density x deriv of the SCF potential)
3275  !Get the total density
3276  CALL qs_rho_get(rho, rho_ao=rho_ao)
3277  DO ispin = 1, nspins
3278  CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
3279  END DO
3280 
3281  CALL get_qs_env(qs_env, pw_env=pw_env)
3282  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
3283  poisson_env=poisson_env)
3284  CALL auxbas_pw_pool%create_pw(vhxc_rspace)
3285 
3286  IF (use_virial) pv_loc = virial%pv_virial
3287 
3288  IF (do_exx) THEN
3289  !Only want response XC contribution, but SCF+response Hartree contribution
3290  DO ispin = 1, nspins
3291  !Hartree
3292  CALL pw_transfer(vh_rspace, vhxc_rspace)
3293  CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
3294  hmat=scrm(ispin), pmat=rho_ao(ispin), &
3295  qs_env=qs_env, calculate_forces=.true.)
3296  !XC
3297  CALL pw_transfer(vxc_rspace(ispin), vhxc_rspace)
3298  CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
3299  hmat=scrm(ispin), pmat=matrix_p_mp2(ispin), &
3300  qs_env=qs_env, calculate_forces=.true.)
3301  IF (do_tau) THEN
3302  CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
3303  hmat=scrm(ispin), pmat=matrix_p_mp2(ispin), &
3304  qs_env=qs_env, calculate_forces=.true., compute_tau=.true.)
3305  END IF
3306  END DO
3307  ELSE
3308  DO ispin = 1, nspins
3309  CALL pw_transfer(vh_rspace, vhxc_rspace)
3310  CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
3311  CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
3312  hmat=scrm(ispin), pmat=rho_ao(ispin), &
3313  qs_env=qs_env, calculate_forces=.true.)
3314  IF (do_tau) THEN
3315  CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
3316  hmat=scrm(ispin), pmat=rho_ao(ispin), &
3317  qs_env=qs_env, calculate_forces=.true., compute_tau=.true.)
3318  END IF
3319  END DO
3320  END IF
3321  CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
3322 
3323  IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3324 
3325  !The admm projection contribution (mp2 + SCF densities). If EXX, then only mp2 density
3326  IF (dft_control%do_admm) THEN
3327  CALL get_admm_env(admm_env, task_list_aux_fit=task_list_aux_fit, rho_aux_fit=rho_aux_fit, &
3328  matrix_s_aux_fit=matrix_s_aux_fit)
3329  CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux)
3330  CALL dbcsr_allocate_matrix_set(scrm_admm, nspins)
3331  DO ispin = 1, nspins
3332  ALLOCATE (scrm_admm(ispin)%matrix)
3333  CALL dbcsr_create(scrm_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
3334  CALL dbcsr_copy(scrm_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
3335  CALL dbcsr_set(scrm_admm(ispin)%matrix, 0.0_dp)
3336  END DO
3337 
3338  IF (use_virial) pv_loc = virial%pv_virial
3339  IF (.NOT. qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
3340  DO ispin = 1, nspins
3341  IF (do_exx) THEN
3342  CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), &
3343  hmat=scrm_admm(ispin), pmat=matrix_p_mp2_admm(ispin), &
3344  qs_env=qs_env, calculate_forces=.true., &
3345  basis_type="AUX_FIT", task_list_external=task_list_aux_fit)
3346  ELSE
3347  CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
3348  CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), &
3349  hmat=scrm_admm(ispin), pmat=rho_ao_aux(ispin), &
3350  qs_env=qs_env, calculate_forces=.true., &
3351  basis_type="AUX_FIT", task_list_external=task_list_aux_fit)
3352  CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
3353  END IF
3354  END DO
3355  END IF
3356  IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3357 
3358  CALL tddft_hfx_matrix(scrm_admm, rho_ao_aux, qs_env, .false., .false.)
3359 
3360  IF (do_exx) THEN
3361  CALL admm_projection_derivative(qs_env, scrm_admm, matrix_p_mp2)
3362  ELSE
3363  CALL admm_projection_derivative(qs_env, scrm_admm, rho_ao)
3364  END IF
3365  END IF
3366 
3367  !The exact-exchange term (mp2 + SCF densities)
3368  xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
3369  hfx_section => section_vals_get_subs_vals(xc_section, "HF")
3370  CALL section_vals_get(hfx_section, explicit=do_hfx)
3371 
3372  IF (do_hfx) THEN
3373  CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
3374  cpassert(n_rep_hf == 1)
3375  IF (use_virial) virial%pv_fock_4c = 0.0_dp
3376 
3377  !In case of EXX, only want to response HFX forces, as the SCF will change according to RI_RPA%HF
3378  IF (do_exx) THEN
3379  IF (dft_control%do_admm) THEN
3380  CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
3381  CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux, rho_ao_kp=dbcsr_work_p)
3382  IF (x_data(1, 1)%do_hfx_ri) THEN
3383 
3384  CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3385  x_data(1, 1)%general_parameter%fraction, &
3386  rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2_admm, &
3387  use_virial=use_virial, resp_only=.true.)
3388  ELSE
3389  CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2_admm, hfx_section, para_env, &
3390  1, use_virial, resp_only=.true.)
3391  END IF
3392  ELSE
3393  DO ispin = 1, nspins
3394  CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
3395  END DO
3396  CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3397  IF (x_data(1, 1)%do_hfx_ri) THEN
3398 
3399  CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3400  x_data(1, 1)%general_parameter%fraction, &
3401  rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2, &
3402  use_virial=use_virial, resp_only=.true.)
3403  ELSE
3404  CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2, hfx_section, para_env, &
3405  1, use_virial, resp_only=.true.)
3406  END IF
3407  DO ispin = 1, nspins
3408  CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
3409  END DO
3410  END IF !admm
3411 
3412  ELSE !No Exx
3413  IF (dft_control%do_admm) THEN
3414  CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
3415  CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux, rho_ao_kp=dbcsr_work_p)
3416  DO ispin = 1, nspins
3417  CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
3418  END DO
3419  IF (x_data(1, 1)%do_hfx_ri) THEN
3420 
3421  CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3422  x_data(1, 1)%general_parameter%fraction, &
3423  rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2_admm, &
3424  use_virial=use_virial, resp_only=.false.)
3425  ELSE
3426  CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2_admm, hfx_section, para_env, &
3427  1, use_virial, resp_only=.false.)
3428  END IF
3429  DO ispin = 1, nspins
3430  CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
3431  END DO
3432  ELSE
3433  CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3434  IF (x_data(1, 1)%do_hfx_ri) THEN
3435 
3436  CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3437  x_data(1, 1)%general_parameter%fraction, &
3438  rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2, &
3439  use_virial=use_virial, resp_only=.false.)
3440  ELSE
3441  CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2, hfx_section, para_env, &
3442  1, use_virial, resp_only=.false.)
3443  END IF
3444  END IF
3445  END IF !do_exx
3446 
3447  IF (use_virial) THEN
3448  virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
3449  virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
3450  END IF
3451  END IF
3452 
3453  !retrieve the SCF density
3454  CALL qs_rho_get(rho, rho_ao=rho_ao)
3455  DO ispin = 1, nspins
3456  CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
3457  END DO
3458 
3459  !From here, we need to do everything twice. Once for the response density, and once for the
3460  !density that is used for the trace Tr[P*F]. The reason is that the former is needed for the
3461  !eventual overlap contribution from matrix_wz
3462  !Only with the mp2 density
3463 
3464  ALLOCATE (current_density(nspins), current_mat_h(nspins), current_density_admm(nspins))
3465  DO idens = 1, 2
3466  DO ispin = 1, nspins
3467  IF (idens == 1) THEN
3468  current_density(ispin)%matrix => matrix_p_f(ispin)%matrix
3469  current_mat_h(ispin)%matrix => scrm(ispin)%matrix
3470  IF (dft_control%do_admm) current_density_admm(ispin)%matrix => matrix_p_f_admm(ispin)%matrix
3471  ELSE
3472  current_density(ispin)%matrix => p_env%p1(ispin)%matrix
3473  current_mat_h(ispin)%matrix => matrix_hz(ispin)%matrix
3474  IF (dft_control%do_admm) current_density_admm(ispin)%matrix => p_env%p1_admm(ispin)%matrix
3475  END IF
3476  END DO
3477 
3478  !The core-denstiy derivative
3479  ALLOCATE (rhoz_r(nspins), rhoz_g(nspins))
3480  DO ispin = 1, nspins
3481  CALL auxbas_pw_pool%create_pw(rhoz_r(ispin))
3482  CALL auxbas_pw_pool%create_pw(rhoz_g(ispin))
3483  END DO
3484  CALL auxbas_pw_pool%create_pw(rhoz_tot_gspace)
3485  CALL auxbas_pw_pool%create_pw(zv_hartree_rspace)
3486  CALL auxbas_pw_pool%create_pw(zv_hartree_gspace)
3487 
3488  CALL pw_zero(rhoz_tot_gspace)
3489  DO ispin = 1, nspins
3490  CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density(ispin)%matrix, &
3491  rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin))
3492  CALL pw_axpy(rhoz_g(ispin), rhoz_tot_gspace)
3493  END DO
3494 
3495  IF (use_virial) THEN
3496 
3497  CALL get_qs_env(qs_env, rho=rho)
3498  CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
3499 
3500  CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
3501 
3502  h_stress(:, :) = 0.0_dp
3503  CALL pw_poisson_solve(poisson_env, &
3504  density=rhoz_tot_gspace, &
3505  ehartree=ehartree, &
3506  vhartree=zv_hartree_gspace, &
3507  h_stress=h_stress, &
3508  aux_density=rho_tot_gspace)
3509 
3510  CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
3511 
3512  !Green contribution
3513  virial%pv_ehartree = virial%pv_ehartree + 2.0_dp*h_stress/real(para_env%num_pe, dp)
3514  virial%pv_virial = virial%pv_virial + 2.0_dp*h_stress/real(para_env%num_pe, dp)
3515 
3516  ELSE
3517  CALL pw_poisson_solve(poisson_env, rhoz_tot_gspace, ehartree, &
3518  zv_hartree_gspace)
3519  END IF
3520 
3521  CALL pw_transfer(zv_hartree_gspace, zv_hartree_rspace)
3522  CALL pw_scale(zv_hartree_rspace, zv_hartree_rspace%pw_grid%dvol)
3523  CALL integrate_v_core_rspace(zv_hartree_rspace, qs_env)
3524 
3525  IF (do_tau) THEN
3526  block
3527  TYPE(pw_c1d_gs_type) :: tauz_g
3528  CALL auxbas_pw_pool%create_pw(tauz_g)
3529  ALLOCATE (tauz_r(nspins))
3530  DO ispin = 1, nspins
3531  CALL auxbas_pw_pool%create_pw(tauz_r(ispin))
3532 
3533  CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density(ispin)%matrix, &
3534  rho=tauz_r(ispin), rho_gspace=tauz_g, compute_tau=.true.)
3535  END DO
3536  CALL auxbas_pw_pool%give_back_pw(tauz_g)
3537  END block
3538  END IF
3539 
3540  !Volume contribution to the virial
3541  IF (use_virial) THEN
3542  !Volume contribution
3543  exc = 0.0_dp
3544  DO ispin = 1, nspins
3545  exc = exc + pw_integral_ab(rhoz_r(ispin), vxc_rspace(ispin))/ &
3546  vxc_rspace(ispin)%pw_grid%dvol
3547  END DO
3548  IF (ASSOCIATED(vtau_rspace)) THEN
3549  DO ispin = 1, nspins
3550  exc = exc + pw_integral_ab(tauz_r(ispin), vtau_rspace(ispin))/ &
3551  vtau_rspace(ispin)%pw_grid%dvol
3552  END DO
3553  END IF
3554  DO i = 1, 3
3555  virial%pv_ehartree(i, i) = virial%pv_ehartree(i, i) - 4.0_dp*ehartree/real(para_env%num_pe, dp)
3556  virial%pv_exc(i, i) = virial%pv_exc(i, i) - exc/real(para_env%num_pe, dp)
3557  virial%pv_virial(i, i) = virial%pv_virial(i, i) - 4.0_dp*ehartree/real(para_env%num_pe, dp) &
3558  - exc/real(para_env%num_pe, dp)
3559  END DO
3560  END IF
3561 
3562  !The xc-kernel term.
3563  IF (dft_control%do_admm) THEN
3564  CALL get_qs_env(qs_env, admm_env=admm_env)
3565  xc_section => admm_env%xc_section_primary
3566  ELSE
3567  xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
3568  END IF
3569 
3570  IF (use_virial) virial%pv_xc = 0.0_dp
3571 
3572  CALL create_kernel(qs_env, &
3573  vxc=v_xc, &
3574  vxc_tau=v_xc_tau, &
3575  rho=rho, &
3576  rho1_r=rhoz_r, &
3577  rho1_g=rhoz_g, &
3578  tau1_r=tauz_r, &
3579  xc_section=xc_section, &
3580  compute_virial=use_virial, &
3581  virial_xc=virial%pv_xc)
3582 
3583  IF (use_virial) THEN
3584  virial%pv_exc = virial%pv_exc + virial%pv_xc
3585  virial%pv_virial = virial%pv_virial + virial%pv_xc
3586 
3587  pv_loc = virial%pv_virial
3588  END IF
3589 
3590  CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3591  DO ispin = 1, nspins
3592  CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
3593  CALL pw_axpy(zv_hartree_rspace, v_xc(ispin))
3594  CALL integrate_v_rspace(qs_env=qs_env, &
3595  v_rspace=v_xc(ispin), &
3596  hmat=current_mat_h(ispin), &
3597  pmat=dbcsr_work_p(ispin, 1), &
3598  calculate_forces=.true.)
3599  CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
3600  END DO
3601  CALL auxbas_pw_pool%give_back_pw(rhoz_tot_gspace)
3602  CALL auxbas_pw_pool%give_back_pw(zv_hartree_rspace)
3603  CALL auxbas_pw_pool%give_back_pw(zv_hartree_gspace)
3604  DEALLOCATE (v_xc)
3605 
3606  IF (do_tau) THEN
3607  DO ispin = 1, nspins
3608  CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
3609  CALL integrate_v_rspace(qs_env=qs_env, &
3610  v_rspace=v_xc_tau(ispin), &
3611  hmat=current_mat_h(ispin), &
3612  pmat=dbcsr_work_p(ispin, 1), &
3613  compute_tau=.true., &
3614  calculate_forces=.true.)
3615  CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
3616  END DO
3617  DEALLOCATE (v_xc_tau)
3618  END IF
3619 
3620  IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3621 
3622  IF (do_hfx) THEN
3623  IF (dft_control%do_admm) THEN
3624  DO ispin = 1, nspins
3625  CALL dbcsr_set(scrm_admm(ispin)%matrix, 0.0_dp)
3626  END DO
3627  CALL qs_rho_get(rho_aux_fit, tau_r_valid=do_tau_admm)
3628 
3629  IF (.NOT. admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
3630  CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
3631  DO ispin = 1, nspins
3632  CALL pw_zero(rhoz_r(ispin))
3633  CALL pw_zero(rhoz_g(ispin))
3634  CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density_admm(ispin)%matrix, &
3635  rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin), &
3636  basis_type="AUX_FIT", task_list_external=task_list_aux_fit)
3637  END DO
3638 
3639  IF (do_tau_admm) THEN
3640  block
3641  TYPE(pw_c1d_gs_type) :: tauz_g
3642  CALL auxbas_pw_pool%create_pw(tauz_g)
3643  DO ispin = 1, nspins
3644  CALL pw_zero(tauz_r(ispin))
3645  CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density(ispin)%matrix, &
3646  rho=tauz_r(ispin), rho_gspace=tauz_g, &
3647  basis_type="AUX_FIT", task_list_external=task_list_aux_fit, &
3648  compute_tau=.true.)
3649  END DO
3650  CALL auxbas_pw_pool%give_back_pw(tauz_g)
3651  END block
3652  END IF
3653 
3654  !Volume contribution to the virial
3655  IF (use_virial) THEN
3656  exc = 0.0_dp
3657  DO ispin = 1, nspins
3658  exc = exc + pw_integral_ab(rhoz_r(ispin), vadmm_rspace(ispin))/ &
3659  vadmm_rspace(ispin)%pw_grid%dvol
3660  END DO
3661  DO i = 1, 3
3662  virial%pv_exc(i, i) = virial%pv_exc(i, i) - exc/real(para_env%num_pe, dp)
3663  virial%pv_virial(i, i) = virial%pv_virial(i, i) - exc/real(para_env%num_pe, dp)
3664  END DO
3665 
3666  virial%pv_xc = 0.0_dp
3667  END IF
3668 
3669  xc_section => admm_env%xc_section_aux
3670  CALL create_kernel(qs_env, v_xc, v_xc_tau, rho_aux_fit, rhoz_r, rhoz_g, tauz_r, xc_section, &
3671  compute_virial=use_virial, virial_xc=virial%pv_xc)
3672 
3673  IF (use_virial) THEN
3674  virial%pv_exc = virial%pv_exc + virial%pv_xc
3675  virial%pv_virial = virial%pv_virial + virial%pv_xc
3676 
3677  pv_loc = virial%pv_virial
3678  END IF
3679 
3680  CALL qs_rho_get(rho_aux_fit, rho_ao_kp=dbcsr_work_p)
3681  DO ispin = 1, nspins
3682  CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
3683  CALL integrate_v_rspace(qs_env=qs_env, &
3684  v_rspace=v_xc(ispin), &
3685  hmat=scrm_admm(ispin), &
3686  pmat=dbcsr_work_p(ispin, 1), &
3687  calculate_forces=.true., &
3688  basis_type="AUX_FIT", &
3689  task_list_external=task_list_aux_fit)
3690  CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
3691  END DO
3692  DEALLOCATE (v_xc)
3693 
3694  IF (do_tau_admm) THEN
3695  DO ispin = 1, nspins
3696  CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
3697  CALL integrate_v_rspace(qs_env=qs_env, &
3698  v_rspace=v_xc_tau(ispin), &
3699  hmat=scrm_admm(ispin), &
3700  pmat=dbcsr_work_p(ispin, 1), &
3701  calculate_forces=.true., &
3702  basis_type="AUX_FIT", &
3703  task_list_external=task_list_aux_fit, &
3704  compute_tau=.true.)
3705  CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
3706  END DO
3707  DEALLOCATE (v_xc_tau)
3708  END IF
3709 
3710  IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3711  END IF
3712 
3713  CALL tddft_hfx_matrix(scrm_admm, current_density_admm, qs_env, .false., .false.)
3714 
3715  CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3716  CALL admm_projection_derivative(qs_env, scrm_admm, dbcsr_work_p(:, 1))
3717 
3718  !If response density, need to get matrix_hz contribution
3719  CALL dbcsr_create(dbcsr_work, template=matrix_s(1)%matrix)
3720  IF (idens == 2) THEN
3721  nao = admm_env%nao_orb
3722  nao_aux = admm_env%nao_aux_fit
3723  DO ispin = 1, nspins
3724  CALL dbcsr_copy(dbcsr_work, matrix_hz(ispin)%matrix)
3725  CALL dbcsr_set(dbcsr_work, 0.0_dp)
3726 
3727  CALL cp_dbcsr_sm_fm_multiply(scrm_admm(ispin)%matrix, admm_env%A, &
3728  admm_env%work_aux_orb, nao)
3729  CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
3730  1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
3731  admm_env%work_orb_orb)
3732  CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_work, keep_sparsity=.true.)
3733  CALL dbcsr_add(matrix_hz(ispin)%matrix, dbcsr_work, 1.0_dp, 1.0_dp)
3734  END DO
3735  END IF
3736 
3737  CALL dbcsr_release(dbcsr_work)
3738  ELSE !no admm
3739 
3740  !Need the contribution to matrix_hz as well
3741  IF (idens == 2) THEN
3742  CALL tddft_hfx_matrix(matrix_hz, current_density, qs_env, .false., .false.)
3743  END IF
3744  END IF !admm
3745  END IF !do_hfx
3746 
3747  DO ispin = 1, nspins
3748  CALL auxbas_pw_pool%give_back_pw(rhoz_r(ispin))
3749  CALL auxbas_pw_pool%give_back_pw(rhoz_g(ispin))
3750  END DO
3751  DEALLOCATE (rhoz_r, rhoz_g)
3752 
3753  IF (do_tau) THEN
3754  DO ispin = 1, nspins
3755  CALL auxbas_pw_pool%give_back_pw(tauz_r(ispin))
3756  END DO
3757  DEALLOCATE (tauz_r)
3758  END IF
3759  END DO !idens
3760  CALL dbcsr_deallocate_matrix_set(scrm_admm)
3761 
3762  DEALLOCATE (current_density, current_mat_h, current_density_admm)
3763  CALL dbcsr_deallocate_matrix_set(scrm)
3764 
3765  !The energy weighted and overlap term. ONLY with the response density
3766  focc = 2.0_dp
3767  IF (nspins == 2) focc = 1.0_dp
3768  CALL get_qs_env(qs_env, mos=mos)
3769  DO ispin = 1, nspins
3770  CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
3771  CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
3772  p_env%w1(ispin)%matrix, focc, nocc)
3773  END DO
3774  IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, p_env%w1(2)%matrix, 1.0_dp, 1.0_dp)
3775 
3776  !Add to it the SCF W matrix, except if EXX (because taken care of by HF response)
3777  IF (.NOT. do_exx) THEN
3778  CALL qs_energies_compute_w(qs_env, calc_forces=.true.)
3779  CALL get_qs_env(qs_env, matrix_w=matrix_w)
3780  CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(1)%matrix, 1.0_dp, 1.0_dp)
3781  IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(2)%matrix, 1.0_dp, 1.0_dp)
3782  END IF
3783 
3784  NULLIFY (scrm)
3785  CALL build_overlap_matrix(qs_env%ks_env, matrix_s=scrm, &
3786  matrix_name="OVERLAP MATRIX", &
3787  basis_type_a="ORB", basis_type_b="ORB", &
3788  sab_nl=sab_orb, calculate_forces=.true., &
3789  matrix_p=p_env%w1(1)%matrix)
3790 
3791  IF (.NOT. do_exx) THEN
3792  CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(1)%matrix, 1.0_dp, -1.0_dp)
3793  IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(2)%matrix, 1.0_dp, -1.0_dp)
3794  DO ispin = 1, nspins
3795  CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
3796  END DO
3797  END IF
3798 
3799  IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, p_env%w1(2)%matrix, 1.0_dp, -1.0_dp)
3800  CALL dbcsr_deallocate_matrix_set(scrm)
3801 
3802  IF (use_virial) virial%pv_calculate = .false.
3803 
3804  !clean-up
3805  CALL auxbas_pw_pool%give_back_pw(vh_rspace)
3806 
3807  DO ispin = 1, nspins
3808  CALL auxbas_pw_pool%give_back_pw(vxc_rspace(ispin))
3809  IF (ASSOCIATED(vtau_rspace)) THEN
3810  CALL auxbas_pw_pool%give_back_pw(vtau_rspace(ispin))
3811  END IF
3812  IF (ASSOCIATED(vadmm_rspace)) THEN
3813  CALL auxbas_pw_pool%give_back_pw(vadmm_rspace(ispin))
3814  END IF
3815  END DO
3816  DEALLOCATE (vxc_rspace)
3817  IF (ASSOCIATED(vtau_rspace)) DEALLOCATE (vtau_rspace)
3818  IF (ASSOCIATED(vadmm_rspace)) DEALLOCATE (vadmm_rspace)
3819 
3820  CALL timestop(handle)
3821 
3822  END SUBROUTINE update_im_time_forces
3823 
3824 ! **************************************************************************************************
3825 !> \brief Iteratively builds the matrix Y = sum_k Y_k until convergence, where
3826 !> Y_k = 1/k*2^n (A/2^n) Y_k-1 + 1/k!*2^n * PR(n) * (A/2^n)^(k-1)
3827 !> n is chosen such that the norm of A is < 1 (and e^A converges fast)
3828 !> PR(n) = e^(A/2^n)*PR(n-1) + PR(n-1)*e^(A/2^n), PR(0) = P*R^T
3829 !> \param Y ...
3830 !> \param A ...
3831 !> \param P ...
3832 !> \param R ...
3833 !> \param filter_eps ...
3834 ! **************************************************************************************************
3835  SUBROUTINE build_y_matrix(Y, A, P, R, filter_eps)
3836 
3837  TYPE(dbcsr_type), INTENT(OUT) :: y
3838  TYPE(dbcsr_type), INTENT(INOUT) :: a, p, r
3839  REAL(dp), INTENT(IN) :: filter_eps
3840 
3841  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_Y_matrix'
3842 
3843  INTEGER :: handle, k, n
3844  REAL(dp) :: norm_scalar, threshold
3845  TYPE(dbcsr_type) :: a2n, exp_a2n, prn, work, work2, yk
3846 
3847  CALL timeset(routinen, handle)
3848 
3849  threshold = 1.0e-16_dp
3850 
3851  !Find n such that norm(A) < 1 and we insure convergence of the exponential
3852  norm_scalar = dbcsr_frobenius_norm(a)
3853 
3854  !checked: result invariant with value of n
3855  n = 1
3856  DO
3857  IF ((norm_scalar/2.0_dp**n) < 1.0_dp) EXIT
3858  n = n + 1
3859  END DO
3860 
3861  !Calculate PR(n) recursively
3862  CALL dbcsr_create(prn, template=a, matrix_type=dbcsr_type_no_symmetry)
3863  CALL dbcsr_create(work, template=a, matrix_type=dbcsr_type_no_symmetry)
3864  CALL dbcsr_multiply('N', 'N', 1.0_dp, p, r, 0.0_dp, work, filter_eps=filter_eps)
3865  CALL dbcsr_create(exp_a2n, template=a, matrix_type=dbcsr_type_no_symmetry)
3866 
3867  DO k = 1, n
3868  CALL matrix_exponential(exp_a2n, a, 1.0_dp, 0.5_dp**k, threshold)
3869  CALL dbcsr_multiply('N', 'N', 1.0_dp, exp_a2n, work, 0.0_dp, prn, filter_eps=filter_eps)
3870  CALL dbcsr_multiply('N', 'N', 1.0_dp, work, exp_a2n, 1.0_dp, prn, filter_eps=filter_eps)
3871  CALL dbcsr_copy(work, prn)
3872  END DO
3873  CALL dbcsr_release(exp_a2n)
3874 
3875  !Calculate Y iteratively, until convergence
3876  CALL dbcsr_create(a2n, template=a, matrix_type=dbcsr_type_no_symmetry)
3877  CALL dbcsr_copy(a2n, a)
3878  CALL dbcsr_scale(a2n, 0.5_dp**n)
3879  CALL dbcsr_create(y, template=a, matrix_type=dbcsr_type_no_symmetry)
3880  CALL dbcsr_create(yk, template=a, matrix_type=dbcsr_type_no_symmetry)
3881  CALL dbcsr_create(work2, template=a, matrix_type=dbcsr_type_no_symmetry)
3882 
3883  !k=1
3884  CALL dbcsr_scale(prn, 0.5_dp**n)
3885  CALL dbcsr_copy(work, prn)
3886  CALL dbcsr_copy(work2, prn)
3887  CALL dbcsr_add(y, prn, 1.0_dp, 1.0_dp)
3888 
3889  k = 1
3890  DO
3891  k = k + 1
3892  CALL dbcsr_multiply('N', 'N', 1.0_dp/real(k, dp), a2n, work, 0.0_dp, yk, filter_eps=filter_eps)
3893  CALL dbcsr_multiply('N', 'N', 1.0_dp/real(k, dp), work2, a2n, 0.0_dp, prn, filter_eps=filter_eps)
3894 
3895  CALL dbcsr_add(yk, prn, 1.0_dp, 1.0_dp)
3896  CALL dbcsr_add(y, yk, 1.0_dp, 1.0_dp)
3897 
3898  IF (dbcsr_frobenius_norm(yk) < threshold) EXIT
3899  CALL dbcsr_copy(work, yk)
3900  CALL dbcsr_copy(work2, prn)
3901  END DO
3902 
3903  CALL dbcsr_release(work)
3904  CALL dbcsr_release(work2)
3905  CALL dbcsr_release(prn)
3906  CALL dbcsr_release(a2n)
3907  CALL dbcsr_release(yk)
3908 
3909  CALL timestop(handle)
3910 
3911  END SUBROUTINE build_y_matrix
3912 
3913 ! **************************************************************************************************
3914 !> \brief Overwrites the "optimal" Laplace quadrature with that of the first step
3915 !> \param tj ...
3916 !> \param wj ...
3917 !> \param tau_tj ...
3918 !> \param tau_wj ...
3919 !> \param weights_cos_tf_t_to_w ...
3920 !> \param weights_cos_tf_w_to_t ...
3921 !> \param do_laplace ...
3922 !> \param do_im_time ...
3923 !> \param num_integ_points ...
3924 !> \param unit_nr ...
3925 !> \param qs_env ...
3926 ! **************************************************************************************************
3927  SUBROUTINE keep_initial_quad(tj, wj, tau_tj, tau_wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, &
3928  do_laplace, do_im_time, num_integ_points, unit_nr, qs_env)
3929 
3930  REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: tj, wj, tau_tj, tau_wj
3931  REAL(dp), ALLOCATABLE, DIMENSION(:, :), &
3932  INTENT(INOUT) :: weights_cos_tf_t_to_w, &
3933  weights_cos_tf_w_to_t
3934  LOGICAL, INTENT(IN) :: do_laplace, do_im_time
3935  INTEGER, INTENT(IN) :: num_integ_points, unit_nr
3936  TYPE(qs_environment_type), POINTER :: qs_env
3937 
3938  INTEGER :: jquad
3939 
3940  IF (do_laplace .OR. do_im_time) THEN
3941  IF (.NOT. ASSOCIATED(qs_env%mp2_env%ri_rpa_im_time%tau_tj)) THEN
3942  ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%tau_tj(0:num_integ_points))
3943  ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%tau_wj(num_integ_points))
3944  qs_env%mp2_env%ri_rpa_im_time%tau_tj(:) = tau_tj(:)
3945  qs_env%mp2_env%ri_rpa_im_time%tau_wj(:) = tau_wj(:)
3946  ELSE
3947  !If weights already stored, we overwrite the new ones
3948  tau_tj(:) = qs_env%mp2_env%ri_rpa_im_time%tau_tj(:)
3949  tau_wj(:) = qs_env%mp2_env%ri_rpa_im_time%tau_wj(:)
3950  END IF
3951  END IF
3952  IF (.NOT. do_laplace) THEN
3953  IF (.NOT. ASSOCIATED(qs_env%mp2_env%ri_rpa_im_time%tj)) THEN
3954  ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%tj(num_integ_points))
3955  ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%wj(num_integ_points))
3956  qs_env%mp2_env%ri_rpa_im_time%tj(:) = tj(:)
3957  qs_env%mp2_env%ri_rpa_im_time%wj(:) = wj(:)
3958  IF (do_im_time) THEN
3959  ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
3960  ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_w_to_t(num_integ_points, num_integ_points))
3961  qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_t_to_w(:, :) = weights_cos_tf_t_to_w(:, :)
3962  qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_w_to_t(:, :) = weights_cos_tf_w_to_t(:, :)
3963  END IF
3964  ELSE
3965  tj(:) = qs_env%mp2_env%ri_rpa_im_time%tj(:)
3966  wj(:) = qs_env%mp2_env%ri_rpa_im_time%wj(:)
3967  IF (do_im_time) THEN
3968  weights_cos_tf_t_to_w(:, :) = qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_t_to_w(:, :)
3969  weights_cos_tf_w_to_t(:, :) = qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_w_to_t(:, :)
3970  END IF
3971  END IF
3972  END IF
3973  IF (unit_nr > 0) THEN
3974  !Printing order same as in mp2_grids.F for consistency
3975  IF (ASSOCIATED(qs_env%mp2_env%ri_rpa_im_time%tj) .AND. (.NOT. do_laplace)) THEN
3976  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
3977  "MINIMAX_INFO| Number of integration points:", num_integ_points
3978  WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") &
3979  "MINIMAX_INFO| Minimax params (freq grid, scaled):", "Weights", "Abscissas"
3980  DO jquad = 1, num_integ_points
3981  WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") wj(jquad), tj(jquad)
3982  END DO
3983  CALL m_flush(unit_nr)
3984  END IF
3985  IF (ASSOCIATED(qs_env%mp2_env%ri_rpa_im_time%tau_tj)) THEN
3986  WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
3987  "MINIMAX_INFO| Number of integration points:", num_integ_points
3988  WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") &
3989  "MINIMAX_INFO| Minimax params (time grid, scaled):", "Weights", "Abscissas"
3990  DO jquad = 1, num_integ_points
3991  WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") tau_wj(jquad), tau_tj(jquad)
3992  END DO
3993  CALL m_flush(unit_nr)
3994  END IF
3995  END IF
3996 
3997  END SUBROUTINE keep_initial_quad
3998 
3999 END MODULE rpa_im_time_force_methods
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
static GRID_HOST_DEVICE int ncoset(const int l)
Number of Cartesian orbitals up to given angular momentum quantum.
Definition: grid_common.h:73
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Contains ADMM methods which require molecular orbitals.
Definition: admm_methods.F:15
subroutine, public admm_projection_derivative(qs_env, matrix_hz, matrix_pz, fval)
Calculate derivatives terms from overlap matrices.
Types and set/get functions for auxiliary density matrix methods.
Definition: admm_types.F:15
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition: admm_types.F:593
All kind of helpful little routines.
Definition: ao_util.F:14
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Definition: ao_util.F:209
Define the atomic kind types and their sub types.
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.
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
integer, save, public bussy2023
Definition: bibliography.F:43
Handles all functions related to the CELL.
Definition: cell_types.F:15
Calculation of the nuclear attraction contribution to the core Hamiltonian <a|erfc|b> :we only calcul...
Definition: core_ae.F:14
subroutine, public build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index)
...
Definition: core_ae.F:84
Calculation of the local pseudopotential contribution to the core Hamiltonian <a|V(local)|b> = <a|Sum...
Definition: core_ppl.F:18
subroutine, public build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, nimages, cell_to_index, basis_type, deltaR)
...
Definition: core_ppl.F:95
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition: core_ppnl.F:15
subroutine, public build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, basis_type, deltaR, matrix_l)
...
Definition: core_ppnl.F:89
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Interface to (sca)lapack for the Cholesky based procedures.
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,...
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, upper_to_full)
used to replace the cholesky decomposition by the inverse
Interface to (sca)lapack for the Cholesky based procedures.
Definition: cp_dbcsr_diag.F:17
subroutine, public cp_dbcsr_power(matrix, exponent, threshold, n_dependent, para_env, blacs_env, verbose, eigenvectors, eigenvalues)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public cp_dbcsr_dist2d_to_dist(dist2d, dist)
Creates a DBCSR distribution from a distribution_2d.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Interface to Minimax-Ewald method for periodic ERI's to be used in CP2K.
subroutine, public cp_eri_mme_update_local_counts(param, para_env, G_count_2c, R_count_2c, GG_count_3c, GR_count_3c, RR_count_3c)
Update local counters to gather statistics on different paths taken in MME algorithm (each Ewald sum ...
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
This is the start of a dbt_api, all publically needed functions are exported here....
Definition: dbt_api.F:17
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
Routines used for Harris functional Kohn-Sham calculation.
Definition: ec_methods.F:15
subroutine, public create_kernel(qs_env, vxc, vxc_tau, rho, rho1_r, rho1_g, tau1_r, xc_section, compute_virial, virial_xc)
Creation of second derivative xc-potential.
Definition: ec_methods.F:88
integer function, public gaussian_gridlevel(gridlevel_info, exponent)
...
Utilities for hfx and admm methods.
subroutine, public tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, external_hfx_sections, external_x_data)
Add the hfx contributions to the Hamiltonian.
Routines to calculate derivatives with respect to basis function origin.
subroutine, public derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_section, para_env, irep, use_virial, adiabatic_rescale_factor, resp_only, external_x_data)
computes four center derivatives for a full basis set and updates the forcesfock_4c arrays....
Routines to calculate EXX in RPA and energy correction methods.
Definition: hfx_exx.F:16
subroutine, public add_exx_to_rhs(rhs, qs_env, ext_hfx_section, x_data, recalc_integrals, do_admm, do_ec, do_exx, reuse_hfx)
Add the EXX contribution to the RHS of the Z-vector equation, namely the HF Hamiltonian.
Definition: hfx_exx.F:325
RI-methods for HFX.
Definition: hfx_ri.F:12
subroutine, public get_force_from_3c_trace(force, t_3c_contr, t_3c_der, atom_of_kind, kind_of, idx_to_at, pref, do_mp2, deriv_dim)
This routines calculates the force contribution from a trace over 3D tensors, i.e....
Definition: hfx_ri.F:3360
subroutine, public get_idx_to_atom(idx_to_at, bsizes_split, bsizes_orig)
a small utility function that returns the atom corresponding to a block of a split tensor
Definition: hfx_ri.F:4106
subroutine, public hfx_ri_update_forces(qs_env, ri_data, nspins, hf_fraction, rho_ao, rho_ao_resp, mos, use_virial, resp_only, rescale_factor)
the general routine that calls the relevant force code
Definition: hfx_ri.F:3036
subroutine, public get_2c_der_force(force, t_2c_contr, t_2c_der, atom_of_kind, kind_of, idx_to_at, pref, do_mp2, do_ovlp)
Update the forces due to the derivative of the a 2-center product d/dR (Q|R)
Definition: hfx_ri.F:3446
Types and set/get functions for HFX.
Definition: hfx_types.F:15
subroutine, public dealloc_containers(DATA, memory_usage)
...
Definition: hfx_types.F:2874
subroutine, public alloc_containers(DATA, bin_size)
...
Definition: hfx_types.F:2906
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_eri_mme
integer, parameter, public ri_rpa_method_gpw
integer, parameter, public do_admm_aux_exch_func_none
integer, parameter, public do_potential_id
integer, parameter, public do_eri_gpw
objects that represent the structure of input sections and the data contained in an input section
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
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
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
Routines useful for iterative matrix calculations.
subroutine, public matrix_exponential(matrix_exp, matrix, omega, alpha, threshold)
...
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public int_8
Definition: kinds.F:54
integer, parameter, public dp
Definition: kinds.F:34
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition: libint_2c_3c.F:14
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition: machine.F:123
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public fourpi
Interface to the message passing library MPI.
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Routines to calculate 2c- and 3c-integrals for RI with GPW.
Definition: mp2_eri_gpw.F:13
subroutine, public cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L)
Cleanup GPW integration for RI-MP2/RI-RPA.
Definition: mp2_eri_gpw.F:1060
subroutine, public virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, potential_parameter, para_env_sub)
Calculates stress tensor contribution from the operator.
Definition: mp2_eri_gpw.F:721
subroutine, public prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub)
Prepares GPW calculation for RI-MP2/RI-RPA.
Definition: mp2_eri_gpw.F:967
subroutine, public calc_potential_gpw(pot_r, rho_g, poisson_env, pot_g, potential_parameter, dvg, no_transfer)
Calculates potential from a given density in g-space.
Definition: mp2_eri_gpw.F:1110
Interface to direct methods for electron repulsion integrals for MP2.
Definition: mp2_eri.F:12
subroutine, public integrate_set_2c(param, potential_parameter, la_min, la_max, lb_min, lb_max, npgfa, npgfb, zeta, zetb, ra, rb, hab, n_hab_a, n_hab_b, offset_hab_a, offset_hab_b, offset_set_a, offset_set_b, sphi_a, sphi_b, sgfa, sgfb, nsgfa, nsgfb, eri_method, pab, force_a, force_b, hdab, hadb, G_count, R_count, do_reflection_a, do_reflection_b)
Integrate set pair and contract with sphi matrix.
Definition: mp2_eri.F:426
Types needed for MP2 calculations.
Definition: mp2_types.F:14
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
basic linear algebra operations for full matrixes
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
container for various plainwaves related things
Definition: pw_env_types.F:14
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
functions related to the poisson solver on regular grids
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type, external_vector)
maps a given wavefunction on the grid
collects routines that calculate density matrices
subroutine, public calculate_whz_matrix(c0vec, hzm, w_matrix, focc, nocc)
Calculate the Wz matrix from the MO eigenvectors, MO eigenvalues, and the MO occupation numbers....
Utility subroutine for qs energy calculation.
subroutine, public qs_energies_compute_w(qs_env, calc_forces)
Refactoring of qs_energies_scf. Moves computation of matrix_w into separate subroutine.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
Some utility functions for the calculation of integrals.
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.
Integrate single or product functions over a potential on a RS grid.
Calculate the interaction radii for the operator matrix calculation.
subroutine, public init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_short)
...
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
Calculation of kinetic energy matrix and forces.
Definition: qs_kinetic.F:15
subroutine, public build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, nderivative)
Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
Definition: qs_kinetic.F:101
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition: qs_ks_methods.F:22
subroutine, public calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density)
...
Calculate the KS reference potentials.
subroutine, public ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress)
calculate the Kohn-Sham reference potential
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition: qs_ks_types.F:567
Type definitiona for linear response calculations.
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
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
Define the neighbor list data types and the corresponding functionality.
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Calculation of overlap matrix, its derivatives and forces.
Definition: qs_overlap.F:19
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition: qs_overlap.F:120
Utility functions for the perturbation calculations.
subroutine, public p_env_psi0_changed(p_env, qs_env)
To be called after the value of psi0 has changed. Recalculates the quantities S_psi0 and m_epsilon.
subroutine, public p_env_create(p_env, qs_env, p1_option, p1_admm_option, orthogonal_orbitals, linres_control)
allocates and initializes the perturbation environment (no setup)
basis types for the calculation of the perturbation of density theory.
subroutine, public p_env_release(p_env)
relases the given p_env (see doc/ReferenceCounting.html)
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
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
Utility methods to build 3-center integral tensors of various types.
subroutine, public distribution_3d_create(dist_3d, dist1, dist2, dist3, nkind, particle_set, mp_comm_3d, own_comm)
Create a 3d distribution.
subroutine, public create_2c_tensor(t2c, dist_1, dist_2, pgrid, sizes_1, sizes_2, order, name)
...
subroutine, public create_tensor_batches(sizes, nbatches, starts_array, ends_array, starts_array_block, ends_array_block)
...
subroutine, public create_3c_tensor(t3c, dist_1, dist_2, dist_3, pgrid, sizes_1, sizes_2, sizes_3, map1, map2, name)
...
Utility methods to build 3-center integral tensors of various types.
Definition: qs_tensors.F:11
subroutine, public calc_3c_virial(work_virial, t3c_trace, pref, qs_env, nl_3c, basis_i, basis_j, basis_k, potential_parameter, der_eps, op_pos)
Calculates the 3c virial contributions on the fly.
Definition: qs_tensors.F:1868
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
subroutine, public build_3c_derivatives(t3c_der_i, t3c_der_k, filter_eps, qs_env, nl_3c, basis_i, basis_j, basis_k, potential_parameter, der_eps, op_pos, do_kpoints, do_hfx_kpoints, bounds_i, bounds_j, bounds_k, RI_range, img_to_RI_cell)
Build 3-center derivative tensors.
Definition: qs_tensors.F:1147
subroutine, public compress_tensor(tensor, blk_indices, compressed, eps, memory)
...
Definition: qs_tensors.F:3666
subroutine, public neighbor_list_3c_destroy(ijk_list)
Destroy 3c neighborlist.
Definition: qs_tensors.F:383
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:3232
subroutine, public decompress_tensor(tensor, blk_indices, compressed, eps)
...
Definition: qs_tensors.F:3794
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:2989
subroutine, public get_tensor_occupancy(tensor, nze, occ)
...
Definition: qs_tensors.F:3885
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
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:3412
pure logical function, public map_gaussian_here(rs_grid, h_inv, ra, offset, group_size, my_pos)
...
Calculate the CPKS equation and the resulting forces.
subroutine, public response_equation_new(qs_env, p_env, cpmos, iounit)
Initializes vectors for MO-coefficient based linear response solver and calculates response density,...
Routines needed for cubic-scaling RPA and SOS-Laplace-MP2 forces.
subroutine, public calc_rpa_loop_forces(force_data, mat_P_omega, t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, num_integ_points, nmo, Eigenval, e_fermi, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, tj, wj, tau_tj, cut_memory, ispin, open_shell, unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
Updates the cubic-scaling RPA contribution to the forces at each quadrature point....
subroutine, public calc_laplace_loop_forces(force_data, mat_P_omega, t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, num_integ_points, nmo, Eigenval, tau_tj, tau_wj, cut_memory, Pspin, Qspin, open_shell, unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
Updates the cubic-scaling SOS-Laplace-MP2 contribution to the forces at each quadrature point.
subroutine, public calc_post_loop_forces(force_data, unit_nr, qs_env)
All the forces that can be calculated after the loop on the Laplace quaradture, using terms collected...
subroutine, public init_im_time_forces(force_data, fm_matrix_PQ, t_3c_M, unit_nr, mp2_env, qs_env)
Initializes and pre-calculates all needed tensors for the forces.
subroutine, public keep_initial_quad(tj, wj, tau_tj, tau_wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, do_laplace, do_im_time, num_integ_points, unit_nr, qs_env)
Overwrites the "optimal" Laplace quadrature with that of the first step.
Types needed for cubic-scaling RPA and SOS-Laplace-MP2 forces.
Routines for low-scaling RPA/GW with imaginary time.
Definition: rpa_im_time.F:13
subroutine, public compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, tau_tj, num_integ_points, nmo, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm_occ_global, mat_dm_virt_global, matrix_s, ispin, Eigenval, e_fermi, eps_filter, memory_info, unit_nr, jquad, do_kpoints_cubic_RPA, do_kpoints_from_Gamma, qs_env, num_cells_dm, index_to_cell_dm, para_env)
...
Definition: rpa_im_time.F:665
Transfers densities from PW to RS grids and potentials from PW to RS.
subroutine, public potential_pw2rs(rs_v, v_rspace, pw_env)
transfers a potential from a pw_grid to a vector of realspace multigrids
types for task lists