(git:858b7a1)
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-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! **************************************************************************************************
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 .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 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) :: 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
436 IF (ASSOCIATED(v_rspace)) THEN
437 CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace, matrix_vxc=matrix_vxc)
438
439 DO i = 1, SIZE(v_rspace)
440 CALL v_rspace(i)%release()
441 END DO
442 DEALLOCATE (v_rspace)
443
444 CALL cp_fm_get_info(mo_coeff(1), nrow_global=nao, ncol_global=nmo)
445
446 DO i = 1, SIZE(matrix_vxc)
447 CALL cp_fm_set_all(fm_xc_ao, 0.0_dp)
448 CALL copy_dbcsr_to_fm(matrix=matrix_vxc(i)%matrix, fm=fm_xc_ao)
449 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 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 fm_xc_ao_mo, mo_coeff(i), 1.0_dp, fm_xc_mo(i))
458
459 END DO
460
461 DO i = 1, SIZE(matrix_vxc)
462 CALL dbcsr_release(matrix_vxc(i)%matrix)
463 DEALLOCATE (matrix_vxc(i)%matrix)
464 END DO
465 DEALLOCATE (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 .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 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:192
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
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:1179
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:205
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1195
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1133
cp_dbcsr_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)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:214
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_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:167
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:238
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_diag
subroutine, public cp_fm_get_diag(matrix, diag)
returns the diagonal elements of a fm
Definition cp_fm_types.F:563
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:1009
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:1472
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:528
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:164
hfx_energy_potential
Routines to calculate HFX energy and potential.
Definition hfx_energy_potential.F:14
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)
computes four center integrals for a full basis set and updates the Kohn-Sham-Matrix and energy....
Definition hfx_energy_potential.F:137
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:735
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:760
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:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
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:23
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, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:518
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)
compute matrix_vxc, defined via the potential created by qs_vxc_create ignores things like tau functi...
Definition qs_ks_utils.F:1175
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)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Definition qs_vxc.F:98
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:586
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
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:763
mp2_types::mp2_type
Definition mp2_types.F:345
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:219
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:136
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