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