(git:3f0d10b)
Loading...
Searching...
No Matches
rpa_rse.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines to compute singles correction to RPA (RSE)
10!> \par History
11!> 08.2019 created [Vladimir Rybkin]
12!> \author Vladimir Rybkin
13! **************************************************************************************************
14MODULE 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 dbcsr_deallocate_matrix_set
31 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
32 USE cp_fm_diag, ONLY: choose_eigv_solver
33 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
34 cp_fm_struct_release,&
35 cp_fm_struct_type
36 USE cp_fm_types, ONLY: cp_fm_create,&
37 cp_fm_get_diag,&
38 cp_fm_get_info,&
39 cp_fm_release,&
40 cp_fm_set_all,&
41 cp_fm_to_fm_submat,&
42 cp_fm_type
43 USE hfx_energy_potential, ONLY: integrate_four_center
44 USE hfx_exx, ONLY: exx_post_hfx,&
45 exx_pre_hfx
46 USE hfx_ri, ONLY: hfx_ri_update_ks
47 USE input_section_types, ONLY: section_vals_get,&
48 section_vals_get_subs_vals,&
49 section_vals_type,&
50 section_vals_val_get
51 USE kinds, ONLY: dp
52 USE message_passing, ONLY: mp_para_env_type
53 USE mp2_types, ONLY: mp2_type
54 USE parallel_gemm_api, ONLY: parallel_gemm
55 USE pw_types, ONLY: pw_r3d_rs_type
56 USE qs_energy_types, ONLY: qs_energy_type
57 USE qs_environment_types, ONLY: get_qs_env,&
58 qs_environment_type
59 USE qs_ks_types, ONLY: qs_ks_env_type
60 USE qs_ks_utils, ONLY: compute_matrix_vxc
61 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
62 USE qs_rho_types, ONLY: qs_rho_get,&
63 qs_rho_type
64 USE qs_vxc, ONLY: qs_vxc_create
65
66!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
67
68#include "./base/base_uses.f90"
69
70 IMPLICIT NONE
71
72 PRIVATE
73
74 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_rse'
75
76 PUBLIC :: rse_energy
77
78CONTAINS
79
80! **************************************************************************************************
81!> \brief Single excitations energy corrections for RPA
82!> \param qs_env ...
83!> \param mp2_env ...
84!> \param para_env ...
85!> \param dft_control ...
86!> \param mo_coeff ...
87!> \param homo ...
88!> \param Eigenval ...
89!> \author Vladimir Rybkin, 08/2019
90! **************************************************************************************************
91 SUBROUTINE rse_energy(qs_env, mp2_env, para_env, dft_control, &
92 mo_coeff, homo, Eigenval)
93 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
94 TYPE(mp2_type), INTENT(INOUT) :: mp2_env
95 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
96 TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
97 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
98 INTEGER, DIMENSION(:), INTENT(IN) :: homo
99 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: eigenval
100
101 CHARACTER(LEN=*), PARAMETER :: routinen = 'rse_energy'
102
103 INTEGER :: handle, i_global, iib, ispin, j_global, &
104 jjb, n_rep_hf, nao, ncol_local, nmo, &
105 nrow_local, nspins
106 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
107 LOGICAL :: do_hfx, hfx_treat_lsd_in_core
108 REAL(kind=dp) :: coeff, corr, rse_corr
109 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: diag_diff
110 TYPE(cp_blacs_env_type), POINTER :: blacs_env
111 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
112 TYPE(cp_fm_type) :: fm_ao, fm_ao_mo
113 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_p_mu_nu, fm_x_mo, fm_xc_mo
114 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_mu_nu, matrix_s, rho_ao
115 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
116 POINTER :: sab_orb
117 TYPE(qs_energy_type), POINTER :: energy
118 TYPE(qs_rho_type), POINTER :: rho
119 TYPE(section_vals_type), POINTER :: hfx_sections, input
120
121 CALL timeset(routinen, handle)
122
123 nspins = dft_control%nspins
124
125 ! Pick the diagonal terms
126 CALL cp_fm_get_info(matrix=mo_coeff(1), &
127 nrow_local=nrow_local, &
128 ncol_local=ncol_local, &
129 row_indices=row_indices, &
130 col_indices=col_indices, &
131 nrow_global=nao, &
132 ncol_global=nmo)
133
134 ! start collecting stuff
135 NULLIFY (input, matrix_s, blacs_env, rho, energy, sab_orb)
136 CALL get_qs_env(qs_env, &
137 input=input, &
138 matrix_s=matrix_s, &
139 blacs_env=blacs_env, &
140 rho=rho, &
141 energy=energy, &
142 sab_orb=sab_orb)
143
144 CALL qs_rho_get(rho, rho_ao=rho_ao)
145
146 ! hfx section
147 NULLIFY (hfx_sections)
148 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
149 CALL section_vals_get(hfx_sections, explicit=do_hfx, n_repetition=n_rep_hf)
150 IF (do_hfx) THEN
151 CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
152 i_rep_section=1)
153 END IF
154
155 ! create work array
156 NULLIFY (mat_mu_nu)
157 CALL dbcsr_allocate_matrix_set(mat_mu_nu, nspins)
158 DO ispin = 1, nspins
159 ALLOCATE (mat_mu_nu(ispin)%matrix)
160 CALL dbcsr_create(matrix=mat_mu_nu(ispin)%matrix, template=matrix_s(1)%matrix, name="T_mu_nu", &
161 matrix_type=dbcsr_type_symmetric)
162 CALL cp_dbcsr_alloc_block_from_nbl(mat_mu_nu(ispin)%matrix, sab_orb)
163 CALL dbcsr_set(mat_mu_nu(ispin)%matrix, 0.0_dp)
164 END DO
165
166 ! Dense (full) matrices
167 ALLOCATE (fm_p_mu_nu(nspins))
168 NULLIFY (fm_struct_tmp)
169 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
170 nrow_global=nao, ncol_global=nao)
171 DO ispin = 1, nspins
172 CALL cp_fm_create(fm_p_mu_nu(ispin), fm_struct_tmp, name="P_mu_nu")
173 CALL cp_fm_set_all(fm_p_mu_nu(ispin), 0.0_dp)
174 END DO
175 CALL cp_fm_create(fm_ao, fm_struct_tmp, name="f_ao")
176 CALL cp_fm_struct_release(fm_struct_tmp)
177 CALL cp_fm_set_all(fm_ao, 0.0_dp)
178 CALL cp_fm_struct_release(fm_struct_tmp)
179
180 NULLIFY (fm_struct_tmp)
181 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
182 nrow_global=nmo, ncol_global=nmo)
183 ALLOCATE (fm_x_mo(nspins), fm_xc_mo(nspins))
184 DO ispin = 1, nspins
185 CALL cp_fm_create(fm_x_mo(ispin), fm_struct_tmp, name="f_X_mo")
186 CALL cp_fm_create(fm_xc_mo(ispin), fm_struct_tmp, name="f_XC_mo")
187 CALL cp_fm_set_all(fm_x_mo(ispin), 0.0_dp)
188 CALL cp_fm_set_all(fm_xc_mo(ispin), 0.0_dp)
189 END DO
190 CALL cp_fm_struct_release(fm_struct_tmp)
191
192 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
193 nrow_global=nmo, ncol_global=nao)
194 CALL cp_fm_create(fm_ao_mo, fm_struct_tmp, name="f_ao_mo")
195 CALL cp_fm_struct_release(fm_struct_tmp)
196 CALL cp_fm_set_all(fm_ao_mo, 0.0_dp)
197
198 !
199 ! Ready with preparations, do the real staff
200 !
201
202 ! Obtain density matrix like quantity
203
204 coeff = 1.0_dp
205 IF (nspins == 1) coeff = 2.0_dp
206 DO ispin = 1, nspins
207 CALL parallel_gemm(transa='N', transb='T', m=nao, n=nao, k=homo(ispin), alpha=coeff, &
208 matrix_a=mo_coeff(ispin), matrix_b=mo_coeff(ispin), &
209 beta=0.0_dp, matrix_c=fm_p_mu_nu(ispin))
210 END DO
211
212 ! Calculate exact exchange contribution
213 CALL exchange_contribution(qs_env, para_env, mo_coeff, &
214 hfx_sections, n_rep_hf, &
215 rho, mat_mu_nu, fm_p_mu_nu, &
216 fm_ao, fm_x_mo, fm_ao_mo)
217
218 ! Calculate DFT exchange-correlation contribution
219 CALL xc_contribution(qs_env, fm_ao, fm_ao_mo, fm_xc_mo, mo_coeff)
220
221 ALLOCATE (diag_diff(nmo))
222 rse_corr = 0.0_dp
223
224 DO ispin = 1, nspins
225 ! Compute the correction matrix: it is stored in fm_X_mo
226 CALL cp_fm_scale_and_add(1.0_dp, fm_x_mo(ispin), -1.0_dp, fm_xc_mo(ispin))
227
228 ! Pick the diagonal terms
229 CALL cp_fm_get_diag(fm_x_mo(ispin), diag_diff)
230
231 ! Compute the correction
232 CALL cp_fm_get_info(matrix=fm_x_mo(ispin), &
233 nrow_local=nrow_local, &
234 ncol_local=ncol_local, &
235 row_indices=row_indices, &
236 col_indices=col_indices)
237
238 corr = 0.0_dp
239
240!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
241!$OMP REDUCTION(+: corr) &
242!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,diag_diff,eigenval,fm_X_mo,homo,ispin)
243 DO jjb = 1, ncol_local
244 j_global = col_indices(jjb)
245 DO iib = 1, nrow_local
246 i_global = row_indices(iib)
247 IF ((i_global <= homo(ispin)) .AND. (j_global > homo(ispin))) THEN
248 corr = corr + fm_x_mo(ispin)%local_data(iib, jjb)**2.0_dp/ &
249 (eigenval(i_global, ispin) - eigenval(j_global, ispin) - diag_diff(i_global) + diag_diff(j_global))
250 END IF
251 END DO
252 END DO
253!$OMP END PARALLEL DO
254
255 rse_corr = rse_corr + corr
256 END DO
257
258 CALL para_env%sum(rse_corr)
259
260 IF (nspins == 1) rse_corr = rse_corr*2.0_dp
261
262 mp2_env%ri_rpa%rse_corr_diag = rse_corr
263
264 CALL non_diag_rse(fm_x_mo, eigenval, homo, para_env, blacs_env, rse_corr)
265
266 IF (nspins == 1) rse_corr = rse_corr*2.0_dp
267
268 mp2_env%ri_rpa%rse_corr = rse_corr
269
270 ! Release staff
271 DEALLOCATE (diag_diff)
272 CALL cp_fm_release(fm_ao)
273 CALL cp_fm_release(fm_ao_mo)
274 CALL cp_fm_release(fm_p_mu_nu)
275 CALL cp_fm_release(fm_x_mo)
276 CALL cp_fm_release(fm_xc_mo)
277 DO ispin = 1, nspins
278 CALL dbcsr_release(mat_mu_nu(ispin)%matrix)
279 DEALLOCATE (mat_mu_nu(ispin)%matrix)
280 END DO
281 DEALLOCATE (mat_mu_nu)
282
283 CALL timestop(handle)
284
285 END SUBROUTINE rse_energy
286
287! **************************************************************************************************
288!> \brief HF exchange occupied-virtual matrix
289!> \param qs_env ...
290!> \param para_env ...
291!> \param mo_coeff ...
292!> \param hfx_sections ...
293!> \param n_rep_hf ...
294!> \param rho_work ...
295!> \param mat_mu_nu ...
296!> \param fm_P_mu_nu ...
297!> \param fm_X_ao ...
298!> \param fm_X_mo ...
299!> \param fm_X_ao_mo ...
300! **************************************************************************************************
301 SUBROUTINE exchange_contribution(qs_env, para_env, mo_coeff, &
302 hfx_sections, n_rep_hf, &
303 rho_work, mat_mu_nu, fm_P_mu_nu, &
304 fm_X_ao, fm_X_mo, fm_X_ao_mo)
305 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
306 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
307 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
308 TYPE(section_vals_type), INTENT(IN), POINTER :: hfx_sections
309 INTEGER, INTENT(IN) :: n_rep_hf
310 TYPE(qs_rho_type), INTENT(IN), POINTER :: rho_work
311 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
312 POINTER :: mat_mu_nu
313 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_p_mu_nu
314 TYPE(cp_fm_type), INTENT(INOUT) :: fm_x_ao
315 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_x_mo
316 TYPE(cp_fm_type), INTENT(IN) :: fm_x_ao_mo
317
318 CHARACTER(LEN=*), PARAMETER :: routinen = 'exchange_contribution'
319
320 INTEGER :: handle, irep, is, nao, nmo, ns
321 LOGICAL :: my_recalc_hfx_integrals
322 REAL(kind=dp) :: ehfx
323 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_mu_nu, rho_work_ao
324 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_2d, rho_ao_2d
325
326 CALL timeset(routinen, handle)
327
328 CALL cp_fm_get_info(mo_coeff(1), nrow_global=nao, ncol_global=nmo)
329
330 CALL qs_rho_get(rho_work, rho_ao=rho_work_ao)
331 ns = SIZE(rho_work_ao)
332 NULLIFY (p_mu_nu)
333 CALL dbcsr_allocate_matrix_set(p_mu_nu, ns)
334 DO is = 1, ns
335 CALL dbcsr_init_p(p_mu_nu(is)%matrix)
336 CALL dbcsr_create(p_mu_nu(is)%matrix, template=rho_work_ao(1)%matrix)
337 CALL dbcsr_copy(p_mu_nu(is)%matrix, rho_work_ao(1)%matrix)
338 CALL dbcsr_set(p_mu_nu(is)%matrix, 0.0_dp)
339 END DO
340
341 my_recalc_hfx_integrals = .true.
342
343 CALL exx_pre_hfx(hfx_sections, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
344 DO is = 1, ns
345 CALL copy_fm_to_dbcsr(fm_p_mu_nu(is), p_mu_nu(1)%matrix, keep_sparsity=.true.)
346
347 CALL dbcsr_set(mat_mu_nu(1)%matrix, 0.0_dp)
348
349 IF (qs_env%mp2_env%ri_rpa%x_data(1, 1)%do_hfx_ri) THEN
350
351 DO irep = 1, n_rep_hf
352 rho_ao_2d(1:ns, 1:1) => p_mu_nu(1:ns)
353 mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns)
354 CALL hfx_ri_update_ks(qs_env, qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data, mat_2d, ehfx, &
355 rho_ao=rho_ao_2d, geometry_did_change=my_recalc_hfx_integrals, nspins=1, &
356 hf_fraction=qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction)
357
358 IF (ns == 2) CALL dbcsr_scale(mat_mu_nu(1)%matrix, 2.0_dp)
359 my_recalc_hfx_integrals = .false.
360 END DO
361
362 ELSE
363
364 DO irep = 1, n_rep_hf
365 rho_ao_2d(1:ns, 1:1) => p_mu_nu(1:ns)
366 mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns)
367 CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
368 para_env, my_recalc_hfx_integrals, irep, .true., &
369 ispin=1)
370
371 my_recalc_hfx_integrals = .false.
372 END DO
373 END IF
374
375 ! copy back to fm
376 CALL cp_fm_set_all(fm_x_ao, 0.0_dp)
377 CALL copy_dbcsr_to_fm(matrix=mat_mu_nu(1)%matrix, fm=fm_x_ao)
378 CALL cp_fm_set_all(fm_x_mo(is), 0.0_dp)
379
380 ! First index
381 CALL parallel_gemm('T', 'N', nmo, nao, nmo, 1.0_dp, &
382 mo_coeff(is), fm_x_ao, 0.0_dp, fm_x_ao_mo)
383
384 ! Second index
385 CALL parallel_gemm('N', 'N', nmo, nmo, nao, 1.0_dp, &
386 fm_x_ao_mo, mo_coeff(is), 1.0_dp, fm_x_mo(is))
387
388 END DO
389 CALL exx_post_hfx(qs_env, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
390
391 ! Release dbcsr objects
392 DO is = 1, SIZE(p_mu_nu)
393 CALL dbcsr_release(p_mu_nu(is)%matrix)
394 DEALLOCATE (p_mu_nu(is)%matrix)
395 END DO
396 DEALLOCATE (p_mu_nu)
397
398 CALL timestop(handle)
399
400 END SUBROUTINE exchange_contribution
401
402! **************************************************************************************************
403!> \brief Exchange-correlation occupied-virtual matrix
404!> \param qs_env ...
405!> \param fm_XC_ao ...
406!> \param fm_XC_ao_mo ...
407!> \param fm_XC_mo ...
408!> \param mo_coeff ...
409! **************************************************************************************************
410 SUBROUTINE xc_contribution(qs_env, fm_XC_ao, fm_XC_ao_mo, fm_XC_mo, mo_coeff)
411 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
412 TYPE(cp_fm_type), INTENT(INOUT) :: fm_xc_ao
413 TYPE(cp_fm_type), INTENT(IN) :: fm_xc_ao_mo
414 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_xc_mo, mo_coeff
415
416 CHARACTER(LEN=*), PARAMETER :: routinen = 'xc_contribution'
417
418 INTEGER :: handle, i, nao, nmo
419 REAL(kind=dp) :: edisp, exc
420 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_vxc
421 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau_rspace, v_rspace
422 TYPE(qs_ks_env_type), POINTER :: ks_env
423 TYPE(qs_rho_type), POINTER :: rho
424 TYPE(section_vals_type), POINTER :: input, xc_section
425
426 CALL timeset(routinen, handle)
427
428 NULLIFY (matrix_vxc, v_rspace, tau_rspace, input, xc_section, ks_env, &
429 rho)
430 CALL get_qs_env(qs_env, input=input, ks_env=ks_env, rho=rho)
431 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
432
433 ! Compute XC matrix in AO basis
434 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, &
435 vxc_rho=v_rspace, vxc_tau=tau_rspace, exc=exc, &
436 edisp=edisp, dispersion_env=qs_env%dispersion_env)
437
438 IF (ASSOCIATED(v_rspace)) THEN
439 ! This temporary belongs to the RSE calculation; the cached XC matrices stay untouched.
440 CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace, matrix_vxc=matrix_vxc)
441 cpassert(SIZE(matrix_vxc, 2) == 1)
442
443 DO i = 1, SIZE(v_rspace)
444 CALL v_rspace(i)%release()
445 END DO
446 DEALLOCATE (v_rspace)
447
448 CALL cp_fm_get_info(mo_coeff(1), nrow_global=nao, ncol_global=nmo)
449
450 DO i = 1, SIZE(matrix_vxc, 1)
451 CALL cp_fm_set_all(fm_xc_ao, 0.0_dp)
452 CALL copy_dbcsr_to_fm(matrix=matrix_vxc(i, 1)%matrix, fm=fm_xc_ao)
453 CALL cp_fm_set_all(fm_xc_mo(i), 0.0_dp)
454
455 ! First index
456 CALL parallel_gemm('T', 'N', nmo, nao, nao, 1.0_dp, &
457 mo_coeff(i), fm_xc_ao, 0.0_dp, fm_xc_ao_mo)
458
459 ! Second index
460 CALL parallel_gemm('N', 'N', nmo, nmo, nao, 1.0_dp, &
461 fm_xc_ao_mo, mo_coeff(i), 1.0_dp, fm_xc_mo(i))
462
463 END DO
464
465 CALL dbcsr_deallocate_matrix_set(matrix_vxc)
466 END IF
467
468 CALL timestop(handle)
469
470 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 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 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
496 REAL(kind=dp) :: corr
497 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 CALL timeset(routinen, handle)
503
504 nmo = SIZE(eigenval, 1)
505 nspins = SIZE(fm_f_mo)
506
507 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 col_indices=col_indices)
514
515!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
516!$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 == 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 rse_corr = 0.0_dp
529
530 DO ispin = 1, nspins
531 IF (homo(ispin) <= 0 .OR. homo(ispin) >= nmo) cycle
532 ! Create the occupied-occupied and virtual-virtual blocks, eigenvectors
533 NULLIFY (fm_struct_tmp)
534 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
535 nrow_global=homo(ispin), ncol_global=homo(ispin))
536 CALL cp_fm_create(fm_f_oo, fm_struct_tmp, name="F_oo")
537 CALL cp_fm_create(fm_o, fm_struct_tmp, name="O")
538 CALL cp_fm_set_all(fm_f_oo, 0.0_dp)
539 CALL cp_fm_set_all(fm_o, 0.0_dp)
540 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 t_firstrow=1, t_firstcol=1)
546 virtual = nmo - homo(ispin)
547 NULLIFY (fm_struct_tmp)
548 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
549 nrow_global=virtual, ncol_global=virtual)
550 CALL cp_fm_create(fm_f_vv, fm_struct_tmp, name="F_vv")
551 CALL cp_fm_create(fm_u, fm_struct_tmp, name="U")
552 CALL cp_fm_set_all(fm_f_vv, 0.0_dp)
553 CALL cp_fm_set_all(fm_u, 0.0_dp)
554 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 t_firstrow=1, t_firstcol=1)
560
561 ! Diagonalize occupied-occupied and virtual-virtual matrices
562 ALLOCATE (eig_o(homo(ispin)))
563 ALLOCATE (eig_v(virtual))
564 eig_v = 0.0_dp
565 eig_o = 0.0_dp
566 CALL choose_eigv_solver(fm_f_oo, fm_o, eig_o)
567 CALL choose_eigv_solver(fm_f_vv, fm_u, eig_v)
568
569 ! Collect the eigenvalues to one array
570 ALLOCATE (eig_semi_can(nmo))
571 eig_semi_can = 0.0_dp
572 eig_semi_can(1:homo(ispin)) = eig_o(:)
573 eig_semi_can(homo(ispin) + 1:nmo) = eig_v(:)
574
575 ! Create occupied-virtual block
576 NULLIFY (fm_struct_tmp)
577 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
578 nrow_global=homo(ispin), ncol_global=virtual)
579 CALL cp_fm_create(fm_f_ov, fm_struct_tmp, name="F_ov")
580 CALL cp_fm_create(fm_tmp, fm_struct_tmp, name="tmp")
581 CALL cp_fm_set_all(fm_f_ov, 0.0_dp)
582 CALL cp_fm_set_all(fm_tmp, 0.0_dp)
583 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 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 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 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 col_indices=col_indices)
602 corr = 0.0_dp
603!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
604!$OMP REDUCTION(+:corr) &
605!$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 rse_corr = rse_corr + corr
617
618 ! Release
619 DEALLOCATE (eig_semi_can)
620 DEALLOCATE (eig_o)
621 DEALLOCATE (eig_v)
622
623 CALL cp_fm_release(fm_f_ov)
624 CALL cp_fm_release(fm_f_oo)
625 CALL cp_fm_release(fm_f_vv)
626 CALL cp_fm_release(fm_u)
627 CALL cp_fm_release(fm_o)
628 CALL cp_fm_release(fm_tmp)
629
630 END DO
631
632 CALL para_env%sum(rse_corr)
633
634 CALL timestop(handle)
635
636 END SUBROUTINE non_diag_rse
637
638END MODULE rpa_rse
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:96
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:104
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
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_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1188
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan, check_plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:344
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:240
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:171
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:262
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1046
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag, n_active)
returns the diagonal elements of a fm
Definition cp_fm_types.F:567
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition cp_fm_types.F:1518
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:536
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
hfx_energy_potential
Routines to calculate HFX energy and potential.
Definition hfx_energy_potential.F:15
hfx_energy_potential::integrate_four_center
subroutine, public integrate_four_center(qs_env, x_data, ks_matrix, ehfx, rho_ao, hfx_section, para_env, geometry_did_change, irep, distribute_fock_matrix, ispin, nspins)
computes four center integrals for a full basis set and updates the Kohn-Sham-Matrix and energy....
Definition hfx_energy_potential.F:145
hfx_exx
Routines to calculate EXX in RPA and energy correction methods.
Definition hfx_exx.F:16
hfx_exx::exx_pre_hfx
subroutine, public exx_pre_hfx(ext_hfx_section, x_data, reuse_hfx)
Prepare the external x_data for integration. Simply change the HFX fraction in case the qs_envx_data ...
Definition hfx_exx.F:739
hfx_exx::exx_post_hfx
subroutine, public exx_post_hfx(qs_env, x_data, reuse_hfx)
Revert back to the proper HFX fraction in case qs_envx_data is reused.
Definition hfx_exx.F:764
hfx_ri
RI-methods for HFX.
Definition hfx_ri.F:12
hfx_ri::hfx_ri_update_ks
subroutine, public hfx_ri_update_ks(qs_env, ri_data, ks_matrix, ehfx, mos, rho_ao, geometry_did_change, nspins, hf_fraction)
...
Definition hfx_ri.F:1041
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
mp2_types
Types needed for MP2 calculations.
Definition mp2_types.F:14
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
pw_types
Definition pw_types.F:24
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_utils
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition qs_ks_utils.F:22
qs_ks_utils::compute_matrix_vxc
subroutine, public compute_matrix_vxc(qs_env, v_rspace, matrix_vxc, gapw_full_basis)
Integrate the scalar XC potential into spin/image-resolved AO matrices. Gamma-only calculations have ...
Definition qs_ks_utils.F:1433
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
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_vxc
Definition qs_vxc.F:16
qs_vxc::qs_vxc_create
subroutine, public qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, just_energy, edisp, dispersion_env, adiabatic_rescale_factor, pw_env_external, native_skala_atom_force, qs_env_external, native_gapw_composite_override, native_skala_defer_to_atom_composite)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Definition qs_vxc.F:120
rpa_rse
Routines to compute singles correction to RPA (RSE).
Definition rpa_rse.F:14
rpa_rse::rse_energy
subroutine, public rse_energy(qs_env, mp2_env, para_env, dft_control, mo_coeff, homo, eigenval)
Single excitations energy corrections for RPA.
Definition rpa_rse.F:93
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:55
cp_control_types::dft_control_type
Definition cp_control_types.F:786
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:774
mp2_types::mp2_type
Definition mp2_types.F:312
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
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