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