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 : USE input_section_types, ONLY: section_vals_get,&
34 : section_vals_get_subs_vals,&
35 : section_vals_type
36 : USE kinds, ONLY: dp
37 : USE message_passing, ONLY: mp_para_env_type
38 : USE parallel_gemm_api, ONLY: parallel_gemm
39 : USE particle_types, ONLY: particle_type
40 : USE pw_env_types, ONLY: pw_env_get,&
41 : pw_env_type
42 : USE pw_methods, ONLY: pw_axpy,&
43 : pw_multiply,&
44 : pw_scale,&
45 : pw_transfer
46 : USE pw_poisson_methods, ONLY: pw_poisson_solve
47 : USE pw_poisson_types, ONLY: pw_poisson_type
48 : USE pw_pool_types, ONLY: pw_pool_p_type
49 : USE pw_types, ONLY: pw_c1d_gs_type,&
50 : pw_r3d_rs_type
51 : USE qs_environment_types, ONLY: get_qs_env,&
52 : qs_environment_type
53 : USE qs_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 realspace_grid_types, ONLY: realspace_grid_desc_p_type,&
62 : realspace_grid_type
63 : #include "./base/base_uses.f90"
64 :
65 : IMPLICIT NONE
66 :
67 : PRIVATE
68 :
69 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_operators'
70 :
71 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
72 : ! number of first derivative components (3: d/dx, d/dy, d/dz)
73 : INTEGER, PARAMETER, PRIVATE :: nderivs = 3
74 : INTEGER, PARAMETER, PRIVATE :: maxspins = 2
75 :
76 : PUBLIC :: tddfpt_apply_energy_diff, tddfpt_apply_coulomb, tddfpt_apply_hfx, &
77 : tddfpt_apply_xc_potential, tddfpt_apply_hfxlr_kernel, tddfpt_apply_hfxsr_kernel
78 :
79 : ! **************************************************************************************************
80 :
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief Apply orbital energy difference term:
85 : !> Aop_evects(spin,state) += KS(spin) * evects(spin,state) -
86 : !> S * evects(spin,state) * diag(evals_occ(spin))
87 : !> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
88 : !> \param evects trial vectors C_{1,i}
89 : !> \param S_evects S * C_{1,i}
90 : !> \param gs_mos molecular orbitals optimised for the ground state (only occupied orbital
91 : !> energies [component %evals_occ] are needed)
92 : !> \param matrix_ks Kohn-Sham matrix
93 : !> \param tddfpt_control ...
94 : !> \par History
95 : !> * 05.2016 initialise all matrix elements in one go [Sergey Chulkov]
96 : !> * 03.2017 renamed from tddfpt_init_energy_diff(), altered prototype [Sergey Chulkov]
97 : !> \note Based on the subroutine p_op_l1() which was originally created by
98 : !> Thomas Chassaing on 08.2002.
99 : ! **************************************************************************************************
100 6814 : SUBROUTINE tddfpt_apply_energy_diff(Aop_evects, evects, S_evects, gs_mos, matrix_ks, tddfpt_control)
101 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects
102 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: 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 : TYPE(tddfpt2_control_type), INTENT(in), POINTER :: tddfpt_control
107 :
108 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_energy_diff'
109 :
110 : INTEGER :: handle, i, ispin, ivect, j, nactive, &
111 : nao, nspins, nvects, spin2
112 6814 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals_active
113 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
114 : TYPE(cp_fm_type) :: hevec
115 :
116 6814 : CALL timeset(routineN, handle)
117 :
118 6814 : nspins = SIZE(evects, 1)
119 6814 : nvects = SIZE(evects, 2)
120 :
121 14644 : DO ispin = 1, SIZE(evects, 1)
122 : CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
123 7830 : nrow_global=nao, ncol_global=nactive)
124 7830 : CALL cp_fm_create(hevec, matrix_struct)
125 23490 : ALLOCATE (evals_active(nactive))
126 87326 : DO i = 1, nactive
127 79496 : j = gs_mos(ispin)%index_active(i)
128 87326 : evals_active(i) = gs_mos(ispin)%evals_occ(j)
129 : END DO
130 :
131 7830 : IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
132 : spin2 = ispin
133 : ELSE
134 96 : spin2 = 2
135 : END IF
136 :
137 28308 : DO ivect = 1, nvects
138 : CALL cp_dbcsr_sm_fm_multiply(matrix_ks(spin2)%matrix, evects(ispin, ivect), &
139 : Aop_evects(ispin, ivect), ncol=nactive, &
140 20478 : alpha=1.0_dp, beta=1.0_dp)
141 :
142 20478 : IF (ASSOCIATED(gs_mos(ispin)%evals_occ_matrix)) THEN
143 : ! orbital energy correction: evals_occ_matrix is not a diagonal matrix
144 : CALL parallel_gemm('N', 'N', nao, nactive, nactive, 1.0_dp, &
145 : S_evects(ispin, ivect), gs_mos(ispin)%evals_occ_matrix, &
146 762 : 0.0_dp, hevec)
147 : ELSE
148 19716 : CALL cp_fm_to_fm(S_evects(ispin, ivect), hevec)
149 19716 : CALL cp_fm_column_scale(hevec, evals_active)
150 : END IF
151 :
152 : ! KS * C1 - S * C1 * occupied_orbital_energies
153 28308 : CALL cp_fm_scale_and_add(1.0_dp, Aop_evects(ispin, ivect), -1.0_dp, hevec)
154 : END DO
155 7830 : DEALLOCATE (evals_active)
156 22474 : CALL cp_fm_release(hevec)
157 : END DO
158 :
159 6814 : CALL timestop(handle)
160 :
161 13628 : END SUBROUTINE tddfpt_apply_energy_diff
162 :
163 : ! **************************************************************************************************
164 : !> \brief Update v_rspace by adding coulomb term.
165 : !> \param A_ia_rspace action of TDDFPT operator on the trial vector expressed in a plane wave
166 : !> representation (modified on exit)
167 : !> \param rho_ia_g response density in reciprocal space for the given trial vector
168 : !> \param local_rho_set ...
169 : !> \param hartree_local ...
170 : !> \param qs_env ...
171 : !> \param sub_env the full sub_environment needed for calculation
172 : !> \param gapw Flag indicating GAPW cacluation
173 : !> \param work_v_gspace work reciprocal-space grid to store Coulomb potential (modified on exit)
174 : !> \param work_v_rspace work real-space grid to store Coulomb potential (modified on exit)
175 : !> \param tddfpt_mgrid ...
176 : !> \par History
177 : !> * 05.2016 compute all coulomb terms in one go [Sergey Chulkov]
178 : !> * 03.2017 proceed excited states sequentially; minimise the number of conversions between
179 : !> DBCSR and FM matrices [Sergey Chulkov]
180 : !> * 06.2018 return the action expressed in the plane wave representation instead of the one
181 : !> in the atomic basis set representation
182 : !> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by
183 : !> Mohamed Fawzi on 10.2002.
184 : ! **************************************************************************************************
185 8078 : SUBROUTINE tddfpt_apply_coulomb(A_ia_rspace, rho_ia_g, local_rho_set, hartree_local, &
186 : qs_env, sub_env, gapw, work_v_gspace, work_v_rspace, tddfpt_mgrid)
187 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT) :: A_ia_rspace
188 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_ia_g
189 : TYPE(local_rho_type), POINTER :: local_rho_set
190 : TYPE(hartree_local_type), POINTER :: hartree_local
191 : TYPE(qs_environment_type), POINTER :: qs_env
192 : TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
193 : LOGICAL, INTENT(IN) :: gapw
194 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: work_v_gspace
195 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: work_v_rspace
196 : LOGICAL, INTENT(IN) :: tddfpt_mgrid
197 :
198 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_coulomb'
199 :
200 : INTEGER :: handle, ispin, nspins
201 : REAL(kind=dp) :: alpha, pair_energy
202 : TYPE(pw_env_type), POINTER :: pw_env
203 : TYPE(pw_poisson_type), POINTER :: poisson_env
204 8078 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: my_pools
205 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
206 8078 : POINTER :: my_rs_descs
207 8078 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: my_rs_grids
208 :
209 8078 : CALL timeset(routineN, handle)
210 :
211 8078 : nspins = SIZE(A_ia_rspace)
212 8078 : pw_env => sub_env%pw_env
213 8078 : IF (tddfpt_mgrid) THEN
214 : CALL pw_env_get(pw_env, poisson_env=poisson_env, rs_grids=my_rs_grids, &
215 86 : rs_descs=my_rs_descs, pw_pools=my_pools)
216 : ELSE
217 7992 : CALL pw_env_get(pw_env, poisson_env=poisson_env)
218 : END IF
219 :
220 8078 : IF (nspins > 1) THEN
221 1822 : alpha = 1.0_dp
222 : ELSE
223 : ! spin-restricted case: alpha == 2 due to singlet state.
224 : ! In case of triplet states alpha == 0, so we should not call this subroutine at all.
225 6256 : alpha = 2.0_dp
226 : END IF
227 :
228 8078 : IF (gapw) THEN
229 1922 : CPASSERT(ASSOCIATED(local_rho_set))
230 1922 : CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rho_ia_g)
231 1922 : IF (ASSOCIATED(local_rho_set%rho0_mpole%rhoz_cneo_s_gs)) THEN
232 0 : CALL pw_axpy(local_rho_set%rho0_mpole%rhoz_cneo_s_gs, rho_ia_g)
233 : END IF
234 : END IF
235 :
236 8078 : CALL pw_poisson_solve(poisson_env, rho_ia_g, pair_energy, work_v_gspace)
237 8078 : CALL pw_transfer(work_v_gspace, work_v_rspace)
238 :
239 : ! (i a || j b) = ( i_alpha a_alpha + i_beta a_beta || j_alpha b_alpha + j_beta b_beta) =
240 : ! tr (Cj_alpha^T * [J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu] * Cb_alpha) +
241 : ! tr (Cj_beta^T * [J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu] * Cb_beta)
242 17978 : DO ispin = 1, nspins
243 17978 : CALL pw_axpy(work_v_rspace, A_ia_rspace(ispin), alpha)
244 : END DO
245 :
246 8078 : IF (gapw) THEN
247 : CALL Vh_1c_gg_integrals(qs_env, pair_energy, &
248 : hartree_local%ecoul_1c, &
249 : local_rho_set, &
250 1922 : sub_env%para_env, tddft=.TRUE., core_2nd=.TRUE.)
251 1922 : CALL pw_scale(work_v_rspace, work_v_rspace%pw_grid%dvol)
252 1922 : IF (tddfpt_mgrid) THEN
253 : CALL integrate_vhg0_rspace(qs_env, work_v_rspace, sub_env%para_env, &
254 : calculate_forces=.FALSE., &
255 : local_rho_set=local_rho_set, my_pools=my_pools, &
256 50 : my_rs_descs=my_rs_descs)
257 : ELSE
258 : CALL integrate_vhg0_rspace(qs_env, work_v_rspace, sub_env%para_env, &
259 : calculate_forces=.FALSE., &
260 1872 : local_rho_set=local_rho_set)
261 : END IF
262 : END IF
263 :
264 8078 : CALL timestop(handle)
265 :
266 8078 : END SUBROUTINE tddfpt_apply_coulomb
267 :
268 : ! **************************************************************************************************
269 : !> \brief Routine for applying fxc potential
270 : !> \param A_ia_rspace action of TDDFPT operator on trial vectors expressed in a plane wave
271 : !> representation (modified on exit)
272 : !> \param fxc_rspace ...
273 : !> \param rho_ia_struct response density for the given trial vector
274 : !> \param is_rks_triplets ...
275 : ! **************************************************************************************************
276 156 : SUBROUTINE tddfpt_apply_xc_potential(A_ia_rspace, fxc_rspace, rho_ia_struct, is_rks_triplets)
277 :
278 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT) :: A_ia_rspace
279 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rspace
280 : TYPE(qs_rho_type), POINTER :: rho_ia_struct
281 : LOGICAL, INTENT(in) :: is_rks_triplets
282 :
283 : INTEGER :: nspins
284 : REAL(KIND=dp) :: alpha
285 156 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r
286 :
287 156 : nspins = SIZE(A_ia_rspace)
288 :
289 156 : alpha = 1.0_dp
290 :
291 156 : CALL qs_rho_get(rho_ia_struct, rho_r=rho1_r)
292 :
293 156 : IF (nspins == 2) THEN
294 0 : CALL pw_multiply(A_ia_rspace(1), fxc_rspace(1), rho1_r(1), alpha)
295 0 : CALL pw_multiply(A_ia_rspace(1), fxc_rspace(2), rho1_r(2), alpha)
296 0 : CALL pw_multiply(A_ia_rspace(2), fxc_rspace(3), rho1_r(2), alpha)
297 0 : CALL pw_multiply(A_ia_rspace(2), fxc_rspace(2), rho1_r(1), alpha)
298 156 : ELSE IF (is_rks_triplets) THEN
299 0 : CALL pw_multiply(A_ia_rspace(1), fxc_rspace(1), rho1_r(1), alpha)
300 0 : CALL pw_multiply(A_ia_rspace(1), fxc_rspace(2), rho1_r(1), -alpha)
301 : ELSE
302 156 : CALL pw_multiply(A_ia_rspace(1), fxc_rspace(1), rho1_r(1), alpha)
303 156 : CALL pw_multiply(A_ia_rspace(1), fxc_rspace(2), rho1_r(1), alpha)
304 : END IF
305 :
306 156 : END SUBROUTINE tddfpt_apply_xc_potential
307 :
308 : ! **************************************************************************************************
309 : !> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
310 : !> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
311 : !> \param evects trial vectors
312 : !> \param gs_mos molecular orbitals optimised for the ground state (only occupied
313 : !> molecular orbitals [component %mos_occ] are needed)
314 : !> \param do_admm perform auxiliary density matrix method calculations
315 : !> \param qs_env Quickstep environment
316 : !> \param work_rho_ia_ao_symm ...
317 : !> \param work_hmat_symm ...
318 : !> \param work_rho_ia_ao_asymm ...
319 : !> \param work_hmat_asymm ...
320 : !> \param wfm_rho_orb ...
321 : !> \par History
322 : !> * 05.2016 compute all exact-exchange terms in one go [Sergey Chulkov]
323 : !> * 03.2017 code related to ADMM correction is now moved to tddfpt_apply_admm_correction()
324 : !> in order to compute this correction within parallel groups [Sergey Chulkov]
325 : !> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by
326 : !> Mohamed Fawzi on 10.2002.
327 : ! **************************************************************************************************
328 1572 : SUBROUTINE tddfpt_apply_hfx(Aop_evects, evects, gs_mos, do_admm, qs_env, &
329 1572 : work_rho_ia_ao_symm, work_hmat_symm, work_rho_ia_ao_asymm, &
330 1572 : work_hmat_asymm, wfm_rho_orb)
331 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: Aop_evects
332 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects
333 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
334 : INTENT(in) :: gs_mos
335 : LOGICAL, INTENT(in) :: do_admm
336 : TYPE(qs_environment_type), POINTER :: qs_env
337 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: work_rho_ia_ao_symm
338 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
339 : TARGET :: work_hmat_symm
340 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: work_rho_ia_ao_asymm
341 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
342 : TARGET :: work_hmat_asymm
343 : TYPE(cp_fm_type), INTENT(IN) :: wfm_rho_orb
344 :
345 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_hfx'
346 :
347 : INTEGER :: handle, ispin, ivect, nao, nao_aux, &
348 : nspins, nvects
349 : INTEGER, DIMENSION(maxspins) :: nactive
350 : LOGICAL :: do_hfx
351 : REAL(kind=dp) :: alpha
352 : TYPE(admm_type), POINTER :: admm_env
353 : TYPE(section_vals_type), POINTER :: hfx_section, input
354 :
355 1572 : CALL timeset(routineN, handle)
356 :
357 : ! Check for hfx section
358 1572 : CALL get_qs_env(qs_env, input=input)
359 1572 : hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
360 1572 : CALL section_vals_get(hfx_section, explicit=do_hfx)
361 :
362 1572 : IF (do_hfx) THEN
363 1572 : nspins = SIZE(evects, 1)
364 1572 : nvects = SIZE(evects, 2)
365 :
366 1572 : IF (SIZE(gs_mos) > 1) THEN
367 98 : alpha = 1.0_dp
368 : ELSE
369 1474 : alpha = 2.0_dp
370 : END IF
371 :
372 1572 : CALL cp_fm_get_info(evects(1, 1), nrow_global=nao)
373 3222 : DO ispin = 1, nspins
374 3222 : CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
375 : END DO
376 :
377 1572 : IF (do_admm) THEN
378 880 : CALL get_qs_env(qs_env, admm_env=admm_env)
379 880 : CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux)
380 : END IF
381 :
382 : !Note: the symmetrized transition density matrix is P = 0.5*(C*evect^T + evect*C^T)
383 : ! in the end, we only want evect*C^T for consistency with the MO formulation of TDDFT
384 : ! therefore, we go in 2 steps: with the symmetric 0.5*(C*evect^T + evect*C^T) and
385 : ! the antisymemtric 0.5*(C*evect^T - evect*C^T)
386 :
387 : ! some stuff from qs_ks_build_kohn_sham_matrix
388 : ! TO DO: add SIC support
389 4552 : DO ivect = 1, nvects
390 6088 : DO ispin = 1, nspins
391 :
392 : !The symmetric density matrix
393 : CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect), &
394 3108 : gs_mos(ispin)%mos_active, 0.0_dp, wfm_rho_orb)
395 : CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, gs_mos(ispin)%mos_active, &
396 3108 : evects(ispin, ivect), 1.0_dp, wfm_rho_orb)
397 :
398 3108 : CALL dbcsr_set(work_hmat_symm(ispin)%matrix, 0.0_dp)
399 6088 : IF (do_admm) THEN
400 : CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, &
401 1622 : wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb)
402 : CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
403 1622 : 0.0_dp, admm_env%work_aux_aux)
404 1622 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao_symm(ispin)%matrix, keep_sparsity=.TRUE.)
405 : ELSE
406 1486 : CALL copy_fm_to_dbcsr(wfm_rho_orb, work_rho_ia_ao_symm(ispin)%matrix, keep_sparsity=.TRUE.)
407 : END IF
408 : END DO
409 :
410 2980 : CALL tddft_hfx_matrix(work_hmat_symm, work_rho_ia_ao_symm, qs_env)
411 :
412 2980 : IF (do_admm) THEN
413 3216 : DO ispin = 1, nspins
414 : CALL cp_dbcsr_sm_fm_multiply(work_hmat_symm(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
415 1622 : ncol=nao, alpha=1.0_dp, beta=0.0_dp)
416 :
417 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
418 1622 : admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb)
419 :
420 : CALL parallel_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, &
421 3216 : gs_mos(ispin)%mos_active, 1.0_dp, Aop_evects(ispin, ivect))
422 : END DO
423 : ELSE
424 2872 : DO ispin = 1, nspins
425 : CALL cp_dbcsr_sm_fm_multiply(work_hmat_symm(ispin)%matrix, gs_mos(ispin)%mos_active, &
426 : Aop_evects(ispin, ivect), ncol=nactive(ispin), &
427 2872 : alpha=alpha, beta=1.0_dp)
428 : END DO
429 : END IF
430 :
431 : !The anti-symmetric density matrix
432 6088 : DO ispin = 1, nspins
433 :
434 : !The symmetric density matrix
435 : CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect), &
436 3108 : gs_mos(ispin)%mos_active, 0.0_dp, wfm_rho_orb)
437 : CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), -0.5_dp, gs_mos(ispin)%mos_active, &
438 3108 : evects(ispin, ivect), 1.0_dp, wfm_rho_orb)
439 :
440 3108 : CALL dbcsr_set(work_hmat_asymm(ispin)%matrix, 0.0_dp)
441 6088 : IF (do_admm) THEN
442 : CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, &
443 1622 : wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb)
444 : CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
445 1622 : 0.0_dp, admm_env%work_aux_aux)
446 1622 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao_asymm(ispin)%matrix, keep_sparsity=.TRUE.)
447 : ELSE
448 1486 : CALL copy_fm_to_dbcsr(wfm_rho_orb, work_rho_ia_ao_asymm(ispin)%matrix, keep_sparsity=.TRUE.)
449 : END IF
450 : END DO
451 :
452 2980 : CALL tddft_hfx_matrix(work_hmat_asymm, work_rho_ia_ao_asymm, qs_env)
453 :
454 4552 : IF (do_admm) THEN
455 3216 : DO ispin = 1, nspins
456 : CALL cp_dbcsr_sm_fm_multiply(work_hmat_asymm(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
457 1622 : ncol=nao, alpha=1.0_dp, beta=0.0_dp)
458 :
459 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
460 1622 : admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb)
461 :
462 : CALL parallel_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, &
463 3216 : gs_mos(ispin)%mos_active, 1.0_dp, Aop_evects(ispin, ivect))
464 : END DO
465 : ELSE
466 2872 : DO ispin = 1, nspins
467 : CALL cp_dbcsr_sm_fm_multiply(work_hmat_asymm(ispin)%matrix, gs_mos(ispin)%mos_active, &
468 : Aop_evects(ispin, ivect), ncol=nactive(ispin), &
469 2872 : alpha=alpha, beta=1.0_dp)
470 : END DO
471 : END IF
472 : END DO
473 : END IF
474 :
475 1572 : CALL timestop(handle)
476 :
477 1572 : END SUBROUTINE tddfpt_apply_hfx
478 :
479 : ! **************************************************************************************************
480 : !> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
481 : !> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
482 : !> \param evects trial vectors
483 : !> \param gs_mos molecular orbitals optimised for the ground state (only occupied
484 : !> molecular orbitals [component %mos_occ] are needed)
485 : !> \param qs_env Quickstep environment
486 : !> \param admm_env ...
487 : !> \param hfx_section ...
488 : !> \param x_data ...
489 : !> \param symmetry ...
490 : !> \param recalc_integrals ...
491 : !> \param work_rho_ia_ao ...
492 : !> \param work_hmat ...
493 : !> \param wfm_rho_orb ...
494 : ! **************************************************************************************************
495 44 : SUBROUTINE tddfpt_apply_hfxsr_kernel(Aop_evects, evects, gs_mos, qs_env, admm_env, &
496 : hfx_section, x_data, symmetry, recalc_integrals, &
497 44 : work_rho_ia_ao, work_hmat, wfm_rho_orb)
498 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: Aop_evects, evects
499 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
500 : INTENT(in) :: gs_mos
501 : TYPE(qs_environment_type), POINTER :: qs_env
502 : TYPE(admm_type), POINTER :: admm_env
503 : TYPE(section_vals_type), POINTER :: hfx_section
504 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
505 : INTEGER, INTENT(IN) :: symmetry
506 : LOGICAL, INTENT(IN) :: recalc_integrals
507 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: work_rho_ia_ao
508 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
509 : TARGET :: work_hmat
510 : TYPE(cp_fm_type), INTENT(IN) :: wfm_rho_orb
511 :
512 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_apply_hfxsr_kernel'
513 :
514 : INTEGER :: handle, ispin, ivect, nao, nao_aux, &
515 : nspins, nvects
516 : INTEGER, DIMENSION(maxspins) :: nactive
517 : LOGICAL :: reint
518 : REAL(kind=dp) :: alpha
519 :
520 44 : CALL timeset(routineN, handle)
521 :
522 44 : nspins = SIZE(evects, 1)
523 44 : nvects = SIZE(evects, 2)
524 :
525 44 : alpha = 2.0_dp
526 44 : IF (nspins > 1) alpha = 1.0_dp
527 :
528 44 : CALL cp_fm_get_info(evects(1, 1), nrow_global=nao)
529 44 : CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux)
530 88 : DO ispin = 1, nspins
531 88 : CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
532 : END DO
533 :
534 44 : reint = recalc_integrals
535 :
536 132 : DO ivect = 1, nvects
537 176 : DO ispin = 1, nspins
538 : CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect), &
539 88 : gs_mos(ispin)%mos_active, 0.0_dp, wfm_rho_orb)
540 : CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp*symmetry, gs_mos(ispin)%mos_active, &
541 88 : evects(ispin, ivect), 1.0_dp, wfm_rho_orb)
542 88 : CALL dbcsr_set(work_hmat(ispin)%matrix, 0.0_dp)
543 : CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, &
544 88 : wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb)
545 : CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
546 88 : 0.0_dp, admm_env%work_aux_aux)
547 176 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao(ispin)%matrix, keep_sparsity=.TRUE.)
548 : END DO
549 :
550 88 : CALL tddft_hfx_matrix(work_hmat, work_rho_ia_ao, qs_env, .FALSE., reint, hfx_section, x_data)
551 88 : reint = .FALSE.
552 :
553 220 : DO ispin = 1, nspins
554 : CALL cp_dbcsr_sm_fm_multiply(work_hmat(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
555 88 : ncol=nao, alpha=1.0_dp, beta=0.0_dp)
556 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
557 88 : admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb)
558 : CALL parallel_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, &
559 176 : gs_mos(ispin)%mos_active, 1.0_dp, Aop_evects(ispin, ivect))
560 : END DO
561 : END DO
562 :
563 44 : CALL timestop(handle)
564 :
565 44 : END SUBROUTINE tddfpt_apply_hfxsr_kernel
566 :
567 : ! **************************************************************************************************
568 : !> \brief ...Calculate the HFXLR kernel contribution by contracting the Lowdin MO coefficients --
569 : !> transition charges with the exchange-type integrals using the sTDA approximation
570 : !> \param qs_env ...
571 : !> \param sub_env ...
572 : !> \param rcut ...
573 : !> \param hfx_scale ...
574 : !> \param work ...
575 : !> \param X ...
576 : !> \param res ... vector AX with A being the sTDA matrix and X the Davidson trial vector of the
577 : !> eigenvalue problem A*X = omega*X
578 : ! **************************************************************************************************
579 72 : SUBROUTINE tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, hfx_scale, work, X, res)
580 :
581 : TYPE(qs_environment_type), POINTER :: qs_env
582 : TYPE(tddfpt_subgroup_env_type) :: sub_env
583 : REAL(KIND=dp), INTENT(IN) :: rcut, hfx_scale
584 : TYPE(tddfpt_work_matrices) :: work
585 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: X
586 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: res
587 :
588 : CHARACTER(len=*), PARAMETER :: routineN = 'tddfpt_apply_hfxlr_kernel'
589 :
590 : INTEGER :: handle, iatom, ispin, jatom, natom, &
591 : nsgf, nspins
592 : INTEGER, DIMENSION(2) :: nactive
593 : REAL(KIND=dp) :: dr, eps_filter, fcut, gabr
594 : REAL(KIND=dp), DIMENSION(3) :: rij
595 72 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pblock
596 : TYPE(cell_type), POINTER :: cell
597 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
598 : TYPE(cp_fm_type) :: cvec
599 72 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: xtransformed
600 : TYPE(cp_fm_type), POINTER :: ct
601 : TYPE(dbcsr_iterator_type) :: iter
602 : TYPE(dbcsr_type) :: pdens
603 : TYPE(dbcsr_type), POINTER :: tempmat
604 : TYPE(mp_para_env_type), POINTER :: para_env
605 72 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
606 :
607 72 : CALL timeset(routineN, handle)
608 :
609 : ! parameters
610 72 : eps_filter = 1.E-08_dp
611 :
612 72 : nspins = SIZE(X)
613 144 : DO ispin = 1, nspins
614 144 : CALL cp_fm_get_info(X(ispin), ncol_global=nactive(ispin))
615 : END DO
616 :
617 72 : para_env => sub_env%para_env
618 :
619 72 : CALL get_qs_env(qs_env, natom=natom, cell=cell, particle_set=particle_set)
620 :
621 : ! calculate Loewdin transformed Davidson trial vector tilde(X)=S^1/2*X
622 : ! and tilde(tilde(X))=S^1/2_A*tilde(X)_A
623 288 : ALLOCATE (xtransformed(nspins))
624 144 : DO ispin = 1, nspins
625 72 : NULLIFY (fmstruct)
626 72 : ct => work%ctransformed(ispin)
627 72 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct)
628 144 : CALL cp_fm_create(matrix=xtransformed(ispin), matrix_struct=fmstruct, name="XTRANSFORMED")
629 : END DO
630 72 : CALL get_lowdin_x(work%shalf, X, xtransformed)
631 :
632 144 : DO ispin = 1, nspins
633 72 : ct => work%ctransformed(ispin)
634 72 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct, nrow_global=nsgf)
635 72 : CALL cp_fm_create(cvec, fmstruct)
636 : !
637 72 : tempmat => work%shalf
638 72 : CALL dbcsr_create(pdens, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
639 : ! P(nu,mu) = SUM_j XT(nu,j)*CT(mu,j)
640 72 : ct => work%ctransformed(ispin)
641 72 : CALL dbcsr_set(pdens, 0.0_dp)
642 : CALL cp_dbcsr_plus_fm_fm_t(pdens, xtransformed(ispin), ct, nactive(ispin), &
643 72 : 1.0_dp, keep_sparsity=.FALSE.)
644 72 : CALL dbcsr_filter(pdens, eps_filter)
645 : ! Apply PP*gab -> PP; gab = gamma_coulomb
646 : ! P(nu,mu) = P(nu,mu)*g(a of nu,b of mu)
647 72 : CALL dbcsr_iterator_start(iter, pdens)
648 396 : DO WHILE (dbcsr_iterator_blocks_left(iter))
649 324 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock)
650 1296 : rij = particle_set(iatom)%r - particle_set(jatom)%r
651 1296 : rij = pbc(rij, cell)
652 1296 : dr = SQRT(SUM(rij(:)**2))
653 324 : gabr = 1._dp/rcut
654 324 : IF (dr < 1.e-6) THEN
655 108 : gabr = 2._dp*gabr/SQRT(3.1415926_dp)
656 : ELSE
657 216 : gabr = ERF(gabr*dr)/dr
658 : fcut = EXP(dr - 4._dp*rcut)
659 216 : fcut = fcut/(fcut + 1._dp)
660 : END IF
661 21924 : pblock = hfx_scale*gabr*pblock
662 : END DO
663 72 : CALL dbcsr_iterator_stop(iter)
664 : ! CV(mu,i) = P(nu,mu)*CT(mu,i)
665 72 : CALL cp_dbcsr_sm_fm_multiply(pdens, ct, cvec, nactive(ispin), 1.0_dp, 0.0_dp)
666 : ! rho(nu,i) = rho(nu,i) + ShalfP(nu,mu)*CV(mu,i)
667 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), &
668 72 : -1.0_dp, 1.0_dp)
669 : !
670 72 : CALL dbcsr_release(pdens)
671 : !
672 288 : CALL cp_fm_release(cvec)
673 : END DO
674 :
675 72 : CALL cp_fm_release(xtransformed)
676 :
677 72 : CALL timestop(handle)
678 :
679 144 : END SUBROUTINE tddfpt_apply_hfxlr_kernel
680 :
681 : ! **************************************************************************************************
682 :
683 : END MODULE qs_tddfpt2_operators
|
