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