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