Line data Source code
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 to compute singles correction to RPA (RSE)
10 : !> \par History
11 : !> 08.2019 created [Vladimir Rybkin]
12 : !> \author Vladimir Rybkin
13 : ! **************************************************************************************************
14 : MODULE rpa_rse
15 :
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
19 : dbcsr_create,&
20 : dbcsr_init_p,&
21 : dbcsr_p_type,&
22 : dbcsr_release,&
23 : dbcsr_scale,&
24 : dbcsr_set,&
25 : dbcsr_type_symmetric
26 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
27 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 : copy_fm_to_dbcsr,&
29 : dbcsr_allocate_matrix_set
30 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
31 : USE cp_fm_diag, ONLY: choose_eigv_solver
32 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 : cp_fm_struct_release,&
34 : cp_fm_struct_type
35 : USE cp_fm_types, ONLY: cp_fm_create,&
36 : cp_fm_get_diag,&
37 : cp_fm_get_info,&
38 : cp_fm_release,&
39 : cp_fm_set_all,&
40 : cp_fm_to_fm_submat,&
41 : cp_fm_type
42 : USE hfx_energy_potential, ONLY: integrate_four_center
43 : USE hfx_exx, ONLY: exx_post_hfx,&
44 : exx_pre_hfx
45 : USE hfx_ri, ONLY: hfx_ri_update_ks
46 : USE input_section_types, ONLY: section_vals_get,&
47 : section_vals_get_subs_vals,&
48 : section_vals_type,&
49 : section_vals_val_get
50 : USE kinds, ONLY: dp
51 : USE message_passing, ONLY: mp_para_env_type
52 : USE mp2_types, ONLY: mp2_type
53 : USE parallel_gemm_api, ONLY: parallel_gemm
54 : USE pw_types, ONLY: pw_r3d_rs_type
55 : USE qs_energy_types, ONLY: qs_energy_type
56 : USE qs_environment_types, ONLY: get_qs_env,&
57 : qs_environment_type
58 : USE qs_ks_types, ONLY: qs_ks_env_type
59 : USE qs_ks_utils, ONLY: compute_matrix_vxc
60 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
61 : USE qs_rho_types, ONLY: qs_rho_get,&
62 : qs_rho_type
63 : USE qs_vxc, ONLY: qs_vxc_create
64 :
65 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
66 :
67 : #include "./base/base_uses.f90"
68 :
69 : IMPLICIT NONE
70 :
71 : PRIVATE
72 :
73 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_rse'
74 :
75 : PUBLIC :: rse_energy
76 :
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief Single excitations energy corrections for RPA
81 : !> \param qs_env ...
82 : !> \param mp2_env ...
83 : !> \param para_env ...
84 : !> \param dft_control ...
85 : !> \param mo_coeff ...
86 : !> \param homo ...
87 : !> \param Eigenval ...
88 : !> \author Vladimir Rybkin, 08/2019
89 : ! **************************************************************************************************
90 6 : SUBROUTINE rse_energy(qs_env, mp2_env, para_env, dft_control, &
91 6 : mo_coeff, homo, Eigenval)
92 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
93 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
94 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
95 : TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
96 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
97 : INTEGER, DIMENSION(:), INTENT(IN) :: homo
98 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Eigenval
99 :
100 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rse_energy'
101 :
102 : INTEGER :: handle, i_global, iiB, ispin, j_global, &
103 : jjB, n_rep_hf, nao, ncol_local, nmo, &
104 : nrow_local, nspins
105 6 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
106 : LOGICAL :: do_hfx, hfx_treat_lsd_in_core
107 : REAL(KIND=dp) :: coeff, corr, rse_corr
108 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag_diff
109 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
110 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
111 : TYPE(cp_fm_type) :: fm_ao, fm_ao_mo
112 6 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_P_mu_nu, fm_X_mo, fm_XC_mo
113 6 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_mu_nu, matrix_s, rho_ao
114 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
115 6 : POINTER :: sab_orb
116 : TYPE(qs_energy_type), POINTER :: energy
117 : TYPE(qs_rho_type), POINTER :: rho
118 : TYPE(section_vals_type), POINTER :: hfx_sections, input
119 :
120 6 : CALL timeset(routineN, handle)
121 :
122 6 : nspins = dft_control%nspins
123 :
124 : ! Pick the diagonal terms
125 : CALL cp_fm_get_info(matrix=mo_coeff(1), &
126 : nrow_local=nrow_local, &
127 : ncol_local=ncol_local, &
128 : row_indices=row_indices, &
129 : col_indices=col_indices, &
130 : nrow_global=nao, &
131 6 : ncol_global=nmo)
132 :
133 : ! start collecting stuff
134 6 : NULLIFY (input, matrix_s, blacs_env, rho, energy, sab_orb)
135 : CALL get_qs_env(qs_env, &
136 : input=input, &
137 : matrix_s=matrix_s, &
138 : blacs_env=blacs_env, &
139 : rho=rho, &
140 : energy=energy, &
141 6 : sab_orb=sab_orb)
142 :
143 6 : CALL qs_rho_get(rho, rho_ao=rho_ao)
144 :
145 : ! hfx section
146 6 : NULLIFY (hfx_sections)
147 6 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
148 6 : CALL section_vals_get(hfx_sections, explicit=do_hfx, n_repetition=n_rep_hf)
149 6 : IF (do_hfx) THEN
150 : CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
151 6 : i_rep_section=1)
152 : END IF
153 :
154 : ! create work array
155 6 : NULLIFY (mat_mu_nu)
156 6 : CALL dbcsr_allocate_matrix_set(mat_mu_nu, nspins)
157 16 : DO ispin = 1, nspins
158 10 : ALLOCATE (mat_mu_nu(ispin)%matrix)
159 : CALL dbcsr_create(matrix=mat_mu_nu(ispin)%matrix, template=matrix_s(1)%matrix, name="T_mu_nu", &
160 10 : matrix_type=dbcsr_type_symmetric)
161 10 : CALL cp_dbcsr_alloc_block_from_nbl(mat_mu_nu(ispin)%matrix, sab_orb)
162 16 : CALL dbcsr_set(mat_mu_nu(ispin)%matrix, 0.0_dp)
163 : END DO
164 :
165 : ! Dense (full) matrices
166 28 : ALLOCATE (fm_P_mu_nu(nspins))
167 6 : NULLIFY (fm_struct_tmp)
168 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
169 6 : nrow_global=nao, ncol_global=nao)
170 16 : DO ispin = 1, nspins
171 10 : CALL cp_fm_create(fm_P_mu_nu(ispin), fm_struct_tmp, name="P_mu_nu")
172 16 : CALL cp_fm_set_all(fm_P_mu_nu(ispin), 0.0_dp)
173 : END DO
174 6 : CALL cp_fm_create(fm_ao, fm_struct_tmp, name="f_ao")
175 6 : CALL cp_fm_struct_release(fm_struct_tmp)
176 6 : CALL cp_fm_set_all(fm_ao, 0.0_dp)
177 6 : CALL cp_fm_struct_release(fm_struct_tmp)
178 :
179 6 : NULLIFY (fm_struct_tmp)
180 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
181 6 : nrow_global=nmo, ncol_global=nmo)
182 50 : ALLOCATE (fm_X_mo(nspins), fm_XC_mo(nspins))
183 16 : DO ispin = 1, nspins
184 10 : CALL cp_fm_create(fm_X_mo(ispin), fm_struct_tmp, name="f_X_mo")
185 10 : CALL cp_fm_create(fm_XC_mo(ispin), fm_struct_tmp, name="f_XC_mo")
186 10 : CALL cp_fm_set_all(fm_X_mo(ispin), 0.0_dp)
187 16 : CALL cp_fm_set_all(fm_XC_mo(ispin), 0.0_dp)
188 : END DO
189 6 : CALL cp_fm_struct_release(fm_struct_tmp)
190 :
191 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
192 6 : nrow_global=nmo, ncol_global=nao)
193 6 : CALL cp_fm_create(fm_ao_mo, fm_struct_tmp, name="f_ao_mo")
194 6 : CALL cp_fm_struct_release(fm_struct_tmp)
195 6 : CALL cp_fm_set_all(fm_ao_mo, 0.0_dp)
196 :
197 : !
198 : ! Ready with preparations, do the real staff
199 : !
200 :
201 : ! Obtain density matrix like quantity
202 :
203 6 : coeff = 1.0_dp
204 6 : IF (nspins == 1) coeff = 2.0_dp
205 16 : DO ispin = 1, nspins
206 : CALL parallel_gemm(transa='N', transb='T', m=nao, n=nao, k=homo(ispin), alpha=coeff, &
207 : matrix_a=mo_coeff(ispin), matrix_b=mo_coeff(ispin), &
208 16 : beta=0.0_dp, matrix_c=fm_P_mu_nu(ispin))
209 : END DO
210 :
211 : ! Calculate exact exchange contribution
212 : CALL exchange_contribution(qs_env, para_env, mo_coeff, &
213 : hfx_sections, n_rep_hf, &
214 : rho, mat_mu_nu, fm_P_mu_nu, &
215 6 : fm_ao, fm_X_mo, fm_ao_mo)
216 :
217 : ! Calculate DFT exchange-correlation contribution
218 6 : CALL xc_contribution(qs_env, fm_ao, fm_ao_mo, fm_XC_mo, mo_coeff)
219 :
220 18 : ALLOCATE (diag_diff(nmo))
221 6 : rse_corr = 0.0_dp
222 :
223 16 : DO ispin = 1, nspins
224 : ! Compute the correction matrix: it is stored in fm_X_mo
225 10 : CALL cp_fm_scale_and_add(1.0_dp, fm_X_mo(ispin), -1.0_dp, fm_XC_mo(ispin))
226 :
227 : ! Pick the diagonal terms
228 10 : CALL cp_fm_get_diag(fm_X_mo(ispin), diag_diff)
229 :
230 : ! Compute the correction
231 : CALL cp_fm_get_info(matrix=fm_X_mo(ispin), &
232 : nrow_local=nrow_local, &
233 : ncol_local=ncol_local, &
234 : row_indices=row_indices, &
235 10 : col_indices=col_indices)
236 :
237 10 : corr = 0.0_dp
238 :
239 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
240 : !$OMP REDUCTION(+: corr) &
241 10 : !$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,diag_diff,eigenval,fm_X_mo,homo,ispin)
242 : DO jjB = 1, ncol_local
243 : j_global = col_indices(jjB)
244 : DO iiB = 1, nrow_local
245 : i_global = row_indices(iiB)
246 : IF ((i_global .LE. homo(ispin)) .AND. (j_global .GT. homo(ispin))) THEN
247 : corr = corr + fm_X_mo(ispin)%local_data(iib, jjb)**2.0_dp/ &
248 : (eigenval(i_global, ispin) - eigenval(j_global, ispin) - diag_diff(i_global) + diag_diff(j_global))
249 : END IF
250 : END DO
251 : END DO
252 : !$OMP END PARALLEL DO
253 :
254 16 : rse_corr = rse_corr + corr
255 : END DO
256 :
257 6 : CALL para_env%sum(rse_corr)
258 :
259 6 : IF (nspins == 1) rse_corr = rse_corr*2.0_dp
260 :
261 6 : mp2_env%ri_rpa%rse_corr_diag = rse_corr
262 :
263 6 : CALL non_diag_rse(fm_X_mo, eigenval, homo, para_env, blacs_env, rse_corr)
264 :
265 6 : IF (nspins == 1) rse_corr = rse_corr*2.0_dp
266 :
267 6 : mp2_env%ri_rpa%rse_corr = rse_corr
268 :
269 : ! Release staff
270 6 : DEALLOCATE (diag_diff)
271 6 : CALL cp_fm_release(fm_ao)
272 6 : CALL cp_fm_release(fm_ao_mo)
273 6 : CALL cp_fm_release(fm_P_mu_nu)
274 6 : CALL cp_fm_release(fm_X_mo)
275 6 : CALL cp_fm_release(fm_XC_mo)
276 16 : DO ispin = 1, nspins
277 10 : CALL dbcsr_release(mat_mu_nu(ispin)%matrix)
278 16 : DEALLOCATE (mat_mu_nu(ispin)%matrix)
279 : END DO
280 6 : DEALLOCATE (mat_mu_nu)
281 :
282 6 : CALL timestop(handle)
283 :
284 30 : END SUBROUTINE rse_energy
285 :
286 : ! **************************************************************************************************
287 : !> \brief HF exchange occupied-virtual matrix
288 : !> \param qs_env ...
289 : !> \param para_env ...
290 : !> \param mo_coeff ...
291 : !> \param hfx_sections ...
292 : !> \param n_rep_hf ...
293 : !> \param rho_work ...
294 : !> \param mat_mu_nu ...
295 : !> \param fm_P_mu_nu ...
296 : !> \param fm_X_ao ...
297 : !> \param fm_X_mo ...
298 : !> \param fm_X_ao_mo ...
299 : ! **************************************************************************************************
300 6 : SUBROUTINE exchange_contribution(qs_env, para_env, mo_coeff, &
301 : hfx_sections, n_rep_hf, &
302 6 : rho_work, mat_mu_nu, fm_P_mu_nu, &
303 6 : fm_X_ao, fm_X_mo, fm_X_ao_mo)
304 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
305 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
306 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
307 : TYPE(section_vals_type), INTENT(IN), POINTER :: hfx_sections
308 : INTEGER, INTENT(IN) :: n_rep_hf
309 : TYPE(qs_rho_type), INTENT(IN), POINTER :: rho_work
310 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
311 : POINTER :: mat_mu_nu
312 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_P_mu_nu
313 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_X_ao
314 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_X_mo
315 : TYPE(cp_fm_type), INTENT(IN) :: fm_X_ao_mo
316 :
317 : CHARACTER(LEN=*), PARAMETER :: routineN = 'exchange_contribution'
318 :
319 : INTEGER :: handle, irep, is, nao, nmo, ns
320 : LOGICAL :: my_recalc_hfx_integrals
321 : REAL(KIND=dp) :: ehfx
322 6 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: P_mu_nu, rho_work_ao
323 6 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_2d, rho_ao_2d
324 :
325 6 : CALL timeset(routineN, handle)
326 :
327 6 : CALL cp_fm_get_info(mo_coeff(1), nrow_global=nao, ncol_global=nmo)
328 :
329 6 : CALL qs_rho_get(rho_work, rho_ao=rho_work_ao)
330 6 : ns = SIZE(rho_work_ao)
331 6 : NULLIFY (P_mu_nu)
332 6 : CALL dbcsr_allocate_matrix_set(P_mu_nu, ns)
333 16 : DO is = 1, ns
334 10 : CALL dbcsr_init_p(P_mu_nu(is)%matrix)
335 10 : CALL dbcsr_create(P_mu_nu(is)%matrix, template=rho_work_ao(1)%matrix)
336 10 : CALL dbcsr_copy(P_mu_nu(is)%matrix, rho_work_ao(1)%matrix)
337 16 : CALL dbcsr_set(P_mu_nu(is)%matrix, 0.0_dp)
338 : END DO
339 :
340 6 : my_recalc_hfx_integrals = .TRUE.
341 :
342 6 : CALL exx_pre_hfx(hfx_sections, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
343 16 : DO is = 1, ns
344 10 : CALL copy_fm_to_dbcsr(fm_P_mu_nu(is), P_mu_nu(1)%matrix, keep_sparsity=.TRUE.)
345 :
346 10 : CALL dbcsr_set(mat_mu_nu(1)%matrix, 0.0_dp)
347 :
348 10 : IF (qs_env%mp2_env%ri_rpa%x_data(1, 1)%do_hfx_ri) THEN
349 :
350 0 : DO irep = 1, n_rep_hf
351 0 : rho_ao_2d(1:ns, 1:1) => P_mu_nu(1:ns)
352 0 : mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns)
353 : CALL hfx_ri_update_ks(qs_env, qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data, mat_2d, ehfx, &
354 : rho_ao=rho_ao_2d, geometry_did_change=my_recalc_hfx_integrals, nspins=1, &
355 0 : hf_fraction=qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction)
356 :
357 0 : IF (ns == 2) CALL dbcsr_scale(mat_mu_nu(1)%matrix, 2.0_dp)
358 0 : my_recalc_hfx_integrals = .FALSE.
359 : END DO
360 :
361 : ELSE
362 :
363 20 : DO irep = 1, n_rep_hf
364 10 : rho_ao_2d(1:ns, 1:1) => P_mu_nu(1:ns)
365 10 : mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns)
366 : CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
367 : para_env, my_recalc_hfx_integrals, irep, .TRUE., &
368 10 : ispin=1)
369 :
370 20 : my_recalc_hfx_integrals = .FALSE.
371 : END DO
372 : END IF
373 :
374 : ! copy back to fm
375 10 : CALL cp_fm_set_all(fm_X_ao, 0.0_dp)
376 10 : CALL copy_dbcsr_to_fm(matrix=mat_mu_nu(1)%matrix, fm=fm_X_ao)
377 10 : CALL cp_fm_set_all(fm_X_mo(is), 0.0_dp)
378 :
379 : ! First index
380 : CALL parallel_gemm('T', 'N', nmo, nao, nmo, 1.0_dp, &
381 10 : mo_coeff(is), fm_X_ao, 0.0_dp, fm_X_ao_mo)
382 :
383 : ! Second index
384 : CALL parallel_gemm('N', 'N', nmo, nmo, nao, 1.0_dp, &
385 16 : fm_X_ao_mo, mo_coeff(is), 1.0_dp, fm_X_mo(is))
386 :
387 : END DO
388 6 : CALL exx_post_hfx(qs_env, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
389 :
390 : ! Release dbcsr objects
391 16 : DO is = 1, SIZE(P_mu_nu)
392 10 : CALL dbcsr_release(P_mu_nu(is)%matrix)
393 16 : DEALLOCATE (P_mu_nu(is)%matrix)
394 : END DO
395 6 : DEALLOCATE (P_mu_nu)
396 :
397 6 : CALL timestop(handle)
398 :
399 6 : END SUBROUTINE exchange_contribution
400 :
401 : ! **************************************************************************************************
402 : !> \brief Exchange-correlation occupied-virtual matrix
403 : !> \param qs_env ...
404 : !> \param fm_XC_ao ...
405 : !> \param fm_XC_ao_mo ...
406 : !> \param fm_XC_mo ...
407 : !> \param mo_coeff ...
408 : ! **************************************************************************************************
409 6 : SUBROUTINE xc_contribution(qs_env, fm_XC_ao, fm_XC_ao_mo, fm_XC_mo, mo_coeff)
410 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
411 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_XC_ao
412 : TYPE(cp_fm_type), INTENT(IN) :: fm_XC_ao_mo
413 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_XC_mo, mo_coeff
414 :
415 : CHARACTER(LEN=*), PARAMETER :: routineN = 'xc_contribution'
416 :
417 : INTEGER :: handle, i, nao, nmo
418 : REAL(KIND=dp) :: exc
419 6 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc
420 6 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau_rspace, v_rspace
421 : TYPE(qs_ks_env_type), POINTER :: ks_env
422 : TYPE(qs_rho_type), POINTER :: rho
423 : TYPE(section_vals_type), POINTER :: input, xc_section
424 :
425 6 : CALL timeset(routineN, handle)
426 :
427 6 : NULLIFY (matrix_vxc, v_rspace, tau_rspace, input, xc_section, ks_env, &
428 6 : rho)
429 6 : CALL get_qs_env(qs_env, matrix_vxc=matrix_vxc, input=input, ks_env=ks_env, rho=rho)
430 6 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
431 :
432 : ! Compute XC matrix in AO basis
433 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, &
434 6 : vxc_rho=v_rspace, vxc_tau=tau_rspace, exc=exc)
435 :
436 6 : IF (ASSOCIATED(v_rspace)) THEN
437 6 : CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace, matrix_vxc=matrix_vxc)
438 :
439 16 : DO i = 1, SIZE(v_rspace)
440 16 : CALL v_rspace(i)%release()
441 : END DO
442 6 : DEALLOCATE (v_rspace)
443 :
444 6 : CALL cp_fm_get_info(mo_coeff(1), nrow_global=nao, ncol_global=nmo)
445 :
446 16 : DO i = 1, SIZE(matrix_vxc)
447 10 : CALL cp_fm_set_all(fm_XC_ao, 0.0_dp)
448 10 : CALL copy_dbcsr_to_fm(matrix=matrix_vxc(i)%matrix, fm=fm_XC_ao)
449 10 : CALL cp_fm_set_all(fm_XC_mo(i), 0.0_dp)
450 :
451 : ! First index
452 : CALL parallel_gemm('T', 'N', nmo, nao, nao, 1.0_dp, &
453 10 : mo_coeff(i), fm_XC_ao, 0.0_dp, fm_XC_ao_mo)
454 :
455 : ! Second index
456 : CALL parallel_gemm('N', 'N', nmo, nmo, nao, 1.0_dp, &
457 16 : fm_XC_ao_mo, mo_coeff(i), 1.0_dp, fm_XC_mo(i))
458 :
459 : END DO
460 :
461 16 : DO i = 1, SIZE(matrix_vxc)
462 10 : CALL dbcsr_release(matrix_vxc(i)%matrix)
463 16 : DEALLOCATE (matrix_vxc(i)%matrix)
464 : END DO
465 6 : DEALLOCATE (matrix_vxc)
466 : END IF
467 :
468 6 : CALL timestop(handle)
469 :
470 6 : END SUBROUTINE xc_contribution
471 :
472 : ! **************************************************************************************************
473 : !> \brief ...
474 : !> \param fm_F_mo ...
475 : !> \param eigenval ...
476 : !> \param homo ...
477 : !> \param para_env ...
478 : !> \param blacs_env ...
479 : !> \param rse_corr ...
480 : ! **************************************************************************************************
481 6 : SUBROUTINE non_diag_rse(fm_F_mo, eigenval, homo, para_env, &
482 : blacs_env, rse_corr)
483 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_F_mo
484 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Eigenval
485 : INTEGER, DIMENSION(:), INTENT(IN) :: homo
486 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
487 : TYPE(cp_blacs_env_type), INTENT(IN), POINTER :: blacs_env
488 : REAL(KIND=dp), INTENT(OUT) :: rse_corr
489 :
490 : CHARACTER(LEN=*), PARAMETER :: routineN = 'non_diag_rse'
491 :
492 : INTEGER :: handle, i_global, iiB, ispin, j_global, &
493 : jjB, ncol_local, nmo, nrow_local, &
494 : nspins, virtual
495 6 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
496 : REAL(KIND=dp) :: corr
497 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig_o, eig_semi_can, eig_v
498 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
499 : TYPE(cp_fm_type) :: fm_F_oo, fm_F_ov, fm_F_vv, fm_O, fm_tmp, &
500 : fm_U
501 :
502 6 : CALL timeset(routineN, handle)
503 :
504 6 : nmo = SIZE(Eigenval, 1)
505 6 : nspins = SIZE(fm_f_mo)
506 :
507 16 : DO ispin = 1, nspins
508 : ! Add eigenvalues on the diagonal
509 : CALL cp_fm_get_info(matrix=fm_F_mo(ispin), &
510 : nrow_local=nrow_local, &
511 : ncol_local=ncol_local, &
512 : row_indices=row_indices, &
513 10 : col_indices=col_indices)
514 :
515 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
516 16 : !$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_F_mo,eigenval,ispin)
517 : DO jjB = 1, ncol_local
518 : j_global = col_indices(jjB)
519 : DO iiB = 1, nrow_local
520 : i_global = row_indices(iiB)
521 : IF (i_global .EQ. j_global) fm_F_mo(ispin)%local_data(iib, jjb) = &
522 : fm_F_mo(ispin)%local_data(iib, jjb) + eigenval(i_global, ispin)
523 : END DO
524 : END DO
525 : !$OMP END PARALLEL DO
526 : END DO
527 :
528 6 : rse_corr = 0.0_dp
529 :
530 16 : DO ispin = 1, nspins
531 10 : IF (homo(ispin) <= 0 .OR. homo(ispin) >= nmo) CYCLE
532 : ! Create the occupied-occupied and virtual-virtual blocks, eigenvectors
533 10 : NULLIFY (fm_struct_tmp)
534 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
535 10 : nrow_global=homo(ispin), ncol_global=homo(ispin))
536 10 : CALL cp_fm_create(fm_F_oo, fm_struct_tmp, name="F_oo")
537 10 : CALL cp_fm_create(fm_O, fm_struct_tmp, name="O")
538 10 : CALL cp_fm_set_all(fm_F_oo, 0.0_dp)
539 10 : CALL cp_fm_set_all(fm_O, 0.0_dp)
540 10 : CALL cp_fm_struct_release(fm_struct_tmp)
541 :
542 : CALL cp_fm_to_fm_submat(msource=fm_F_mo(ispin), mtarget=fm_F_oo, &
543 : nrow=homo(ispin), ncol=homo(ispin), &
544 : s_firstrow=1, s_firstcol=1, &
545 10 : t_firstrow=1, t_firstcol=1)
546 10 : virtual = nmo - homo(ispin)
547 10 : NULLIFY (fm_struct_tmp)
548 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
549 10 : nrow_global=virtual, ncol_global=virtual)
550 10 : CALL cp_fm_create(fm_F_vv, fm_struct_tmp, name="F_vv")
551 10 : CALL cp_fm_create(fm_U, fm_struct_tmp, name="U")
552 10 : CALL cp_fm_set_all(fm_F_vv, 0.0_dp)
553 10 : CALL cp_fm_set_all(fm_U, 0.0_dp)
554 10 : CALL cp_fm_struct_release(fm_struct_tmp)
555 :
556 : CALL cp_fm_to_fm_submat(msource=fm_F_mo(ispin), mtarget=fm_F_vv, &
557 : nrow=virtual, ncol=virtual, &
558 : s_firstrow=homo(ispin) + 1, s_firstcol=homo(ispin) + 1, &
559 10 : t_firstrow=1, t_firstcol=1)
560 :
561 : ! Diagonalize occupied-occupied and virtual-virtual matrices
562 30 : ALLOCATE (eig_o(homo(ispin)))
563 30 : ALLOCATE (eig_v(virtual))
564 88 : eig_v = 0.0_dp
565 24 : eig_o = 0.0_dp
566 10 : CALL choose_eigv_solver(fm_F_oo, fm_O, eig_o)
567 10 : CALL choose_eigv_solver(fm_F_vv, fm_U, eig_v)
568 :
569 : ! Collect the eigenvalues to one array
570 30 : ALLOCATE (eig_semi_can(nmo))
571 102 : eig_semi_can = 0.0_dp
572 24 : eig_semi_can(1:homo(ispin)) = eig_o(:)
573 88 : eig_semi_can(homo(ispin) + 1:nmo) = eig_v(:)
574 :
575 : ! Create occupied-virtual block
576 10 : NULLIFY (fm_struct_tmp)
577 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
578 10 : nrow_global=homo(ispin), ncol_global=virtual)
579 10 : CALL cp_fm_create(fm_F_ov, fm_struct_tmp, name="F_ov")
580 10 : CALL cp_fm_create(fm_tmp, fm_struct_tmp, name="tmp")
581 10 : CALL cp_fm_set_all(fm_F_ov, 0.0_dp)
582 10 : CALL cp_fm_set_all(fm_tmp, 0.0_dp)
583 10 : CALL cp_fm_struct_release(fm_struct_tmp)
584 :
585 : CALL cp_fm_to_fm_submat(msource=fm_F_mo(ispin), mtarget=fm_F_ov, &
586 : nrow=homo(ispin), ncol=virtual, &
587 : s_firstrow=1, s_firstcol=homo(ispin) + 1, &
588 10 : t_firstrow=1, t_firstcol=1)
589 :
590 : CALL parallel_gemm(transa='T', transb='N', m=homo(ispin), n=virtual, k=homo(ispin), alpha=1.0_dp, &
591 10 : matrix_a=fm_O, matrix_b=fm_F_ov, beta=0.0_dp, matrix_c=fm_tmp)
592 :
593 : CALL parallel_gemm(transa='N', transb='N', m=homo(ispin), n=virtual, k=virtual, alpha=1.0_dp, &
594 10 : matrix_a=fm_tmp, matrix_b=fm_U, beta=0.0_dp, matrix_c=fm_F_ov)
595 :
596 : ! Compute the correction
597 : CALL cp_fm_get_info(matrix=fm_F_ov, &
598 : nrow_local=nrow_local, &
599 : ncol_local=ncol_local, &
600 : row_indices=row_indices, &
601 10 : col_indices=col_indices)
602 10 : corr = 0.0_dp
603 : !$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
604 : !$OMP REDUCTION(+:corr) &
605 10 : !$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_F_ov,eig_semi_can,homo,ispin)
606 : DO jjB = 1, ncol_local
607 : j_global = col_indices(jjB)
608 : DO iiB = 1, nrow_local
609 : i_global = row_indices(iiB)
610 : corr = corr + fm_F_ov%local_data(iib, jjb)**2.0_dp/ &
611 : (eig_semi_can(i_global) - eig_semi_can(j_global + homo(ispin)))
612 : END DO
613 : END DO
614 : !$OMP END PARALLEL DO
615 :
616 10 : rse_corr = rse_corr + corr
617 :
618 : ! Release
619 10 : DEALLOCATE (eig_semi_can)
620 10 : DEALLOCATE (eig_o)
621 10 : DEALLOCATE (eig_v)
622 :
623 10 : CALL cp_fm_release(fm_F_ov)
624 10 : CALL cp_fm_release(fm_F_oo)
625 10 : CALL cp_fm_release(fm_F_vv)
626 10 : CALL cp_fm_release(fm_U)
627 10 : CALL cp_fm_release(fm_O)
628 56 : CALL cp_fm_release(fm_tmp)
629 :
630 : END DO
631 :
632 6 : CALL para_env%sum(rse_corr)
633 :
634 6 : CALL timestop(handle)
635 :
636 12 : END SUBROUTINE non_diag_rse
637 :
638 : END MODULE rpa_rse
|