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