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