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