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