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