(git:6344619)
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 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
77CONTAINS
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 SUBROUTINE rse_energy(qs_env, mp2_env, para_env, dft_control, &
91 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 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 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 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_p_mu_nu, fm_x_mo, fm_xc_mo
113 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_mu_nu, matrix_s, rho_ao
114 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
115 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 CALL timeset(routinen, handle)
121
122 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 ncol_global=nmo)
132
133 ! start collecting stuff
134 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 sab_orb=sab_orb)
142
143 CALL qs_rho_get(rho, rho_ao=rho_ao)
144
145 ! hfx section
146 NULLIFY (hfx_sections)
147 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
148 CALL section_vals_get(hfx_sections, explicit=do_hfx, n_repetition=n_rep_hf)
149 IF (do_hfx) THEN
150 CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
151 i_rep_section=1)
152 END IF
153
154 ! create work array
155 NULLIFY (mat_mu_nu)
156 CALL dbcsr_allocate_matrix_set(mat_mu_nu, nspins)
157 DO ispin = 1, nspins
158 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 matrix_type=dbcsr_type_symmetric)
161 CALL cp_dbcsr_alloc_block_from_nbl(mat_mu_nu(ispin)%matrix, sab_orb)
162 CALL dbcsr_set(mat_mu_nu(ispin)%matrix, 0.0_dp)
163 END DO
164
165 ! Dense (full) matrices
166 ALLOCATE (fm_p_mu_nu(nspins))
167 NULLIFY (fm_struct_tmp)
168 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
169 nrow_global=nao, ncol_global=nao)
170 DO ispin = 1, nspins
171 CALL cp_fm_create(fm_p_mu_nu(ispin), fm_struct_tmp, name="P_mu_nu")
172 CALL cp_fm_set_all(fm_p_mu_nu(ispin), 0.0_dp)
173 END DO
174 CALL cp_fm_create(fm_ao, fm_struct_tmp, name="f_ao")
175 CALL cp_fm_struct_release(fm_struct_tmp)
176 CALL cp_fm_set_all(fm_ao, 0.0_dp)
177 CALL cp_fm_struct_release(fm_struct_tmp)
178
179 NULLIFY (fm_struct_tmp)
180 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
181 nrow_global=nmo, ncol_global=nmo)
182 ALLOCATE (fm_x_mo(nspins), fm_xc_mo(nspins))
183 DO ispin = 1, nspins
184 CALL cp_fm_create(fm_x_mo(ispin), fm_struct_tmp, name="f_X_mo")
185 CALL cp_fm_create(fm_xc_mo(ispin), fm_struct_tmp, name="f_XC_mo")
186 CALL cp_fm_set_all(fm_x_mo(ispin), 0.0_dp)
187 CALL cp_fm_set_all(fm_xc_mo(ispin), 0.0_dp)
188 END DO
189 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 nrow_global=nmo, ncol_global=nao)
193 CALL cp_fm_create(fm_ao_mo, fm_struct_tmp, name="f_ao_mo")
194 CALL cp_fm_struct_release(fm_struct_tmp)
195 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 coeff = 1.0_dp
204 IF (nspins == 1) coeff = 2.0_dp
205 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 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 fm_ao, fm_x_mo, fm_ao_mo)
216
217 ! Calculate DFT exchange-correlation contribution
218 CALL xc_contribution(qs_env, fm_ao, fm_ao_mo, fm_xc_mo, mo_coeff)
219
220 ALLOCATE (diag_diff(nmo))
221 rse_corr = 0.0_dp
222
223 DO ispin = 1, nspins
224 ! Compute the correction matrix: it is stored in fm_X_mo
225 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 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 col_indices=col_indices)
236
237 corr = 0.0_dp
238
239!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
240!$OMP REDUCTION(+: corr) &
241!$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 <= homo(ispin)) .AND. (j_global > 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 rse_corr = rse_corr + corr
255 END DO
256
257 CALL para_env%sum(rse_corr)
258
259 IF (nspins == 1) rse_corr = rse_corr*2.0_dp
260
261 mp2_env%ri_rpa%rse_corr_diag = rse_corr
262
263 CALL non_diag_rse(fm_x_mo, eigenval, homo, para_env, blacs_env, rse_corr)
264
265 IF (nspins == 1) rse_corr = rse_corr*2.0_dp
266
267 mp2_env%ri_rpa%rse_corr = rse_corr
268
269 ! Release staff
270 DEALLOCATE (diag_diff)
271 CALL cp_fm_release(fm_ao)
272 CALL cp_fm_release(fm_ao_mo)
273 CALL cp_fm_release(fm_p_mu_nu)
274 CALL cp_fm_release(fm_x_mo)
275 CALL cp_fm_release(fm_xc_mo)
276 DO ispin = 1, nspins
277 CALL dbcsr_release(mat_mu_nu(ispin)%matrix)
278 DEALLOCATE (mat_mu_nu(ispin)%matrix)
279 END DO
280 DEALLOCATE (mat_mu_nu)
281
282 CALL timestop(handle)
283
284 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 SUBROUTINE exchange_contribution(qs_env, para_env, mo_coeff, &
301 hfx_sections, n_rep_hf, &
302 rho_work, mat_mu_nu, fm_P_mu_nu, &
303 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_mu_nu, rho_work_ao
323 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_2d, rho_ao_2d
324
325 CALL timeset(routinen, handle)
326
327 CALL cp_fm_get_info(mo_coeff(1), nrow_global=nao, ncol_global=nmo)
328
329 CALL qs_rho_get(rho_work, rho_ao=rho_work_ao)
330 ns = SIZE(rho_work_ao)
331 NULLIFY (p_mu_nu)
332 CALL dbcsr_allocate_matrix_set(p_mu_nu, ns)
333 DO is = 1, ns
334 CALL dbcsr_init_p(p_mu_nu(is)%matrix)
335 CALL dbcsr_create(p_mu_nu(is)%matrix, template=rho_work_ao(1)%matrix)
336 CALL dbcsr_copy(p_mu_nu(is)%matrix, rho_work_ao(1)%matrix)
337 CALL dbcsr_set(p_mu_nu(is)%matrix, 0.0_dp)
338 END DO
339
340 my_recalc_hfx_integrals = .true.
341
342 CALL exx_pre_hfx(hfx_sections, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
343 DO is = 1, ns
344 CALL copy_fm_to_dbcsr(fm_p_mu_nu(is), p_mu_nu(1)%matrix, keep_sparsity=.true.)
345
346 CALL dbcsr_set(mat_mu_nu(1)%matrix, 0.0_dp)
347
348 IF (qs_env%mp2_env%ri_rpa%x_data(1, 1)%do_hfx_ri) THEN
349
350 DO irep = 1, n_rep_hf
351 rho_ao_2d(1:ns, 1:1) => p_mu_nu(1:ns)
352 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 hf_fraction=qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction)
356
357 IF (ns == 2) CALL dbcsr_scale(mat_mu_nu(1)%matrix, 2.0_dp)
358 my_recalc_hfx_integrals = .false.
359 END DO
360
361 ELSE
362
363 DO irep = 1, n_rep_hf
364 rho_ao_2d(1:ns, 1:1) => p_mu_nu(1:ns)
365 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 ispin=1)
369
370 my_recalc_hfx_integrals = .false.
371 END DO
372 END IF
373
374 ! copy back to fm
375 CALL cp_fm_set_all(fm_x_ao, 0.0_dp)
376 CALL copy_dbcsr_to_fm(matrix=mat_mu_nu(1)%matrix, fm=fm_x_ao)
377 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 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 fm_x_ao_mo, mo_coeff(is), 1.0_dp, fm_x_mo(is))
386
387 END DO
388 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 DO is = 1, SIZE(p_mu_nu)
392 CALL dbcsr_release(p_mu_nu(is)%matrix)
393 DEALLOCATE (p_mu_nu(is)%matrix)
394 END DO
395 DEALLOCATE (p_mu_nu)
396
397 CALL timestop(handle)
398
399 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 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) :: edisp, exc
419 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc
420 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 CALL timeset(routinen, handle)
426
427 NULLIFY (matrix_vxc, v_rspace, tau_rspace, input, xc_section, ks_env, &
428 rho)
429 CALL get_qs_env(qs_env, matrix_vxc=matrix_vxc, input=input, ks_env=ks_env, rho=rho)
430 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 vxc_rho=v_rspace, vxc_tau=tau_rspace, exc=exc, &
435 edisp=edisp, dispersion_env=qs_env%dispersion_env)
436
437 IF (ASSOCIATED(v_rspace)) THEN
438 CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace, matrix_vxc=matrix_vxc)
439
440 DO i = 1, SIZE(v_rspace)
441 CALL v_rspace(i)%release()
442 END DO
443 DEALLOCATE (v_rspace)
444
445 CALL cp_fm_get_info(mo_coeff(1), nrow_global=nao, ncol_global=nmo)
446
447 DO i = 1, SIZE(matrix_vxc)
448 CALL cp_fm_set_all(fm_xc_ao, 0.0_dp)
449 CALL copy_dbcsr_to_fm(matrix=matrix_vxc(i)%matrix, fm=fm_xc_ao)
450 CALL cp_fm_set_all(fm_xc_mo(i), 0.0_dp)
451
452 ! First index
453 CALL parallel_gemm('T', 'N', nmo, nao, nao, 1.0_dp, &
454 mo_coeff(i), fm_xc_ao, 0.0_dp, fm_xc_ao_mo)
455
456 ! Second index
457 CALL parallel_gemm('N', 'N', nmo, nmo, nao, 1.0_dp, &
458 fm_xc_ao_mo, mo_coeff(i), 1.0_dp, fm_xc_mo(i))
459
460 END DO
461
462 DO i = 1, SIZE(matrix_vxc)
463 CALL dbcsr_release(matrix_vxc(i)%matrix)
464 DEALLOCATE (matrix_vxc(i)%matrix)
465 END DO
466 DEALLOCATE (matrix_vxc)
467 END IF
468
469 CALL timestop(handle)
470
471 END SUBROUTINE xc_contribution
472
473! **************************************************************************************************
474!> \brief ...
475!> \param fm_F_mo ...
476!> \param eigenval ...
477!> \param homo ...
478!> \param para_env ...
479!> \param blacs_env ...
480!> \param rse_corr ...
481! **************************************************************************************************
482 SUBROUTINE non_diag_rse(fm_F_mo, eigenval, homo, para_env, &
483 blacs_env, rse_corr)
484 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_f_mo
485 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: eigenval
486 INTEGER, DIMENSION(:), INTENT(IN) :: homo
487 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
488 TYPE(cp_blacs_env_type), INTENT(IN), POINTER :: blacs_env
489 REAL(kind=dp), INTENT(OUT) :: rse_corr
490
491 CHARACTER(LEN=*), PARAMETER :: routinen = 'non_diag_rse'
492
493 INTEGER :: handle, i_global, iib, ispin, j_global, &
494 jjb, ncol_local, nmo, nrow_local, &
495 nspins, virtual
496 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
497 REAL(kind=dp) :: corr
498 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eig_o, eig_semi_can, eig_v
499 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
500 TYPE(cp_fm_type) :: fm_f_oo, fm_f_ov, fm_f_vv, fm_o, fm_tmp, &
501 fm_u
502
503 CALL timeset(routinen, handle)
504
505 nmo = SIZE(eigenval, 1)
506 nspins = SIZE(fm_f_mo)
507
508 DO ispin = 1, nspins
509 ! Add eigenvalues on the diagonal
510 CALL cp_fm_get_info(matrix=fm_f_mo(ispin), &
511 nrow_local=nrow_local, &
512 ncol_local=ncol_local, &
513 row_indices=row_indices, &
514 col_indices=col_indices)
515
516!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
517!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_F_mo,eigenval,ispin)
518 DO jjb = 1, ncol_local
519 j_global = col_indices(jjb)
520 DO iib = 1, nrow_local
521 i_global = row_indices(iib)
522 IF (i_global == j_global) fm_f_mo(ispin)%local_data(iib, jjb) = &
523 fm_f_mo(ispin)%local_data(iib, jjb) + eigenval(i_global, ispin)
524 END DO
525 END DO
526!$OMP END PARALLEL DO
527 END DO
528
529 rse_corr = 0.0_dp
530
531 DO ispin = 1, nspins
532 IF (homo(ispin) <= 0 .OR. homo(ispin) >= nmo) cycle
533 ! Create the occupied-occupied and virtual-virtual blocks, eigenvectors
534 NULLIFY (fm_struct_tmp)
535 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
536 nrow_global=homo(ispin), ncol_global=homo(ispin))
537 CALL cp_fm_create(fm_f_oo, fm_struct_tmp, name="F_oo")
538 CALL cp_fm_create(fm_o, fm_struct_tmp, name="O")
539 CALL cp_fm_set_all(fm_f_oo, 0.0_dp)
540 CALL cp_fm_set_all(fm_o, 0.0_dp)
541 CALL cp_fm_struct_release(fm_struct_tmp)
542
543 CALL cp_fm_to_fm_submat(msource=fm_f_mo(ispin), mtarget=fm_f_oo, &
544 nrow=homo(ispin), ncol=homo(ispin), &
545 s_firstrow=1, s_firstcol=1, &
546 t_firstrow=1, t_firstcol=1)
547 virtual = nmo - homo(ispin)
548 NULLIFY (fm_struct_tmp)
549 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
550 nrow_global=virtual, ncol_global=virtual)
551 CALL cp_fm_create(fm_f_vv, fm_struct_tmp, name="F_vv")
552 CALL cp_fm_create(fm_u, fm_struct_tmp, name="U")
553 CALL cp_fm_set_all(fm_f_vv, 0.0_dp)
554 CALL cp_fm_set_all(fm_u, 0.0_dp)
555 CALL cp_fm_struct_release(fm_struct_tmp)
556
557 CALL cp_fm_to_fm_submat(msource=fm_f_mo(ispin), mtarget=fm_f_vv, &
558 nrow=virtual, ncol=virtual, &
559 s_firstrow=homo(ispin) + 1, s_firstcol=homo(ispin) + 1, &
560 t_firstrow=1, t_firstcol=1)
561
562 ! Diagonalize occupied-occupied and virtual-virtual matrices
563 ALLOCATE (eig_o(homo(ispin)))
564 ALLOCATE (eig_v(virtual))
565 eig_v = 0.0_dp
566 eig_o = 0.0_dp
567 CALL choose_eigv_solver(fm_f_oo, fm_o, eig_o)
568 CALL choose_eigv_solver(fm_f_vv, fm_u, eig_v)
569
570 ! Collect the eigenvalues to one array
571 ALLOCATE (eig_semi_can(nmo))
572 eig_semi_can = 0.0_dp
573 eig_semi_can(1:homo(ispin)) = eig_o(:)
574 eig_semi_can(homo(ispin) + 1:nmo) = eig_v(:)
575
576 ! Create occupied-virtual block
577 NULLIFY (fm_struct_tmp)
578 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
579 nrow_global=homo(ispin), ncol_global=virtual)
580 CALL cp_fm_create(fm_f_ov, fm_struct_tmp, name="F_ov")
581 CALL cp_fm_create(fm_tmp, fm_struct_tmp, name="tmp")
582 CALL cp_fm_set_all(fm_f_ov, 0.0_dp)
583 CALL cp_fm_set_all(fm_tmp, 0.0_dp)
584 CALL cp_fm_struct_release(fm_struct_tmp)
585
586 CALL cp_fm_to_fm_submat(msource=fm_f_mo(ispin), mtarget=fm_f_ov, &
587 nrow=homo(ispin), ncol=virtual, &
588 s_firstrow=1, s_firstcol=homo(ispin) + 1, &
589 t_firstrow=1, t_firstcol=1)
590
591 CALL parallel_gemm(transa='T', transb='N', m=homo(ispin), n=virtual, k=homo(ispin), alpha=1.0_dp, &
592 matrix_a=fm_o, matrix_b=fm_f_ov, beta=0.0_dp, matrix_c=fm_tmp)
593
594 CALL parallel_gemm(transa='N', transb='N', m=homo(ispin), n=virtual, k=virtual, alpha=1.0_dp, &
595 matrix_a=fm_tmp, matrix_b=fm_u, beta=0.0_dp, matrix_c=fm_f_ov)
596
597 ! Compute the correction
598 CALL cp_fm_get_info(matrix=fm_f_ov, &
599 nrow_local=nrow_local, &
600 ncol_local=ncol_local, &
601 row_indices=row_indices, &
602 col_indices=col_indices)
603 corr = 0.0_dp
604!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
605!$OMP REDUCTION(+:corr) &
606!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_F_ov,eig_semi_can,homo,ispin)
607 DO jjb = 1, ncol_local
608 j_global = col_indices(jjb)
609 DO iib = 1, nrow_local
610 i_global = row_indices(iib)
611 corr = corr + fm_f_ov%local_data(iib, jjb)**2.0_dp/ &
612 (eig_semi_can(i_global) - eig_semi_can(j_global + homo(ispin)))
613 END DO
614 END DO
615!$OMP END PARALLEL DO
616
617 rse_corr = rse_corr + corr
618
619 ! Release
620 DEALLOCATE (eig_semi_can)
621 DEALLOCATE (eig_o)
622 DEALLOCATE (eig_v)
623
624 CALL cp_fm_release(fm_f_ov)
625 CALL cp_fm_release(fm_f_oo)
626 CALL cp_fm_release(fm_f_vv)
627 CALL cp_fm_release(fm_u)
628 CALL cp_fm_release(fm_o)
629 CALL cp_fm_release(fm_tmp)
630
631 END DO
632
633 CALL para_env%sum(rse_corr)
634
635 CALL timestop(handle)
636
637 END SUBROUTINE non_diag_rse
638
639END 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:95
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
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:233
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:129
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:1045
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:566
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:1500
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:535
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:166
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:150
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)
compute matrix_vxc, defined via the potential created by qs_vxc_create ignores things like tau functi...
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:92
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: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:115
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:768
mp2_types::mp2_type
Definition mp2_types.F:309
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