Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Auxiliary routines needed for RPA-AXK
10 : !> given blacs_env to another
11 : !> \par History
12 : !> 09.2016 created [Vladimir Rybkin]
13 : !> 03.2019 Renamed [Frederick Stein]
14 : !> 03.2019 Moved Functions from rpa_ri_gpw.F [Frederick Stein]
15 : !> \author Vladimir Rybkin
16 : ! **************************************************************************************************
17 : MODULE rpa_axk
18 : USE atomic_kind_types, ONLY: atomic_kind_type
19 : USE cell_types, ONLY: cell_type
20 : USE cp_control_types, ONLY: dft_control_type
21 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
22 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
23 : cp_fm_scale
24 : USE cp_fm_diag, ONLY: choose_eigv_solver
25 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
26 : cp_fm_struct_release,&
27 : cp_fm_struct_type
28 : USE cp_fm_types, ONLY: cp_fm_create,&
29 : cp_fm_get_info,&
30 : cp_fm_release,&
31 : cp_fm_set_all,&
32 : cp_fm_to_fm,&
33 : cp_fm_to_fm_submat_general,&
34 : cp_fm_type
35 : USE dbcsr_api, ONLY: &
36 : dbcsr_copy, dbcsr_create, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
37 : dbcsr_set, dbcsr_trace, dbcsr_type, dbcsr_type_no_symmetry
38 : USE hfx_energy_potential, ONLY: integrate_four_center
39 : USE hfx_ri, ONLY: hfx_ri_update_ks
40 : USE hfx_types, ONLY: hfx_create,&
41 : hfx_release,&
42 : hfx_type
43 : USE input_section_types, ONLY: section_vals_get,&
44 : section_vals_get_subs_vals,&
45 : section_vals_type
46 : USE kinds, ONLY: dp
47 : USE message_passing, ONLY: mp_para_env_type
48 : USE mp2_types, ONLY: mp2_type
49 : USE parallel_gemm_api, ONLY: parallel_gemm
50 : USE particle_types, ONLY: particle_type
51 : USE qs_energy_types, ONLY: qs_energy_type
52 : USE qs_environment_types, ONLY: get_qs_env,&
53 : qs_environment_type
54 : USE qs_kind_types, ONLY: qs_kind_type
55 : USE qs_subsys_types, ONLY: qs_subsys_get,&
56 : qs_subsys_type
57 : USE rpa_communication, ONLY: gamma_fm_to_dbcsr
58 : USE rpa_util, ONLY: calc_fm_mat_S_rpa,&
59 : remove_scaling_factor_rpa
60 : USE scf_control_types, ONLY: scf_control_type
61 : USE util, ONLY: get_limit
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 :
66 : PRIVATE
67 :
68 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_axk'
69 :
70 : PUBLIC :: compute_axk_ener
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief Main driver for RPA-AXK energies
76 : !> \param qs_env ...
77 : !> \param fm_mat_Q ...
78 : !> \param fm_mat_Q_gemm ...
79 : !> \param dimen_RI ...
80 : !> \param dimen_ia ...
81 : !> \param para_env_sub ...
82 : !> \param para_env_RPA ...
83 : !> \param eig ...
84 : !> \param fm_mat_S ...
85 : !> \param homo ...
86 : !> \param virtual ...
87 : !> \param omega ...
88 : !> \param mp2_env ...
89 : !> \param mat_munu ...
90 : !> \param unit_nr ...
91 : !> \param e_axk_corr ... AXK energy correctrion for a quadrature point
92 : !> \author Vladimir Rybkin, 07/2016
93 : ! **************************************************************************************************
94 14 : SUBROUTINE compute_axk_ener(qs_env, fm_mat_Q, fm_mat_Q_gemm, dimen_RI, dimen_ia, &
95 : para_env_sub, para_env_RPA, &
96 14 : eig, fm_mat_S, homo, virtual, omega, &
97 : mp2_env, mat_munu, unit_nr, e_axk_corr)
98 : TYPE(qs_environment_type), POINTER :: qs_env
99 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_Q, fm_mat_Q_gemm
100 : INTEGER, INTENT(IN) :: dimen_RI, dimen_ia
101 : TYPE(mp_para_env_type), POINTER :: para_env_sub, para_env_RPA
102 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: Eig
103 : TYPE(cp_fm_type), INTENT(IN) :: fm_mat_S
104 : INTEGER, INTENT(IN) :: homo, virtual
105 : REAL(KIND=dp), INTENT(IN) :: omega
106 : TYPE(mp2_type) :: mp2_env
107 : TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu
108 : INTEGER, INTENT(IN) :: unit_nr
109 : REAL(KIND=dp), INTENT(INOUT) :: e_axk_corr
110 :
111 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_axk_ener'
112 : REAL(KIND=dp), PARAMETER :: thresh = 0.0000001_dp
113 :
114 : INTEGER :: color_sub, handle, iib, iitmp(2), kkb, L_counter, my_group_L_end, &
115 : my_group_L_size, my_group_L_start, ncol_local, ngroup
116 14 : INTEGER, DIMENSION(:), POINTER :: col_indices
117 : REAL(KIND=dp) :: eps_filter, trace_corr
118 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval
119 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
120 : TYPE(cp_fm_type) :: fm_mat_Gamma_3, fm_mat_Q_tmp, &
121 : fm_mat_R_half, fm_mat_R_half_gemm, &
122 : fm_mat_U
123 14 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_Gamma_3, dbcsr_Gamma_inu_P, &
124 14 : dbcsr_Gamma_munu_P
125 : TYPE(dbcsr_type), POINTER :: mo_coeff_o, mo_coeff_v
126 :
127 14 : CALL timeset(routineN, handle)
128 :
129 : ! Eigenvalues
130 42 : ALLOCATE (eigenval(dimen_RI))
131 926 : eigenval = 0.0_dp
132 : ! create the R_half and U matrices with a different blacs env similar to Q
133 : ! and a tmp_Q needed for diagonalization
134 :
135 14 : NULLIFY (fm_struct)
136 :
137 : CALL cp_fm_get_info(matrix=fm_mat_Q, &
138 14 : matrix_struct=fm_struct)
139 14 : CALL cp_fm_create(fm_mat_U, fm_struct, name="fm_mat_U")
140 14 : CALL cp_fm_create(fm_mat_R_half, fm_struct, name="fm_mat_R_half")
141 14 : CALL cp_fm_create(fm_mat_Q_tmp, fm_struct, name="fm_mat_Q_tmp")
142 14 : CALL cp_fm_set_all(matrix=fm_mat_Q_tmp, alpha=0.0_dp)
143 14 : CALL cp_fm_set_all(matrix=fm_mat_U, alpha=0.0_dp)
144 14 : CALL cp_fm_set_all(matrix=fm_mat_R_half, alpha=0.0_dp)
145 :
146 : ! Copy Q to Q_tmp
147 14 : CALL cp_fm_to_fm(fm_mat_Q, fm_mat_Q_tmp)
148 :
149 14 : CALL cp_fm_scale(0.50_dp, fm_mat_Q_tmp)
150 : ! Diagonalize Q
151 14 : CALL choose_eigv_solver(fm_mat_Q_tmp, fm_mat_U, eigenval)
152 :
153 : !Calculate diagonal matrix for R_half
154 :
155 : ! U*diag stored in U, whereas eigenvectors are in fm_mat_Q_tmp
156 : !CALL cp_fm_to_fm(fm_mat_Q_tmp, fm_mat_U)
157 14 : CALL cp_fm_to_fm(fm_mat_U, fm_mat_Q_tmp)
158 :
159 : ! Manipulate eigenvalues to get diagonal matrix
160 926 : DO iib = 1, dimen_RI
161 926 : IF (ABS(eigenval(iib)) .GE. thresh) THEN
162 : eigenval(iib) = &
163 : SQRT((1.0_dp/(eigenval(iib)**2))*LOG(1.0_dp + eigenval(iib)) &
164 738 : - 1.0_dp/(eigenval(iib)*(eigenval(iib) + 1.0_dp)))
165 : ELSE
166 174 : eigenval(iib) = 0.707_dp
167 : END IF
168 : END DO
169 :
170 14 : CALL cp_fm_column_scale(fm_mat_U, eigenval)
171 :
172 : ! Release memory
173 14 : DEALLOCATE (eigenval)
174 :
175 : ! Get R_half by multiplication
176 : CALL parallel_gemm(transa="N", transb="T", m=dimen_RI, n=dimen_RI, k=dimen_RI, alpha=1.0_dp, &
177 : matrix_a=fm_mat_U, matrix_b=fm_mat_Q_tmp, beta=0.0_dp, &
178 14 : matrix_c=fm_mat_R_half)
179 :
180 : ! get info of fm_mat_S and initialize Gamma_3
181 14 : NULLIFY (fm_struct)
182 14 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=fm_mat_S%matrix_struct, nrow_global=dimen_ia, ncol_global=dimen_RI)
183 14 : CALL cp_fm_create(fm_mat_Gamma_3, fm_struct)
184 14 : CALL cp_fm_struct_release(fm_struct)
185 14 : CALL cp_fm_set_all(matrix=fm_mat_Gamma_3, alpha=0.0_dp)
186 : CALL cp_fm_get_info(matrix=fm_mat_S, &
187 : ncol_local=ncol_local, &
188 14 : col_indices=col_indices)
189 :
190 : ! Update G with a new value of Omega: in practice, it is G*S
191 :
192 : ! Here eig are orbital energies, don't confuse with eigenval, which are eigenvalues of Q!
193 :
194 : ! Scale fm_work_iaP
195 : CALL calc_fm_mat_S_rpa(fm_mat_S, .TRUE., virtual, eig, &
196 14 : homo, omega, 0.0_dp)
197 :
198 : ! Redistribute fm_mat_R_half for "rectangular" multiplication: ia*P P*P
199 14 : CALL cp_fm_create(fm_mat_R_half_gemm, fm_mat_Q_gemm%matrix_struct)
200 14 : CALL cp_fm_set_all(matrix=fm_mat_R_half_gemm, alpha=0.0_dp)
201 :
202 : CALL cp_fm_to_fm_submat_general(fm_mat_R_half, fm_mat_R_half_gemm, dimen_RI, dimen_RI, 1, 1, 1, 1, &
203 14 : fm_mat_R_half%matrix_struct%context)
204 :
205 : ! Calculate Gamma_3: Gamma_3 = G*S*R^(1/2) = G*S*R^(1/2) )
206 : CALL parallel_gemm(transa="T", transb="N", m=dimen_ia, n=dimen_RI, k=dimen_RI, alpha=1.0_dp, &
207 : matrix_a=fm_mat_S, matrix_b=fm_mat_R_half_gemm, beta=0.0_dp, &
208 14 : matrix_c=fm_mat_Gamma_3)
209 :
210 : ! Remove extra factor from S after the multiplication
211 14 : CALL remove_scaling_factor_rpa(fm_mat_S, virtual, eig, homo, omega)
212 :
213 : ! Release full matrix stuff
214 14 : CALL cp_fm_release(fm_mat_Q_tmp)
215 14 : CALL cp_fm_release(fm_mat_U)
216 14 : CALL cp_fm_release(fm_mat_R_half)
217 14 : CALL cp_fm_release(fm_mat_R_half_gemm)
218 :
219 : ! Retrieve mo coefficients in dbcsr format
220 14 : NULLIFY (mo_coeff_o, mo_coeff_v)
221 14 : mo_coeff_o => mp2_env%ri_rpa%mo_coeff_o
222 14 : mo_coeff_v => mp2_env%ri_rpa%mo_coeff_v
223 :
224 : ! Get aux sizes
225 14 : ngroup = para_env_RPA%num_pe/para_env_sub%num_pe
226 :
227 14 : color_sub = para_env_RPA%mepos/para_env_sub%num_pe
228 :
229 14 : iitmp = get_limit(dimen_RI, ngroup, color_sub)
230 14 : my_group_L_start = iitmp(1)
231 14 : my_group_L_end = iitmp(2)
232 14 : my_group_L_size = iitmp(2) - iitmp(1) + 1
233 :
234 : ! Copy Gamma_ia_P^3 to dbcsr matrix set
235 : CALL gamma_fm_to_dbcsr(fm_mat_Gamma_3, dbcsr_Gamma_3, para_env_RPA, para_env_sub, &
236 : homo, virtual, mo_coeff_o, ngroup, my_group_L_start, &
237 14 : my_group_L_end, my_group_L_size)
238 :
239 : ! Create more dbcsr matrices
240 :
241 14 : NULLIFY (dbcsr_Gamma_inu_P)
242 : !CALL dbcsr_allocate_matrix_set(dbcsr_Gamma_inu_P, ncol_local)
243 14 : CALL dbcsr_allocate_matrix_set(dbcsr_Gamma_inu_P, my_group_L_size)
244 14 : NULLIFY (dbcsr_Gamma_munu_P)
245 : !CALL dbcsr_allocate_matrix_set(dbcsr_Gamma_munu_P, ncol_local)
246 14 : CALL dbcsr_allocate_matrix_set(dbcsr_Gamma_munu_P, my_group_L_size)
247 14 : eps_filter = mp2_env%mp2_gpw%eps_filter
248 :
249 14 : L_counter = 0
250 470 : DO kkb = my_group_L_start, my_group_L_end
251 456 : L_counter = L_counter + 1
252 : ! One-index transformed Gamma_3
253 456 : ALLOCATE (dbcsr_Gamma_inu_P(L_counter)%matrix)
254 456 : CALL dbcsr_init_p(dbcsr_Gamma_inu_P(L_counter)%matrix)
255 456 : CALL dbcsr_create(dbcsr_Gamma_inu_P(L_counter)%matrix, template=mo_coeff_o)
256 456 : CALL dbcsr_copy(dbcsr_Gamma_inu_P(L_counter)%matrix, mo_coeff_o)
257 456 : CALL dbcsr_set(dbcsr_Gamma_inu_P(L_counter)%matrix, 0.0_dp)
258 : ! Init Gamma_3 in AO basis
259 456 : ALLOCATE (dbcsr_Gamma_munu_P(L_counter)%matrix)
260 456 : CALL dbcsr_init_p(dbcsr_Gamma_munu_P(L_counter)%matrix)
261 : CALL dbcsr_create(dbcsr_Gamma_munu_P(L_counter)%matrix, template=mat_munu%matrix, &
262 456 : matrix_type=dbcsr_type_no_symmetry)
263 456 : CALL dbcsr_copy(dbcsr_Gamma_munu_P(L_counter)%matrix, mat_munu%matrix)
264 470 : CALL dbcsr_set(dbcsr_Gamma_munu_P(L_counter)%matrix, 0.0_dp)
265 : END DO
266 :
267 : !! Loup over auxiliary basis functions: multiplication
268 : L_counter = 0
269 470 : DO kkb = my_group_L_start, my_group_L_end
270 456 : L_counter = L_counter + 1
271 : ! Do dbcsr multiplication: transform the virtual index
272 : CALL dbcsr_multiply("N", "T", 1.0_dp, mo_coeff_v, dbcsr_Gamma_3(L_counter)%matrix, &
273 456 : 0.0_dp, dbcsr_Gamma_inu_P(L_counter)%matrix, filter_eps=eps_filter)
274 :
275 : !Do dbcsr multiplication: transform the occupied index
276 : CALL dbcsr_multiply("N", "T", 1.0_dp, dbcsr_Gamma_inu_P(L_counter)%matrix, mo_coeff_o, &
277 456 : 0.0_dp, dbcsr_Gamma_munu_P(L_counter)%matrix, filter_eps=eps_filter)
278 : !
279 470 : CALL dbcsr_trace(dbcsr_Gamma_munu_P(L_counter)%matrix, trace_corr)
280 : END DO
281 :
282 : ! Gamma_3 not needed anymore
283 : L_counter = 0
284 470 : DO kkb = my_group_L_start, my_group_L_end
285 456 : L_counter = L_counter + 1
286 456 : CALL dbcsr_release(dbcsr_Gamma_3(L_counter)%matrix)
287 470 : DEALLOCATE (dbcsr_Gamma_3(L_counter)%matrix)
288 : END DO
289 14 : DEALLOCATE (dbcsr_Gamma_3)
290 :
291 : ! Contract DM with exchange integrals
292 : !CALL integrate_exchange(qs_env, dbcsr_Gamma_munu_P, mat_munu, para_env_sub, ncol_local, eps_filter, e_axk_corr)
293 : CALL integrate_exchange(qs_env, dbcsr_Gamma_munu_P, mat_munu, para_env_sub, my_group_L_size, eps_filter, e_axk_corr, &
294 14 : my_group_L_start, my_group_L_end)
295 :
296 : !CALL para_env_RPA%sum(e_axk_corr)
297 :
298 : ! Print AXK correlation energy to the file
299 21 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T68,F25.14,A4)') 'AXK correlation energy for a quadrature point:', &
300 14 : e_axk_corr, ' a.u.'
301 :
302 : L_counter = 0
303 470 : DO kkb = my_group_L_start, my_group_L_end
304 456 : L_counter = L_counter + 1
305 456 : CALL dbcsr_release(dbcsr_Gamma_inu_P(L_counter)%matrix)
306 456 : CALL dbcsr_release(dbcsr_Gamma_munu_P(L_counter)%matrix)
307 456 : DEALLOCATE (dbcsr_Gamma_inu_P(L_counter)%matrix)
308 470 : DEALLOCATE (dbcsr_Gamma_munu_P(L_counter)%matrix)
309 : END DO
310 14 : DEALLOCATE (dbcsr_Gamma_inu_P)
311 14 : DEALLOCATE (dbcsr_Gamma_munu_P)
312 :
313 14 : CALL timestop(handle)
314 :
315 28 : END SUBROUTINE compute_axk_ener
316 :
317 : ! **************************************************************************************************
318 : !> \brief Contract RPA-AXK density matrix with HF exchange integrals and evaluate the correction
319 : !> \param qs_env ...
320 : !> \param dbcsr_Gamma_munu_P ... AXK density matrix in AO basis to be contracted
321 : !> \param mat_munu ...
322 : !> \param para_env_sub ...
323 : !> \param P_stack_size ...
324 : !> \param eps_filter ...
325 : !> \param axk_corr ... The AXK energy correction
326 : !> \param my_group_L_start ...
327 : !> \param my_group_L_end ...
328 : !> \author Vladimir Rybkin, 08/2016
329 : ! **************************************************************************************************
330 28 : SUBROUTINE integrate_exchange(qs_env, dbcsr_Gamma_munu_P, mat_munu, para_env_sub, P_stack_size, &
331 : eps_filter, axk_corr, &
332 : my_group_L_start, my_group_L_end)
333 : TYPE(qs_environment_type), POINTER :: qs_env
334 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dbcsr_Gamma_munu_P
335 : TYPE(dbcsr_p_type), INTENT(IN) :: mat_munu
336 : TYPE(mp_para_env_type), POINTER :: para_env_sub
337 : INTEGER, INTENT(INOUT) :: P_stack_size
338 : REAL(KIND=dp), INTENT(IN) :: eps_filter
339 : REAL(KIND=dp), INTENT(OUT) :: axk_corr
340 : INTEGER, INTENT(IN) :: my_group_L_start, my_group_L_end
341 :
342 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_exchange'
343 :
344 : INTEGER :: aux, handle, irep, kkb, n_rep_hf, ns
345 : LOGICAL :: my_recalc_hfx_integrals
346 : REAL(KIND=dp) :: e_axk_P, ehfx
347 14 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_work_ao
348 14 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_2d, rho_ao_2d
349 14 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
350 : TYPE(qs_energy_type), POINTER :: energy
351 : TYPE(section_vals_type), POINTER :: hfx_sections
352 :
353 14 : CALL timeset(routineN, handle)
354 :
355 : ! Get qs environment
356 14 : NULLIFY (energy)
357 : CALL get_qs_env(qs_env, &
358 14 : energy=energy)
359 :
360 : ! hfx section
361 14 : CALL hfx_create_subgroup(qs_env, para_env_sub, hfx_sections, x_data, n_rep_hf)
362 :
363 : ! create a working rho environment
364 14 : NULLIFY (rho_work_ao)
365 14 : CALL dbcsr_allocate_matrix_set(rho_work_ao, 1)
366 14 : ALLOCATE (rho_work_ao(1)%matrix)
367 14 : CALL dbcsr_init_p(rho_work_ao(1)%matrix)
368 14 : CALL dbcsr_create(rho_work_ao(1)%matrix, template=mat_munu%matrix)
369 :
370 : ! For the first aux function in the group we recalculate integrals, but only for the first
371 14 : my_recalc_hfx_integrals = .TRUE.
372 :
373 14 : NULLIFY (mat_2d)
374 14 : CALL dbcsr_allocate_matrix_set(mat_2d, 1, 1)
375 14 : ALLOCATE (mat_2d(1, 1)%matrix)
376 14 : CALL dbcsr_init_p(mat_2d(1, 1)%matrix)
377 : CALL dbcsr_create(mat_2d(1, 1)%matrix, template=mat_munu%matrix, &
378 14 : matrix_type=dbcsr_type_no_symmetry)
379 14 : CALL dbcsr_copy(mat_2d(1, 1)%matrix, mat_munu%matrix)
380 :
381 : ! The loop over auxiliary basis functions
382 14 : axk_corr = 0.0_dp
383 : !DO aux = 1, P_stack_size
384 : P_stack_size = P_stack_size
385 14 : aux = 0
386 470 : DO kkb = my_group_L_start, my_group_L_end
387 456 : aux = aux + 1
388 :
389 456 : CALL dbcsr_copy(rho_work_ao(1)%matrix, dbcsr_Gamma_munu_P(aux)%matrix)
390 :
391 912 : DO irep = 1, n_rep_hf
392 456 : ns = SIZE(rho_work_ao)
393 456 : rho_ao_2d(1:ns, 1:1) => rho_work_ao(1:ns)
394 :
395 456 : CALL dbcsr_set(mat_2d(1, 1)%matrix, 0.0_dp)
396 :
397 912 : IF (x_data(irep, 1)%do_hfx_ri) THEN
398 : CALL hfx_ri_update_ks(qs_env, x_data(irep, 1)%ri_data, mat_2d, ehfx, &
399 : rho_ao=rho_ao_2d, geometry_did_change=my_recalc_hfx_integrals, &
400 0 : nspins=ns, hf_fraction=x_data(irep, 1)%general_parameter%fraction)
401 :
402 : ELSE
403 : CALL integrate_four_center(qs_env, x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
404 : para_env_sub, my_recalc_hfx_integrals, irep, .TRUE., &
405 456 : ispin=1)
406 : END IF
407 : END DO
408 :
409 456 : my_recalc_hfx_integrals = .FALSE.
410 : ! One more dbcsr multiplication and trace
411 : CALL dbcsr_multiply("T", "N", 1.0_dp, mat_2d(1, 1)%matrix, rho_work_ao(1)%matrix, &
412 456 : 0.0_dp, dbcsr_Gamma_munu_P(aux)%matrix, filter_eps=eps_filter)
413 456 : CALL dbcsr_trace(dbcsr_Gamma_munu_P(aux)%matrix, e_axk_p)
414 470 : axk_corr = axk_corr + e_axk_P
415 : END DO
416 :
417 14 : CALL dbcsr_release(mat_2d(1, 1)%matrix)
418 : ! release rho stuff
419 14 : CALL dbcsr_release(mat_2d(1, 1)%matrix)
420 14 : DEALLOCATE (mat_2d(1, 1)%matrix)
421 14 : DEALLOCATE (mat_2d)
422 14 : CALL dbcsr_release(rho_work_ao(1)%matrix)
423 14 : DEALLOCATE (rho_work_ao(1)%matrix)
424 14 : DEALLOCATE (rho_work_ao)
425 14 : CALL hfx_release(x_data)
426 :
427 14 : CALL timestop(handle)
428 :
429 14 : END SUBROUTINE integrate_exchange
430 :
431 : ! **************************************************************************************************
432 : !> \brief ... Initializes x_data on a subgroup
433 : !> \param qs_env ...
434 : !> \param para_env_sub ...
435 : !> \param hfx_section ...
436 : !> \param x_data ...
437 : !> \param n_rep_hf ...
438 : !> \author Vladimir Rybkin
439 : ! **************************************************************************************************
440 28 : SUBROUTINE hfx_create_subgroup(qs_env, para_env_sub, hfx_section, x_data, n_rep_hf)
441 : TYPE(qs_environment_type), POINTER :: qs_env
442 : TYPE(mp_para_env_type), POINTER :: para_env_sub
443 : TYPE(section_vals_type), POINTER :: hfx_section
444 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
445 : INTEGER, INTENT(OUT) :: n_rep_hf
446 :
447 : CHARACTER(LEN=*), PARAMETER :: routineN = 'hfx_create_subgroup'
448 :
449 : INTEGER :: handle, nelectron_total
450 : LOGICAL :: do_hfx
451 14 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
452 : TYPE(cell_type), POINTER :: my_cell
453 : TYPE(dft_control_type), POINTER :: dft_control
454 14 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
455 14 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
456 : TYPE(qs_subsys_type), POINTER :: subsys
457 : TYPE(scf_control_type), POINTER :: scf_control
458 : TYPE(section_vals_type), POINTER :: input
459 :
460 14 : CALL timeset(routineN, handle)
461 :
462 14 : NULLIFY (my_cell, atomic_kind_set, particle_set, dft_control, x_data, qs_kind_set, scf_control)
463 :
464 : CALL get_qs_env(qs_env, &
465 : subsys=subsys, &
466 : input=input, &
467 : scf_control=scf_control, &
468 14 : nelectron_total=nelectron_total)
469 :
470 : CALL qs_subsys_get(subsys, &
471 : cell=my_cell, &
472 : atomic_kind_set=atomic_kind_set, &
473 : qs_kind_set=qs_kind_set, &
474 14 : particle_set=particle_set)
475 :
476 : do_hfx = .TRUE.
477 14 : hfx_section => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
478 : !hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
479 14 : CALL section_vals_get(hfx_section, explicit=do_hfx, n_repetition=n_rep_hf)
480 14 : CALL get_qs_env(qs_env, dft_control=dft_control)
481 :
482 14 : IF (do_hfx) THEN
483 : ! Retrieve particle_set and atomic_kind_set
484 : CALL hfx_create(x_data, para_env_sub, hfx_section, atomic_kind_set, &
485 : qs_kind_set, particle_set, dft_control, my_cell, orb_basis='ORB', &
486 14 : nelectron_total=nelectron_total)
487 : END IF
488 :
489 14 : CALL timestop(handle)
490 :
491 14 : END SUBROUTINE hfx_create_subgroup
492 :
493 : END MODULE rpa_axk
|
