(git:e7e05ae)
qs_tddfpt2_fhxc.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
13  USE cp_dbcsr_operations, ONLY: copy_fm_to_dbcsr,&
14  cp_dbcsr_plus_fm_fm_t,&
15  cp_dbcsr_sm_fm_multiply
16  USE cp_fm_types, ONLY: cp_fm_create,&
17  cp_fm_get_info,&
18  cp_fm_release,&
19  cp_fm_type
20  USE dbcsr_api, ONLY: &
21  dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_info, &
22  dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_symmetric
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 !> \par History
91 !> * 06.2016 created [Sergey Chulkov]
92 !> * 03.2017 refactored [Sergey Chulkov]
93 !> * 04.2019 refactored [JHU]
94 ! **************************************************************************************************
95  SUBROUTINE fhxc_kernel(Aop_evects, evects, is_rks_triplets, &
96  do_hfx, do_admm, qs_env, kernel_env, kernel_env_admm_aux, &
97  sub_env, work_matrices, admm_symm, admm_xc_correction, do_lrigpw)
98  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: aop_evects, evects
99  LOGICAL, INTENT(in) :: is_rks_triplets, do_hfx, do_admm
100  TYPE(qs_environment_type), POINTER :: qs_env
101  TYPE(full_kernel_env_type), POINTER :: kernel_env, kernel_env_admm_aux
102  TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
103  TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices
104  LOGICAL, INTENT(in) :: admm_symm, admm_xc_correction, do_lrigpw
105 
106  CHARACTER(LEN=*), PARAMETER :: routinen = 'fhxc_kernel'
107 
108  CHARACTER(LEN=default_string_length) :: basis_type
109  INTEGER :: handle, ikind, ispin, ivect, nao, &
110  nao_aux, nkind, nspins, nvects
111  INTEGER, DIMENSION(:), POINTER :: blk_sizes
112  INTEGER, DIMENSION(maxspins) :: nactive
113  LOGICAL :: gapw, gapw_xc
114  TYPE(admm_type), POINTER :: admm_env
115  TYPE(cp_fm_type) :: work_aux_orb, work_orb_orb
116  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: a_xc_munu_sub, rho_ia_ao, &
117  rho_ia_ao_aux_fit
118  TYPE(dbcsr_type), POINTER :: dbwork
119  TYPE(dft_control_type), POINTER :: dft_control
120  TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
121  TYPE(mp_para_env_type), POINTER :: para_env
122  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_ia_g, rho_ia_g_aux_fit
123  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
124  TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: v_rspace_sub
125  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_ia_r, rho_ia_r_aux_fit
126  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set, rho_atom_set
127  TYPE(task_list_type), POINTER :: task_list
128 
129  CALL timeset(routinen, handle)
130 
131  nspins = SIZE(evects, 1)
132  nvects = SIZE(evects, 2)
133  IF (do_admm) THEN
134  cpassert(do_hfx)
135  cpassert(ASSOCIATED(sub_env%admm_A))
136  END IF
137  CALL get_qs_env(qs_env, dft_control=dft_control)
138 
139  gapw = dft_control%qs_control%gapw
140  gapw_xc = dft_control%qs_control%gapw_xc
141 
142  CALL cp_fm_get_info(evects(1, 1), nrow_global=nao)
143  DO ispin = 1, nspins
144  CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
145  END DO
146 
147  CALL qs_rho_get(work_matrices%rho_orb_struct_sub, rho_ao=rho_ia_ao, &
148  rho_g=rho_ia_g, rho_r=rho_ia_r)
149  IF (do_hfx .AND. do_admm) THEN
150  CALL get_qs_env(qs_env, admm_env=admm_env)
151  CALL qs_rho_get(work_matrices%rho_aux_fit_struct_sub, &
152  rho_ao=rho_ia_ao_aux_fit, rho_g=rho_ia_g_aux_fit, &
153  rho_r=rho_ia_r_aux_fit)
154  END IF
155 
156  DO ivect = 1, nvects
157  IF (ALLOCATED(work_matrices%evects_sub)) THEN
158  IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct)) THEN
159  DO ispin = 1, nspins
160  CALL dbcsr_set(rho_ia_ao(ispin)%matrix, 0.0_dp)
161  CALL cp_dbcsr_plus_fm_fm_t(rho_ia_ao(ispin)%matrix, &
162  matrix_v=sub_env%mos_occ(ispin), &
163  matrix_g=work_matrices%evects_sub(ispin, ivect), &
164  ncol=nactive(ispin), symmetry_mode=1)
165  END DO
166  ELSE
167  ! skip trial vectors which are assigned to different parallel groups
168  cycle
169  END IF
170  ELSE
171  DO ispin = 1, nspins
172  CALL dbcsr_set(rho_ia_ao(ispin)%matrix, 0.0_dp)
173  CALL cp_dbcsr_plus_fm_fm_t(rho_ia_ao(ispin)%matrix, &
174  matrix_v=sub_env%mos_occ(ispin), &
175  matrix_g=evects(ispin, ivect), &
176  ncol=nactive(ispin), symmetry_mode=1)
177  END DO
178  END IF
179 
180  IF (do_lrigpw) THEN
181  CALL qs_rho_update_tddfpt(work_matrices%rho_orb_struct_sub, qs_env, &
182  pw_env_external=sub_env%pw_env, &
183  task_list_external=sub_env%task_list_orb, &
184  para_env_external=sub_env%para_env, &
185  tddfpt_lri_env=kernel_env%lri_env, &
186  tddfpt_lri_density=kernel_env%lri_density)
187  ELSEIF (dft_control%qs_control%lrigpw .OR. &
188  dft_control%qs_control%rigpw) THEN
189  CALL qs_rho_update_tddfpt(work_matrices%rho_orb_struct_sub, qs_env, &
190  pw_env_external=sub_env%pw_env, &
191  task_list_external=sub_env%task_list_orb, &
192  para_env_external=sub_env%para_env)
193  ELSE
194  IF (gapw) THEN
195  CALL qs_rho_update_rho(work_matrices%rho_orb_struct_sub, qs_env, &
196  local_rho_set=work_matrices%local_rho_set, &
197  pw_env_external=sub_env%pw_env, &
198  task_list_external=sub_env%task_list_orb_soft, &
199  para_env_external=sub_env%para_env)
200  CALL prepare_gapw_den(qs_env, work_matrices%local_rho_set, &
201  do_rho0=(.NOT. is_rks_triplets))
202  ELSEIF (gapw_xc) THEN
203  CALL qs_rho_update_rho(work_matrices%rho_orb_struct_sub, qs_env, &
204  rho_xc_external=work_matrices%rho_xc_struct_sub, &
205  local_rho_set=work_matrices%local_rho_set, &
206  pw_env_external=sub_env%pw_env, &
207  task_list_external=sub_env%task_list_orb, &
208  task_list_external_soft=sub_env%task_list_orb_soft, &
209  para_env_external=sub_env%para_env)
210  CALL prepare_gapw_den(qs_env, work_matrices%local_rho_set, do_rho0=.false.)
211  ELSE
212  CALL qs_rho_update_rho(work_matrices%rho_orb_struct_sub, qs_env, &
213  pw_env_external=sub_env%pw_env, &
214  task_list_external=sub_env%task_list_orb, &
215  para_env_external=sub_env%para_env)
216  END IF
217  END IF
218 
219  DO ispin = 1, nspins
220  CALL dbcsr_set(work_matrices%A_ia_munu_sub(ispin)%matrix, 0.0_dp)
221  END DO
222 
223  ! electron-hole exchange-correlation interaction
224  DO ispin = 1, nspins
225  CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin))
226  END DO
227 
228  ! C_x d^{2}E_{x}^{DFT}[\rho] / d\rho^2
229  ! + C_{HF} d^{2}E_{x, ADMM}^{DFT}[\rho] / d\rho^2 in case of ADMM calculation
230  IF (gapw_xc) THEN
231  IF (kernel_env%do_exck) THEN
232  cpabort("NYA")
233  ELSE
234  CALL tddfpt_apply_xc(a_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, &
235  rho_ia_struct=work_matrices%rho_xc_struct_sub, &
236  is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
237  work_v_xc=work_matrices%wpw_rspace_sub, &
238  work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
239  END IF
240  DO ispin = 1, nspins
241  CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), &
242  work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol)
243  CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
244  hmat=work_matrices%A_ia_munu_sub(ispin), &
245  qs_env=qs_env, calculate_forces=.false., gapw=gapw_xc, &
246  pw_env_external=sub_env%pw_env, &
247  task_list_external=sub_env%task_list_orb_soft)
248  CALL pw_zero(work_matrices%A_ia_rspace_sub(ispin))
249  END DO
250  ELSE
251  IF (kernel_env%do_exck) THEN
252  CALL tddfpt_apply_xc_potential(work_matrices%A_ia_rspace_sub, work_matrices%fxc_rspace_sub, &
253  work_matrices%rho_orb_struct_sub, is_rks_triplets)
254  ELSE
255  CALL tddfpt_apply_xc(a_ia_rspace=work_matrices%A_ia_rspace_sub, kernel_env=kernel_env, &
256  rho_ia_struct=work_matrices%rho_orb_struct_sub, &
257  is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
258  work_v_xc=work_matrices%wpw_rspace_sub, &
259  work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
260  END IF
261  END IF
262  IF (gapw .OR. gapw_xc) THEN
263  rho_atom_set => sub_env%local_rho_set%rho_atom_set
264  rho1_atom_set => work_matrices%local_rho_set%rho_atom_set
265  CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, kernel_env%xc_section, &
266  sub_env%para_env, do_tddfpt2=.true., do_triplet=is_rks_triplets)
267  END IF
268 
269  ! ADMM correction
270  IF (do_admm .AND. admm_xc_correction) THEN
271  IF (dft_control%admm_control%aux_exch_func /= do_admm_aux_exch_func_none) THEN
272  CALL tddfpt_construct_aux_fit_density(rho_orb_struct=work_matrices%rho_orb_struct_sub, &
273  rho_aux_fit_struct=work_matrices%rho_aux_fit_struct_sub, &
274  local_rho_set=work_matrices%local_rho_set_admm, &
275  qs_env=qs_env, sub_env=sub_env, &
276  wfm_rho_orb=work_matrices%rho_ao_orb_fm_sub, &
277  wfm_rho_aux_fit=work_matrices%rho_ao_aux_fit_fm_sub, &
278  wfm_aux_orb=work_matrices%wfm_aux_orb_sub)
279  ! - C_{HF} d^{2}E_{x, ADMM}^{DFT}[\hat{\rho}] / d\hat{\rho}^2
280  IF (admm_symm) THEN
281  CALL dbcsr_get_info(rho_ia_ao_aux_fit(1)%matrix, row_blk_size=blk_sizes)
282  ALLOCATE (a_xc_munu_sub(nspins))
283  DO ispin = 1, nspins
284  ALLOCATE (a_xc_munu_sub(ispin)%matrix)
285  CALL dbcsr_create(matrix=a_xc_munu_sub(ispin)%matrix, name="ADMM_XC", &
286  dist=sub_env%dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
287  row_blk_size=blk_sizes, col_blk_size=blk_sizes, nze=0)
288  CALL cp_dbcsr_alloc_block_from_nbl(a_xc_munu_sub(ispin)%matrix, sub_env%sab_aux_fit)
289  CALL dbcsr_set(a_xc_munu_sub(ispin)%matrix, 0.0_dp)
290  END DO
291 
292  CALL pw_env_get(sub_env%pw_env, auxbas_pw_pool=auxbas_pw_pool)
293  ALLOCATE (v_rspace_sub(nspins))
294  DO ispin = 1, nspins
295  CALL auxbas_pw_pool%create_pw(v_rspace_sub(ispin))
296  CALL pw_zero(v_rspace_sub(ispin))
297  END DO
298 
299  IF (admm_env%do_gapw) THEN
300  basis_type = "AUX_FIT_SOFT"
301  task_list => sub_env%task_list_aux_fit_soft
302  ELSE
303  basis_type = "AUX_FIT"
304  task_list => sub_env%task_list_aux_fit
305  END IF
306 
307  CALL tddfpt_apply_xc(a_ia_rspace=v_rspace_sub, &
308  kernel_env=kernel_env_admm_aux, &
309  rho_ia_struct=work_matrices%rho_aux_fit_struct_sub, &
310  is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
311  work_v_xc=work_matrices%wpw_rspace_sub, &
312  work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
313  DO ispin = 1, nspins
314  CALL pw_scale(v_rspace_sub(ispin), v_rspace_sub(ispin)%pw_grid%dvol)
315  CALL integrate_v_rspace(v_rspace=v_rspace_sub(ispin), &
316  hmat=a_xc_munu_sub(ispin), &
317  qs_env=qs_env, calculate_forces=.false., &
318  pw_env_external=sub_env%pw_env, &
319  basis_type=basis_type, &
320  task_list_external=task_list)
321  END DO
322  IF (admm_env%do_gapw) THEN
323  rho_atom_set => sub_env%local_rho_set_admm%rho_atom_set
324  rho1_atom_set => work_matrices%local_rho_set_admm%rho_atom_set
325  CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, &
326  kernel_env_admm_aux%xc_section, &
327  sub_env%para_env, do_tddfpt2=.true., &
328  do_triplet=is_rks_triplets, &
329  kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
330  CALL update_ks_atom(qs_env, a_xc_munu_sub, rho_ia_ao_aux_fit, forces=.false., tddft=.true., &
331  rho_atom_external=rho1_atom_set, &
332  kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
333  oce_external=admm_env%admm_gapw_env%oce, &
334  sab_external=sub_env%sab_aux_fit)
335  END IF
336  ALLOCATE (dbwork)
337  CALL dbcsr_create(dbwork, template=work_matrices%A_ia_munu_sub(1)%matrix)
338  CALL cp_fm_create(work_aux_orb, &
339  matrix_struct=work_matrices%wfm_aux_orb_sub%matrix_struct)
340  CALL cp_fm_create(work_orb_orb, &
341  matrix_struct=work_matrices%rho_ao_orb_fm_sub%matrix_struct)
342  CALL cp_fm_get_info(work_aux_orb, nrow_global=nao_aux, ncol_global=nao)
343  DO ispin = 1, nspins
344  CALL cp_dbcsr_sm_fm_multiply(a_xc_munu_sub(ispin)%matrix, sub_env%admm_A, &
345  work_aux_orb, nao)
346  CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, sub_env%admm_A, &
347  work_aux_orb, 0.0_dp, work_orb_orb)
348  CALL dbcsr_copy(dbwork, work_matrices%A_ia_munu_sub(1)%matrix)
349  CALL dbcsr_set(dbwork, 0.0_dp)
350  CALL copy_fm_to_dbcsr(work_orb_orb, dbwork, keep_sparsity=.true.)
351  CALL dbcsr_add(work_matrices%A_ia_munu_sub(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
352  END DO
353  CALL dbcsr_release(dbwork)
354  DEALLOCATE (dbwork)
355  DO ispin = 1, nspins
356  CALL auxbas_pw_pool%give_back_pw(v_rspace_sub(ispin))
357  END DO
358  DEALLOCATE (v_rspace_sub)
359  CALL cp_fm_release(work_aux_orb)
360  CALL cp_fm_release(work_orb_orb)
361  DO ispin = 1, nspins
362  CALL dbcsr_deallocate_matrix(a_xc_munu_sub(ispin)%matrix)
363  END DO
364  DEALLOCATE (a_xc_munu_sub)
365  ELSE
366  CALL tddfpt_apply_xc(a_ia_rspace=work_matrices%A_ia_rspace_sub, &
367  kernel_env=kernel_env_admm_aux, &
368  rho_ia_struct=work_matrices%rho_aux_fit_struct_sub, &
369  is_rks_triplets=is_rks_triplets, pw_env=sub_env%pw_env, &
370  work_v_xc=work_matrices%wpw_rspace_sub, &
371  work_v_xc_tau=work_matrices%wpw_tau_rspace_sub)
372  IF (admm_env%do_gapw) THEN
373  cpwarn("GAPW/ADMM needs symmetric ADMM kernel")
374  cpabort("GAPW/ADMM@TDDFT")
375  END IF
376  END IF
377  END IF
378  END IF
379 
380  ! electron-hole Coulomb interaction
381  IF (.NOT. is_rks_triplets) THEN
382  ! a sum J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu can be computed by solving
383  ! the Poisson equation for combined density (rho_{ia,alpha} + rho_{ia,beta}) .
384  ! The following action will destroy reciprocal-space grid in spin-unrestricted case.
385  DO ispin = 2, nspins
386  CALL pw_axpy(rho_ia_g(ispin), rho_ia_g(1))
387  END DO
388  CALL tddfpt_apply_coulomb(a_ia_rspace=work_matrices%A_ia_rspace_sub, &
389  rho_ia_g=rho_ia_g(1), &
390  local_rho_set=work_matrices%local_rho_set, &
391  hartree_local=work_matrices%hartree_local, &
392  qs_env=qs_env, sub_env=sub_env, gapw=gapw, &
393  work_v_gspace=work_matrices%wpw_gspace_sub(1), &
394  work_v_rspace=work_matrices%wpw_rspace_sub(1))
395  END IF
396 
397  ! convert from the plane-wave representation into the Gaussian basis set representation
398  DO ispin = 1, nspins
399  IF (.NOT. do_lrigpw) THEN
400  CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), &
401  work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol)
402 
403  IF (gapw) THEN
404  CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
405  hmat=work_matrices%A_ia_munu_sub(ispin), &
406  qs_env=qs_env, calculate_forces=.false., gapw=gapw, &
407  pw_env_external=sub_env%pw_env, &
408  task_list_external=sub_env%task_list_orb_soft)
409  ELSEIF (gapw_xc) THEN
410  IF (.NOT. is_rks_triplets) THEN
411  CALL integrate_v_rspace(v_rspace=work_matrices%A_ia_rspace_sub(ispin), &
412  hmat=work_matrices%A_ia_munu_sub(ispin), &
413  qs_env=qs_env, calculate_forces=.false., gapw=.false., &
414  pw_env_external=sub_env%pw_env, task_list_external=sub_env%task_list_orb)
415  END IF
416  ELSE
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  pw_env_external=sub_env%pw_env, task_list_external=sub_env%task_list_orb)
421  END IF
422  ELSE ! for full kernel using lri
423  CALL pw_scale(work_matrices%A_ia_rspace_sub(ispin), &
424  work_matrices%A_ia_rspace_sub(ispin)%pw_grid%dvol)
425  lri_v_int => kernel_env%lri_density%lri_coefs(ispin)%lri_kinds
426  CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
427  DO ikind = 1, nkind
428  lri_v_int(ikind)%v_int = 0.0_dp
429  END DO
430  CALL integrate_v_rspace_one_center(work_matrices%A_ia_rspace_sub(ispin), &
431  qs_env, lri_v_int, .false., "P_LRI_AUX")
432  DO ikind = 1, nkind
433  CALL para_env%sum(lri_v_int(ikind)%v_int)
434  END DO
435  END IF ! for full kernel using lri
436  END DO
437 
438  ! local atom contributions
439  IF (.NOT. do_lrigpw) THEN
440  IF (gapw .OR. gapw_xc) THEN
441  ! rho_ia_ao will not be touched
442  CALL update_ks_atom(qs_env, work_matrices%A_ia_munu_sub, rho_ia_ao, forces=.false., &
443  rho_atom_external=work_matrices%local_rho_set%rho_atom_set, &
444  tddft=.true.)
445  END IF
446  END IF
447 
448  ! calculate Coulomb contribution to response vector for lrigpw !
449  ! this is restricting lri to Coulomb only at the moment !
450  IF (do_lrigpw .AND. (.NOT. is_rks_triplets)) THEN !
451  CALL tddfpt2_lri_amat(qs_env, sub_env, kernel_env%lri_env, lri_v_int, work_matrices%A_ia_munu_sub)
452  END IF
453 
454  IF (ALLOCATED(work_matrices%evects_sub)) THEN
455  DO ispin = 1, nspins
456  CALL cp_dbcsr_sm_fm_multiply(work_matrices%A_ia_munu_sub(ispin)%matrix, &
457  sub_env%mos_occ(ispin), &
458  work_matrices%Aop_evects_sub(ispin, ivect), &
459  ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
460  END DO
461  ELSE
462  DO ispin = 1, nspins
463  CALL cp_dbcsr_sm_fm_multiply(work_matrices%A_ia_munu_sub(ispin)%matrix, &
464  sub_env%mos_occ(ispin), &
465  aop_evects(ispin, ivect), &
466  ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
467  END DO
468  END IF
469  END DO
470 
471  CALL timestop(handle)
472 
473  END SUBROUTINE fhxc_kernel
474 
475 ! **************************************************************************************************
476 !> \brief Compute action matrix-vector products with the sTDA Kernel
477 !> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
478 !> \param evects TDDFPT trial vectors
479 !> \param is_rks_triplets indicates that a triplet excited states calculation using
480 !> spin-unpolarised molecular orbitals has been requested
481 !> \param qs_env Quickstep environment
482 !> \param stda_control control parameters for sTDA kernel
483 !> \param stda_env ...
484 !> \param sub_env parallel (sub)group environment
485 !> \param work_matrices collection of work matrices (modified on exit)
486 !> \par History
487 !> * 04.2019 initial version [JHU]
488 ! **************************************************************************************************
489  SUBROUTINE stda_kernel(Aop_evects, evects, is_rks_triplets, &
490  qs_env, stda_control, stda_env, &
491  sub_env, work_matrices)
492 
493  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: aop_evects, evects
494  LOGICAL, INTENT(in) :: is_rks_triplets
495  TYPE(qs_environment_type), POINTER :: qs_env
496  TYPE(stda_control_type) :: stda_control
497  TYPE(stda_env_type) :: stda_env
498  TYPE(tddfpt_subgroup_env_type) :: sub_env
499  TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices
500 
501  CHARACTER(LEN=*), PARAMETER :: routinen = 'stda_kernel'
502 
503  INTEGER :: handle, ivect, nvects
504 
505  CALL timeset(routinen, handle)
506 
507  nvects = SIZE(evects, 2)
508 
509  DO ivect = 1, nvects
510  IF (ALLOCATED(work_matrices%evects_sub)) THEN
511  IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct)) THEN
512  CALL stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, work_matrices, &
513  is_rks_triplets, work_matrices%evects_sub(:, ivect), &
514  work_matrices%Aop_evects_sub(:, ivect))
515  ELSE
516  ! skip trial vectors which are assigned to different parallel groups
517  cycle
518  END IF
519  ELSE
520  CALL stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, work_matrices, &
521  is_rks_triplets, evects(:, ivect), aop_evects(:, ivect))
522  END IF
523  END DO
524 
525  CALL timestop(handle)
526 
527  END SUBROUTINE stda_kernel
528 
529 ! **************************************************************************************************
530 
531 END MODULE qs_tddfpt2_fhxc
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition: admm_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
Definition: cp_dbcsr_operations.F:887
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition: cp_dbcsr_operations.F:118
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_admm_aux_exch_func_none
integer, parameter, public do_admm_aux_exch_func_none
Definition: input_constants.F:863
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition: lri_environment_types.F:18
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_methods
Definition: pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_gapw_densities
Definition: qs_gapw_densities.F:9
qs_gapw_densities::prepare_gapw_den
subroutine, public prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external)
...
Definition: qs_gapw_densities.F:49
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition: qs_integrate_potential.F:22
qs_kernel_types
Definition: qs_kernel_types.F:8
qs_ks_atom
routines that build the Kohn-Sham matrix contributions coming from local atomic densities
Definition: qs_ks_atom.F:12
qs_ks_atom::update_ks_atom
subroutine, public update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, kind_set_external, oce_external, sab_external, kscale, kintegral, kforce, fscale)
The correction to the KS matrix due to the GAPW local terms to the hartree and XC contributions is he...
Definition: qs_ks_atom.F:110
qs_rho_atom_types
Definition: qs_rho_atom_types.F:8
qs_rho_methods
methods of the rho structure (defined in qs_rho_types)
Definition: qs_rho_methods.F:15
qs_rho_methods::qs_rho_update_tddfpt
subroutine, public qs_rho_update_tddfpt(rho_struct, qs_env, pw_env_external, task_list_external, para_env_external, tddfpt_lri_env, tddfpt_lri_density)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g
Definition: qs_rho_methods.F:670
qs_rho_methods::qs_rho_update_rho
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
Definition: qs_rho_methods.F:375
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
qs_tddfpt2_densities
Definition: qs_tddfpt2_densities.F:8
qs_tddfpt2_densities::tddfpt_construct_aux_fit_density
subroutine, public tddfpt_construct_aux_fit_density(rho_orb_struct, rho_aux_fit_struct, local_rho_set, qs_env, sub_env, wfm_rho_orb, wfm_rho_aux_fit, wfm_aux_orb)
Project a charge density expressed in primary basis set into the auxiliary basis set.
Definition: qs_tddfpt2_densities.F:160
qs_tddfpt2_fhxc
Definition: qs_tddfpt2_fhxc.F:8
qs_tddfpt2_fhxc::stda_kernel
subroutine, public stda_kernel(Aop_evects, evects, is_rks_triplets, qs_env, stda_control, stda_env, sub_env, work_matrices)
Compute action matrix-vector products with the sTDA Kernel.
Definition: qs_tddfpt2_fhxc.F:492
qs_tddfpt2_fhxc::fhxc_kernel
subroutine, public fhxc_kernel(Aop_evects, evects, is_rks_triplets, do_hfx, do_admm, qs_env, kernel_env, kernel_env_admm_aux, sub_env, work_matrices, admm_symm, admm_xc_correction, do_lrigpw)
Compute action matrix-vector products with the FHxc Kernel.
Definition: qs_tddfpt2_fhxc.F:98
qs_tddfpt2_lri_utils
Definition: qs_tddfpt2_lri_utils.F:8
qs_tddfpt2_lri_utils::tddfpt2_lri_amat
subroutine, public tddfpt2_lri_amat(qs_env, sub_env, lri_env, lri_v_int, A_ia_munu_sub)
Calculate contribution to response vector for LRI.
Definition: qs_tddfpt2_lri_utils.F:251
qs_tddfpt2_operators
Definition: qs_tddfpt2_operators.F:8
qs_tddfpt2_operators::tddfpt_apply_xc_potential
subroutine, public tddfpt_apply_xc_potential(A_ia_rspace, fxc_rspace, rho_ia_struct, is_rks_triplets)
Routine for applying fxc potential.
Definition: qs_tddfpt2_operators.F:284
qs_tddfpt2_operators::tddfpt_apply_xc
subroutine, public tddfpt_apply_xc(A_ia_rspace, kernel_env, rho_ia_struct, is_rks_triplets, pw_env, work_v_xc, work_v_xc_tau)
Driver routine for applying fxc (analyic vs. finite difference for testing.
Definition: qs_tddfpt2_operators.F:248
qs_tddfpt2_operators::tddfpt_apply_coulomb
subroutine, public tddfpt_apply_coulomb(A_ia_rspace, rho_ia_g, local_rho_set, hartree_local, qs_env, sub_env, gapw, work_v_gspace, work_v_rspace)
Update v_rspace by adding coulomb term.
Definition: qs_tddfpt2_operators.F:172
qs_tddfpt2_stda_types
Simplified Tamm Dancoff approach (sTDA).
Definition: qs_tddfpt2_stda_types.F:10
qs_tddfpt2_stda_utils
Simplified Tamm Dancoff approach (sTDA).
Definition: qs_tddfpt2_stda_utils.F:10
qs_tddfpt2_stda_utils::stda_calculate_kernel
subroutine, public stda_calculate_kernel(qs_env, stda_control, stda_env, sub_env, work, is_rks_triplets, X, res)
...Calculate the sTDA kernel contribution by contracting the Lowdin MO coefficients – transition char...
Definition: qs_tddfpt2_stda_utils.F:508
qs_tddfpt2_subgroups
Definition: qs_tddfpt2_subgroups.F:8
qs_tddfpt2_types
Definition: qs_tddfpt2_types.F:8
qs_vxc_atom
routines that build the integrals of the Vxc potential calculated for the atomic density in the basis...
Definition: qs_vxc_atom.F:12
qs_vxc_atom::calculate_xc_2nd_deriv_atom
subroutine, public calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, do_tddft, do_tddfpt2, do_triplet, kind_set_external)
...
Definition: qs_vxc_atom.F:446
task_list_types
types for task lists
Definition: task_list_types.F:14