(git:ed6f26b)
Loading...
Searching...
No Matches
qs_tddfpt2_operators.F
Go to the documentation of this file.
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
9 USE admm_types, ONLY: admm_type
10 USE cell_types, ONLY: cell_type,&
11 pbc
12 USE cp_dbcsr_api, ONLY: &
15 dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
22 USE cp_fm_types, ONLY: cp_fm_create,&
30 USE hfx_types, ONLY: hfx_type
34 USE kinds, ONLY: dp
38 USE pw_env_types, ONLY: pw_env_get,&
40 USE pw_methods, ONLY: pw_axpy,&
42 pw_scale,&
47 USE pw_pool_types, ONLY: pw_pool_p_type,&
49 USE pw_types, ONLY: pw_c1d_gs_type,&
56 USE qs_rho_types, ONLY: qs_rho_get,&
68#include "./base/base_uses.f90"
69
70 IMPLICIT NONE
71
72 PRIVATE
73
74 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_operators'
75
76 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
77 ! number of first derivative components (3: d/dx, d/dy, d/dz)
78 INTEGER, PARAMETER, PRIVATE :: nderivs = 3
79 INTEGER, PARAMETER, PRIVATE :: maxspins = 2
80
83
84! **************************************************************************************************
85
86CONTAINS
87
88! **************************************************************************************************
89!> \brief Apply orbital energy difference term:
90!> Aop_evects(spin,state) += KS(spin) * evects(spin,state) -
91!> S * evects(spin,state) * diag(evals_occ(spin))
92!> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
93!> \param evects trial vectors C_{1,i}
94!> \param S_evects S * C_{1,i}
95!> \param gs_mos molecular orbitals optimised for the ground state (only occupied orbital
96!> energies [component %evals_occ] are needed)
97!> \param matrix_ks Kohn-Sham matrix
98!> \par History
99!> * 05.2016 initialise all matrix elements in one go [Sergey Chulkov]
100!> * 03.2017 renamed from tddfpt_init_energy_diff(), altered prototype [Sergey Chulkov]
101!> \note Based on the subroutine p_op_l1() which was originally created by
102!> Thomas Chassaing on 08.2002.
103! **************************************************************************************************
104 SUBROUTINE tddfpt_apply_energy_diff(Aop_evects, evects, S_evects, gs_mos, matrix_ks)
105 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: aop_evects
106 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects, s_evects
107 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
108 INTENT(in) :: gs_mos
109 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_ks
110
111 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_apply_energy_diff'
112
113 INTEGER :: handle, ispin, ivect, nactive, nao, &
114 nspins, nvects
115 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
116 TYPE(cp_fm_type) :: hevec
117
118 CALL timeset(routinen, handle)
119
120 nspins = SIZE(evects, 1)
121 nvects = SIZE(evects, 2)
122
123 DO ispin = 1, nspins
124 CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
125 nrow_global=nao, ncol_global=nactive)
126 CALL cp_fm_create(hevec, matrix_struct)
127
128 DO ivect = 1, nvects
129 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, evects(ispin, ivect), &
130 aop_evects(ispin, ivect), ncol=nactive, &
131 alpha=1.0_dp, beta=1.0_dp)
132
133 IF (ASSOCIATED(gs_mos(ispin)%evals_occ_matrix)) THEN
134 ! orbital energy correction: evals_occ_matrix is not a diagonal matrix
135 CALL parallel_gemm('N', 'N', nao, nactive, nactive, 1.0_dp, &
136 s_evects(ispin, ivect), gs_mos(ispin)%evals_occ_matrix, &
137 0.0_dp, hevec)
138 ELSE
139 CALL cp_fm_to_fm(s_evects(ispin, ivect), hevec)
140 CALL cp_fm_column_scale(hevec, gs_mos(ispin)%evals_occ)
141 END IF
142
143 ! KS * C1 - S * C1 * occupied_orbital_energies
144 CALL cp_fm_scale_and_add(1.0_dp, aop_evects(ispin, ivect), -1.0_dp, hevec)
145 END DO
146 CALL cp_fm_release(hevec)
147 END DO
148
149 CALL timestop(handle)
150
151 END SUBROUTINE tddfpt_apply_energy_diff
152
153! **************************************************************************************************
154!> \brief Update v_rspace by adding coulomb term.
155!> \param A_ia_rspace action of TDDFPT operator on the trial vector expressed in a plane wave
156!> representation (modified on exit)
157!> \param rho_ia_g response density in reciprocal space for the given trial vector
158!> \param local_rho_set ...
159!> \param hartree_local ...
160!> \param qs_env ...
161!> \param sub_env the full sub_environment needed for calculation
162!> \param gapw Flag indicating GAPW cacluation
163!> \param work_v_gspace work reciprocal-space grid to store Coulomb potential (modified on exit)
164!> \param work_v_rspace work real-space grid to store Coulomb potential (modified on exit)
165!> \param tddfpt_mgrid ...
166!> \par History
167!> * 05.2016 compute all coulomb terms in one go [Sergey Chulkov]
168!> * 03.2017 proceed excited states sequentially; minimise the number of conversions between
169!> DBCSR and FM matrices [Sergey Chulkov]
170!> * 06.2018 return the action expressed in the plane wave representation instead of the one
171!> in the atomic basis set representation
172!> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by
173!> Mohamed Fawzi on 10.2002.
174! **************************************************************************************************
175 SUBROUTINE tddfpt_apply_coulomb(A_ia_rspace, rho_ia_g, local_rho_set, hartree_local, &
176 qs_env, sub_env, gapw, work_v_gspace, work_v_rspace, tddfpt_mgrid)
177 TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT) :: a_ia_rspace
178 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_ia_g
179 TYPE(local_rho_type), POINTER :: local_rho_set
180 TYPE(hartree_local_type), POINTER :: hartree_local
181 TYPE(qs_environment_type), POINTER :: qs_env
182 TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
183 LOGICAL, INTENT(IN) :: gapw
184 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: work_v_gspace
185 TYPE(pw_r3d_rs_type), INTENT(INOUT) :: work_v_rspace
186 LOGICAL, INTENT(IN) :: tddfpt_mgrid
187
188 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_apply_coulomb'
189
190 INTEGER :: handle, ispin, nspins
191 REAL(kind=dp) :: alpha, pair_energy
192 TYPE(pw_env_type), POINTER :: pw_env
193 TYPE(pw_poisson_type), POINTER :: poisson_env
194 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: my_pools
195 TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
196 POINTER :: my_rs_descs
197 TYPE(realspace_grid_type), DIMENSION(:), POINTER :: my_rs_grids
198
199 CALL timeset(routinen, handle)
200
201 nspins = SIZE(a_ia_rspace)
202 pw_env => sub_env%pw_env
203 IF (tddfpt_mgrid) THEN
204 CALL pw_env_get(pw_env, poisson_env=poisson_env, rs_grids=my_rs_grids, &
205 rs_descs=my_rs_descs, pw_pools=my_pools)
206 ELSE
207 CALL pw_env_get(pw_env, poisson_env=poisson_env)
208 END IF
209
210 IF (nspins > 1) THEN
211 alpha = 1.0_dp
212 ELSE
213 ! spin-restricted case: alpha == 2 due to singlet state.
214 ! In case of triplet states alpha == 0, so we should not call this subroutine at all.
215 alpha = 2.0_dp
216 END IF
217
218 IF (gapw) THEN
219 cpassert(ASSOCIATED(local_rho_set))
220 CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rho_ia_g)
221 END IF
222
223 CALL pw_poisson_solve(poisson_env, rho_ia_g, pair_energy, work_v_gspace)
224 CALL pw_transfer(work_v_gspace, work_v_rspace)
225
226 ! (i a || j b) = ( i_alpha a_alpha + i_beta a_beta || j_alpha b_alpha + j_beta b_beta) =
227 ! tr (Cj_alpha^T * [J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu] * Cb_alpha) +
228 ! tr (Cj_beta^T * [J_i{alpha}a{alpha}_munu + J_i{beta}a{beta}_munu] * Cb_beta)
229 DO ispin = 1, nspins
230 CALL pw_axpy(work_v_rspace, a_ia_rspace(ispin), alpha)
231 END DO
232
233 IF (gapw) THEN
234 CALL vh_1c_gg_integrals(qs_env, pair_energy, &
235 hartree_local%ecoul_1c, &
236 local_rho_set, &
237 sub_env%para_env, tddft=.true., core_2nd=.true.)
238 CALL pw_scale(work_v_rspace, work_v_rspace%pw_grid%dvol)
239 IF (tddfpt_mgrid) THEN
240 CALL integrate_vhg0_rspace(qs_env, work_v_rspace, sub_env%para_env, &
241 calculate_forces=.false., &
242 local_rho_set=local_rho_set, my_pools=my_pools, &
243 my_rs_descs=my_rs_descs)
244 ELSE
245 CALL integrate_vhg0_rspace(qs_env, work_v_rspace, sub_env%para_env, &
246 calculate_forces=.false., &
247 local_rho_set=local_rho_set)
248 END IF
249 END IF
250
251 CALL timestop(handle)
252
253 END SUBROUTINE tddfpt_apply_coulomb
254
255! **************************************************************************************************
256!> \brief Driver routine for applying fxc (analyic vs. finite difference for testing
257!> \param A_ia_rspace action of TDDFPT operator on trial vectors expressed in a plane wave
258!> representation (modified on exit)
259!> \param kernel_env kernel environment
260!> \param rho_ia_struct response density for the given trial vector
261!> \param is_rks_triplets indicates that the triplet excited states calculation using
262!> spin-unpolarised molecular orbitals has been requested
263!> \param pw_env plain wave environment
264!> \param work_v_xc work real-space grid to store the gradient of the exchange-correlation
265!> potential with respect to the response density (modified on exit)
266!> \param work_v_xc_tau ...
267! **************************************************************************************************
268 SUBROUTINE tddfpt_apply_xc(A_ia_rspace, kernel_env, rho_ia_struct, is_rks_triplets, &
269 pw_env, work_v_xc, work_v_xc_tau)
270
271 TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT) :: a_ia_rspace
272 TYPE(full_kernel_env_type), INTENT(IN) :: kernel_env
273 TYPE(qs_rho_type), POINTER :: rho_ia_struct
274 LOGICAL, INTENT(in) :: is_rks_triplets
275 TYPE(pw_env_type), POINTER :: pw_env
276 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: work_v_xc, work_v_xc_tau
277
278 INTEGER :: ispin, nspins
279
280 nspins = SIZE(a_ia_rspace)
281
282 IF (kernel_env%deriv2_analytic) THEN
283 CALL tddfpt_apply_xc_analytic(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
284 pw_env, work_v_xc, work_v_xc_tau)
285 ELSE
286 CALL tddfpt_apply_xc_fd(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
287 pw_env, work_v_xc, work_v_xc_tau)
288 END IF
289
290 DO ispin = 1, nspins
291 ! pw2 = pw2 + alpha * pw1
292 CALL pw_axpy(work_v_xc(ispin), a_ia_rspace(ispin), kernel_env%alpha)
293 END DO
294
295 END SUBROUTINE tddfpt_apply_xc
296
297! **************************************************************************************************
298!> \brief Routine for applying fxc potential
299!> \param A_ia_rspace action of TDDFPT operator on trial vectors expressed in a plane wave
300!> representation (modified on exit)
301!> \param fxc_rspace ...
302!> \param rho_ia_struct response density for the given trial vector
303!> \param is_rks_triplets ...
304! **************************************************************************************************
305 SUBROUTINE tddfpt_apply_xc_potential(A_ia_rspace, fxc_rspace, rho_ia_struct, is_rks_triplets)
306
307 TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(INOUT) :: a_ia_rspace
308 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rspace
309 TYPE(qs_rho_type), POINTER :: rho_ia_struct
310 LOGICAL, INTENT(in) :: is_rks_triplets
311
312 INTEGER :: nspins
313 REAL(kind=dp) :: alpha
314 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r
315
316 nspins = SIZE(a_ia_rspace)
317
318 alpha = 1.0_dp
319
320 CALL qs_rho_get(rho_ia_struct, rho_r=rho1_r)
321
322 IF (nspins == 2) THEN
323 CALL pw_multiply(a_ia_rspace(1), fxc_rspace(1), rho1_r(1), alpha)
324 CALL pw_multiply(a_ia_rspace(1), fxc_rspace(2), rho1_r(2), alpha)
325 CALL pw_multiply(a_ia_rspace(2), fxc_rspace(3), rho1_r(2), alpha)
326 CALL pw_multiply(a_ia_rspace(2), fxc_rspace(2), rho1_r(1), alpha)
327 ELSE IF (is_rks_triplets) THEN
328 CALL pw_multiply(a_ia_rspace(1), fxc_rspace(1), rho1_r(1), alpha)
329 CALL pw_multiply(a_ia_rspace(1), fxc_rspace(2), rho1_r(1), -alpha)
330 ELSE
331 CALL pw_multiply(a_ia_rspace(1), fxc_rspace(1), rho1_r(1), alpha)
332 CALL pw_multiply(a_ia_rspace(1), fxc_rspace(2), rho1_r(1), alpha)
333 END IF
334
335 END SUBROUTINE tddfpt_apply_xc_potential
336
337! **************************************************************************************************
338!> \brief Update A_ia_munu by adding exchange-correlation term.
339!> \param kernel_env kernel environment
340!> \param rho_ia_struct response density for the given trial vector
341!> \param is_rks_triplets indicates that the triplet excited states calculation using
342!> spin-unpolarised molecular orbitals has been requested
343!> \param nspins ...
344!> \param pw_env plain wave environment
345!> \param work_v_xc work real-space grid to store the gradient of the exchange-correlation
346!> potential with respect to the response density (modified on exit)
347!> \param work_v_xc_tau ...
348!> \par History
349!> * 05.2016 compute all kernel terms in one go [Sergey Chulkov]
350!> * 03.2017 proceed excited states sequentially; minimise the number of conversions between
351!> DBCSR and FM matrices [Sergey Chulkov]
352!> * 06.2018 return the action expressed in the plane wave representation instead of the one
353!> in the atomic basis set representation
354!> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by
355!> Mohamed Fawzi on 10.2002.
356! **************************************************************************************************
357 SUBROUTINE tddfpt_apply_xc_analytic(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
358 pw_env, work_v_xc, work_v_xc_tau)
359 TYPE(full_kernel_env_type), INTENT(in) :: kernel_env
360 TYPE(qs_rho_type), POINTER :: rho_ia_struct
361 LOGICAL, INTENT(in) :: is_rks_triplets
362 INTEGER, INTENT(in) :: nspins
363 TYPE(pw_env_type), POINTER :: pw_env
364 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: work_v_xc, work_v_xc_tau
365
366 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_apply_xc_analytic'
367
368 INTEGER :: handle, ispin
369 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_ia_g, rho_ia_g2
370 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
371 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_ia_r, rho_ia_r2, tau_ia_r, tau_ia_r2
372
373 CALL timeset(routinen, handle)
374
375 CALL qs_rho_get(rho_ia_struct, rho_g=rho_ia_g, rho_r=rho_ia_r, tau_r=tau_ia_r)
376 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
377
378 IF (debug_this_module) THEN
379 cpassert(SIZE(rho_ia_g) == nspins)
380 cpassert(SIZE(rho_ia_r) == nspins)
381 cpassert((.NOT. ASSOCIATED(tau_ia_r)) .OR. SIZE(tau_ia_r) == nspins)
382 cpassert((.NOT. is_rks_triplets) .OR. nspins == 1)
383 END IF
384
385 NULLIFY (tau_ia_r2)
386 IF (is_rks_triplets) THEN
387 ALLOCATE (rho_ia_r2(2))
388 ALLOCATE (rho_ia_g2(2))
389 rho_ia_r2(1) = rho_ia_r(1)
390 rho_ia_r2(2) = rho_ia_r(1)
391 rho_ia_g2(1) = rho_ia_g(1)
392 rho_ia_g2(2) = rho_ia_g(1)
393
394 IF (ASSOCIATED(tau_ia_r)) THEN
395 ALLOCATE (tau_ia_r2(2))
396 tau_ia_r2(1) = tau_ia_r(1)
397 tau_ia_r2(2) = tau_ia_r(1)
398 END IF
399 ELSE
400 rho_ia_r2 => rho_ia_r
401 rho_ia_g2 => rho_ia_g
402
403 tau_ia_r2 => tau_ia_r
404 END IF
405
406 DO ispin = 1, nspins
407 CALL pw_zero(work_v_xc(ispin))
408 IF (ASSOCIATED(work_v_xc_tau)) CALL pw_zero(work_v_xc_tau(ispin))
409 END DO
410
411 CALL xc_rho_set_update(rho_set=kernel_env%xc_rho1_set, rho_r=rho_ia_r2, rho_g=rho_ia_g2, tau=tau_ia_r2, &
412 needs=kernel_env%xc_rho1_cflags, xc_deriv_method_id=kernel_env%deriv_method_id, &
413 xc_rho_smooth_id=kernel_env%rho_smooth_id, pw_pool=auxbas_pw_pool)
414
415 CALL xc_calc_2nd_deriv_analytical(v_xc=work_v_xc, v_xc_tau=work_v_xc_tau, deriv_set=kernel_env%xc_deriv_set, &
416 rho_set=kernel_env%xc_rho_set, &
417 rho1_set=kernel_env%xc_rho1_set, pw_pool=auxbas_pw_pool, &
418 xc_section=kernel_env%xc_section, gapw=.false., tddfpt_fac=kernel_env%beta)
419
420 IF (is_rks_triplets) THEN
421 DEALLOCATE (rho_ia_r2)
422 DEALLOCATE (rho_ia_g2)
423 IF (ASSOCIATED(tau_ia_r2)) DEALLOCATE (tau_ia_r2)
424 END IF
425
426 CALL timestop(handle)
427
428 END SUBROUTINE tddfpt_apply_xc_analytic
429
430! **************************************************************************************************
431!> \brief Update A_ia_munu by adding exchange-correlation term using finite difference methods.
432!> \param kernel_env kernel environment
433!> \param rho_ia_struct response density for the given trial vector
434!> \param is_rks_triplets indicates that the triplet excited states calculation using
435!> spin-unpolarised molecular orbitals has been requested
436!> \param nspins ...
437!> \param pw_env plain wave environment
438!> \param work_v_xc work real-space grid to store the gradient of the exchange-correlation
439!> potential with respect to the response density (modified on exit)
440!> \param work_v_xc_tau ...
441! **************************************************************************************************
442 SUBROUTINE tddfpt_apply_xc_fd(kernel_env, rho_ia_struct, is_rks_triplets, nspins, &
443 pw_env, work_v_xc, work_v_xc_tau)
444 TYPE(full_kernel_env_type), INTENT(in) :: kernel_env
445 TYPE(qs_rho_type), POINTER :: rho_ia_struct
446 LOGICAL, INTENT(in) :: is_rks_triplets
447 INTEGER, INTENT(in) :: nspins
448 TYPE(pw_env_type), POINTER :: pw_env
449 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: work_v_xc, work_v_xc_tau
450
451 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_apply_xc_fd'
452
453 INTEGER :: handle, ispin
454 LOGICAL :: lsd, singlet, triplet
455 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g
456 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
457 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, tau1_r
458 TYPE(xc_rho_set_type), POINTER :: rho_set
459
460 CALL timeset(routinen, handle)
461
462 CALL qs_rho_get(rho_ia_struct, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r)
463 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
464 DO ispin = 1, nspins
465 CALL pw_zero(work_v_xc(ispin))
466 END DO
467 rho_set => kernel_env%xc_rho_set
468
469 singlet = .false.
470 triplet = .false.
471 lsd = .false.
472 IF (nspins == 1 .AND. .NOT. is_rks_triplets) THEN
473 singlet = .true.
474 ELSE IF (nspins == 1 .AND. is_rks_triplets) THEN
475 triplet = .true.
476 ELSE IF (nspins == 2) THEN
477 lsd = .true.
478 ELSE
479 cpabort("illegal options")
480 END IF
481
482 IF (ASSOCIATED(tau1_r)) THEN
483 DO ispin = 1, nspins
484 CALL pw_zero(work_v_xc_tau(ispin))
485 END DO
486 END IF
487
488 CALL xc_calc_2nd_deriv_numerical(work_v_xc, work_v_xc_tau, rho_set, rho1_r, rho1_g, tau1_r, &
489 auxbas_pw_pool, kernel_env%xc_section, &
490 is_rks_triplets)
491
492 CALL timestop(handle)
493
494 END SUBROUTINE tddfpt_apply_xc_fd
495
496! **************************************************************************************************
497!> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
498!> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
499!> \param evects trial vectors
500!> \param gs_mos molecular orbitals optimised for the ground state (only occupied
501!> molecular orbitals [component %mos_occ] are needed)
502!> \param do_admm perform auxiliary density matrix method calculations
503!> \param qs_env Quickstep environment
504!> \param work_rho_ia_ao_symm ...
505!> \param work_hmat_symm ...
506!> \param work_rho_ia_ao_asymm ...
507!> \param work_hmat_asymm ...
508!> \param wfm_rho_orb ...
509!> \par History
510!> * 05.2016 compute all exact-exchange terms in one go [Sergey Chulkov]
511!> * 03.2017 code related to ADMM correction is now moved to tddfpt_apply_admm_correction()
512!> in order to compute this correction within parallel groups [Sergey Chulkov]
513!> \note Based on the subroutine kpp1_calc_k_p_p1() which was originally created by
514!> Mohamed Fawzi on 10.2002.
515! **************************************************************************************************
516 SUBROUTINE tddfpt_apply_hfx(Aop_evects, evects, gs_mos, do_admm, qs_env, &
517 work_rho_ia_ao_symm, work_hmat_symm, work_rho_ia_ao_asymm, &
518 work_hmat_asymm, wfm_rho_orb)
519 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: aop_evects
520 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects
521 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
522 INTENT(in) :: gs_mos
523 LOGICAL, INTENT(in) :: do_admm
524 TYPE(qs_environment_type), POINTER :: qs_env
525 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: work_rho_ia_ao_symm
526 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
527 TARGET :: work_hmat_symm
528 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: work_rho_ia_ao_asymm
529 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
530 TARGET :: work_hmat_asymm
531 TYPE(cp_fm_type), INTENT(IN) :: wfm_rho_orb
532
533 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_apply_hfx'
534
535 INTEGER :: handle, ispin, ivect, nao, nao_aux, &
536 nspins, nvects
537 INTEGER, DIMENSION(maxspins) :: nactive
538 LOGICAL :: do_hfx
539 REAL(kind=dp) :: alpha
540 TYPE(admm_type), POINTER :: admm_env
541 TYPE(section_vals_type), POINTER :: hfx_section, input
542
543 CALL timeset(routinen, handle)
544
545 ! Check for hfx section
546 CALL get_qs_env(qs_env, input=input)
547 hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
548 CALL section_vals_get(hfx_section, explicit=do_hfx)
549
550 IF (do_hfx) THEN
551 nspins = SIZE(evects, 1)
552 nvects = SIZE(evects, 2)
553
554 IF (nspins > 1) THEN
555 alpha = 1.0_dp
556 ELSE
557 alpha = 2.0_dp
558 END IF
559
560 CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
561 DO ispin = 1, nspins
562 CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
563 END DO
564
565 IF (do_admm) THEN
566 CALL get_qs_env(qs_env, admm_env=admm_env)
567 CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux)
568 END IF
569
570 !Note: the symmetrized transition density matrix is P = 0.5*(C*evect^T + evect*C^T)
571 ! in the end, we only want evect*C^T for consistency with the MO formulation of TDDFT
572 ! therefore, we go in 2 steps: with the symmetric 0.5*(C*evect^T + evect*C^T) and
573 ! the antisymemtric 0.5*(C*evect^T - evect*C^T)
574
575 ! some stuff from qs_ks_build_kohn_sham_matrix
576 ! TO DO: add SIC support
577 DO ivect = 1, nvects
578 DO ispin = 1, nspins
579
580 !The symmetric density matrix
581 CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect), &
582 gs_mos(ispin)%mos_occ, 0.0_dp, wfm_rho_orb)
583 CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, gs_mos(ispin)%mos_occ, &
584 evects(ispin, ivect), 1.0_dp, wfm_rho_orb)
585
586 CALL dbcsr_set(work_hmat_symm(ispin)%matrix, 0.0_dp)
587 IF (do_admm) THEN
588 CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, &
589 wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb)
590 CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
591 0.0_dp, admm_env%work_aux_aux)
592 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao_symm(ispin)%matrix, keep_sparsity=.true.)
593 ELSE
594 CALL copy_fm_to_dbcsr(wfm_rho_orb, work_rho_ia_ao_symm(ispin)%matrix, keep_sparsity=.true.)
595 END IF
596 END DO
597
598 CALL tddft_hfx_matrix(work_hmat_symm, work_rho_ia_ao_symm, qs_env)
599
600 IF (do_admm) THEN
601 DO ispin = 1, nspins
602 CALL cp_dbcsr_sm_fm_multiply(work_hmat_symm(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
603 ncol=nao, alpha=1.0_dp, beta=0.0_dp)
604
605 CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
606 admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb)
607
608 CALL parallel_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, &
609 gs_mos(ispin)%mos_occ, 1.0_dp, aop_evects(ispin, ivect))
610 END DO
611 ELSE
612 DO ispin = 1, nspins
613 CALL cp_dbcsr_sm_fm_multiply(work_hmat_symm(ispin)%matrix, gs_mos(ispin)%mos_occ, &
614 aop_evects(ispin, ivect), ncol=nactive(ispin), &
615 alpha=alpha, beta=1.0_dp)
616 END DO
617 END IF
618
619 !The anti-symmetric density matrix
620 DO ispin = 1, nspins
621
622 !The symmetric density matrix
623 CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect), &
624 gs_mos(ispin)%mos_occ, 0.0_dp, wfm_rho_orb)
625 CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), -0.5_dp, gs_mos(ispin)%mos_occ, &
626 evects(ispin, ivect), 1.0_dp, wfm_rho_orb)
627
628 CALL dbcsr_set(work_hmat_asymm(ispin)%matrix, 0.0_dp)
629 IF (do_admm) THEN
630 CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, &
631 wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb)
632 CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
633 0.0_dp, admm_env%work_aux_aux)
634 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao_asymm(ispin)%matrix, keep_sparsity=.true.)
635 ELSE
636 CALL copy_fm_to_dbcsr(wfm_rho_orb, work_rho_ia_ao_asymm(ispin)%matrix, keep_sparsity=.true.)
637 END IF
638 END DO
639
640 CALL tddft_hfx_matrix(work_hmat_asymm, work_rho_ia_ao_asymm, qs_env)
641
642 IF (do_admm) THEN
643 DO ispin = 1, nspins
644 CALL cp_dbcsr_sm_fm_multiply(work_hmat_asymm(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
645 ncol=nao, alpha=1.0_dp, beta=0.0_dp)
646
647 CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
648 admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb)
649
650 CALL parallel_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, &
651 gs_mos(ispin)%mos_occ, 1.0_dp, aop_evects(ispin, ivect))
652 END DO
653 ELSE
654 DO ispin = 1, nspins
655 CALL cp_dbcsr_sm_fm_multiply(work_hmat_asymm(ispin)%matrix, gs_mos(ispin)%mos_occ, &
656 aop_evects(ispin, ivect), ncol=nactive(ispin), &
657 alpha=alpha, beta=1.0_dp)
658 END DO
659 END IF
660 END DO
661 END IF
662
663 CALL timestop(handle)
664
665 END SUBROUTINE tddfpt_apply_hfx
666
667! **************************************************************************************************
668!> \brief Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
669!> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
670!> \param evects trial vectors
671!> \param gs_mos molecular orbitals optimised for the ground state (only occupied
672!> molecular orbitals [component %mos_occ] are needed)
673!> \param qs_env Quickstep environment
674!> \param admm_env ...
675!> \param hfx_section ...
676!> \param x_data ...
677!> \param symmetry ...
678!> \param recalc_integrals ...
679!> \param work_rho_ia_ao ...
680!> \param work_hmat ...
681!> \param wfm_rho_orb ...
682! **************************************************************************************************
683 SUBROUTINE tddfpt_apply_hfxsr_kernel(Aop_evects, evects, gs_mos, qs_env, admm_env, &
684 hfx_section, x_data, symmetry, recalc_integrals, &
685 work_rho_ia_ao, work_hmat, wfm_rho_orb)
686 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: aop_evects, evects
687 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
688 INTENT(in) :: gs_mos
689 TYPE(qs_environment_type), POINTER :: qs_env
690 TYPE(admm_type), POINTER :: admm_env
691 TYPE(section_vals_type), POINTER :: hfx_section
692 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
693 INTEGER, INTENT(IN) :: symmetry
694 LOGICAL, INTENT(IN) :: recalc_integrals
695 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: work_rho_ia_ao
696 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
697 TARGET :: work_hmat
698 TYPE(cp_fm_type), INTENT(IN) :: wfm_rho_orb
699
700 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_apply_hfxsr_kernel'
701
702 INTEGER :: handle, ispin, ivect, nao, nao_aux, &
703 nspins, nvects
704 INTEGER, DIMENSION(maxspins) :: nactive
705 LOGICAL :: reint
706 REAL(kind=dp) :: alpha
707
708 CALL timeset(routinen, handle)
709
710 nspins = SIZE(evects, 1)
711 nvects = SIZE(evects, 2)
712
713 alpha = 2.0_dp
714 IF (nspins > 1) alpha = 1.0_dp
715
716 CALL cp_fm_get_info(gs_mos(1)%mos_occ, nrow_global=nao)
717 CALL cp_fm_get_info(admm_env%A, nrow_global=nao_aux)
718 DO ispin = 1, nspins
719 CALL cp_fm_get_info(evects(ispin, 1), ncol_global=nactive(ispin))
720 END DO
721
722 reint = recalc_integrals
723
724 DO ivect = 1, nvects
725 DO ispin = 1, nspins
726 CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp, evects(ispin, ivect), &
727 gs_mos(ispin)%mos_occ, 0.0_dp, wfm_rho_orb)
728 CALL parallel_gemm('N', 'T', nao, nao, nactive(ispin), 0.5_dp*symmetry, gs_mos(ispin)%mos_occ, &
729 evects(ispin, ivect), 1.0_dp, wfm_rho_orb)
730 CALL dbcsr_set(work_hmat(ispin)%matrix, 0.0_dp)
731 CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, &
732 wfm_rho_orb, 0.0_dp, admm_env%work_aux_orb)
733 CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
734 0.0_dp, admm_env%work_aux_aux)
735 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_rho_ia_ao(ispin)%matrix, keep_sparsity=.true.)
736 END DO
737
738 CALL tddft_hfx_matrix(work_hmat, work_rho_ia_ao, qs_env, .false., reint, hfx_section, x_data)
739 reint = .false.
740
741 DO ispin = 1, nspins
742 CALL cp_dbcsr_sm_fm_multiply(work_hmat(ispin)%matrix, admm_env%A, admm_env%work_aux_orb, &
743 ncol=nao, alpha=1.0_dp, beta=0.0_dp)
744 CALL parallel_gemm('T', 'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, &
745 admm_env%work_aux_orb, 0.0_dp, wfm_rho_orb)
746 CALL parallel_gemm('N', 'N', nao, nactive(ispin), nao, alpha, wfm_rho_orb, &
747 gs_mos(ispin)%mos_occ, 1.0_dp, aop_evects(ispin, ivect))
748 END DO
749 END DO
750
751 CALL timestop(handle)
752
753 END SUBROUTINE tddfpt_apply_hfxsr_kernel
754
755! **************************************************************************************************
756!> \brief ...Calculate the HFXLR kernel contribution by contracting the Lowdin MO coefficients --
757!> transition charges with the exchange-type integrals using the sTDA approximation
758!> \param qs_env ...
759!> \param sub_env ...
760!> \param rcut ...
761!> \param hfx_scale ...
762!> \param work ...
763!> \param X ...
764!> \param res ... vector AX with A being the sTDA matrix and X the Davidson trial vector of the
765!> eigenvalue problem A*X = omega*X
766! **************************************************************************************************
767 SUBROUTINE tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, hfx_scale, work, X, res)
768
769 TYPE(qs_environment_type), POINTER :: qs_env
770 TYPE(tddfpt_subgroup_env_type) :: sub_env
771 REAL(kind=dp), INTENT(IN) :: rcut, hfx_scale
772 TYPE(tddfpt_work_matrices) :: work
773 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: x
774 TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: res
775
776 CHARACTER(len=*), PARAMETER :: routinen = 'tddfpt_apply_hfxlr_kernel'
777
778 INTEGER :: handle, iatom, ispin, jatom, natom, &
779 nsgf, nspins
780 INTEGER, DIMENSION(2) :: nactive
781 REAL(kind=dp) :: dr, eps_filter, fcut, gabr
782 REAL(kind=dp), DIMENSION(3) :: rij
783 REAL(kind=dp), DIMENSION(:, :), POINTER :: pblock
784 TYPE(cell_type), POINTER :: cell
785 TYPE(cp_fm_struct_type), POINTER :: fmstruct
786 TYPE(cp_fm_type) :: cvec
787 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: xtransformed
788 TYPE(cp_fm_type), POINTER :: ct
789 TYPE(dbcsr_iterator_type) :: iter
790 TYPE(dbcsr_type) :: pdens
791 TYPE(dbcsr_type), POINTER :: tempmat
792 TYPE(mp_para_env_type), POINTER :: para_env
793 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
794
795 CALL timeset(routinen, handle)
796
797 ! parameters
798 eps_filter = 1.e-08_dp
799
800 nspins = SIZE(x)
801 DO ispin = 1, nspins
802 CALL cp_fm_get_info(x(ispin), ncol_global=nactive(ispin))
803 END DO
804
805 para_env => sub_env%para_env
806
807 CALL get_qs_env(qs_env, natom=natom, cell=cell, particle_set=particle_set)
808
809 ! calculate Loewdin transformed Davidson trial vector tilde(X)=S^1/2*X
810 ! and tilde(tilde(X))=S^1/2_A*tilde(X)_A
811 ALLOCATE (xtransformed(nspins))
812 DO ispin = 1, nspins
813 NULLIFY (fmstruct)
814 ct => work%ctransformed(ispin)
815 CALL cp_fm_get_info(ct, matrix_struct=fmstruct)
816 CALL cp_fm_create(matrix=xtransformed(ispin), matrix_struct=fmstruct, name="XTRANSFORMED")
817 END DO
818 CALL get_lowdin_x(work%shalf, x, xtransformed)
819
820 DO ispin = 1, nspins
821 ct => work%ctransformed(ispin)
822 CALL cp_fm_get_info(ct, matrix_struct=fmstruct, nrow_global=nsgf)
823 CALL cp_fm_create(cvec, fmstruct)
824 !
825 tempmat => work%shalf
826 CALL dbcsr_create(pdens, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
827 ! P(nu,mu) = SUM_j XT(nu,j)*CT(mu,j)
828 ct => work%ctransformed(ispin)
829 CALL dbcsr_set(pdens, 0.0_dp)
830 CALL cp_dbcsr_plus_fm_fm_t(pdens, xtransformed(ispin), ct, nactive(ispin), &
831 1.0_dp, keep_sparsity=.false.)
832 CALL dbcsr_filter(pdens, eps_filter)
833 ! Apply PP*gab -> PP; gab = gamma_coulomb
834 ! P(nu,mu) = P(nu,mu)*g(a of nu,b of mu)
835 CALL dbcsr_iterator_start(iter, pdens)
836 DO WHILE (dbcsr_iterator_blocks_left(iter))
837 CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock)
838 rij = particle_set(iatom)%r - particle_set(jatom)%r
839 rij = pbc(rij, cell)
840 dr = sqrt(sum(rij(:)**2))
841 gabr = 1._dp/rcut
842 IF (dr < 1.e-6) THEN
843 gabr = 2._dp*gabr/sqrt(3.1415926_dp)
844 ELSE
845 gabr = erf(gabr*dr)/dr
846 fcut = exp(dr - 4._dp*rcut)
847 fcut = fcut/(fcut + 1._dp)
848 END IF
849 pblock = hfx_scale*gabr*pblock
850 END DO
851 CALL dbcsr_iterator_stop(iter)
852 ! CV(mu,i) = P(nu,mu)*CT(mu,i)
853 CALL cp_dbcsr_sm_fm_multiply(pdens, ct, cvec, nactive(ispin), 1.0_dp, 0.0_dp)
854 ! rho(nu,i) = rho(nu,i) + ShalfP(nu,mu)*CV(mu,i)
855 CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, res(ispin), nactive(ispin), &
856 -1.0_dp, 1.0_dp)
857 !
858 CALL dbcsr_release(pdens)
859 !
860 CALL cp_fm_release(cvec)
861 END DO
862
863 CALL cp_fm_release(xtransformed)
864
865 CALL timestop(handle)
866
867 END SUBROUTINE tddfpt_apply_hfxlr_kernel
868
869! **************************************************************************************************
870
871END MODULE qs_tddfpt2_operators
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset)
...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_release(matrix)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
basic linear algebra operations for full matrices
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
represent the structure of a full matrix
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
subroutine, public vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set, para_env, tddft, local_rho_set_2nd, core_2nd)
Calculates one center GAPW Hartree energies and matrix elements Hartree potentials are input Takes po...
Utilities for hfx and admm methods.
subroutine, public tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, external_hfx_sections, external_x_data, external_para_env)
Add the hfx contributions to the Hamiltonian.
Types and set/get functions for HFX.
Definition hfx_types.F:15
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Interface to the message passing library MPI.
basic linear algebra operations for full matrixes
Define the data structure for the particle information.
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
functions related to the poisson solver on regular grids
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs)
Get the QUICKSTEP environment.
subroutine, public integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, local_rho_set_2nd, atener, kforce, my_pools, my_rs_descs)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
subroutine, public tddfpt_apply_hfx(aop_evects, evects, gs_mos, do_admm, qs_env, work_rho_ia_ao_symm, work_hmat_symm, work_rho_ia_ao_asymm, work_hmat_asymm, wfm_rho_orb)
Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
subroutine, public tddfpt_apply_coulomb(a_ia_rspace, rho_ia_g, local_rho_set, hartree_local, qs_env, sub_env, gapw, work_v_gspace, work_v_rspace, tddfpt_mgrid)
Update v_rspace by adding coulomb term.
subroutine, public tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, hfx_scale, work, x, res)
...Calculate the HFXLR kernel contribution by contracting the Lowdin MO coefficients – transition cha...
subroutine, public tddfpt_apply_energy_diff(aop_evects, evects, s_evects, gs_mos, matrix_ks)
Apply orbital energy difference term: Aop_evects(spin,state) += KS(spin) * evects(spin,...
subroutine, public tddfpt_apply_xc(a_ia_rspace, kernel_env, rho_ia_struct, is_rks_triplets, pw_env, work_v_xc, work_v_xc_tau)
Driver routine for applying fxc (analyic vs. finite difference for testing.
subroutine, public tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, admm_env, hfx_section, x_data, symmetry, recalc_integrals, work_rho_ia_ao, work_hmat, wfm_rho_orb)
Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
subroutine, public tddfpt_apply_xc_potential(a_ia_rspace, fxc_rspace, rho_ia_struct, is_rks_triplets)
Routine for applying fxc potential.
Simplified Tamm Dancoff approach (sTDA).
subroutine, public get_lowdin_x(shalf, xvec, xt)
Calculate Lowdin transformed Davidson trial vector X shalf (dbcsr), xvec, xt (fm) are defined in the ...
contains the structure
subroutine, public xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, xc_deriv_method_id, xc_rho_smooth_id, pw_pool)
updates the given rho set with the density given by rho_r (and rho_g). The rho set will contain the c...
Exchange and Correlation functional calculations.
Definition xc.F:17
subroutine, public xc_calc_2nd_deriv_analytical(v_xc, v_xc_tau, deriv_set, rho_set, rho1_set, pw_pool, xc_section, gapw, vxg, tddfpt_fac, compute_virial, virial_xc)
Calculates the second derivative of E_xc at rho in the direction rho1 (if you see the second derivati...
Definition xc.F:2640
subroutine, public xc_calc_2nd_deriv_numerical(v_xc, v_tau, rho_set, rho1_r, rho1_g, tau1_r, pw_pool, xc_section, do_triplet, calc_virial, virial_xc, deriv_set)
calculates 2nd derivative numerically
Definition xc.F:1653
stores some data used in wavefunction fitting
Definition admm_types.F:120
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
keeps the information about the structure of a full matrix
represent a full matrix
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F:509
stores all the informations relevant to an mpi environment
contained for different pw related things
environment for the poisson solver
to create arrays of pools
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Collection of variables required to evaluate adiabatic TDDFPT kernel.
keeps the density in various representations, keeping track of which ones are valid.
Ground state molecular orbitals.
Set of temporary ("work") matrices.
represent a density, with all the representation and data needed to perform a functional evaluation