(git:c9ed8e0)
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 nmo ...
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, nmo, 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, INTENT(IN) :: nmo
99 INTEGER, DIMENSION(:), INTENT(IN) :: homo
100 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: eigenval
101
102 CHARACTER(LEN=*), PARAMETER :: routinen = 'rse_energy'
103
104 INTEGER :: dimen, handle, i_global, iib, ispin, &
105 j_global, jjb, n_rep_hf, ncol_local, &
106 nrow_local, nspins
107 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
108 LOGICAL :: do_hfx, hfx_treat_lsd_in_core
109 REAL(kind=dp) :: coeff, corr, rse_corr
110 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: diag_diff
111 TYPE(cp_blacs_env_type), POINTER :: blacs_env
112 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
113 TYPE(cp_fm_type) :: fm_ao, fm_ao_mo
114 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_p_mu_nu, fm_x_mo, fm_xc_mo
115 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_mu_nu, matrix_s, rho_ao
116 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
117 POINTER :: sab_orb
118 TYPE(qs_energy_type), POINTER :: energy
119 TYPE(qs_rho_type), POINTER :: rho
120 TYPE(section_vals_type), POINTER :: hfx_sections, input
121
122 CALL timeset(routinen, handle)
123
124 nspins = dft_control%nspins
125
126 ! Pick the diagonal terms
127 CALL cp_fm_get_info(matrix=mo_coeff(1), &
128 nrow_local=nrow_local, &
129 ncol_local=ncol_local, &
130 row_indices=row_indices, &
131 col_indices=col_indices)
132
133 ! start collecting stuff
134 dimen = nmo
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=dimen, ncol_global=dimen)
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_struct_release(fm_struct_tmp)
176
177 NULLIFY (fm_struct_tmp)
178 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
179 nrow_global=dimen, ncol_global=dimen)
180 ALLOCATE (fm_x_mo(nspins), fm_xc_mo(nspins))
181 DO ispin = 1, nspins
182 CALL cp_fm_create(fm_x_mo(ispin), fm_struct_tmp, name="f_X_mo")
183 CALL cp_fm_create(fm_xc_mo(ispin), fm_struct_tmp, name="f_XC_mo")
184 CALL cp_fm_set_all(fm_x_mo(ispin), 0.0_dp)
185 CALL cp_fm_set_all(fm_xc_mo(ispin), 0.0_dp)
186 END DO
187 CALL cp_fm_struct_release(fm_struct_tmp)
188
189 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
190 nrow_global=dimen, ncol_global=dimen)
191 CALL cp_fm_create(fm_ao, fm_struct_tmp, name="f_ao")
192 CALL cp_fm_struct_release(fm_struct_tmp)
193 CALL cp_fm_set_all(fm_ao, 0.0_dp)
194 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
195 nrow_global=dimen, ncol_global=dimen)
196 CALL cp_fm_create(fm_ao_mo, fm_struct_tmp, name="f_ao_mo")
197 CALL cp_fm_struct_release(fm_struct_tmp)
198 CALL cp_fm_set_all(fm_ao_mo, 0.0_dp)
199
200 !
201 ! Ready with preparations, do the real staff
202 !
203
204 ! Obtain density matrix like quantity
205
206 coeff = 1.0_dp
207 IF (nspins == 1) coeff = 2.0_dp
208 DO ispin = 1, nspins
209 CALL parallel_gemm(transa='N', transb='T', m=dimen, n=dimen, k=homo(ispin), alpha=coeff, &
210 matrix_a=mo_coeff(ispin), matrix_b=mo_coeff(ispin), &
211 beta=0.0_dp, matrix_c=fm_p_mu_nu(ispin))
212 END DO
213
214 ! Calculate exact exchange contribution
215 CALL exchange_contribution(qs_env, para_env, dimen, mo_coeff, &
216 hfx_sections, n_rep_hf, &
217 rho, mat_mu_nu, fm_p_mu_nu, &
218 fm_ao, fm_x_mo, fm_ao_mo)
219
220 ! Calculate DFT exchange-correlation contribution
221 CALL xc_contribution(qs_env, fm_ao, fm_ao_mo, fm_xc_mo, mo_coeff, dimen)
222
223 ALLOCATE (diag_diff(dimen))
224 rse_corr = 0.0_dp
225
226 DO ispin = 1, nspins
227 ! Compute the correction matrix: it is stored in fm_X_mo
228 CALL cp_fm_scale_and_add(1.0_dp, fm_x_mo(ispin), -1.0_dp, fm_xc_mo(ispin))
229
230 ! Pick the diagonal terms
231 CALL cp_fm_get_diag(fm_x_mo(ispin), diag_diff)
232
233 ! Compute the correction
234 CALL cp_fm_get_info(matrix=fm_x_mo(ispin), &
235 nrow_local=nrow_local, &
236 ncol_local=ncol_local, &
237 row_indices=row_indices, &
238 col_indices=col_indices)
239
240 corr = 0.0_dp
241
242!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
243!$OMP REDUCTION(+: corr) &
244!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,diag_diff,eigenval,fm_X_mo,homo,ispin)
245 DO jjb = 1, ncol_local
246 j_global = col_indices(jjb)
247 DO iib = 1, nrow_local
248 i_global = row_indices(iib)
249 IF ((i_global .LE. homo(ispin)) .AND. (j_global .GT. homo(ispin))) THEN
250 corr = corr + fm_x_mo(ispin)%local_data(iib, jjb)**2.0_dp/ &
251 (eigenval(i_global, ispin) - eigenval(j_global, ispin) - diag_diff(i_global) + diag_diff(j_global))
252 END IF
253 END DO
254 END DO
255!$OMP END PARALLEL DO
256
257 rse_corr = rse_corr + corr
258 END DO
259
260 CALL para_env%sum(rse_corr)
261
262 IF (nspins == 1) rse_corr = rse_corr*2.0_dp
263
264 mp2_env%ri_rpa%rse_corr_diag = rse_corr
265
266 CALL non_diag_rse(fm_x_mo, eigenval, dimen, homo, para_env, blacs_env, rse_corr)
267
268 IF (nspins == 1) rse_corr = rse_corr*2.0_dp
269
270 mp2_env%ri_rpa%rse_corr = rse_corr
271
272 ! Release staff
273 DEALLOCATE (diag_diff)
274 CALL cp_fm_release(fm_ao)
275 CALL cp_fm_release(fm_ao_mo)
276 CALL cp_fm_release(fm_p_mu_nu)
277 CALL cp_fm_release(fm_x_mo)
278 CALL cp_fm_release(fm_xc_mo)
279 DO ispin = 1, nspins
280 CALL dbcsr_release(mat_mu_nu(ispin)%matrix)
281 DEALLOCATE (mat_mu_nu(ispin)%matrix)
282 END DO
283 DEALLOCATE (mat_mu_nu)
284
285 CALL timestop(handle)
286
287 END SUBROUTINE rse_energy
288
289! **************************************************************************************************
290!> \brief HF exchange occupied-virtual matrix
291!> \param qs_env ...
292!> \param para_env ...
293!> \param dimen ...
294!> \param mo_coeff ...
295!> \param hfx_sections ...
296!> \param n_rep_hf ...
297!> \param rho_work ...
298!> \param mat_mu_nu ...
299!> \param fm_P_mu_nu ...
300!> \param fm_X_ao ...
301!> \param fm_X_mo ...
302!> \param fm_X_ao_mo ...
303! **************************************************************************************************
304 SUBROUTINE exchange_contribution(qs_env, para_env, dimen, mo_coeff, &
305 hfx_sections, n_rep_hf, &
306 rho_work, mat_mu_nu, fm_P_mu_nu, &
307 fm_X_ao, fm_X_mo, fm_X_ao_mo)
308 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
309 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
310 INTEGER, INTENT(IN) :: dimen
311 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
312 TYPE(section_vals_type), INTENT(IN), POINTER :: hfx_sections
313 INTEGER, INTENT(IN) :: n_rep_hf
314 TYPE(qs_rho_type), INTENT(IN), POINTER :: rho_work
315 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
316 POINTER :: mat_mu_nu
317 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_p_mu_nu
318 TYPE(cp_fm_type), INTENT(INOUT) :: fm_x_ao
319 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_x_mo
320 TYPE(cp_fm_type), INTENT(IN) :: fm_x_ao_mo
321
322 CHARACTER(LEN=*), PARAMETER :: routinen = 'exchange_contribution'
323
324 INTEGER :: handle, irep, is, ns
325 LOGICAL :: my_recalc_hfx_integrals
326 REAL(kind=dp) :: ehfx
327 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_mu_nu, rho_work_ao
328 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_2d, rho_ao_2d
329
330 CALL timeset(routinen, handle)
331
332 CALL qs_rho_get(rho_work, rho_ao=rho_work_ao)
333 ns = SIZE(rho_work_ao)
334 NULLIFY (p_mu_nu)
335 CALL dbcsr_allocate_matrix_set(p_mu_nu, ns)
336 DO is = 1, ns
337 CALL dbcsr_init_p(p_mu_nu(is)%matrix)
338 CALL dbcsr_create(p_mu_nu(is)%matrix, template=rho_work_ao(1)%matrix)
339 CALL dbcsr_copy(p_mu_nu(is)%matrix, rho_work_ao(1)%matrix)
340 CALL dbcsr_set(p_mu_nu(is)%matrix, 0.0_dp)
341 END DO
342
343 my_recalc_hfx_integrals = .true.
344
345 CALL exx_pre_hfx(hfx_sections, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
346 DO is = 1, ns
347 CALL copy_fm_to_dbcsr(fm_p_mu_nu(is), p_mu_nu(1)%matrix, keep_sparsity=.true.)
348
349 CALL dbcsr_set(mat_mu_nu(1)%matrix, 0.0_dp)
350
351 IF (qs_env%mp2_env%ri_rpa%x_data(1, 1)%do_hfx_ri) THEN
352
353 DO irep = 1, n_rep_hf
354 rho_ao_2d(1:ns, 1:1) => p_mu_nu(1:ns)
355 mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns)
356 CALL hfx_ri_update_ks(qs_env, qs_env%mp2_env%ri_rpa%x_data(irep, 1)%ri_data, mat_2d, ehfx, &
357 rho_ao=rho_ao_2d, geometry_did_change=my_recalc_hfx_integrals, nspins=1, &
358 hf_fraction=qs_env%mp2_env%ri_rpa%x_data(irep, 1)%general_parameter%fraction)
359
360 IF (ns == 2) CALL dbcsr_scale(mat_mu_nu(1)%matrix, 2.0_dp)
361 my_recalc_hfx_integrals = .false.
362 END DO
363
364 ELSE
365
366 DO irep = 1, n_rep_hf
367 rho_ao_2d(1:ns, 1:1) => p_mu_nu(1:ns)
368 mat_2d(1:ns, 1:1) => mat_mu_nu(1:ns)
369 CALL integrate_four_center(qs_env, qs_env%mp2_env%ri_rpa%x_data, mat_2d, ehfx, rho_ao_2d, hfx_sections, &
370 para_env, my_recalc_hfx_integrals, irep, .true., &
371 ispin=1)
372
373 my_recalc_hfx_integrals = .false.
374 END DO
375 END IF
376
377 ! copy back to fm
378 CALL cp_fm_set_all(fm_x_ao, 0.0_dp)
379 CALL copy_dbcsr_to_fm(matrix=mat_mu_nu(1)%matrix, fm=fm_x_ao)
380 CALL cp_fm_set_all(fm_x_mo(is), 0.0_dp)
381
382 ! First index
383 CALL parallel_gemm('T', 'N', dimen, dimen, dimen, 1.0_dp, &
384 mo_coeff(is), fm_x_ao, 0.0_dp, fm_x_ao_mo)
385
386 ! Second index
387 CALL parallel_gemm('N', 'N', dimen, dimen, dimen, 1.0_dp, &
388 fm_x_ao_mo, mo_coeff(is), 1.0_dp, fm_x_mo(is))
389
390 END DO
391 CALL exx_post_hfx(qs_env, qs_env%mp2_env%ri_rpa%x_data, qs_env%mp2_env%ri_rpa%reuse_hfx)
392
393 ! Release dbcsr objects
394 DO is = 1, SIZE(p_mu_nu)
395 CALL dbcsr_release(p_mu_nu(is)%matrix)
396 DEALLOCATE (p_mu_nu(is)%matrix)
397 END DO
398 DEALLOCATE (p_mu_nu)
399
400 CALL timestop(handle)
401
402 END SUBROUTINE exchange_contribution
403
404! **************************************************************************************************
405!> \brief Exchange-correlation occupied-virtual matrix
406!> \param qs_env ...
407!> \param fm_XC_ao ...
408!> \param fm_XC_ao_mo ...
409!> \param fm_XC_mo ...
410!> \param mo_coeff ...
411!> \param dimen ...
412! **************************************************************************************************
413 SUBROUTINE xc_contribution(qs_env, fm_XC_ao, fm_XC_ao_mo, fm_XC_mo, mo_coeff, dimen)
414 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
415 TYPE(cp_fm_type), INTENT(INOUT) :: fm_xc_ao
416 TYPE(cp_fm_type), INTENT(IN) :: fm_xc_ao_mo
417 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_xc_mo, mo_coeff
418 INTEGER, INTENT(IN) :: dimen
419
420 CHARACTER(LEN=*), PARAMETER :: routinen = 'xc_contribution'
421
422 INTEGER :: handle, i
423 REAL(kind=dp) :: exc
424 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc
425 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: tau_rspace, v_rspace
426 TYPE(qs_ks_env_type), POINTER :: ks_env
427 TYPE(qs_rho_type), POINTER :: rho
428 TYPE(section_vals_type), POINTER :: input, xc_section
429
430 CALL timeset(routinen, handle)
431
432 NULLIFY (matrix_vxc, v_rspace, tau_rspace, input, xc_section, ks_env, &
433 rho)
434 CALL get_qs_env(qs_env, matrix_vxc=matrix_vxc, input=input, ks_env=ks_env, rho=rho)
435 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
436
437 ! Compute XC matrix in AO basis
438 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, &
439 vxc_rho=v_rspace, vxc_tau=tau_rspace, exc=exc)
440
441 IF (ASSOCIATED(v_rspace)) THEN
442 CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace, matrix_vxc=matrix_vxc)
443
444 DO i = 1, SIZE(v_rspace)
445 CALL v_rspace(i)%release()
446 END DO
447 DEALLOCATE (v_rspace)
448
449 DO i = 1, SIZE(matrix_vxc)
450 CALL cp_fm_set_all(fm_xc_ao, 0.0_dp)
451 CALL copy_dbcsr_to_fm(matrix=matrix_vxc(i)%matrix, fm=fm_xc_ao)
452 CALL cp_fm_set_all(fm_xc_mo(i), 0.0_dp)
453
454 ! First index
455 CALL parallel_gemm('T', 'N', dimen, dimen, dimen, 1.0_dp, &
456 mo_coeff(i), fm_xc_ao, 0.0_dp, fm_xc_ao_mo)
457
458 ! Second index
459 CALL parallel_gemm('N', 'N', dimen, dimen, dimen, 1.0_dp, &
460 fm_xc_ao_mo, mo_coeff(i), 1.0_dp, fm_xc_mo(i))
461
462 END DO
463
464 DO i = 1, SIZE(matrix_vxc)
465 CALL dbcsr_release(matrix_vxc(i)%matrix)
466 DEALLOCATE (matrix_vxc(i)%matrix)
467 END DO
468 DEALLOCATE (matrix_vxc)
469 END IF
470
471 CALL timestop(handle)
472
473 END SUBROUTINE xc_contribution
474
475! **************************************************************************************************
476!> \brief ...
477!> \param fm_F_mo ...
478!> \param eigenval ...
479!> \param dimen ...
480!> \param homo ...
481!> \param para_env ...
482!> \param blacs_env ...
483!> \param rse_corr ...
484! **************************************************************************************************
485 SUBROUTINE non_diag_rse(fm_F_mo, eigenval, dimen, homo, para_env, &
486 blacs_env, rse_corr)
487 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_f_mo
488 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: eigenval
489 INTEGER, INTENT(IN) :: dimen
490 INTEGER, DIMENSION(:), INTENT(IN) :: homo
491 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
492 TYPE(cp_blacs_env_type), INTENT(IN), POINTER :: blacs_env
493 REAL(kind=dp), INTENT(OUT) :: rse_corr
494
495 CHARACTER(LEN=*), PARAMETER :: routinen = 'non_diag_rse'
496
497 INTEGER :: handle, i_global, iib, ispin, j_global, &
498 jjb, ncol_local, nrow_local, nspins, &
499 virtual
500 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
501 REAL(kind=dp) :: corr
502 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eig_o, eig_semi_can, eig_v
503 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
504 TYPE(cp_fm_type) :: fm_f_oo, fm_f_ov, fm_f_vv, fm_o, fm_tmp, &
505 fm_u
506
507 CALL timeset(routinen, handle)
508
509 nspins = SIZE(fm_f_mo)
510
511 DO ispin = 1, nspins
512 ! Add eigenvalues on the diagonal
513 CALL cp_fm_get_info(matrix=fm_f_mo(ispin), &
514 nrow_local=nrow_local, &
515 ncol_local=ncol_local, &
516 row_indices=row_indices, &
517 col_indices=col_indices)
518
519!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
520!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_F_mo,eigenval,ispin)
521 DO jjb = 1, ncol_local
522 j_global = col_indices(jjb)
523 DO iib = 1, nrow_local
524 i_global = row_indices(iib)
525 IF (i_global .EQ. j_global) fm_f_mo(ispin)%local_data(iib, jjb) = &
526 fm_f_mo(ispin)%local_data(iib, jjb) + eigenval(i_global, ispin)
527 END DO
528 END DO
529!$OMP END PARALLEL DO
530 END DO
531
532 rse_corr = 0.0_dp
533
534 DO ispin = 1, nspins
535 IF (homo(ispin) <= 0 .OR. homo(ispin) >= dimen) cycle
536 ! Create the occupied-occupied and virtual-virtual blocks, eigenvectors
537 NULLIFY (fm_struct_tmp)
538 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
539 nrow_global=homo(ispin), ncol_global=homo(ispin))
540 CALL cp_fm_create(fm_f_oo, fm_struct_tmp, name="F_oo")
541 CALL cp_fm_create(fm_o, fm_struct_tmp, name="O")
542 CALL cp_fm_set_all(fm_f_oo, 0.0_dp)
543 CALL cp_fm_set_all(fm_o, 0.0_dp)
544 CALL cp_fm_struct_release(fm_struct_tmp)
545
546 CALL cp_fm_to_fm_submat(msource=fm_f_mo(ispin), mtarget=fm_f_oo, &
547 nrow=homo(ispin), ncol=homo(ispin), &
548 s_firstrow=1, s_firstcol=1, &
549 t_firstrow=1, t_firstcol=1)
550 virtual = dimen - homo(ispin)
551 NULLIFY (fm_struct_tmp)
552 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
553 nrow_global=virtual, ncol_global=virtual)
554 CALL cp_fm_create(fm_f_vv, fm_struct_tmp, name="F_vv")
555 CALL cp_fm_create(fm_u, fm_struct_tmp, name="U")
556 CALL cp_fm_set_all(fm_f_vv, 0.0_dp)
557 CALL cp_fm_set_all(fm_u, 0.0_dp)
558 CALL cp_fm_struct_release(fm_struct_tmp)
559
560 CALL cp_fm_to_fm_submat(msource=fm_f_mo(ispin), mtarget=fm_f_vv, &
561 nrow=virtual, ncol=virtual, &
562 s_firstrow=homo(ispin) + 1, s_firstcol=homo(ispin) + 1, &
563 t_firstrow=1, t_firstcol=1)
564
565 ! Diagonalize occupied-occupied and virtual-virtual matrices
566 ALLOCATE (eig_o(homo(ispin)))
567 ALLOCATE (eig_v(virtual))
568 eig_v = 0.0_dp
569 eig_o = 0.0_dp
570 CALL choose_eigv_solver(fm_f_oo, fm_o, eig_o)
571 CALL choose_eigv_solver(fm_f_vv, fm_u, eig_v)
572
573 ! Collect the eigenvalues to one array
574 ALLOCATE (eig_semi_can(dimen))
575 eig_semi_can = 0.0_dp
576 eig_semi_can(1:homo(ispin)) = eig_o(:)
577 eig_semi_can(homo(ispin) + 1:dimen) = eig_v(:)
578
579 ! Create occupied-virtual block
580 NULLIFY (fm_struct_tmp)
581 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
582 nrow_global=homo(ispin), ncol_global=virtual)
583 CALL cp_fm_create(fm_f_ov, fm_struct_tmp, name="F_ov")
584 CALL cp_fm_create(fm_tmp, fm_struct_tmp, name="tmp")
585 CALL cp_fm_set_all(fm_f_ov, 0.0_dp)
586 CALL cp_fm_set_all(fm_tmp, 0.0_dp)
587 CALL cp_fm_struct_release(fm_struct_tmp)
588
589 CALL cp_fm_to_fm_submat(msource=fm_f_mo(ispin), mtarget=fm_f_ov, &
590 nrow=homo(ispin), ncol=virtual, &
591 s_firstrow=1, s_firstcol=homo(ispin) + 1, &
592 t_firstrow=1, t_firstcol=1)
593
594 CALL parallel_gemm(transa='T', transb='N', m=homo(ispin), n=virtual, k=homo(ispin), alpha=1.0_dp, &
595 matrix_a=fm_o, matrix_b=fm_f_ov, beta=0.0_dp, matrix_c=fm_tmp)
596
597 CALL parallel_gemm(transa='N', transb='N', m=homo(ispin), n=virtual, k=virtual, alpha=1.0_dp, &
598 matrix_a=fm_tmp, matrix_b=fm_u, beta=0.0_dp, matrix_c=fm_f_ov)
599
600 ! Compute the correction
601 CALL cp_fm_get_info(matrix=fm_f_ov, &
602 nrow_local=nrow_local, &
603 ncol_local=ncol_local, &
604 row_indices=row_indices, &
605 col_indices=col_indices)
606 corr = 0.0_dp
607!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(jjB,iiB,i_global,j_global) &
608!$OMP REDUCTION(+:corr) &
609!$OMP SHARED(ncol_local,nrow_local,col_indices,row_indices,fm_F_ov,eig_semi_can,homo,ispin)
610 DO jjb = 1, ncol_local
611 j_global = col_indices(jjb)
612 DO iib = 1, nrow_local
613 i_global = row_indices(iib)
614 corr = corr + fm_f_ov%local_data(iib, jjb)**2.0_dp/ &
615 (eig_semi_can(i_global) - eig_semi_can(j_global + homo(ispin)))
616 END DO
617 END DO
618!$OMP END PARALLEL DO
619
620 rse_corr = rse_corr + corr
621
622 ! Release
623 DEALLOCATE (eig_semi_can)
624 DEALLOCATE (eig_o)
625 DEALLOCATE (eig_v)
626
627 CALL cp_fm_release(fm_f_ov)
628 CALL cp_fm_release(fm_f_oo)
629 CALL cp_fm_release(fm_f_vv)
630 CALL cp_fm_release(fm_u)
631 CALL cp_fm_release(fm_o)
632 CALL cp_fm_release(fm_tmp)
633
634 END DO
635
636 CALL para_env%sum(rse_corr)
637
638 CALL timestop(handle)
639
640 END SUBROUTINE non_diag_rse
641
642END 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, nmo, 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: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