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 : MODULE qs_tddfpt2_bse_utils
9 : USE cp_blacs_env, ONLY: cp_blacs_env_type
10 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
11 : dbcsr_p_type,&
12 : dbcsr_release,&
13 : dbcsr_set,&
14 : dbcsr_type
15 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
16 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
17 : copy_fm_to_dbcsr,&
18 : cp_dbcsr_sm_fm_multiply
19 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
20 : cp_fm_matvec,&
21 : cp_fm_scale_and_add
22 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
23 : cp_fm_struct_release,&
24 : cp_fm_struct_type
25 : USE cp_fm_types, ONLY: &
26 : cp_fm_create, cp_fm_get_info, cp_fm_release, cp_fm_set_all, cp_fm_set_element, &
27 : cp_fm_to_fm, cp_fm_to_fm_submat, cp_fm_type, cp_fm_vectorssum, cp_fm_write_formatted
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_get_default_io_unit,&
30 : cp_logger_type
31 : USE exstates_types, ONLY: excited_energy_type
32 : USE kinds, ONLY: dp
33 : USE message_passing, ONLY: mp_para_env_type
34 : USE parallel_gemm_api, ONLY: parallel_gemm
35 : USE qs_environment_types, ONLY: get_qs_env,&
36 : qs_environment_type
37 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
38 : USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos
39 : #include "./base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_bse_utils'
46 :
47 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
48 : ! number of first derivative components (3: d/dx, d/dy, d/dz)
49 : INTEGER, PARAMETER, PRIVATE :: nderivs = 3
50 : INTEGER, PARAMETER, PRIVATE :: maxspins = 2
51 :
52 : PUBLIC:: tddfpt_apply_bse_debug, tddfpt_apply_bse
53 : PUBLIC:: zeroth_order_gw
54 :
55 : CONTAINS
56 : ! **************************************************************************************************
57 : !> \brief ...
58 : !> \param qs_env ...
59 : !> \param Aop_evects ...
60 : !> \param evects ...
61 : !> \param S_evects ...
62 : !> \param gs_mos ...
63 : !> \param matrix_s ...
64 : !> \param matrix_ks ...
65 : ! **************************************************************************************************
66 20 : SUBROUTINE zeroth_order_gw(qs_env, Aop_evects, evects, S_evects, gs_mos, matrix_s, matrix_ks)
67 : TYPE(qs_environment_type), POINTER :: qs_env
68 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects
69 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects, S_evects
70 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
71 : INTENT(in) :: gs_mos
72 : TYPE(dbcsr_type), INTENT(in), POINTER :: matrix_s
73 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_ks
74 :
75 : CHARACTER(LEN=*), PARAMETER :: routineN = 'zeroth_order_gw'
76 :
77 : INTEGER :: handle, i, ispin, ivect, j, nactive, &
78 : nao, nmo, nspins, nvects, occ, virt
79 20 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: gw_occ, gw_virt
80 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
81 : TYPE(cp_fm_struct_type), POINTER :: fmstruct, matrix_struct
82 : TYPE(cp_fm_type) :: fms, hevec, matrixtmp, matrixtmp2, &
83 : matrixtmp3, Sweighted_vect, &
84 : weighted_vect
85 : TYPE(dbcsr_type) :: matrixf
86 : TYPE(excited_energy_type), POINTER :: ex_env
87 : TYPE(mp_para_env_type), POINTER :: para_env
88 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
89 20 : POINTER :: sab_orb
90 :
91 20 : CALL timeset(routineN, handle)
92 :
93 20 : NULLIFY (ex_env, sab_orb)
94 20 : CALL get_qs_env(qs_env, exstate_env=ex_env, sab_orb=sab_orb)
95 :
96 20 : nspins = SIZE(matrix_ks, 1)
97 20 : nspins = SIZE(evects, 1)
98 20 : nvects = SIZE(evects, 2)
99 :
100 40 : DO ispin = 1, nspins
101 :
102 20 : CPASSERT(.NOT. ASSOCIATED(gs_mos(ispin)%evals_occ_matrix))
103 :
104 20 : CALL dbcsr_create(matrixf, template=matrix_s)
105 20 : nmo = SIZE(ex_env%gw_eigen)
106 : CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
107 20 : nrow_global=nao, ncol_global=nactive)
108 20 : NULLIFY (blacs_env, para_env)
109 20 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
110 :
111 20 : occ = SIZE(gs_mos(ispin)%evals_occ)
112 20 : nactive = gs_mos(ispin)%nmo_active
113 20 : nmo = SIZE(ex_env%gw_eigen)
114 20 : virt = SIZE(gs_mos(ispin)%evals_virt)
115 20 : NULLIFY (fmstruct)
116 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
117 20 : context=blacs_env, nrow_global=virt, ncol_global=virt)
118 20 : CALL cp_fm_create(matrixtmp, fmstruct)
119 20 : CALL cp_fm_struct_release(fmstruct)
120 :
121 20 : NULLIFY (fmstruct)
122 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
123 20 : context=blacs_env, nrow_global=virt, ncol_global=nao)
124 20 : CALL cp_fm_create(matrixtmp2, fmstruct)
125 20 : CALL cp_fm_struct_release(fmstruct)
126 :
127 20 : NULLIFY (fmstruct)
128 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
129 20 : context=blacs_env, nrow_global=nao, ncol_global=nao)
130 20 : CALL cp_fm_create(matrixtmp3, fmstruct)
131 20 : CALL cp_fm_create(fms, fmstruct)
132 20 : CALL cp_fm_struct_release(fmstruct)
133 20 : CALL cp_dbcsr_alloc_block_from_nbl(matrixf, sab_orb)
134 :
135 : !--add virt eigenvalues
136 20 : CALL dbcsr_set(matrixf, 0.0_dp)
137 20 : CALL cp_fm_create(weighted_vect, gs_mos(ispin)%mos_virt%matrix_struct)
138 20 : CALL cp_fm_create(Sweighted_vect, gs_mos(ispin)%mos_virt%matrix_struct)
139 20 : CALL cp_fm_to_fm(gs_mos(ispin)%mos_virt, weighted_vect)
140 20 : CALL copy_dbcsr_to_fm(matrix_s, fms)
141 :
142 60 : ALLOCATE (gw_virt(virt))
143 60 : ALLOCATE (gw_occ(nactive))
144 400 : gw_virt(1:virt) = ex_env%gw_eigen(occ + 1:nmo)
145 100 : DO i = 1, nactive
146 80 : j = gs_mos(ispin)%index_active(i)
147 100 : gw_occ(i) = ex_env%gw_eigen(j)
148 : END DO
149 :
150 20 : CALL cp_fm_set_all(matrixtmp, 0.0_dp)
151 400 : DO i = 1, virt
152 400 : CALL cp_fm_set_element(matrixtmp, i, i, gw_virt(i))
153 : END DO
154 20 : DEALLOCATE (gw_virt)
155 20 : CALL parallel_gemm('N', 'N', nao, virt, nao, 1.0_dp, fms, weighted_vect, 0.0_dp, Sweighted_vect)
156 20 : CALL parallel_gemm('N', 'T', virt, nao, virt, 1.0_dp, matrixtmp, Sweighted_vect, 0.0_dp, matrixtmp2)
157 20 : CALL parallel_gemm('N', 'N', nao, nao, virt, 1.0_dp, Sweighted_vect, matrixtmp2, 0.0_dp, matrixtmp3)
158 20 : CALL copy_fm_to_dbcsr(matrixtmp3, matrixf)
159 :
160 20 : CALL cp_fm_release(weighted_vect)
161 20 : CALL cp_fm_release(Sweighted_vect)
162 20 : CALL cp_fm_release(fmS)
163 : !--add occ eigenvalues
164 : CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
165 20 : nrow_global=nao, ncol_global=nactive)
166 20 : CALL cp_fm_create(hevec, matrix_struct)
167 :
168 212 : DO ivect = 1, nvects
169 : CALL cp_dbcsr_sm_fm_multiply(matrixf, evects(ispin, ivect), &
170 : Aop_evects(ispin, ivect), ncol=nactive, &
171 192 : alpha=1.0_dp, beta=1.0_dp)
172 :
173 192 : CALL cp_fm_to_fm(S_evects(ispin, ivect), hevec)
174 192 : CALL cp_fm_column_scale(hevec, gw_occ)
175 :
176 212 : CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect), -1.0_dp, hevec)
177 : END DO !ivect
178 20 : DEALLOCATE (gw_occ)
179 :
180 20 : CALL cp_fm_release(matrixtmp)
181 20 : CALL cp_fm_release(matrixtmp2)
182 20 : CALL cp_fm_release(matrixtmp3)
183 :
184 20 : CALL dbcsr_release(matrixf)
185 160 : CALL cp_fm_release(hevec)
186 : END DO !ispin
187 :
188 20 : virt = SIZE(Aop_evects, 2)
189 20 : CALL timestop(handle)
190 :
191 40 : END SUBROUTINE zeroth_order_gw
192 :
193 : ! **************************************************************************************************
194 : !> \brief Update action of TDDFPT operator on trial vectors by adding BSE W term.
195 : !> \brief debug version
196 : !> \param Aop_evects ...
197 : !> \param evects ...
198 : !> \param gs_mos ...
199 : !> \param qs_env ...
200 : ! **************************************************************************************************
201 0 : SUBROUTINE tddfpt_apply_bse_debug(Aop_evects, evects, gs_mos, qs_env)
202 :
203 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects
204 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects
205 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
206 : INTENT(in) :: gs_mos
207 : TYPE(qs_environment_type), POINTER :: qs_env
208 :
209 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_bse_debug'
210 :
211 : INTEGER :: a_nao_col, a_virt_col, b_nao_col, c_virt_col, handle, i_occ_row, i_row_global, &
212 : ii, iounit, ispin, ivect, j_col_global, j_occ_row, jj, k_occ_col, mu_col_global, nao, &
213 : ncol_block, ncol_block_bse, ncol_block_cs, ncol_local, ncol_local_bse, ncol_local_cs, &
214 : nrow_block, nrow_block_bse, nrow_block_cs, nrow_local, nrow_local_bse, nrow_local_cs, &
215 : nspins, nvects, nvirt
216 : INTEGER, DIMENSION(2) :: nactive
217 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, col_indices_bse, &
218 0 : col_indices_cs, row_indices, &
219 0 : row_indices_bse, row_indices_cs
220 : REAL(KIND=dp) :: alpha
221 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
222 0 : POINTER :: my_block, my_bse_w_matrix_MO, my_CSvirt
223 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
224 : TYPE(cp_fm_struct_type), POINTER :: fmstruct, matrix_struct
225 : TYPE(cp_fm_type) :: CSvirt, fms, WXaoao, WXmat2, WXvirtao
226 : TYPE(cp_fm_type), POINTER :: bse_w_matrix_MO
227 : TYPE(cp_logger_type), POINTER :: logger
228 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
229 : TYPE(excited_energy_type), POINTER :: ex_env
230 : TYPE(mp_para_env_type), POINTER :: para_env
231 :
232 0 : NULLIFY (logger) !get output_unit
233 0 : logger => cp_get_default_logger()
234 0 : iounit = cp_logger_get_default_io_unit(logger)
235 :
236 0 : CALL timeset(routineN, handle)
237 :
238 0 : nspins = SIZE(evects, 1)
239 0 : nvects = SIZE(evects, 2)
240 : IF (nspins > 1) THEN
241 : alpha = 1.0_dp
242 : ELSE
243 : alpha = 2.0_dp
244 : END IF
245 0 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
246 0 : DO ispin = 1, nspins
247 0 : CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
248 : END DO
249 :
250 0 : NULLIFY (ex_env, para_env, blacs_env, matrix_s)
251 : CALL get_qs_env(qs_env, exstate_env=ex_env, para_env=para_env, blacs_env=blacs_env, &
252 0 : matrix_s=matrix_s)
253 :
254 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
255 0 : context=blacs_env, nrow_global=nao, ncol_global=nao)
256 0 : CALL cp_fm_create(fms, fmstruct)
257 0 : CALL cp_fm_struct_release(fmstruct)
258 0 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, fms)
259 :
260 0 : NULLIFY (bse_w_matrix_MO)
261 0 : bse_w_matrix_MO => ex_env%bse_w_matrix_MO(1, 1)
262 :
263 0 : DO ivect = 1, nvects
264 0 : DO ispin = 1, nspins
265 0 : NULLIFY (matrix_struct, fmstruct)
266 : CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
267 0 : nrow_global=nao, ncol_global=nactive(ispin))
268 0 : nvirt = SIZE(gs_mos(ispin)%evals_virt)
269 :
270 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
271 0 : context=blacs_env, nrow_global=nvirt, ncol_global=nao)
272 0 : CALL cp_fm_create(CSvirt, fmstruct)
273 0 : CALL cp_fm_struct_release(fmstruct)
274 :
275 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
276 : context=blacs_env, nrow_global=nactive(ispin)*nactive(ispin), &
277 0 : ncol_global=nvirt*nao)
278 0 : CALL cp_fm_create(WXvirtao, fmstruct)
279 0 : CALL cp_fm_struct_release(fmstruct)
280 0 : CALL cp_fm_set_all(WXvirtao, 0.0_dp)
281 :
282 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
283 : context=blacs_env, nrow_global=nactive(ispin)*nactive(ispin), &
284 0 : ncol_global=nao*nao)
285 0 : CALL cp_fm_create(WXaoao, fmstruct)
286 0 : CALL cp_fm_struct_release(fmstruct)
287 0 : CALL cp_fm_set_all(WXaoao, 0.0_dp)
288 :
289 0 : CALL parallel_gemm('T', 'N', nvirt, nao, nao, 1.0_dp, gs_mos(ispin)%mos_virt, fms, 0.0_dp, CSvirt)
290 0 : NULLIFY (row_indices, col_indices)
291 : CALL cp_fm_get_info(matrix=WXvirtao, nrow_local=nrow_local, ncol_local=ncol_local, &
292 : row_indices=row_indices, col_indices=col_indices, &
293 0 : nrow_block=nrow_block, ncol_block=ncol_block, local_data=my_block)
294 : CALL cp_fm_get_info(matrix=CSvirt, nrow_local=nrow_local_cs, ncol_local=ncol_local_cs, &
295 : row_indices=row_indices_cs, col_indices=col_indices_cs, &
296 0 : nrow_block=nrow_block_cs, ncol_block=ncol_block_cs, local_data=my_CSvirt)
297 : CALL cp_fm_get_info(matrix=bse_w_matrix_MO, nrow_local=nrow_local_bse, ncol_local=ncol_local_bse, &
298 : row_indices=row_indices_bse, col_indices=col_indices_bse, &
299 0 : nrow_block=nrow_block_bse, ncol_block=ncol_block_bse, local_data=my_bse_w_matrix_MO)
300 :
301 0 : CALL cp_fm_set_all(WXvirtao, 0.0_dp)
302 :
303 0 : DO ii = 1, nrow_local
304 0 : i_row_global = row_indices(ii)
305 0 : DO jj = 1, ncol_local
306 0 : j_col_global = col_indices(jj)
307 :
308 0 : i_occ_row = (i_row_global - 1)/nactive(ispin) + 1
309 0 : j_occ_row = MOD(i_row_global - 1, nactive(ispin)) + 1
310 0 : a_virt_col = (j_col_global - 1)/nao + 1
311 0 : b_nao_col = MOD(j_col_global - 1, nao) + 1
312 :
313 0 : DO c_virt_col = 1, nvirt
314 0 : mu_col_global = (a_virt_col - 1)*nvirt + c_virt_col
315 :
316 : WXvirtao%local_data(i_row_global, j_col_global) = WXvirtao%local_data(i_row_global, j_col_global) + &
317 0 : bse_w_matrix_MO%local_data(i_row_global, mu_col_global)*CSvirt%local_data(c_virt_col, b_nao_col)
318 :
319 : END DO
320 : END DO
321 : END DO
322 :
323 0 : NULLIFY (row_indices, col_indices) ! redefine indices
324 : CALL cp_fm_get_info(matrix=WXaoao, nrow_local=nrow_local, ncol_local=ncol_local, &
325 : row_indices=row_indices, col_indices=col_indices, &
326 0 : nrow_block=nrow_block, ncol_block=ncol_block)
327 :
328 0 : CALL cp_fm_set_all(WXaoao, 0.0_dp)
329 0 : DO ii = 1, nrow_local
330 0 : i_row_global = row_indices(ii)
331 0 : DO jj = 1, ncol_local
332 0 : j_col_global = col_indices(jj)
333 :
334 0 : i_occ_row = (i_row_global - 1)/nactive(ispin) + 1
335 0 : j_occ_row = MOD(i_row_global - 1, nactive(ispin)) + 1
336 0 : a_nao_col = (j_col_global - 1)/nao + 1
337 0 : b_nao_col = MOD(j_col_global - 1, nao) + 1
338 :
339 0 : DO k_occ_col = 1, nvirt
340 0 : mu_col_global = (k_occ_col - 1)*nao + a_nao_col
341 :
342 : WXaoao%local_data(i_row_global, j_col_global) = WXaoao%local_data(i_row_global, j_col_global) + &
343 0 : WXvirtao%local_data(i_row_global, mu_col_global)*CSvirt%local_data(k_occ_col, b_nao_col)
344 :
345 : END DO
346 : END DO
347 : END DO
348 :
349 0 : CALL cp_fm_release(WXvirtao)
350 0 : CALL cp_fm_release(CSvirt)
351 :
352 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
353 0 : context=blacs_env, nrow_global=nao, ncol_global=nactive(ispin))
354 0 : CALL cp_fm_create(WXmat2, fmstruct)
355 0 : CALL cp_fm_struct_release(fmstruct)
356 0 : CALL cp_fm_set_all(WXmat2, 0.0_dp)
357 :
358 0 : DO ii = 1, nrow_local
359 0 : i_row_global = row_indices(ii)
360 0 : DO jj = 1, ncol_local
361 0 : j_col_global = col_indices(jj)
362 :
363 0 : i_occ_row = (i_row_global - 1)/nactive(ispin) + 1
364 0 : j_occ_row = MOD(i_row_global - 1, nactive(ispin)) + 1
365 0 : a_nao_col = (j_col_global - 1)/nao + 1
366 0 : b_nao_col = MOD(j_col_global - 1, nao) + 1
367 :
368 : WXmat2%local_data(a_nao_col, i_occ_row) = WXmat2%local_data(a_nao_col, i_occ_row) + &
369 0 : WXaoao%local_data(i_row_global, j_col_global)*evects(ispin, ivect)%local_data(b_nao_col, j_occ_row)
370 : END DO
371 : END DO
372 :
373 0 : CALL cp_fm_release(WXaoao)
374 :
375 0 : IF (iounit > 0) THEN
376 0 : CALL cp_fm_write_formatted(WXmat2, iounit, "WXmat2")
377 : END IF
378 :
379 0 : CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect), -1.0_dp, WXmat2)
380 :
381 0 : CALL cp_fm_release(WXmat2)
382 :
383 : END DO! ispin
384 : END DO !ivect
385 :
386 0 : CALL cp_fm_release(fms)
387 :
388 0 : CALL timestop(handle)
389 :
390 0 : END SUBROUTINE tddfpt_apply_bse_debug
391 :
392 : ! **************************************************************************************************
393 : !> \brief ...
394 : !> \param Aop_evects ...
395 : !> \param evects ...
396 : !> \param gs_mos ...
397 : !> \param qs_env ...
398 : !> \param S_evects ...
399 : ! **************************************************************************************************
400 36 : SUBROUTINE tddfpt_apply_bse(Aop_evects, evects, gs_mos, qs_env, S_evects)
401 :
402 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects
403 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects
404 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
405 : INTENT(in) :: gs_mos
406 : TYPE(qs_environment_type), POINTER :: qs_env
407 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: S_evects
408 :
409 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_bse'
410 :
411 : INTEGER :: handle, ispin, ivect, nao, nspins, &
412 : nvects, nvirt
413 : INTEGER, DIMENSION(2) :: nactive
414 : REAL(KIND=dp) :: alpha
415 36 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries, rho
416 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
417 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
418 : TYPE(cp_fm_type) :: CSvirt, evects_mo, evects_unsplit, fms, &
419 : rhomu, rhosplit
420 : TYPE(cp_fm_type), POINTER :: bse_a_matrix_MO
421 36 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
422 : TYPE(excited_energy_type), POINTER :: ex_env
423 : TYPE(mp_para_env_type), POINTER :: para_env
424 :
425 36 : CALL timeset(routineN, handle)
426 :
427 36 : nspins = SIZE(evects, 1)
428 36 : nvects = SIZE(evects, 2)
429 : IF (nspins > 1) THEN
430 : alpha = 1.0_dp
431 : ELSE
432 : alpha = 2.0_dp
433 : END IF
434 36 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
435 72 : DO ispin = 1, nspins
436 72 : CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
437 : END DO
438 :
439 36 : NULLIFY (ex_env, para_env, blacs_env, matrix_s)
440 : CALL get_qs_env(qs_env, exstate_env=ex_env, para_env=para_env, blacs_env=blacs_env, &
441 36 : matrix_s=matrix_s)
442 :
443 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
444 36 : context=blacs_env, nrow_global=nao, ncol_global=nao)
445 36 : CALL cp_fm_create(fms, fmstruct)
446 36 : CALL cp_fm_struct_release(fmstruct)
447 36 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, fms)
448 :
449 36 : NULLIFY (bse_a_matrix_MO)
450 36 : bse_a_matrix_MO => ex_env%bse_a_matrix_MO(1, 1)
451 36 : para_env => bse_a_matrix_MO%matrix_struct%para_env
452 :
453 380 : DO ivect = 1, nvects
454 724 : DO ispin = 1, nspins
455 344 : NULLIFY (fmstruct)
456 : CALL cp_fm_get_info(matrix=evects(ispin, 1), &
457 344 : nrow_global=nao, ncol_global=nactive(ispin))
458 344 : nvirt = SIZE(gs_mos(ispin)%evals_virt)
459 :
460 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
461 344 : context=blacs_env, nrow_global=nvirt, ncol_global=nao)
462 344 : CALL cp_fm_create(CSvirt, fmstruct)
463 344 : CALL cp_fm_struct_release(fmstruct)
464 : CALL parallel_gemm('T', 'N', nvirt, nao, nao, 1.0_dp, gs_mos(ispin)%mos_virt, &
465 344 : fms, 0.0_dp, CSvirt)
466 :
467 344 : NULLIFY (fmstruct)
468 : CALL cp_fm_struct_create(fmstruct, para_env=para_env, context=blacs_env, &
469 344 : nrow_global=nvirt, ncol_global=nactive(ispin))
470 344 : CALL cp_fm_create(evects_mo, fmstruct)
471 344 : CALL cp_fm_struct_release(fmstruct)
472 344 : NULLIFY (fmstruct)
473 : CALL cp_fm_struct_create(fmstruct, para_env=para_env, context=blacs_env, &
474 344 : nrow_global=nvirt, ncol_global=nactive(ispin))
475 344 : CALL cp_fm_create(rhosplit, fmstruct)
476 344 : CALL cp_fm_struct_release(fmstruct)
477 344 : NULLIFY (fmstruct)
478 : CALL cp_fm_struct_create(fmstruct, para_env=para_env, context=blacs_env, &
479 344 : nrow_global=nao, ncol_global=nactive(ispin))
480 344 : CALL cp_fm_create(rhomu, fmstruct)
481 344 : CALL cp_fm_struct_release(fmstruct)
482 344 : NULLIFY (fmstruct)
483 : CALL cp_fm_struct_create(fmstruct, para_env=para_env, context=blacs_env, &
484 344 : nrow_global=nactive(ispin)*nvirt, ncol_global=1)
485 344 : CALL cp_fm_create(evects_unsplit, fmstruct)
486 344 : CALL cp_fm_struct_release(fmstruct)
487 :
488 : ! get X_jb
489 : CALL parallel_gemm("T", "N", nvirt, nactive(ispin), nao, 1.0_dp, gs_mos(ispin)%mos_virt, &
490 344 : S_evects(ispin, ivect), 0.0_dp, evects_mo)
491 : ! rearrange X_jb
492 344 : CALL contract_bse(evects_mo, evects_unsplit, nactive(ispin), nvirt, eigvec_entries)
493 : ! contract A_iajb X_jb
494 1032 : ALLOCATE (rho(nvirt*nactive(ispin)))
495 344 : CALL cp_fm_matvec(bse_a_matrix_MO, eigvec_entries, rho, 1.0_dp, 0.0_dp)
496 : ! rearrange rho_ia
497 344 : CALL split_bse(rho, rhosplit, nvirt)
498 : ! get rho_imu
499 : CALL parallel_gemm("T", "N", nao, nactive(ispin), nvirt, 1.0_dp, CSvirt, &
500 344 : rhosplit, 0.0_dp, rhomu)
501 :
502 344 : CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect), 1.0_dp, rhomu)
503 :
504 344 : CALL cp_fm_release(rhomu)
505 344 : CALL cp_fm_release(CSvirt)
506 344 : CALL cp_fm_release(evects_mo)
507 344 : CALL cp_fm_release(rhosplit)
508 344 : CALL cp_fm_release(evects_unsplit)
509 344 : DEALLOCATE (rho)
510 2752 : DEALLOCATE (eigvec_entries)
511 :
512 : END DO! ispin
513 : END DO !ivect
514 :
515 36 : CALL cp_fm_release(fms)
516 :
517 36 : CALL timestop(handle)
518 :
519 72 : END SUBROUTINE tddfpt_apply_bse
520 :
521 : ! **************************************************************************************************
522 : !> \brief ...
523 : !> \param evects_mo ...
524 : !> \param evects_unsplit ...
525 : !> \param nactive ...
526 : !> \param nvirt ...
527 : !> \param eigvec_entries ...
528 : ! **************************************************************************************************
529 344 : SUBROUTINE contract_bse(evects_mo, evects_unsplit, nactive, nvirt, eigvec_entries)
530 : TYPE(cp_fm_type), INTENT(IN) :: evects_mo
531 : TYPE(cp_fm_type), INTENT(INOUT) :: evects_unsplit
532 : INTEGER :: nactive, nvirt
533 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigvec_entries
534 :
535 : CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_bse'
536 :
537 : INTEGER :: handle, jj
538 :
539 344 : CALL timeset(routineN, handle)
540 :
541 1032 : ALLOCATE (eigvec_entries(nvirt*nactive))
542 1720 : DO jj = 1, nactive
543 : CALL cp_fm_to_fm_submat(msource=evects_mo, mtarget=evects_unsplit, &
544 : nrow=nvirt, ncol=1, s_firstrow=1, s_firstcol=jj, &
545 1720 : t_firstrow=(jj - 1)*nvirt + 1, t_firstcol=1)
546 : END DO
547 344 : CALL cp_fm_vectorssum(evects_unsplit, eigvec_entries, "R")
548 :
549 344 : CALL timestop(handle)
550 344 : END SUBROUTINE contract_bse
551 :
552 : ! **************************************************************************************************
553 : !> \brief ...
554 : !> \param rho ...
555 : !> \param rhosplit ...
556 : !> \param nvirt ...
557 : ! **************************************************************************************************
558 344 : SUBROUTINE split_bse(rho, rhosplit, nvirt)
559 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
560 : INTENT(IN) :: rho
561 : TYPE(cp_fm_type), INTENT(INOUT) :: rhosplit
562 : INTEGER, INTENT(IN) :: nvirt
563 :
564 : CHARACTER(LEN=*), PARAMETER :: routineN = 'split_bse'
565 :
566 : INTEGER :: handle, i_occ_row, ii, j_occ_row
567 :
568 344 : CALL timeset(routineN, handle)
569 :
570 344 : CALL cp_fm_set_all(rhosplit, 0.0_dp)
571 :
572 26488 : DO ii = 1, SIZE(rho)
573 26144 : i_occ_row = (ii - 1)/nvirt + 1
574 26144 : j_occ_row = MOD(ii - 1, nvirt) + 1
575 26488 : CALL cp_fm_set_element(rhosplit, j_occ_row, i_occ_row, rho(ii))
576 : END DO
577 :
578 344 : CALL timestop(handle)
579 344 : END SUBROUTINE split_bse
580 :
581 : END MODULE qs_tddfpt2_bse_utils
|
