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