(git:a535875)
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-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines needed for cubic-scaling RPA and SOS-Laplace-MP2 forces
10!> \author Augustin Bussy
11! **************************************************************************************************
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 cp_blacs_env, ONLY: cp_blacs_env_type
26 USE cp_control_types, ONLY: dft_control_type
27 USE cp_dbcsr_api, ONLY: &
28 dbcsr_add, dbcsr_clear, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, &
29 dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, &
30 dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
31 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
32 dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
33 dbcsr_type_no_symmetry, dbcsr_type_symmetric
34 USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
35 cp_dbcsr_cholesky_invert
36 USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
37 dbcsr_frobenius_norm
38 USE cp_dbcsr_diag, ONLY: cp_dbcsr_power
39 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
40 copy_fm_to_dbcsr,&
41 cp_dbcsr_dist2d_to_dist,&
42 cp_dbcsr_sm_fm_multiply,&
43 dbcsr_allocate_matrix_set,&
44 dbcsr_deallocate_matrix_set
45 USE cp_eri_mme_interface, ONLY: cp_eri_mme_update_local_counts
46 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
47 cp_fm_struct_release,&
48 cp_fm_struct_type
49 USE cp_fm_types, ONLY: cp_fm_create,&
50 cp_fm_release,&
51 cp_fm_set_all,&
52 cp_fm_to_fm,&
53 cp_fm_type
54 USE dbt_api, ONLY: &
55 dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_clear, dbt_contract, &
56 dbt_copy, dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, dbt_destroy, &
57 dbt_filter, dbt_get_info, dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, &
58 dbt_pgrid_type, dbt_scale, dbt_type
59 USE distribution_2d_types, ONLY: distribution_2d_type
60 USE ec_methods, ONLY: create_kernel
61 USE gaussian_gridlevels, ONLY: gaussian_gridlevel
62 USE hfx_admm_utils, ONLY: tddft_hfx_matrix
63 USE hfx_derivatives, ONLY: derivatives_four_center
64 USE hfx_exx, ONLY: add_exx_to_rhs
65 USE hfx_ri, ONLY: get_2c_der_force,&
66 get_force_from_3c_trace,&
67 get_idx_to_atom,&
68 hfx_ri_update_forces
69 USE hfx_types, ONLY: alloc_containers,&
70 block_ind_type,&
71 dealloc_containers,&
72 hfx_compression_type,&
73 hfx_type
74 USE input_constants, ONLY: do_admm_aux_exch_func_none,&
75 do_eri_gpw,&
76 do_eri_mme,&
77 do_potential_id,&
78 ri_rpa_method_gpw
79 USE input_section_types, ONLY: section_vals_get,&
80 section_vals_get_subs_vals,&
81 section_vals_type,&
82 section_vals_val_get
83 USE iterate_matrix, ONLY: matrix_exponential
84 USE kinds, ONLY: dp,&
85 int_8
86 USE libint_2c_3c, ONLY: libint_potential_type
87 USE machine, ONLY: m_flush,&
88 m_walltime
89 USE mathconstants, ONLY: fourpi
90 USE message_passing, ONLY: mp_cart_type,&
91 mp_para_env_release,&
92 mp_para_env_type
93 USE mp2_eri, ONLY: integrate_set_2c
94 USE mp2_eri_gpw, ONLY: calc_potential_gpw,&
95 cleanup_gpw,&
96 prepare_gpw,&
97 virial_gpw_potential
98 USE mp2_types, ONLY: mp2_type
99 USE orbital_pointers, ONLY: ncoset
100 USE parallel_gemm_api, ONLY: parallel_gemm
101 USE particle_methods, ONLY: get_particle_set
102 USE particle_types, ONLY: particle_type
103 USE pw_env_types, ONLY: pw_env_get,&
104 pw_env_type
105 USE pw_methods, ONLY: pw_axpy,&
106 pw_copy,&
107 pw_integral_ab,&
108 pw_scale,&
109 pw_transfer,&
110 pw_zero
111 USE pw_poisson_methods, ONLY: pw_poisson_solve
112 USE pw_poisson_types, ONLY: pw_poisson_type
113 USE pw_pool_types, ONLY: pw_pool_type
114 USE pw_types, ONLY: pw_c1d_gs_type,&
115 pw_r3d_rs_type
116 USE qs_collocate_density, ONLY: calculate_rho_elec,&
117 collocate_function
118 USE qs_core_matrices, ONLY: core_matrices,&
119 kinetic_energy_matrix
120 USE qs_density_matrices, ONLY: calculate_whz_matrix
121 USE qs_environment_types, ONLY: get_qs_env,&
122 qs_environment_type,&
123 set_qs_env
124 USE qs_force_types, ONLY: qs_force_type
125 USE qs_integral_utils, ONLY: basis_set_list_setup
126 USE qs_integrate_potential, ONLY: integrate_pgf_product,&
127 integrate_v_core_rspace,&
128 integrate_v_rspace
129 USE qs_interactions, ONLY: init_interaction_radii_orb_basis
130 USE qs_kind_types, ONLY: qs_kind_type
131 USE qs_ks_methods, ONLY: calc_rho_tot_gspace
132 USE qs_ks_reference, ONLY: ks_ref_potential
133 USE qs_ks_types, ONLY: set_ks_env
134 USE qs_linres_types, ONLY: linres_control_type
135 USE qs_matrix_w, ONLY: compute_matrix_w
136 USE qs_mo_types, ONLY: get_mo_set,&
137 mo_set_type
138 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
139 release_neighbor_list_sets
140 USE qs_overlap, ONLY: build_overlap_matrix
141 USE qs_p_env_methods, ONLY: p_env_create,&
142 p_env_psi0_changed
143 USE qs_p_env_types, ONLY: p_env_release,&
144 qs_p_env_type
145 USE qs_rho_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, uplo_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, uplo_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(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(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, uplo_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 with 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 IF (cut_memory > 0) CALL dbt_batched_contract_init(t_kbkt)
2108 DO j_mem = 1, cut_memory
2109 jbounds(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
2110
2111 bounds_cpy(:, 1) = [1, bounds_3c(1)]
2112 bounds_cpy(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
2113 bounds_cpy(:, 3) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
2114 CALL dbt_copy(t_3c_sparse, t_3c_7, bounds=bounds_cpy)
2115
2116 CALL dbt_batched_contract_init(t_dm_virt)
2117 DO k_mem = 1, n_mem_ri
2118 bounds_2c(:, 1) = [batch_start_ri(k_mem), batch_end_ri(k_mem)]
2119 bounds_2c(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
2120
2121 CALL timeset(routinen//"_3c_dm", handle2)
2122
2123 !Calculate (mu nu| P) * D_occ * D_virt
2124 !Note: technically need M_occ*D_virt + M_virt*D_occ, but it is equivalent to 2*M_occ*D_virt
2125 CALL dbt_contract(2.0_dp, t_3c_4, t_dm_virt, 0.0_dp, t_3c_5, &
2126 contract_1=[3], notcontract_1=[1, 2], &
2127 contract_2=[1], notcontract_2=[2], &
2128 map_1=[1, 2], map_2=[3], filter_eps=eps_filter, &
2129 bounds_2=bounds_2c, bounds_3=jbounds, flop=flop, unit_nr=unit_nr_dbcsr)
2130 dbcsr_nflop = dbcsr_nflop + flop
2131
2132 CALL get_tensor_occupancy(t_3c_5, nze, occ)
2133 nze_ddint = nze_ddint + nze
2134 occ_ddint = occ_ddint + occ
2135
2136 ! Skip the expensive KBK^T contraction when the intermediate block is empty
2137 IF (nze == 0) THEN
2138 CALL dbt_clear(t_3c_5)
2139 cycle
2140 END IF
2141
2142 CALL dbt_copy(t_3c_5, t_3c_6, move_data=.true.)
2143 CALL timestop(handle2)
2144
2145 !Calculate the contraction of the above with K*B*K^T
2146 CALL timeset(routinen//"_3c_KBK", handle2)
2147 CALL dbt_contract(1.0_dp, t_kbkt, t_3c_6, 0.0_dp, t_3c_7, &
2148 contract_1=[2], notcontract_1=[1], &
2149 contract_2=[1], notcontract_2=[2, 3], &
2150 map_1=[1], map_2=[2, 3], &
2151 retain_sparsity=.true., flop=flop, unit_nr=unit_nr_dbcsr)
2152 dbcsr_nflop = dbcsr_nflop + flop
2153 CALL timestop(handle2)
2154 CALL dbt_copy(t_3c_7, t_3c_8, summation=.true.)
2155
2156 END DO !k_mem
2157 CALL dbt_batched_contract_finalize(t_dm_virt)
2158 END DO !j_mem
2159 IF (cut_memory > 0) CALL dbt_batched_contract_finalize(t_kbkt)
2160
2161 CALL dbt_copy(t_3c_8, t_3c_help_1, move_data=.true.)
2162
2163 pref = 1.0_dp*fac
2164 DO k_mem = 1, cut_memory
2165 DO i_xyz = 1, 3
2166 CALL dbt_clear(force_data%t_3c_der_RI(i_xyz))
2167 CALL decompress_tensor(force_data%t_3c_der_RI(i_xyz), force_data%t_3c_der_RI_ind(k_mem, i_xyz)%ind, &
2168 force_data%t_3c_der_RI_comp(k_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
2169 END DO
2170 CALL get_force_from_3c_trace(force, t_3c_help_1, force_data%t_3c_der_RI, atom_of_kind, kind_of, &
2171 force_data%idx_to_at_RI, pref, do_mp2=.true., deriv_dim=1)
2172 END DO
2173
2174 IF (use_virial) THEN
2175 CALL dbt_copy(t_3c_help_1, t_3c_help_2)
2176 CALL dbt_scale(t_3c_help_2, pref)
2177 CALL dbt_copy(t_3c_help_2, force_data%t_3c_virial_split, summation=.true., move_data=.true.)
2178 END IF
2179
2180 CALL dbt_copy(t_3c_help_1, t_3c_help_2)
2181 CALL dbt_copy(t_3c_help_1, t_3c_help_2, order=[1, 3, 2], move_data=.true., summation=.true.)
2182 DO k_mem = 1, cut_memory
2183 DO i_xyz = 1, 3
2184 CALL dbt_clear(force_data%t_3c_der_AO(i_xyz))
2185 CALL decompress_tensor(force_data%t_3c_der_AO(i_xyz), force_data%t_3c_der_AO_ind(k_mem, i_xyz)%ind, &
2186 force_data%t_3c_der_AO_comp(k_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
2187 END DO
2188 CALL get_force_from_3c_trace(force, t_3c_help_2, force_data%t_3c_der_AO, atom_of_kind, kind_of, &
2189 force_data%idx_to_at_AO, pref, do_mp2=.true., deriv_dim=3)
2190 END DO
2191
2192 CALL dbt_clear(t_3c_help_2)
2193 END DO !i_mem
2194 CALL dbt_batched_contract_finalize(t_r_occ)
2195 CALL dbt_batched_contract_finalize(t_r_virt)
2196
2197 DO k_mem = 1, n_mem_ri
2198 DO i_mem = 1, cut_memory
2199 CALL dealloc_containers(store_3c(k_mem, i_mem), dummy_int)
2200 END DO
2201 END DO
2202 DEALLOCATE (store_3c, blk_indices)
2203
2204 CALL timestop(handle)
2205
2206 END SUBROUTINE perform_3c_ops
2207
2208! **************************************************************************************************
2209!> \brief All the forces that can be calculated after the loop on the Laplace quaradture, using
2210!> terms collected during the said loop. This inludes the z-vector equation and its reponse
2211!> forces, as well as the force coming from the trace with the derivative of the KS matrix
2212!> \param force_data ...
2213!> \param unit_nr ...
2214!> \param qs_env ...
2215! **************************************************************************************************
2216 SUBROUTINE calc_post_loop_forces(force_data, unit_nr, qs_env)
2217
2218 TYPE(im_time_force_type), INTENT(INOUT) :: force_data
2219 INTEGER, INTENT(IN) :: unit_nr
2220 TYPE(qs_environment_type), POINTER :: qs_env
2221
2222 CHARACTER(len=*), PARAMETER :: routinen = 'calc_post_loop_forces'
2223
2224 INTEGER :: handle, ispin, nao, nao_aux, nocc, nspins
2225 LOGICAL :: do_exx
2226 REAL(dp) :: focc
2227 TYPE(admm_type), POINTER :: admm_env
2228 TYPE(cp_fm_struct_type), POINTER :: fm_struct
2229 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: cpmos, mo_occ
2230 TYPE(cp_fm_type), POINTER :: mo_coeff
2231 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_p_work, matrix_p_mp2, &
2232 matrix_p_mp2_admm, matrix_s, &
2233 matrix_s_aux, work_admm, yp_admm
2234 TYPE(dft_control_type), POINTER :: dft_control
2235 TYPE(linres_control_type), POINTER :: linres_control
2236 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2237 TYPE(qs_p_env_type), POINTER :: p_env
2238 TYPE(section_vals_type), POINTER :: hfx_section, lr_section
2239
2240 NULLIFY (linres_control, p_env, dft_control, matrix_s, mos, mo_coeff, fm_struct, lr_section, &
2241 dbcsr_p_work, yp_admm, matrix_p_mp2, admm_env, work_admm, matrix_s_aux, matrix_p_mp2_admm)
2242
2243 CALL timeset(routinen, handle)
2244
2245 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_s=matrix_s, mos=mos)
2246 nspins = dft_control%nspins
2247
2248 ! Setting up for the z-vector equation
2249
2250 ! Initialize linres_control
2251 lr_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%LOW_SCALING%CPHF")
2252
2253 ALLOCATE (linres_control)
2254 CALL section_vals_val_get(lr_section, "MAX_ITER", i_val=linres_control%max_iter)
2255 CALL section_vals_val_get(lr_section, "EPS_CONV", r_val=linres_control%eps)
2256 CALL section_vals_val_get(lr_section, "PRECONDITIONER", i_val=linres_control%preconditioner_type)
2257 CALL section_vals_val_get(lr_section, "ENERGY_GAP", r_val=linres_control%energy_gap)
2258
2259 linres_control%do_kernel = .true.
2260 linres_control%lr_triplet = .false.
2261 linres_control%converged = .false.
2262 linres_control%eps_filter = qs_env%mp2_env%ri_rpa_im_time%eps_filter
2263
2264 CALL set_qs_env(qs_env, linres_control=linres_control)
2265
2266 IF (unit_nr > 0) THEN
2267 WRITE (unit_nr, *)
2268 WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations'
2269 WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Convergence threshold:', linres_control%eps
2270 WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Maximum number of iterations: ', linres_control%max_iter
2271 END IF
2272
2273 ALLOCATE (p_env)
2274 CALL p_env_create(p_env, qs_env, orthogonal_orbitals=.true., linres_control=linres_control)
2275 CALL p_env_psi0_changed(p_env, qs_env)
2276
2277 ! Matrix allocation
2278 CALL dbcsr_allocate_matrix_set(p_env%p1, nspins)
2279 CALL dbcsr_allocate_matrix_set(p_env%w1, nspins)
2280 CALL dbcsr_allocate_matrix_set(dbcsr_p_work, nspins)
2281 DO ispin = 1, nspins
2282 ALLOCATE (p_env%p1(ispin)%matrix, p_env%w1(ispin)%matrix, dbcsr_p_work(ispin)%matrix)
2283 CALL dbcsr_create(matrix=p_env%p1(ispin)%matrix, template=matrix_s(1)%matrix)
2284 CALL dbcsr_create(matrix=p_env%w1(ispin)%matrix, template=matrix_s(1)%matrix)
2285 CALL dbcsr_create(matrix=dbcsr_p_work(ispin)%matrix, template=matrix_s(1)%matrix)
2286 CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_s(1)%matrix)
2287 CALL dbcsr_copy(p_env%w1(ispin)%matrix, matrix_s(1)%matrix)
2288 CALL dbcsr_copy(dbcsr_p_work(ispin)%matrix, matrix_s(1)%matrix)
2289 CALL dbcsr_set(p_env%p1(ispin)%matrix, 0.0_dp)
2290 CALL dbcsr_set(p_env%w1(ispin)%matrix, 0.0_dp)
2291 CALL dbcsr_set(dbcsr_p_work(ispin)%matrix, 0.0_dp)
2292 END DO
2293
2294 IF (dft_control%do_admm) THEN
2295 CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux)
2296 CALL dbcsr_allocate_matrix_set(p_env%p1_admm, nspins)
2297 CALL dbcsr_allocate_matrix_set(work_admm, nspins)
2298 DO ispin = 1, nspins
2299 ALLOCATE (p_env%p1_admm(ispin)%matrix, work_admm(ispin)%matrix)
2300 CALL dbcsr_create(p_env%p1_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2301 CALL dbcsr_copy(p_env%p1_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2302 CALL dbcsr_set(p_env%p1_admm(ispin)%matrix, 0.0_dp)
2303 CALL dbcsr_create(work_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2304 CALL dbcsr_copy(work_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2305 CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp)
2306 END DO
2307 END IF
2308
2309 ! Preparing the RHS of the z-vector equation
2310 CALL prepare_for_response(force_data, qs_env)
2311 ALLOCATE (cpmos(nspins), mo_occ(nspins))
2312 DO ispin = 1, nspins
2313 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, homo=nocc)
2314 NULLIFY (fm_struct)
2315 CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, &
2316 template_fmstruct=mo_coeff%matrix_struct)
2317 CALL cp_fm_create(cpmos(ispin), fm_struct)
2318 CALL cp_fm_set_all(cpmos(ispin), 0.0_dp)
2319 CALL cp_fm_create(mo_occ(ispin), fm_struct)
2320 CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc)
2321 CALL cp_fm_struct_release(fm_struct)
2322 END DO
2323
2324 ! in case of EXX, need to add the HF Hamiltonian to the RHS of the Z-vector equation
2325 ! Strategy: we take the ks_matrix, remove the current xc contribution, and then add the RPA HF one
2326 do_exx = .false.
2327 IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
2328 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
2329 CALL section_vals_get(hfx_section, explicit=do_exx)
2330 END IF
2331
2332 IF (do_exx) THEN
2333 CALL add_exx_to_rhs(rhs=force_data%sum_O_tau, &
2334 qs_env=qs_env, &
2335 ext_hfx_section=hfx_section, &
2336 x_data=qs_env%mp2_env%ri_rpa%x_data, &
2337 recalc_integrals=.false., &
2338 do_admm=qs_env%mp2_env%ri_rpa%do_admm, &
2339 do_exx=do_exx, &
2340 reuse_hfx=qs_env%mp2_env%ri_rpa%reuse_hfx)
2341 END IF
2342
2343 focc = 2.0_dp
2344 IF (nspins == 1) focc = 4.0_dp
2345 DO ispin = 1, nspins
2346 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
2347 CALL cp_dbcsr_sm_fm_multiply(force_data%sum_O_tau(ispin)%matrix, mo_occ(ispin), &
2348 cpmos(ispin), nocc, &
2349 alpha=focc, beta=0.0_dp)
2350 END DO
2351
2352 ! The z-vector equation and associated forces
2353 CALL response_equation_new(qs_env, p_env, cpmos, unit_nr)
2354
2355 ! Save the mp2 density matrix
2356 CALL get_qs_env(qs_env, matrix_p_mp2=matrix_p_mp2)
2357 IF (ASSOCIATED(matrix_p_mp2)) CALL dbcsr_deallocate_matrix_set(matrix_p_mp2)
2358 DO ispin = 1, nspins
2359 CALL dbcsr_copy(dbcsr_p_work(ispin)%matrix, p_env%p1(ispin)%matrix)
2360 CALL dbcsr_add(dbcsr_p_work(ispin)%matrix, force_data%sum_YP_tau(ispin)%matrix, 1.0_dp, 1.0_dp)
2361 END DO
2362 CALL set_ks_env(qs_env%ks_env, matrix_p_mp2=dbcsr_p_work)
2363
2364 IF (dft_control%do_admm) THEN
2365 CALL dbcsr_allocate_matrix_set(yp_admm, nspins)
2366 CALL get_qs_env(qs_env, matrix_p_mp2_admm=matrix_p_mp2_admm, admm_env=admm_env)
2367 nao = admm_env%nao_orb
2368 nao_aux = admm_env%nao_aux_fit
2369 IF (ASSOCIATED(matrix_p_mp2_admm)) CALL dbcsr_deallocate_matrix_set(matrix_p_mp2_admm)
2370 DO ispin = 1, nspins
2371
2372 !sum_YP_tau in the auxiliary basis
2373 CALL copy_dbcsr_to_fm(force_data%sum_YP_tau(ispin)%matrix, admm_env%work_orb_orb)
2374 CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, admm_env%work_orb_orb, &
2375 0.0_dp, admm_env%work_aux_orb)
2376 CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
2377 0.0_dp, admm_env%work_aux_aux)
2378 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_admm(ispin)%matrix, keep_sparsity=.true.)
2379
2380 !save the admm representation od sum_YP_tau
2381 ALLOCATE (yp_admm(ispin)%matrix)
2382 CALL dbcsr_create(yp_admm(ispin)%matrix, template=work_admm(ispin)%matrix)
2383 CALL dbcsr_copy(yp_admm(ispin)%matrix, work_admm(ispin)%matrix)
2384
2385 CALL dbcsr_add(work_admm(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
2386
2387 END DO
2388 CALL set_ks_env(qs_env%ks_env, matrix_p_mp2_admm=work_admm)
2389 END IF
2390
2391 !Calculate the response force and the force from the trace with F
2392 CALL update_im_time_forces(p_env, force_data%sum_O_tau, force_data%sum_YP_tau, yp_admm, qs_env)
2393
2394 !clean-up
2395 IF (dft_control%do_admm) CALL dbcsr_deallocate_matrix_set(yp_admm)
2396
2397 CALL cp_fm_release(cpmos)
2398 CALL cp_fm_release(mo_occ)
2399 CALL p_env_release(p_env)
2400 DEALLOCATE (p_env)
2401
2402 CALL timestop(handle)
2403
2404 END SUBROUTINE calc_post_loop_forces
2405
2406! **************************************************************************************************
2407!> \brief Prepares the RHS of the z-vector equation. Apply the xc and HFX kernel on the previously
2408!> stored sum_YP_tau density, and add it to the final force_data%sum_O_tau quantity
2409!> \param force_data ...
2410!> \param qs_env ...
2411! **************************************************************************************************
2412 SUBROUTINE prepare_for_response(force_data, qs_env)
2413
2414 TYPE(im_time_force_type), INTENT(INOUT) :: force_data
2415 TYPE(qs_environment_type), POINTER :: qs_env
2416
2417 CHARACTER(LEN=*), PARAMETER :: routinen = 'prepare_for_response'
2418
2419 INTEGER :: handle, ispin, nao, nao_aux, nspins
2420 LOGICAL :: do_hfx, do_tau, do_tau_admm
2421 REAL(dp) :: ehartree
2422 TYPE(admm_type), POINTER :: admm_env
2423 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_p_work, ker_tau_admm, matrix_s, &
2424 matrix_s_aux, work_admm
2425 TYPE(dbcsr_type) :: dbcsr_work
2426 TYPE(dft_control_type), POINTER :: dft_control
2427 TYPE(pw_c1d_gs_type) :: rhoz_tot_gspace, zv_hartree_gspace
2428 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rhoz_g
2429 TYPE(pw_env_type), POINTER :: pw_env
2430 TYPE(pw_poisson_type), POINTER :: poisson_env
2431 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2432 TYPE(pw_r3d_rs_type) :: zv_hartree_rspace
2433 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rhoz_r, tauz_r, v_xc, v_xc_tau
2434 TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
2435 TYPE(section_vals_type), POINTER :: hfx_section, xc_section
2436 TYPE(task_list_type), POINTER :: task_list_aux_fit
2437
2438 NULLIFY (pw_env, rhoz_r, rhoz_g, tauz_r, v_xc, v_xc_tau, &
2439 poisson_env, auxbas_pw_pool, dft_control, admm_env, xc_section, rho, rho_aux_fit, &
2440 task_list_aux_fit, ker_tau_admm, work_admm, dbcsr_p_work, matrix_s, hfx_section)
2441
2442 CALL timeset(routinen, handle)
2443
2444 CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env, rho=rho, matrix_s=matrix_s)
2445 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
2446 nspins = dft_control%nspins
2447
2448 CALL dbcsr_allocate_matrix_set(dbcsr_p_work, nspins)
2449 DO ispin = 1, nspins
2450 ALLOCATE (dbcsr_p_work(ispin)%matrix)
2451 CALL dbcsr_create(matrix=dbcsr_p_work(ispin)%matrix, template=matrix_s(1)%matrix)
2452 CALL dbcsr_copy(dbcsr_p_work(ispin)%matrix, matrix_s(1)%matrix)
2453 CALL dbcsr_set(dbcsr_p_work(ispin)%matrix, 0.0_dp)
2454 END DO
2455
2456 !Apply the kernel on the density saved in force_data%sum_YP_tau
2457 ALLOCATE (rhoz_r(nspins), rhoz_g(nspins))
2458 DO ispin = 1, nspins
2459 CALL auxbas_pw_pool%create_pw(rhoz_r(ispin))
2460 CALL auxbas_pw_pool%create_pw(rhoz_g(ispin))
2461 END DO
2462 CALL auxbas_pw_pool%create_pw(rhoz_tot_gspace)
2463 CALL auxbas_pw_pool%create_pw(zv_hartree_rspace)
2464 CALL auxbas_pw_pool%create_pw(zv_hartree_gspace)
2465
2466 CALL pw_zero(rhoz_tot_gspace)
2467 DO ispin = 1, nspins
2468 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=force_data%sum_YP_tau(ispin)%matrix, &
2469 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin))
2470 CALL pw_axpy(rhoz_g(ispin), rhoz_tot_gspace)
2471 END DO
2472
2473 CALL pw_poisson_solve(poisson_env, rhoz_tot_gspace, ehartree, &
2474 zv_hartree_gspace)
2475
2476 CALL pw_transfer(zv_hartree_gspace, zv_hartree_rspace)
2477 CALL pw_scale(zv_hartree_rspace, zv_hartree_rspace%pw_grid%dvol)
2478
2479 CALL qs_rho_get(rho, tau_r_valid=do_tau)
2480 IF (do_tau) THEN
2481 block
2482 TYPE(pw_c1d_gs_type) :: tauz_g
2483 ALLOCATE (tauz_r(nspins))
2484 CALL auxbas_pw_pool%create_pw(tauz_g)
2485 DO ispin = 1, nspins
2486 CALL auxbas_pw_pool%create_pw(tauz_r(ispin))
2487
2488 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=force_data%sum_YP_tau(ispin)%matrix, &
2489 rho=tauz_r(ispin), rho_gspace=tauz_g, compute_tau=.true.)
2490 END DO
2491 CALL auxbas_pw_pool%give_back_pw(tauz_g)
2492 END block
2493 END IF
2494
2495 IF (dft_control%do_admm) THEN
2496 CALL get_qs_env(qs_env, admm_env=admm_env)
2497 xc_section => admm_env%xc_section_primary
2498 ELSE
2499 xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
2500 END IF
2501
2502 !Primary XC kernel
2503 CALL create_kernel(qs_env, v_xc, v_xc_tau, rho, rhoz_r, rhoz_g, tauz_r, xc_section)
2504
2505 DO ispin = 1, nspins
2506 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2507 CALL pw_axpy(zv_hartree_rspace, v_xc(ispin))
2508 CALL integrate_v_rspace(qs_env=qs_env, &
2509 v_rspace=v_xc(ispin), &
2510 hmat=dbcsr_p_work(ispin), &
2511 calculate_forces=.false.)
2512
2513 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2514 END DO
2515 CALL auxbas_pw_pool%give_back_pw(rhoz_tot_gspace)
2516 CALL auxbas_pw_pool%give_back_pw(zv_hartree_rspace)
2517 CALL auxbas_pw_pool%give_back_pw(zv_hartree_gspace)
2518 DEALLOCATE (v_xc)
2519
2520 IF (do_tau) THEN
2521 DO ispin = 1, nspins
2522 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2523 CALL integrate_v_rspace(qs_env=qs_env, &
2524 v_rspace=v_xc_tau(ispin), &
2525 hmat=dbcsr_p_work(ispin), &
2526 compute_tau=.true., &
2527 calculate_forces=.false.)
2528 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2529 END DO
2530 DEALLOCATE (v_xc_tau)
2531 END IF
2532
2533 !Auxiliary xc kernel (admm)
2534 IF (dft_control%do_admm) THEN
2535 CALL get_qs_env(qs_env, admm_env=admm_env)
2536 CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux, &
2537 task_list_aux_fit=task_list_aux_fit, rho_aux_fit=rho_aux_fit)
2538
2539 CALL qs_rho_get(rho_aux_fit, tau_r_valid=do_tau_admm)
2540
2541 CALL dbcsr_allocate_matrix_set(work_admm, nspins)
2542 CALL dbcsr_allocate_matrix_set(ker_tau_admm, nspins)
2543 DO ispin = 1, nspins
2544 ALLOCATE (work_admm(ispin)%matrix, ker_tau_admm(ispin)%matrix)
2545 CALL dbcsr_create(work_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2546 CALL dbcsr_copy(work_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2547 CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp)
2548 CALL dbcsr_create(ker_tau_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2549 CALL dbcsr_copy(ker_tau_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2550 CALL dbcsr_set(ker_tau_admm(ispin)%matrix, 0.0_dp)
2551 END DO
2552
2553 !get the density in the auxuliary density
2554 cpassert(ASSOCIATED(admm_env%work_orb_orb))
2555 cpassert(ASSOCIATED(admm_env%work_aux_orb))
2556 cpassert(ASSOCIATED(admm_env%work_aux_aux))
2557 nao = admm_env%nao_orb
2558 nao_aux = admm_env%nao_aux_fit
2559 DO ispin = 1, nspins
2560 CALL copy_dbcsr_to_fm(force_data%sum_YP_tau(ispin)%matrix, admm_env%work_orb_orb)
2561 CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, admm_env%work_orb_orb, &
2562 0.0_dp, admm_env%work_aux_orb)
2563 CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
2564 0.0_dp, admm_env%work_aux_aux)
2565 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, ker_tau_admm(ispin)%matrix, keep_sparsity=.true.)
2566 END DO
2567
2568 IF (.NOT. qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
2569 DO ispin = 1, nspins
2570 CALL pw_zero(rhoz_r(ispin))
2571 CALL pw_zero(rhoz_g(ispin))
2572 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=ker_tau_admm(ispin)%matrix, &
2573 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin), &
2574 basis_type="AUX_FIT", task_list_external=task_list_aux_fit)
2575 END DO
2576
2577 IF (do_tau_admm) THEN
2578 block
2579 TYPE(pw_c1d_gs_type) :: tauz_g
2580 CALL auxbas_pw_pool%create_pw(tauz_g)
2581 DO ispin = 1, nspins
2582 CALL pw_zero(tauz_r(ispin))
2583 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=ker_tau_admm(ispin)%matrix, &
2584 rho=tauz_r(ispin), rho_gspace=tauz_g, &
2585 basis_type="AUX_FIT", task_list_external=task_list_aux_fit, &
2586 compute_tau=.true.)
2587 END DO
2588 CALL auxbas_pw_pool%give_back_pw(tauz_g)
2589 END block
2590 END IF
2591
2592 xc_section => admm_env%xc_section_aux
2593 CALL create_kernel(qs_env, v_xc, v_xc_tau, rho_aux_fit, rhoz_r, rhoz_g, tauz_r, xc_section)
2594
2595 DO ispin = 1, nspins
2596 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2597 CALL integrate_v_rspace(qs_env=qs_env, &
2598 v_rspace=v_xc(ispin), &
2599 hmat=work_admm(ispin), &
2600 calculate_forces=.false., &
2601 basis_type="AUX_FIT", &
2602 task_list_external=task_list_aux_fit)
2603 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2604 END DO
2605 DEALLOCATE (v_xc)
2606
2607 IF (do_tau_admm) THEN
2608 DO ispin = 1, nspins
2609 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2610 CALL integrate_v_rspace(qs_env=qs_env, &
2611 v_rspace=v_xc_tau(ispin), &
2612 hmat=work_admm(ispin), &
2613 calculate_forces=.false., &
2614 basis_type="AUX_FIT", &
2615 task_list_external=task_list_aux_fit, &
2616 compute_tau=.true.)
2617 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2618 END DO
2619 DEALLOCATE (v_xc_tau)
2620 END IF
2621 END IF !admm
2622 END IF
2623
2624 DO ispin = 1, nspins
2625 CALL auxbas_pw_pool%give_back_pw(rhoz_r(ispin))
2626 CALL auxbas_pw_pool%give_back_pw(rhoz_g(ispin))
2627 END DO
2628 DEALLOCATE (rhoz_r, rhoz_g)
2629
2630 IF (do_tau) THEN
2631 DO ispin = 1, nspins
2632 CALL auxbas_pw_pool%give_back_pw(tauz_r(ispin))
2633 END DO
2634 DEALLOCATE (tauz_r)
2635 END IF
2636
2637 !HFX kernel
2638 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
2639 CALL section_vals_get(hfx_section, explicit=do_hfx)
2640 IF (do_hfx) THEN
2641 IF (dft_control%do_admm) THEN
2642 CALL tddft_hfx_matrix(work_admm, ker_tau_admm, qs_env, .false., .false.)
2643
2644 !Going back to primary basis
2645 CALL dbcsr_create(dbcsr_work, template=dbcsr_p_work(1)%matrix)
2646 CALL dbcsr_copy(dbcsr_work, dbcsr_p_work(1)%matrix)
2647 CALL dbcsr_set(dbcsr_work, 0.0_dp)
2648 DO ispin = 1, nspins
2649 CALL copy_dbcsr_to_fm(work_admm(ispin)%matrix, admm_env%work_aux_aux)
2650 CALL parallel_gemm('N', 'N', nao_aux, nao, nao_aux, 1.0_dp, admm_env%work_aux_aux, admm_env%A, &
2651 0.0_dp, admm_env%work_aux_orb)
2652 CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, admm_env%work_aux_orb, &
2653 0.0_dp, admm_env%work_orb_orb)
2654 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_work, keep_sparsity=.true.)
2655 CALL dbcsr_add(dbcsr_p_work(ispin)%matrix, dbcsr_work, 1.0_dp, 1.0_dp)
2656 END DO
2657 CALL dbcsr_release(dbcsr_work)
2658 CALL dbcsr_deallocate_matrix_set(ker_tau_admm)
2659 ELSE
2660 CALL tddft_hfx_matrix(dbcsr_p_work, force_data%sum_YP_tau, qs_env, .false., .false.)
2661 END IF
2662 END IF
2663
2664 DO ispin = 1, nspins
2665 CALL dbcsr_add(force_data%sum_O_tau(ispin)%matrix, dbcsr_p_work(ispin)%matrix, 1.0_dp, 1.0_dp)
2666 END DO
2667
2668 CALL dbcsr_deallocate_matrix_set(dbcsr_p_work)
2669 CALL dbcsr_deallocate_matrix_set(work_admm)
2670
2671 CALL timestop(handle)
2672
2673 END SUBROUTINE prepare_for_response
2674
2675! **************************************************************************************************
2676!> \brief Calculate the force and virial due to the (P|Q) GPW integral derivatives
2677!> \param G_PQ ...
2678!> \param force ...
2679!> \param h_stress ...
2680!> \param use_virial ...
2681!> \param mp2_env ...
2682!> \param qs_env ...
2683! **************************************************************************************************
2684 SUBROUTINE get_2c_gpw_forces(G_PQ, force, h_stress, use_virial, mp2_env, qs_env)
2685
2686 TYPE(dbcsr_type), INTENT(INOUT) :: g_pq
2687 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2688 REAL(dp), DIMENSION(3, 3), INTENT(INOUT) :: h_stress
2689 LOGICAL, INTENT(IN) :: use_virial
2690 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
2691 TYPE(qs_environment_type), POINTER :: qs_env
2692
2693 CHARACTER(len=*), PARAMETER :: routinen = 'get_2c_gpw_forces'
2694
2695 INTEGER :: atom_a, color, handle, i, i_ri, i_xyz, iatom, igrid_level, ikind, ipgf, iset, j, &
2696 j_ri, jatom, lb_ri, n_ri, natom, ncoa, ncoms, nkind, nproc, nseta, o1, offset, pdims(2), &
2697 sgfa, ub_ri
2698 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, iproc_map, kind_of, &
2699 sizes_ri
2700 INTEGER, DIMENSION(:), POINTER :: col_dist, la_max, la_min, npgfa, nsgfa, &
2701 row_dist
2702 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, pgrid
2703 LOGICAL :: found, one_proc_group
2704 REAL(dp) :: cutoff_old, radius, relative_cutoff_old
2705 REAL(dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old, wf_vector
2706 REAL(dp), DIMENSION(3) :: force_a, force_b, ra
2707 REAL(dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
2708 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_tmp, i_ab, pab, pblock, sphi_a, zeta
2709 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2710 TYPE(cell_type), POINTER :: cell
2711 TYPE(dbcsr_distribution_type) :: dbcsr_dist
2712 TYPE(dbcsr_type) :: tmp_g_pq
2713 TYPE(dft_control_type), POINTER :: dft_control
2714 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
2715 DIMENSION(:), TARGET :: basis_set_ri_aux
2716 TYPE(gto_basis_set_type), POINTER :: basis_set_a
2717 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2718 TYPE(mp_para_env_type), POINTER :: para_env, para_env_ext
2719 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2720 POINTER :: sab_orb
2721 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2722 TYPE(pw_c1d_gs_type) :: dvg(3), pot_g, rho_g, rho_g_copy
2723 TYPE(pw_env_type), POINTER :: pw_env_ext
2724 TYPE(pw_poisson_type), POINTER :: poisson_env
2725 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2726 TYPE(pw_r3d_rs_type) :: psi_l, rho_r
2727 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2728 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
2729 TYPE(task_list_type), POINTER :: task_list_ext
2730
2731 NULLIFY (sab_orb, task_list_ext, particle_set, qs_kind_set, dft_control, pw_env_ext, auxbas_pw_pool, &
2732 poisson_env, atomic_kind_set, para_env, cell, rs_v, mos, basis_set_a)
2733
2734 CALL timeset(routinen, handle)
2735
2736 CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env, sab_orb=sab_orb, &
2737 natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, particle_set=particle_set, &
2738 mos=mos, cell=cell, atomic_kind_set=atomic_kind_set)
2739
2740 !The idea is to use GPW to compute the integrals and derivatives. Because the potential needs
2741 !to be calculated for each phi_j (column) of all AO pairs, and because that is expensive, we want
2742 !to minimize the amount of time we do that. Therefore, we work with a special distribution, where
2743 !each column of the resulting DBCSR matrix is mapped to a sub-communicator.
2744
2745 !Try to get the optimal pdims (we want a grid that is flat: many cols, few rows)
2746 IF (para_env%num_pe <= natom) THEN
2747 pdims(1) = 1
2748 pdims(2) = para_env%num_pe
2749 ELSE
2750 DO i = natom, 1, -1
2751 IF (modulo(para_env%num_pe, i) == 0) THEN
2752 pdims(1) = para_env%num_pe/i
2753 pdims(2) = i
2754 EXIT
2755 END IF
2756 END DO
2757 END IF
2758
2759 ALLOCATE (row_dist(natom), col_dist(natom))
2760 DO iatom = 1, natom
2761 row_dist(iatom) = modulo(iatom, pdims(1))
2762 END DO
2763 DO jatom = 1, natom
2764 col_dist(jatom) = modulo(jatom, pdims(2))
2765 END DO
2766
2767 ALLOCATE (pgrid(0:pdims(1) - 1, 0:pdims(2) - 1))
2768 nproc = 0
2769 DO i = 0, pdims(1) - 1
2770 DO j = 0, pdims(2) - 1
2771 pgrid(i, j) = nproc
2772 nproc = nproc + 1
2773 END DO
2774 END DO
2775
2776 CALL dbcsr_distribution_new(dbcsr_dist, group=para_env%get_handle(), pgrid=pgrid, row_dist=row_dist, col_dist=col_dist)
2777
2778 !The temporary DBCSR integrals and derivatives matrices in this flat distribution
2779 CALL dbcsr_create(tmp_g_pq, template=g_pq, matrix_type=dbcsr_type_no_symmetry, dist=dbcsr_dist)
2780 CALL dbcsr_complete_redistribute(g_pq, tmp_g_pq)
2781
2782 ALLOCATE (basis_set_ri_aux(nkind), sizes_ri(natom))
2783 CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
2784 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ri, basis=basis_set_ri_aux)
2785 n_ri = sum(sizes_ri)
2786
2787 one_proc_group = mp2_env%mp2_num_proc == 1
2788 ALLOCATE (para_env_ext)
2789 IF (one_proc_group) THEN
2790 !one subgroup per proc
2791 CALL para_env_ext%from_split(para_env, para_env%mepos)
2792 ELSE
2793 !Split the communicator accross the columns of the matrix
2794 ncoms = min(pdims(2), para_env%num_pe/mp2_env%mp2_num_proc)
2795 DO i = 0, pdims(1) - 1
2796 DO j = 0, pdims(2) - 1
2797 IF (pgrid(i, j) == para_env%mepos) color = modulo(j + 1, ncoms)
2798 END DO
2799 END DO
2800 CALL para_env_ext%from_split(para_env, color)
2801 END IF
2802
2803 !sab_orb and task_list_ext are essentially dummies
2804 CALL prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_ext, pw_env_ext, &
2805 auxbas_pw_pool, poisson_env, task_list_ext, rho_r, rho_g, pot_g, psi_l, sab_orb)
2806
2807 IF (use_virial) THEN
2808 CALL auxbas_pw_pool%create_pw(rho_g_copy)
2809 DO i_xyz = 1, 3
2810 CALL auxbas_pw_pool%create_pw(dvg(i_xyz))
2811 END DO
2812 END IF
2813
2814 ALLOCATE (wf_vector(n_ri))
2815
2816 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
2817
2818 ALLOCATE (iproc_map(natom))
2819
2820 !Loop over the atomic blocks
2821 DO jatom = 1, natom
2822
2823 !Only calculate if on the correct sub-communicator/proc
2824 IF (one_proc_group) THEN
2825 iproc_map = 0
2826 DO iatom = 1, natom
2827 IF (pgrid(row_dist(iatom), col_dist(jatom)) == para_env%mepos) iproc_map(iatom) = 1
2828 END DO
2829 IF (.NOT. any(iproc_map == 1)) cycle
2830 ELSE
2831 IF (.NOT. modulo(col_dist(jatom) + 1, ncoms) == color) cycle
2832 END IF
2833
2834 lb_ri = sum(sizes_ri(1:jatom - 1))
2835 ub_ri = lb_ri + sizes_ri(jatom)
2836 DO j_ri = lb_ri + 1, ub_ri
2837
2838 wf_vector = 0.0_dp
2839 wf_vector(j_ri) = 1.0_dp
2840
2841 CALL collocate_function(wf_vector, psi_l, rho_g, atomic_kind_set, qs_kind_set, cell, &
2842 particle_set, pw_env_ext, dft_control%qs_control%eps_rho_rspace, &
2843 basis_type="RI_AUX")
2844
2845 IF (use_virial) THEN
2846 CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, mp2_env%potential_parameter, dvg)
2847
2848 wf_vector = 0.0_dp
2849 DO iatom = 1, natom
2850 !only compute if i,j atom pair on correct proc
2851 IF (one_proc_group) THEN
2852 IF (.NOT. iproc_map(iatom) == 1) cycle
2853 END IF
2854
2855 CALL dbcsr_get_block_p(tmp_g_pq, iatom, jatom, pblock, found)
2856 IF (.NOT. found) cycle
2857
2858 i_ri = sum(sizes_ri(1:iatom - 1))
2859 wf_vector(i_ri + 1:i_ri + sizes_ri(iatom)) = pblock(:, j_ri - lb_ri)
2860 END DO
2861
2862 CALL pw_copy(rho_g, rho_g_copy)
2863 CALL collocate_function(wf_vector, psi_l, rho_g, atomic_kind_set, qs_kind_set, cell, &
2864 particle_set, pw_env_ext, dft_control%qs_control%eps_rho_rspace, &
2865 basis_type="RI_AUX")
2866
2867 CALL calc_potential_gpw(psi_l, rho_g, poisson_env, pot_g, mp2_env%potential_parameter, &
2868 no_transfer=.true.)
2869 CALL virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, &
2870 mp2_env%potential_parameter, para_env_ext)
2871 ELSE
2872 CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, mp2_env%potential_parameter)
2873 END IF
2874
2875 NULLIFY (rs_v)
2876 CALL pw_env_get(pw_env_ext, rs_grids=rs_v)
2877 CALL potential_pw2rs(rs_v, rho_r, pw_env_ext)
2878
2879 DO iatom = 1, natom
2880
2881 !only compute if i,j atom pair on correct proc
2882 IF (one_proc_group) THEN
2883 IF (.NOT. iproc_map(iatom) == 1) cycle
2884 END IF
2885
2886 force_a(:) = 0.0_dp
2887 force_b(:) = 0.0_dp
2888 IF (use_virial) THEN
2889 my_virial_a = 0.0_dp
2890 my_virial_b = 0.0_dp
2891 END IF
2892
2893 ikind = kind_of(iatom)
2894 atom_a = atom_of_kind(iatom)
2895
2896 basis_set_a => basis_set_ri_aux(ikind)%gto_basis_set
2897 first_sgfa => basis_set_a%first_sgf
2898 la_max => basis_set_a%lmax
2899 la_min => basis_set_a%lmin
2900 nseta = basis_set_a%nset
2901 nsgfa => basis_set_a%nsgf_set
2902 sphi_a => basis_set_a%sphi
2903 zeta => basis_set_a%zet
2904 npgfa => basis_set_a%npgf
2905
2906 ra(:) = pbc(particle_set(iatom)%r, cell)
2907
2908 CALL dbcsr_get_block_p(tmp_g_pq, iatom, jatom, pblock, found)
2909 IF (.NOT. found) cycle
2910
2911 offset = 0
2912 DO iset = 1, nseta
2913 ncoa = npgfa(iset)*ncoset(la_max(iset))
2914 sgfa = first_sgfa(1, iset)
2915
2916 ALLOCATE (h_tmp(ncoa, 1)); h_tmp = 0.0_dp
2917 ALLOCATE (i_ab(nsgfa(iset), 1)); i_ab = 0.0_dp
2918 ALLOCATE (pab(ncoa, 1)); pab = 0.0_dp
2919
2920 i_ab(1:nsgfa(iset), 1) = 2.0_dp*pblock(offset + 1:offset + nsgfa(iset), j_ri - lb_ri)
2921 CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
2922 i_ab(1, 1), nsgfa(iset), 0.0_dp, pab(1, 1), ncoa)
2923
2924 igrid_level = gaussian_gridlevel(pw_env_ext%gridlevel_info, minval(zeta(:, iset)))
2925
2926 ! The last three parameters are used to check whether a given function is within the own range.
2927 ! Here, it is always the case, so let's enforce it because mod(0, 1)==0
2928 IF (map_gaussian_here(rs_v(igrid_level), cell%h_inv, ra, 0, 1, 0)) THEN
2929 DO ipgf = 1, npgfa(iset)
2930 o1 = (ipgf - 1)*ncoset(la_max(iset))
2931 igrid_level = gaussian_gridlevel(pw_env_ext%gridlevel_info, zeta(ipgf, iset))
2932
2933 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
2934 lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
2935 zetp=zeta(ipgf, iset), &
2936 eps=dft_control%qs_control%eps_gvg_rspace, &
2937 prefactor=1.0_dp, cutoff=1.0_dp)
2938
2939 CALL integrate_pgf_product( &
2940 la_max=la_max(iset), zeta=zeta(ipgf, iset), la_min=la_min(iset), &
2941 lb_max=0, zetb=0.0_dp, lb_min=0, &
2942 ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
2943 rsgrid=rs_v(igrid_level), &
2944 hab=h_tmp, pab=pab, &
2945 o1=o1, &
2946 o2=0, &
2947 radius=radius, &
2948 calculate_forces=.true., &
2949 force_a=force_a, force_b=force_b, &
2950 use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
2951
2952 END DO
2953
2954 END IF
2955
2956 offset = offset + nsgfa(iset)
2957 DEALLOCATE (pab, h_tmp, i_ab)
2958 END DO !iset
2959
2960 force(ikind)%mp2_non_sep(:, atom_a) = force(ikind)%mp2_non_sep(:, atom_a) + force_a + force_b
2961 IF (use_virial) h_stress = h_stress + my_virial_a + my_virial_b
2962
2963 END DO !iatom
2964 END DO !j_RI
2965 END DO !jatom
2966
2967 IF (use_virial) THEN
2968 CALL auxbas_pw_pool%give_back_pw(rho_g_copy)
2969 DO i_xyz = 1, 3
2970 CALL auxbas_pw_pool%give_back_pw(dvg(i_xyz))
2971 END DO
2972 END IF
2973
2974 CALL cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_ext, pw_env_ext, &
2975 task_list_ext, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_l)
2976
2977 CALL dbcsr_release(tmp_g_pq)
2978 CALL dbcsr_distribution_release(dbcsr_dist)
2979 DEALLOCATE (col_dist, row_dist, pgrid)
2980
2981 CALL mp_para_env_release(para_env_ext)
2982
2983 CALL timestop(handle)
2984
2985 END SUBROUTINE get_2c_gpw_forces
2986
2987! **************************************************************************************************
2988!> \brief Calculate the forces due to the (P|Q) MME integral derivatives
2989!> \param G_PQ ...
2990!> \param force ...
2991!> \param mp2_env ...
2992!> \param qs_env ...
2993! **************************************************************************************************
2994 SUBROUTINE get_2c_mme_forces(G_PQ, force, mp2_env, qs_env)
2995
2996 TYPE(dbcsr_type), INTENT(INOUT) :: g_pq
2997 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2998 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
2999 TYPE(qs_environment_type), POINTER :: qs_env
3000
3001 CHARACTER(len=*), PARAMETER :: routinen = 'get_2c_mme_forces'
3002
3003 INTEGER :: atom_a, atom_b, g_count, handle, i_xyz, iatom, ikind, iset, jatom, jkind, jset, &
3004 natom, nkind, nseta, nsetb, offset_hab_a, offset_hab_b, r_count, sgfa, sgfb
3005 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
3006 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
3007 npgfb, nsgfa, nsgfb
3008 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
3009 LOGICAL :: found
3010 REAL(dp) :: new_force, pref
3011 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: hab
3012 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: hdab
3013 REAL(dp), DIMENSION(:, :), POINTER :: pblock
3014 REAL(kind=dp), DIMENSION(3) :: ra, rb
3015 REAL(kind=dp), DIMENSION(:, :), POINTER :: sphi_a, sphi_b, zeta, zetb
3016 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3017 TYPE(cell_type), POINTER :: cell
3018 TYPE(dbcsr_iterator_type) :: iter
3019 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
3020 DIMENSION(:), TARGET :: basis_set_ri_aux
3021 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
3022 TYPE(mp_para_env_type), POINTER :: para_env
3023 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3024 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3025
3026 NULLIFY (qs_kind_set, basis_set_a, basis_set_b, pblock, particle_set, &
3027 cell, la_max, la_min, lb_min, npgfa, lb_max, npgfb, nsgfa, &
3028 nsgfb, first_sgfa, first_sgfb, sphi_a, sphi_b, zeta, zetb, &
3029 atomic_kind_set, para_env)
3030
3031 CALL timeset(routinen, handle)
3032
3033 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind, particle_set=particle_set, &
3034 cell=cell, atomic_kind_set=atomic_kind_set, natom=natom, para_env=para_env)
3035
3036 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
3037
3038 ALLOCATE (basis_set_ri_aux(nkind))
3039 CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
3040
3041 g_count = 0; r_count = 0
3042
3043 CALL dbcsr_iterator_start(iter, g_pq)
3044 DO WHILE (dbcsr_iterator_blocks_left(iter))
3045
3046 CALL dbcsr_iterator_next_block(iter, row=iatom, column=jatom)
3047 CALL dbcsr_get_block_p(g_pq, iatom, jatom, pblock, found)
3048 IF (.NOT. found) cycle
3049 IF (iatom > jatom) cycle
3050 pref = 2.0_dp
3051 IF (iatom == jatom) pref = 1.0_dp
3052
3053 ikind = kind_of(iatom)
3054 jkind = kind_of(jatom)
3055
3056 atom_a = atom_of_kind(iatom)
3057 atom_b = atom_of_kind(jatom)
3058
3059 basis_set_a => basis_set_ri_aux(ikind)%gto_basis_set
3060 first_sgfa => basis_set_a%first_sgf
3061 la_max => basis_set_a%lmax
3062 la_min => basis_set_a%lmin
3063 nseta = basis_set_a%nset
3064 nsgfa => basis_set_a%nsgf_set
3065 sphi_a => basis_set_a%sphi
3066 zeta => basis_set_a%zet
3067 npgfa => basis_set_a%npgf
3068
3069 basis_set_b => basis_set_ri_aux(jkind)%gto_basis_set
3070 first_sgfb => basis_set_b%first_sgf
3071 lb_max => basis_set_b%lmax
3072 lb_min => basis_set_b%lmin
3073 nsetb = basis_set_b%nset
3074 nsgfb => basis_set_b%nsgf_set
3075 sphi_b => basis_set_b%sphi
3076 zetb => basis_set_b%zet
3077 npgfb => basis_set_b%npgf
3078
3079 ra(:) = pbc(particle_set(iatom)%r, cell)
3080 rb(:) = pbc(particle_set(jatom)%r, cell)
3081
3082 ALLOCATE (hab(basis_set_a%nsgf, basis_set_b%nsgf))
3083 ALLOCATE (hdab(3, basis_set_a%nsgf, basis_set_b%nsgf))
3084 hab(:, :) = 0.0_dp
3085 hdab(:, :, :) = 0.0_dp
3086
3087 offset_hab_a = 0
3088 DO iset = 1, nseta
3089 sgfa = first_sgfa(1, iset)
3090
3091 offset_hab_b = 0
3092 DO jset = 1, nsetb
3093 sgfb = first_sgfb(1, jset)
3094
3095 CALL integrate_set_2c(mp2_env%eri_mme_param%par, mp2_env%potential_parameter, la_min(iset), &
3096 la_max(iset), lb_min(jset), lb_max(jset), npgfa(iset), npgfb(jset), &
3097 zeta(:, iset), zetb(:, jset), ra, rb, hab, nsgfa(iset), nsgfb(jset), &
3098 offset_hab_a, offset_hab_b, 0, 0, sphi_a, sphi_b, sgfa, sgfb, &
3099 nsgfa(iset), nsgfb(jset), do_eri_mme, hdab=hdab, &
3100 g_count=g_count, r_count=r_count)
3101
3102 offset_hab_b = offset_hab_b + nsgfb(jset)
3103 END DO
3104 offset_hab_a = offset_hab_a + nsgfa(iset)
3105 END DO
3106
3107 DO i_xyz = 1, 3
3108 new_force = pref*sum(pblock(:, :)*hdab(i_xyz, :, :))
3109 force(ikind)%mp2_non_sep(i_xyz, atom_a) = force(ikind)%mp2_non_sep(i_xyz, atom_a) + new_force
3110 force(jkind)%mp2_non_sep(i_xyz, atom_b) = force(jkind)%mp2_non_sep(i_xyz, atom_b) - new_force
3111 END DO
3112
3113 DEALLOCATE (hab, hdab)
3114 END DO
3115 CALL dbcsr_iterator_stop(iter)
3116
3117 CALL cp_eri_mme_update_local_counts(mp2_env%eri_mme_param, para_env, g_count_2c=g_count, r_count_2c=r_count)
3118
3119 CALL timestop(handle)
3120
3121 END SUBROUTINE get_2c_mme_forces
3122
3123! **************************************************************************************************
3124!> \brief This routines gather all the force updates due to the response density and the trace with F
3125!> Also update the forces due to the SCF density for XC and exact exchange
3126!> \param p_env the p_env coming from the response calculation
3127!> \param matrix_hz the matrix going into the RHS of the response equation
3128!> \param matrix_p_F the density matrix with which we evaluate Trace[P*F]
3129!> \param matrix_p_F_admm ...
3130!> \param qs_env ...
3131!> \note very much inspired from the response_force routine in response_solver.F, especially for virial
3132! **************************************************************************************************
3133 SUBROUTINE update_im_time_forces(p_env, matrix_hz, matrix_p_F, matrix_p_F_admm, qs_env)
3134
3135 TYPE(qs_p_env_type), POINTER :: p_env
3136 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz, matrix_p_f, matrix_p_f_admm
3137 TYPE(qs_environment_type), POINTER :: qs_env
3138
3139 CHARACTER(len=*), PARAMETER :: routinen = 'update_im_time_forces'
3140
3141 INTEGER :: handle, i, idens, ispin, n_rep_hf, nao, &
3142 nao_aux, nder, nimages, nocc, nspins
3143 LOGICAL :: do_exx, do_hfx, do_tau, do_tau_admm, &
3144 use_virial
3145 REAL(dp) :: dummy_real1, dummy_real2, ehartree, exc, &
3146 focc
3147 REAL(dp), DIMENSION(3, 3) :: h_stress, pv_loc
3148 TYPE(admm_type), POINTER :: admm_env
3149 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3150 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: current_density, current_density_admm, &
3151 current_mat_h, matrix_p_mp2, matrix_p_mp2_admm, matrix_s, matrix_s_aux_fit, matrix_w, &
3152 rho_ao, rho_ao_aux, scrm, scrm_admm
3153 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: dbcsr_work_h, dbcsr_work_p, mpa2
3154 TYPE(dbcsr_type) :: dbcsr_work
3155 TYPE(dft_control_type), POINTER :: dft_control
3156 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
3157 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3158 TYPE(mp_para_env_type), POINTER :: para_env
3159 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3160 POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
3161 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3162 TYPE(pw_c1d_gs_type) :: rho_tot_gspace, rhoz_tot_gspace, &
3163 zv_hartree_gspace
3164 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rhoz_g
3165 TYPE(pw_env_type), POINTER :: pw_env
3166 TYPE(pw_poisson_type), POINTER :: poisson_env
3167 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
3168 TYPE(pw_r3d_rs_type) :: vh_rspace, vhxc_rspace, zv_hartree_rspace
3169 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rhoz_r, tauz_r, v_xc, v_xc_tau, &
3170 vadmm_rspace, vtau_rspace, vxc_rspace
3171 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3172 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3173 TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
3174 TYPE(section_vals_type), POINTER :: hfx_section, xc_section
3175 TYPE(task_list_type), POINTER :: task_list_aux_fit
3176 TYPE(virial_type), POINTER :: virial
3177
3178 NULLIFY (scrm, rho, dft_control, matrix_p_mp2, matrix_s, &
3179 matrix_p_mp2_admm, admm_env, sab_orb, dbcsr_work_p, &
3180 dbcsr_work_h, sac_ae, sac_ppl, sap_ppnl, force, virial, &
3181 qs_kind_set, atomic_kind_set, particle_set, pw_env, poisson_env, &
3182 auxbas_pw_pool, task_list_aux_fit, matrix_s_aux_fit, scrm_admm, &
3183 rho_aux_fit, rho_ao_aux, x_data, hfx_section, xc_section, &
3184 para_env, rhoz_g, rhoz_r, tauz_r, v_xc, v_xc_tau, vxc_rspace, &
3185 vtau_rspace, vadmm_rspace, rho_ao, matrix_w)
3186
3187 CALL timeset(routinen, handle)
3188
3189 CALL get_qs_env(qs_env, rho=rho, dft_control=dft_control, matrix_s=matrix_s, admm_env=admm_env, &
3190 sab_orb=sab_orb, sac_ae=sac_ae, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl, force=force, &
3191 virial=virial, particle_set=particle_set, qs_kind_set=qs_kind_set, &
3192 atomic_kind_set=atomic_kind_set, x_data=x_data, para_env=para_env)
3193 nspins = dft_control%nspins
3194
3195 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
3196 IF (use_virial) virial%pv_calculate = .true.
3197
3198 !Whether we replace the force/energy of SCF XC with HF in RPA
3199 do_exx = .false.
3200 IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
3201 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
3202 CALL section_vals_get(hfx_section, explicit=do_exx)
3203 END IF
3204
3205 !Get the mp2 density matrix which is p_env%p1 + matrix_p_F
3206 CALL get_qs_env(qs_env, matrix_p_mp2=matrix_p_mp2, matrix_p_mp2_admm=matrix_p_mp2_admm)
3207
3208 !The kinetic term (only response density)
3209 NULLIFY (scrm)
3210 mpa2(1:nspins, 1:1) => matrix_p_mp2(1:nspins)
3211 CALL kinetic_energy_matrix(qs_env, matrix_t=scrm, matrix_p=mpa2, &
3212 matrix_name="KINETIC ENERGY MATRIX", &
3213 basis_type="ORB", &
3214 sab_orb=sab_orb, calculate_forces=.true.)
3215 CALL dbcsr_deallocate_matrix_set(scrm)
3216
3217 !The pseudo-potential terms (only reponse density)
3218 CALL dbcsr_allocate_matrix_set(scrm, nspins)
3219 DO ispin = 1, nspins
3220 ALLOCATE (scrm(ispin)%matrix)
3221 CALL dbcsr_create(scrm(ispin)%matrix, template=matrix_s(1)%matrix)
3222 CALL dbcsr_copy(scrm(ispin)%matrix, matrix_s(1)%matrix)
3223 CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
3224 END DO
3225
3226 nder = 1
3227 nimages = 1
3228 ALLOCATE (dbcsr_work_p(nspins, 1), dbcsr_work_h(nspins, 1))
3229 DO ispin = 1, nspins
3230 dbcsr_work_p(ispin, 1)%matrix => matrix_p_mp2(ispin)%matrix
3231 dbcsr_work_h(ispin, 1)%matrix => scrm(ispin)%matrix
3232 END DO
3233
3234 CALL core_matrices(qs_env, dbcsr_work_h, dbcsr_work_p, .true., nder)
3235
3236 DEALLOCATE (dbcsr_work_p, dbcsr_work_h)
3237
3238 IF (use_virial) THEN
3239 h_stress = 0.0_dp
3240 virial%pv_xc = 0.0_dp
3241 NULLIFY (vxc_rspace, vtau_rspace, vadmm_rspace)
3242 CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, &
3243 dummy_real1, dummy_real2, h_stress)
3244 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
3245 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
3246 IF (.NOT. do_exx) THEN
3247 !if RPA EXX, then do not consider XC virial (replaced by RPA%HF virial)
3248 virial%pv_exc = virial%pv_exc - virial%pv_xc
3249 virial%pv_virial = virial%pv_virial - virial%pv_xc
3250 END IF
3251 ELSE
3252 CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, dummy_real1, dummy_real2)
3253 END IF
3254 do_tau = ASSOCIATED(vtau_rspace)
3255
3256 !Core forces from the SCF
3257 CALL integrate_v_core_rspace(vh_rspace, qs_env)
3258
3259 !The Hartree-xc potential term, P*dVHxc (mp2 + SCF density x deriv of the SCF potential)
3260 !Get the total density
3261 CALL qs_rho_get(rho, rho_ao=rho_ao)
3262 DO ispin = 1, nspins
3263 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
3264 END DO
3265
3266 CALL get_qs_env(qs_env, pw_env=pw_env)
3267 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
3268 poisson_env=poisson_env)
3269 CALL auxbas_pw_pool%create_pw(vhxc_rspace)
3270
3271 IF (use_virial) pv_loc = virial%pv_virial
3272
3273 IF (do_exx) THEN
3274 !Only want response XC contribution, but SCF+response Hartree contribution
3275 DO ispin = 1, nspins
3276 !Hartree
3277 CALL pw_transfer(vh_rspace, vhxc_rspace)
3278 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
3279 hmat=scrm(ispin), pmat=rho_ao(ispin), &
3280 qs_env=qs_env, calculate_forces=.true.)
3281 !XC
3282 CALL pw_transfer(vxc_rspace(ispin), vhxc_rspace)
3283 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
3284 hmat=scrm(ispin), pmat=matrix_p_mp2(ispin), &
3285 qs_env=qs_env, calculate_forces=.true.)
3286 IF (do_tau) THEN
3287 CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
3288 hmat=scrm(ispin), pmat=matrix_p_mp2(ispin), &
3289 qs_env=qs_env, calculate_forces=.true., compute_tau=.true.)
3290 END IF
3291 END DO
3292 ELSE
3293 DO ispin = 1, nspins
3294 CALL pw_transfer(vh_rspace, vhxc_rspace)
3295 CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
3296 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
3297 hmat=scrm(ispin), pmat=rho_ao(ispin), &
3298 qs_env=qs_env, calculate_forces=.true.)
3299 IF (do_tau) THEN
3300 CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
3301 hmat=scrm(ispin), pmat=rho_ao(ispin), &
3302 qs_env=qs_env, calculate_forces=.true., compute_tau=.true.)
3303 END IF
3304 END DO
3305 END IF
3306 CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
3307
3308 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3309
3310 !The admm projection contribution (mp2 + SCF densities). If EXX, then only mp2 density
3311 IF (dft_control%do_admm) THEN
3312 CALL get_admm_env(admm_env, task_list_aux_fit=task_list_aux_fit, rho_aux_fit=rho_aux_fit, &
3313 matrix_s_aux_fit=matrix_s_aux_fit)
3314 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux)
3315 CALL dbcsr_allocate_matrix_set(scrm_admm, nspins)
3316 DO ispin = 1, nspins
3317 ALLOCATE (scrm_admm(ispin)%matrix)
3318 CALL dbcsr_create(scrm_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
3319 CALL dbcsr_copy(scrm_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
3320 CALL dbcsr_set(scrm_admm(ispin)%matrix, 0.0_dp)
3321 END DO
3322
3323 IF (use_virial) pv_loc = virial%pv_virial
3324 IF (.NOT. qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
3325 DO ispin = 1, nspins
3326 IF (do_exx) THEN
3327 CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), &
3328 hmat=scrm_admm(ispin), pmat=matrix_p_mp2_admm(ispin), &
3329 qs_env=qs_env, calculate_forces=.true., &
3330 basis_type="AUX_FIT", task_list_external=task_list_aux_fit)
3331 ELSE
3332 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
3333 CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), &
3334 hmat=scrm_admm(ispin), pmat=rho_ao_aux(ispin), &
3335 qs_env=qs_env, calculate_forces=.true., &
3336 basis_type="AUX_FIT", task_list_external=task_list_aux_fit)
3337 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
3338 END IF
3339 END DO
3340 END IF
3341 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3342
3343 CALL tddft_hfx_matrix(scrm_admm, rho_ao_aux, qs_env, .false., .false.)
3344
3345 IF (do_exx) THEN
3346 CALL admm_projection_derivative(qs_env, scrm_admm, matrix_p_mp2)
3347 ELSE
3348 CALL admm_projection_derivative(qs_env, scrm_admm, rho_ao)
3349 END IF
3350 END IF
3351
3352 !The exact-exchange term (mp2 + SCF densities)
3353 xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
3354 hfx_section => section_vals_get_subs_vals(xc_section, "HF")
3355 CALL section_vals_get(hfx_section, explicit=do_hfx)
3356
3357 IF (do_hfx) THEN
3358 CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
3359 cpassert(n_rep_hf == 1)
3360 IF (use_virial) virial%pv_fock_4c = 0.0_dp
3361
3362 !In case of EXX, only want to response HFX forces, as the SCF will change according to RI_RPA%HF
3363 IF (do_exx) THEN
3364 IF (dft_control%do_admm) THEN
3365 CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
3366 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux, rho_ao_kp=dbcsr_work_p)
3367 IF (x_data(1, 1)%do_hfx_ri) THEN
3368
3369 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3370 x_data(1, 1)%general_parameter%fraction, &
3371 rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2_admm, &
3372 use_virial=use_virial, resp_only=.true.)
3373 ELSE
3374 CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2_admm, hfx_section, para_env, &
3375 1, use_virial, resp_only=.true.)
3376 END IF
3377 ELSE
3378 DO ispin = 1, nspins
3379 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
3380 END DO
3381 CALL qs_rho_get(rho, 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, &
3387 use_virial=use_virial, resp_only=.true.)
3388 ELSE
3389 CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2, hfx_section, para_env, &
3390 1, use_virial, resp_only=.true.)
3391 END IF
3392 DO ispin = 1, nspins
3393 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
3394 END DO
3395 END IF !admm
3396
3397 ELSE !No Exx
3398 IF (dft_control%do_admm) THEN
3399 CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
3400 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux, rho_ao_kp=dbcsr_work_p)
3401 DO ispin = 1, nspins
3402 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
3403 END DO
3404 IF (x_data(1, 1)%do_hfx_ri) THEN
3405
3406 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3407 x_data(1, 1)%general_parameter%fraction, &
3408 rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2_admm, &
3409 use_virial=use_virial, resp_only=.false.)
3410 ELSE
3411 CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2_admm, hfx_section, para_env, &
3412 1, use_virial, resp_only=.false.)
3413 END IF
3414 DO ispin = 1, nspins
3415 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
3416 END DO
3417 ELSE
3418 CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
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, &
3424 use_virial=use_virial, resp_only=.false.)
3425 ELSE
3426 CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2, hfx_section, para_env, &
3427 1, use_virial, resp_only=.false.)
3428 END IF
3429 END IF
3430 END IF !do_exx
3431
3432 IF (use_virial) THEN
3433 virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
3434 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
3435 END IF
3436 END IF
3437
3438 !retrieve the SCF density
3439 CALL qs_rho_get(rho, rho_ao=rho_ao)
3440 DO ispin = 1, nspins
3441 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
3442 END DO
3443
3444 !From here, we need to do everything twice. Once for the response density, and once for the
3445 !density that is used for the trace Tr[P*F]. The reason is that the former is needed for the
3446 !eventual overlap contribution from matrix_wz
3447 !Only with the mp2 density
3448
3449 ALLOCATE (current_density(nspins), current_mat_h(nspins), current_density_admm(nspins))
3450 DO idens = 1, 2
3451 DO ispin = 1, nspins
3452 IF (idens == 1) THEN
3453 current_density(ispin)%matrix => matrix_p_f(ispin)%matrix
3454 current_mat_h(ispin)%matrix => scrm(ispin)%matrix
3455 IF (dft_control%do_admm) current_density_admm(ispin)%matrix => matrix_p_f_admm(ispin)%matrix
3456 ELSE
3457 current_density(ispin)%matrix => p_env%p1(ispin)%matrix
3458 current_mat_h(ispin)%matrix => matrix_hz(ispin)%matrix
3459 IF (dft_control%do_admm) current_density_admm(ispin)%matrix => p_env%p1_admm(ispin)%matrix
3460 END IF
3461 END DO
3462
3463 !The core-denstiy derivative
3464 ALLOCATE (rhoz_r(nspins), rhoz_g(nspins))
3465 DO ispin = 1, nspins
3466 CALL auxbas_pw_pool%create_pw(rhoz_r(ispin))
3467 CALL auxbas_pw_pool%create_pw(rhoz_g(ispin))
3468 END DO
3469 CALL auxbas_pw_pool%create_pw(rhoz_tot_gspace)
3470 CALL auxbas_pw_pool%create_pw(zv_hartree_rspace)
3471 CALL auxbas_pw_pool%create_pw(zv_hartree_gspace)
3472
3473 CALL pw_zero(rhoz_tot_gspace)
3474 DO ispin = 1, nspins
3475 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density(ispin)%matrix, &
3476 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin))
3477 CALL pw_axpy(rhoz_g(ispin), rhoz_tot_gspace)
3478 END DO
3479
3480 IF (use_virial) THEN
3481
3482 CALL get_qs_env(qs_env, rho=rho)
3483 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
3484
3485 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
3486
3487 h_stress(:, :) = 0.0_dp
3488 CALL pw_poisson_solve(poisson_env, &
3489 density=rhoz_tot_gspace, &
3490 ehartree=ehartree, &
3491 vhartree=zv_hartree_gspace, &
3492 h_stress=h_stress, &
3493 aux_density=rho_tot_gspace)
3494
3495 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
3496
3497 !Green contribution
3498 virial%pv_ehartree = virial%pv_ehartree + 2.0_dp*h_stress/real(para_env%num_pe, dp)
3499 virial%pv_virial = virial%pv_virial + 2.0_dp*h_stress/real(para_env%num_pe, dp)
3500
3501 ELSE
3502 CALL pw_poisson_solve(poisson_env, rhoz_tot_gspace, ehartree, &
3503 zv_hartree_gspace)
3504 END IF
3505
3506 CALL pw_transfer(zv_hartree_gspace, zv_hartree_rspace)
3507 CALL pw_scale(zv_hartree_rspace, zv_hartree_rspace%pw_grid%dvol)
3508 CALL integrate_v_core_rspace(zv_hartree_rspace, qs_env)
3509
3510 IF (do_tau) THEN
3511 block
3512 TYPE(pw_c1d_gs_type) :: tauz_g
3513 CALL auxbas_pw_pool%create_pw(tauz_g)
3514 ALLOCATE (tauz_r(nspins))
3515 DO ispin = 1, nspins
3516 CALL auxbas_pw_pool%create_pw(tauz_r(ispin))
3517
3518 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density(ispin)%matrix, &
3519 rho=tauz_r(ispin), rho_gspace=tauz_g, compute_tau=.true.)
3520 END DO
3521 CALL auxbas_pw_pool%give_back_pw(tauz_g)
3522 END block
3523 END IF
3524
3525 !Volume contribution to the virial
3526 IF (use_virial) THEN
3527 !Volume contribution
3528 exc = 0.0_dp
3529 DO ispin = 1, nspins
3530 exc = exc + pw_integral_ab(rhoz_r(ispin), vxc_rspace(ispin))/ &
3531 vxc_rspace(ispin)%pw_grid%dvol
3532 END DO
3533 IF (ASSOCIATED(vtau_rspace)) THEN
3534 DO ispin = 1, nspins
3535 exc = exc + pw_integral_ab(tauz_r(ispin), vtau_rspace(ispin))/ &
3536 vtau_rspace(ispin)%pw_grid%dvol
3537 END DO
3538 END IF
3539 DO i = 1, 3
3540 virial%pv_ehartree(i, i) = virial%pv_ehartree(i, i) - 4.0_dp*ehartree/real(para_env%num_pe, dp)
3541 virial%pv_exc(i, i) = virial%pv_exc(i, i) - exc/real(para_env%num_pe, dp)
3542 virial%pv_virial(i, i) = virial%pv_virial(i, i) - 4.0_dp*ehartree/real(para_env%num_pe, dp) &
3543 - exc/real(para_env%num_pe, dp)
3544 END DO
3545 END IF
3546
3547 !The xc-kernel term.
3548 IF (dft_control%do_admm) THEN
3549 CALL get_qs_env(qs_env, admm_env=admm_env)
3550 xc_section => admm_env%xc_section_primary
3551 ELSE
3552 xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
3553 END IF
3554
3555 IF (use_virial) virial%pv_xc = 0.0_dp
3556
3557 CALL create_kernel(qs_env, &
3558 vxc=v_xc, &
3559 vxc_tau=v_xc_tau, &
3560 rho=rho, &
3561 rho1_r=rhoz_r, &
3562 rho1_g=rhoz_g, &
3563 tau1_r=tauz_r, &
3564 xc_section=xc_section, &
3565 compute_virial=use_virial, &
3566 virial_xc=virial%pv_xc)
3567
3568 IF (use_virial) THEN
3569 virial%pv_exc = virial%pv_exc + virial%pv_xc
3570 virial%pv_virial = virial%pv_virial + virial%pv_xc
3571
3572 pv_loc = virial%pv_virial
3573 END IF
3574
3575 CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3576 DO ispin = 1, nspins
3577 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
3578 CALL pw_axpy(zv_hartree_rspace, v_xc(ispin))
3579 CALL integrate_v_rspace(qs_env=qs_env, &
3580 v_rspace=v_xc(ispin), &
3581 hmat=current_mat_h(ispin), &
3582 pmat=dbcsr_work_p(ispin, 1), &
3583 calculate_forces=.true.)
3584 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
3585 END DO
3586 CALL auxbas_pw_pool%give_back_pw(rhoz_tot_gspace)
3587 CALL auxbas_pw_pool%give_back_pw(zv_hartree_rspace)
3588 CALL auxbas_pw_pool%give_back_pw(zv_hartree_gspace)
3589 DEALLOCATE (v_xc)
3590
3591 IF (do_tau) THEN
3592 DO ispin = 1, nspins
3593 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
3594 CALL integrate_v_rspace(qs_env=qs_env, &
3595 v_rspace=v_xc_tau(ispin), &
3596 hmat=current_mat_h(ispin), &
3597 pmat=dbcsr_work_p(ispin, 1), &
3598 compute_tau=.true., &
3599 calculate_forces=.true.)
3600 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
3601 END DO
3602 DEALLOCATE (v_xc_tau)
3603 END IF
3604
3605 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3606
3607 IF (do_hfx) THEN
3608 IF (dft_control%do_admm) THEN
3609 DO ispin = 1, nspins
3610 CALL dbcsr_set(scrm_admm(ispin)%matrix, 0.0_dp)
3611 END DO
3612 CALL qs_rho_get(rho_aux_fit, tau_r_valid=do_tau_admm)
3613
3614 IF (.NOT. admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
3615 CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
3616 DO ispin = 1, nspins
3617 CALL pw_zero(rhoz_r(ispin))
3618 CALL pw_zero(rhoz_g(ispin))
3619 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density_admm(ispin)%matrix, &
3620 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin), &
3621 basis_type="AUX_FIT", task_list_external=task_list_aux_fit)
3622 END DO
3623
3624 IF (do_tau_admm) THEN
3625 block
3626 TYPE(pw_c1d_gs_type) :: tauz_g
3627 CALL auxbas_pw_pool%create_pw(tauz_g)
3628 DO ispin = 1, nspins
3629 CALL pw_zero(tauz_r(ispin))
3630 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density(ispin)%matrix, &
3631 rho=tauz_r(ispin), rho_gspace=tauz_g, &
3632 basis_type="AUX_FIT", task_list_external=task_list_aux_fit, &
3633 compute_tau=.true.)
3634 END DO
3635 CALL auxbas_pw_pool%give_back_pw(tauz_g)
3636 END block
3637 END IF
3638
3639 !Volume contribution to the virial
3640 IF (use_virial) THEN
3641 exc = 0.0_dp
3642 DO ispin = 1, nspins
3643 exc = exc + pw_integral_ab(rhoz_r(ispin), vadmm_rspace(ispin))/ &
3644 vadmm_rspace(ispin)%pw_grid%dvol
3645 END DO
3646 DO i = 1, 3
3647 virial%pv_exc(i, i) = virial%pv_exc(i, i) - exc/real(para_env%num_pe, dp)
3648 virial%pv_virial(i, i) = virial%pv_virial(i, i) - exc/real(para_env%num_pe, dp)
3649 END DO
3650
3651 virial%pv_xc = 0.0_dp
3652 END IF
3653
3654 xc_section => admm_env%xc_section_aux
3655 CALL create_kernel(qs_env, v_xc, v_xc_tau, rho_aux_fit, rhoz_r, rhoz_g, tauz_r, xc_section, &
3656 compute_virial=use_virial, virial_xc=virial%pv_xc)
3657
3658 IF (use_virial) THEN
3659 virial%pv_exc = virial%pv_exc + virial%pv_xc
3660 virial%pv_virial = virial%pv_virial + virial%pv_xc
3661
3662 pv_loc = virial%pv_virial
3663 END IF
3664
3665 CALL qs_rho_get(rho_aux_fit, rho_ao_kp=dbcsr_work_p)
3666 DO ispin = 1, nspins
3667 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
3668 CALL integrate_v_rspace(qs_env=qs_env, &
3669 v_rspace=v_xc(ispin), &
3670 hmat=scrm_admm(ispin), &
3671 pmat=dbcsr_work_p(ispin, 1), &
3672 calculate_forces=.true., &
3673 basis_type="AUX_FIT", &
3674 task_list_external=task_list_aux_fit)
3675 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
3676 END DO
3677 DEALLOCATE (v_xc)
3678
3679 IF (do_tau_admm) THEN
3680 DO ispin = 1, nspins
3681 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
3682 CALL integrate_v_rspace(qs_env=qs_env, &
3683 v_rspace=v_xc_tau(ispin), &
3684 hmat=scrm_admm(ispin), &
3685 pmat=dbcsr_work_p(ispin, 1), &
3686 calculate_forces=.true., &
3687 basis_type="AUX_FIT", &
3688 task_list_external=task_list_aux_fit, &
3689 compute_tau=.true.)
3690 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
3691 END DO
3692 DEALLOCATE (v_xc_tau)
3693 END IF
3694
3695 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3696 END IF
3697
3698 CALL tddft_hfx_matrix(scrm_admm, current_density_admm, qs_env, .false., .false.)
3699
3700 CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3701 CALL admm_projection_derivative(qs_env, scrm_admm, dbcsr_work_p(:, 1))
3702
3703 !If response density, need to get matrix_hz contribution
3704 CALL dbcsr_create(dbcsr_work, template=matrix_s(1)%matrix)
3705 IF (idens == 2) THEN
3706 nao = admm_env%nao_orb
3707 nao_aux = admm_env%nao_aux_fit
3708 DO ispin = 1, nspins
3709 CALL dbcsr_copy(dbcsr_work, matrix_hz(ispin)%matrix)
3710 CALL dbcsr_set(dbcsr_work, 0.0_dp)
3711
3712 CALL cp_dbcsr_sm_fm_multiply(scrm_admm(ispin)%matrix, admm_env%A, &
3713 admm_env%work_aux_orb, nao)
3714 CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
3715 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
3716 admm_env%work_orb_orb)
3717 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_work, keep_sparsity=.true.)
3718 CALL dbcsr_add(matrix_hz(ispin)%matrix, dbcsr_work, 1.0_dp, 1.0_dp)
3719 END DO
3720 END IF
3721
3722 CALL dbcsr_release(dbcsr_work)
3723 ELSE !no admm
3724
3725 !Need the contribution to matrix_hz as well
3726 IF (idens == 2) THEN
3727 CALL tddft_hfx_matrix(matrix_hz, current_density, qs_env, .false., .false.)
3728 END IF
3729 END IF !admm
3730 END IF !do_hfx
3731
3732 DO ispin = 1, nspins
3733 CALL auxbas_pw_pool%give_back_pw(rhoz_r(ispin))
3734 CALL auxbas_pw_pool%give_back_pw(rhoz_g(ispin))
3735 END DO
3736 DEALLOCATE (rhoz_r, rhoz_g)
3737
3738 IF (do_tau) THEN
3739 DO ispin = 1, nspins
3740 CALL auxbas_pw_pool%give_back_pw(tauz_r(ispin))
3741 END DO
3742 DEALLOCATE (tauz_r)
3743 END IF
3744 END DO !idens
3745 CALL dbcsr_deallocate_matrix_set(scrm_admm)
3746
3747 DEALLOCATE (current_density, current_mat_h, current_density_admm)
3748 CALL dbcsr_deallocate_matrix_set(scrm)
3749
3750 !The energy weighted and overlap term. ONLY with the response density
3751 focc = 2.0_dp
3752 IF (nspins == 2) focc = 1.0_dp
3753 CALL get_qs_env(qs_env, mos=mos)
3754 DO ispin = 1, nspins
3755 CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
3756 CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
3757 p_env%w1(ispin)%matrix, focc, nocc)
3758 END DO
3759 IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, p_env%w1(2)%matrix, 1.0_dp, 1.0_dp)
3760
3761 !Add to it the SCF W matrix, except if EXX (because taken care of by HF response)
3762 IF (.NOT. do_exx) THEN
3763 CALL compute_matrix_w(qs_env, calc_forces=.true.)
3764 CALL get_qs_env(qs_env, matrix_w=matrix_w)
3765 CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(1)%matrix, 1.0_dp, 1.0_dp)
3766 IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(2)%matrix, 1.0_dp, 1.0_dp)
3767 END IF
3768
3769 NULLIFY (scrm)
3770 CALL build_overlap_matrix(qs_env%ks_env, matrix_s=scrm, &
3771 matrix_name="OVERLAP MATRIX", &
3772 basis_type_a="ORB", basis_type_b="ORB", &
3773 sab_nl=sab_orb, calculate_forces=.true., &
3774 matrix_p=p_env%w1(1)%matrix)
3775
3776 IF (.NOT. do_exx) THEN
3777 CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(1)%matrix, 1.0_dp, -1.0_dp)
3778 IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(2)%matrix, 1.0_dp, -1.0_dp)
3779 DO ispin = 1, nspins
3780 CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
3781 END DO
3782 END IF
3783
3784 IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, p_env%w1(2)%matrix, 1.0_dp, -1.0_dp)
3785 CALL dbcsr_deallocate_matrix_set(scrm)
3786
3787 IF (use_virial) virial%pv_calculate = .false.
3788
3789 !clean-up
3790 CALL auxbas_pw_pool%give_back_pw(vh_rspace)
3791
3792 DO ispin = 1, nspins
3793 CALL auxbas_pw_pool%give_back_pw(vxc_rspace(ispin))
3794 IF (ASSOCIATED(vtau_rspace)) THEN
3795 CALL auxbas_pw_pool%give_back_pw(vtau_rspace(ispin))
3796 END IF
3797 IF (ASSOCIATED(vadmm_rspace)) THEN
3798 CALL auxbas_pw_pool%give_back_pw(vadmm_rspace(ispin))
3799 END IF
3800 END DO
3801 DEALLOCATE (vxc_rspace)
3802 IF (ASSOCIATED(vtau_rspace)) DEALLOCATE (vtau_rspace)
3803 IF (ASSOCIATED(vadmm_rspace)) DEALLOCATE (vadmm_rspace)
3804
3805 CALL timestop(handle)
3806
3807 END SUBROUTINE update_im_time_forces
3808
3809! **************************************************************************************************
3810!> \brief Iteratively builds the matrix Y = sum_k Y_k until convergence, where
3811!> Y_k = 1/k*2^n (A/2^n) Y_k-1 + 1/k!*2^n * PR(n) * (A/2^n)^(k-1)
3812!> n is chosen such that the norm of A is < 1 (and e^A converges fast)
3813!> PR(n) = e^(A/2^n)*PR(n-1) + PR(n-1)*e^(A/2^n), PR(0) = P*R^T
3814!> \param Y ...
3815!> \param A ...
3816!> \param P ...
3817!> \param R ...
3818!> \param filter_eps ...
3819! **************************************************************************************************
3820 SUBROUTINE build_y_matrix(Y, A, P, R, filter_eps)
3821
3822 TYPE(dbcsr_type), INTENT(OUT) :: y
3823 TYPE(dbcsr_type), INTENT(INOUT) :: a, p, r
3824 REAL(dp), INTENT(IN) :: filter_eps
3825
3826 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_Y_matrix'
3827
3828 INTEGER :: handle, k, n
3829 REAL(dp) :: norm_scalar, threshold
3830 TYPE(dbcsr_type) :: a2n, exp_a2n, prn, work, work2, yk
3831
3832 CALL timeset(routinen, handle)
3833
3834 threshold = 1.0e-16_dp
3835
3836 !Find n such that norm(A) < 1 and we insure convergence of the exponential
3837 norm_scalar = dbcsr_frobenius_norm(a)
3838
3839 !checked: result invariant with value of n
3840 n = 1
3841 DO
3842 IF ((norm_scalar/2.0_dp**n) < 1.0_dp) EXIT
3843 n = n + 1
3844 END DO
3845
3846 !Calculate PR(n) recursively
3847 CALL dbcsr_create(prn, template=a, matrix_type=dbcsr_type_no_symmetry)
3848 CALL dbcsr_create(work, template=a, matrix_type=dbcsr_type_no_symmetry)
3849 CALL dbcsr_multiply('N', 'N', 1.0_dp, p, r, 0.0_dp, work, filter_eps=filter_eps)
3850 CALL dbcsr_create(exp_a2n, template=a, matrix_type=dbcsr_type_no_symmetry)
3851
3852 DO k = 1, n
3853 CALL matrix_exponential(exp_a2n, a, 1.0_dp, 0.5_dp**k, threshold)
3854 CALL dbcsr_multiply('N', 'N', 1.0_dp, exp_a2n, work, 0.0_dp, prn, filter_eps=filter_eps)
3855 CALL dbcsr_multiply('N', 'N', 1.0_dp, work, exp_a2n, 1.0_dp, prn, filter_eps=filter_eps)
3856 CALL dbcsr_copy(work, prn)
3857 END DO
3858 CALL dbcsr_release(exp_a2n)
3859
3860 !Calculate Y iteratively, until convergence
3861 CALL dbcsr_create(a2n, template=a, matrix_type=dbcsr_type_no_symmetry)
3862 CALL dbcsr_copy(a2n, a)
3863 CALL dbcsr_scale(a2n, 0.5_dp**n)
3864 CALL dbcsr_create(y, template=a, matrix_type=dbcsr_type_no_symmetry)
3865 CALL dbcsr_create(yk, template=a, matrix_type=dbcsr_type_no_symmetry)
3866 CALL dbcsr_create(work2, template=a, matrix_type=dbcsr_type_no_symmetry)
3867
3868 !k=1
3869 CALL dbcsr_scale(prn, 0.5_dp**n)
3870 CALL dbcsr_copy(work, prn)
3871 CALL dbcsr_copy(work2, prn)
3872 CALL dbcsr_add(y, prn, 1.0_dp, 1.0_dp)
3873
3874 k = 1
3875 DO
3876 k = k + 1
3877 CALL dbcsr_multiply('N', 'N', 1.0_dp/real(k, dp), a2n, work, 0.0_dp, yk, filter_eps=filter_eps)
3878 CALL dbcsr_multiply('N', 'N', 1.0_dp/real(k, dp), work2, a2n, 0.0_dp, prn, filter_eps=filter_eps)
3879
3880 CALL dbcsr_add(yk, prn, 1.0_dp, 1.0_dp)
3881 CALL dbcsr_add(y, yk, 1.0_dp, 1.0_dp)
3882
3883 IF (dbcsr_frobenius_norm(yk) < threshold) EXIT
3884 CALL dbcsr_copy(work, yk)
3885 CALL dbcsr_copy(work2, prn)
3886 END DO
3887
3888 CALL dbcsr_release(work)
3889 CALL dbcsr_release(work2)
3890 CALL dbcsr_release(prn)
3891 CALL dbcsr_release(a2n)
3892 CALL dbcsr_release(yk)
3893
3894 CALL timestop(handle)
3895
3896 END SUBROUTINE build_y_matrix
3897
3898! **************************************************************************************************
3899!> \brief Overwrites the "optimal" Laplace quadrature with that of the first step
3900!> \param tj ...
3901!> \param wj ...
3902!> \param tau_tj ...
3903!> \param tau_wj ...
3904!> \param weights_cos_tf_t_to_w ...
3905!> \param weights_cos_tf_w_to_t ...
3906!> \param do_laplace ...
3907!> \param do_im_time ...
3908!> \param num_integ_points ...
3909!> \param unit_nr ...
3910!> \param qs_env ...
3911! **************************************************************************************************
3912 SUBROUTINE keep_initial_quad(tj, wj, tau_tj, tau_wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, &
3913 do_laplace, do_im_time, num_integ_points, unit_nr, qs_env)
3914
3915 REAL(dp), ALLOCATABLE, DIMENSION(:), INTENT(INOUT) :: tj, wj, tau_tj, tau_wj
3916 REAL(dp), ALLOCATABLE, DIMENSION(:, :), &
3917 INTENT(INOUT) :: weights_cos_tf_t_to_w, &
3918 weights_cos_tf_w_to_t
3919 LOGICAL, INTENT(IN) :: do_laplace, do_im_time
3920 INTEGER, INTENT(IN) :: num_integ_points, unit_nr
3921 TYPE(qs_environment_type), POINTER :: qs_env
3922
3923 INTEGER :: jquad
3924
3925 IF (do_laplace .OR. do_im_time) THEN
3926 IF (.NOT. ASSOCIATED(qs_env%mp2_env%ri_rpa_im_time%tau_tj)) THEN
3927 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%tau_tj(num_integ_points))
3928 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%tau_wj(num_integ_points))
3929 qs_env%mp2_env%ri_rpa_im_time%tau_tj(:) = tau_tj(:)
3930 qs_env%mp2_env%ri_rpa_im_time%tau_wj(:) = tau_wj(:)
3931 ELSE
3932 !If weights already stored, we overwrite the new ones
3933 tau_tj(:) = qs_env%mp2_env%ri_rpa_im_time%tau_tj(:)
3934 tau_wj(:) = qs_env%mp2_env%ri_rpa_im_time%tau_wj(:)
3935 END IF
3936 END IF
3937 IF (.NOT. do_laplace) THEN
3938 IF (.NOT. ASSOCIATED(qs_env%mp2_env%ri_rpa_im_time%tj)) THEN
3939 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%tj(num_integ_points))
3940 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%wj(num_integ_points))
3941 qs_env%mp2_env%ri_rpa_im_time%tj(:) = tj(:)
3942 qs_env%mp2_env%ri_rpa_im_time%wj(:) = wj(:)
3943 IF (do_im_time) THEN
3944 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
3945 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_w_to_t(num_integ_points, num_integ_points))
3946 qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_t_to_w(:, :) = weights_cos_tf_t_to_w(:, :)
3947 qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_w_to_t(:, :) = weights_cos_tf_w_to_t(:, :)
3948 END IF
3949 ELSE
3950 tj(:) = qs_env%mp2_env%ri_rpa_im_time%tj(:)
3951 wj(:) = qs_env%mp2_env%ri_rpa_im_time%wj(:)
3952 IF (do_im_time) THEN
3953 weights_cos_tf_t_to_w(:, :) = qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_t_to_w(:, :)
3954 weights_cos_tf_w_to_t(:, :) = qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_w_to_t(:, :)
3955 END IF
3956 END IF
3957 END IF
3958 IF (unit_nr > 0) THEN
3959 !Printing order same as in mp2_grids.F for consistency
3960 IF (ASSOCIATED(qs_env%mp2_env%ri_rpa_im_time%tj) .AND. (.NOT. do_laplace)) THEN
3961 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
3962 "MINIMAX_INFO| Number of integration points:", num_integ_points
3963 WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") &
3964 "MINIMAX_INFO| Minimax params (freq grid, scaled):", "Weights", "Abscissas"
3965 DO jquad = 1, num_integ_points
3966 WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") wj(jquad), tj(jquad)
3967 END DO
3968 CALL m_flush(unit_nr)
3969 END IF
3970 IF (ASSOCIATED(qs_env%mp2_env%ri_rpa_im_time%tau_tj)) THEN
3971 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
3972 "MINIMAX_INFO| Number of integration points:", num_integ_points
3973 WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") &
3974 "MINIMAX_INFO| Minimax params (time grid, scaled):", "Weights", "Abscissas"
3975 DO jquad = 1, num_integ_points
3976 WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") tau_wj(jquad), tau_tj(jquad)
3977 END DO
3978 CALL m_flush(unit_nr)
3979 END IF
3980 END IF
3981
3982 END SUBROUTINE keep_initial_quad
3983
3984END 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:76
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:120
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:214
cell_types::pbc
Definition cell_types.F:97
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:84
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
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:21
bibliography::bussy2023
integer, save, public bussy2023
Definition bibliography.F:36
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
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_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_distribution_release
subroutine, public dbcsr_distribution_release(dist)
...
Definition cp_dbcsr_api.F:642
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1178
cp_dbcsr_api::dbcsr_distribution_new
subroutine, public dbcsr_distribution_new(dist, template, group, pgrid, row_dist, col_dist, reuse_arrays)
...
Definition cp_dbcsr_api.F:615
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
Definition cp_dbcsr_api.F:969
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:943
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1040
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:692
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1086
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1002
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1194
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1132
cp_dbcsr_api::dbcsr_complete_redistribute
subroutine, public dbcsr_complete_redistribute(matrix, redist)
...
Definition cp_dbcsr_api.F:319
cp_dbcsr_api::dbcsr_clear
subroutine, public dbcsr_clear(matrix)
...
Definition cp_dbcsr_api.F:304
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:253
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, uplo_to_full)
used to replace the cholesky decomposition by the inverse
Definition cp_dbcsr_cholesky.F:114
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_add_on_diag
subroutine, public dbcsr_add_on_diag(matrix, alpha)
Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
Definition cp_dbcsr_contrib.F:331
cp_dbcsr_contrib::dbcsr_frobenius_norm
real(dp) function, public dbcsr_frobenius_norm(matrix)
Compute the frobenius norm of a dbcsr matrix.
Definition cp_dbcsr_contrib.F:127
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:279
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:552
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:419
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:214
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:111
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:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:169
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:556
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, external_para_env)
Add the hfx contributions to the Hamiltonian.
Definition hfx_admm_utils.F:2487
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, nspins)
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:3368
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:4238
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:3044
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:3454
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:2910
hfx_types::dealloc_containers
subroutine, public dealloc_containers(data, memory_usage)
...
Definition hfx_types.F:2878
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:1259
input_constants::ri_rpa_method_gpw
integer, parameter, public ri_rpa_method_gpw
Definition input_constants.F:1064
input_constants::do_admm_aux_exch_func_none
integer, parameter, public do_admm_aux_exch_func_none
Definition input_constants.F:874
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition input_constants.F:812
input_constants::do_eri_gpw
integer, parameter, public do_eri_gpw
Definition input_constants.F:1259
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:731
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:704
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:1047
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:2108
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:136
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:153
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:1958
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:962
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:716
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:1105
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:1055
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:99
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:1715
qs_collocate_density::collocate_function
subroutine, public collocate_function(vector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, eps_rho_rspace, basis_type)
maps a given function on the grid
Definition qs_collocate_density.F:2920
qs_core_matrices
Calculation of the core Hamiltonian integral matrix <a|H|b> over Cartesian Gaussian-type functions.
Definition qs_core_matrices.F:69
qs_core_matrices::kinetic_energy_matrix
subroutine, public kinetic_energy_matrix(qs_env, matrixkp_t, matrix_t, matrix_p, matrix_name, calculate_forces, nderivative, sab_orb, eps_filter, basis_type, debug_forces, debug_stress)
Calculate kinetic energy matrix and possible relativistic correction.
Definition qs_core_matrices.F:396
qs_core_matrices::core_matrices
subroutine, public core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, ec_env, dcdr_env, ec_env_matrices, basis_type, debug_forces, debug_stress, atcore)
...
Definition qs_core_matrices.F:142
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:407
qs_environment_types
Definition qs_environment_types.F:14
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_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, harris_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, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:526
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, rhoz_cneo_set, 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, harris_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, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1098
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:439
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
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:956
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_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition qs_ks_types.F:578
qs_linres_types
Type definitiona for linear response calculations.
Definition qs_linres_types.F:12
qs_matrix_w
Utility subroutine for qs energy calculation.
Definition qs_matrix_w.F:14
qs_matrix_w::compute_matrix_w
subroutine, public compute_matrix_w(qs_env, calc_forces)
Refactoring of qs_energies_scf. Moves computation of matrix_w into separate subroutine.
Definition qs_matrix_w.F:62
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:399
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:1095
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:479
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:130
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:3156
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:1626
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:926
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:144
qs_tensors::compress_tensor
subroutine, public compress_tensor(tensor, blk_indices, compressed, eps, memory)
...
Definition qs_tensors.F:3405
qs_tensors::neighbor_list_3c_destroy
subroutine, public neighbor_list_3c_destroy(ijk_list)
Destroy 3c neighborlist.
Definition qs_tensors.F:381
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:2977
qs_tensors::decompress_tensor
subroutine, public decompress_tensor(tensor, blk_indices, compressed, eps)
...
Definition qs_tensors.F:3533
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:2735
qs_tensors::get_tensor_occupancy
subroutine, public get_tensor_occupancy(tensor, nze, occ)
...
Definition qs_tensors.F:3624
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:280
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:2144
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:560
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:2217
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:3914
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:664
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:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
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:628
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
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:115
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:368
hfx_types::hfx_compression_type
Definition hfx_types.F:352
hfx_types::hfx_type
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F:510
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
libint_2c_3c::libint_potential_type
Definition libint_2c_3c.F:67
message_passing::mp_cart_type
Definition message_passing.F:698
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:721
mp2_types::mp2_type
Definition mp2_types.F:345
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:123
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:220
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:177
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:47
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:136
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