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