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