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