Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Static Coulomb-hole screened-exchange (COHSEX) quasiparticle energies and MOs in the
10 : !> O(N^4) (RI_RPA%GW) GW code for molecules
11 : !>
12 : !> Symbols
13 : !> μ,ν atomic orbitals (AOs)
14 : !> P,Q RI auxiliary functions (after the RI-basis reduction, dimen_RI)
15 : !> i / a occupied / virtual Kohn-Sham MOs, p any MO, m,n MOs in the window of
16 : !> corrected MOs [homo-CORR_MOS_OCC+1, homo+CORR_MOS_VIRT]
17 : !> ϵ_n DFT eigenvalue, ϵ_F Fermi level (between HOMO and LUMO)
18 : !> B^P_nm = Σ_Q (P|Q)^(-1/2) (Q|nm) RI 3-center integrals (fm_mat_S_gw)
19 : !>
20 : !> Static screened Coulomb interaction (RPA, ω = 0), closed shell:
21 : !> Q_PQ(0) = 4 Σ_ia B^P_ia B^Q_ia / (ϵ_a - ϵ_i)
22 : !> W^c_PQ(0) = [1 + Q(0)]^(-1)_PQ - δ_PQ (-δ_PQ removes the bare Coulomb part)
23 : !> W^c_mp,pn = Σ_PQ B^P_mp W^c_PQ(0) B^Q_pn
24 : !>
25 : !> static COHSEX self-energy = Coulomb hole (COH) + screened exchange (SEX):
26 : !> Σ^COH_mn = 1/2 Σ_p W^c_mp,pn
27 : !> Σ^SEX_mn = - Σ_i^occ W^c_mi,in (+ bare exchange Σ^x_mn)
28 : !> Σ^c_mn = Σ^COH_mn + Σ^SEX_mn = 1/2 Σ_p sgn(ϵ_p - ϵ_F) W^c_mp,pn
29 : !> Σ^x_mn = -Σ_i^occ Σ_P B^P_mi B^P_in (RI_SIGMA_X; else from HFX)
30 : !>
31 : !> Quasiparticle Hamiltonian in the window:
32 : !> H^QP_mn = ϵ_n δ_mn + Σ^x_mn - v^xc_mn + α_mn Σ^c_mn
33 : !> α_mn = α_IP = COHSEX_SCALE_OCC (m,n occupied), α_EA = COHSEX_SCALE_VIRT (m,n virtual), H^QP_ia = 0
34 : !> ϵ_n^DG = H^QP_nn (COHSEX DIAGONAL)
35 : !> H^QP U = U ϵ^CN, occ-occ and virt-virt block separately (COHSEX QP_UPDATE)
36 : !> C^QP_μn = Σ_m C_μm U_mn (COHSEX QP_UPDATE MOs)
37 : !> COHSEX DIAGONAL prints ϵ^DG only; COHSEX QP_UPDATE also diagonalizes H^QP, prints ϵ^CN
38 : !> and the canonical MOs.
39 : !>
40 : !> \par History
41 : !> 09.2026 created [Ritaj Tyagi]
42 : ! **************************************************************************************************
43 : MODULE rpa_gw_static_cohsex
44 : USE atomic_kind_types, ONLY: atomic_kind_type
45 : USE bibliography, ONLY: Hedin1965,&
46 : Hybertsen1986,&
47 : Tyagi2024,&
48 : Tyagi2026,&
49 : cite_reference
50 : USE cell_types, ONLY: cell_type
51 : USE cp_control_types, ONLY: dft_control_type
52 : USE cp_dbcsr_api, ONLY: dbcsr_release_p,&
53 : dbcsr_type
54 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
55 : USE cp_fm_basic_linalg, ONLY: cp_fm_add_on_diag,&
56 : cp_fm_column_scale
57 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
58 : cp_fm_cholesky_invert
59 : USE cp_fm_redistribution, ONLY: cp_fm_permute
60 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
61 : cp_fm_struct_release,&
62 : cp_fm_struct_type
63 : USE cp_fm_types, ONLY: cp_fm_create,&
64 : cp_fm_get_info,&
65 : cp_fm_get_submatrix,&
66 : cp_fm_release,&
67 : cp_fm_set_submatrix,&
68 : cp_fm_to_fm,&
69 : cp_fm_type
70 : USE cp_log_handling, ONLY: cp_get_default_logger,&
71 : cp_logger_type
72 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
73 : cp_print_key_should_output,&
74 : cp_print_key_unit_nr
75 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
76 : USE input_constants, ONLY: cohsex_qp_update
77 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
78 : section_vals_type,&
79 : section_vals_val_get
80 : USE kinds, ONLY: default_path_length,&
81 : dp
82 : USE mathlib, ONLY: diamat_all
83 : USE message_passing, ONLY: mp_para_env_type
84 : USE mp2_types, ONLY: mp2_type
85 : USE parallel_gemm_api, ONLY: parallel_gemm
86 : USE particle_list_types, ONLY: particle_list_type
87 : USE particle_types, ONLY: particle_type
88 : USE physcon, ONLY: evolt
89 : USE pw_env_types, ONLY: pw_env_get,&
90 : pw_env_type
91 : USE pw_pool_types, ONLY: pw_pool_type
92 : USE pw_types, ONLY: pw_c1d_gs_type,&
93 : pw_r3d_rs_type
94 : USE qs_collocate_density, ONLY: calculate_wavefunction
95 : USE qs_environment_types, ONLY: get_qs_env,&
96 : qs_environment_type
97 : USE qs_kind_types, ONLY: qs_kind_type
98 : USE qs_mo_io, ONLY: write_mo_set_to_output_unit
99 : USE qs_mo_types, ONLY: allocate_mo_set,&
100 : deallocate_mo_set,&
101 : init_mo_set,&
102 : mo_set_type
103 : USE qs_subsys_types, ONLY: qs_subsys_get,&
104 : qs_subsys_type
105 : #include "./base/base_uses.f90"
106 :
107 : IMPLICIT NONE
108 :
109 : PRIVATE
110 :
111 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rpa_gw_static_cohsex'
112 :
113 : PUBLIC :: static_cohsex_qp_energies, rotate_bse_slab, get_static_cohsex_mo_coeff
114 :
115 : CONTAINS
116 :
117 : ! **************************************************************************************************
118 : !> \brief Compute and print the diagonal (ϵ^DG) static COHSEX quasiparticle energies; with
119 : !> COHSEX QP_UPDATE also the canonical (ϵ^CN) energies and the canonical MOs
120 : !> \param fm_mat_Q on entry the spin-summed 1 + Q_PQ(ω=0); on exit W^c_PQ(0) = [1 + Q(0)]^(-1)_PQ - δ_PQ
121 : !> \param fm_mat_S_gw per-spin B^P_mp
122 : !> \param fm_mat_S_gw_work ...
123 : !> \param dimen_RI ...
124 : !> \param nmo ...
125 : !> \param homo ...
126 : !> \param gw_corr_lev_occ ...
127 : !> \param gw_corr_lev_virt ...
128 : !> \param Eigenval ...
129 : !> \param mo_coeff ...
130 : !> \param unit_nr ...
131 : !> \param qs_env ...
132 : !> \param mp2_env ...
133 : !> \param eigenval_dg ϵ_n^DG of all MOs (scissor outside the window), nmo x nspins
134 : !> \param eigenval_cn ϵ_n^CN of all MOs (scissor outside the window), nmo x nspins, QP_UPDATE only
135 : !> \param u_qp_all U_mn of all MOs (δ_mn outside the window), nmo x nmo x nspins, QP_UPDATE only
136 : ! **************************************************************************************************
137 16 : SUBROUTINE static_cohsex_qp_energies(fm_mat_Q, fm_mat_S_gw, fm_mat_S_gw_work, dimen_RI, nmo, homo, &
138 16 : gw_corr_lev_occ, gw_corr_lev_virt, Eigenval, mo_coeff, unit_nr, &
139 8 : qs_env, mp2_env, eigenval_dg, eigenval_cn, u_qp_all)
140 :
141 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_mat_Q
142 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mat_S_gw, fm_mat_S_gw_work
143 : INTEGER, INTENT(IN) :: dimen_RI, nmo
144 : INTEGER, DIMENSION(:), INTENT(IN) :: homo, gw_corr_lev_occ, gw_corr_lev_virt
145 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Eigenval
146 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
147 : INTEGER, INTENT(IN) :: unit_nr
148 : TYPE(qs_environment_type), POINTER :: qs_env
149 : TYPE(mp2_type) :: mp2_env
150 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
151 : OPTIONAL :: eigenval_dg, eigenval_cn
152 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT), &
153 : OPTIONAL :: u_qp_all
154 :
155 : CHARACTER(LEN=*), PARAMETER :: routineN = 'static_cohsex_qp_energies'
156 :
157 : CHARACTER(LEN=10) :: spin_label
158 : CHARACTER(LEN=3) :: occ_vir
159 : CHARACTER(LEN=40) :: cohsex_label
160 : INTEGER :: handle, ispin, m, n, n_mo_win, no, &
161 : nspins, nv, p, win_start
162 : LOGICAL :: qp_update
163 : REAL(KIND=dp) :: shift_occ, shift_virt
164 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: e_can, e_diag, e_qp, scale_vec
165 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: h_qp, sigma_c, sigma_x, &
166 8 : sigma_x_minus_vxc, u_qp, vxc
167 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_B_re, fm_struct_win
168 : TYPE(cp_fm_type) :: fm_B_re, fm_sigma, fm_WB_re, fm_work
169 : TYPE(mp_para_env_type), POINTER :: para_env
170 :
171 8 : CALL timeset(routineN, handle)
172 :
173 8 : nspins = SIZE(fm_mat_S_gw)
174 : ! COHSEX QP_UPDATE: diagonalize H^QP (canonical energies and MOs); DIAGONAL: ϵ^DG only
175 8 : qp_update = (mp2_env%ri_g0w0%cohsex_flavor == cohsex_qp_update)
176 :
177 : ! static COHSEX; with QP_UPDATE the quasiparticle update; scaled self-energy if requested
178 8 : CALL cite_reference(Hedin1965)
179 8 : CALL cite_reference(Hybertsen1986)
180 8 : IF (qp_update) CALL cite_reference(Tyagi2024)
181 8 : IF (mp2_env%ri_g0w0%cohsex_scale_occ /= 1.0_dp .OR. mp2_env%ri_g0w0%cohsex_scale_virt /= 1.0_dp) THEN
182 6 : CALL cite_reference(Tyagi2026)
183 : END IF
184 8 : CALL cp_fm_get_info(mo_coeff(1), para_env=para_env)
185 :
186 : ! ==============================================================================
187 : ! 1. Static correlation part of the screened Coulomb interaction
188 : ! ε_PQ(0) = δ_PQ + Q_PQ(0)
189 : ! W^c_PQ(0) = ε^(-1)_PQ(0) - δ_PQ
190 : ! ==============================================================================
191 8 : CALL cp_fm_cholesky_decompose(fm_mat_Q, n=dimen_RI)
192 8 : CALL cp_fm_create(fm_work, fm_mat_Q%matrix_struct)
193 8 : CALL cp_fm_cholesky_invert(fm_mat_Q, n=dimen_RI, work=fm_work)
194 8 : CALL cp_fm_release(fm_work)
195 8 : CALL cp_fm_add_on_diag(fm_mat_Q, -1.0_dp, n_active=dimen_RI)
196 :
197 16 : DO ispin = 1, nspins
198 :
199 8 : no = gw_corr_lev_occ(ispin)
200 8 : nv = gw_corr_lev_virt(ispin)
201 8 : n_mo_win = no + nv
202 8 : win_start = homo(ispin) - no + 1
203 :
204 : ! ==============================================================================
205 : ! 2. Screened 3-center integrals, weighted with the COHSEX occupation factor
206 : ! (W^c B)^P_mp = Σ_Q W^c_PQ(0) B^Q_mp s_p, s_p = 1/2 sgn(ϵ_p - ϵ_F)
207 : ! s_p = -1/2 for p occupied, s_p = +1/2 for p virtual (COH)
208 : ! ==============================================================================
209 : CALL parallel_gemm(transa="N", transb="N", m=dimen_RI, n=n_mo_win*nmo, k=dimen_RI, &
210 : alpha=1.0_dp, matrix_a=fm_mat_Q, matrix_b=fm_mat_S_gw(ispin), &
211 8 : beta=0.0_dp, matrix_c=fm_mat_S_gw_work(ispin))
212 24 : ALLOCATE (scale_vec(n_mo_win*nmo))
213 162 : DO m = 1, n_mo_win
214 3704 : DO p = 1, nmo
215 3696 : IF (p <= homo(ispin)) THEN
216 616 : scale_vec((m - 1)*nmo + p) = -0.5_dp
217 : ELSE
218 2926 : scale_vec((m - 1)*nmo + p) = 0.5_dp
219 : END IF
220 : END DO
221 : END DO
222 8 : CALL cp_fm_column_scale(fm_mat_S_gw_work(ispin), scale_vec)
223 8 : DEALLOCATE (scale_vec)
224 :
225 : ! ==============================================================================
226 : ! 3. Reorder the compound index for O(N^4) contractions over (p,P)
227 : ! B'_(p,P),m = B^P_mp
228 : ! (W^c B)'_(p,P),m = (W^c B)^P_mp
229 : ! row index (p-1)*dimen_RI + P, P fast: the occupied p (p <= homo) are the
230 : ! first homo*dimen_RI rows. cp_fm_permute maps the input indices
231 : ! (P,1 ; m,p) with reordering [4,1,3,2] to (p,P ; m,1).
232 : ! ==============================================================================
233 8 : NULLIFY (fm_struct_B_re)
234 : CALL cp_fm_struct_create(fm_struct_B_re, template_fmstruct=fm_mat_S_gw(ispin)%matrix_struct, &
235 8 : nrow_global=nmo*dimen_RI, ncol_global=n_mo_win)
236 8 : CALL cp_fm_create(fm_B_re, fm_struct_B_re, name="static_cohsex_B_reordered", set_zero=.TRUE.)
237 8 : CALL cp_fm_create(fm_WB_re, fm_struct_B_re, name="static_cohsex_WB_reordered", set_zero=.TRUE.)
238 8 : CALL cp_fm_struct_release(fm_struct_B_re)
239 8 : CALL cp_fm_permute(fm_B_re, fm_mat_S_gw(ispin), 1.0_dp, 1, nmo, dimen_RI, 1, [4, 1, 3, 2])
240 8 : CALL cp_fm_permute(fm_WB_re, fm_mat_S_gw_work(ispin), 1.0_dp, 1, nmo, dimen_RI, 1, [4, 1, 3, 2])
241 :
242 8 : NULLIFY (fm_struct_win)
243 : CALL cp_fm_struct_create(fm_struct_win, template_fmstruct=fm_mat_S_gw(ispin)%matrix_struct, &
244 8 : nrow_global=n_mo_win, ncol_global=n_mo_win)
245 8 : CALL cp_fm_create(fm_sigma, fm_struct_win)
246 8 : CALL cp_fm_struct_release(fm_struct_win)
247 :
248 48 : ALLOCATE (sigma_c(n_mo_win, n_mo_win), sigma_x(n_mo_win, n_mo_win))
249 :
250 : ! ==============================================================================
251 : ! 4. Correlation self-energy
252 : ! Σ^c_mn = Σ_p s_p Σ_PQ B^P_mp W^c_PQ(0) B^Q_pn
253 : ! = Σ_(p,P) (W^c B)'_(p,P),m B'_(p,P),n
254 : ! ==============================================================================
255 : CALL parallel_gemm(transa="T", transb="N", m=n_mo_win, n=n_mo_win, k=nmo*dimen_RI, &
256 : alpha=1.0_dp, matrix_a=fm_WB_re, matrix_b=fm_B_re, &
257 8 : beta=0.0_dp, matrix_c=fm_sigma)
258 8 : CALL cp_fm_get_submatrix(fm_sigma, sigma_c)
259 :
260 : ! ==============================================================================
261 : ! 5. RI exchange self-energy (RI_SIGMA_X)
262 : ! Σ^x_mn = -Σ_i^occ Σ_P B^P_mi B^P_in = -Σ_(i,P) B'_(i,P),m B'_(i,P),n
263 : ! ==============================================================================
264 8 : sigma_x(:, :) = 0.0_dp
265 8 : IF (mp2_env%ri_g0w0%do_ri_Sigma_x) THEN
266 : CALL parallel_gemm(transa="T", transb="N", m=n_mo_win, n=n_mo_win, k=homo(ispin)*dimen_RI, &
267 : alpha=-1.0_dp, matrix_a=fm_B_re, matrix_b=fm_B_re, &
268 6 : beta=0.0_dp, matrix_c=fm_sigma)
269 6 : CALL cp_fm_get_submatrix(fm_sigma, sigma_x)
270 : END IF
271 :
272 8 : CALL cp_fm_release(fm_sigma)
273 8 : CALL cp_fm_release(fm_B_re)
274 8 : CALL cp_fm_release(fm_WB_re)
275 :
276 : ! ==============================================================================
277 : ! 6. v^xc and Σ^x - v^xc in the MO window
278 : ! v^xc_mn = Σ_μν C_μm v^xc_μν C_νn
279 : ! (Σ^x - v^xc)_mn = Σ_μν C_μm (Σ^x - v^xc)_μν C_νn
280 : ! The AO matrix (Σ^x - v^xc)_μν holds the HFX exchange Σ^x_μν without
281 : ! RI_SIGMA_X and only -v^xc_μν with RI_SIGMA_X, hence in both cases
282 : ! Σ^x_mn = Σ^x,RI_mn + (Σ^x - v^xc)_mn + v^xc_mn (Σ^x,RI = 0 without RI)
283 : ! ==============================================================================
284 48 : ALLOCATE (vxc(n_mo_win, n_mo_win), sigma_x_minus_vxc(n_mo_win, n_mo_win))
285 : CALL ao_to_mo_window(mp2_env%ri_g0w0%matrix_vxc(ispin)%matrix, mo_coeff(ispin), &
286 8 : win_start, n_mo_win, vxc)
287 : CALL ao_to_mo_window(mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix, mo_coeff(ispin), &
288 8 : win_start, n_mo_win, sigma_x_minus_vxc)
289 3464 : sigma_x(:, :) = sigma_x(:, :) + sigma_x_minus_vxc(:, :) + vxc(:, :)
290 :
291 : ! ==============================================================================
292 : ! 7. static COHSEX quasiparticle Hamiltonian in the window
293 : ! H^QP_ij = ϵ_i δ_ij + Σ^x_ij - v^xc_ij + α_IP Σ^c_ij (occ-occ)
294 : ! H^QP_ab = ϵ_a δ_ab + Σ^x_ab - v^xc_ab + α_EA Σ^c_ab (virt-virt)
295 : ! H^QP_ia = 0 (occ-virt)
296 : ! symmetrized, H^QP = (H^QP + H^QP^T)/2
297 : ! ϵ_n^DG = H^QP_nn (diagonal static COHSEX)
298 : ! ==============================================================================
299 32 : ALLOCATE (h_qp(n_mo_win, n_mo_win))
300 8 : h_qp(:, :) = 0.0_dp
301 : h_qp(1:no, 1:no) = sigma_x(1:no, 1:no) - vxc(1:no, 1:no) + &
302 152 : mp2_env%ri_g0w0%cohsex_scale_occ*sigma_c(1:no, 1:no)
303 : h_qp(no + 1:, no + 1:) = sigma_x(no + 1:, no + 1:) - vxc(no + 1:, no + 1:) + &
304 2348 : mp2_env%ri_g0w0%cohsex_scale_virt*sigma_c(no + 1:, no + 1:)
305 6920 : h_qp(:, :) = 0.5_dp*(h_qp(:, :) + TRANSPOSE(h_qp(:, :)))
306 162 : DO n = 1, n_mo_win
307 162 : h_qp(n, n) = h_qp(n, n) + Eigenval(win_start + n - 1, ispin)
308 : END DO
309 :
310 24 : ALLOCATE (e_diag(n_mo_win))
311 162 : DO n = 1, n_mo_win
312 162 : e_diag(n) = h_qp(n, n)
313 : END DO
314 :
315 : ! ϵ^DG of all MOs for a subsequent BSE, scissor shift of the DFT eigenvalue outside the window
316 8 : IF (PRESENT(eigenval_dg)) THEN
317 48 : DO n = 1, nmo
318 48 : IF (n < win_start) THEN
319 0 : eigenval_dg(n, ispin) = Eigenval(n, ispin) + e_diag(1) - Eigenval(win_start, ispin)
320 46 : ELSE IF (n >= win_start + n_mo_win) THEN
321 : eigenval_dg(n, ispin) = Eigenval(n, ispin) + e_diag(n_mo_win) - &
322 0 : Eigenval(win_start + n_mo_win - 1, ispin)
323 : ELSE
324 46 : eigenval_dg(n, ispin) = e_diag(n - win_start + 1)
325 : END IF
326 : END DO
327 : END IF
328 :
329 8 : IF (qp_update) THEN
330 : ! ==============================================================================
331 : ! 8. Canonical static COHSEX (QP_UPDATE): diagonalize the occ-occ and virt-virt blocks
332 : ! Σ_m H^QP_nm U_mk = U_nk ϵ_k^CN
333 : ! ==============================================================================
334 8 : ALLOCATE (e_can(n_mo_win))
335 4 : e_can(:) = 0.0_dp
336 4 : IF (para_env%is_source()) THEN
337 2 : IF (no > 0) CALL diamat_all(h_qp(1:no, 1:no), e_can(1:no))
338 2 : IF (nv > 0) CALL diamat_all(h_qp(no + 1:, no + 1:), e_can(no + 1:))
339 : END IF
340 4 : CALL para_env%bcast(h_qp, para_env%source)
341 4 : CALL para_env%bcast(e_can, para_env%source)
342 :
343 : ! ==============================================================================
344 : ! 9. Energies and MO transformation for all nmo MOs
345 : ! inside the window: ϵ_n^QP = ϵ_n^CN, U_mn from step 8
346 : ! outside the window: U_mn = δ_mn and scissor shift of the DFT eigenvalue
347 : ! ϵ_n^QP = ϵ_n + (ϵ^CN_first - ϵ_first) below the window (occupied)
348 : ! ϵ_n^QP = ϵ_n + (ϵ^CN_last - ϵ_last) above the window (virtual)
349 : ! ==============================================================================
350 24 : ALLOCATE (e_qp(nmo), u_qp(nmo, nmo))
351 4 : shift_occ = e_can(1) - Eigenval(win_start, ispin)
352 4 : shift_virt = e_can(n_mo_win) - Eigenval(win_start + n_mo_win - 1, ispin)
353 96 : DO n = 1, nmo
354 96 : IF (n < win_start) THEN
355 0 : e_qp(n) = Eigenval(n, ispin) + shift_occ
356 92 : ELSE IF (n >= win_start + n_mo_win) THEN
357 0 : e_qp(n) = Eigenval(n, ispin) + shift_virt
358 : ELSE
359 92 : e_qp(n) = e_can(n - win_start + 1)
360 : END IF
361 : END DO
362 4 : u_qp(:, :) = 0.0_dp
363 96 : DO n = 1, nmo
364 96 : u_qp(n, n) = 1.0_dp
365 : END DO
366 2212 : u_qp(win_start:win_start + n_mo_win - 1, win_start:win_start + n_mo_win - 1) = h_qp(:, :)
367 :
368 : ! for a subsequent BSE: ϵ^CN and U of all MOs
369 50 : IF (PRESENT(eigenval_cn)) eigenval_cn(:, ispin) = e_qp(:)
370 1108 : IF (PRESENT(u_qp_all)) u_qp_all(:, :, ispin) = u_qp(:, :)
371 : END IF
372 :
373 : ! ==============================================================================
374 : ! 10. Print per MO: ϵ_n, Σ^c_nn (unscaled), Σ^x_nn, v^xc_nn, ϵ_n^DG and, with QP_UPDATE,
375 : ! ϵ_n^CN; then HOMO, LUMO and gap. ϵ^CN is sorted in each block, i.e.
376 : ! ϵ^CN_HOMO = max over the occ block and ϵ^CN_LUMO = min over the virt block
377 : ! ==============================================================================
378 8 : IF (unit_nr > 0) THEN
379 4 : spin_label = ""
380 4 : IF (nspins == 2 .AND. ispin == 1) spin_label = "Alpha "
381 4 : IF (nspins == 2 .AND. ispin == 2) spin_label = "Beta "
382 4 : WRITE (unit_nr, "(/,T3,A)") REPEAT("-", 77)
383 : ! "static COHSEX" for closed shell, "Alpha static COHSEX" / "Beta static COHSEX" otherwise
384 4 : cohsex_label = TRIM(ADJUSTL(TRIM(spin_label)//" static COHSEX"))
385 4 : WRITE (unit_nr, "(T3,A)") TRIM(cohsex_label)//" quasiparticle energies"
386 4 : WRITE (unit_nr, "(T3,A,F6.3,A,F6.3)") "COHSEX_SCALE_OCC =", mp2_env%ri_g0w0%cohsex_scale_occ, &
387 8 : ", COHSEX_SCALE_VIRT =", mp2_env%ri_g0w0%cohsex_scale_virt
388 4 : WRITE (unit_nr, "(T3,A)") REPEAT("-", 77)
389 4 : IF (qp_update) THEN
390 : WRITE (unit_nr, "(T3,A)") "DG: diagonal of H^QP, CN: canonical, eigenvalues of the "// &
391 2 : "occ-occ and virt-virt blocks"
392 2 : WRITE (unit_nr, "(/,T3,A,T16,A,T27,A,T38,A,T49,A,T60,A,T71,A)") "MO", "e_n^DFT", "Sigma^c_nn", &
393 4 : "Sigma^x_nn", "v^xc_nn", "e_n^DG", "e_n^CN"
394 2 : WRITE (unit_nr, "(T3,A,T16,A,T27,A,T38,A,T49,A,T60,A,T71,A)") "", "(eV)", "(eV)", "(eV)", &
395 4 : "(eV)", "(eV)", "(eV)"
396 : ELSE
397 2 : WRITE (unit_nr, "(T3,A)") "DG: diagonal of H^QP"
398 2 : WRITE (unit_nr, "(/,T3,A,T16,A,T27,A,T38,A,T49,A,T60,A)") "MO", "e_n^DFT", "Sigma^c_nn", &
399 4 : "Sigma^x_nn", "v^xc_nn", "e_n^DG"
400 2 : WRITE (unit_nr, "(T3,A,T16,A,T27,A,T38,A,T49,A,T60,A)") "", "(eV)", "(eV)", "(eV)", &
401 4 : "(eV)", "(eV)"
402 : END IF
403 81 : DO n = 1, n_mo_win
404 77 : IF (n <= no) THEN
405 15 : occ_vir = "occ"
406 : ELSE
407 62 : occ_vir = "vir"
408 : END IF
409 81 : IF (qp_update) THEN
410 46 : WRITE (unit_nr, "(T3,I5,A,A3,A,6F11.3)") win_start + n - 1, " (", occ_vir, ")", &
411 46 : Eigenval(win_start + n - 1, ispin)*evolt, sigma_c(n, n)*evolt, sigma_x(n, n)*evolt, &
412 92 : vxc(n, n)*evolt, e_diag(n)*evolt, e_can(n)*evolt
413 : ELSE
414 31 : WRITE (unit_nr, "(T3,I5,A,A3,A,5F11.3)") win_start + n - 1, " (", occ_vir, ")", &
415 31 : Eigenval(win_start + n - 1, ispin)*evolt, sigma_c(n, n)*evolt, sigma_x(n, n)*evolt, &
416 62 : vxc(n, n)*evolt, e_diag(n)*evolt
417 : END IF
418 : END DO
419 4 : IF (no > 0 .AND. nv > 0) THEN
420 4 : WRITE (unit_nr, "(/,T3,A,T66,F15.3)") TRIM(cohsex_label)//" HOMO, diagonal (eV)", &
421 23 : MAXVAL(e_diag(1:no))*evolt
422 4 : WRITE (unit_nr, "(T3,A,T66,F15.3)") TRIM(cohsex_label)//" LUMO, diagonal (eV)", &
423 70 : MINVAL(e_diag(no + 1:))*evolt
424 4 : WRITE (unit_nr, "(T3,A,T66,F15.3)") TRIM(cohsex_label)//" HOMO-LUMO gap, diagonal (eV)", &
425 85 : (MINVAL(e_diag(no + 1:)) - MAXVAL(e_diag(1:no)))*evolt
426 4 : IF (qp_update) THEN
427 2 : WRITE (unit_nr, "(T3,A,T66,F15.3)") TRIM(cohsex_label)//" HOMO, canonical (eV)", &
428 12 : MAXVAL(e_can(1:no))*evolt
429 2 : WRITE (unit_nr, "(T3,A,T66,F15.3)") TRIM(cohsex_label)//" LUMO, canonical (eV)", &
430 42 : MINVAL(e_can(no + 1:))*evolt
431 2 : WRITE (unit_nr, "(T3,A,T66,F15.3)") TRIM(cohsex_label)//" HOMO-LUMO gap, canonical (eV)", &
432 50 : (MINVAL(e_can(no + 1:)) - MAXVAL(e_can(1:no)))*evolt
433 : END IF
434 : END IF
435 4 : WRITE (unit_nr, "(T3,A)") REPEAT("-", 77)
436 : END IF
437 :
438 8 : IF (qp_update) THEN
439 : ! ==============================================================================
440 : ! 11. Canonical static COHSEX MOs C^QP = C U via DFT%PRINT%MO and DFT%PRINT%MO_CUBES
441 : ! ==============================================================================
442 4 : CALL print_static_cohsex_mos(qs_env, mo_coeff(ispin), u_qp, e_qp, homo(ispin), nspins, ispin)
443 4 : DEALLOCATE (e_can, e_qp, u_qp)
444 : END IF
445 :
446 24 : DEALLOCATE (sigma_c, sigma_x, sigma_x_minus_vxc, vxc, h_qp, e_diag)
447 :
448 : END DO
449 :
450 16 : DO ispin = 1, nspins
451 8 : CALL dbcsr_release_p(mp2_env%ri_g0w0%matrix_vxc(ispin)%matrix)
452 16 : CALL dbcsr_release_p(mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc(ispin)%matrix)
453 : END DO
454 8 : DEALLOCATE (mp2_env%ri_g0w0%matrix_vxc, mp2_env%ri_g0w0%matrix_sigma_x_minus_vxc)
455 :
456 8 : CALL timestop(handle)
457 :
458 16 : END SUBROUTINE static_cohsex_qp_energies
459 :
460 : ! **************************************************************************************************
461 : !> \brief MO-window block of an AO matrix
462 : !> A_mn = Σ_μν C_μm A_μν C_νn, m,n = win_start, ..., win_start + n_mo_win - 1
463 : !> \param matrix_ao ...
464 : !> \param mo_coeff ...
465 : !> \param win_start ...
466 : !> \param n_mo_win ...
467 : !> \param mat_mo ...
468 : ! **************************************************************************************************
469 16 : SUBROUTINE ao_to_mo_window(matrix_ao, mo_coeff, win_start, n_mo_win, mat_mo)
470 : TYPE(dbcsr_type), POINTER :: matrix_ao
471 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
472 : INTEGER, INTENT(IN) :: win_start, n_mo_win
473 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: mat_mo
474 :
475 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ao_to_mo_window'
476 :
477 : INTEGER :: handle, nao
478 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
479 : TYPE(cp_fm_type) :: fm_ao, fm_ao_mo, fm_mo_mo
480 :
481 16 : CALL timeset(routineN, handle)
482 :
483 16 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao)
484 :
485 16 : NULLIFY (fm_struct)
486 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=mo_coeff%matrix_struct, &
487 16 : nrow_global=nao, ncol_global=nao)
488 16 : CALL cp_fm_create(fm_ao, fm_struct)
489 16 : CALL cp_fm_struct_release(fm_struct)
490 16 : CALL copy_dbcsr_to_fm(matrix_ao, fm_ao)
491 :
492 : ! (A C)_μn = Σ_ν A_μν C_νn, only the window columns n of C
493 :
494 16 : NULLIFY (fm_struct)
495 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=mo_coeff%matrix_struct, &
496 16 : nrow_global=nao, ncol_global=n_mo_win)
497 16 : CALL cp_fm_create(fm_ao_mo, fm_struct)
498 16 : CALL cp_fm_struct_release(fm_struct)
499 : CALL parallel_gemm(transa="N", transb="N", m=nao, n=n_mo_win, k=nao, alpha=1.0_dp, &
500 : matrix_a=fm_ao, matrix_b=mo_coeff, beta=0.0_dp, matrix_c=fm_ao_mo, &
501 16 : b_first_col=win_start)
502 16 : CALL cp_fm_release(fm_ao)
503 :
504 : ! A_mn = Σ_μ C_μm (A C)_μn
505 16 : NULLIFY (fm_struct)
506 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=mo_coeff%matrix_struct, &
507 16 : nrow_global=n_mo_win, ncol_global=n_mo_win)
508 16 : CALL cp_fm_create(fm_mo_mo, fm_struct)
509 16 : CALL cp_fm_struct_release(fm_struct)
510 : CALL parallel_gemm(transa="T", transb="N", m=n_mo_win, n=n_mo_win, k=nao, alpha=1.0_dp, &
511 : matrix_a=mo_coeff, matrix_b=fm_ao_mo, beta=0.0_dp, matrix_c=fm_mo_mo, &
512 16 : a_first_col=win_start)
513 16 : CALL cp_fm_release(fm_ao_mo)
514 :
515 16 : CALL cp_fm_get_submatrix(fm_mo_mo, mat_mo)
516 16 : CALL cp_fm_release(fm_mo_mo)
517 :
518 16 : CALL timestop(handle)
519 :
520 16 : END SUBROUTINE ao_to_mo_window
521 :
522 : ! **************************************************************************************************
523 : !> \brief Print the canonical static COHSEX MOs with the DFT%PRINT%MO and DFT%PRINT%MO_CUBES
524 : !> \param qs_env ...
525 : !> \param mo_coeff ...
526 : !> \param u_qp ...
527 : !> \param e_qp ...
528 : !> \param homo ...
529 : !> \param nspins ...
530 : !> \param ispin ...
531 : ! **************************************************************************************************
532 4 : SUBROUTINE print_static_cohsex_mos(qs_env, mo_coeff, u_qp, e_qp, homo, nspins, ispin)
533 : TYPE(qs_environment_type), POINTER :: qs_env
534 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
535 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: u_qp
536 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: e_qp
537 : INTEGER, INTENT(IN) :: homo, nspins, ispin
538 :
539 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_static_cohsex_mos'
540 :
541 : CHARACTER(LEN=default_path_length) :: fname, title
542 : INTEGER :: first_mo, handle, ivector, last_mo, nao, &
543 : nelectron, nhomo, nlumo, nmo, unit_nr
544 4 : INTEGER, DIMENSION(:), POINTER :: stride
545 : LOGICAL :: print_coeff, print_cubes, &
546 : print_energies, print_mo, print_occ
547 : REAL(KIND=dp) :: maxocc
548 4 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
549 : TYPE(cell_type), POINTER :: cell
550 : TYPE(cp_fm_type) :: fm_qp_mo
551 : TYPE(cp_logger_type), POINTER :: logger
552 : TYPE(dft_control_type), POINTER :: dft_control
553 : TYPE(mo_set_type) :: qp_mo_set
554 : TYPE(particle_list_type), POINTER :: particles
555 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
556 : TYPE(pw_c1d_gs_type) :: wf_g
557 : TYPE(pw_env_type), POINTER :: pw_env
558 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
559 : TYPE(pw_r3d_rs_type) :: wf_r
560 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
561 : TYPE(qs_subsys_type), POINTER :: subsys
562 : TYPE(section_vals_type), POINTER :: dft_section, input
563 :
564 4 : CALL timeset(routineN, handle)
565 :
566 4 : CALL get_qs_env(qs_env, input=input)
567 4 : dft_section => section_vals_get_subs_vals(input, "DFT")
568 4 : logger => cp_get_default_logger()
569 :
570 4 : CALL section_vals_val_get(dft_section, "PRINT%MO%ENERGIES", l_val=print_energies)
571 4 : CALL section_vals_val_get(dft_section, "PRINT%MO%COEFFICIENTS", l_val=print_coeff)
572 4 : CALL section_vals_val_get(dft_section, "PRINT%MO%OCCUPATION_NUMBERS", l_val=print_occ)
573 : print_mo = (cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO") /= 0) .AND. &
574 4 : (print_energies .OR. print_occ .OR. print_coeff)
575 4 : print_cubes = (cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO_CUBES") /= 0)
576 :
577 4 : IF (print_mo .OR. print_cubes) THEN
578 :
579 0 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_global=nmo)
580 :
581 : ! C^QP_μn = Σ_m C_μm U_mn
582 0 : CALL get_static_cohsex_mo_coeff(mo_coeff, u_qp, fm_qp_mo)
583 :
584 : ! MO energies ϵ_n^QP, occupations f_n (maxocc for n <= homo) and coefficients C^QP_μn
585 0 : IF (print_mo) THEN
586 0 : IF (nspins == 1) THEN
587 0 : maxocc = 2.0_dp
588 : ELSE
589 0 : maxocc = 1.0_dp
590 : END IF
591 0 : nelectron = NINT(maxocc*REAL(homo, KIND=dp))
592 : CALL allocate_mo_set(qp_mo_set, nao=nao, nmo=nmo, nelectron=nelectron, &
593 : n_el_f=REAL(nelectron, KIND=dp), maxocc=maxocc, &
594 0 : flexible_electron_count=0.0_dp)
595 0 : CALL init_mo_set(qp_mo_set, fm_struct=fm_qp_mo%matrix_struct, name="static_cohsex")
596 0 : qp_mo_set%eigenvalues(:) = e_qp(:)
597 0 : CALL cp_fm_to_fm(fm_qp_mo, qp_mo_set%mo_coeff)
598 :
599 0 : CALL get_qs_env(qs_env, particle_set=particle_set, qs_kind_set=qs_kind_set)
600 0 : IF (nspins == 1) THEN
601 : CALL write_mo_set_to_output_unit(qp_mo_set, qs_kind_set, particle_set, dft_section, &
602 0 : 4, 0, final_mos=.TRUE.)
603 0 : ELSE IF (ispin == 1) THEN
604 : CALL write_mo_set_to_output_unit(qp_mo_set, qs_kind_set, particle_set, dft_section, &
605 0 : 4, 0, final_mos=.TRUE., spin="ALPHA")
606 : ELSE
607 : CALL write_mo_set_to_output_unit(qp_mo_set, qs_kind_set, particle_set, dft_section, &
608 0 : 4, 0, final_mos=.TRUE., spin="BETA")
609 : END IF
610 0 : CALL deallocate_mo_set(qp_mo_set)
611 : END IF
612 :
613 : ! cube files of ψ_n(r) = Σ_μ C^QP_μn ϕ_μ(r) for n = homo-NHOMO+1, ..., homo+NLUMO
614 0 : IF (print_cubes) THEN
615 0 : CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%NHOMO", i_val=nhomo)
616 0 : CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%NLUMO", i_val=nlumo)
617 0 : CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%STRIDE", i_vals=stride)
618 0 : IF (nhomo == -1) nhomo = homo
619 0 : IF (nlumo == -1) nlumo = nmo - homo
620 0 : first_mo = MAX(1, homo - nhomo + 1)
621 0 : last_mo = MIN(nmo, homo + nlumo)
622 :
623 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
624 : cell=cell, dft_control=dft_control, particle_set=particle_set, &
625 0 : pw_env=pw_env, subsys=subsys)
626 0 : CALL qs_subsys_get(subsys, particles=particles)
627 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
628 0 : CALL auxbas_pw_pool%create_pw(wf_r)
629 0 : CALL auxbas_pw_pool%create_pw(wf_g)
630 :
631 0 : DO ivector = first_mo, last_mo
632 : CALL calculate_wavefunction(fm_qp_mo, ivector, wf_r, wf_g, atomic_kind_set, &
633 0 : qs_kind_set, cell, dft_control, particle_set, pw_env)
634 0 : WRITE (fname, '(A,I5.5,A,I1)') "static_cohsex_WFN_", ivector, "_", ispin
635 : unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%MO_CUBES", &
636 : extension=".cube", middle_name=TRIM(fname), &
637 0 : file_position="REWIND", log_filename=.FALSE.)
638 0 : WRITE (title, *) "static COHSEX WAVEFUNCTION ", ivector, " spin ", ispin, &
639 0 : " i.e. HOMO ", ivector - homo
640 0 : CALL cp_pw_to_cube(wf_r, unit_nr, title, particles, stride=stride)
641 0 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%MO_CUBES")
642 : END DO
643 :
644 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
645 0 : CALL auxbas_pw_pool%give_back_pw(wf_g)
646 : END IF
647 :
648 0 : CALL cp_fm_release(fm_qp_mo)
649 : END IF
650 :
651 4 : CALL timestop(handle)
652 :
653 4 : END SUBROUTINE print_static_cohsex_mos
654 :
655 : ! **************************************************************************************************
656 : !> \brief Canonical static COHSEX MO coefficients
657 : !> C^QP_μn = Σ_m C_μm U_mn
658 : !> \param mo_coeff ...
659 : !> \param u_qp ...
660 : !> \param fm_qp_mo ...
661 : ! **************************************************************************************************
662 2 : SUBROUTINE get_static_cohsex_mo_coeff(mo_coeff, u_qp, fm_qp_mo)
663 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
664 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: u_qp
665 : TYPE(cp_fm_type), INTENT(OUT) :: fm_qp_mo
666 :
667 : INTEGER :: nao, nmo
668 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
669 : TYPE(cp_fm_type) :: fm_u
670 :
671 2 : CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_global=nmo)
672 2 : NULLIFY (fm_struct)
673 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=mo_coeff%matrix_struct, &
674 2 : nrow_global=nmo, ncol_global=nmo)
675 2 : CALL cp_fm_create(fm_u, fm_struct)
676 2 : CALL cp_fm_struct_release(fm_struct)
677 2 : CALL cp_fm_set_submatrix(fm_u, u_qp)
678 2 : CALL cp_fm_create(fm_qp_mo, mo_coeff%matrix_struct)
679 : CALL parallel_gemm(transa="N", transb="N", m=nao, n=nmo, k=nmo, alpha=1.0_dp, &
680 2 : matrix_a=mo_coeff, matrix_b=fm_u, beta=0.0_dp, matrix_c=fm_qp_mo)
681 2 : CALL cp_fm_release(fm_u)
682 :
683 2 : END SUBROUTINE get_static_cohsex_mo_coeff
684 :
685 : ! **************************************************************************************************
686 : !> \brief Rotate an RI 3-center slab to the canonical static COHSEX MOs
687 : !> \param fm_slab ...
688 : !> \param dimen_RI ...
689 : !> \param n_outer ...
690 : !> \param n_inner ...
691 : !> \param u_outer ...
692 : !> \param u_inner ...
693 : ! **************************************************************************************************
694 6 : SUBROUTINE rotate_bse_slab(fm_slab, dimen_RI, n_outer, n_inner, u_outer, u_inner)
695 : TYPE(cp_fm_type), INTENT(IN) :: fm_slab
696 : INTEGER, INTENT(IN) :: dimen_RI, n_outer, n_inner
697 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: u_outer, u_inner
698 :
699 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rotate_bse_slab'
700 :
701 : INTEGER :: handle
702 :
703 6 : CALL timeset(routineN, handle)
704 :
705 : ! (a) inner index: (P ; o,i) -> (P,o ; i), P fast, reordering [3,1,2,4] and back [3,2,1,4]
706 : CALL rotate_one_index(fm_slab, dimen_RI, n_outer, n_inner, u_inner, [3, 1, 2, 4], [3, 2, 1, 4], &
707 6 : n_inner)
708 : ! (b) outer index: (P ; o,i) -> (P,i ; o), P fast, reordering [4,1,2,3] and back [3,2,4,1]
709 : CALL rotate_one_index(fm_slab, dimen_RI, n_inner, n_outer, u_outer, [4, 1, 2, 3], [3, 2, 4, 1], &
710 6 : n_inner)
711 :
712 6 : CALL timestop(handle)
713 :
714 6 : END SUBROUTINE rotate_bse_slab
715 :
716 : ! **************************************************************************************************
717 : !> \brief rotate one MO index of the compound column index
718 : !> \param fm_slab ...
719 : !> \param dimen_RI ...
720 : !> \param n_spect ...
721 : !> \param n_rot ...
722 : !> \param u_rot ...
723 : !> \param reord_to ...
724 : !> \param reord_back ...
725 : !> \param n_inner size ...
726 : ! **************************************************************************************************
727 12 : SUBROUTINE rotate_one_index(fm_slab, dimen_RI, n_spect, n_rot, u_rot, reord_to, reord_back, &
728 : n_inner)
729 : TYPE(cp_fm_type), INTENT(IN) :: fm_slab
730 : INTEGER, INTENT(IN) :: dimen_RI, n_spect, n_rot
731 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: u_rot
732 : INTEGER, DIMENSION(4) :: reord_to, reord_back
733 : INTEGER, INTENT(IN) :: n_inner
734 :
735 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
736 : TYPE(cp_fm_type) :: fm_slab_out, fm_t, fm_t_rot, fm_u
737 :
738 12 : NULLIFY (fm_struct)
739 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=fm_slab%matrix_struct, &
740 12 : nrow_global=n_spect*dimen_RI, ncol_global=n_rot)
741 12 : CALL cp_fm_create(fm_t, fm_struct, name="static_cohsex_slab_reordered", set_zero=.TRUE.)
742 12 : CALL cp_fm_create(fm_t_rot, fm_struct, name="static_cohsex_slab_rotated", set_zero=.TRUE.)
743 12 : CALL cp_fm_struct_release(fm_struct)
744 12 : CALL cp_fm_permute(fm_t, fm_slab, 1.0_dp, 1, n_inner, dimen_RI, n_rot, reord_to)
745 :
746 12 : NULLIFY (fm_struct)
747 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=fm_slab%matrix_struct, &
748 12 : nrow_global=n_rot, ncol_global=n_rot)
749 12 : CALL cp_fm_create(fm_u, fm_struct, name="static_cohsex_U")
750 12 : CALL cp_fm_struct_release(fm_struct)
751 12 : CALL cp_fm_set_submatrix(fm_u, u_rot)
752 : CALL parallel_gemm(transa="N", transb="N", m=n_spect*dimen_RI, n=n_rot, k=n_rot, alpha=1.0_dp, &
753 12 : matrix_a=fm_t, matrix_b=fm_u, beta=0.0_dp, matrix_c=fm_t_rot)
754 12 : CALL cp_fm_release(fm_u)
755 12 : CALL cp_fm_release(fm_t)
756 :
757 12 : CALL cp_fm_create(fm_slab_out, fm_slab%matrix_struct, name="static_cohsex_slab_back", set_zero=.TRUE.)
758 12 : CALL cp_fm_permute(fm_slab_out, fm_t_rot, 1.0_dp, dimen_RI, n_rot, 1, n_inner, reord_back)
759 12 : CALL cp_fm_to_fm(fm_slab_out, fm_slab)
760 12 : CALL cp_fm_release(fm_slab_out)
761 12 : CALL cp_fm_release(fm_t_rot)
762 :
763 12 : END SUBROUTINE rotate_one_index
764 :
765 : END MODULE rpa_gw_static_cohsex
|
