(git:f8a708f)
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, 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
1377 !calculate (1 + Q(w))^-1 - 1 for the given freq.
1378 !Always take spin alpha (get 2*alpha in closed shell, and alpha+beta in open-shell)
1379 CALL dbcsr_copy(dbcsr_work_symm, mat_p_omega(iquad, 1)%matrix)
1380 CALL dbt_copy_matrix_to_tensor(dbcsr_work_symm, t_2c_tmp)
1381 CALL dbt_copy(t_2c_tmp, t_2c_ri, move_data=.true.)
1382
1383 CALL dbt_contract(1.0_dp, t_2c_ri, force_data%t_2c_K, 0.0_dp, t_2c_ri_2, &
1384 contract_1=[2], notcontract_1=[1], &
1385 contract_2=[1], notcontract_2=[2], &
1386 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1387 flop=flop, unit_nr=unit_nr_dbcsr)
1388 dbcsr_nflop = dbcsr_nflop + flop
1389 CALL dbt_contract(1.0_dp, force_data%t_2c_K, t_2c_ri_2, 0.0_dp, t_2c_ri, &
1390 contract_1=[1], notcontract_1=[2], &
1391 contract_2=[1], notcontract_2=[2], &
1392 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1393 flop=flop, unit_nr=unit_nr_dbcsr)
1394 CALL dbt_copy(t_2c_ri, t_2c_tmp, move_data=.true.)
1395 CALL dbt_copy_tensor_to_matrix(t_2c_tmp, dbcsr_work_symm)
1396 CALL dbcsr_add_on_diag(dbcsr_work_symm, 1.0_dp)
1397
1398 CALL cp_dbcsr_cholesky_decompose(dbcsr_work_symm, para_env=para_env, blacs_env=blacs_env)
1399 CALL cp_dbcsr_cholesky_invert(dbcsr_work_symm, para_env=para_env, blacs_env=blacs_env, uplo_to_full=.true.)
1400
1401 CALL dbcsr_add_on_diag(dbcsr_work_symm, -1.0_dp)
1402
1403 DO jquad = 1, num_integ_points
1404
1405 !the P matrix to time.
1406 weight = grid%cosine_frequency_to_time(jquad, iquad)
1407 IF (open_shell) THEN
1408 IF (ispin == 1) THEN
1409 !mat_P_omega contains the sum of alpha and beta spin => we only want alpha
1410 CALL dbcsr_add(mat_p_tau(jquad)%matrix, mat_p_omega(iquad, 1)%matrix, 1.0_dp, weight)
1411 CALL dbcsr_add(mat_p_tau(jquad)%matrix, mat_p_omega(iquad, 2)%matrix, 1.0_dp, -weight)
1412 ELSE
1413 CALL dbcsr_add(mat_p_tau(jquad)%matrix, mat_p_omega(iquad, 2)%matrix, 1.0_dp, weight)
1414 END IF
1415 ELSE
1416 !factor 0.5 because originam matrix Q is scaled by 2 in RPA (spin)
1417 weight = 0.5_dp*weight
1418 CALL dbcsr_add(mat_p_tau(jquad)%matrix, mat_p_omega(iquad, 1)%matrix, 1.0_dp, weight)
1419 END IF
1420
1421 ! The energy derivative uses the weighted adjoint of the forward transform.
1422 weight = grid%cosine_time_to_frequency(iquad, jquad)* &
1423 grid%frequency_weights(iquad)
1424 CALL dbt_copy_matrix_to_tensor(dbcsr_work_symm, t_2c_tmp)
1425 CALL dbt_scale(t_2c_tmp, weight)
1426 CALL dbt_copy(t_2c_tmp, t_b(jquad), summation=.true., move_data=.true.)
1427 END DO
1428 END DO
1429 CALL dbt_destroy(t_2c_tmp)
1430 CALL dbcsr_release(dbcsr_work_symm)
1431 CALL dbt_clear(t_2c_ri)
1432 CALL dbt_clear(t_2c_ri_2)
1433
1434 !Pre-define the sparsity of t_3c_4 as a function of the derivatives
1435 occ_der_ao = 0; nze_der_ao = 0
1436 occ_der_ri = 0; nze_der_ri = 0
1437 DO i_xyz = 1, 3
1438 DO i_mem = 1, cut_memory
1439 CALL decompress_tensor(force_data%t_3c_der_RI(i_xyz), force_data%t_3c_der_RI_ind(i_mem, i_xyz)%ind, &
1440 force_data%t_3c_der_RI_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
1441 CALL get_tensor_occupancy(force_data%t_3c_der_RI(i_xyz), nze, occ)
1442 occ_der_ri = occ_der_ri + occ
1443 nze_der_ri = nze_der_ri + nze
1444 CALL dbt_copy(force_data%t_3c_der_RI(i_xyz), t_3c_sparse, summation=.true., move_data=.true.)
1445
1446 CALL decompress_tensor(force_data%t_3c_der_AO(i_xyz), force_data%t_3c_der_AO_ind(i_mem, i_xyz)%ind, &
1447 force_data%t_3c_der_AO_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
1448 CALL get_tensor_occupancy(force_data%t_3c_der_AO(i_xyz), nze, occ)
1449 occ_der_ao = occ_der_ao + occ
1450 nze_der_ao = nze_der_ao + nze
1451 CALL dbt_copy(force_data%t_3c_der_AO(i_xyz), t_3c_sparse, order=[1, 3, 2], summation=.true.)
1452 CALL dbt_copy(force_data%t_3c_der_AO(i_xyz), t_3c_sparse, summation=.true., move_data=.true.)
1453 END DO
1454 END DO
1455 occ_der_ri = occ_der_ri/3.0_dp
1456 occ_der_ao = occ_der_ao/3.0_dp
1457 nze_der_ri = nze_der_ri/3
1458 nze_der_ao = nze_der_ao/3
1459
1460 CALL dbcsr_create(r_occ, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1461 CALL dbcsr_create(r_virt, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1462 CALL dbcsr_create(dbcsr_work_symm, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_symmetric)
1463 CALL dbcsr_create(dbcsr_work1, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1464 CALL dbcsr_create(dbcsr_work2, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1465 CALL dbcsr_create(dbcsr_work3, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1466 CALL dbcsr_create(exp_occ, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1467 CALL dbcsr_create(exp_virt, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1468 IF (use_virial) CALL dbcsr_create(virial_ovlp, template=dbcsr_work1)
1469
1470 CALL dbt_batched_contract_init(t_3c_0, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
1471 CALL dbt_batched_contract_init(t_3c_1, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
1472 CALL dbt_batched_contract_init(t_3c_3, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1473 CALL dbt_batched_contract_init(t_m_occ, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1474 CALL dbt_batched_contract_init(t_m_virt, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1475
1476 CALL dbt_batched_contract_init(t_3c_ints, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1477 CALL dbt_batched_contract_init(t_3c_work, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1478
1479 CALL dbt_batched_contract_init(t_3c_4, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1480 batch_range_3=mc_ranges)
1481 CALL dbt_batched_contract_init(t_3c_5, 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_6, 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_7, 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_8, 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_sparse, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1490 batch_range_3=mc_ranges)
1491
1492 fac = 1.0_dp/fourpi*mp2_env%ri_rpa%scale_rpa
1493 IF (open_shell) fac = 0.5_dp*fac
1494
1495 work_virial = 0.0_dp
1496 work_virial_ovlp = 0.0_dp
1497 DO jquad = 1, num_integ_points
1498 tau = grid%imaginary_time(jquad)
1499 occ_ddint = 0; nze_ddint = 0
1500
1501 CALL para_env%sync()
1502 t1 = m_walltime()
1503
1504 !Deal with the force contributions where there is no explicit 3-center quantities, i.e. the
1505 !forces due to the metric and potential derivatives
1506 CALL dbt_create(mat_p_tau(jquad)%matrix, t_2c_tmp)
1507 CALL dbt_copy_matrix_to_tensor(mat_p_tau(jquad)%matrix, t_2c_tmp)
1508 CALL dbt_copy(t_2c_tmp, t_p, move_data=.true.)
1509 CALL dbt_filter(t_p, eps_filter)
1510 CALL dbt_destroy(t_2c_tmp)
1511
1512 CALL perform_2c_ops(force, t_kbkt, force_data, fac, t_b(jquad), t_p, t_2c_ri, t_2c_ri_2, &
1513 use_virial, atom_of_kind, kind_of, eps_filter, dbcsr_nflop, unit_nr_dbcsr)
1514 CALL get_tensor_occupancy(t_kbkt, nze_kbk, occ_kbk)
1515
1516 !Calculate the pseudo-density matrix in tensor form. There are a few useless arguments for SOS-MP2
1517 CALL compute_mat_dm_global(grid, nmo, fm_mo_coeff(ispin), homo(ispin), propagator, &
1518 matrix_s, ispin, eigenval(:, ispin), e_fermi, eps_filter, &
1519 mp2_env%ri_rpa_im_time%memory_info, unit_nr, &
1520 jquad, .false., .false., qs_env, dummy_int, dummy_ptr, para_env)
1521
1522 CALL dbt_create(propagator(propagator_sector_occupied, jquad, 1)%matrix, t_2c_tmp)
1523 CALL dbt_copy_matrix_to_tensor(propagator(propagator_sector_occupied, jquad, 1)%matrix, t_2c_tmp)
1524 CALL dbt_copy(t_2c_tmp, t_dm_occ, move_data=.true.)
1525 CALL dbt_filter(t_dm_occ, eps_filter)
1526 CALL dbt_destroy(t_2c_tmp)
1527
1528 CALL dbt_create(propagator(propagator_sector_virtual, jquad, 1)%matrix, t_2c_tmp)
1529 CALL dbt_copy_matrix_to_tensor(propagator(propagator_sector_virtual, jquad, 1)%matrix, t_2c_tmp)
1530 CALL dbt_copy(t_2c_tmp, t_dm_virt, move_data=.true.)
1531 CALL dbt_filter(t_dm_virt, eps_filter)
1532 CALL dbt_destroy(t_2c_tmp)
1533
1534 !Deal with the 3-center quantities.
1535 CALL perform_3c_ops(force, t_r_occ, t_r_virt, force_data, fac, cut_memory, n_mem_ri, &
1536 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, &
1537 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, &
1538 t_3c_ints, t_3c_work, starts_array_mc, ends_array_mc, batch_start_ri, &
1539 batch_end_ri, t_3c_o_compressed, t_3c_o_ind, use_virial, &
1540 atom_of_kind, kind_of, eps_filter, occ_ddint, nze_ddint, dbcsr_nflop, &
1541 unit_nr_dbcsr, mp2_env)
1542
1543 CALL timeset(routinen//"_dbcsr", handle2)
1544 !We go back to DBCSR matrices from now on
1545 !Note: R matrices are in fact symmetric, but use a normal type for convenience
1546 CALL dbt_create(matrix_s(1)%matrix, t_2c_tmp)
1547 CALL dbt_copy(t_r_occ, t_2c_tmp, move_data=.true.)
1548 CALL dbt_copy_tensor_to_matrix(t_2c_tmp, r_occ)
1549
1550 CALL dbt_copy(t_r_virt, t_2c_tmp, move_data=.true.)
1551 CALL dbt_copy_tensor_to_matrix(t_2c_tmp, r_virt)
1552
1553 !Iteratively calculate the Y1 and Y2 matrices
1554 CALL dbcsr_copy(dbcsr_work_symm, matrix_ks(ispin)%matrix)
1555 CALL dbcsr_add(dbcsr_work_symm, matrix_s(1)%matrix, 1.0_dp, -e_fermi)
1556 CALL dbcsr_multiply('N', 'N', tau, force_data%P_occ(ispin)%matrix, &
1557 dbcsr_work_symm, 0.0_dp, dbcsr_work1)
1558 CALL build_y_matrix(y_1, dbcsr_work1, force_data%P_occ(ispin)%matrix, r_virt, eps_filter)
1559 CALL matrix_exponential(exp_occ, dbcsr_work1, 1.0_dp, 1.0_dp, eps_filter)
1560
1561 CALL dbcsr_multiply('N', 'N', -tau, force_data%P_virt(ispin)%matrix, &
1562 dbcsr_work_symm, 0.0_dp, dbcsr_work1)
1563 CALL build_y_matrix(y_2, dbcsr_work1, force_data%P_virt(ispin)%matrix, r_occ, eps_filter)
1564 CALL matrix_exponential(exp_virt, dbcsr_work1, 1.0_dp, 1.0_dp, eps_filter)
1565
1566 !The force contribution coming from [-S^-1*(e^-tau*P_virt*F)^T*R_occ*S^-1
1567 ! +tau*S^-1*Y_2^T*F*S^-1] * der_S
1568 !as well as -tau*e_fermi*Y_1*P^occ + tau*e_fermi*Y_2*P^virt
1569 CALL dbcsr_multiply('N', 'N', 1.0_dp, r_occ, force_data%inv_ovlp, 0.0_dp, dbcsr_work1)
1570 CALL dbcsr_multiply('T', 'N', 1.0_dp, exp_virt, dbcsr_work1, 0.0_dp, dbcsr_work3)
1571 CALL dbcsr_multiply('N', 'N', 1.0_dp, force_data%inv_ovlp, dbcsr_work3, 0.0_dp, dbcsr_work2)
1572
1573 CALL dbcsr_multiply('N', 'T', tau, force_data%inv_ovlp, y_2, 0.0_dp, dbcsr_work3)
1574 CALL dbcsr_multiply('N', 'N', 1.0_dp, dbcsr_work3, dbcsr_work_symm, 0.0_dp, dbcsr_work1)
1575 CALL dbcsr_multiply('N', 'N', -1.0_dp, dbcsr_work1, force_data%inv_ovlp, 1.0_dp, dbcsr_work2)
1576
1577 CALL dbcsr_multiply('N', 'T', tau*e_fermi, force_data%P_occ(ispin)%matrix, y_1, 1.0_dp, dbcsr_work2)
1578 CALL dbcsr_multiply('N', 'T', -tau*e_fermi, force_data%P_virt(ispin)%matrix, y_2, 1.0_dp, dbcsr_work2)
1579
1580 CALL dbt_copy_matrix_to_tensor(dbcsr_work2, t_2c_tmp)
1581 CALL dbt_copy(t_2c_tmp, t_2c_ao, move_data=.true.)
1582
1583 pref = -1.0_dp*fac
1584 CALL get_2c_der_force(force, t_2c_ao, force_data%t_2c_der_ovlp, atom_of_kind, &
1585 kind_of, force_data%idx_to_at_AO, pref, do_ovlp=.true.)
1586
1587 IF (use_virial) CALL dbcsr_add(virial_ovlp, dbcsr_work2, 1.0_dp, pref)
1588
1589 !The final contribution from Tr[(tau*Y_1*P_occ - tau*Y_2*P_virt) * der_F]
1590 CALL dbcsr_multiply('N', 'N', fac*tau, y_1, force_data%P_occ(ispin)%matrix, 1.0_dp, &
1591 force_data%sum_YP_tau(ispin)%matrix, retain_sparsity=.true.)
1592 CALL dbcsr_multiply('N', 'N', -fac*tau, y_2, force_data%P_virt(ispin)%matrix, 1.0_dp, &
1593 force_data%sum_YP_tau(ispin)%matrix, retain_sparsity=.true.)
1594
1595 spin_fac = 0.5_dp*fac
1596 IF (open_shell) spin_fac = 2.0_dp*spin_fac
1597 !Build-up the RHS of the response equation.
1598 CALL dbcsr_multiply('N', 'N', 1.0_dp*spin_fac, r_virt, exp_occ, 1.0_dp, &
1599 force_data%sum_O_tau(ispin)%matrix, retain_sparsity=.true.)
1600 CALL dbcsr_multiply('N', 'N', -1.0_dp*spin_fac, r_occ, exp_virt, 1.0_dp, &
1601 force_data%sum_O_tau(ispin)%matrix, retain_sparsity=.true.)
1602 CALL dbcsr_multiply('N', 'N', tau*spin_fac, dbcsr_work_symm, y_1, 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_2, 1.0_dp, &
1605 force_data%sum_O_tau(ispin)%matrix, retain_sparsity=.true.)
1606
1607 CALL timestop(handle2)
1608
1609 !Print some info
1610 CALL para_env%sync()
1611 t2 = m_walltime()
1612 dbcsr_time = dbcsr_time + t2 - t1
1613
1614 IF (unit_nr > 0) THEN
1615 WRITE (unit_nr, '(/T3,A,1X,I3,A)') &
1616 'RPA_LOW_SCALING_INFO| Info for time point', jquad, ' (gradients)'
1617 WRITE (unit_nr, '(T6,A,T56,F25.6)') &
1618 'Time:', t2 - t1
1619 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1620 'Occupancy of 3c AO derivs:', real(nze_der_ao, dp), '/', occ_der_ao*100, '%'
1621 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1622 'Occupancy of 3c RI derivs:', real(nze_der_ri, dp), '/', occ_der_ri*100, '%'
1623 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1624 'Occupancy of the Docc * Dvirt * 3c-int tensor', real(nze_ddint, dp), '/', occ_ddint*100, '%'
1625 WRITE (unit_nr, '(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1626 'Occupancy of KBK^T 2c-tensor:', real(nze_kbk, dp), '/', occ_kbk*100, '%'
1627 CALL m_flush(unit_nr)
1628 END IF
1629
1630 !intermediate clean-up
1631 CALL dbcsr_release(y_1)
1632 CALL dbcsr_release(y_2)
1633 CALL dbt_destroy(t_2c_tmp)
1634
1635 END DO !jquad
1636
1637 CALL dbt_batched_contract_finalize(t_3c_0)
1638 CALL dbt_batched_contract_finalize(t_3c_1)
1639 CALL dbt_batched_contract_finalize(t_3c_3)
1640 CALL dbt_batched_contract_finalize(t_m_occ)
1641 CALL dbt_batched_contract_finalize(t_m_virt)
1642
1643 CALL dbt_batched_contract_finalize(t_3c_ints)
1644 CALL dbt_batched_contract_finalize(t_3c_work)
1645
1646 CALL dbt_batched_contract_finalize(t_3c_4)
1647 CALL dbt_batched_contract_finalize(t_3c_5)
1648 CALL dbt_batched_contract_finalize(t_3c_6)
1649 CALL dbt_batched_contract_finalize(t_3c_7)
1650 CALL dbt_batched_contract_finalize(t_3c_8)
1651 CALL dbt_batched_contract_finalize(t_3c_sparse)
1652
1653 !Calculate the 2c and 3c contributions to the virial
1654 IF (use_virial) THEN
1655 CALL dbt_copy(force_data%t_3c_virial_split, force_data%t_3c_virial, move_data=.true.)
1656 CALL calc_3c_virial(work_virial, force_data%t_3c_virial, 1.0_dp, qs_env, force_data%nl_3c, &
1657 basis_set_ri_aux, basis_set_ao, basis_set_ao, mp2_env%ri_metric, &
1658 der_eps=mp2_env%ri_rpa_im_time%eps_filter, op_pos=1)
1659
1660 CALL calc_2c_virial(work_virial, force_data%RI_virial_met, 1.0_dp, qs_env, force_data%nl_2c_met, &
1661 basis_set_ri_aux, basis_set_ri_aux, mp2_env%ri_metric)
1662 CALL dbcsr_clear(force_data%RI_virial_met)
1663
1664 IF (.NOT. force_data%do_periodic) THEN
1665 CALL calc_2c_virial(work_virial, force_data%RI_virial_pot, 1.0_dp, qs_env, force_data%nl_2c_pot, &
1666 basis_set_ri_aux, basis_set_ri_aux, mp2_env%potential_parameter)
1667 CALL dbcsr_clear(force_data%RI_virial_pot)
1668 END IF
1669
1670 identity_pot%potential_type = do_potential_id
1671 CALL calc_2c_virial(work_virial_ovlp, virial_ovlp, 1.0_dp, qs_env, force_data%nl_2c_ovlp, &
1672 basis_set_ao, basis_set_ao, identity_pot)
1673 CALL dbcsr_release(virial_ovlp)
1674
1675 DO k_xyz = 1, 3
1676 DO j_xyz = 1, 3
1677 DO i_xyz = 1, 3
1678 virial%pv_mp2(i_xyz, j_xyz) = virial%pv_mp2(i_xyz, j_xyz) &
1679 - work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1680 virial%pv_overlap(i_xyz, j_xyz) = virial%pv_overlap(i_xyz, j_xyz) &
1681 - work_virial_ovlp(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1682 virial%pv_virial(i_xyz, j_xyz) = virial%pv_virial(i_xyz, j_xyz) &
1683 - work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz) &
1684 - work_virial_ovlp(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1685 END DO
1686 END DO
1687 END DO
1688 END IF
1689
1690 !Calculate the periodic contributions of (P|Q) to the force and the virial
1691 work_virial = 0.0_dp
1692 IF (force_data%do_periodic) THEN
1693 IF (mp2_env%eri_method == do_eri_gpw) THEN
1694 CALL get_2c_gpw_forces(force_data%G_PQ, force, work_virial, use_virial, mp2_env, qs_env)
1695 ELSE IF (mp2_env%eri_method == do_eri_mme) THEN
1696 CALL get_2c_mme_forces(force_data%G_PQ, force, mp2_env, qs_env)
1697 IF (use_virial) cpabort("Stress tensor not available with MME intrgrals")
1698 ELSE
1699 cpabort("Periodic case not possible with OS integrals")
1700 END IF
1701 CALL dbcsr_clear(force_data%G_PQ)
1702 END IF
1703
1704 IF (use_virial) THEN
1705 virial%pv_mp2 = virial%pv_mp2 + work_virial
1706 virial%pv_virial = virial%pv_virial + work_virial
1707 virial%pv_calculate = .false.
1708
1709 DO ibasis = 1, SIZE(basis_set_ao)
1710 orb_basis => basis_set_ao(ibasis)%gto_basis_set
1711 CALL init_interaction_radii_orb_basis(orb_basis, eps_pgf_orb_old)
1712 ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
1713 CALL init_interaction_radii_orb_basis(ri_basis, eps_pgf_orb_old)
1714 END DO
1715 END IF
1716
1717 !clean-up
1718 IF (ASSOCIATED(dummy_ptr)) DEALLOCATE (dummy_ptr)
1719 DO jquad = 1, num_integ_points
1720 CALL dbt_destroy(t_b(jquad))
1721 END DO
1722 CALL dbt_destroy(t_p)
1723 CALL dbt_destroy(t_3c_0)
1724 CALL dbt_destroy(t_3c_1)
1725 CALL dbt_destroy(t_3c_3)
1726 CALL dbt_destroy(t_3c_4)
1727 CALL dbt_destroy(t_3c_5)
1728 CALL dbt_destroy(t_3c_6)
1729 CALL dbt_destroy(t_3c_7)
1730 CALL dbt_destroy(t_3c_8)
1731 CALL dbt_destroy(t_3c_sparse)
1732 CALL dbt_destroy(t_3c_help_1)
1733 CALL dbt_destroy(t_3c_help_2)
1734 CALL dbt_destroy(t_3c_ints)
1735 CALL dbt_destroy(t_3c_work)
1736 CALL dbt_destroy(t_r_occ)
1737 CALL dbt_destroy(t_r_virt)
1738 CALL dbt_destroy(t_dm_occ)
1739 CALL dbt_destroy(t_dm_virt)
1740 CALL dbt_destroy(t_kbkt)
1741 CALL dbt_destroy(t_m_occ)
1742 CALL dbt_destroy(t_m_virt)
1743 CALL dbcsr_release(r_occ)
1744 CALL dbcsr_release(r_virt)
1745 CALL dbcsr_release(dbcsr_work_symm)
1746 CALL dbcsr_release(dbcsr_work1)
1747 CALL dbcsr_release(dbcsr_work2)
1748 CALL dbcsr_release(dbcsr_work3)
1749 CALL dbcsr_release(exp_occ)
1750 CALL dbcsr_release(exp_virt)
1751
1752 CALL dbt_destroy(t_2c_ri)
1753 CALL dbt_destroy(t_2c_ri_2)
1754 CALL dbt_destroy(t_2c_ao)
1755 CALL dbcsr_deallocate_matrix_set(propagator)
1756 CALL dbcsr_deallocate_matrix_set(mat_p_tau)
1757
1758 CALL timestop(handle)
1759
1760 END SUBROUTINE calc_rpa_loop_forces
1761
1762! **************************************************************************************************
1763!> \brief This subroutines performs the 2c tensor operations that are common accros low-scaling RPA
1764!> and SOS-MP2, including forces and virial
1765!> \param force ...
1766!> \param t_KBKT returns the 2c tensor product of K*B*K^T
1767!> \param force_data ...
1768!> \param fac ...
1769!> \param t_B depending on RPA or SOS-MP2, t_B contains (1 + Q)^-1 - 1 or simply Q, respectively
1770!> \param t_P ...
1771!> \param t_2c_RI ...
1772!> \param t_2c_RI_2 ...
1773!> \param use_virial ...
1774!> \param atom_of_kind ...
1775!> \param kind_of ...
1776!> \param eps_filter ...
1777!> \param dbcsr_nflop ...
1778!> \param unit_nr_dbcsr ...
1779! **************************************************************************************************
1780 SUBROUTINE perform_2c_ops(force, t_KBKT, force_data, fac, t_B, t_P, t_2c_RI, t_2c_RI_2, use_virial, &
1781 atom_of_kind, kind_of, eps_filter, dbcsr_nflop, unit_nr_dbcsr)
1782
1783 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1784 TYPE(dbt_type), INTENT(INOUT) :: t_kbkt
1785 TYPE(im_time_force_type), INTENT(INOUT) :: force_data
1786 REAL(dp), INTENT(IN) :: fac
1787 TYPE(dbt_type), INTENT(INOUT) :: t_b, t_p, t_2c_ri, t_2c_ri_2
1788 LOGICAL, INTENT(IN) :: use_virial
1789 INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind, kind_of
1790 REAL(dp), INTENT(IN) :: eps_filter
1791 INTEGER(int_8), INTENT(INOUT) :: dbcsr_nflop
1792 INTEGER, INTENT(IN) :: unit_nr_dbcsr
1793
1794 CHARACTER(LEN=*), PARAMETER :: routinen = 'perform_2c_ops'
1795
1796 INTEGER :: handle
1797 INTEGER(int_8) :: flop
1798 REAL(dp) :: pref
1799 TYPE(dbt_type) :: t_2c_tmp, t_2c_virial
1800
1801 CALL timeset(routinen, handle)
1802
1803 IF (use_virial) CALL dbt_create(force_data%RI_virial_pot, t_2c_virial)
1804
1805 !P^T*K*B + P*K*B^T (note we calculate and save K*B*K^T for later, and P=P^T)
1806 CALL dbt_contract(1.0_dp, force_data%t_2c_K, t_b, 0.0_dp, t_2c_ri, &
1807 contract_1=[2], notcontract_1=[1], &
1808 contract_2=[1], notcontract_2=[2], &
1809 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1810 flop=flop, unit_nr=unit_nr_dbcsr)
1811 dbcsr_nflop = dbcsr_nflop + flop
1812
1813 CALL dbt_contract(1.0_dp, t_2c_ri, force_data%t_2c_K, 0.0_dp, t_kbkt, &
1814 contract_1=[2], notcontract_1=[1], &
1815 contract_2=[2], notcontract_2=[1], &
1816 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1817 flop=flop, unit_nr=unit_nr_dbcsr)
1818 dbcsr_nflop = dbcsr_nflop + flop
1819
1820 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
1821 contract_1=[2], notcontract_1=[1], &
1822 contract_2=[1], notcontract_2=[2], &
1823 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1824 flop=flop, unit_nr=unit_nr_dbcsr)
1825 dbcsr_nflop = dbcsr_nflop + flop
1826 CALL dbt_clear(t_2c_ri)
1827 !t_2c_RI_2 currently holds 2*P^T*K*B = P^T*K*B + P*K*B^T (because of symmetry)
1828
1829 !For the metric contribution, we need S^-1*(P^T*K*B + P*K*B^T)*K^T
1830 CALL dbt_contract(1.0_dp, force_data%t_2c_inv_metric, t_2c_ri_2, 0.0_dp, t_2c_ri, &
1831 contract_1=[2], notcontract_1=[1], &
1832 contract_2=[1], notcontract_2=[2], &
1833 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1834 flop=flop, unit_nr=unit_nr_dbcsr)
1835 dbcsr_nflop = dbcsr_nflop + flop
1836
1837 CALL dbt_contract(1.0_dp, t_2c_ri, force_data%t_2c_K, 0.0_dp, t_2c_ri_2, &
1838 contract_1=[2], notcontract_1=[1], &
1839 contract_2=[2], notcontract_2=[1], &
1840 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1841 flop=flop, unit_nr=unit_nr_dbcsr)
1842 dbcsr_nflop = dbcsr_nflop + flop
1843
1844 !Here we do the trace for the force
1845 pref = -1.0_dp*fac
1846 CALL get_2c_der_force(force, t_2c_ri_2, force_data%t_2c_der_metric, atom_of_kind, &
1847 kind_of, force_data%idx_to_at_RI, pref, do_mp2=.true.)
1848 IF (use_virial) THEN
1849 CALL dbt_copy(t_2c_ri_2, t_2c_virial)
1850 CALL dbt_scale(t_2c_virial, pref)
1851 CALL dbt_copy_tensor_to_matrix(t_2c_virial, force_data%RI_virial_met, summation=.true.)
1852 CALL dbt_clear(t_2c_virial)
1853 END IF
1854
1855 !For the potential contribution, we need S^-1*(P^T*K*B + P*K*B^T)*V^-0.5
1856 !some of it is still in t_2c_RI: ( S^-1*(P^T*K*B + P*K*B^T) )
1857 CALL dbt_contract(1.0_dp, t_2c_ri, force_data%t_2c_pot_msqrt, 0.0_dp, t_2c_ri_2, &
1858 contract_1=[2], notcontract_1=[1], &
1859 contract_2=[1], notcontract_2=[2], &
1860 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1861 flop=flop, unit_nr=unit_nr_dbcsr)
1862 dbcsr_nflop = dbcsr_nflop + flop
1863
1864 !Here we do the trace for the force. In the periodic case, we store the matrix in G_PQ for later
1865 pref = 0.5_dp*fac
1866 IF (force_data%do_periodic) THEN
1867 CALL dbt_scale(t_2c_ri_2, pref)
1868 CALL dbt_create(force_data%G_PQ, t_2c_tmp)
1869 CALL dbt_copy(t_2c_ri_2, t_2c_tmp, move_data=.true.)
1870 CALL dbt_copy_tensor_to_matrix(t_2c_tmp, force_data%G_PQ, summation=.true.)
1871 CALL dbt_destroy(t_2c_tmp)
1872 ELSE
1873 CALL get_2c_der_force(force, t_2c_ri_2, force_data%t_2c_der_pot, atom_of_kind, &
1874 kind_of, force_data%idx_to_at_RI, pref, do_mp2=.true.)
1875
1876 IF (use_virial) THEN
1877 CALL dbt_copy(t_2c_ri_2, t_2c_virial)
1878 CALL dbt_scale(t_2c_virial, pref)
1879 CALL dbt_copy_tensor_to_matrix(t_2c_virial, force_data%RI_virial_pot, summation=.true.)
1880 CALL dbt_clear(t_2c_virial)
1881 END IF
1882 END IF
1883
1884 CALL dbt_clear(t_2c_ri)
1885 CALL dbt_clear(t_2c_ri_2)
1886
1887 IF (use_virial) CALL dbt_destroy(t_2c_virial)
1888
1889 CALL timestop(handle)
1890
1891 END SUBROUTINE perform_2c_ops
1892
1893! **************************************************************************************************
1894!> \brief This subroutines performs the 3c tensor operations that are common accros low-scaling RPA
1895!> and SOS-MP2, including forces and virial
1896!> \param force ...
1897!> \param t_R_occ ...
1898!> \param t_R_virt ...
1899!> \param force_data ...
1900!> \param fac ...
1901!> \param cut_memory ...
1902!> \param n_mem_RI ...
1903!> \param t_KBKT ...
1904!> \param t_dm_occ ...
1905!> \param t_dm_virt ...
1906!> \param t_3c_O ...
1907!> \param t_3c_M ...
1908!> \param t_M_occ ...
1909!> \param t_M_virt ...
1910!> \param t_3c_0 ...
1911!> \param t_3c_1 ...
1912!> \param t_3c_3 ...
1913!> \param t_3c_4 ...
1914!> \param t_3c_5 ...
1915!> \param t_3c_6 ...
1916!> \param t_3c_7 ...
1917!> \param t_3c_8 ...
1918!> \param t_3c_sparse ...
1919!> \param t_3c_help_1 ...
1920!> \param t_3c_help_2 ...
1921!> \param t_3c_ints ...
1922!> \param t_3c_work ...
1923!> \param starts_array_mc ...
1924!> \param ends_array_mc ...
1925!> \param batch_start_RI ...
1926!> \param batch_end_RI ...
1927!> \param t_3c_O_compressed ...
1928!> \param t_3c_O_ind ...
1929!> \param use_virial ...
1930!> \param atom_of_kind ...
1931!> \param kind_of ...
1932!> \param eps_filter ...
1933!> \param occ_ddint ...
1934!> \param nze_ddint ...
1935!> \param dbcsr_nflop ...
1936!> \param unit_nr_dbcsr ...
1937!> \param mp2_env ...
1938! **************************************************************************************************
1939 SUBROUTINE perform_3c_ops(force, t_R_occ, t_R_virt, force_data, fac, cut_memory, n_mem_RI, &
1940 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, &
1941 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, &
1942 t_3c_ints, t_3c_work, starts_array_mc, ends_array_mc, batch_start_RI, &
1943 batch_end_RI, t_3c_O_compressed, t_3c_O_ind, use_virial, &
1944 atom_of_kind, kind_of, eps_filter, occ_ddint, nze_ddint, dbcsr_nflop, &
1945 unit_nr_dbcsr, mp2_env)
1946
1947 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1948 TYPE(dbt_type), INTENT(INOUT) :: t_r_occ, t_r_virt
1949 TYPE(im_time_force_type), INTENT(INOUT) :: force_data
1950 REAL(dp), INTENT(IN) :: fac
1951 INTEGER, INTENT(IN) :: cut_memory, n_mem_ri
1952 TYPE(dbt_type), INTENT(INOUT) :: t_kbkt, t_dm_occ, t_dm_virt, t_3c_o, t_3c_m, t_m_occ, &
1953 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, &
1954 t_3c_help_1, t_3c_help_2, t_3c_ints, t_3c_work
1955 INTEGER, DIMENSION(:), INTENT(IN) :: starts_array_mc, ends_array_mc, &
1956 batch_start_ri, batch_end_ri
1957 TYPE(hfx_compression_type), DIMENSION(:) :: t_3c_o_compressed
1958 TYPE(block_ind_type), DIMENSION(:), INTENT(INOUT) :: t_3c_o_ind
1959 LOGICAL, INTENT(IN) :: use_virial
1960 INTEGER, DIMENSION(:), INTENT(IN) :: atom_of_kind, kind_of
1961 REAL(dp), INTENT(IN) :: eps_filter
1962 REAL(dp), INTENT(INOUT) :: occ_ddint
1963 INTEGER(int_8), INTENT(INOUT) :: nze_ddint, dbcsr_nflop
1964 INTEGER, INTENT(IN) :: unit_nr_dbcsr
1965 TYPE(mp2_type) :: mp2_env
1966
1967 CHARACTER(LEN=*), PARAMETER :: routinen = 'perform_3c_ops'
1968
1969 INTEGER :: dummy_int, handle, handle2, i_mem, &
1970 i_xyz, j_mem, k_mem
1971 INTEGER(int_8) :: flop, nze
1972 INTEGER, DIMENSION(2, 1) :: ibounds, jbounds, kbounds
1973 INTEGER, DIMENSION(2, 2) :: bounds_2c
1974 INTEGER, DIMENSION(2, 3) :: bounds_cpy
1975 INTEGER, DIMENSION(3) :: bounds_3c
1976 REAL(dp) :: memory, occ, pref
1977 TYPE(block_ind_type), ALLOCATABLE, DIMENSION(:, :) :: blk_indices
1978 TYPE(hfx_compression_type), ALLOCATABLE, &
1979 DIMENSION(:, :) :: store_3c
1980
1981 CALL timeset(routinen, handle)
1982
1983 CALL dbt_get_info(t_3c_m, nfull_total=bounds_3c)
1984
1985 !Pre-compute and compress KBK^T * (pq|R)
1986 ALLOCATE (store_3c(n_mem_ri, cut_memory))
1987 ALLOCATE (blk_indices(n_mem_ri, cut_memory))
1988 memory = 0.0_dp
1989 CALL timeset(routinen//"_pre_3c", handle2)
1990 !temporarily build the full int 3c tensor
1991 CALL dbt_copy(t_3c_o, t_3c_0)
1992 DO i_mem = 1, cut_memory
1993 CALL decompress_tensor(t_3c_o, t_3c_o_ind(i_mem)%ind, t_3c_o_compressed(i_mem), &
1994 mp2_env%ri_rpa_im_time%eps_compress)
1995 CALL dbt_copy(t_3c_o, t_3c_ints)
1996 CALL dbt_copy(t_3c_o, t_3c_0, move_data=.true., summation=.true.)
1997
1998 DO k_mem = 1, n_mem_ri
1999 kbounds(:, 1) = [batch_start_ri(k_mem), batch_end_ri(k_mem)]
2000
2001 CALL alloc_containers(store_3c(k_mem, i_mem), 1)
2002
2003 !contract with KBK^T over the RI index and store
2004 CALL dbt_batched_contract_init(t_kbkt)
2005 CALL dbt_contract(1.0_dp, t_kbkt, t_3c_ints, 0.0_dp, t_3c_work, &
2006 contract_1=[2], notcontract_1=[1], &
2007 contract_2=[1], notcontract_2=[2, 3], &
2008 map_1=[1], map_2=[2, 3], filter_eps=eps_filter, &
2009 bounds_2=kbounds, flop=flop, unit_nr=unit_nr_dbcsr)
2010 CALL dbt_batched_contract_finalize(t_kbkt)
2011 dbcsr_nflop = dbcsr_nflop + flop
2012
2013 CALL dbt_copy(t_3c_work, t_3c_m, move_data=.true.)
2014 CALL compress_tensor(t_3c_m, blk_indices(k_mem, i_mem)%ind, store_3c(k_mem, i_mem), &
2015 mp2_env%ri_rpa_im_time%eps_compress, memory)
2016 END DO
2017 END DO !i_mem
2018 CALL dbt_clear(t_3c_m)
2019 CALL dbt_copy(t_3c_m, t_3c_ints)
2020 CALL timestop(handle2)
2021
2022 CALL dbt_batched_contract_init(t_r_occ)
2023 CALL dbt_batched_contract_init(t_r_virt)
2024 DO i_mem = 1, cut_memory
2025 ibounds(:, 1) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
2026
2027 !Compute the matrices M (integrals in t_3c_0)
2028 CALL timeset(routinen//"_3c_M", handle2)
2029 CALL dbt_batched_contract_init(t_dm_occ)
2030 CALL dbt_contract(1.0_dp, t_3c_0, t_dm_occ, 0.0_dp, t_3c_1, &
2031 contract_1=[3], notcontract_1=[1, 2], &
2032 contract_2=[1], notcontract_2=[2], &
2033 map_1=[1, 2], map_2=[3], filter_eps=eps_filter, &
2034 bounds_3=ibounds, flop=flop, unit_nr=unit_nr_dbcsr)
2035 dbcsr_nflop = dbcsr_nflop + flop
2036 CALL dbt_batched_contract_finalize(t_dm_occ)
2037 CALL dbt_copy(t_3c_1, t_m_occ, order=[1, 3, 2], move_data=.true.)
2038
2039 CALL dbt_batched_contract_init(t_dm_virt)
2040 CALL dbt_contract(1.0_dp, t_3c_0, t_dm_virt, 0.0_dp, t_3c_1, &
2041 contract_1=[3], notcontract_1=[1, 2], &
2042 contract_2=[1], notcontract_2=[2], &
2043 map_1=[1, 2], map_2=[3], filter_eps=eps_filter, &
2044 bounds_3=ibounds, flop=flop, unit_nr=unit_nr_dbcsr)
2045 dbcsr_nflop = dbcsr_nflop + flop
2046 CALL dbt_batched_contract_finalize(t_dm_virt)
2047 CALL dbt_copy(t_3c_1, t_m_virt, order=[1, 3, 2], move_data=.true.)
2048 CALL timestop(handle2)
2049
2050 !Compute the R matrices
2051 CALL timeset(routinen//"_3c_R", handle2)
2052 DO k_mem = 1, n_mem_ri
2053 CALL decompress_tensor(t_3c_m, blk_indices(k_mem, i_mem)%ind, store_3c(k_mem, i_mem), &
2054 mp2_env%ri_rpa_im_time%eps_compress)
2055 CALL dbt_copy(t_3c_m, t_3c_3, move_data=.true.)
2056
2057 CALL dbt_contract(1.0_dp, t_m_occ, t_3c_3, 1.0_dp, t_r_occ, &
2058 contract_1=[1, 2], notcontract_1=[3], &
2059 contract_2=[1, 2], notcontract_2=[3], &
2060 map_1=[1], map_2=[2], filter_eps=eps_filter, &
2061 flop=flop, unit_nr=unit_nr_dbcsr)
2062 dbcsr_nflop = dbcsr_nflop + flop
2063
2064 CALL dbt_contract(1.0_dp, t_m_virt, t_3c_3, 1.0_dp, t_r_virt, &
2065 contract_1=[1, 2], notcontract_1=[3], &
2066 contract_2=[1, 2], notcontract_2=[3], &
2067 map_1=[1], map_2=[2], filter_eps=eps_filter, &
2068 flop=flop, unit_nr=unit_nr_dbcsr)
2069 dbcsr_nflop = dbcsr_nflop + flop
2070 END DO
2071 CALL dbt_copy(t_3c_m, t_3c_3)
2072 CALL dbt_copy(t_3c_m, t_m_virt)
2073 CALL timestop(handle2)
2074
2075 CALL dbt_copy(t_m_occ, t_3c_4, move_data=.true.)
2076
2077 IF (cut_memory > 0) CALL dbt_batched_contract_init(t_kbkt)
2078 DO j_mem = 1, cut_memory
2079 jbounds(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
2080
2081 bounds_cpy(:, 1) = [1, bounds_3c(1)]
2082 bounds_cpy(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
2083 bounds_cpy(:, 3) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
2084 CALL dbt_copy(t_3c_sparse, t_3c_7, bounds=bounds_cpy)
2085
2086 CALL dbt_batched_contract_init(t_dm_virt)
2087 DO k_mem = 1, n_mem_ri
2088 bounds_2c(:, 1) = [batch_start_ri(k_mem), batch_end_ri(k_mem)]
2089 bounds_2c(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
2090
2091 CALL timeset(routinen//"_3c_dm", handle2)
2092
2093 !Calculate (mu nu| P) * D_occ * D_virt
2094 !Note: technically need M_occ*D_virt + M_virt*D_occ, but it is equivalent to 2*M_occ*D_virt
2095 CALL dbt_contract(2.0_dp, t_3c_4, t_dm_virt, 0.0_dp, t_3c_5, &
2096 contract_1=[3], notcontract_1=[1, 2], &
2097 contract_2=[1], notcontract_2=[2], &
2098 map_1=[1, 2], map_2=[3], filter_eps=eps_filter, &
2099 bounds_2=bounds_2c, bounds_3=jbounds, flop=flop, unit_nr=unit_nr_dbcsr)
2100 dbcsr_nflop = dbcsr_nflop + flop
2101
2102 CALL get_tensor_occupancy(t_3c_5, nze, occ)
2103 nze_ddint = nze_ddint + nze
2104 occ_ddint = occ_ddint + occ
2105
2106 ! Skip the expensive KBK^T contraction when the intermediate block is empty
2107 IF (nze == 0) THEN
2108 CALL dbt_clear(t_3c_5)
2109 cycle
2110 END IF
2111
2112 CALL dbt_copy(t_3c_5, t_3c_6, move_data=.true.)
2113 CALL timestop(handle2)
2114
2115 !Calculate the contraction of the above with K*B*K^T
2116 CALL timeset(routinen//"_3c_KBK", handle2)
2117 CALL dbt_contract(1.0_dp, t_kbkt, t_3c_6, 0.0_dp, t_3c_7, &
2118 contract_1=[2], notcontract_1=[1], &
2119 contract_2=[1], notcontract_2=[2, 3], &
2120 map_1=[1], map_2=[2, 3], &
2121 retain_sparsity=.true., flop=flop, unit_nr=unit_nr_dbcsr)
2122 dbcsr_nflop = dbcsr_nflop + flop
2123 CALL timestop(handle2)
2124 CALL dbt_copy(t_3c_7, t_3c_8, summation=.true.)
2125
2126 END DO !k_mem
2127 CALL dbt_batched_contract_finalize(t_dm_virt)
2128 END DO !j_mem
2129 IF (cut_memory > 0) CALL dbt_batched_contract_finalize(t_kbkt)
2130
2131 CALL dbt_copy(t_3c_8, t_3c_help_1, move_data=.true.)
2132
2133 pref = 1.0_dp*fac
2134 DO k_mem = 1, cut_memory
2135 DO i_xyz = 1, 3
2136 CALL dbt_clear(force_data%t_3c_der_RI(i_xyz))
2137 CALL decompress_tensor(force_data%t_3c_der_RI(i_xyz), force_data%t_3c_der_RI_ind(k_mem, i_xyz)%ind, &
2138 force_data%t_3c_der_RI_comp(k_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
2139 END DO
2140 CALL get_force_from_3c_trace(force, t_3c_help_1, force_data%t_3c_der_RI, atom_of_kind, kind_of, &
2141 force_data%idx_to_at_RI, pref, do_mp2=.true., deriv_dim=1)
2142 END DO
2143
2144 IF (use_virial) THEN
2145 CALL dbt_copy(t_3c_help_1, t_3c_help_2)
2146 CALL dbt_scale(t_3c_help_2, pref)
2147 CALL dbt_copy(t_3c_help_2, force_data%t_3c_virial_split, summation=.true., move_data=.true.)
2148 END IF
2149
2150 CALL dbt_copy(t_3c_help_1, t_3c_help_2)
2151 CALL dbt_copy(t_3c_help_1, t_3c_help_2, order=[1, 3, 2], move_data=.true., summation=.true.)
2152 DO k_mem = 1, cut_memory
2153 DO i_xyz = 1, 3
2154 CALL dbt_clear(force_data%t_3c_der_AO(i_xyz))
2155 CALL decompress_tensor(force_data%t_3c_der_AO(i_xyz), force_data%t_3c_der_AO_ind(k_mem, i_xyz)%ind, &
2156 force_data%t_3c_der_AO_comp(k_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
2157 END DO
2158 CALL get_force_from_3c_trace(force, t_3c_help_2, force_data%t_3c_der_AO, atom_of_kind, kind_of, &
2159 force_data%idx_to_at_AO, pref, do_mp2=.true., deriv_dim=3)
2160 END DO
2161
2162 CALL dbt_clear(t_3c_help_2)
2163 END DO !i_mem
2164 CALL dbt_batched_contract_finalize(t_r_occ)
2165 CALL dbt_batched_contract_finalize(t_r_virt)
2166
2167 DO k_mem = 1, n_mem_ri
2168 DO i_mem = 1, cut_memory
2169 CALL dealloc_containers(store_3c(k_mem, i_mem), dummy_int)
2170 END DO
2171 END DO
2172 DEALLOCATE (store_3c, blk_indices)
2173
2174 CALL timestop(handle)
2175
2176 END SUBROUTINE perform_3c_ops
2177
2178! **************************************************************************************************
2179!> \brief All the forces that can be calculated after the loop on the Laplace quaradture, using
2180!> terms collected during the said loop. This inludes the z-vector equation and its reponse
2181!> forces, as well as the force coming from the trace with the derivative of the KS matrix
2182!> \param force_data ...
2183!> \param unit_nr ...
2184!> \param qs_env ...
2185! **************************************************************************************************
2186 SUBROUTINE calc_post_loop_forces(force_data, unit_nr, qs_env)
2187
2188 TYPE(im_time_force_type), INTENT(INOUT) :: force_data
2189 INTEGER, INTENT(IN) :: unit_nr
2190 TYPE(qs_environment_type), POINTER :: qs_env
2191
2192 CHARACTER(len=*), PARAMETER :: routinen = 'calc_post_loop_forces'
2193
2194 INTEGER :: handle, ispin, nao, nao_aux, nocc, nspins
2195 LOGICAL :: do_exx
2196 REAL(dp) :: focc
2197 TYPE(admm_type), POINTER :: admm_env
2198 TYPE(cp_fm_struct_type), POINTER :: fm_struct
2199 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: cpmos, mo_occ
2200 TYPE(cp_fm_type), POINTER :: mo_coeff
2201 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_p_work, matrix_p_mp2, &
2202 matrix_p_mp2_admm, matrix_s, &
2203 matrix_s_aux, work_admm, yp_admm
2204 TYPE(dft_control_type), POINTER :: dft_control
2205 TYPE(linres_control_type), POINTER :: linres_control
2206 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2207 TYPE(qs_p_env_type), POINTER :: p_env
2208 TYPE(section_vals_type), POINTER :: hfx_section, lr_section
2209
2210 NULLIFY (linres_control, p_env, dft_control, matrix_s, mos, mo_coeff, fm_struct, lr_section, &
2211 dbcsr_p_work, yp_admm, matrix_p_mp2, admm_env, work_admm, matrix_s_aux, matrix_p_mp2_admm)
2212
2213 CALL timeset(routinen, handle)
2214
2215 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_s=matrix_s, mos=mos)
2216 nspins = dft_control%nspins
2217
2218 ! Setting up for the z-vector equation
2219
2220 ! Initialize linres_control
2221 lr_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%LOW_SCALING%CPHF")
2222
2223 ALLOCATE (linres_control)
2224 CALL section_vals_val_get(lr_section, "MAX_ITER", i_val=linres_control%max_iter)
2225 CALL section_vals_val_get(lr_section, "EPS_CONV", r_val=linres_control%eps)
2226 CALL section_vals_val_get(lr_section, "PRECONDITIONER", i_val=linres_control%preconditioner_type)
2227 CALL section_vals_val_get(lr_section, "ENERGY_GAP", r_val=linres_control%energy_gap)
2228
2229 linres_control%do_kernel = .true.
2230 linres_control%lr_triplet = .false.
2231 linres_control%converged = .false.
2232 linres_control%eps_filter = qs_env%mp2_env%ri_rpa_im_time%eps_filter
2233
2234 CALL set_qs_env(qs_env, linres_control=linres_control)
2235
2236 IF (unit_nr > 0) THEN
2237 WRITE (unit_nr, *)
2238 WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations'
2239 WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Convergence threshold:', linres_control%eps
2240 WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Maximum number of iterations: ', linres_control%max_iter
2241 END IF
2242
2243 ALLOCATE (p_env)
2244 CALL p_env_create(p_env, qs_env, orthogonal_orbitals=.true., linres_control=linres_control)
2245 CALL p_env_psi0_changed(p_env, qs_env)
2246
2247 ! Matrix allocation
2248 CALL dbcsr_allocate_matrix_set(p_env%p1, nspins)
2249 CALL dbcsr_allocate_matrix_set(p_env%w1, nspins)
2250 CALL dbcsr_allocate_matrix_set(dbcsr_p_work, nspins)
2251 DO ispin = 1, nspins
2252 ALLOCATE (p_env%p1(ispin)%matrix, p_env%w1(ispin)%matrix, dbcsr_p_work(ispin)%matrix)
2253 CALL dbcsr_create(matrix=p_env%p1(ispin)%matrix, template=matrix_s(1)%matrix)
2254 CALL dbcsr_create(matrix=p_env%w1(ispin)%matrix, template=matrix_s(1)%matrix)
2255 CALL dbcsr_create(matrix=dbcsr_p_work(ispin)%matrix, template=matrix_s(1)%matrix)
2256 CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_s(1)%matrix)
2257 CALL dbcsr_copy(p_env%w1(ispin)%matrix, matrix_s(1)%matrix)
2258 CALL dbcsr_copy(dbcsr_p_work(ispin)%matrix, matrix_s(1)%matrix)
2259 CALL dbcsr_set(p_env%p1(ispin)%matrix, 0.0_dp)
2260 CALL dbcsr_set(p_env%w1(ispin)%matrix, 0.0_dp)
2261 CALL dbcsr_set(dbcsr_p_work(ispin)%matrix, 0.0_dp)
2262 END DO
2263
2264 IF (dft_control%do_admm) THEN
2265 CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux)
2266 CALL dbcsr_allocate_matrix_set(p_env%p1_admm, nspins)
2267 CALL dbcsr_allocate_matrix_set(work_admm, nspins)
2268 DO ispin = 1, nspins
2269 ALLOCATE (p_env%p1_admm(ispin)%matrix, work_admm(ispin)%matrix)
2270 CALL dbcsr_create(p_env%p1_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2271 CALL dbcsr_copy(p_env%p1_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2272 CALL dbcsr_set(p_env%p1_admm(ispin)%matrix, 0.0_dp)
2273 CALL dbcsr_create(work_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2274 CALL dbcsr_copy(work_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2275 CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp)
2276 END DO
2277 END IF
2278
2279 ! Preparing the RHS of the z-vector equation
2280 CALL prepare_for_response(force_data, qs_env)
2281 ALLOCATE (cpmos(nspins), mo_occ(nspins))
2282 DO ispin = 1, nspins
2283 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, homo=nocc)
2284 NULLIFY (fm_struct)
2285 CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, &
2286 template_fmstruct=mo_coeff%matrix_struct)
2287 CALL cp_fm_create(cpmos(ispin), fm_struct)
2288 CALL cp_fm_set_all(cpmos(ispin), 0.0_dp)
2289 CALL cp_fm_create(mo_occ(ispin), fm_struct)
2290 CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc)
2291 CALL cp_fm_struct_release(fm_struct)
2292 END DO
2293
2294 ! in case of EXX, need to add the HF Hamiltonian to the RHS of the Z-vector equation
2295 ! Strategy: we take the ks_matrix, remove the current xc contribution, and then add the RPA HF one
2296 do_exx = .false.
2297 IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
2298 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
2299 CALL section_vals_get(hfx_section, explicit=do_exx)
2300 END IF
2301
2302 IF (do_exx) THEN
2303 CALL add_exx_to_rhs(rhs=force_data%sum_O_tau, &
2304 qs_env=qs_env, &
2305 ext_hfx_section=hfx_section, &
2306 x_data=qs_env%mp2_env%ri_rpa%x_data, &
2307 recalc_integrals=.false., &
2308 do_admm=qs_env%mp2_env%ri_rpa%do_admm, &
2309 do_exx=do_exx, &
2310 reuse_hfx=qs_env%mp2_env%ri_rpa%reuse_hfx)
2311 END IF
2312
2313 focc = 2.0_dp
2314 IF (nspins == 1) focc = 4.0_dp
2315 DO ispin = 1, nspins
2316 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
2317 CALL cp_dbcsr_sm_fm_multiply(force_data%sum_O_tau(ispin)%matrix, mo_occ(ispin), &
2318 cpmos(ispin), nocc, &
2319 alpha=focc, beta=0.0_dp)
2320 END DO
2321
2322 ! The z-vector equation and associated forces
2323 CALL response_equation_new(qs_env, p_env, cpmos, unit_nr)
2324
2325 ! Save the mp2 density matrix
2326 CALL get_qs_env(qs_env, matrix_p_mp2=matrix_p_mp2)
2327 IF (ASSOCIATED(matrix_p_mp2)) CALL dbcsr_deallocate_matrix_set(matrix_p_mp2)
2328 DO ispin = 1, nspins
2329 CALL dbcsr_copy(dbcsr_p_work(ispin)%matrix, p_env%p1(ispin)%matrix)
2330 CALL dbcsr_add(dbcsr_p_work(ispin)%matrix, force_data%sum_YP_tau(ispin)%matrix, 1.0_dp, 1.0_dp)
2331 END DO
2332 CALL set_ks_env(qs_env%ks_env, matrix_p_mp2=dbcsr_p_work)
2333
2334 IF (dft_control%do_admm) THEN
2335 CALL dbcsr_allocate_matrix_set(yp_admm, nspins)
2336 CALL get_qs_env(qs_env, matrix_p_mp2_admm=matrix_p_mp2_admm, admm_env=admm_env)
2337 nao = admm_env%nao_orb
2338 nao_aux = admm_env%nao_aux_fit
2339 IF (ASSOCIATED(matrix_p_mp2_admm)) CALL dbcsr_deallocate_matrix_set(matrix_p_mp2_admm)
2340 DO ispin = 1, nspins
2341
2342 !sum_YP_tau in the auxiliary basis
2343 CALL copy_dbcsr_to_fm(force_data%sum_YP_tau(ispin)%matrix, admm_env%work_orb_orb)
2344 CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, admm_env%work_orb_orb, &
2345 0.0_dp, admm_env%work_aux_orb)
2346 CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
2347 0.0_dp, admm_env%work_aux_aux)
2348 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_admm(ispin)%matrix, keep_sparsity=.true.)
2349
2350 !save the admm representation od sum_YP_tau
2351 ALLOCATE (yp_admm(ispin)%matrix)
2352 CALL dbcsr_create(yp_admm(ispin)%matrix, template=work_admm(ispin)%matrix)
2353 CALL dbcsr_copy(yp_admm(ispin)%matrix, work_admm(ispin)%matrix)
2354
2355 CALL dbcsr_add(work_admm(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
2356
2357 END DO
2358 CALL set_ks_env(qs_env%ks_env, matrix_p_mp2_admm=work_admm)
2359 END IF
2360
2361 !Calculate the response force and the force from the trace with F
2362 CALL update_im_time_forces(p_env, force_data%sum_O_tau, force_data%sum_YP_tau, yp_admm, qs_env)
2363
2364 !clean-up
2365 IF (dft_control%do_admm) CALL dbcsr_deallocate_matrix_set(yp_admm)
2366
2367 CALL cp_fm_release(cpmos)
2368 CALL cp_fm_release(mo_occ)
2369 CALL p_env_release(p_env)
2370 DEALLOCATE (p_env)
2371
2372 CALL timestop(handle)
2373
2374 END SUBROUTINE calc_post_loop_forces
2375
2376! **************************************************************************************************
2377!> \brief Prepares the RHS of the z-vector equation. Apply the xc and HFX kernel on the previously
2378!> stored sum_YP_tau density, and add it to the final force_data%sum_O_tau quantity
2379!> \param force_data ...
2380!> \param qs_env ...
2381! **************************************************************************************************
2382 SUBROUTINE prepare_for_response(force_data, qs_env)
2383
2384 TYPE(im_time_force_type), INTENT(INOUT) :: force_data
2385 TYPE(qs_environment_type), POINTER :: qs_env
2386
2387 CHARACTER(LEN=*), PARAMETER :: routinen = 'prepare_for_response'
2388
2389 INTEGER :: handle, ispin, nao, nao_aux, nspins
2390 LOGICAL :: do_hfx, do_tau, do_tau_admm
2391 REAL(dp) :: ehartree
2392 TYPE(admm_type), POINTER :: admm_env
2393 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_p_work, ker_tau_admm, matrix_s, &
2394 matrix_s_aux, work_admm
2395 TYPE(dbcsr_type) :: dbcsr_work
2396 TYPE(dft_control_type), POINTER :: dft_control
2397 TYPE(pw_c1d_gs_type) :: rhoz_tot_gspace, zv_hartree_gspace
2398 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rhoz_g
2399 TYPE(pw_env_type), POINTER :: pw_env
2400 TYPE(pw_poisson_type), POINTER :: poisson_env
2401 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2402 TYPE(pw_r3d_rs_type) :: zv_hartree_rspace
2403 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rhoz_r, tauz_r, v_xc, v_xc_tau
2404 TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit, rhoz
2405 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set
2406 TYPE(section_vals_type), POINTER :: hfx_section, xc_section
2407 TYPE(task_list_type), POINTER :: task_list_aux_fit
2408
2409 NULLIFY (pw_env, rhoz_r, rhoz_g, tauz_r, v_xc, v_xc_tau, &
2410 poisson_env, auxbas_pw_pool, dft_control, admm_env, xc_section, rho, rho_aux_fit, &
2411 task_list_aux_fit, ker_tau_admm, work_admm, dbcsr_p_work, matrix_s, hfx_section)
2412 NULLIFY (rho0_atom_set, rho1_atom_set)
2413
2414 CALL timeset(routinen, handle)
2415
2416 CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env, rho=rho, matrix_s=matrix_s)
2417 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
2418 nspins = dft_control%nspins
2419
2420 CALL dbcsr_allocate_matrix_set(dbcsr_p_work, nspins)
2421 DO ispin = 1, nspins
2422 ALLOCATE (dbcsr_p_work(ispin)%matrix)
2423 CALL dbcsr_create(matrix=dbcsr_p_work(ispin)%matrix, template=matrix_s(1)%matrix)
2424 CALL dbcsr_copy(dbcsr_p_work(ispin)%matrix, matrix_s(1)%matrix)
2425 CALL dbcsr_set(dbcsr_p_work(ispin)%matrix, 0.0_dp)
2426 END DO
2427
2428 !Apply the kernel on the density saved in force_data%sum_YP_tau
2429 ALLOCATE (rhoz_r(nspins), rhoz_g(nspins))
2430 DO ispin = 1, nspins
2431 CALL auxbas_pw_pool%create_pw(rhoz_r(ispin))
2432 CALL auxbas_pw_pool%create_pw(rhoz_g(ispin))
2433 END DO
2434 CALL auxbas_pw_pool%create_pw(rhoz_tot_gspace)
2435 CALL auxbas_pw_pool%create_pw(zv_hartree_rspace)
2436 CALL auxbas_pw_pool%create_pw(zv_hartree_gspace)
2437
2438 CALL pw_zero(rhoz_tot_gspace)
2439 DO ispin = 1, nspins
2440 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=force_data%sum_YP_tau(ispin)%matrix, &
2441 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin))
2442 CALL pw_axpy(rhoz_g(ispin), rhoz_tot_gspace)
2443 END DO
2444
2445 CALL pw_poisson_solve(poisson_env, rhoz_tot_gspace, ehartree, &
2446 zv_hartree_gspace)
2447
2448 CALL pw_transfer(zv_hartree_gspace, zv_hartree_rspace)
2449 CALL pw_scale(zv_hartree_rspace, zv_hartree_rspace%pw_grid%dvol)
2450
2451 CALL qs_rho_get(rho, tau_r_valid=do_tau)
2452 IF (do_tau) THEN
2453 block
2454 TYPE(pw_c1d_gs_type) :: tauz_g
2455 ALLOCATE (tauz_r(nspins))
2456 CALL auxbas_pw_pool%create_pw(tauz_g)
2457 DO ispin = 1, nspins
2458 CALL auxbas_pw_pool%create_pw(tauz_r(ispin))
2459
2460 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=force_data%sum_YP_tau(ispin)%matrix, &
2461 rho=tauz_r(ispin), rho_gspace=tauz_g, compute_tau=.true.)
2462 END DO
2463 CALL auxbas_pw_pool%give_back_pw(tauz_g)
2464 END block
2465 END IF
2466
2467 IF (dft_control%do_admm) THEN
2468 CALL get_qs_env(qs_env, admm_env=admm_env)
2469 xc_section => admm_env%xc_section_primary
2470 ELSE
2471 xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
2472 END IF
2473
2474 !Primary XC kernel
2475 ALLOCATE (rhoz)
2476 CALL qs_rho_create(rhoz)
2477 IF (ASSOCIATED(rhoz_r)) THEN
2478 CALL qs_rho_set(rhoz, rho_r=rhoz_r, rho_r_valid=.true.)
2479 END IF
2480 IF (ASSOCIATED(rhoz_g)) THEN
2481 CALL qs_rho_set(rhoz, rho_g=rhoz_g, rho_g_valid=.true.)
2482 END IF
2483 IF (ASSOCIATED(tauz_r)) THEN
2484 CALL qs_rho_set(rhoz, tau_r=tauz_r, tau_r_valid=.true.)
2485 END IF
2486 !
2487 CALL qs_fxc_create(qs_env, rho, rhoz, rho0_atom_set, xc_section, .false., &
2488 v_xc, v_xc_tau, rho1_atom_set, &
2489 dispersion_env=qs_env%dispersion_env)
2490 !
2491 DEALLOCATE (rhoz)
2492
2493 DO ispin = 1, nspins
2494 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2495 CALL pw_axpy(zv_hartree_rspace, v_xc(ispin))
2496 CALL integrate_v_rspace(qs_env=qs_env, &
2497 v_rspace=v_xc(ispin), &
2498 hmat=dbcsr_p_work(ispin), &
2499 calculate_forces=.false.)
2500
2501 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2502 END DO
2503 CALL auxbas_pw_pool%give_back_pw(rhoz_tot_gspace)
2504 CALL auxbas_pw_pool%give_back_pw(zv_hartree_rspace)
2505 CALL auxbas_pw_pool%give_back_pw(zv_hartree_gspace)
2506 DEALLOCATE (v_xc)
2507
2508 IF (do_tau) THEN
2509 DO ispin = 1, nspins
2510 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2511 CALL integrate_v_rspace(qs_env=qs_env, &
2512 v_rspace=v_xc_tau(ispin), &
2513 hmat=dbcsr_p_work(ispin), &
2514 compute_tau=.true., &
2515 calculate_forces=.false.)
2516 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2517 END DO
2518 DEALLOCATE (v_xc_tau)
2519 END IF
2520
2521 !Auxiliary xc kernel (admm)
2522 IF (dft_control%do_admm) THEN
2523 CALL get_qs_env(qs_env, admm_env=admm_env)
2524 CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux, &
2525 task_list_aux_fit=task_list_aux_fit, rho_aux_fit=rho_aux_fit)
2526
2527 CALL qs_rho_get(rho_aux_fit, tau_r_valid=do_tau_admm)
2528
2529 CALL dbcsr_allocate_matrix_set(work_admm, nspins)
2530 CALL dbcsr_allocate_matrix_set(ker_tau_admm, nspins)
2531 DO ispin = 1, nspins
2532 ALLOCATE (work_admm(ispin)%matrix, ker_tau_admm(ispin)%matrix)
2533 CALL dbcsr_create(work_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2534 CALL dbcsr_copy(work_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2535 CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp)
2536 CALL dbcsr_create(ker_tau_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2537 CALL dbcsr_copy(ker_tau_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2538 CALL dbcsr_set(ker_tau_admm(ispin)%matrix, 0.0_dp)
2539 END DO
2540
2541 !get the density in the auxuliary density
2542 cpassert(ASSOCIATED(admm_env%work_orb_orb))
2543 cpassert(ASSOCIATED(admm_env%work_aux_orb))
2544 cpassert(ASSOCIATED(admm_env%work_aux_aux))
2545 nao = admm_env%nao_orb
2546 nao_aux = admm_env%nao_aux_fit
2547 DO ispin = 1, nspins
2548 CALL copy_dbcsr_to_fm(force_data%sum_YP_tau(ispin)%matrix, admm_env%work_orb_orb)
2549 CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, admm_env%work_orb_orb, &
2550 0.0_dp, admm_env%work_aux_orb)
2551 CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
2552 0.0_dp, admm_env%work_aux_aux)
2553 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, ker_tau_admm(ispin)%matrix, keep_sparsity=.true.)
2554 END DO
2555
2556 IF (.NOT. qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
2557 DO ispin = 1, nspins
2558 CALL pw_zero(rhoz_r(ispin))
2559 CALL pw_zero(rhoz_g(ispin))
2560 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=ker_tau_admm(ispin)%matrix, &
2561 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin), &
2562 basis_type="AUX_FIT", task_list_external=task_list_aux_fit)
2563 END DO
2564
2565 IF (do_tau_admm) THEN
2566 block
2567 TYPE(pw_c1d_gs_type) :: tauz_g
2568 CALL auxbas_pw_pool%create_pw(tauz_g)
2569 DO ispin = 1, nspins
2570 CALL pw_zero(tauz_r(ispin))
2571 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=ker_tau_admm(ispin)%matrix, &
2572 rho=tauz_r(ispin), rho_gspace=tauz_g, &
2573 basis_type="AUX_FIT", task_list_external=task_list_aux_fit, &
2574 compute_tau=.true.)
2575 END DO
2576 CALL auxbas_pw_pool%give_back_pw(tauz_g)
2577 END block
2578 END IF
2579
2580 xc_section => admm_env%xc_section_aux
2581 ALLOCATE (rhoz)
2582 CALL qs_rho_create(rhoz)
2583 IF (ASSOCIATED(rhoz_r)) THEN
2584 CALL qs_rho_set(rhoz, rho_r=rhoz_r, rho_r_valid=.true.)
2585 END IF
2586 IF (ASSOCIATED(rhoz_g)) THEN
2587 CALL qs_rho_set(rhoz, rho_g=rhoz_g, rho_g_valid=.true.)
2588 END IF
2589 IF (ASSOCIATED(tauz_r)) THEN
2590 CALL qs_rho_set(rhoz, tau_r=tauz_r, tau_r_valid=.true.)
2591 END IF
2592 !
2593 CALL qs_fxc_create(qs_env, rho_aux_fit, rhoz, rho0_atom_set, xc_section, .false., &
2594 v_xc, v_xc_tau, rho1_atom_set)
2595 !
2596 DEALLOCATE (rhoz)
2597
2598 DO ispin = 1, nspins
2599 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2600 CALL integrate_v_rspace(qs_env=qs_env, &
2601 v_rspace=v_xc(ispin), &
2602 hmat=work_admm(ispin), &
2603 calculate_forces=.false., &
2604 basis_type="AUX_FIT", &
2605 task_list_external=task_list_aux_fit)
2606 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2607 END DO
2608 DEALLOCATE (v_xc)
2609
2610 IF (do_tau_admm) THEN
2611 DO ispin = 1, nspins
2612 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2613 CALL integrate_v_rspace(qs_env=qs_env, &
2614 v_rspace=v_xc_tau(ispin), &
2615 hmat=work_admm(ispin), &
2616 calculate_forces=.false., &
2617 basis_type="AUX_FIT", &
2618 task_list_external=task_list_aux_fit, &
2619 compute_tau=.true.)
2620 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2621 END DO
2622 DEALLOCATE (v_xc_tau)
2623 END IF
2624 END IF !admm
2625 END IF
2626
2627 DO ispin = 1, nspins
2628 CALL auxbas_pw_pool%give_back_pw(rhoz_r(ispin))
2629 CALL auxbas_pw_pool%give_back_pw(rhoz_g(ispin))
2630 END DO
2631 DEALLOCATE (rhoz_r, rhoz_g)
2632
2633 IF (do_tau) THEN
2634 DO ispin = 1, nspins
2635 CALL auxbas_pw_pool%give_back_pw(tauz_r(ispin))
2636 END DO
2637 DEALLOCATE (tauz_r)
2638 END IF
2639
2640 !HFX kernel
2641 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
2642 CALL section_vals_get(hfx_section, explicit=do_hfx)
2643 IF (do_hfx) THEN
2644 IF (dft_control%do_admm) THEN
2645 CALL tddft_hfx_matrix(work_admm, ker_tau_admm, qs_env, .false., .false.)
2646
2647 !Going back to primary basis
2648 CALL dbcsr_create(dbcsr_work, template=dbcsr_p_work(1)%matrix)
2649 CALL dbcsr_copy(dbcsr_work, dbcsr_p_work(1)%matrix)
2650 CALL dbcsr_set(dbcsr_work, 0.0_dp)
2651 DO ispin = 1, nspins
2652 CALL copy_dbcsr_to_fm(work_admm(ispin)%matrix, admm_env%work_aux_aux)
2653 CALL parallel_gemm('N', 'N', nao_aux, nao, nao_aux, 1.0_dp, admm_env%work_aux_aux, admm_env%A, &
2654 0.0_dp, admm_env%work_aux_orb)
2655 CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, admm_env%work_aux_orb, &
2656 0.0_dp, admm_env%work_orb_orb)
2657 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_work, keep_sparsity=.true.)
2658 CALL dbcsr_add(dbcsr_p_work(ispin)%matrix, dbcsr_work, 1.0_dp, 1.0_dp)
2659 END DO
2660 CALL dbcsr_release(dbcsr_work)
2661 CALL dbcsr_deallocate_matrix_set(ker_tau_admm)
2662 ELSE
2663 CALL tddft_hfx_matrix(dbcsr_p_work, force_data%sum_YP_tau, qs_env, .false., .false.)
2664 END IF
2665 END IF
2666
2667 DO ispin = 1, nspins
2668 CALL dbcsr_add(force_data%sum_O_tau(ispin)%matrix, dbcsr_p_work(ispin)%matrix, 1.0_dp, 1.0_dp)
2669 END DO
2670
2671 CALL dbcsr_deallocate_matrix_set(dbcsr_p_work)
2672 CALL dbcsr_deallocate_matrix_set(work_admm)
2673
2674 CALL timestop(handle)
2675
2676 END SUBROUTINE prepare_for_response
2677
2678! **************************************************************************************************
2679!> \brief Calculate the force and virial due to the (P|Q) GPW integral derivatives
2680!> \param G_PQ ...
2681!> \param force ...
2682!> \param h_stress ...
2683!> \param use_virial ...
2684!> \param mp2_env ...
2685!> \param qs_env ...
2686! **************************************************************************************************
2687 SUBROUTINE get_2c_gpw_forces(G_PQ, force, h_stress, use_virial, mp2_env, qs_env)
2688
2689 TYPE(dbcsr_type), INTENT(INOUT) :: g_pq
2690 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2691 REAL(dp), DIMENSION(3, 3), INTENT(INOUT) :: h_stress
2692 LOGICAL, INTENT(IN) :: use_virial
2693 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
2694 TYPE(qs_environment_type), POINTER :: qs_env
2695
2696 CHARACTER(len=*), PARAMETER :: routinen = 'get_2c_gpw_forces'
2697
2698 INTEGER :: atom_a, color, handle, i, i_ri, i_xyz, iatom, igrid_level, ikind, ipgf, iset, j, &
2699 j_ri, jatom, lb_ri, n_ri, natom, ncoa, ncoms, nkind, nproc, nseta, o1, offset, pdims(2), &
2700 sgfa, ub_ri
2701 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, iproc_map, kind_of, &
2702 sizes_ri
2703 INTEGER, DIMENSION(:), POINTER :: col_dist, la_max, la_min, npgfa, nsgfa, &
2704 row_dist
2705 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, pgrid
2706 LOGICAL :: found, one_proc_group
2707 REAL(dp) :: cutoff_old, radius, relative_cutoff_old
2708 REAL(dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old, wf_vector
2709 REAL(dp), DIMENSION(3) :: force_a, force_b, ra
2710 REAL(dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
2711 REAL(kind=dp), DIMENSION(:, :), POINTER :: h_tmp, i_ab, pab, pblock, sphi_a, zeta
2712 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2713 TYPE(cell_type), POINTER :: cell
2714 TYPE(dbcsr_distribution_type) :: dbcsr_dist
2715 TYPE(dbcsr_type) :: tmp_g_pq
2716 TYPE(dft_control_type), POINTER :: dft_control
2717 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
2718 DIMENSION(:), TARGET :: basis_set_ri_aux
2719 TYPE(gto_basis_set_type), POINTER :: basis_set_a
2720 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2721 TYPE(mp_para_env_type), POINTER :: para_env, para_env_ext
2722 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2723 POINTER :: sab_orb
2724 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2725 TYPE(pw_c1d_gs_type) :: dvg(3), pot_g, rho_g, rho_g_copy
2726 TYPE(pw_env_type), POINTER :: pw_env_ext
2727 TYPE(pw_poisson_type), POINTER :: poisson_env
2728 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2729 TYPE(pw_r3d_rs_type) :: psi_l, rho_r
2730 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2731 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
2732 TYPE(task_list_type), POINTER :: task_list_ext
2733
2734 NULLIFY (sab_orb, task_list_ext, particle_set, qs_kind_set, dft_control, pw_env_ext, auxbas_pw_pool, &
2735 poisson_env, atomic_kind_set, para_env, cell, rs_v, mos, basis_set_a)
2736
2737 CALL timeset(routinen, handle)
2738
2739 CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env, sab_orb=sab_orb, &
2740 natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, particle_set=particle_set, &
2741 mos=mos, cell=cell, atomic_kind_set=atomic_kind_set)
2742
2743 !The idea is to use GPW to compute the integrals and derivatives. Because the potential needs
2744 !to be calculated for each phi_j (column) of all AO pairs, and because that is expensive, we want
2745 !to minimize the amount of time we do that. Therefore, we work with a special distribution, where
2746 !each column of the resulting DBCSR matrix is mapped to a sub-communicator.
2747
2748 !Try to get the optimal pdims (we want a grid that is flat: many cols, few rows)
2749 IF (para_env%num_pe <= natom) THEN
2750 pdims(1) = 1
2751 pdims(2) = para_env%num_pe
2752 ELSE
2753 DO i = natom, 1, -1
2754 IF (modulo(para_env%num_pe, i) == 0) THEN
2755 pdims(1) = para_env%num_pe/i
2756 pdims(2) = i
2757 EXIT
2758 END IF
2759 END DO
2760 END IF
2761
2762 ALLOCATE (row_dist(natom), col_dist(natom))
2763 DO iatom = 1, natom
2764 row_dist(iatom) = modulo(iatom, pdims(1))
2765 END DO
2766 DO jatom = 1, natom
2767 col_dist(jatom) = modulo(jatom, pdims(2))
2768 END DO
2769
2770 ALLOCATE (pgrid(0:pdims(1) - 1, 0:pdims(2) - 1))
2771 nproc = 0
2772 DO i = 0, pdims(1) - 1
2773 DO j = 0, pdims(2) - 1
2774 pgrid(i, j) = nproc
2775 nproc = nproc + 1
2776 END DO
2777 END DO
2778
2779 CALL dbcsr_distribution_new(dbcsr_dist, group=para_env%get_handle(), pgrid=pgrid, row_dist=row_dist, col_dist=col_dist)
2780
2781 !The temporary DBCSR integrals and derivatives matrices in this flat distribution
2782 CALL dbcsr_create(tmp_g_pq, template=g_pq, matrix_type=dbcsr_type_no_symmetry, dist=dbcsr_dist)
2783 CALL dbcsr_complete_redistribute(g_pq, tmp_g_pq)
2784
2785 ALLOCATE (basis_set_ri_aux(nkind), sizes_ri(natom))
2786 CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
2787 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ri, basis=basis_set_ri_aux)
2788 n_ri = sum(sizes_ri)
2789
2790 one_proc_group = mp2_env%mp2_num_proc == 1
2791 ALLOCATE (para_env_ext)
2792 IF (one_proc_group) THEN
2793 !one subgroup per proc
2794 CALL para_env_ext%from_split(para_env, para_env%mepos)
2795 ELSE
2796 !Split the communicator accross the columns of the matrix
2797 ncoms = min(pdims(2), para_env%num_pe/mp2_env%mp2_num_proc)
2798 DO i = 0, pdims(1) - 1
2799 DO j = 0, pdims(2) - 1
2800 IF (pgrid(i, j) == para_env%mepos) color = modulo(j + 1, ncoms)
2801 END DO
2802 END DO
2803 CALL para_env_ext%from_split(para_env, color)
2804 END IF
2805
2806 !sab_orb and task_list_ext are essentially dummies
2807 CALL prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_ext, pw_env_ext, &
2808 auxbas_pw_pool, poisson_env, task_list_ext, rho_r, rho_g, pot_g, psi_l, sab_orb)
2809
2810 IF (use_virial) THEN
2811 CALL auxbas_pw_pool%create_pw(rho_g_copy)
2812 DO i_xyz = 1, 3
2813 CALL auxbas_pw_pool%create_pw(dvg(i_xyz))
2814 END DO
2815 END IF
2816
2817 ALLOCATE (wf_vector(n_ri))
2818
2819 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
2820
2821 ALLOCATE (iproc_map(natom))
2822
2823 !Loop over the atomic blocks
2824 DO jatom = 1, natom
2825
2826 !Only calculate if on the correct sub-communicator/proc
2827 IF (one_proc_group) THEN
2828 iproc_map = 0
2829 DO iatom = 1, natom
2830 IF (pgrid(row_dist(iatom), col_dist(jatom)) == para_env%mepos) iproc_map(iatom) = 1
2831 END DO
2832 IF (.NOT. any(iproc_map == 1)) cycle
2833 ELSE
2834 IF (.NOT. modulo(col_dist(jatom) + 1, ncoms) == color) cycle
2835 END IF
2836
2837 lb_ri = sum(sizes_ri(1:jatom - 1))
2838 ub_ri = lb_ri + sizes_ri(jatom)
2839 DO j_ri = lb_ri + 1, ub_ri
2840
2841 wf_vector = 0.0_dp
2842 wf_vector(j_ri) = 1.0_dp
2843
2844 CALL collocate_function(wf_vector, psi_l, rho_g, atomic_kind_set, qs_kind_set, cell, &
2845 particle_set, pw_env_ext, dft_control%qs_control%eps_rho_rspace, &
2846 basis_type="RI_AUX")
2847
2848 IF (use_virial) THEN
2849 CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, mp2_env%potential_parameter, dvg)
2850
2851 wf_vector = 0.0_dp
2852 DO iatom = 1, natom
2853 !only compute if i,j atom pair on correct proc
2854 IF (one_proc_group) THEN
2855 IF (.NOT. iproc_map(iatom) == 1) cycle
2856 END IF
2857
2858 CALL dbcsr_get_block_p(tmp_g_pq, iatom, jatom, pblock, found)
2859 IF (.NOT. found) cycle
2860
2861 i_ri = sum(sizes_ri(1:iatom - 1))
2862 wf_vector(i_ri + 1:i_ri + sizes_ri(iatom)) = pblock(:, j_ri - lb_ri)
2863 END DO
2864
2865 CALL pw_copy(rho_g, rho_g_copy)
2866 CALL collocate_function(wf_vector, psi_l, rho_g, atomic_kind_set, qs_kind_set, cell, &
2867 particle_set, pw_env_ext, dft_control%qs_control%eps_rho_rspace, &
2868 basis_type="RI_AUX")
2869
2870 CALL calc_potential_gpw(psi_l, rho_g, poisson_env, pot_g, mp2_env%potential_parameter, &
2871 no_transfer=.true.)
2872 CALL virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, &
2873 mp2_env%potential_parameter, para_env_ext)
2874 ELSE
2875 CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, mp2_env%potential_parameter)
2876 END IF
2877
2878 NULLIFY (rs_v)
2879 CALL pw_env_get(pw_env_ext, rs_grids=rs_v)
2880 CALL potential_pw2rs(rs_v, rho_r, pw_env_ext)
2881
2882 DO iatom = 1, natom
2883
2884 !only compute if i,j atom pair on correct proc
2885 IF (one_proc_group) THEN
2886 IF (.NOT. iproc_map(iatom) == 1) cycle
2887 END IF
2888
2889 force_a(:) = 0.0_dp
2890 force_b(:) = 0.0_dp
2891 IF (use_virial) THEN
2892 my_virial_a = 0.0_dp
2893 my_virial_b = 0.0_dp
2894 END IF
2895
2896 ikind = kind_of(iatom)
2897 atom_a = atom_of_kind(iatom)
2898
2899 basis_set_a => basis_set_ri_aux(ikind)%gto_basis_set
2900 first_sgfa => basis_set_a%first_sgf
2901 la_max => basis_set_a%lmax
2902 la_min => basis_set_a%lmin
2903 nseta = basis_set_a%nset
2904 nsgfa => basis_set_a%nsgf_set
2905 sphi_a => basis_set_a%sphi
2906 zeta => basis_set_a%zet
2907 npgfa => basis_set_a%npgf
2908
2909 ra(:) = pbc(particle_set(iatom)%r, cell)
2910
2911 CALL dbcsr_get_block_p(tmp_g_pq, iatom, jatom, pblock, found)
2912 IF (.NOT. found) cycle
2913
2914 offset = 0
2915 DO iset = 1, nseta
2916 ncoa = npgfa(iset)*ncoset(la_max(iset))
2917 sgfa = first_sgfa(1, iset)
2918
2919 ALLOCATE (h_tmp(ncoa, 1)); h_tmp = 0.0_dp
2920 ALLOCATE (i_ab(nsgfa(iset), 1)); i_ab = 0.0_dp
2921 ALLOCATE (pab(ncoa, 1)); pab = 0.0_dp
2922
2923 i_ab(1:nsgfa(iset), 1) = 2.0_dp*pblock(offset + 1:offset + nsgfa(iset), j_ri - lb_ri)
2924 CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
2925 i_ab(1, 1), nsgfa(iset), 0.0_dp, pab(1, 1), ncoa)
2926
2927 igrid_level = gaussian_gridlevel(pw_env_ext%gridlevel_info, minval(zeta(:, iset)))
2928
2929 ! The last three parameters are used to check whether a given function is within the own range.
2930 ! Here, it is always the case, so let's enforce it because mod(0, 1)==0
2931 IF (map_gaussian_here(rs_v(igrid_level), cell%h_inv, ra, 0, 1, 0)) THEN
2932 DO ipgf = 1, npgfa(iset)
2933 o1 = (ipgf - 1)*ncoset(la_max(iset))
2934 igrid_level = gaussian_gridlevel(pw_env_ext%gridlevel_info, zeta(ipgf, iset))
2935
2936 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
2937 lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
2938 zetp=zeta(ipgf, iset), &
2939 eps=dft_control%qs_control%eps_gvg_rspace, &
2940 prefactor=1.0_dp, cutoff=1.0_dp)
2941
2942 CALL integrate_pgf_product( &
2943 la_max=la_max(iset), zeta=zeta(ipgf, iset), la_min=la_min(iset), &
2944 lb_max=0, zetb=0.0_dp, lb_min=0, &
2945 ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
2946 rsgrid=rs_v(igrid_level), &
2947 hab=h_tmp, pab=pab, &
2948 o1=o1, &
2949 o2=0, &
2950 radius=radius, &
2951 calculate_forces=.true., &
2952 force_a=force_a, force_b=force_b, &
2953 use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
2954
2955 END DO
2956
2957 END IF
2958
2959 offset = offset + nsgfa(iset)
2960 DEALLOCATE (pab, h_tmp, i_ab)
2961 END DO !iset
2962
2963 force(ikind)%mp2_non_sep(:, atom_a) = force(ikind)%mp2_non_sep(:, atom_a) + force_a + force_b
2964 IF (use_virial) h_stress = h_stress + my_virial_a + my_virial_b
2965
2966 END DO !iatom
2967 END DO !j_RI
2968 END DO !jatom
2969
2970 IF (use_virial) THEN
2971 CALL auxbas_pw_pool%give_back_pw(rho_g_copy)
2972 DO i_xyz = 1, 3
2973 CALL auxbas_pw_pool%give_back_pw(dvg(i_xyz))
2974 END DO
2975 END IF
2976
2977 CALL cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_ext, pw_env_ext, &
2978 task_list_ext, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_l)
2979
2980 CALL dbcsr_release(tmp_g_pq)
2981 CALL dbcsr_distribution_release(dbcsr_dist)
2982 DEALLOCATE (col_dist, row_dist, pgrid)
2983
2984 CALL mp_para_env_release(para_env_ext)
2985
2986 CALL timestop(handle)
2987
2988 END SUBROUTINE get_2c_gpw_forces
2989
2990! **************************************************************************************************
2991!> \brief Calculate the forces due to the (P|Q) MME integral derivatives
2992!> \param G_PQ ...
2993!> \param force ...
2994!> \param mp2_env ...
2995!> \param qs_env ...
2996! **************************************************************************************************
2997 SUBROUTINE get_2c_mme_forces(G_PQ, force, mp2_env, qs_env)
2998
2999 TYPE(dbcsr_type), INTENT(INOUT) :: g_pq
3000 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3001 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
3002 TYPE(qs_environment_type), POINTER :: qs_env
3003
3004 CHARACTER(len=*), PARAMETER :: routinen = 'get_2c_mme_forces'
3005
3006 INTEGER :: atom_a, atom_b, g_count, handle, i_xyz, iatom, ikind, iset, jatom, jkind, jset, &
3007 natom, nkind, nseta, nsetb, offset_hab_a, offset_hab_b, r_count, sgfa, sgfb
3008 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
3009 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
3010 npgfb, nsgfa, nsgfb
3011 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
3012 LOGICAL :: found
3013 REAL(dp) :: new_force, pref
3014 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: hab
3015 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: hdab
3016 REAL(dp), DIMENSION(:, :), POINTER :: pblock
3017 REAL(kind=dp), DIMENSION(3) :: ra, rb
3018 REAL(kind=dp), DIMENSION(:, :), POINTER :: sphi_a, sphi_b, zeta, zetb
3019 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3020 TYPE(cell_type), POINTER :: cell
3021 TYPE(dbcsr_iterator_type) :: iter
3022 TYPE(gto_basis_set_p_type), ALLOCATABLE, &
3023 DIMENSION(:), TARGET :: basis_set_ri_aux
3024 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
3025 TYPE(mp_para_env_type), POINTER :: para_env
3026 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3027 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3028
3029 NULLIFY (qs_kind_set, basis_set_a, basis_set_b, pblock, particle_set, &
3030 cell, la_max, la_min, lb_min, npgfa, lb_max, npgfb, nsgfa, &
3031 nsgfb, first_sgfa, first_sgfb, sphi_a, sphi_b, zeta, zetb, &
3032 atomic_kind_set, para_env)
3033
3034 CALL timeset(routinen, handle)
3035
3036 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind, particle_set=particle_set, &
3037 cell=cell, atomic_kind_set=atomic_kind_set, natom=natom, para_env=para_env)
3038
3039 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
3040
3041 ALLOCATE (basis_set_ri_aux(nkind))
3042 CALL basis_set_list_setup(basis_set_ri_aux, "RI_AUX", qs_kind_set)
3043
3044 g_count = 0; r_count = 0
3045
3046 CALL dbcsr_iterator_start(iter, g_pq)
3047 DO WHILE (dbcsr_iterator_blocks_left(iter))
3048
3049 CALL dbcsr_iterator_next_block(iter, row=iatom, column=jatom)
3050 CALL dbcsr_get_block_p(g_pq, iatom, jatom, pblock, found)
3051 IF (.NOT. found) cycle
3052 IF (iatom > jatom) cycle
3053 pref = 2.0_dp
3054 IF (iatom == jatom) pref = 1.0_dp
3055
3056 ikind = kind_of(iatom)
3057 jkind = kind_of(jatom)
3058
3059 atom_a = atom_of_kind(iatom)
3060 atom_b = atom_of_kind(jatom)
3061
3062 basis_set_a => basis_set_ri_aux(ikind)%gto_basis_set
3063 first_sgfa => basis_set_a%first_sgf
3064 la_max => basis_set_a%lmax
3065 la_min => basis_set_a%lmin
3066 nseta = basis_set_a%nset
3067 nsgfa => basis_set_a%nsgf_set
3068 sphi_a => basis_set_a%sphi
3069 zeta => basis_set_a%zet
3070 npgfa => basis_set_a%npgf
3071
3072 basis_set_b => basis_set_ri_aux(jkind)%gto_basis_set
3073 first_sgfb => basis_set_b%first_sgf
3074 lb_max => basis_set_b%lmax
3075 lb_min => basis_set_b%lmin
3076 nsetb = basis_set_b%nset
3077 nsgfb => basis_set_b%nsgf_set
3078 sphi_b => basis_set_b%sphi
3079 zetb => basis_set_b%zet
3080 npgfb => basis_set_b%npgf
3081
3082 ra(:) = pbc(particle_set(iatom)%r, cell)
3083 rb(:) = pbc(particle_set(jatom)%r, cell)
3084
3085 ALLOCATE (hab(basis_set_a%nsgf, basis_set_b%nsgf))
3086 ALLOCATE (hdab(3, basis_set_a%nsgf, basis_set_b%nsgf))
3087 hab(:, :) = 0.0_dp
3088 hdab(:, :, :) = 0.0_dp
3089
3090 offset_hab_a = 0
3091 DO iset = 1, nseta
3092 sgfa = first_sgfa(1, iset)
3093
3094 offset_hab_b = 0
3095 DO jset = 1, nsetb
3096 sgfb = first_sgfb(1, jset)
3097
3098 CALL integrate_set_2c(mp2_env%eri_mme_param%par, mp2_env%potential_parameter, la_min(iset), &
3099 la_max(iset), lb_min(jset), lb_max(jset), npgfa(iset), npgfb(jset), &
3100 zeta(:, iset), zetb(:, jset), ra, rb, hab, nsgfa(iset), nsgfb(jset), &
3101 offset_hab_a, offset_hab_b, 0, 0, sphi_a, sphi_b, sgfa, sgfb, &
3102 nsgfa(iset), nsgfb(jset), do_eri_mme, hdab=hdab, &
3103 g_count=g_count, r_count=r_count)
3104
3105 offset_hab_b = offset_hab_b + nsgfb(jset)
3106 END DO
3107 offset_hab_a = offset_hab_a + nsgfa(iset)
3108 END DO
3109
3110 DO i_xyz = 1, 3
3111 new_force = pref*sum(pblock(:, :)*hdab(i_xyz, :, :))
3112 force(ikind)%mp2_non_sep(i_xyz, atom_a) = force(ikind)%mp2_non_sep(i_xyz, atom_a) + new_force
3113 force(jkind)%mp2_non_sep(i_xyz, atom_b) = force(jkind)%mp2_non_sep(i_xyz, atom_b) - new_force
3114 END DO
3115
3116 DEALLOCATE (hab, hdab)
3117 END DO
3118 CALL dbcsr_iterator_stop(iter)
3119
3120 CALL cp_eri_mme_update_local_counts(mp2_env%eri_mme_param, para_env, g_count_2c=g_count, r_count_2c=r_count)
3121
3122 CALL timestop(handle)
3123
3124 END SUBROUTINE get_2c_mme_forces
3125
3126! **************************************************************************************************
3127!> \brief This routines gather all the force updates due to the response density and the trace with F
3128!> Also update the forces due to the SCF density for XC and exact exchange
3129!> \param p_env the p_env coming from the response calculation
3130!> \param matrix_hz the matrix going into the RHS of the response equation
3131!> \param matrix_p_F the density matrix with which we evaluate Trace[P*F]
3132!> \param matrix_p_F_admm ...
3133!> \param qs_env ...
3134!> \note very much inspired from the response_force routine in response_solver.F, especially for virial
3135! **************************************************************************************************
3136 SUBROUTINE update_im_time_forces(p_env, matrix_hz, matrix_p_F, matrix_p_F_admm, qs_env)
3137
3138 TYPE(qs_p_env_type), POINTER :: p_env
3139 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz, matrix_p_f, matrix_p_f_admm
3140 TYPE(qs_environment_type), POINTER :: qs_env
3141
3142 CHARACTER(len=*), PARAMETER :: routinen = 'update_im_time_forces'
3143
3144 INTEGER :: handle, i, idens, ispin, n_rep_hf, nao, &
3145 nao_aux, nder, nimages, nocc, nspins
3146 LOGICAL :: do_exx, do_hfx, do_tau, do_tau_admm, &
3147 use_virial
3148 REAL(dp) :: dummy_real1, dummy_real2, ehartree, exc, &
3149 focc
3150 REAL(dp), DIMENSION(3, 3) :: h_stress, pv_loc
3151 TYPE(admm_type), POINTER :: admm_env
3152 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3153 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: current_density, current_density_admm, &
3154 current_mat_h, matrix_p_mp2, matrix_p_mp2_admm, matrix_s, matrix_s_aux_fit, matrix_w, &
3155 rho_ao, rho_ao_aux, scrm, scrm_admm
3156 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: dbcsr_work_h, dbcsr_work_p, mpa2
3157 TYPE(dbcsr_type) :: dbcsr_work
3158 TYPE(dft_control_type), POINTER :: dft_control
3159 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
3160 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3161 TYPE(mp_para_env_type), POINTER :: para_env
3162 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3163 POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
3164 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3165 TYPE(pw_c1d_gs_type) :: rho_tot_gspace, rhoz_tot_gspace, &
3166 zv_hartree_gspace
3167 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rhoz_g
3168 TYPE(pw_env_type), POINTER :: pw_env
3169 TYPE(pw_poisson_type), POINTER :: poisson_env
3170 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
3171 TYPE(pw_r3d_rs_type) :: vh_rspace, vhxc_rspace, zv_hartree_rspace
3172 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rhoz_r, tauz_r, v_xc, v_xc_tau, &
3173 vadmm_rspace, vtau_rspace, vxc_rspace
3174 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
3175 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3176 TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit, rhoz
3177 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set
3178 TYPE(section_vals_type), POINTER :: hfx_section, xc_section
3179 TYPE(task_list_type), POINTER :: task_list_aux_fit
3180 TYPE(virial_type), POINTER :: virial
3181
3182 NULLIFY (scrm, rho, dft_control, matrix_p_mp2, matrix_s, &
3183 matrix_p_mp2_admm, admm_env, sab_orb, dbcsr_work_p, &
3184 dbcsr_work_h, sac_ae, sac_ppl, sap_ppnl, force, virial, &
3185 qs_kind_set, atomic_kind_set, particle_set, pw_env, poisson_env, &
3186 auxbas_pw_pool, task_list_aux_fit, matrix_s_aux_fit, scrm_admm, &
3187 rho_aux_fit, rho_ao_aux, x_data, hfx_section, xc_section, &
3188 para_env, rhoz_g, rhoz_r, tauz_r, v_xc, v_xc_tau, vxc_rspace, &
3189 vtau_rspace, vadmm_rspace, rho_ao, matrix_w)
3190 NULLIFY (rho0_atom_set, rho1_atom_set)
3191
3192 CALL timeset(routinen, handle)
3193
3194 CALL get_qs_env(qs_env, rho=rho, dft_control=dft_control, matrix_s=matrix_s, admm_env=admm_env, &
3195 sab_orb=sab_orb, sac_ae=sac_ae, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl, force=force, &
3196 virial=virial, particle_set=particle_set, qs_kind_set=qs_kind_set, &
3197 atomic_kind_set=atomic_kind_set, x_data=x_data, para_env=para_env)
3198 nspins = dft_control%nspins
3199
3200 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
3201 IF (use_virial) virial%pv_calculate = .true.
3202
3203 !Whether we replace the force/energy of SCF XC with HF in RPA
3204 do_exx = .false.
3205 IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
3206 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
3207 CALL section_vals_get(hfx_section, explicit=do_exx)
3208 END IF
3209
3210 !Get the mp2 density matrix which is p_env%p1 + matrix_p_F
3211 CALL get_qs_env(qs_env, matrix_p_mp2=matrix_p_mp2, matrix_p_mp2_admm=matrix_p_mp2_admm)
3212
3213 !The kinetic term (only response density)
3214 NULLIFY (scrm)
3215 mpa2(1:nspins, 1:1) => matrix_p_mp2(1:nspins)
3216 CALL kinetic_energy_matrix(qs_env, matrix_t=scrm, matrix_p=mpa2, &
3217 matrix_name="KINETIC ENERGY MATRIX", &
3218 basis_type="ORB", &
3219 sab_orb=sab_orb, calculate_forces=.true.)
3220 CALL dbcsr_deallocate_matrix_set(scrm)
3221
3222 !The pseudo-potential terms (only reponse density)
3223 CALL dbcsr_allocate_matrix_set(scrm, nspins)
3224 DO ispin = 1, nspins
3225 ALLOCATE (scrm(ispin)%matrix)
3226 CALL dbcsr_create(scrm(ispin)%matrix, template=matrix_s(1)%matrix)
3227 CALL dbcsr_copy(scrm(ispin)%matrix, matrix_s(1)%matrix)
3228 CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
3229 END DO
3230
3231 nder = 1
3232 nimages = 1
3233 ALLOCATE (dbcsr_work_p(nspins, 1), dbcsr_work_h(nspins, 1))
3234 DO ispin = 1, nspins
3235 dbcsr_work_p(ispin, 1)%matrix => matrix_p_mp2(ispin)%matrix
3236 dbcsr_work_h(ispin, 1)%matrix => scrm(ispin)%matrix
3237 END DO
3238
3239 CALL core_matrices(qs_env, dbcsr_work_h, dbcsr_work_p, .true., nder)
3240
3241 DEALLOCATE (dbcsr_work_p, dbcsr_work_h)
3242
3243 IF (use_virial) THEN
3244 h_stress = 0.0_dp
3245 virial%pv_xc = 0.0_dp
3246 NULLIFY (vxc_rspace, vtau_rspace, vadmm_rspace)
3247 CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, &
3248 dummy_real1, dummy_real2, h_stress)
3249 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
3250 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
3251 IF (.NOT. do_exx) THEN
3252 !if RPA EXX, then do not consider XC virial (replaced by RPA%HF virial)
3253 virial%pv_exc = virial%pv_exc - virial%pv_xc
3254 virial%pv_virial = virial%pv_virial - virial%pv_xc
3255 END IF
3256 ELSE
3257 CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, dummy_real1, dummy_real2)
3258 END IF
3259 do_tau = ASSOCIATED(vtau_rspace)
3260
3261 !Core forces from the SCF
3262 CALL integrate_v_core_rspace(vh_rspace, qs_env)
3263
3264 !The Hartree-xc potential term, P*dVHxc (mp2 + SCF density x deriv of the SCF potential)
3265 !Get the total density
3266 CALL qs_rho_get(rho, rho_ao=rho_ao)
3267 DO ispin = 1, nspins
3268 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
3269 END DO
3270
3271 CALL get_qs_env(qs_env, pw_env=pw_env)
3272 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
3273 poisson_env=poisson_env)
3274 CALL auxbas_pw_pool%create_pw(vhxc_rspace)
3275
3276 IF (use_virial) pv_loc = virial%pv_virial
3277
3278 IF (do_exx) THEN
3279 !Only want response XC contribution, but SCF+response Hartree contribution
3280 DO ispin = 1, nspins
3281 !Hartree
3282 CALL pw_transfer(vh_rspace, vhxc_rspace)
3283 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
3284 hmat=scrm(ispin), pmat=rho_ao(ispin), &
3285 qs_env=qs_env, calculate_forces=.true.)
3286 !XC
3287 CALL pw_transfer(vxc_rspace(ispin), vhxc_rspace)
3288 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
3289 hmat=scrm(ispin), pmat=matrix_p_mp2(ispin), &
3290 qs_env=qs_env, calculate_forces=.true.)
3291 IF (do_tau) THEN
3292 CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
3293 hmat=scrm(ispin), pmat=matrix_p_mp2(ispin), &
3294 qs_env=qs_env, calculate_forces=.true., compute_tau=.true.)
3295 END IF
3296 END DO
3297 ELSE
3298 DO ispin = 1, nspins
3299 CALL pw_transfer(vh_rspace, vhxc_rspace)
3300 CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
3301 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
3302 hmat=scrm(ispin), pmat=rho_ao(ispin), &
3303 qs_env=qs_env, calculate_forces=.true.)
3304 IF (do_tau) THEN
3305 CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
3306 hmat=scrm(ispin), pmat=rho_ao(ispin), &
3307 qs_env=qs_env, calculate_forces=.true., compute_tau=.true.)
3308 END IF
3309 END DO
3310 END IF
3311 CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
3312
3313 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3314
3315 !The admm projection contribution (mp2 + SCF densities). If EXX, then only mp2 density
3316 IF (dft_control%do_admm) THEN
3317 CALL get_admm_env(admm_env, task_list_aux_fit=task_list_aux_fit, rho_aux_fit=rho_aux_fit, &
3318 matrix_s_aux_fit=matrix_s_aux_fit)
3319 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux)
3320 CALL dbcsr_allocate_matrix_set(scrm_admm, nspins)
3321 DO ispin = 1, nspins
3322 ALLOCATE (scrm_admm(ispin)%matrix)
3323 CALL dbcsr_create(scrm_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
3324 CALL dbcsr_copy(scrm_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
3325 CALL dbcsr_set(scrm_admm(ispin)%matrix, 0.0_dp)
3326 END DO
3327
3328 IF (use_virial) pv_loc = virial%pv_virial
3329 IF (.NOT. qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
3330 DO ispin = 1, nspins
3331 IF (do_exx) THEN
3332 CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), &
3333 hmat=scrm_admm(ispin), pmat=matrix_p_mp2_admm(ispin), &
3334 qs_env=qs_env, calculate_forces=.true., &
3335 basis_type="AUX_FIT", task_list_external=task_list_aux_fit)
3336 ELSE
3337 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
3338 CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), &
3339 hmat=scrm_admm(ispin), pmat=rho_ao_aux(ispin), &
3340 qs_env=qs_env, calculate_forces=.true., &
3341 basis_type="AUX_FIT", task_list_external=task_list_aux_fit)
3342 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
3343 END IF
3344 END DO
3345 END IF
3346 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3347
3348 CALL tddft_hfx_matrix(scrm_admm, rho_ao_aux, qs_env, .false., .false.)
3349
3350 IF (do_exx) THEN
3351 CALL admm_projection_derivative(qs_env, scrm_admm, matrix_p_mp2)
3352 ELSE
3353 CALL admm_projection_derivative(qs_env, scrm_admm, rho_ao)
3354 END IF
3355 END IF
3356
3357 !The exact-exchange term (mp2 + SCF densities)
3358 xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
3359 hfx_section => section_vals_get_subs_vals(xc_section, "HF")
3360 CALL section_vals_get(hfx_section, explicit=do_hfx)
3361
3362 IF (do_hfx) THEN
3363 CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
3364 cpassert(n_rep_hf == 1)
3365 IF (use_virial) virial%pv_fock_4c = 0.0_dp
3366
3367 !In case of EXX, only want to response HFX forces, as the SCF will change according to RI_RPA%HF
3368 IF (do_exx) THEN
3369 IF (dft_control%do_admm) THEN
3370 CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
3371 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux, rho_ao_kp=dbcsr_work_p)
3372 IF (x_data(1, 1)%do_hfx_ri) THEN
3373
3374 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3375 x_data(1, 1)%general_parameter%fraction, &
3376 rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2_admm, &
3377 use_virial=use_virial, resp_only=.true.)
3378 ELSE
3379 CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2_admm, hfx_section, para_env, &
3380 1, use_virial, resp_only=.true.)
3381 END IF
3382 ELSE
3383 DO ispin = 1, nspins
3384 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
3385 END DO
3386 CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3387 IF (x_data(1, 1)%do_hfx_ri) THEN
3388
3389 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3390 x_data(1, 1)%general_parameter%fraction, &
3391 rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2, &
3392 use_virial=use_virial, resp_only=.true.)
3393 ELSE
3394 CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2, hfx_section, para_env, &
3395 1, use_virial, resp_only=.true.)
3396 END IF
3397 DO ispin = 1, nspins
3398 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
3399 END DO
3400 END IF !admm
3401
3402 ELSE !No Exx
3403 IF (dft_control%do_admm) THEN
3404 CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
3405 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux, rho_ao_kp=dbcsr_work_p)
3406 DO ispin = 1, nspins
3407 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
3408 END DO
3409 IF (x_data(1, 1)%do_hfx_ri) THEN
3410
3411 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3412 x_data(1, 1)%general_parameter%fraction, &
3413 rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2_admm, &
3414 use_virial=use_virial, resp_only=.false.)
3415 ELSE
3416 CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2_admm, hfx_section, para_env, &
3417 1, use_virial, resp_only=.false.)
3418 END IF
3419 DO ispin = 1, nspins
3420 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
3421 END DO
3422 ELSE
3423 CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3424 IF (x_data(1, 1)%do_hfx_ri) THEN
3425
3426 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3427 x_data(1, 1)%general_parameter%fraction, &
3428 rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2, &
3429 use_virial=use_virial, resp_only=.false.)
3430 ELSE
3431 CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2, hfx_section, para_env, &
3432 1, use_virial, resp_only=.false.)
3433 END IF
3434 END IF
3435 END IF !do_exx
3436
3437 IF (use_virial) THEN
3438 virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
3439 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
3440 END IF
3441 END IF
3442
3443 !retrieve the SCF density
3444 CALL qs_rho_get(rho, rho_ao=rho_ao)
3445 DO ispin = 1, nspins
3446 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
3447 END DO
3448
3449 !From here, we need to do everything twice. Once for the response density, and once for the
3450 !density that is used for the trace Tr[P*F]. The reason is that the former is needed for the
3451 !eventual overlap contribution from matrix_wz
3452 !Only with the mp2 density
3453
3454 ALLOCATE (current_density(nspins), current_mat_h(nspins), current_density_admm(nspins))
3455 DO idens = 1, 2
3456 DO ispin = 1, nspins
3457 IF (idens == 1) THEN
3458 current_density(ispin)%matrix => matrix_p_f(ispin)%matrix
3459 current_mat_h(ispin)%matrix => scrm(ispin)%matrix
3460 IF (dft_control%do_admm) current_density_admm(ispin)%matrix => matrix_p_f_admm(ispin)%matrix
3461 ELSE
3462 current_density(ispin)%matrix => p_env%p1(ispin)%matrix
3463 current_mat_h(ispin)%matrix => matrix_hz(ispin)%matrix
3464 IF (dft_control%do_admm) current_density_admm(ispin)%matrix => p_env%p1_admm(ispin)%matrix
3465 END IF
3466 END DO
3467
3468 !The core-denstiy derivative
3469 ALLOCATE (rhoz_r(nspins), rhoz_g(nspins))
3470 DO ispin = 1, nspins
3471 CALL auxbas_pw_pool%create_pw(rhoz_r(ispin))
3472 CALL auxbas_pw_pool%create_pw(rhoz_g(ispin))
3473 END DO
3474 CALL auxbas_pw_pool%create_pw(rhoz_tot_gspace)
3475 CALL auxbas_pw_pool%create_pw(zv_hartree_rspace)
3476 CALL auxbas_pw_pool%create_pw(zv_hartree_gspace)
3477
3478 CALL pw_zero(rhoz_tot_gspace)
3479 DO ispin = 1, nspins
3480 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density(ispin)%matrix, &
3481 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin))
3482 CALL pw_axpy(rhoz_g(ispin), rhoz_tot_gspace)
3483 END DO
3484
3485 IF (use_virial) THEN
3486
3487 CALL get_qs_env(qs_env, rho=rho)
3488 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
3489
3490 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
3491
3492 h_stress(:, :) = 0.0_dp
3493 CALL pw_poisson_solve(poisson_env, &
3494 density=rhoz_tot_gspace, &
3495 ehartree=ehartree, &
3496 vhartree=zv_hartree_gspace, &
3497 h_stress=h_stress, &
3498 aux_density=rho_tot_gspace)
3499
3500 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
3501
3502 !Green contribution
3503 virial%pv_ehartree = virial%pv_ehartree + 2.0_dp*h_stress/real(para_env%num_pe, dp)
3504 virial%pv_virial = virial%pv_virial + 2.0_dp*h_stress/real(para_env%num_pe, dp)
3505
3506 ELSE
3507 CALL pw_poisson_solve(poisson_env, rhoz_tot_gspace, ehartree, &
3508 zv_hartree_gspace)
3509 END IF
3510
3511 CALL pw_transfer(zv_hartree_gspace, zv_hartree_rspace)
3512 CALL pw_scale(zv_hartree_rspace, zv_hartree_rspace%pw_grid%dvol)
3513 CALL integrate_v_core_rspace(zv_hartree_rspace, qs_env)
3514
3515 IF (do_tau) THEN
3516 block
3517 TYPE(pw_c1d_gs_type) :: tauz_g
3518 CALL auxbas_pw_pool%create_pw(tauz_g)
3519 ALLOCATE (tauz_r(nspins))
3520 DO ispin = 1, nspins
3521 CALL auxbas_pw_pool%create_pw(tauz_r(ispin))
3522
3523 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density(ispin)%matrix, &
3524 rho=tauz_r(ispin), rho_gspace=tauz_g, compute_tau=.true.)
3525 END DO
3526 CALL auxbas_pw_pool%give_back_pw(tauz_g)
3527 END block
3528 END IF
3529
3530 !Volume contribution to the virial
3531 IF (use_virial) THEN
3532 !Volume contribution
3533 exc = 0.0_dp
3534 DO ispin = 1, nspins
3535 exc = exc + pw_integral_ab(rhoz_r(ispin), vxc_rspace(ispin))/ &
3536 vxc_rspace(ispin)%pw_grid%dvol
3537 END DO
3538 IF (ASSOCIATED(vtau_rspace)) THEN
3539 DO ispin = 1, nspins
3540 exc = exc + pw_integral_ab(tauz_r(ispin), vtau_rspace(ispin))/ &
3541 vtau_rspace(ispin)%pw_grid%dvol
3542 END DO
3543 END IF
3544 DO i = 1, 3
3545 virial%pv_ehartree(i, i) = virial%pv_ehartree(i, i) - 4.0_dp*ehartree/real(para_env%num_pe, dp)
3546 virial%pv_exc(i, i) = virial%pv_exc(i, i) - exc/real(para_env%num_pe, dp)
3547 virial%pv_virial(i, i) = virial%pv_virial(i, i) - 4.0_dp*ehartree/real(para_env%num_pe, dp) &
3548 - exc/real(para_env%num_pe, dp)
3549 END DO
3550 END IF
3551
3552 !The xc-kernel term.
3553 IF (dft_control%do_admm) THEN
3554 CALL get_qs_env(qs_env, admm_env=admm_env)
3555 xc_section => admm_env%xc_section_primary
3556 ELSE
3557 xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
3558 END IF
3559
3560 IF (use_virial) virial%pv_xc = 0.0_dp
3561
3562 ALLOCATE (rhoz)
3563 CALL qs_rho_create(rhoz)
3564 IF (ASSOCIATED(rhoz_r)) THEN
3565 CALL qs_rho_set(rhoz, rho_r=rhoz_r, rho_r_valid=.true.)
3566 END IF
3567 IF (ASSOCIATED(rhoz_g)) THEN
3568 CALL qs_rho_set(rhoz, rho_g=rhoz_g, rho_g_valid=.true.)
3569 END IF
3570 IF (ASSOCIATED(tauz_r)) THEN
3571 CALL qs_rho_set(rhoz, tau_r=tauz_r, tau_r_valid=.true.)
3572 END IF
3573 !
3574 CALL qs_fxc_create(qs_env, rho, rhoz, rho0_atom_set, xc_section, .false., &
3575 v_xc, v_xc_tau, rho1_atom_set, &
3576 dispersion_env=qs_env%dispersion_env, &
3577 compute_virial=use_virial, virial_xc=virial%pv_xc)
3578 !
3579 DEALLOCATE (rhoz)
3580
3581 IF (use_virial) THEN
3582 virial%pv_exc = virial%pv_exc + virial%pv_xc
3583 virial%pv_virial = virial%pv_virial + virial%pv_xc
3584
3585 pv_loc = virial%pv_virial
3586 END IF
3587
3588 CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3589 DO ispin = 1, nspins
3590 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
3591 CALL pw_axpy(zv_hartree_rspace, v_xc(ispin))
3592 CALL integrate_v_rspace(qs_env=qs_env, &
3593 v_rspace=v_xc(ispin), &
3594 hmat=current_mat_h(ispin), &
3595 pmat=dbcsr_work_p(ispin, 1), &
3596 calculate_forces=.true.)
3597 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
3598 END DO
3599 CALL auxbas_pw_pool%give_back_pw(rhoz_tot_gspace)
3600 CALL auxbas_pw_pool%give_back_pw(zv_hartree_rspace)
3601 CALL auxbas_pw_pool%give_back_pw(zv_hartree_gspace)
3602 DEALLOCATE (v_xc)
3603
3604 IF (do_tau) THEN
3605 DO ispin = 1, nspins
3606 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
3607 CALL integrate_v_rspace(qs_env=qs_env, &
3608 v_rspace=v_xc_tau(ispin), &
3609 hmat=current_mat_h(ispin), &
3610 pmat=dbcsr_work_p(ispin, 1), &
3611 compute_tau=.true., &
3612 calculate_forces=.true.)
3613 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
3614 END DO
3615 DEALLOCATE (v_xc_tau)
3616 END IF
3617
3618 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3619
3620 IF (do_hfx) THEN
3621 IF (dft_control%do_admm) THEN
3622 DO ispin = 1, nspins
3623 CALL dbcsr_set(scrm_admm(ispin)%matrix, 0.0_dp)
3624 END DO
3625 CALL qs_rho_get(rho_aux_fit, tau_r_valid=do_tau_admm)
3626
3627 IF (.NOT. admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
3628 CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
3629 DO ispin = 1, nspins
3630 CALL pw_zero(rhoz_r(ispin))
3631 CALL pw_zero(rhoz_g(ispin))
3632 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density_admm(ispin)%matrix, &
3633 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin), &
3634 basis_type="AUX_FIT", task_list_external=task_list_aux_fit)
3635 END DO
3636
3637 IF (do_tau_admm) THEN
3638 block
3639 TYPE(pw_c1d_gs_type) :: tauz_g
3640 CALL auxbas_pw_pool%create_pw(tauz_g)
3641 DO ispin = 1, nspins
3642 CALL pw_zero(tauz_r(ispin))
3643 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density(ispin)%matrix, &
3644 rho=tauz_r(ispin), rho_gspace=tauz_g, &
3645 basis_type="AUX_FIT", task_list_external=task_list_aux_fit, &
3646 compute_tau=.true.)
3647 END DO
3648 CALL auxbas_pw_pool%give_back_pw(tauz_g)
3649 END block
3650 END IF
3651
3652 !Volume contribution to the virial
3653 IF (use_virial) THEN
3654 exc = 0.0_dp
3655 DO ispin = 1, nspins
3656 exc = exc + pw_integral_ab(rhoz_r(ispin), vadmm_rspace(ispin))/ &
3657 vadmm_rspace(ispin)%pw_grid%dvol
3658 END DO
3659 DO i = 1, 3
3660 virial%pv_exc(i, i) = virial%pv_exc(i, i) - exc/real(para_env%num_pe, dp)
3661 virial%pv_virial(i, i) = virial%pv_virial(i, i) - exc/real(para_env%num_pe, dp)
3662 END DO
3663
3664 virial%pv_xc = 0.0_dp
3665 END IF
3666
3667 xc_section => admm_env%xc_section_aux
3668 ALLOCATE (rhoz)
3669 CALL qs_rho_create(rhoz)
3670 IF (ASSOCIATED(rhoz_r)) THEN
3671 CALL qs_rho_set(rhoz, rho_r=rhoz_r, rho_r_valid=.true.)
3672 END IF
3673 IF (ASSOCIATED(rhoz_g)) THEN
3674 CALL qs_rho_set(rhoz, rho_g=rhoz_g, rho_g_valid=.true.)
3675 END IF
3676 IF (ASSOCIATED(tauz_r)) THEN
3677 CALL qs_rho_set(rhoz, tau_r=tauz_r, tau_r_valid=.true.)
3678 END IF
3679 !
3680 CALL qs_fxc_create(qs_env, rho_aux_fit, rhoz, rho0_atom_set, xc_section, .false., &
3681 v_xc, v_xc_tau, rho1_atom_set, &
3682 compute_virial=use_virial, virial_xc=virial%pv_xc)
3683 !
3684 DEALLOCATE (rhoz)
3685
3686 IF (use_virial) THEN
3687 virial%pv_exc = virial%pv_exc + virial%pv_xc
3688 virial%pv_virial = virial%pv_virial + virial%pv_xc
3689
3690 pv_loc = virial%pv_virial
3691 END IF
3692
3693 CALL qs_rho_get(rho_aux_fit, rho_ao_kp=dbcsr_work_p)
3694 DO ispin = 1, nspins
3695 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
3696 CALL integrate_v_rspace(qs_env=qs_env, &
3697 v_rspace=v_xc(ispin), &
3698 hmat=scrm_admm(ispin), &
3699 pmat=dbcsr_work_p(ispin, 1), &
3700 calculate_forces=.true., &
3701 basis_type="AUX_FIT", &
3702 task_list_external=task_list_aux_fit)
3703 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
3704 END DO
3705 DEALLOCATE (v_xc)
3706
3707 IF (do_tau_admm) THEN
3708 DO ispin = 1, nspins
3709 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
3710 CALL integrate_v_rspace(qs_env=qs_env, &
3711 v_rspace=v_xc_tau(ispin), &
3712 hmat=scrm_admm(ispin), &
3713 pmat=dbcsr_work_p(ispin, 1), &
3714 calculate_forces=.true., &
3715 basis_type="AUX_FIT", &
3716 task_list_external=task_list_aux_fit, &
3717 compute_tau=.true.)
3718 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
3719 END DO
3720 DEALLOCATE (v_xc_tau)
3721 END IF
3722
3723 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3724 END IF
3725
3726 CALL tddft_hfx_matrix(scrm_admm, current_density_admm, qs_env, .false., .false.)
3727
3728 CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3729 CALL admm_projection_derivative(qs_env, scrm_admm, dbcsr_work_p(:, 1))
3730
3731 !If response density, need to get matrix_hz contribution
3732 CALL dbcsr_create(dbcsr_work, template=matrix_s(1)%matrix)
3733 IF (idens == 2) THEN
3734 nao = admm_env%nao_orb
3735 nao_aux = admm_env%nao_aux_fit
3736 DO ispin = 1, nspins
3737 CALL dbcsr_copy(dbcsr_work, matrix_hz(ispin)%matrix)
3738 CALL dbcsr_set(dbcsr_work, 0.0_dp)
3739
3740 CALL cp_dbcsr_sm_fm_multiply(scrm_admm(ispin)%matrix, admm_env%A, &
3741 admm_env%work_aux_orb, nao)
3742 CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
3743 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
3744 admm_env%work_orb_orb)
3745 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_work, keep_sparsity=.true.)
3746 CALL dbcsr_add(matrix_hz(ispin)%matrix, dbcsr_work, 1.0_dp, 1.0_dp)
3747 END DO
3748 END IF
3749
3750 CALL dbcsr_release(dbcsr_work)
3751 ELSE !no admm
3752
3753 !Need the contribution to matrix_hz as well
3754 IF (idens == 2) THEN
3755 CALL tddft_hfx_matrix(matrix_hz, current_density, qs_env, .false., .false.)
3756 END IF
3757 END IF !admm
3758 END IF !do_hfx
3759
3760 DO ispin = 1, nspins
3761 CALL auxbas_pw_pool%give_back_pw(rhoz_r(ispin))
3762 CALL auxbas_pw_pool%give_back_pw(rhoz_g(ispin))
3763 END DO
3764 DEALLOCATE (rhoz_r, rhoz_g)
3765
3766 IF (do_tau) THEN
3767 DO ispin = 1, nspins
3768 CALL auxbas_pw_pool%give_back_pw(tauz_r(ispin))
3769 END DO
3770 DEALLOCATE (tauz_r)
3771 END IF
3772 END DO !idens
3773 CALL dbcsr_deallocate_matrix_set(scrm_admm)
3774
3775 DEALLOCATE (current_density, current_mat_h, current_density_admm)
3776 CALL dbcsr_deallocate_matrix_set(scrm)
3777
3778 !The energy weighted and overlap term. ONLY with the response density
3779 focc = 2.0_dp
3780 IF (nspins == 2) focc = 1.0_dp
3781 CALL get_qs_env(qs_env, mos=mos)
3782 DO ispin = 1, nspins
3783 CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
3784 CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
3785 p_env%w1(ispin)%matrix, focc, nocc)
3786 END DO
3787 IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, p_env%w1(2)%matrix, 1.0_dp, 1.0_dp)
3788
3789 !Add to it the SCF W matrix, except if EXX (because taken care of by HF response)
3790 IF (.NOT. do_exx) THEN
3791 CALL compute_matrix_w(qs_env, calc_forces=.true.)
3792 CALL get_qs_env(qs_env, matrix_w=matrix_w)
3793 CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(1)%matrix, 1.0_dp, 1.0_dp)
3794 IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(2)%matrix, 1.0_dp, 1.0_dp)
3795 END IF
3796
3797 NULLIFY (scrm)
3798 CALL build_overlap_matrix(qs_env%ks_env, matrix_s=scrm, &
3799 matrix_name="OVERLAP MATRIX", &
3800 basis_type_a="ORB", basis_type_b="ORB", &
3801 sab_nl=sab_orb, calculate_forces=.true., &
3802 matrix_p=p_env%w1(1)%matrix)
3803
3804 IF (.NOT. do_exx) THEN
3805 CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(1)%matrix, 1.0_dp, -1.0_dp)
3806 IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(2)%matrix, 1.0_dp, -1.0_dp)
3807 DO ispin = 1, nspins
3808 CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
3809 END DO
3810 END IF
3811
3812 IF (nspins == 2) CALL dbcsr_add(p_env%w1(1)%matrix, p_env%w1(2)%matrix, 1.0_dp, -1.0_dp)
3813 CALL dbcsr_deallocate_matrix_set(scrm)
3814
3815 IF (use_virial) virial%pv_calculate = .false.
3816
3817 !clean-up
3818 CALL auxbas_pw_pool%give_back_pw(vh_rspace)
3819
3820 DO ispin = 1, nspins
3821 CALL auxbas_pw_pool%give_back_pw(vxc_rspace(ispin))
3822 IF (ASSOCIATED(vtau_rspace)) THEN
3823 CALL auxbas_pw_pool%give_back_pw(vtau_rspace(ispin))
3824 END IF
3825 IF (ASSOCIATED(vadmm_rspace)) THEN
3826 CALL auxbas_pw_pool%give_back_pw(vadmm_rspace(ispin))
3827 END IF
3828 END DO
3829 DEALLOCATE (vxc_rspace)
3830 IF (ASSOCIATED(vtau_rspace)) DEALLOCATE (vtau_rspace)
3831 IF (ASSOCIATED(vadmm_rspace)) DEALLOCATE (vadmm_rspace)
3832
3833 CALL timestop(handle)
3834
3835 END SUBROUTINE update_im_time_forces
3836
3837! **************************************************************************************************
3838!> \brief Iteratively builds the matrix Y = sum_k Y_k until convergence, where
3839!> Y_k = 1/k*2^n (A/2^n) Y_k-1 + 1/k!*2^n * PR(n) * (A/2^n)^(k-1)
3840!> n is chosen such that the norm of A is < 1 (and e^A converges fast)
3841!> PR(n) = e^(A/2^n)*PR(n-1) + PR(n-1)*e^(A/2^n), PR(0) = P*R^T
3842!> \param Y ...
3843!> \param A ...
3844!> \param P ...
3845!> \param R ...
3846!> \param filter_eps ...
3847! **************************************************************************************************
3848 SUBROUTINE build_y_matrix(Y, A, P, R, filter_eps)
3849
3850 TYPE(dbcsr_type), INTENT(OUT) :: y
3851 TYPE(dbcsr_type), INTENT(INOUT) :: a, p, r
3852 REAL(dp), INTENT(IN) :: filter_eps
3853
3854 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_Y_matrix'
3855
3856 INTEGER :: handle, k, n
3857 REAL(dp) :: norm_scalar, threshold
3858 TYPE(dbcsr_type) :: a2n, exp_a2n, prn, work, work2, yk
3859
3860 CALL timeset(routinen, handle)
3861
3862 threshold = 1.0e-16_dp
3863
3864 !Find n such that norm(A) < 1 and we insure convergence of the exponential
3865 norm_scalar = dbcsr_frobenius_norm(a)
3866
3867 !checked: result invariant with value of n
3868 n = 1
3869 DO
3870 IF ((norm_scalar/2.0_dp**n) < 1.0_dp) EXIT
3871 n = n + 1
3872 END DO
3873
3874 !Calculate PR(n) recursively
3875 CALL dbcsr_create(prn, template=a, matrix_type=dbcsr_type_no_symmetry)
3876 CALL dbcsr_create(work, template=a, matrix_type=dbcsr_type_no_symmetry)
3877 CALL dbcsr_multiply('N', 'N', 1.0_dp, p, r, 0.0_dp, work, filter_eps=filter_eps)
3878 CALL dbcsr_create(exp_a2n, template=a, matrix_type=dbcsr_type_no_symmetry)
3879
3880 DO k = 1, n
3881 CALL matrix_exponential(exp_a2n, a, 1.0_dp, 0.5_dp**k, threshold)
3882 CALL dbcsr_multiply('N', 'N', 1.0_dp, exp_a2n, work, 0.0_dp, prn, filter_eps=filter_eps)
3883 CALL dbcsr_multiply('N', 'N', 1.0_dp, work, exp_a2n, 1.0_dp, prn, filter_eps=filter_eps)
3884 CALL dbcsr_copy(work, prn)
3885 END DO
3886 CALL dbcsr_release(exp_a2n)
3887
3888 !Calculate Y iteratively, until convergence
3889 CALL dbcsr_create(a2n, template=a, matrix_type=dbcsr_type_no_symmetry)
3890 CALL dbcsr_copy(a2n, a)
3891 CALL dbcsr_scale(a2n, 0.5_dp**n)
3892 CALL dbcsr_create(y, template=a, matrix_type=dbcsr_type_no_symmetry)
3893 CALL dbcsr_create(yk, template=a, matrix_type=dbcsr_type_no_symmetry)
3894 CALL dbcsr_create(work2, template=a, matrix_type=dbcsr_type_no_symmetry)
3895
3896 !k=1
3897 CALL dbcsr_scale(prn, 0.5_dp**n)
3898 CALL dbcsr_copy(work, prn)
3899 CALL dbcsr_copy(work2, prn)
3900 CALL dbcsr_add(y, prn, 1.0_dp, 1.0_dp)
3901
3902 k = 1
3903 DO
3904 k = k + 1
3905 CALL dbcsr_multiply('N', 'N', 1.0_dp/real(k, dp), a2n, work, 0.0_dp, yk, filter_eps=filter_eps)
3906 CALL dbcsr_multiply('N', 'N', 1.0_dp/real(k, dp), work2, a2n, 0.0_dp, prn, filter_eps=filter_eps)
3907
3908 CALL dbcsr_add(yk, prn, 1.0_dp, 1.0_dp)
3909 CALL dbcsr_add(y, yk, 1.0_dp, 1.0_dp)
3910
3911 IF (dbcsr_frobenius_norm(yk) < threshold) EXIT
3912 CALL dbcsr_copy(work, yk)
3913 CALL dbcsr_copy(work2, prn)
3914 END DO
3915
3916 CALL dbcsr_release(work)
3917 CALL dbcsr_release(work2)
3918 CALL dbcsr_release(prn)
3919 CALL dbcsr_release(a2n)
3920 CALL dbcsr_release(yk)
3921
3922 CALL timestop(handle)
3923
3924 END SUBROUTINE build_y_matrix
3925
3926! **************************************************************************************************
3927!> \brief Overwrites the "optimal" Laplace quadrature with that of the first step
3928!> \param grid ...
3929!> \param do_laplace ...
3930!> \param do_im_time ...
3931!> \param unit_nr ...
3932!> \param qs_env ...
3933! **************************************************************************************************
3934 SUBROUTINE keep_initial_quad(grid, do_laplace, do_im_time, unit_nr, qs_env)
3935
3936 TYPE(time_frequency_grid_type), INTENT(INOUT) :: grid
3937 LOGICAL, INTENT(IN) :: do_laplace, do_im_time
3938 INTEGER, INTENT(IN) :: unit_nr
3939 TYPE(qs_environment_type), POINTER :: qs_env
3940
3941 INTEGER :: jquad
3942
3943 IF (do_laplace .OR. do_im_time) THEN
3944 IF (.NOT. ALLOCATED(qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%imaginary_time)) THEN
3945 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%imaginary_time, source=grid%imaginary_time)
3946 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%time_weights_at_zero_frequency, &
3947 source=grid%time_weights_at_zero_frequency)
3948 ELSE
3949 !If weights already stored, we overwrite the new ones
3950 grid%imaginary_time(:) = qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%imaginary_time(:)
3951 grid%time_weights_at_zero_frequency(:) = &
3952 qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%time_weights_at_zero_frequency(:)
3953 END IF
3954 END IF
3955 IF (.NOT. do_laplace) THEN
3956 IF (.NOT. ALLOCATED(qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%frequency)) THEN
3957 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%frequency, source=grid%frequency)
3958 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%frequency_weights, &
3959 source=grid%frequency_weights)
3960 IF (do_im_time) THEN
3961 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%cosine_time_to_frequency, &
3962 source=grid%cosine_time_to_frequency)
3963 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%cosine_frequency_to_time, &
3964 source=grid%cosine_frequency_to_time)
3965 END IF
3966 ELSE
3967 grid%frequency(:) = qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%frequency(:)
3968 grid%frequency_weights(:) = qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%frequency_weights(:)
3969 IF (do_im_time) THEN
3970 grid%cosine_time_to_frequency(:, :) = &
3971 qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%cosine_time_to_frequency(:, :)
3972 grid%cosine_frequency_to_time(:, :) = &
3973 qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%cosine_frequency_to_time(:, :)
3974 END IF
3975 END IF
3976 END IF
3977 IF (unit_nr > 0) THEN
3978 !Printing order same as in mp2_grids.F for consistency
3979 IF (ALLOCATED(qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%frequency) .AND. (.NOT. do_laplace)) THEN
3980 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
3981 "MINIMAX_INFO| Number of integration points:", SIZE(grid%frequency)
3982 WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") &
3983 "MINIMAX_INFO| Minimax params (freq grid, scaled):", "Weights", "Abscissas"
3984 DO jquad = 1, SIZE(grid%frequency)
3985 WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") &
3986 grid%frequency_weights(jquad), grid%frequency(jquad)
3987 END DO
3988 CALL m_flush(unit_nr)
3989 END IF
3990 IF (ALLOCATED(qs_env%mp2_env%ri_rpa_im_time%time_frequency_grid%imaginary_time)) THEN
3991 WRITE (unit=unit_nr, fmt="(T3,A,T75,i6)") &
3992 "MINIMAX_INFO| Number of integration points:", SIZE(grid%imaginary_time)
3993 WRITE (unit=unit_nr, fmt="(T3,A,T54,A,T72,A)") &
3994 "MINIMAX_INFO| Minimax params (time grid, scaled):", "Weights", "Abscissas"
3995 DO jquad = 1, SIZE(grid%imaginary_time)
3996 WRITE (unit=unit_nr, fmt="(T41,F20.10,F20.10)") &
3997 grid%time_weights_at_zero_frequency(jquad), grid%imaginary_time(jquad)
3998 END DO
3999 CALL m_flush(unit_nr)
4000 END IF
4001 END IF
4002
4003 END SUBROUTINE keep_initial_quad
4004
4005END 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:82
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:126
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:105
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:96
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:104
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F: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:3163
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:216
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:336
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:129
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:940
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:344
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:807
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:240
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:536
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:167
dbt_api
This is the start of a dbt_api, all publically needed functions are exported here....
Definition dbt_api.F:17
distribution_2d_types
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
Definition distribution_2d_types.F:15
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:115
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:2937
hfx_types::dealloc_containers
subroutine, public dealloc_containers(data, memory_usage)
...
Definition hfx_types.F:2905
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:1338
input_constants::ri_rpa_method_gpw
integer, parameter, public ri_rpa_method_gpw
Definition input_constants.F:1130
input_constants::do_admm_aux_exch_func_none
integer, parameter, public do_admm_aux_exch_func_none
Definition input_constants.F:922
input_constants::do_potential_id
integer, parameter, public do_potential_id
Definition input_constants.F:856
input_constants::do_eri_gpw
integer, parameter, public do_eri_gpw
Definition input_constants.F:1338
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:127
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:144
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:2169
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:1713
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:2924
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:544
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:1134
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:454
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:1530
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:440
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:3535
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:3626
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:2187
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:3935
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:76
rpa_im_time::propagator_sector_occupied
integer, parameter, public propagator_sector_occupied
Definition rpa_im_time.F:75
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:666
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:786
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:116
distribution_2d_types::distribution_2d_type
distributes pairs on a 2d grid of processors
Definition distribution_2d_types.F:62
hfx_types::block_ind_type
Definition hfx_types.F:322
hfx_types::hfx_compression_type
Definition hfx_types.F:306
hfx_types::hfx_type
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F:472
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:751
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:774
mp2_types::mp2_type
Definition mp2_types.F:315
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:223
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:30
virial_types::virial_type
Definition virial_types.F:26