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_fhxc
9 : USE admm_types, ONLY: admm_type
10 : USE cp_control_types, ONLY: dft_control_type,&
11 : stda_control_type
12 : USE cp_dbcsr_api, ONLY: &
13 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_info, &
14 : dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_symmetric
15 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
16 : USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,&
17 : cp_dbcsr_plus_fm_fm_t,&
18 : cp_dbcsr_sm_fm_multiply
19 : USE cp_fm_types, ONLY: cp_fm_create,&
20 : cp_fm_get_info,&
21 : cp_fm_release,&
22 : cp_fm_type
23 : USE input_constants, ONLY: do_admm_aux_exch_func_none,&
24 : no_sf_tddfpt,&
25 : tddfpt_sf_col,&
26 : tddfpt_sf_noncol
27 : USE kinds, ONLY: default_string_length,&
28 : dp
29 : USE lri_environment_types, ONLY: lri_kind_type
30 : USE message_passing, ONLY: mp_para_env_type
31 : USE parallel_gemm_api, ONLY: parallel_gemm
32 : USE pw_env_types, ONLY: pw_env_get
33 : USE pw_methods, ONLY: pw_axpy,&
34 : pw_scale,&
35 : pw_zero
36 : USE pw_pool_types, ONLY: pw_pool_type
37 : USE pw_types, ONLY: pw_c1d_gs_type,&
38 : pw_r3d_rs_type
39 : USE qs_environment_types, ONLY: get_qs_env,&
40 : qs_environment_type
41 : USE qs_fxc, ONLY: qs_fxc_apply
42 : USE qs_gapw_densities, ONLY: prepare_gapw_den
43 : USE qs_integrate_potential, ONLY: integrate_v_rspace,&
44 : integrate_v_rspace_one_center
45 : USE qs_kernel_types, ONLY: full_kernel_env_type
46 : USE qs_ks_atom, ONLY: update_ks_atom
47 : USE qs_rho_atom_types, ONLY: rho_atom_type
48 : USE qs_rho_methods, ONLY: qs_rho_update_rho,&
49 : qs_rho_update_tddfpt
50 : USE qs_rho_types, ONLY: qs_rho_get
51 : USE qs_tddfpt2_densities, ONLY: tddfpt_construct_aux_fit_density
52 : USE qs_tddfpt2_lri_utils, ONLY: tddfpt2_lri_Amat
53 : USE qs_tddfpt2_operators, ONLY: tddfpt_apply_coulomb,&
54 : tddfpt_apply_xc_potential
55 : USE qs_tddfpt2_stda_types, ONLY: stda_env_type
56 : USE qs_tddfpt2_stda_utils, ONLY: stda_calculate_kernel
57 : USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
58 : USE qs_tddfpt2_types, ONLY: tddfpt_work_matrices
59 : USE task_list_types, ONLY: task_list_type
60 : #include "./base/base_uses.f90"
61 :
62 : IMPLICIT NONE
63 :
64 : PRIVATE
65 :
66 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_fhxc'
67 :
68 : INTEGER, PARAMETER, PRIVATE :: maxspins = 2
69 :
70 : PUBLIC :: fhxc_kernel, stda_kernel
71 :
72 : ! **************************************************************************************************
73 :
74 : CONTAINS
75 :
76 : ! **************************************************************************************************
77 : !> \brief Compute action matrix-vector products with the FHxc Kernel
78 : !> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
79 : !> \param evects TDDFPT trial vectors
80 : !> \param is_rks_triplets indicates that a triplet excited states calculation using
81 : !> spin-unpolarised molecular orbitals has been requested
82 : !> \param do_hfx flag that activates computation of exact-exchange terms
83 : !> \param do_admm ...
84 : !> \param qs_env Quickstep environment
85 : !> \param kernel_env kernel environment
86 : !> \param kernel_env_admm_aux kernel environment for ADMM correction
87 : !> \param sub_env parallel (sub)group environment
88 : !> \param work_matrices collection of work matrices (modified on exit)
89 : !> \param admm_symm use symmetric definition of ADMM kernel correction
90 : !> \param admm_xc_correction use ADMM XC kernel correction
91 : !> \param do_lrigpw ...
92 : !> \param tddfpt_mgrid ...
93 : !> \par History
94 : !> * 06.2016 created [Sergey Chulkov]
95 : !> * 03.2017 refactored [Sergey Chulkov]
96 : !> * 04.2019 refactored [JHU]
97 : ! **************************************************************************************************
98 4420 : SUBROUTINE fhxc_kernel(Aop_evects, evects, is_rks_triplets, &
99 : do_hfx, do_admm, qs_env, kernel_env, kernel_env_admm_aux, &
100 : sub_env, work_matrices, admm_symm, admm_xc_correction, do_lrigpw, &
101 : tddfpt_mgrid)
102 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects
103 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects
104 : LOGICAL, INTENT(in) :: is_rks_triplets, do_hfx, do_admm
105 : TYPE(qs_environment_type), POINTER :: qs_env
106 : TYPE(full_kernel_env_type), POINTER :: kernel_env, kernel_env_admm_aux
107 : TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
108 : TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices
109 : LOGICAL, INTENT(in) :: admm_symm, admm_xc_correction, &
110 : do_lrigpw, tddfpt_mgrid
111 :
112 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fhxc_kernel'
113 :
114 : CHARACTER(LEN=default_string_length) :: basis_type
115 : INTEGER :: handle, ikind, ispin, ivect, nao, &
116 : nao_aux, nkind, nspins, nvects, &
117 : spinflip
118 4420 : INTEGER, DIMENSION(:), POINTER :: blk_sizes
119 : INTEGER, DIMENSION(maxspins) :: nactive
120 : LOGICAL :: do_noncol, do_onecenter, do_scale, gapw, &
121 : gapw_xc
122 : TYPE(admm_type), POINTER :: admm_env
123 : TYPE(cp_fm_type) :: work_aux_orb, work_orb_orb
124 4420 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: A_xc_munu_sub, rho_ia_ao, &
125 4420 : rho_ia_ao_aux_fit
126 : TYPE(dbcsr_type), POINTER :: dbwork
127 : TYPE(dft_control_type), POINTER :: dft_control
128 4420 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
129 : TYPE(mp_para_env_type), POINTER :: para_env
130 4420 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_ia_g, rho_ia_g_aux_fit
131 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
132 4420 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: V_rspace_sub
133 4420 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_ia_r, rho_ia_r_aux_fit
134 : TYPE(pw_r3d_rs_type), POINTER :: weights
135 4420 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set, rho_atom_set
136 : TYPE(task_list_type), POINTER :: task_list
137 :
138 4420 : CALL timeset(routineN, handle)
139 :
140 4420 : nspins = SIZE(evects, 1)
141 4420 : nvects = SIZE(evects, 2)
142 4420 : IF (do_admm) THEN
143 880 : CPASSERT(do_hfx)
144 880 : CPASSERT(ASSOCIATED(sub_env%admm_A))
145 : END IF
146 4420 : CALL get_qs_env(qs_env, dft_control=dft_control)
147 :
148 4420 : gapw = dft_control%qs_control%gapw
149 4420 : gapw_xc = dft_control%qs_control%gapw_xc
150 4420 : spinflip = dft_control%tddfpt2_control%spinflip
151 :
152 4420 : do_noncol = spinflip == tddfpt_sf_noncol
153 4420 : do_scale = is_rks_triplets
154 4420 : do_onecenter = gapw .OR. gapw_xc
155 :
156 4420 : CALL cp_fm_get_info(evects(1, 1), nrow_global=nao)
157 9604 : DO ispin = 1, nspins
158 9604 : CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
159 : END DO
160 :
161 : CALL qs_rho_get(work_matrices%rho_orb_struct_sub, rho_ao=rho_ia_ao, &
162 4420 : rho_g=rho_ia_g, rho_r=rho_ia_r)
163 4420 : IF (do_hfx .AND. do_admm) THEN
164 880 : CALL get_qs_env(qs_env, admm_env=admm_env)
165 : CALL qs_rho_get(work_matrices%rho_aux_fit_struct_sub, &
166 : rho_ao=rho_ia_ao_aux_fit, rho_g=rho_ia_g_aux_fit, &
167 880 : rho_r=rho_ia_r_aux_fit)
168 : END IF
169 :
170 4420 : NULLIFY (weights)
171 4420 : CALL get_qs_env(qs_env, xcint_weights=weights)
172 :
173 14108 : DO ivect = 1, nvects
174 :
175 : ! Transform TDDFT vectors to AO space and store them into rho_ia_ao
176 9688 : IF (ALLOCATED(work_matrices%evects_sub)) THEN
177 16 : IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct)) THEN
178 16 : DO ispin = 1, nspins
179 8 : CALL dbcsr_set(rho_ia_ao(ispin)%matrix, 0.0_dp)
180 : CALL cp_dbcsr_plus_fm_fm_t(rho_ia_ao(ispin)%matrix, &
181 : matrix_v=sub_env%mos_active(ispin), &
182 : matrix_g=work_matrices%evects_sub(ispin, ivect), &
183 16 : ncol=nactive(ispin), symmetry_mode=1)
184 : END DO
185 : ELSE
186 : ! skip trial vectors which are assigned to different parallel groups
187 : CYCLE
188 : END IF
189 : ELSE
190 21288 : DO ispin = 1, nspins
191 11616 : CALL dbcsr_set(rho_ia_ao(ispin)%matrix, 0.0_dp)
192 : CALL cp_dbcsr_plus_fm_fm_t(rho_ia_ao(ispin)%matrix, &
193 : matrix_v=sub_env%mos_active(ispin), &
194 : matrix_g=evects(ispin, ivect), &
195 21288 : ncol=nactive(ispin), symmetry_mode=1)
196 : END DO
197 : END IF
198 :
199 9680 : IF (do_lrigpw) THEN
200 : CALL qs_rho_update_tddfpt(work_matrices%rho_orb_struct_sub, qs_env, &
201 : pw_env_external=sub_env%pw_env, &
202 : task_list_external=sub_env%task_list_orb, &
203 : para_env_external=sub_env%para_env, &
204 : tddfpt_lri_env=kernel_env%lri_env, &
205 172 : tddfpt_lri_density=kernel_env%lri_density)
206 9508 : ELSE IF (dft_control%qs_control%lrigpw .OR. &
207 : dft_control%qs_control%rigpw) THEN
208 : CALL qs_rho_update_tddfpt(work_matrices%rho_orb_struct_sub, qs_env, &
209 : pw_env_external=sub_env%pw_env, &
210 : task_list_external=sub_env%task_list_orb, &
211 0 : para_env_external=sub_env%para_env)
212 : ELSE
213 9508 : IF (gapw) THEN
214 : CALL qs_rho_update_rho(work_matrices%rho_orb_struct_sub, qs_env, &
215 : local_rho_set=work_matrices%local_rho_set, &
216 : pw_env_external=sub_env%pw_env, &
217 : task_list_external=sub_env%task_list_orb_soft, &
218 2232 : para_env_external=sub_env%para_env)
219 : CALL prepare_gapw_den(qs_env, work_matrices%local_rho_set, &
220 2232 : do_rho0=(.NOT. is_rks_triplets), pw_env_sub=sub_env%pw_env)
221 7276 : ELSE IF (gapw_xc) THEN
222 : CALL qs_rho_update_rho(work_matrices%rho_orb_struct_sub, qs_env, &
223 : rho_xc_external=work_matrices%rho_xc_struct_sub, &
224 : local_rho_set=work_matrices%local_rho_set, &
225 : pw_env_external=sub_env%pw_env, &
226 : task_list_external=sub_env%task_list_orb, &
227 : task_list_external_soft=sub_env%task_list_orb_soft, &
228 442 : para_env_external=sub_env%para_env)
229 : CALL prepare_gapw_den(qs_env, work_matrices%local_rho_set, do_rho0=.FALSE., &
230 442 : pw_env_sub=sub_env%pw_env)
231 : ELSE
232 : CALL qs_rho_update_rho(work_matrices%rho_orb_struct_sub, qs_env, &
233 : pw_env_external=sub_env%pw_env, &
234 : task_list_external=sub_env%task_list_orb, &
235 6834 : para_env_external=sub_env%para_env)
236 : END IF
237 : END IF
238 :
239 21304 : DO ispin = 1, nspins
240 21304 : CALL dbcsr_set(work_matrices%A_ia_munu_sub(ispin)%matrix, 0.0_dp)
241 : END DO
242 :
243 : ! electron-hole exchange-correlation interaction
244 21304 : DO ispin = 1, nspins
245 21304 : CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin))
246 : END DO
247 :
248 : ! Skip kernel if collinear xc-kernel for spin-flip is requested
249 9680 : IF (spinflip /= tddfpt_sf_col) THEN
250 9562 : IF ((.NOT. dft_control%tddfpt2_control%do_bse) .AND. (.NOT. dft_control%tddfpt2_control%do_bse_w_only)) THEN
251 9218 : IF ((.NOT. dft_control%tddfpt2_control%do_bse_gw_only)) THEN
252 : ! C_x d^{2}E_{x}^{DFT}[\rho] / d\rho^2
253 : ! + C_{HF} d^{2}E_{x, ADMM}^{DFT}[\rho] / d\rho^2 in case of ADMM calculation
254 20380 : DO ispin = 1, nspins
255 11162 : CALL pw_zero(work_matrices%wpw_rspace_sub(ispin))
256 20380 : IF (ASSOCIATED(work_matrices%wpw_tau_rspace_sub)) CALL pw_zero(work_matrices%wpw_tau_rspace_sub(ispin))
257 : END DO
258 9218 : IF (do_onecenter) THEN
259 2674 : rho_atom_set => sub_env%local_rho_set%rho_atom_set
260 2674 : rho1_atom_set => work_matrices%local_rho_set%rho_atom_set
261 : END IF
262 : !
263 9218 : IF (gapw_xc) THEN
264 442 : IF (kernel_env%do_exck) THEN
265 0 : CPABORT("NYA")
266 : ELSE
267 : CALL qs_fxc_apply(qs_env, kernel_env%xc_deriv_set, kernel_env%xc_rho_set, &
268 : work_matrices%rho_xc_struct_sub, rho_atom_set, &
269 : kernel_env%xc_section, do_onecenter, work_matrices%wpw_rspace_sub, &
270 : work_matrices%wpw_tau_rspace_sub, rho1_atom_set, &
271 : do_scale=do_scale, is_triplet=is_rks_triplets, spinflip=do_noncol, &
272 442 : pw_env_ext=sub_env%pw_env, para_env_external=sub_env%para_env)
273 884 : DO ispin = 1, nspins
274 : CALL pw_axpy(work_matrices%wpw_rspace_sub(ispin), &
275 884 : work_matrices%A_ia_rspace_sub(ispin), kernel_env%alpha)
276 : END DO
277 :
278 : END IF
279 884 : DO ispin = 1, nspins
280 : CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), &
281 442 : work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol)
282 : CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
283 : hmat=work_matrices%A_ia_munu_sub(ispin), &
284 : qs_env=qs_env, calculate_forces=.FALSE., gapw=gapw_xc, &
285 : pw_env_external=sub_env%pw_env, &
286 442 : task_list_external=sub_env%task_list_orb_soft)
287 884 : CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin))
288 : END DO
289 : ELSE
290 8776 : IF (kernel_env%do_exck) THEN
291 : CALL tddfpt_apply_xc_potential(work_matrices%A_ia_rspace_sub, work_matrices%fxc_rspace_sub, &
292 156 : work_matrices%rho_orb_struct_sub, is_rks_triplets)
293 : ELSE
294 : CALL qs_fxc_apply(qs_env, kernel_env%xc_deriv_set, kernel_env%xc_rho_set, &
295 : work_matrices%rho_orb_struct_sub, rho_atom_set, &
296 : kernel_env%xc_section, do_onecenter, work_matrices%wpw_rspace_sub, &
297 : work_matrices%wpw_tau_rspace_sub, rho1_atom_set, &
298 : do_scale=do_scale, is_triplet=is_rks_triplets, spinflip=do_noncol, &
299 8620 : pw_env_ext=sub_env%pw_env, para_env_external=sub_env%para_env)
300 19184 : DO ispin = 1, nspins
301 : CALL pw_axpy(work_matrices%wpw_rspace_sub(ispin), &
302 19184 : work_matrices%A_ia_rspace_sub(ispin), kernel_env%alpha)
303 : END DO
304 :
305 : END IF
306 : END IF
307 :
308 : END IF ! do_bse
309 : END IF ! do_bse
310 : END IF ! spin-flip
311 :
312 : ! ADMM correction
313 : IF ((.NOT. dft_control%tddfpt2_control%do_bse) .AND. (.NOT. dft_control%tddfpt2_control%do_bse_w_only) &
314 9680 : .AND. (.NOT. dft_control%tddfpt2_control%do_bse_gw_only)) THEN
315 9336 : IF (do_admm .AND. admm_xc_correction) THEN
316 1370 : IF (dft_control%admm_control%aux_exch_func /= do_admm_aux_exch_func_none) THEN
317 : CALL tddfpt_construct_aux_fit_density(rho_orb_struct=work_matrices%rho_orb_struct_sub, &
318 : rho_aux_fit_struct=work_matrices%rho_aux_fit_struct_sub, &
319 : local_rho_set=work_matrices%local_rho_set_admm, &
320 : qs_env=qs_env, sub_env=sub_env, &
321 : wfm_rho_orb=work_matrices%rho_ao_orb_fm_sub, &
322 : wfm_rho_aux_fit=work_matrices%rho_ao_aux_fit_fm_sub, &
323 900 : wfm_aux_orb=work_matrices%wfm_aux_orb_sub)
324 : ! - C_{HF} d^{2}E_{x, ADMM}^{DFT}[\hat{\rho}] / d\hat{\rho}^2
325 900 : IF (admm_symm) THEN
326 900 : CALL dbcsr_get_info(rho_ia_ao_aux_fit(1)%matrix, row_blk_size=blk_sizes)
327 3600 : ALLOCATE (A_xc_munu_sub(nspins))
328 1800 : DO ispin = 1, nspins
329 900 : ALLOCATE (A_xc_munu_sub(ispin)%matrix)
330 : CALL dbcsr_create(matrix=A_xc_munu_sub(ispin)%matrix, name="ADMM_XC", &
331 : dist=sub_env%dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
332 900 : row_blk_size=blk_sizes, col_blk_size=blk_sizes)
333 900 : CALL cp_dbcsr_alloc_block_from_nbl(A_xc_munu_sub(ispin)%matrix, sub_env%sab_aux_fit)
334 1800 : CALL dbcsr_set(A_xc_munu_sub(ispin)%matrix, 0.0_dp)
335 : END DO
336 :
337 900 : CALL pw_env_get(sub_env%pw_env, auxbas_pw_pool=auxbas_pw_pool)
338 4500 : ALLOCATE (V_rspace_sub(nspins))
339 1800 : DO ispin = 1, nspins
340 900 : CALL auxbas_pw_pool%create_pw(V_rspace_sub(ispin))
341 1800 : CALL pw_zero(V_rspace_sub(ispin))
342 : END DO
343 :
344 900 : IF (admm_env%do_gapw) THEN
345 190 : basis_type = "AUX_FIT_SOFT"
346 190 : task_list => sub_env%task_list_aux_fit_soft
347 : ELSE
348 710 : basis_type = "AUX_FIT"
349 710 : task_list => sub_env%task_list_aux_fit
350 : END IF
351 :
352 900 : IF (admm_env%do_gapw) THEN
353 190 : rho_atom_set => sub_env%local_rho_set_admm%rho_atom_set
354 190 : rho1_atom_set => work_matrices%local_rho_set_admm%rho_atom_set
355 : END IF
356 1800 : DO ispin = 1, nspins
357 900 : CALL pw_zero(work_matrices%wpw_rspace_sub(ispin))
358 1800 : IF (ASSOCIATED(work_matrices%wpw_tau_rspace_sub)) THEN
359 900 : CALL pw_zero(work_matrices%wpw_tau_rspace_sub(ispin))
360 : END IF
361 : END DO
362 :
363 : CALL qs_fxc_apply(qs_env, kernel_env_admm_aux%xc_deriv_set, kernel_env_admm_aux%xc_rho_set, &
364 : work_matrices%rho_aux_fit_struct_sub, rho_atom_set, &
365 : kernel_env_admm_aux%xc_section, admm_env%do_gapw, &
366 : work_matrices%wpw_rspace_sub, work_matrices%wpw_tau_rspace_sub, rho1_atom_set, &
367 : do_scale=do_scale, is_triplet=is_rks_triplets, spinflip=do_noncol, &
368 : pw_env_ext=sub_env%pw_env, kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
369 900 : para_env_external=sub_env%para_env)
370 1800 : DO ispin = 1, nspins
371 : CALL pw_axpy(work_matrices%wpw_rspace_sub(ispin), V_rspace_sub(ispin), &
372 1800 : kernel_env_admm_aux%alpha)
373 : END DO
374 :
375 1800 : DO ispin = 1, nspins
376 900 : CALL pw_scale(V_rspace_sub(ispin), V_rspace_sub(ispin)%pw_grid%dvol)
377 : CALL integrate_v_rspace(v_rspace=V_rspace_sub(ispin), &
378 : hmat=A_xc_munu_sub(ispin), &
379 : qs_env=qs_env, calculate_forces=.FALSE., &
380 : pw_env_external=sub_env%pw_env, &
381 : basis_type=basis_type, &
382 1800 : task_list_external=task_list)
383 : END DO
384 900 : IF (admm_env%do_gapw) THEN
385 : CALL update_ks_atom(qs_env, A_xc_munu_sub, rho_ia_ao_aux_fit, forces=.FALSE., tddft=.TRUE., &
386 : rho_atom_external=rho1_atom_set, &
387 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
388 : oce_external=admm_env%admm_gapw_env%oce, &
389 190 : sab_external=sub_env%sab_aux_fit)
390 : END IF
391 900 : ALLOCATE (dbwork)
392 900 : CALL dbcsr_create(dbwork, template=work_matrices%A_ia_munu_sub(1)%matrix)
393 : CALL cp_fm_create(work_aux_orb, &
394 900 : matrix_struct=work_matrices%wfm_aux_orb_sub%matrix_struct)
395 : CALL cp_fm_create(work_orb_orb, &
396 900 : matrix_struct=work_matrices%rho_ao_orb_fm_sub%matrix_struct)
397 900 : CALL cp_fm_get_info(work_aux_orb, nrow_global=nao_aux, ncol_global=nao)
398 1800 : DO ispin = 1, nspins
399 : CALL cp_dbcsr_sm_fm_multiply(A_xc_munu_sub(ispin)%matrix, sub_env%admm_A, &
400 900 : work_aux_orb, nao)
401 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, sub_env%admm_A, &
402 900 : work_aux_orb, 0.0_dp, work_orb_orb)
403 900 : CALL dbcsr_copy(dbwork, work_matrices%A_ia_munu_sub(1)%matrix)
404 900 : CALL dbcsr_set(dbwork, 0.0_dp)
405 900 : CALL copy_fm_to_dbcsr(work_orb_orb, dbwork, keep_sparsity=.TRUE.)
406 1800 : CALL dbcsr_add(work_matrices%A_ia_munu_sub(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
407 : END DO
408 900 : CALL dbcsr_release(dbwork)
409 900 : DEALLOCATE (dbwork)
410 1800 : DO ispin = 1, nspins
411 1800 : CALL auxbas_pw_pool%give_back_pw(V_rspace_sub(ispin))
412 : END DO
413 900 : DEALLOCATE (V_rspace_sub)
414 900 : CALL cp_fm_release(work_aux_orb)
415 900 : CALL cp_fm_release(work_orb_orb)
416 1800 : DO ispin = 1, nspins
417 1800 : CALL dbcsr_deallocate_matrix(A_xc_munu_sub(ispin)%matrix)
418 : END DO
419 1800 : DEALLOCATE (A_xc_munu_sub)
420 : ELSE
421 0 : IF (admm_env%do_gapw) THEN
422 0 : CPWARN("GAPW/ADMM needs symmetric ADMM kernel")
423 0 : CPABORT("GAPW/ADMM@TDDFT")
424 : END IF
425 0 : DO ispin = 1, nspins
426 0 : CALL pw_zero(work_matrices%wpw_rspace_sub(ispin))
427 0 : IF (ASSOCIATED(work_matrices%wpw_tau_rspace_sub)) THEN
428 0 : CALL pw_zero(work_matrices%wpw_tau_rspace_sub(ispin))
429 : END IF
430 : END DO
431 :
432 : CALL qs_fxc_apply(qs_env, kernel_env_admm_aux%xc_deriv_set, kernel_env_admm_aux%xc_rho_set, &
433 : work_matrices%rho_aux_fit_struct_sub, rho_atom_set, &
434 : kernel_env_admm_aux%xc_section, admm_env%do_gapw, &
435 : work_matrices%wpw_rspace_sub, work_matrices%wpw_tau_rspace_sub, rho1_atom_set, &
436 : do_scale=do_scale, is_triplet=is_rks_triplets, spinflip=do_noncol, &
437 0 : pw_env_ext=sub_env%pw_env)
438 :
439 0 : DO ispin = 1, nspins
440 : CALL pw_axpy(work_matrices%wpw_rspace_sub(ispin), &
441 0 : work_matrices%A_ia_rspace_sub(ispin), kernel_env_admm_aux%alpha)
442 : END DO
443 :
444 : END IF
445 : END IF
446 : END IF
447 : END IF
448 :
449 : ! electron-hole Coulomb interaction
450 14100 : IF (.NOT. dft_control%tddfpt2_control%do_bse_w_only) THEN
451 9528 : IF (.NOT. dft_control%tddfpt2_control%do_bse_gw_only) THEN
452 9528 : IF ((.NOT. is_rks_triplets) .AND. (spinflip == no_sf_tddfpt)) THEN
453 : ! a sum J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu can be computed by solving
454 : ! the Poisson equation for combined density (rho_{ia,alpha} + rho_{ia,beta}) .
455 : ! The following action will destroy reciprocal-space grid in spin-unrestricted case.
456 10222 : DO ispin = 2, nspins
457 10222 : CALL pw_axpy(rho_ia_g(ispin), rho_ia_g(1))
458 : END DO
459 : CALL tddfpt_apply_coulomb(A_ia_rspace=work_matrices%A_ia_rspace_sub, &
460 : rho_ia_g=rho_ia_g(1), &
461 : local_rho_set=work_matrices%local_rho_set, &
462 : hartree_local=work_matrices%hartree_local, &
463 : qs_env=qs_env, sub_env=sub_env, gapw=gapw, &
464 : work_v_gspace=work_matrices%wpw_gspace_sub(1), &
465 : work_v_rspace=work_matrices%wpw_rspace_sub(1), &
466 8278 : tddfpt_mgrid=tddfpt_mgrid)
467 : END IF
468 :
469 : ! convert from the plane-wave representation into the Gaussian basis set representation
470 21000 : DO ispin = 1, nspins
471 21000 : IF (.NOT. do_lrigpw) THEN
472 : CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), &
473 11300 : work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol)
474 :
475 11300 : IF (gapw) THEN
476 : CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
477 : hmat=work_matrices%A_ia_munu_sub(ispin), &
478 : qs_env=qs_env, calculate_forces=.FALSE., gapw=gapw, &
479 : pw_env_external=sub_env%pw_env, &
480 2352 : task_list_external=sub_env%task_list_orb_soft)
481 8948 : ELSE IF (gapw_xc) THEN
482 442 : IF (.NOT. is_rks_triplets) THEN
483 : CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
484 : hmat=work_matrices%A_ia_munu_sub(ispin), &
485 : qs_env=qs_env, calculate_forces=.FALSE., gapw=.FALSE., &
486 442 : pw_env_external=sub_env%pw_env, task_list_external=sub_env%task_list_orb)
487 : END IF
488 : ELSE
489 : CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
490 : hmat=work_matrices%A_ia_munu_sub(ispin), &
491 : qs_env=qs_env, calculate_forces=.FALSE., gapw=.FALSE., &
492 8506 : pw_env_external=sub_env%pw_env, task_list_external=sub_env%task_list_orb)
493 : END IF
494 : ELSE ! for full kernel using lri
495 : CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), &
496 172 : work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol)
497 172 : lri_v_int => kernel_env%lri_density%lri_coefs(ispin)%lri_kinds
498 172 : CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
499 516 : DO ikind = 1, nkind
500 102304 : lri_v_int(ikind)%v_int = 0.0_dp
501 : END DO
502 : CALL integrate_v_rspace_one_center(work_matrices%A_ia_rspace_sub(ispin), &
503 172 : qs_env, lri_v_int, .FALSE., "P_LRI_AUX")
504 516 : DO ikind = 1, nkind
505 204092 : CALL para_env%sum(lri_v_int(ikind)%v_int)
506 : END DO
507 : END IF ! for full kernel using lri
508 : END DO
509 :
510 : ! local atom contributions
511 9528 : IF (.NOT. do_lrigpw) THEN
512 9356 : IF (gapw .OR. gapw_xc) THEN
513 : ! rho_ia_ao will not be touched
514 : CALL update_ks_atom(qs_env, work_matrices%A_ia_munu_sub, rho_ia_ao, forces=.FALSE., &
515 : rho_atom_external=work_matrices%local_rho_set%rho_atom_set, &
516 2674 : tddft=.TRUE.)
517 : END IF
518 : END IF
519 :
520 : ! calculate Coulomb contribution to response vector for lrigpw !
521 : ! this is restricting lri to Coulomb only at the moment !
522 9528 : IF (do_lrigpw .AND. (.NOT. is_rks_triplets)) THEN !
523 172 : CALL tddfpt2_lri_Amat(qs_env, sub_env, kernel_env%lri_env, lri_v_int, work_matrices%A_ia_munu_sub)
524 : END IF
525 :
526 21000 : DO ispin = 1, nspins
527 21000 : IF (ALLOCATED(work_matrices%evects_sub)) THEN
528 : CALL cp_dbcsr_sm_fm_multiply(work_matrices%A_ia_munu_sub(ispin)%matrix, &
529 : sub_env%mos_active(ispin), &
530 : work_matrices%Aop_evects_sub(ispin, ivect), &
531 8 : ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
532 : ELSE
533 : CALL cp_dbcsr_sm_fm_multiply(work_matrices%A_ia_munu_sub(ispin)%matrix, &
534 : sub_env%mos_active(ispin), &
535 : Aop_evects(ispin, ivect), &
536 11464 : ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
537 : END IF
538 : END DO
539 : END IF
540 : END IF
541 : END DO
542 :
543 4420 : CALL timestop(handle)
544 :
545 8840 : END SUBROUTINE fhxc_kernel
546 :
547 : ! **************************************************************************************************
548 : !> \brief Compute action matrix-vector products with the sTDA Kernel
549 : !> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
550 : !> \param evects TDDFPT trial vectors
551 : !> \param is_rks_triplets indicates that a triplet excited states calculation using
552 : !> spin-unpolarised molecular orbitals has been requested
553 : !> \param qs_env Quickstep environment
554 : !> \param stda_control control parameters for sTDA kernel
555 : !> \param stda_env ...
556 : !> \param sub_env parallel (sub)group environment
557 : !> \param work_matrices collection of work matrices (modified on exit)
558 : !> \par History
559 : !> * 04.2019 initial version [JHU]
560 : ! **************************************************************************************************
561 2512 : SUBROUTINE stda_kernel(Aop_evects, evects, is_rks_triplets, &
562 : qs_env, stda_control, stda_env, &
563 : sub_env, work_matrices)
564 :
565 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects
566 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects
567 : LOGICAL, INTENT(in) :: is_rks_triplets
568 : TYPE(qs_environment_type), POINTER :: qs_env
569 : TYPE(stda_control_type) :: stda_control
570 : TYPE(stda_env_type) :: stda_env
571 : TYPE(tddfpt_subgroup_env_type) :: sub_env
572 : TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices
573 :
574 : CHARACTER(LEN=*), PARAMETER :: routineN = 'stda_kernel'
575 :
576 : INTEGER :: handle, ivect, nvects
577 :
578 2512 : CALL timeset(routineN, handle)
579 :
580 2512 : nvects = SIZE(evects, 2)
581 :
582 9980 : DO ivect = 1, nvects
583 9980 : IF (ALLOCATED(work_matrices%evects_sub)) THEN
584 0 : IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct)) THEN
585 : CALL stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, work_matrices, &
586 : is_rks_triplets, work_matrices%evects_sub(:, ivect), &
587 0 : work_matrices%Aop_evects_sub(:, ivect))
588 : ELSE
589 : ! skip trial vectors which are assigned to different parallel groups
590 : CYCLE
591 : END IF
592 : ELSE
593 : CALL stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, work_matrices, &
594 7468 : is_rks_triplets, evects(:, ivect), Aop_evects(:, ivect))
595 : END IF
596 : END DO
597 :
598 2512 : CALL timestop(handle)
599 :
600 2512 : END SUBROUTINE stda_kernel
601 :
602 : ! **************************************************************************************************
603 :
604 : END MODULE qs_tddfpt2_fhxc
|
