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