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