Line data Source code
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 : 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_scale_and_add
21 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
22 : cp_fm_struct_release,&
23 : cp_fm_struct_type
24 : USE cp_fm_types, ONLY: cp_fm_create,&
25 : cp_fm_get_info,&
26 : cp_fm_release,&
27 : cp_fm_set_all,&
28 : cp_fm_set_element,&
29 : cp_fm_to_fm,&
30 : cp_fm_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
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 0 : 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 0 : 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 0 : POINTER :: sab_orb
90 :
91 0 : CALL timeset(routineN, handle)
92 :
93 0 : NULLIFY (ex_env, sab_orb)
94 0 : CALL get_qs_env(qs_env, exstate_env=ex_env, sab_orb=sab_orb)
95 :
96 0 : nspins = SIZE(matrix_ks, 1)
97 0 : nspins = SIZE(evects, 1)
98 0 : nvects = SIZE(evects, 2)
99 :
100 0 : DO ispin = 1, nspins
101 :
102 0 : CPASSERT(.NOT. ASSOCIATED(gs_mos(ispin)%evals_occ_matrix))
103 :
104 0 : CALL dbcsr_create(matrixf, template=matrix_s)
105 0 : nmo = SIZE(ex_env%gw_eigen)
106 : CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
107 0 : nrow_global=nao, ncol_global=nactive)
108 0 : NULLIFY (blacs_env, para_env)
109 0 : CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
110 :
111 0 : occ = SIZE(gs_mos(ispin)%evals_occ)
112 0 : nactive = gs_mos(ispin)%nmo_active
113 0 : nmo = SIZE(ex_env%gw_eigen)
114 0 : virt = SIZE(gs_mos(ispin)%evals_virt)
115 0 : NULLIFY (fmstruct)
116 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
117 0 : context=blacs_env, nrow_global=virt, ncol_global=virt)
118 0 : CALL cp_fm_create(matrixtmp, fmstruct)
119 0 : CALL cp_fm_struct_release(fmstruct)
120 :
121 0 : NULLIFY (fmstruct)
122 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
123 0 : context=blacs_env, nrow_global=virt, ncol_global=nao)
124 0 : CALL cp_fm_create(matrixtmp2, fmstruct)
125 0 : CALL cp_fm_struct_release(fmstruct)
126 :
127 0 : NULLIFY (fmstruct)
128 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
129 0 : context=blacs_env, nrow_global=nao, ncol_global=nao)
130 0 : CALL cp_fm_create(matrixtmp3, fmstruct)
131 0 : CALL cp_fm_create(fms, fmstruct)
132 0 : CALL cp_fm_struct_release(fmstruct)
133 0 : CALL cp_dbcsr_alloc_block_from_nbl(matrixf, sab_orb)
134 :
135 : !--add virt eigenvalues
136 0 : CALL dbcsr_set(matrixf, 0.0_dp)
137 0 : CALL cp_fm_create(weighted_vect, gs_mos(ispin)%mos_virt%matrix_struct)
138 0 : CALL cp_fm_create(Sweighted_vect, gs_mos(ispin)%mos_virt%matrix_struct)
139 0 : CALL cp_fm_to_fm(gs_mos(ispin)%mos_virt, weighted_vect)
140 0 : CALL copy_dbcsr_to_fm(matrix_s, fms)
141 :
142 0 : ALLOCATE (gw_virt(virt))
143 0 : ALLOCATE (gw_occ(nactive))
144 0 : gw_virt(1:virt) = ex_env%gw_eigen(occ + 1:nmo)
145 0 : DO i = 1, nactive
146 0 : j = gs_mos(ispin)%index_active(i)
147 0 : gw_occ(i) = ex_env%gw_eigen(j)
148 : END DO
149 :
150 0 : CALL cp_fm_set_all(matrixtmp, 0.0_dp)
151 0 : DO i = 1, virt
152 0 : CALL cp_fm_set_element(matrixtmp, i, i, gw_virt(i))
153 : END DO
154 0 : DEALLOCATE (gw_virt)
155 0 : CALL parallel_gemm('N', 'N', nao, virt, nao, 1.0_dp, fms, weighted_vect, 0.0_dp, Sweighted_vect)
156 0 : CALL parallel_gemm('N', 'T', virt, nao, virt, 1.0_dp, matrixtmp, Sweighted_vect, 0.0_dp, matrixtmp2)
157 0 : CALL parallel_gemm('N', 'N', nao, nao, virt, 1.0_dp, Sweighted_vect, matrixtmp2, 0.0_dp, matrixtmp3)
158 0 : CALL copy_fm_to_dbcsr(matrixtmp3, matrixf)
159 :
160 0 : CALL cp_fm_release(weighted_vect)
161 0 : CALL cp_fm_release(Sweighted_vect)
162 0 : CALL cp_fm_release(fmS)
163 : !--add occ eigenvalues
164 : CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
165 0 : nrow_global=nao, ncol_global=nactive)
166 0 : CALL cp_fm_create(hevec, matrix_struct)
167 :
168 0 : DO ivect = 1, nvects
169 : CALL cp_dbcsr_sm_fm_multiply(matrixf, evects(ispin, ivect), &
170 : Aop_evects(ispin, ivect), ncol=nactive, &
171 0 : alpha=1.0_dp, beta=1.0_dp)
172 :
173 0 : CALL cp_fm_to_fm(S_evects(ispin, ivect), hevec)
174 0 : CALL cp_fm_column_scale(hevec, gw_occ)
175 :
176 0 : CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect), -1.0_dp, hevec)
177 : END DO !ivect
178 0 : DEALLOCATE (gw_occ)
179 :
180 0 : CALL cp_fm_release(matrixtmp)
181 0 : CALL cp_fm_release(matrixtmp2)
182 0 : CALL cp_fm_release(matrixtmp3)
183 :
184 0 : CALL dbcsr_release(matrixf)
185 0 : CALL cp_fm_release(hevec)
186 : END DO !ispin
187 :
188 0 : virt = SIZE(Aop_evects, 2)
189 0 : CALL timestop(handle)
190 :
191 0 : END SUBROUTINE zeroth_order_gw
192 :
193 : ! **************************************************************************************************
194 : !> \brief Update action of TDDFPT operator on trial vectors by adding BSE W term.
195 : !> \param Aop_evects ...
196 : !> \param evects ...
197 : !> \param gs_mos ...
198 : !> \param qs_env ...
199 : !> \note Based on the subroutine tddfpt_apply_hfx() which was originally created by
200 : !> Mohamed Fawzi on 10.2002.
201 : ! **************************************************************************************************
202 0 : SUBROUTINE tddfpt_apply_bse(Aop_evects, evects, gs_mos, qs_env)
203 :
204 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects
205 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects
206 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
207 : INTENT(in) :: gs_mos
208 : TYPE(qs_environment_type), POINTER :: qs_env
209 :
210 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_bse'
211 :
212 : INTEGER :: a_nao_col, a_virt_col, b_nao_col, c_virt_col, handle, i_occ_row, i_row_global, &
213 : ii, ispin, ivect, j_col_global, j_occ_row, jj, k_occ_col, mu_col_global, nao, ncol_block, &
214 : ncol_local, nrow_block, nrow_local, nspins, nvects, nvirt
215 : INTEGER, DIMENSION(2) :: nactive
216 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
217 : REAL(KIND=dp) :: alpha
218 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
219 : TYPE(cp_fm_struct_type), POINTER :: fmstruct, matrix_struct
220 : TYPE(cp_fm_type) :: CSvirt, fms, WXaoao, WXmat2, WXvirtao
221 : TYPE(cp_fm_type), POINTER :: bse_w_matrix_MO
222 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
223 : TYPE(excited_energy_type), POINTER :: ex_env
224 : TYPE(mp_para_env_type), POINTER :: para_env
225 :
226 0 : CALL timeset(routineN, handle)
227 :
228 0 : nspins = SIZE(evects, 1)
229 0 : nvects = SIZE(evects, 2)
230 : IF (nspins > 1) THEN
231 : alpha = 1.0_dp
232 : ELSE
233 : alpha = 2.0_dp
234 : END IF
235 0 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
236 0 : DO ispin = 1, nspins
237 0 : CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
238 : END DO
239 :
240 0 : NULLIFY (ex_env, para_env, blacs_env, matrix_s)
241 : CALL get_qs_env(qs_env, exstate_env=ex_env, para_env=para_env, blacs_env=blacs_env, &
242 0 : matrix_s=matrix_s)
243 :
244 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
245 0 : context=blacs_env, nrow_global=nao, ncol_global=nao)
246 0 : CALL cp_fm_create(fms, fmstruct)
247 0 : CALL cp_fm_struct_release(fmstruct)
248 0 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, fms)
249 :
250 0 : NULLIFY (bse_w_matrix_MO)
251 0 : bse_w_matrix_MO => ex_env%bse_w_matrix_MO(1, 1)
252 :
253 0 : DO ivect = 1, nvects
254 0 : DO ispin = 1, nspins
255 0 : NULLIFY (matrix_struct, fmstruct)
256 : CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
257 0 : nrow_global=nao, ncol_global=nactive(ispin))
258 0 : nvirt = SIZE(gs_mos(ispin)%evals_virt)
259 :
260 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
261 0 : context=blacs_env, nrow_global=nvirt, ncol_global=nao)
262 0 : CALL cp_fm_create(CSvirt, fmstruct)
263 0 : CALL cp_fm_struct_release(fmstruct)
264 :
265 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
266 : context=blacs_env, nrow_global=nactive(ispin)*nactive(ispin), &
267 0 : ncol_global=nvirt*nao)
268 0 : CALL cp_fm_create(WXvirtao, fmstruct)
269 0 : CALL cp_fm_struct_release(fmstruct)
270 0 : CALL cp_fm_set_all(WXvirtao, 0.0_dp)
271 :
272 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
273 : context=blacs_env, nrow_global=nactive(ispin)*nactive(ispin), &
274 0 : ncol_global=nao*nao)
275 0 : CALL cp_fm_create(WXaoao, fmstruct)
276 0 : CALL cp_fm_struct_release(fmstruct)
277 0 : CALL cp_fm_set_all(WXaoao, 0.0_dp)
278 :
279 0 : CALL parallel_gemm('T', 'N', nvirt, nao, nao, 1.0_dp, gs_mos(ispin)%mos_virt, fms, 0.0_dp, CSvirt)
280 0 : NULLIFY (row_indices, col_indices)
281 : CALL cp_fm_get_info(matrix=WXvirtao, nrow_local=nrow_local, ncol_local=ncol_local, &
282 : row_indices=row_indices, col_indices=col_indices, &
283 0 : nrow_block=nrow_block, ncol_block=ncol_block)
284 :
285 0 : CALL cp_fm_set_all(WXvirtao, 0.0_dp)
286 0 : DO ii = 1, nrow_local
287 0 : i_row_global = row_indices(ii)
288 0 : DO jj = 1, ncol_local
289 0 : j_col_global = col_indices(jj)
290 :
291 0 : i_occ_row = (i_row_global - 1)/nactive(ispin) + 1
292 0 : j_occ_row = MOD(i_row_global - 1, nactive(ispin)) + 1
293 :
294 0 : a_virt_col = (j_col_global - 1)/nao + 1
295 0 : b_nao_col = MOD(j_col_global - 1, nao) + 1
296 :
297 0 : DO c_virt_col = 1, nvirt
298 0 : mu_col_global = (a_virt_col - 1)*nvirt + c_virt_col
299 :
300 : WXvirtao%local_data(i_row_global, j_col_global) = WXvirtao%local_data(i_row_global, j_col_global) + &
301 0 : bse_w_matrix_MO%local_data(i_row_global, mu_col_global)*CSvirt%local_data(c_virt_col, b_nao_col)
302 :
303 : END DO
304 : END DO
305 : END DO
306 :
307 0 : NULLIFY (row_indices, col_indices) ! redefine indices
308 : CALL cp_fm_get_info(matrix=WXaoao, nrow_local=nrow_local, ncol_local=ncol_local, &
309 : row_indices=row_indices, col_indices=col_indices, &
310 0 : nrow_block=nrow_block, ncol_block=ncol_block)
311 :
312 0 : CALL cp_fm_set_all(WXaoao, 0.0_dp)
313 0 : DO ii = 1, nrow_local
314 0 : i_row_global = row_indices(ii)
315 0 : DO jj = 1, ncol_local
316 0 : j_col_global = col_indices(jj)
317 :
318 0 : i_occ_row = (i_row_global - 1)/nactive(ispin) + 1
319 0 : j_occ_row = MOD(i_row_global - 1, nactive(ispin)) + 1
320 :
321 0 : a_nao_col = (j_col_global - 1)/nao + 1
322 0 : b_nao_col = MOD(j_col_global - 1, nao) + 1
323 :
324 0 : DO k_occ_col = 1, nvirt
325 0 : mu_col_global = (k_occ_col - 1)*nao + a_nao_col
326 :
327 : WXaoao%local_data(i_row_global, j_col_global) = WXaoao%local_data(i_row_global, j_col_global) + &
328 0 : WXvirtao%local_data(i_row_global, mu_col_global)*CSvirt%local_data(k_occ_col, b_nao_col)
329 :
330 : END DO
331 : END DO
332 : END DO
333 :
334 0 : CALL cp_fm_release(WXvirtao)
335 0 : CALL cp_fm_release(CSvirt)
336 :
337 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, &
338 0 : context=blacs_env, nrow_global=nao, ncol_global=nactive(ispin))
339 0 : CALL cp_fm_create(WXmat2, fmstruct)
340 0 : CALL cp_fm_struct_release(fmstruct)
341 0 : CALL cp_fm_set_all(WXmat2, 0.0_dp)
342 :
343 0 : DO ii = 1, nrow_local
344 0 : i_row_global = row_indices(ii)
345 0 : DO jj = 1, ncol_local
346 0 : j_col_global = col_indices(jj)
347 :
348 0 : i_occ_row = (i_row_global - 1)/nactive(ispin) + 1
349 0 : j_occ_row = MOD(i_row_global - 1, nactive(ispin)) + 1
350 0 : a_nao_col = (j_col_global - 1)/nao + 1
351 0 : b_nao_col = MOD(j_col_global - 1, nao) + 1
352 0 : IF (a_nao_col == b_nao_col) THEN
353 : WXmat2%local_data(b_nao_col, i_occ_row) = WXmat2%local_data(b_nao_col, i_occ_row) + &
354 0 : WXaoao%local_data(i_row_global, j_col_global)*evects(ispin, ivect)%local_data(a_nao_col, j_occ_row)
355 : END IF
356 0 : IF (a_nao_col /= b_nao_col) THEN
357 : WXmat2%local_data(a_nao_col, i_occ_row) = WXmat2%local_data(a_nao_col, i_occ_row) + &
358 0 : WXaoao%local_data(i_row_global, j_col_global)*evects(ispin, ivect)%local_data(b_nao_col, j_occ_row)
359 : END IF
360 : END DO
361 : END DO
362 :
363 0 : CALL cp_fm_release(WXaoao)
364 :
365 0 : CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect), -1.0_dp, WXmat2)
366 :
367 0 : CALL cp_fm_release(WXmat2)
368 :
369 : END DO! ispin
370 : END DO !ivect
371 :
372 0 : CALL cp_fm_release(fms)
373 :
374 0 : CALL timestop(handle)
375 :
376 0 : END SUBROUTINE tddfpt_apply_bse
377 : END MODULE qs_tddfpt2_bse_utils
|
