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_operators
9 : USE admm_types, ONLY: admm_type
10 : USE cell_types, ONLY: cell_type,&
11 : pbc
12 : USE cp_dbcsr_api, ONLY: &
13 : dbcsr_create, dbcsr_filter, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
14 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
15 : dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
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_basic_linalg, ONLY: cp_fm_column_scale,&
20 : cp_fm_scale_and_add
21 : USE cp_fm_struct, ONLY: cp_fm_struct_type
22 : USE cp_fm_types, ONLY: cp_fm_create,&
23 : cp_fm_get_info,&
24 : cp_fm_release,&
25 : cp_fm_to_fm,&
26 : cp_fm_type
27 : USE hartree_local_methods, ONLY: Vh_1c_gg_integrals
28 : USE hartree_local_types, ONLY: hartree_local_type
29 : USE hfx_admm_utils, ONLY: tddft_hfx_matrix
30 : USE hfx_types, ONLY: hfx_type
31 : USE input_section_types, ONLY: section_vals_get,&
32 : section_vals_get_subs_vals,&
33 : section_vals_type
34 : USE kinds, ONLY: dp
35 : USE message_passing, ONLY: mp_para_env_type
36 : USE parallel_gemm_api, ONLY: parallel_gemm
37 : USE particle_types, ONLY: particle_type
38 : USE pw_env_types, ONLY: pw_env_get,&
39 : pw_env_type
40 : USE pw_methods, ONLY: pw_axpy,&
41 : pw_multiply,&
42 : pw_scale,&
43 : pw_transfer,&
44 : pw_zero
45 : USE pw_poisson_methods, ONLY: pw_poisson_solve
46 : USE pw_poisson_types, ONLY: pw_poisson_type
47 : USE pw_pool_types, ONLY: pw_pool_p_type,&
48 : pw_pool_type
49 : USE pw_types, ONLY: pw_c1d_gs_type,&
50 : pw_r3d_rs_type
51 : USE qs_environment_types, ONLY: get_qs_env,&
52 : qs_environment_type
53 : USE qs_kernel_types, ONLY: full_kernel_env_type
54 : USE qs_local_rho_types, ONLY: local_rho_type
55 : USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace
56 : USE qs_rho_types, ONLY: qs_rho_get,&
57 : qs_rho_type
58 : USE qs_tddfpt2_stda_utils, ONLY: get_lowdin_x
59 : USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
60 : USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos,&
61 : tddfpt_work_matrices
62 : USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,&
63 : realspace_grid_type
64 : USE xc, ONLY: xc_calc_2nd_deriv_analytical,&
65 : xc_calc_2nd_deriv_numerical
66 : USE xc_rho_set_types, ONLY: xc_rho_set_type,&
67 : xc_rho_set_update
68 : #include "./base/base_uses.f90"
69 :
70 : IMPLICIT NONE
71 :
72 : PRIVATE
73 :
74 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_operators'
75 :
76 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
77 : ! number of first derivative components (3: d/dx, d/dy, d/dz)
78 : INTEGER, PARAMETER, PRIVATE :: nderivs = 3
79 : INTEGER, PARAMETER, PRIVATE :: maxspins = 2
80 :
81 : PUBLIC :: tddfpt_apply_energy_diff, tddfpt_apply_coulomb, tddfpt_apply_xc, tddfpt_apply_hfx, &
82 : tddfpt_apply_xc_potential, tddfpt_apply_hfxlr_kernel, tddfpt_apply_hfxsr_kernel
83 :
84 : ! **************************************************************************************************
85 :
86 : CONTAINS
87 :
88 : ! **************************************************************************************************
89 : !> \brief Apply orbital energy difference term:
90 : !> Aop_evects(spin,state) += KS(spin) * evects(spin,state) -
91 : !> S * evects(spin,state) * diag(evals_occ(spin))
92 : !> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
93 : !> \param evects trial vectors C_{1,i}
94 : !> \param S_evects S * C_{1,i}
95 : !> \param gs_mos molecular orbitals optimised for the ground state (only occupied orbital
96 : !> energies [component %evals_occ] are needed)
97 : !> \param matrix_ks Kohn-Sham matrix
98 : !> \par History
99 : !> * 05.2016 initialise all matrix elements in one go [Sergey Chulkov]
100 : !> * 03.2017 renamed from tddfpt_init_energy_diff(), altered prototype [Sergey Chulkov]
101 : !> \note Based on the subroutine p_op_l1() which was originally created by
102 : !> Thomas Chassaing on 08.2002.
103 : ! **************************************************************************************************
104 5218 : SUBROUTINE tddfpt_apply_energy_diff(Aop_evects, evects, S_evects, gs_mos, matrix_ks)
105 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects
106 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects, S_evects
107 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
108 : INTENT(in) :: gs_mos
109 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_ks
110 :
111 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_energy_diff'
112 :
113 : INTEGER :: handle, ispin, ivect, nactive, nao, &
114 : nspins, nvects
115 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
116 : TYPE(cp_fm_type) :: hevec
117 :
118 5218 : CALL timeset(routineN, handle)
119 :
120 5218 : nspins = SIZE(evects, 1)
121 5218 : nvects = SIZE(evects, 2)
122 :
123 11114 : DO ispin = 1, nspins
124 : CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
125 5896 : nrow_global=nao, ncol_global=nactive)
126 5896 : CALL cp_fm_create(hevec, matrix_struct)
127 :
128 21254 : DO ivect = 1, nvects
129 : CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, evects(ispin, ivect), &
130 : Aop_evects(ispin, ivect), ncol=nactive, &
131 15358 : alpha=1.0_dp, beta=1.0_dp)
132 :
133 15358 : IF (ASSOCIATED(gs_mos(ispin)%evals_occ_matrix)) THEN
134 : ! orbital energy correction: evals_occ_matrix is not a diagonal matrix
135 : CALL parallel_gemm('N', 'N', nao, nactive, nactive, 1.0_dp, &
136 : S_evects(ispin, ivect), gs_mos(ispin)%evals_occ_matrix, &
137 756 : 0.0_dp, hevec)
138 : ELSE
139 14602 : CALL cp_fm_to_fm(S_evects(ispin, ivect), hevec)
140 14602 : CALL cp_fm_column_scale(hevec, gs_mos(ispin)%evals_occ)
141 : END IF
142 :
143 : ! KS * C1 - S * C1 * occupied_orbital_energies
144 21254 : CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect), -1.0_dp, hevec)
145 : END DO
146 17010 : CALL cp_fm_release(hevec)
147 : END DO
148 :
149 5218 : CALL timestop(handle)
150 :
151 5218 : END SUBROUTINE tddfpt_apply_energy_diff
152 :
153 : ! **************************************************************************************************
154 : !> \brief Update v_rspace by adding coulomb term.
155 : !> \param A_ia_rspace action of TDDFPT operator on the trial vector expressed in a plane wave
156 : !> representation (modified on exit)
157 : !> \param rho_ia_g response density in reciprocal space for the given trial vector
158 : !> \param local_rho_set ...
159 : !> \param hartree_local ...
160 : !> \param qs_env ...
161 : !> \param sub_env the full sub_environment needed for calculation
162 : !> \param gapw Flag indicating GAPW cacluation
163 : !> \param work_v_gspace work reciprocal-space grid to store Coulomb potential (modified on exit)
164 : !> \param work_v_rspace work real-space grid to store Coulomb potential (modified on exit)
165 : !> \param tddfpt_mgrid ...
166 : !> \par History
167 : !> * 05.2016 compute all coulomb terms in one go [Sergey Chulkov]
168 : !> * 03.2017 proceed excited states sequentially; minimise the number of conversions between
169 : !> DBCSR and FM matrices [Sergey Chulkov]
170 : !> * 06.2018 return the action expressed in the plane wave representation instead of the one
171 : !> in the atomic basis set representation
172 : !> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by
173 : !> Mohamed Fawzi on 10.2002.
174 : ! **************************************************************************************************
175 6070 : SUBROUTINE tddfpt_apply_coulomb(A_ia_rspace, rho_ia_g, local_rho_set, hartree_local, &
176 : qs_env, sub_env, gapw, work_v_gspace, work_v_rspace, tddfpt_mgrid)
177 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT) :: A_ia_rspace
178 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_ia_g
179 : TYPE(local_rho_type), POINTER :: local_rho_set
180 : TYPE(hartree_local_type), POINTER :: hartree_local
181 : TYPE(qs_environment_type), POINTER :: qs_env
182 : TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
183 : LOGICAL, INTENT(IN) :: gapw
184 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: work_v_gspace
185 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: work_v_rspace
186 : LOGICAL, INTENT(IN) :: tddfpt_mgrid
187 :
188 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_coulomb'
189 :
190 : INTEGER :: handle, ispin, nspins
191 : REAL(kind=dp) :: alpha, pair_energy
192 : TYPE(pw_env_type), POINTER :: pw_env
193 : TYPE(pw_poisson_type), POINTER :: poisson_env
194 6070 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: my_pools
195 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
196 6070 : POINTER :: my_rs_descs
197 6070 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: my_rs_grids
198 :
199 6070 : CALL timeset(routineN, handle)
200 :
201 6070 : nspins = SIZE(A_ia_rspace)
202 6070 : pw_env => sub_env%pw_env
203 6070 : IF (tddfpt_mgrid) THEN
204 : CALL pw_env_get(pw_env, poisson_env=poisson_env, rs_grids=my_rs_grids, &
205 86 : rs_descs=my_rs_descs, pw_pools=my_pools)
206 : ELSE
207 5984 : CALL pw_env_get(pw_env, poisson_env=poisson_env)
208 : END IF
209 :
210 6070 : IF (nspins > 1) THEN
211 1358 : alpha = 1.0_dp
212 : ELSE
213 : ! spin-restricted case: alpha == 2 due to singlet state.
214 : ! In case of triplet states alpha == 0, so we should not call this subroutine at all.
215 4712 : alpha = 2.0_dp
216 : END IF
217 :
218 6070 : IF (gapw) THEN
219 808 : CPASSERT(ASSOCIATED(local_rho_set))
220 808 : CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rho_ia_g)
221 : END IF
222 :
223 6070 : CALL pw_poisson_solve(poisson_env, rho_ia_g, pair_energy, work_v_gspace)
224 6070 : CALL pw_transfer(work_v_gspace, work_v_rspace)
225 :
226 : ! (i a || j b) = ( i_alpha a_alpha + i_beta a_beta || j_alpha b_alpha + j_beta b_beta) =
227 : ! tr (Cj_alpha^T * [J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu] * Cb_alpha) +
228 : ! tr (Cj_beta^T * [J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu] * Cb_beta)
229 13498 : DO ispin = 1, nspins
230 13498 : CALL pw_axpy(work_v_rspace, A_ia_rspace(ispin), alpha)
231 : END DO
232 :
233 6070 : IF (gapw) THEN
234 : CALL Vh_1c_gg_integrals(qs_env, pair_energy, &
235 : hartree_local%ecoul_1c, &
236 : local_rho_set, &
237 808 : sub_env%para_env, tddft=.TRUE., core_2nd=.TRUE.)
238 808 : CALL pw_scale(work_v_rspace, work_v_rspace%pw_grid%dvol)
239 808 : IF (tddfpt_mgrid) THEN
240 : CALL integrate_vhg0_rspace(qs_env, work_v_rspace, sub_env%para_env, &
241 : calculate_forces=.FALSE., &
242 : local_rho_set=local_rho_set, my_pools=my_pools, &
243 50 : my_rs_descs=my_rs_descs)
244 : ELSE
245 : CALL integrate_vhg0_rspace(qs_env, work_v_rspace, sub_env%para_env, &
246 : calculate_forces=.FALSE., &
247 758 : local_rho_set=local_rho_set)
248 : END IF
249 : END IF
250 :
251 6070 : CALL timestop(handle)
252 :
253 6070 : END SUBROUTINE tddfpt_apply_coulomb
254 :
255 : ! **************************************************************************************************
256 : !> \brief Driver routine for applying fxc (analyic vs. finite difference for testing
257 : !> \param A_ia_rspace action of TDDFPT operator on trial vectors expressed in a plane wave
258 : !> representation (modified on exit)
259 : !> \param kernel_env kernel environment
260 : !> \param rho_ia_struct response density for the given trial vector
261 : !> \param is_rks_triplets indicates that the triplet excited states calculation using
262 : !> spin-unpolarised molecular orbitals has been requested
263 : !> \param pw_env plain wave environment
264 : !> \param work_v_xc work real-space grid to store the gradient of the exchange-correlation
265 : !> potential with respect to the response density (modified on exit)
266 : !> \param work_v_xc_tau ...
267 : ! **************************************************************************************************
268 7224 : SUBROUTINE tddfpt_apply_xc(A_ia_rspace, kernel_env, rho_ia_struct, is_rks_triplets, &
269 : pw_env, work_v_xc, work_v_xc_tau)
270 :
271 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT) :: A_ia_rspace
272 : TYPE(full_kernel_env_type), INTENT(IN) :: kernel_env
273 : TYPE(qs_rho_type), POINTER :: rho_ia_struct
274 : LOGICAL, INTENT(in) :: is_rks_triplets
275 : TYPE(pw_env_type), POINTER :: pw_env
276 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: work_v_xc, work_v_xc_tau
277 :
278 : INTEGER :: ispin, nspins
279 :
280 7224 : nspins = SIZE(A_ia_rspace)
281 :
282 7224 : IF (kernel_env%deriv2_analytic) THEN
283 : CALL tddfpt_apply_xc_analytic(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
284 7184 : pw_env, work_v_xc, work_v_xc_tau)
285 : ELSE
286 : CALL tddfpt_apply_xc_fd(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
287 40 : pw_env, work_v_xc, work_v_xc_tau)
288 : END IF
289 :
290 15806 : DO ispin = 1, nspins
291 : ! pw2 = pw2 + alpha * pw1
292 15806 : CALL pw_axpy(work_v_xc(ispin), A_ia_rspace(ispin), kernel_env%alpha)
293 : END DO
294 :
295 7224 : END SUBROUTINE tddfpt_apply_xc
296 :
297 : ! **************************************************************************************************
298 : !> \brief Routine for applying fxc potential
299 : !> \param A_ia_rspace action of TDDFPT operator on trial vectors expressed in a plane wave
300 : !> representation (modified on exit)
301 : !> \param fxc_rspace ...
302 : !> \param rho_ia_struct response density for the given trial vector
303 : !> \param is_rks_triplets ...
304 : ! **************************************************************************************************
305 156 : SUBROUTINE tddfpt_apply_xc_potential(A_ia_rspace, fxc_rspace, rho_ia_struct, is_rks_triplets)
306 :
307 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT) :: A_ia_rspace
308 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rspace
309 : TYPE(qs_rho_type), POINTER :: rho_ia_struct
310 : LOGICAL, INTENT(in) :: is_rks_triplets
311 :
312 : INTEGER :: nspins
313 : REAL(KIND=dp) :: alpha
314 156 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r
315 :
316 156 : nspins = SIZE(A_ia_rspace)
317 :
318 156 : alpha = 1.0_dp
319 :
320 156 : CALL qs_rho_get(rho_ia_struct, rho_r=rho1_r)
321 :
322 156 : IF (nspins == 2) THEN
323 0 : CALL pw_multiply(A_ia_rspace(1), fxc_rspace(1), rho1_r(1), alpha)
324 0 : CALL pw_multiply(A_ia_rspace(1), fxc_rspace(2), rho1_r(2), alpha)
325 0 : CALL pw_multiply(A_ia_rspace(2), fxc_rspace(3), rho1_r(2), alpha)
326 0 : CALL pw_multiply(A_ia_rspace(2), fxc_rspace(2), rho1_r(1), alpha)
327 156 : ELSE IF (is_rks_triplets) THEN
328 0 : CALL pw_multiply(A_ia_rspace(1), fxc_rspace(1), rho1_r(1), alpha)
329 0 : CALL pw_multiply(A_ia_rspace(1), fxc_rspace(2), rho1_r(1), -alpha)
330 : ELSE
331 156 : CALL pw_multiply(A_ia_rspace(1), fxc_rspace(1), rho1_r(1), alpha)
332 156 : CALL pw_multiply(A_ia_rspace(1), fxc_rspace(2), rho1_r(1), alpha)
333 : END IF
334 :
335 156 : END SUBROUTINE tddfpt_apply_xc_potential
336 :
337 : ! **************************************************************************************************
338 : !> \brief Update A_ia_munu by adding exchange-correlation term.
339 : !> \param kernel_env kernel environment
340 : !> \param rho_ia_struct response density for the given trial vector
341 : !> \param is_rks_triplets indicates that the triplet excited states calculation using
342 : !> spin-unpolarised molecular orbitals has been requested
343 : !> \param nspins ...
344 : !> \param pw_env plain wave environment
345 : !> \param work_v_xc work real-space grid to store the gradient of the exchange-correlation
346 : !> potential with respect to the response density (modified on exit)
347 : !> \param work_v_xc_tau ...
348 : !> \par History
349 : !> * 05.2016 compute all kernel terms in one go [Sergey Chulkov]
350 : !> * 03.2017 proceed excited states sequentially; minimise the number of conversions between
351 : !> DBCSR and FM matrices [Sergey Chulkov]
352 : !> * 06.2018 return the action expressed in the plane wave representation instead of the one
353 : !> in the atomic basis set representation
354 : !> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by
355 : !> Mohamed Fawzi on 10.2002.
356 : ! **************************************************************************************************
357 7184 : SUBROUTINE tddfpt_apply_xc_analytic(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
358 : pw_env, work_v_xc, work_v_xc_tau)
359 : TYPE(full_kernel_env_type), INTENT(in) :: kernel_env
360 : TYPE(qs_rho_type), POINTER :: rho_ia_struct
361 : LOGICAL, INTENT(in) :: is_rks_triplets
362 : INTEGER, INTENT(in) :: nspins
363 : TYPE(pw_env_type), POINTER :: pw_env
364 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: work_v_xc, work_v_xc_tau
365 :
366 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_xc_analytic'
367 :
368 : INTEGER :: handle, ispin
369 7184 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_ia_g, rho_ia_g2
370 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
371 7184 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_ia_r, rho_ia_r2, tau_ia_r, tau_ia_r2
372 :
373 7184 : CALL timeset(routineN, handle)
374 :
375 7184 : CALL qs_rho_get(rho_ia_struct, rho_g=rho_ia_g, rho_r=rho_ia_r, tau_r=tau_ia_r)
376 7184 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
377 :
378 : IF (debug_this_module) THEN
379 : CPASSERT(SIZE(rho_ia_g) == nspins)
380 : CPASSERT(SIZE(rho_ia_r) == nspins)
381 : CPASSERT((.NOT. ASSOCIATED(tau_ia_r)) .OR. SIZE(tau_ia_r) == nspins)
382 : CPASSERT((.NOT. is_rks_triplets) .OR. nspins == 1)
383 : END IF
384 :
385 7184 : NULLIFY (tau_ia_r2)
386 7184 : IF (is_rks_triplets) THEN
387 1512 : ALLOCATE (rho_ia_r2(2))
388 1512 : ALLOCATE (rho_ia_g2(2))
389 504 : rho_ia_r2(1) = rho_ia_r(1)
390 504 : rho_ia_r2(2) = rho_ia_r(1)
391 504 : rho_ia_g2(1) = rho_ia_g(1)
392 504 : rho_ia_g2(2) = rho_ia_g(1)
393 :
394 504 : IF (ASSOCIATED(tau_ia_r)) THEN
395 0 : ALLOCATE (tau_ia_r2(2))
396 0 : tau_ia_r2(1) = tau_ia_r(1)
397 0 : tau_ia_r2(2) = tau_ia_r(1)
398 : END IF
399 : ELSE
400 6680 : rho_ia_r2 => rho_ia_r
401 6680 : rho_ia_g2 => rho_ia_g
402 :
403 6680 : tau_ia_r2 => tau_ia_r
404 : END IF
405 :
406 15706 : DO ispin = 1, nspins
407 8522 : CALL pw_zero(work_v_xc(ispin))
408 15706 : IF (ASSOCIATED(work_v_xc_tau)) CALL pw_zero(work_v_xc_tau(ispin))
409 : END DO
410 :
411 : CALL xc_rho_set_update(rho_set=kernel_env%xc_rho1_set, rho_r=rho_ia_r2, rho_g=rho_ia_g2, tau=tau_ia_r2, &
412 : needs=kernel_env%xc_rho1_cflags, xc_deriv_method_id=kernel_env%deriv_method_id, &
413 7184 : xc_rho_smooth_id=kernel_env%rho_smooth_id, pw_pool=auxbas_pw_pool)
414 :
415 : CALL xc_calc_2nd_deriv_analytical(v_xc=work_v_xc, v_xc_tau=work_v_xc_tau, deriv_set=kernel_env%xc_deriv_set, &
416 : rho_set=kernel_env%xc_rho_set, &
417 : rho1_set=kernel_env%xc_rho1_set, pw_pool=auxbas_pw_pool, &
418 7184 : xc_section=kernel_env%xc_section, gapw=.FALSE., tddfpt_fac=kernel_env%beta)
419 :
420 7184 : IF (is_rks_triplets) THEN
421 504 : DEALLOCATE (rho_ia_r2)
422 504 : DEALLOCATE (rho_ia_g2)
423 504 : IF (ASSOCIATED(tau_ia_r2)) DEALLOCATE (tau_ia_r2)
424 : END IF
425 :
426 7184 : CALL timestop(handle)
427 :
428 7184 : END SUBROUTINE tddfpt_apply_xc_analytic
429 :
430 : ! **************************************************************************************************
431 : !> \brief Update A_ia_munu by adding exchange-correlation term using finite difference methods.
432 : !> \param kernel_env kernel environment
433 : !> \param rho_ia_struct response density for the given trial vector
434 : !> \param is_rks_triplets indicates that the triplet excited states calculation using
435 : !> spin-unpolarised molecular orbitals has been requested
436 : !> \param nspins ...
437 : !> \param pw_env plain wave environment
438 : !> \param work_v_xc work real-space grid to store the gradient of the exchange-correlation
439 : !> potential with respect to the response density (modified on exit)
440 : !> \param work_v_xc_tau ...
441 : ! **************************************************************************************************
442 40 : SUBROUTINE tddfpt_apply_xc_fd(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
443 : pw_env, work_v_xc, work_v_xc_tau)
444 : TYPE(full_kernel_env_type), INTENT(in) :: kernel_env
445 : TYPE(qs_rho_type), POINTER :: rho_ia_struct
446 : LOGICAL, INTENT(in) :: is_rks_triplets
447 : INTEGER, INTENT(in) :: nspins
448 : TYPE(pw_env_type), POINTER :: pw_env
449 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: work_v_xc, work_v_xc_tau
450 :
451 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_xc_fd'
452 :
453 : INTEGER :: handle, ispin
454 : LOGICAL :: lsd, singlet, triplet
455 40 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g
456 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
457 40 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r
458 : TYPE(xc_rho_set_type), POINTER :: rho_set
459 :
460 40 : CALL timeset(routineN, handle)
461 :
462 40 : CALL qs_rho_get(rho_ia_struct, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
463 40 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
464 100 : DO ispin = 1, nspins
465 100 : CALL pw_zero(work_v_xc(ispin))
466 : END DO
467 40 : rho_set => kernel_env%xc_rho_set
468 :
469 40 : singlet = .FALSE.
470 40 : triplet = .FALSE.
471 40 : lsd = .FALSE.
472 40 : IF (nspins == 1 .AND. .NOT. is_rks_triplets) THEN
473 40 : singlet = .TRUE.
474 40 : ELSE IF (nspins == 1 .AND. is_rks_triplets) THEN
475 40 : triplet = .TRUE.
476 20 : ELSE IF (nspins == 2) THEN
477 40 : lsd = .TRUE.
478 : ELSE
479 0 : CPABORT("illegal options")
480 : END IF
481 :
482 40 : IF (ASSOCIATED(tau1_r)) THEN
483 0 : DO ispin = 1, nspins
484 0 : CALL pw_zero(work_v_xc_tau(ispin))
485 : END DO
486 : END IF
487 :
488 : CALL xc_calc_2nd_deriv_numerical(work_v_xc, work_v_xc_tau, rho_set, rho1_r, rho1_g, tau1_r, &
489 : auxbas_pw_pool, kernel_env%xc_section, &
490 40 : is_rks_triplets)
491 :
492 40 : CALL timestop(handle)
493 :
494 40 : END SUBROUTINE tddfpt_apply_xc_fd
495 :
496 : ! **************************************************************************************************
497 : !> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
498 : !> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
499 : !> \param evects trial vectors
500 : !> \param gs_mos molecular orbitals optimised for the ground state (only occupied
501 : !> molecular orbitals [component %mos_occ] are needed)
502 : !> \param do_admm perform auxiliary density matrix method calculations
503 : !> \param qs_env Quickstep environment
504 : !> \param work_rho_ia_ao_symm ...
505 : !> \param work_hmat_symm ...
506 : !> \param work_rho_ia_ao_asymm ...
507 : !> \param work_hmat_asymm ...
508 : !> \param wfm_rho_orb ...
509 : !> \par History
510 : !> * 05.2016 compute all exact-exchange terms in one go [Sergey Chulkov]
511 : !> * 03.2017 code related to ADMM correction is now moved to tddfpt_apply_admm_correction()
512 : !> in order to compute this correction within parallel groups [Sergey Chulkov]
513 : !> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by
514 : !> Mohamed Fawzi on 10.2002.
515 : ! **************************************************************************************************
516 1118 : SUBROUTINE tddfpt_apply_hfx(Aop_evects, evects, gs_mos, do_admm, qs_env, &
517 1118 : work_rho_ia_ao_symm, work_hmat_symm, work_rho_ia_ao_asymm, &
518 1118 : work_hmat_asymm, wfm_rho_orb)
519 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects
520 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects
521 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
522 : INTENT(in) :: gs_mos
523 : LOGICAL, INTENT(in) :: do_admm
524 : TYPE(qs_environment_type), POINTER :: qs_env
525 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: work_rho_ia_ao_symm
526 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
527 : TARGET :: work_hmat_symm
528 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: work_rho_ia_ao_asymm
529 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
530 : TARGET :: work_hmat_asymm
531 : TYPE(cp_fm_type), INTENT(IN) :: wfm_rho_orb
532 :
533 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_hfx'
534 :
535 : INTEGER :: handle, ispin, ivect, nao, nao_aux, &
536 : nspins, nvects
537 : INTEGER, DIMENSION(maxspins) :: nactive
538 : LOGICAL :: do_hfx
539 : REAL(kind=dp) :: alpha
540 : TYPE(admm_type), POINTER :: admm_env
541 : TYPE(section_vals_type), POINTER :: hfx_section, input
542 :
543 1118 : CALL timeset(routineN, handle)
544 :
545 : ! Check for hfx section
546 1118 : CALL get_qs_env(qs_env, input=input)
547 1118 : hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
548 1118 : CALL section_vals_get(hfx_section, explicit=do_hfx)
549 :
550 1118 : IF (do_hfx) THEN
551 1118 : nspins = SIZE(evects, 1)
552 1118 : nvects = SIZE(evects, 2)
553 :
554 1118 : IF (nspins > 1) THEN
555 78 : alpha = 1.0_dp
556 : ELSE
557 1040 : alpha = 2.0_dp
558 : END IF
559 :
560 1118 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
561 2314 : DO ispin = 1, nspins
562 2314 : CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
563 : END DO
564 :
565 1118 : IF (do_admm) THEN
566 612 : CALL get_qs_env(qs_env, admm_env=admm_env)
567 612 : CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux)
568 : END IF
569 :
570 : !Note: the symmetrized transition density matrix is P = 0.5*(C*evect^T + evect*C^T)
571 : ! in the end, we only want evect*C^T for consistency with the MO formulation of TDDFT
572 : ! therefore, we go in 2 steps: with the symmetric 0.5*(C*evect^T + evect*C^T) and
573 : ! the antisymemtric 0.5*(C*evect^T - evect*C^T)
574 :
575 : ! some stuff from qs_ks_build_kohn_sham_matrix
576 : ! TO DO: add SIC support
577 3426 : DO ivect = 1, nvects
578 4744 : DO ispin = 1, nspins
579 :
580 : !The symmetric density matrix
581 : CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect), &
582 2436 : gs_mos(ispin)%mos_occ, 0.0_dp, wfm_rho_orb)
583 : CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, gs_mos(ispin)%mos_occ, &
584 2436 : evects(ispin, ivect), 1.0_dp, wfm_rho_orb)
585 :
586 2436 : CALL dbcsr_set(work_hmat_symm(ispin)%matrix, 0.0_dp)
587 4744 : IF (do_admm) THEN
588 : CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, &
589 1300 : wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb)
590 : CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
591 1300 : 0.0_dp, admm_env%work_aux_aux)
592 1300 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao_symm(ispin)%matrix, keep_sparsity=.TRUE.)
593 : ELSE
594 1136 : CALL copy_fm_to_dbcsr(wfm_rho_orb, work_rho_ia_ao_symm(ispin)%matrix, keep_sparsity=.TRUE.)
595 : END IF
596 : END DO
597 :
598 2308 : CALL tddft_hfx_matrix(work_hmat_symm, work_rho_ia_ao_symm, qs_env)
599 :
600 2308 : IF (do_admm) THEN
601 2572 : DO ispin = 1, nspins
602 : CALL cp_dbcsr_sm_fm_multiply(work_hmat_symm(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
603 1300 : ncol=nao, alpha=1.0_dp, beta=0.0_dp)
604 :
605 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
606 1300 : admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb)
607 :
608 : CALL parallel_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, &
609 2572 : gs_mos(ispin)%mos_occ, 1.0_dp, Aop_evects(ispin, ivect))
610 : END DO
611 : ELSE
612 2172 : DO ispin = 1, nspins
613 : CALL cp_dbcsr_sm_fm_multiply(work_hmat_symm(ispin)%matrix, gs_mos(ispin)%mos_occ, &
614 : Aop_evects(ispin, ivect), ncol=nactive(ispin), &
615 2172 : alpha=alpha, beta=1.0_dp)
616 : END DO
617 : END IF
618 :
619 : !The anti-symmetric density matrix
620 4744 : DO ispin = 1, nspins
621 :
622 : !The symmetric density matrix
623 : CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect), &
624 2436 : gs_mos(ispin)%mos_occ, 0.0_dp, wfm_rho_orb)
625 : CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), -0.5_dp, gs_mos(ispin)%mos_occ, &
626 2436 : evects(ispin, ivect), 1.0_dp, wfm_rho_orb)
627 :
628 2436 : CALL dbcsr_set(work_hmat_asymm(ispin)%matrix, 0.0_dp)
629 4744 : IF (do_admm) THEN
630 : CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, &
631 1300 : wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb)
632 : CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
633 1300 : 0.0_dp, admm_env%work_aux_aux)
634 1300 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao_asymm(ispin)%matrix, keep_sparsity=.TRUE.)
635 : ELSE
636 1136 : CALL copy_fm_to_dbcsr(wfm_rho_orb, work_rho_ia_ao_asymm(ispin)%matrix, keep_sparsity=.TRUE.)
637 : END IF
638 : END DO
639 :
640 2308 : CALL tddft_hfx_matrix(work_hmat_asymm, work_rho_ia_ao_asymm, qs_env)
641 :
642 3426 : IF (do_admm) THEN
643 2572 : DO ispin = 1, nspins
644 : CALL cp_dbcsr_sm_fm_multiply(work_hmat_asymm(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
645 1300 : ncol=nao, alpha=1.0_dp, beta=0.0_dp)
646 :
647 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
648 1300 : admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb)
649 :
650 : CALL parallel_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, &
651 2572 : gs_mos(ispin)%mos_occ, 1.0_dp, Aop_evects(ispin, ivect))
652 : END DO
653 : ELSE
654 2172 : DO ispin = 1, nspins
655 : CALL cp_dbcsr_sm_fm_multiply(work_hmat_asymm(ispin)%matrix, gs_mos(ispin)%mos_occ, &
656 : Aop_evects(ispin, ivect), ncol=nactive(ispin), &
657 2172 : alpha=alpha, beta=1.0_dp)
658 : END DO
659 : END IF
660 : END DO
661 : END IF
662 :
663 1118 : CALL timestop(handle)
664 :
665 1118 : END SUBROUTINE tddfpt_apply_hfx
666 :
667 : ! **************************************************************************************************
668 : !> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
669 : !> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
670 : !> \param evects trial vectors
671 : !> \param gs_mos molecular orbitals optimised for the ground state (only occupied
672 : !> molecular orbitals [component %mos_occ] are needed)
673 : !> \param qs_env Quickstep environment
674 : !> \param admm_env ...
675 : !> \param hfx_section ...
676 : !> \param x_data ...
677 : !> \param symmetry ...
678 : !> \param recalc_integrals ...
679 : !> \param work_rho_ia_ao ...
680 : !> \param work_hmat ...
681 : !> \param wfm_rho_orb ...
682 : ! **************************************************************************************************
683 44 : SUBROUTINE tddfpt_apply_hfxsr_kernel(Aop_evects, evects, gs_mos, qs_env, admm_env, &
684 : hfx_section, x_data, symmetry, recalc_integrals, &
685 44 : work_rho_ia_ao, work_hmat, wfm_rho_orb)
686 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: Aop_evects, evects
687 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
688 : INTENT(in) :: gs_mos
689 : TYPE(qs_environment_type), POINTER :: qs_env
690 : TYPE(admm_type), POINTER :: admm_env
691 : TYPE(section_vals_type), POINTER :: hfx_section
692 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
693 : INTEGER, INTENT(IN) :: symmetry
694 : LOGICAL, INTENT(IN) :: recalc_integrals
695 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: work_rho_ia_ao
696 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
697 : TARGET :: work_hmat
698 : TYPE(cp_fm_type), INTENT(IN) :: wfm_rho_orb
699 :
700 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_hfxsr_kernel'
701 :
702 : INTEGER :: handle, ispin, ivect, nao, nao_aux, &
703 : nspins, nvects
704 : INTEGER, DIMENSION(maxspins) :: nactive
705 : LOGICAL :: reint
706 : REAL(kind=dp) :: alpha
707 :
708 44 : CALL timeset(routineN, handle)
709 :
710 44 : nspins = SIZE(evects, 1)
711 44 : nvects = SIZE(evects, 2)
712 :
713 44 : alpha = 2.0_dp
714 44 : IF (nspins > 1) alpha = 1.0_dp
715 :
716 44 : CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
717 44 : CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux)
718 88 : DO ispin = 1, nspins
719 88 : CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
720 : END DO
721 :
722 44 : reint = recalc_integrals
723 :
724 132 : DO ivect = 1, nvects
725 176 : DO ispin = 1, nspins
726 : CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect), &
727 88 : gs_mos(ispin)%mos_occ, 0.0_dp, wfm_rho_orb)
728 : CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp*symmetry, gs_mos(ispin)%mos_occ, &
729 88 : evects(ispin, ivect), 1.0_dp, wfm_rho_orb)
730 88 : CALL dbcsr_set(work_hmat(ispin)%matrix, 0.0_dp)
731 : CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, &
732 88 : wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb)
733 : CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
734 88 : 0.0_dp, admm_env%work_aux_aux)
735 176 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao(ispin)%matrix, keep_sparsity=.TRUE.)
736 : END DO
737 :
738 88 : CALL tddft_hfx_matrix(work_hmat, work_rho_ia_ao, qs_env, .FALSE., reint, hfx_section, x_data)
739 88 : reint = .FALSE.
740 :
741 220 : DO ispin = 1, nspins
742 : CALL cp_dbcsr_sm_fm_multiply(work_hmat(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
743 88 : ncol=nao, alpha=1.0_dp, beta=0.0_dp)
744 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
745 88 : admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb)
746 : CALL parallel_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, &
747 176 : gs_mos(ispin)%mos_occ, 1.0_dp, Aop_evects(ispin, ivect))
748 : END DO
749 : END DO
750 :
751 44 : CALL timestop(handle)
752 :
753 44 : END SUBROUTINE tddfpt_apply_hfxsr_kernel
754 :
755 : ! **************************************************************************************************
756 : !> \brief ...Calculate the HFXLR kernel contribution by contracting the Lowdin MO coefficients --
757 : !> transition charges with the exchange-type integrals using the sTDA approximation
758 : !> \param qs_env ...
759 : !> \param sub_env ...
760 : !> \param rcut ...
761 : !> \param hfx_scale ...
762 : !> \param work ...
763 : !> \param X ...
764 : !> \param res ... vector AX with A being the sTDA matrix and X the Davidson trial vector of the
765 : !> eigenvalue problem A*X = omega*X
766 : ! **************************************************************************************************
767 72 : SUBROUTINE tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, hfx_scale, work, X, res)
768 :
769 : TYPE(qs_environment_type), POINTER :: qs_env
770 : TYPE(tddfpt_subgroup_env_type) :: sub_env
771 : REAL(KIND=dp), INTENT(IN) :: rcut, hfx_scale
772 : TYPE(tddfpt_work_matrices) :: work
773 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: X
774 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: res
775 :
776 : CHARACTER(len=*), PARAMETER :: routineN = 'tddfpt_apply_hfxlr_kernel'
777 :
778 : INTEGER :: handle, iatom, ispin, jatom, natom, &
779 : nsgf, nspins
780 : INTEGER, DIMENSION(2) :: nactive
781 : REAL(KIND=dp) :: dr, eps_filter, fcut, gabr
782 : REAL(KIND=dp), DIMENSION(3) :: rij
783 72 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock
784 : TYPE(cell_type), POINTER :: cell
785 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
786 : TYPE(cp_fm_type) :: cvec
787 72 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: xtransformed
788 : TYPE(cp_fm_type), POINTER :: ct
789 : TYPE(dbcsr_iterator_type) :: iter
790 : TYPE(dbcsr_type) :: pdens
791 : TYPE(dbcsr_type), POINTER :: tempmat
792 : TYPE(mp_para_env_type), POINTER :: para_env
793 72 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
794 :
795 72 : CALL timeset(routineN, handle)
796 :
797 : ! parameters
798 72 : eps_filter = 1.E-08_dp
799 :
800 72 : nspins = SIZE(X)
801 144 : DO ispin = 1, nspins
802 144 : CALL cp_fm_get_info(X(ispin), ncol_global=nactive(ispin))
803 : END DO
804 :
805 72 : para_env => sub_env%para_env
806 :
807 72 : CALL get_qs_env(qs_env, natom=natom, cell=cell, particle_set=particle_set)
808 :
809 : ! calculate Loewdin transformed Davidson trial vector tilde(X)=S^1/2*X
810 : ! and tilde(tilde(X))=S^1/2_A*tilde(X)_A
811 288 : ALLOCATE (xtransformed(nspins))
812 144 : DO ispin = 1, nspins
813 72 : NULLIFY (fmstruct)
814 72 : ct => work%ctransformed(ispin)
815 72 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct)
816 144 : CALL cp_fm_create(matrix=xtransformed(ispin), matrix_struct=fmstruct, name="XTRANSFORMED")
817 : END DO
818 72 : CALL get_lowdin_x(work%shalf, X, xtransformed)
819 :
820 144 : DO ispin = 1, nspins
821 72 : ct => work%ctransformed(ispin)
822 72 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct, nrow_global=nsgf)
823 72 : CALL cp_fm_create(cvec, fmstruct)
824 : !
825 72 : tempmat => work%shalf
826 72 : CALL dbcsr_create(pdens, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
827 : ! P(nu,mu) = SUM_j XT(nu,j)*CT(mu,j)
828 72 : ct => work%ctransformed(ispin)
829 72 : CALL dbcsr_set(pdens, 0.0_dp)
830 : CALL cp_dbcsr_plus_fm_fm_t(pdens, xtransformed(ispin), ct, nactive(ispin), &
831 72 : 1.0_dp, keep_sparsity=.FALSE.)
832 72 : CALL dbcsr_filter(pdens, eps_filter)
833 : ! Apply PP*gab -> PP; gab = gamma_coulomb
834 : ! P(nu,mu) = P(nu,mu)*g(a of nu,b of mu)
835 72 : CALL dbcsr_iterator_start(iter, pdens)
836 396 : DO WHILE (dbcsr_iterator_blocks_left(iter))
837 324 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock)
838 1296 : rij = particle_set(iatom)%r - particle_set(jatom)%r
839 1296 : rij = pbc(rij, cell)
840 1296 : dr = SQRT(SUM(rij(:)**2))
841 324 : gabr = 1._dp/rcut
842 324 : IF (dr < 1.e-6) THEN
843 108 : gabr = 2._dp*gabr/SQRT(3.1415926_dp)
844 : ELSE
845 216 : gabr = ERF(gabr*dr)/dr
846 : fcut = EXP(dr - 4._dp*rcut)
847 216 : fcut = fcut/(fcut + 1._dp)
848 : END IF
849 21924 : pblock = hfx_scale*gabr*pblock
850 : END DO
851 72 : CALL dbcsr_iterator_stop(iter)
852 : ! CV(mu,i) = P(nu,mu)*CT(mu,i)
853 72 : CALL cp_dbcsr_sm_fm_multiply(pdens, ct, cvec, nactive(ispin), 1.0_dp, 0.0_dp)
854 : ! rho(nu,i) = rho(nu,i) + ShalfP(nu,mu)*CV(mu,i)
855 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), &
856 72 : -1.0_dp, 1.0_dp)
857 : !
858 72 : CALL dbcsr_release(pdens)
859 : !
860 288 : CALL cp_fm_release(cvec)
861 : END DO
862 :
863 72 : CALL cp_fm_release(xtransformed)
864 :
865 72 : CALL timestop(handle)
866 :
867 144 : END SUBROUTINE tddfpt_apply_hfxlr_kernel
868 :
869 : ! **************************************************************************************************
870 :
871 : END MODULE qs_tddfpt2_operators
|
