(git:9a40378)
Loading...
Searching...
No Matches
qs_tddfpt2_forces.F
Go to the documentation of this file.
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
8MODULE qs_tddfpt2_forces
9 USE admm_types, ONLY: admm_type,&
10 get_admm_env
11 USE atomic_kind_types, ONLY: atomic_kind_type,&
12 get_atomic_kind,&
13 get_atomic_kind_set
14 USE bibliography, ONLY: hehn2022,&
15 hehn2024,&
16 sertcan2024,&
17 cite_reference
18 USE cp_control_types, ONLY: dft_control_type,&
19 tddfpt2_control_type
20 USE cp_dbcsr_api, ONLY: &
21 dbcsr_add, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_p_type, &
22 dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric
23 USE cp_dbcsr_contrib, ONLY: dbcsr_dot
24 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
25 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
26 copy_fm_to_dbcsr,&
27 cp_dbcsr_plus_fm_fm_t,&
28 cp_dbcsr_sm_fm_multiply,&
29 dbcsr_allocate_matrix_set,&
30 dbcsr_deallocate_matrix_set
31 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
32 cp_fm_struct_release,&
33 cp_fm_struct_type
34 USE cp_fm_types, ONLY: cp_fm_copy_general,&
35 cp_fm_create,&
36 cp_fm_get_info,&
37 cp_fm_release,&
38 cp_fm_set_all,&
39 cp_fm_type
40 USE cp_log_handling, ONLY: cp_get_default_logger,&
41 cp_logger_get_default_unit_nr,&
42 cp_logger_type
43 USE exstates_types, ONLY: excited_energy_type,&
44 exstate_potential_release
45 USE hartree_local_methods, ONLY: vh_1c_gg_integrals,&
46 init_coulomb_local
47 USE hartree_local_types, ONLY: hartree_local_create,&
48 hartree_local_release,&
49 hartree_local_type
50 USE hfx_energy_potential, ONLY: integrate_four_center
51 USE hfx_ri, ONLY: hfx_ri_update_ks
52 USE hfx_types, ONLY: hfx_type
53 USE input_constants, ONLY: do_admm_aux_exch_func_none,&
54 no_sf_tddfpt,&
55 oe_shift,&
56 tddfpt_kernel_full,&
57 tddfpt_kernel_none,&
58 tddfpt_kernel_stda
59 USE input_section_types, ONLY: section_get_lval,&
60 section_vals_get,&
61 section_vals_get_subs_vals,&
62 section_vals_type,&
63 section_vals_val_get
64 USE kinds, ONLY: default_string_length,&
65 dp
66 USE message_passing, ONLY: mp_para_env_type
67 USE mulliken, ONLY: ao_charges
68 USE parallel_gemm_api, ONLY: parallel_gemm
69 USE particle_types, ONLY: particle_type
70 USE pw_env_types, ONLY: pw_env_get,&
71 pw_env_type
72 USE pw_methods, ONLY: pw_axpy,&
73 pw_scale,&
74 pw_transfer,&
75 pw_zero
76 USE pw_poisson_methods, ONLY: pw_poisson_solve
77 USE pw_poisson_types, ONLY: pw_poisson_type
78 USE pw_pool_types, ONLY: pw_pool_type
79 USE pw_types, ONLY: pw_c1d_gs_type,&
80 pw_r3d_rs_type
81 USE qs_collocate_density, ONLY: calculate_rho_elec
82 USE qs_density_matrices, ONLY: calculate_wx_matrix,&
83 calculate_xwx_matrix
84 USE qs_environment_types, ONLY: get_qs_env,&
85 qs_environment_type,&
86 set_qs_env
87 USE qs_force_types, ONLY: allocate_qs_force,&
88 deallocate_qs_force,&
89 qs_force_type,&
90 sum_qs_force,&
91 total_qs_force,&
92 zero_qs_force
93 USE qs_fxc, ONLY: qs_fxc_analytic,&
94 qs_fxc_fdiff
95 USE qs_gapw_densities, ONLY: prepare_gapw_den
96 USE qs_integrate_potential, ONLY: integrate_v_rspace
97 USE qs_kernel_types, ONLY: kernel_env_type
98 USE qs_kind_types, ONLY: get_qs_kind,&
99 get_qs_kind_set,&
100 qs_kind_type
101 USE qs_ks_atom, ONLY: update_ks_atom
102 USE qs_ks_reference, ONLY: ks_ref_potential,&
103 ks_ref_potential_atom
104 USE qs_ks_types, ONLY: qs_ks_env_type
105 USE qs_local_rho_types, ONLY: local_rho_set_create,&
106 local_rho_set_release,&
107 local_rho_type
108 USE qs_mo_types, ONLY: get_mo_set,&
109 mo_set_type
110 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
111 USE qs_oce_types, ONLY: oce_matrix_type
112 USE qs_overlap, ONLY: build_overlap_matrix
113 USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace,&
114 rho0_s_grid_create
115 USE qs_rho0_methods, ONLY: init_rho0
116 USE qs_rho0_types, ONLY: get_rho0_mpole
117 USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals,&
118 calculate_rho_atom_coeff
119 USE qs_rho_atom_types, ONLY: rho_atom_type
120 USE qs_rho_types, ONLY: qs_rho_create,&
121 qs_rho_get,&
122 qs_rho_set,&
123 qs_rho_type
124 USE qs_tddfpt2_fhxc_forces, ONLY: fhxc_force,&
125 stda_force
126 USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
127 USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos,&
128 tddfpt_work_matrices
129 USE qs_vxc_atom, ONLY: calculate_xc_2nd_deriv_atom
130 USE task_list_types, ONLY: task_list_type
131 USE xtb_ehess, ONLY: xtb_coulomb_hessian
132 USE xtb_types, ONLY: get_xtb_atom_param,&
133 xtb_atom_type
134#include "./base/base_uses.f90"
135
136 IMPLICIT NONE
137
138 PRIVATE
139
140 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_forces'
141
142 PUBLIC :: tddfpt_forces_main
143
144! **************************************************************************************************
145
146CONTAINS
147
148! **************************************************************************************************
149!> \brief Perform TDDFPT gradient calculation. This routine calculates the response vector R of Eq. 49
150!> in J. Chem. Theory Comput. 2022, 18, 4186−4202 (https://doi.org/10.1021/acs.jctc.2c00144)
151!> in ex_env%cpmos and a few contributions to the gradient.
152!> \param qs_env Quickstep environment
153!> \param gs_mos ...
154!> \param ex_env Holds: Response vector ex_env%cpmos = R
155!> Difference density ex_env%matrix_pe = T
156!> Matrix ex_env%matrix_hz = H_munu[T]
157!> \param kernel_env ...
158!> \param sub_env ...
159!> \param work_matrices ...
160!> \par History
161!> * 10.2022 created JHU
162! **************************************************************************************************
163 SUBROUTINE tddfpt_forces_main(qs_env, gs_mos, ex_env, kernel_env, sub_env, work_matrices)
164 TYPE(qs_environment_type), POINTER :: qs_env
165 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
166 POINTER :: gs_mos
167 TYPE(excited_energy_type), POINTER :: ex_env
168 TYPE(kernel_env_type) :: kernel_env
169 TYPE(tddfpt_subgroup_env_type) :: sub_env
170 TYPE(tddfpt_work_matrices) :: work_matrices
171
172 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_forces_main'
173
174 INTEGER :: handle, ispin, nspins, spin
175 LOGICAL :: do_sf
176 TYPE(admm_type), POINTER :: admm_env
177 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
178 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pe_asymm, matrix_pe_symm, &
179 matrix_s, matrix_s_aux_fit
180 TYPE(dft_control_type), POINTER :: dft_control
181 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
182
183 CALL timeset(routinen, handle)
184
185 CALL get_qs_env(qs_env, dft_control=dft_control)
186
187 CALL cite_reference(hehn2022)
188 CALL cite_reference(hehn2024)
189 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) CALL cite_reference(sertcan2024)
190
191 nspins = dft_control%nspins
192 tddfpt_control => dft_control%tddfpt2_control
193 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
194 do_sf = .false.
195 ELSE
196 do_sf = .true.
197 END IF
198
199 ! disable RES-TDDFPT for now
200 DO ispin = 1, nspins
201 IF (gs_mos(ispin)%nmo_occ /= gs_mos(ispin)%nmo_active) THEN
202 CALL cp_abort(__location__, "RES-TDDFPT Forces NYA")
203 END IF
204 END DO
205
206 ! rhs of linres equation
207 IF (ASSOCIATED(ex_env%cpmos)) THEN
208 DO ispin = 1, SIZE(ex_env%cpmos)
209 CALL cp_fm_release(ex_env%cpmos(ispin))
210 END DO
211 DEALLOCATE (ex_env%cpmos)
212 END IF
213 ALLOCATE (ex_env%cpmos(nspins))
214 ! Create and initialize rectangular matrices of nao*occ dimension for alpha and beta R vectors
215 ! for the Z-vector equation system: AZ=-R
216 DO ispin = 1, nspins
217 CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=matrix_struct)
218 CALL cp_fm_create(ex_env%cpmos(ispin), matrix_struct)
219 CALL cp_fm_set_all(ex_env%cpmos(ispin), 0.0_dp)
220 END DO
221 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
222 NULLIFY (matrix_pe_asymm, matrix_pe_symm)
223
224 ! Build difference density matrix_pe = X*X^T - (C*X^T*S*X*C^T + (C*X^T*S*X*C^T)^T)/2
225 !
226 CALL dbcsr_allocate_matrix_set(ex_env%matrix_pe, nspins)
227 CALL dbcsr_allocate_matrix_set(matrix_pe_symm, nspins)
228 CALL dbcsr_allocate_matrix_set(matrix_pe_asymm, nspins)
229 DO ispin = 1, nspins
230
231 ! Initialize matrix_pe as a sparse matrix with zeros
232 ALLOCATE (ex_env%matrix_pe(ispin)%matrix)
233 CALL dbcsr_create(ex_env%matrix_pe(ispin)%matrix, template=matrix_s(1)%matrix)
234 CALL dbcsr_copy(ex_env%matrix_pe(ispin)%matrix, matrix_s(1)%matrix)
235 CALL dbcsr_set(ex_env%matrix_pe(ispin)%matrix, 0.0_dp)
236
237 ALLOCATE (matrix_pe_symm(ispin)%matrix)
238 CALL dbcsr_create(matrix_pe_symm(ispin)%matrix, template=matrix_s(1)%matrix)
239 CALL dbcsr_copy(matrix_pe_symm(ispin)%matrix, ex_env%matrix_pe(ispin)%matrix)
240
241 ALLOCATE (matrix_pe_asymm(ispin)%matrix)
242 CALL dbcsr_create(matrix_pe_asymm(ispin)%matrix, template=matrix_s(1)%matrix, &
243 matrix_type=dbcsr_type_antisymmetric)
244 CALL dbcsr_complete_redistribute(ex_env%matrix_pe(ispin)%matrix, matrix_pe_asymm(ispin)%matrix)
245
246 IF (do_sf) THEN
247 spin = 1
248 ELSE
249 spin = ispin
250 END IF
251 ! Add difference density to matrix_pe
252 CALL tddfpt_resvec1(ex_env%evect(spin), gs_mos(spin)%mos_active, &
253 matrix_s(1)%matrix, ex_env%matrix_pe(ispin)%matrix, ispin, do_sf)
254 END DO
255 !
256 ! Project the difference density into auxiliary basis for ADMM
257 IF (dft_control%do_admm) THEN
258 CALL get_qs_env(qs_env, admm_env=admm_env)
259 CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
260 CALL dbcsr_allocate_matrix_set(ex_env%matrix_pe_admm, nspins)
261 DO ispin = 1, nspins
262 ALLOCATE (ex_env%matrix_pe_admm(ispin)%matrix)
263 CALL dbcsr_create(ex_env%matrix_pe_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
264 CALL dbcsr_copy(ex_env%matrix_pe_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
265 CALL dbcsr_set(ex_env%matrix_pe_admm(ispin)%matrix, 0.0_dp)
266 CALL tddfpt_resvec1_admm(ex_env%matrix_pe(ispin)%matrix, &
267 admm_env, ex_env%matrix_pe_admm(ispin)%matrix)
268 END DO
269 END IF
270 !
271 CALL dbcsr_allocate_matrix_set(ex_env%matrix_hz, nspins)
272 DO ispin = 1, nspins
273 ALLOCATE (ex_env%matrix_hz(ispin)%matrix)
274 CALL dbcsr_create(ex_env%matrix_hz(ispin)%matrix, template=matrix_s(1)%matrix)
275 CALL dbcsr_copy(ex_env%matrix_hz(ispin)%matrix, matrix_s(1)%matrix)
276 CALL dbcsr_set(ex_env%matrix_hz(ispin)%matrix, 0.0_dp)
277 END DO
278 ! Calculate first term of R vector: H_{\mu i\sigma}[T]
279 IF (dft_control%qs_control%xtb) THEN
280 CALL tddfpt_resvec2_xtb(qs_env, ex_env%matrix_pe, gs_mos, ex_env%matrix_hz, ex_env%cpmos)
281 ELSE
282 CALL tddfpt_resvec2(qs_env, ex_env%matrix_pe, ex_env%matrix_pe_admm, &
283 gs_mos, ex_env%matrix_hz, ex_env%cpmos)
284 END IF
285 !
286 CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1, SIZE(ex_env%evect, 1))
287 CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_asymm, SIZE(ex_env%evect, 1))
288 DO ispin = 1, SIZE(ex_env%evect, 1)
289 ALLOCATE (ex_env%matrix_px1(ispin)%matrix)
290 CALL dbcsr_create(ex_env%matrix_px1(ispin)%matrix, template=matrix_s(1)%matrix)
291 CALL dbcsr_copy(ex_env%matrix_px1(ispin)%matrix, matrix_s(1)%matrix)
292 CALL dbcsr_set(ex_env%matrix_px1(ispin)%matrix, 0.0_dp)
293
294 ALLOCATE (ex_env%matrix_px1_asymm(ispin)%matrix)
295 CALL dbcsr_create(ex_env%matrix_px1_asymm(ispin)%matrix, template=matrix_s(1)%matrix, &
296 matrix_type=dbcsr_type_antisymmetric)
297 CALL dbcsr_complete_redistribute(ex_env%matrix_px1(ispin)%matrix, ex_env%matrix_px1_asymm(ispin)%matrix)
298 END DO
299 ! Kernel ADMM
300 IF (tddfpt_control%do_admm) THEN
301 CALL get_qs_env(qs_env, admm_env=admm_env)
302 CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
303 CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_admm, SIZE(ex_env%evect, 1))
304 CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_admm_asymm, SIZE(ex_env%evect, 1))
305 DO ispin = 1, SIZE(ex_env%evect, 1)
306 ALLOCATE (ex_env%matrix_px1_admm(ispin)%matrix)
307 CALL dbcsr_create(ex_env%matrix_px1_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
308 CALL dbcsr_copy(ex_env%matrix_px1_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
309 CALL dbcsr_set(ex_env%matrix_px1_admm(ispin)%matrix, 0.0_dp)
310
311 ALLOCATE (ex_env%matrix_px1_admm_asymm(ispin)%matrix)
312 CALL dbcsr_create(ex_env%matrix_px1_admm_asymm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix, &
313 matrix_type=dbcsr_type_antisymmetric)
314 CALL dbcsr_complete_redistribute(ex_env%matrix_px1_admm(ispin)%matrix, &
315 ex_env%matrix_px1_admm_asymm(ispin)%matrix)
316 END DO
317 END IF
318 ! TDA forces. Calculates and adds all missing terms for the response vector, Eq. 49.
319 CALL tddfpt_forces(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices)
320 ! Rotate res vector cpmos into original frame of occupied orbitals.
321 CALL tddfpt_resvec3(qs_env, ex_env%cpmos, work_matrices)
322
323 CALL dbcsr_deallocate_matrix_set(matrix_pe_symm)
324 CALL dbcsr_deallocate_matrix_set(matrix_pe_asymm)
325
326 CALL timestop(handle)
327
328 END SUBROUTINE tddfpt_forces_main
329
330! **************************************************************************************************
331!> \brief Calculate direct tddft forces
332!> \param qs_env ...
333!> \param ex_env ...
334!> \param gs_mos ...
335!> \param kernel_env ...
336!> \param sub_env ...
337!> \param work_matrices ...
338!> \par History
339!> * 01.2020 screated [JGH]
340! **************************************************************************************************
341 SUBROUTINE tddfpt_forces(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices)
342
343 TYPE(qs_environment_type), POINTER :: qs_env
344 TYPE(excited_energy_type), POINTER :: ex_env
345 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
346 POINTER :: gs_mos
347 TYPE(kernel_env_type), INTENT(IN) :: kernel_env
348 TYPE(tddfpt_subgroup_env_type) :: sub_env
349 TYPE(tddfpt_work_matrices) :: work_matrices
350
351 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_forces'
352
353 INTEGER :: handle
354 INTEGER, ALLOCATABLE, DIMENSION(:) :: natom_of_kind
355 LOGICAL :: debug_forces
356 REAL(kind=dp) :: ehartree, exc
357 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
358 TYPE(dft_control_type), POINTER :: dft_control
359 TYPE(qs_force_type), DIMENSION(:), POINTER :: ks_force, td_force
360
361 CALL timeset(routinen, handle)
362
363 ! for extended debug output
364 debug_forces = ex_env%debug_forces
365 ! prepare force array
366 CALL get_qs_env(qs_env, dft_control=dft_control, force=ks_force, &
367 atomic_kind_set=atomic_kind_set)
368 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
369 NULLIFY (td_force)
370 CALL allocate_qs_force(td_force, natom_of_kind)
371 DEALLOCATE (natom_of_kind)
372 CALL zero_qs_force(td_force)
373 CALL set_qs_env(qs_env, force=td_force)
374 !
375 IF (dft_control%qs_control%xtb) THEN
376 CALL tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, &
377 work_matrices, debug_forces)
378 ELSE
379 !
380 CALL exstate_potential_release(ex_env)
381 ! Build the values of hartree, fock and exchange-correlation potential on the grid
382 CALL ks_ref_potential(qs_env, ex_env%vh_rspace, ex_env%vxc_rspace, &
383 ex_env%vtau_rspace, ex_env%vadmm_rspace, ehartree, exc)
384 CALL ks_ref_potential_atom(qs_env, ex_env%local_rho_set, ex_env%local_rho_set_admm, &
385 ex_env%vh_rspace)
386 CALL tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, &
387 work_matrices, debug_forces)
388 END IF
389 !
390 ! add TD and KS forces
391 CALL get_qs_env(qs_env, force=td_force)
392 CALL sum_qs_force(ks_force, td_force)
393 CALL set_qs_env(qs_env, force=ks_force)
394 CALL deallocate_qs_force(td_force)
395 !
396 CALL timestop(handle)
397
398 END SUBROUTINE tddfpt_forces
399
400! **************************************************************************************************
401!> \brief Calculate direct tddft forces.
402!> J. Chem. Theory Comput. 2022, 18, 7, 4186–4202 (https://doi.org/10.1021/acs.jctc.2c00144)
403!> \param qs_env ...
404!> \param ex_env Holds on exit
405!> cpmos = R, Response vector, Eq. 49.
406!> matrix_pe = T, Difference density, Eq. 44.
407!> matrix_wx1 = CK[D^X]X^T, Third term of Eq. 51.
408!> matrix_wz = CX^T(\omegaS - K)XC^T, Last term of Eq. 51.
409!> \param gs_mos ...
410!> \param kernel_env ...
411!> \param sub_env ...
412!> \param work_matrices ...
413!> \param debug_forces ...
414!> \par History
415!> * 01.2020 screated [JGH]
416! **************************************************************************************************
417 SUBROUTINE tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, &
418 debug_forces)
419
420 TYPE(qs_environment_type), POINTER :: qs_env
421 TYPE(excited_energy_type), POINTER :: ex_env
422 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
423 POINTER :: gs_mos
424 TYPE(kernel_env_type), INTENT(IN) :: kernel_env
425 TYPE(tddfpt_subgroup_env_type) :: sub_env
426 TYPE(tddfpt_work_matrices) :: work_matrices
427 LOGICAL :: debug_forces
428
429 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_force_direct'
430
431 INTEGER :: handle, iounit, ispin, nact, natom, &
432 nspins, spin
433 LOGICAL :: do_sf
434 REAL(kind=dp) :: evalue
435 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ftot1, ftot2
436 REAL(kind=dp), DIMENSION(3) :: fodeb
437 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
438 TYPE(cp_fm_type), DIMENSION(:), POINTER :: evect
439 TYPE(cp_logger_type), POINTER :: logger
440 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_wx1, &
441 matrix_wz, scrm
442 TYPE(dft_control_type), POINTER :: dft_control
443 TYPE(mp_para_env_type), POINTER :: para_env
444 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
445 POINTER :: sab_orb
446 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
447 TYPE(qs_ks_env_type), POINTER :: ks_env
448 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
449
450 CALL timeset(routinen, handle)
451
452 logger => cp_get_default_logger()
453 IF (logger%para_env%is_source()) THEN
454 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
455 ELSE
456 iounit = -1
457 END IF
458
459 evect => ex_env%evect
460
461 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, para_env=para_env, &
462 sab_orb=sab_orb, dft_control=dft_control, force=force)
463 NULLIFY (tddfpt_control)
464 tddfpt_control => dft_control%tddfpt2_control
465 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
466 do_sf = .false.
467 ELSE
468 do_sf = .true.
469 END IF
470 nspins = dft_control%nspins
471
472 IF (debug_forces) THEN
473 CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
474 ALLOCATE (ftot1(3, natom))
475 CALL total_qs_force(ftot1, force, atomic_kind_set)
476 END IF
477
478 ! Build last terms of the response vector, Eq. 49, and third term of Lambda_munu, Eq. 51.
479 ! the response vector is in ex_env%cpmos and Lambda is in ex_env%matrix_wx1
480 CALL tddfpt_kernel_force(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, debug_forces)
481
482 ! Overlap matrix, build the Lambda matrix, Eq. 51.
483 NULLIFY (matrix_wx1, matrix_wz)
484 CALL dbcsr_allocate_matrix_set(matrix_wz, nspins)
485 matrix_wx1 => ex_env%matrix_wx1
486 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, matrix_ks=matrix_ks)
487 DO ispin = 1, nspins
488 IF (do_sf) THEN
489 spin = 1
490 ELSE
491 spin = ispin
492 END IF
493 ! Create and initialize the Lambda matrix as a sparse matrix
494 ALLOCATE (matrix_wz(ispin)%matrix)
495 CALL dbcsr_create(matrix=matrix_wz(ispin)%matrix, template=matrix_s(1)%matrix)
496 CALL cp_dbcsr_alloc_block_from_nbl(matrix_wz(ispin)%matrix, sab_orb)
497 CALL dbcsr_set(matrix_wz(ispin)%matrix, 0.0_dp)
498 ! For spin-flip excitations only the beta component of the Lambda matrix
499 ! contains the excitation energy term
500 IF (.NOT. (do_sf .AND. (ispin == 1))) THEN
501 CALL cp_fm_get_info(evect(spin), ncol_global=nact)
502 CALL cp_dbcsr_plus_fm_fm_t(matrix_wz(ispin)%matrix, matrix_v=evect(spin), ncol=nact)
503 evalue = ex_env%evalue
504 IF (tddfpt_control%oe_corr == oe_shift) THEN
505 evalue = ex_env%evalue - tddfpt_control%ev_shift
506 END IF
507 CALL dbcsr_scale(matrix_wz(ispin)%matrix, evalue)
508 END IF
509 ! For spin-flip excitations only the alpha component of the Lambda matrix
510 ! contains the occupied MO energy term
511 IF (.NOT. (do_sf .AND. (ispin == 2))) THEN
512 CALL calculate_wx_matrix(gs_mos(ispin)%mos_active, evect(spin), matrix_ks(ispin)%matrix, &
513 matrix_wz(ispin)%matrix)
514 END IF
515 END DO
516 IF (nspins == 2) THEN
517 CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
518 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
519 END IF
520 NULLIFY (scrm)
521 IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
522 ! Calculate the force contribution of matrix_xz into the force in ks_env.
523 ! force%overlap = Tr(dS*matrix_wz), last term of Eq. 51.
524 CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
525 matrix_name="OVERLAP MATRIX", &
526 basis_type_a="ORB", basis_type_b="ORB", &
527 sab_nl=sab_orb, calculate_forces=.true., &
528 matrix_p=matrix_wz(1)%matrix)
529 CALL dbcsr_deallocate_matrix_set(scrm)
530 CALL dbcsr_deallocate_matrix_set(matrix_wz)
531 IF (debug_forces) THEN
532 fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
533 CALL para_env%sum(fodeb)
534 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wx*dS ", fodeb
535 END IF
536
537 ! Overlap matrix. Build a part of the first term of Lamda, Eq. 51, corresponding to
538 ! the second term of Eq. 41. matrix_wz = C*X^T*(omega*S - K)*X*C^T
539 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, matrix_ks=matrix_ks)
540 NULLIFY (matrix_wz)
541 CALL dbcsr_allocate_matrix_set(matrix_wz, nspins)
542 DO ispin = 1, nspins
543 ALLOCATE (matrix_wz(ispin)%matrix)
544 CALL dbcsr_create(matrix=matrix_wz(ispin)%matrix, template=matrix_s(1)%matrix)
545 CALL cp_dbcsr_alloc_block_from_nbl(matrix_wz(ispin)%matrix, sab_orb)
546 CALL dbcsr_set(matrix_wz(ispin)%matrix, 0.0_dp)
547 ! For spin-flip excitations only the alpha component of Lambda has contributions
548 ! of this term, so skip beta
549 IF (.NOT. (do_sf .AND. (ispin == 2))) THEN
550 evalue = ex_env%evalue
551 IF (tddfpt_control%oe_corr == oe_shift) THEN
552 evalue = ex_env%evalue - tddfpt_control%ev_shift
553 END IF
554 IF (do_sf) THEN
555 spin = 2
556 ELSE
557 spin = ispin
558 END IF
559 CALL calculate_xwx_matrix(gs_mos(ispin)%mos_active, evect(ispin), matrix_s(1)%matrix, &
560 matrix_ks(spin)%matrix, matrix_wz(ispin)%matrix, evalue)
561 END IF
562 END DO
563 IF (nspins == 2) THEN
564 CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
565 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
566 END IF
567 NULLIFY (scrm)
568 IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
569 CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
570 matrix_name="OVERLAP MATRIX", &
571 basis_type_a="ORB", basis_type_b="ORB", &
572 sab_nl=sab_orb, calculate_forces=.true., &
573 matrix_p=matrix_wz(1)%matrix)
574 CALL dbcsr_deallocate_matrix_set(scrm)
575 CALL dbcsr_deallocate_matrix_set(matrix_wz)
576 IF (debug_forces) THEN
577 fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
578 CALL para_env%sum(fodeb)
579 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: xWx*dS ", fodeb
580 END IF
581
582 ! Compute force contribution of the first term of Eq. 41 in the first term of Eq. 51
583 ! that was calculated in tddfpt_kernel_force,
584 ! force%overlap = 0.5C*H[T]*C^T
585 IF (ASSOCIATED(matrix_wx1)) THEN
586 IF (nspins == 2 .AND. .NOT. do_sf) THEN
587 CALL dbcsr_add(matrix_wx1(1)%matrix, matrix_wx1(2)%matrix, &
588 alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
589 ELSE IF (nspins == 2 .AND. do_sf) THEN
590 CALL dbcsr_add(matrix_wx1(1)%matrix, matrix_wx1(2)%matrix, &
591 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
592 END IF
593 NULLIFY (scrm)
594 IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
595 CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
596 matrix_name="OVERLAP MATRIX", &
597 basis_type_a="ORB", basis_type_b="ORB", &
598 sab_nl=sab_orb, calculate_forces=.true., &
599 matrix_p=matrix_wx1(1)%matrix)
600 CALL dbcsr_deallocate_matrix_set(scrm)
601 IF (debug_forces) THEN
602 fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
603 CALL para_env%sum(fodeb)
604 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: D^XKP*dS ", fodeb
605 END IF
606 END IF
607
608 IF (debug_forces) THEN
609 ALLOCATE (ftot2(3, natom))
610 CALL total_qs_force(ftot2, force, atomic_kind_set)
611 fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
612 CALL para_env%sum(fodeb)
613 IF (iounit > 0) WRITE (iounit, "(T3,A,T30,3F16.8)") "DEBUG:: Excitation Force", fodeb
614 DEALLOCATE (ftot1, ftot2)
615 END IF
616
617 CALL timestop(handle)
618
619 END SUBROUTINE tddfpt_force_direct
620
621! **************************************************************************************************
622!> \brief Build the spin difference density,
623!> matrix_pe = matrix_pe + X*X^T - (C*X^T*S*X*C^T + (C*X^T*S*X*C^T)^T)/2
624!> \param evect ...
625!> \param mos_active ...
626!> \param matrix_s ...
627!> \param matrix_pe ...
628!> \param spin ...
629!> \param do_sf ...
630! **************************************************************************************************
631 SUBROUTINE tddfpt_resvec1(evect, mos_active, matrix_s, matrix_pe, spin, do_sf)
632
633 TYPE(cp_fm_type), INTENT(IN) :: evect, mos_active
634 TYPE(dbcsr_type), POINTER :: matrix_s, matrix_pe
635 INTEGER, INTENT(IN) :: spin
636 LOGICAL, INTENT(IN) :: do_sf
637
638 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_resvec1'
639
640 INTEGER :: handle, iounit, nact, nao, norb
641 REAL(kind=dp) :: tmp
642 TYPE(cp_fm_struct_type), POINTER :: fmstruct, fmstruct2
643 TYPE(cp_fm_type) :: cxmat, xxmat
644 TYPE(cp_logger_type), POINTER :: logger
645
646 CALL timeset(routinen, handle)
647 CALL cp_fm_get_info(mos_active, nrow_global=nao, ncol_global=norb)
648 CALL cp_fm_get_info(evect, nrow_global=nao, ncol_global=nact)
649 cpassert(norb == nact)
650
651 ! matrix_pe = X*X^T
652 IF (.NOT. do_sf .OR. (do_sf .AND. (spin == 2))) THEN
653 CALL cp_dbcsr_plus_fm_fm_t(matrix_pe, matrix_v=evect, ncol=norb)
654 END IF
655
656 ! matrix_pe = matrix_pe - (C*X^T*S*X*C^T + (C*X^T*S*X*C^T)^T)/2
657 IF (.NOT. do_sf .OR. (do_sf .AND. (spin == 1))) THEN
658 CALL cp_fm_get_info(evect, matrix_struct=fmstruct)
659 NULLIFY (fmstruct2)
660 CALL cp_fm_struct_create(fmstruct=fmstruct2, template_fmstruct=fmstruct, &
661 nrow_global=norb, ncol_global=norb)
662 CALL cp_fm_create(xxmat, matrix_struct=fmstruct2)
663 CALL cp_fm_struct_release(fmstruct2)
664 CALL cp_fm_create(cxmat, matrix_struct=fmstruct)
665 ! S*X
666 CALL cp_dbcsr_sm_fm_multiply(matrix_s, evect, cxmat, norb, alpha=1.0_dp, beta=0.0_dp)
667 ! (S*X)^T*X
668 CALL parallel_gemm('T', 'N', norb, norb, nao, 1.0_dp, cxmat, evect, 0.0_dp, xxmat)
669 ! C*X^T*S*X
670 CALL parallel_gemm('N', 'N', nao, norb, norb, 1.0_dp, mos_active, xxmat, 0.0_dp, cxmat)
671 CALL cp_fm_release(xxmat)
672 ! matrix_pe = matrix_pe - (C*(C^T*X^T*S*X)^T + C^T*(C^T*X^T*S*X))/2
673 CALL cp_dbcsr_plus_fm_fm_t(matrix_pe, matrix_v=mos_active, matrix_g=cxmat, &
674 ncol=norb, alpha=-1.0_dp, symmetry_mode=1)
675 CALL cp_fm_release(cxmat)
676 END IF
677 !
678 ! Test for Tr(Pe*S)=0
679 CALL dbcsr_dot(matrix_pe, matrix_s, tmp)
680 IF (.NOT. do_sf) THEN
681 IF (abs(tmp) > 1.e-08_dp) THEN
682 logger => cp_get_default_logger()
683 IF (logger%para_env%is_source()) THEN
684 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
685 ELSE
686 iounit = -1
687 END IF
688 cpwarn("Electron count of excitation density matrix is non-zero.")
689 IF (iounit > 0) THEN
690 WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", tmp
691 WRITE (iounit, "(T2,A,/)") repeat("*", 79)
692 END IF
693 END IF
694 ELSE IF (spin == 1) THEN
695 IF (abs(tmp + 1) > 1.e-08_dp) THEN
696 logger => cp_get_default_logger()
697 IF (logger%para_env%is_source()) THEN
698 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
699 ELSE
700 iounit = -1
701 END IF
702 cpwarn("Count of occupied occupation number change is not -1.")
703 IF (iounit > 0) THEN
704 WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", tmp
705 WRITE (iounit, "(T2,A,/)") repeat("*", 79)
706 END IF
707 END IF
708 ELSE IF (spin == 2) THEN
709 IF (abs(tmp - 1) > 1.e-08_dp) THEN
710 logger => cp_get_default_logger()
711 IF (logger%para_env%is_source()) THEN
712 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
713 ELSE
714 iounit = -1
715 END IF
716 cpwarn("Count of unoccupied occupation number change is not 1.")
717 IF (iounit > 0) THEN
718 WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", tmp
719 WRITE (iounit, "(T2,A,/)") repeat("*", 79)
720 END IF
721 END IF
722 END IF
723 !
724
725 CALL timestop(handle)
726
727 END SUBROUTINE tddfpt_resvec1
728
729! **************************************************************************************************
730!> \brief PA = A * P * A(T)
731!> \param matrix_pe ...
732!> \param admm_env ...
733!> \param matrix_pe_admm ...
734! **************************************************************************************************
735 SUBROUTINE tddfpt_resvec1_admm(matrix_pe, admm_env, matrix_pe_admm)
736
737 TYPE(dbcsr_type), POINTER :: matrix_pe
738 TYPE(admm_type), POINTER :: admm_env
739 TYPE(dbcsr_type), POINTER :: matrix_pe_admm
740
741 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_resvec1_admm'
742
743 INTEGER :: handle, nao, nao_aux
744
745 CALL timeset(routinen, handle)
746 !
747 nao_aux = admm_env%nao_aux_fit
748 nao = admm_env%nao_orb
749 !
750 CALL copy_dbcsr_to_fm(matrix_pe, admm_env%work_orb_orb)
751 CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
752 1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
753 admm_env%work_aux_orb)
754 CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
755 1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
756 admm_env%work_aux_aux)
757 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, matrix_pe_admm, keep_sparsity=.true.)
758 !
759 CALL timestop(handle)
760
761 END SUBROUTINE tddfpt_resvec1_admm
762
763! **************************************************************************************************
764!> \brief Calculates the action of the H operator as in the first term of equation 49 in
765!> https://doi.org/10.1021/acs.jctc.2c00144 (J. Chem. Theory Comput. 2022, 18, 4186−4202)
766!> cpmos = H_{\mu i\sigma}[matrix_pe]
767!> \param qs_env ...
768!> \param matrix_pe Input square matrix with the size of the number of atomic orbitals squared nao^2
769!> \param matrix_pe_admm ...
770!> \param gs_mos ...
771!> \param matrix_hz Holds H_{\mu\nu\sigma}[matrix_pe] on exit
772!> \param cpmos ...
773! **************************************************************************************************
774 SUBROUTINE tddfpt_resvec2(qs_env, matrix_pe, matrix_pe_admm, gs_mos, matrix_hz, cpmos)
775
776 TYPE(qs_environment_type), POINTER :: qs_env
777 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pe, matrix_pe_admm
778 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
779 POINTER :: gs_mos
780 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz
781 TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: cpmos
782
783 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_resvec2'
784
785 CHARACTER(LEN=default_string_length) :: basis_type
786 INTEGER :: handle, iounit, ispin, mspin, n_rep_hf, &
787 nao, nao_aux, natom, norb, nspins
788 LOGICAL :: deriv2_analytic, distribute_fock_matrix, &
789 do_hfx, gapw, gapw_xc, &
790 hfx_treat_lsd_in_core, &
791 s_mstruct_changed
792 REAL(kind=dp) :: eh1, focc, rhotot, thartree
793 REAL(kind=dp), DIMENSION(2) :: total_rho
794 REAL(kind=dp), DIMENSION(:), POINTER :: qlm_tot
795 TYPE(admm_type), POINTER :: admm_env
796 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
797 TYPE(cp_fm_type), POINTER :: mos
798 TYPE(cp_logger_type), POINTER :: logger
799 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: msaux
800 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mhz, mpe
801 TYPE(dbcsr_type), POINTER :: dbwork
802 TYPE(dft_control_type), POINTER :: dft_control
803 TYPE(hartree_local_type), POINTER :: hartree_local
804 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
805 TYPE(local_rho_type), POINTER :: local_rho_set, local_rho_set_admm
806 TYPE(mp_para_env_type), POINTER :: para_env
807 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
808 POINTER :: sab, sab_aux_fit
809 TYPE(oce_matrix_type), POINTER :: oce
810 TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
811 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_g_aux, rhoz_g_aux, trho_g, &
812 trho_xc_g
813 TYPE(pw_env_type), POINTER :: pw_env
814 TYPE(pw_poisson_type), POINTER :: poisson_env
815 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
816 TYPE(pw_r3d_rs_type) :: v_hartree_rspace
817 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_r_aux, rhoz_r_aux, tau_r, &
818 trho_r, trho_xc_r, v_xc, v_xc_tau
819 TYPE(pw_r3d_rs_type), POINTER :: weights
820 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
821 TYPE(qs_ks_env_type), POINTER :: ks_env
822 TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit, rho_xc, rhoz_aux, trho
823 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set, rho_atom_set
824 TYPE(section_vals_type), POINTER :: hfx_section, input, xc_section
825 TYPE(task_list_type), POINTER :: task_list
826
827 CALL timeset(routinen, handle)
828
829 NULLIFY (pw_env)
830 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, ks_env=ks_env, &
831 dft_control=dft_control, para_env=para_env)
832 cpassert(ASSOCIATED(pw_env))
833 nspins = dft_control%nspins
834 gapw = dft_control%qs_control%gapw
835 gapw_xc = dft_control%qs_control%gapw_xc
836
837 cpassert(.NOT. dft_control%tddfpt2_control%do_exck)
838 cpassert(.NOT. dft_control%tddfpt2_control%do_hfxsr)
839 cpassert(.NOT. dft_control%tddfpt2_control%do_hfxlr)
840
841 NULLIFY (auxbas_pw_pool, poisson_env)
842 ! gets the tmp grids
843 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
844 poisson_env=poisson_env)
845
846 CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
847 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
848 CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
849
850 ALLOCATE (trho_r(nspins), trho_g(nspins))
851 DO ispin = 1, nspins
852 CALL auxbas_pw_pool%create_pw(trho_r(ispin))
853 CALL auxbas_pw_pool%create_pw(trho_g(ispin))
854 END DO
855 IF (gapw_xc) THEN
856 ALLOCATE (trho_xc_r(nspins), trho_xc_g(nspins))
857 DO ispin = 1, nspins
858 CALL auxbas_pw_pool%create_pw(trho_xc_r(ispin))
859 CALL auxbas_pw_pool%create_pw(trho_xc_g(ispin))
860 END DO
861 END IF
862
863 ! GAPW/GAPW_XC initializations
864 NULLIFY (hartree_local, local_rho_set)
865 IF (gapw) THEN
866 CALL get_qs_env(qs_env, &
867 atomic_kind_set=atomic_kind_set, &
868 natom=natom, &
869 qs_kind_set=qs_kind_set)
870 CALL local_rho_set_create(local_rho_set)
871 CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
872 qs_kind_set, dft_control, para_env)
873 CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control, &
874 zcore=0.0_dp)
875 CALL rho0_s_grid_create(pw_env, local_rho_set%rho0_mpole)
876 CALL hartree_local_create(hartree_local)
877 CALL init_coulomb_local(hartree_local, natom)
878 ELSEIF (gapw_xc) THEN
879 CALL get_qs_env(qs_env, &
880 atomic_kind_set=atomic_kind_set, &
881 qs_kind_set=qs_kind_set)
882 CALL local_rho_set_create(local_rho_set)
883 CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
884 qs_kind_set, dft_control, para_env)
885 END IF
886
887 total_rho = 0.0_dp
888 CALL pw_zero(rho_tot_gspace)
889 DO ispin = 1, nspins
890 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_pe(ispin)%matrix, &
891 rho=trho_r(ispin), &
892 rho_gspace=trho_g(ispin), &
893 soft_valid=gapw, &
894 total_rho=total_rho(ispin))
895 CALL pw_axpy(trho_g(ispin), rho_tot_gspace)
896 IF (gapw_xc) THEN
897 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_pe(ispin)%matrix, &
898 rho=trho_xc_r(ispin), &
899 rho_gspace=trho_xc_g(ispin), &
900 soft_valid=gapw_xc, &
901 total_rho=rhotot)
902 END IF
903 END DO
904
905 ! GAPW o GAPW_XC require the calculation of hard and soft local densities
906 IF (gapw .OR. gapw_xc) THEN
907 CALL get_qs_env(qs_env=qs_env, oce=oce, sab_orb=sab)
908 CALL calculate_rho_atom_coeff(qs_env, matrix_pe, local_rho_set%rho_atom_set, &
909 qs_kind_set, oce, sab, para_env)
910 CALL prepare_gapw_den(qs_env, local_rho_set, do_rho0=gapw)
911 END IF
912 rhotot = sum(total_rho)
913 IF (gapw) THEN
914 CALL get_rho0_mpole(local_rho_set%rho0_mpole, qlm_tot=qlm_tot)
915 rhotot = rhotot + local_rho_set%rho0_mpole%total_rho0_h
916 CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rho_tot_gspace)
917 IF (ASSOCIATED(local_rho_set%rho0_mpole%rhoz_cneo_s_gs)) THEN
918 CALL pw_axpy(local_rho_set%rho0_mpole%rhoz_cneo_s_gs, rho_tot_gspace)
919 END IF
920 END IF
921
922 IF (abs(rhotot) > 1.e-05_dp) THEN
923 logger => cp_get_default_logger()
924 IF (logger%para_env%is_source()) THEN
925 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
926 ELSE
927 iounit = -1
928 END IF
929 cpwarn("Real space electron count of excitation density is non-zero.")
930 IF (iounit > 0) THEN
931 WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", rhotot
932 WRITE (iounit, "(T2,A,/)") repeat("*", 79)
933 END IF
934 END IF
935
936 ! calculate associated hartree potential
937 CALL pw_poisson_solve(poisson_env, rho_tot_gspace, thartree, &
938 v_hartree_gspace)
939 CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
940 CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
941 IF (gapw) THEN
942 CALL vh_1c_gg_integrals(qs_env, thartree, hartree_local%ecoul_1c, &
943 local_rho_set, para_env, tddft=.true.)
944 CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, &
945 calculate_forces=.false., &
946 local_rho_set=local_rho_set)
947 END IF
948
949 ! Fxc*drho term
950 CALL get_qs_env(qs_env, rho=rho)
951 CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
952 !
953 CALL get_qs_env(qs_env, input=input)
954 IF (dft_control%do_admm) THEN
955 CALL get_qs_env(qs_env, admm_env=admm_env)
956 xc_section => admm_env%xc_section_primary
957 ELSE
958 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
959 END IF
960 !
961 deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
962 IF (deriv2_analytic) THEN
963 NULLIFY (v_xc, v_xc_tau, tau_r)
964 CALL get_qs_env(qs_env=qs_env, xcint_weights=weights)
965 IF (gapw_xc) THEN
966 CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
967 CALL qs_fxc_analytic(rho_xc, trho_xc_r, tau_r, xc_section, weights, auxbas_pw_pool, &
968 .false., v_xc, v_xc_tau)
969 ELSE
970 CALL qs_fxc_analytic(rho, trho_r, tau_r, xc_section, weights, auxbas_pw_pool, &
971 .false., v_xc, v_xc_tau)
972 END IF
973 IF (gapw .OR. gapw_xc) THEN
974 CALL get_qs_env(qs_env, rho_atom_set=rho_atom_set)
975 rho1_atom_set => local_rho_set%rho_atom_set
976 CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
977 do_triplet=.false.)
978 END IF
979 ELSE
980 cpabort("NYA 00006")
981 NULLIFY (v_xc, trho)
982 ALLOCATE (trho)
983 CALL qs_rho_create(trho)
984 CALL qs_rho_set(trho, rho_r=trho_r, rho_g=trho_g)
985 CALL qs_fxc_fdiff(ks_env, rho, trho, xc_section, 6, .false., v_xc, v_xc_tau)
986 DEALLOCATE (trho)
987 END IF
988
989 DO ispin = 1, nspins
990 CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
991 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
992 END DO
993 IF (gapw_xc) THEN
994 DO ispin = 1, nspins
995 CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_hartree_rspace, &
996 hmat=matrix_hz(ispin), &
997 calculate_forces=.false.)
998 CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
999 hmat=matrix_hz(ispin), &
1000 gapw=gapw_xc, calculate_forces=.false.)
1001 END DO
1002 ELSE
1003 ! vtot = v_xc(ispin) + v_hartree
1004 DO ispin = 1, nspins
1005 CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
1006 hmat=matrix_hz(ispin), &
1007 gapw=gapw, calculate_forces=.false.)
1008 CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_hartree_rspace, &
1009 hmat=matrix_hz(ispin), &
1010 gapw=gapw, calculate_forces=.false.)
1011 END DO
1012 END IF
1013 IF (gapw .OR. gapw_xc) THEN
1014 mhz(1:nspins, 1:1) => matrix_hz(1:nspins)
1015 mpe(1:nspins, 1:1) => matrix_pe(1:nspins)
1016 CALL update_ks_atom(qs_env, mhz, mpe, forces=.false., &
1017 rho_atom_external=local_rho_set%rho_atom_set)
1018 END IF
1019
1020 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
1021 CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
1022 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1023 DO ispin = 1, nspins
1024 CALL auxbas_pw_pool%give_back_pw(trho_r(ispin))
1025 CALL auxbas_pw_pool%give_back_pw(trho_g(ispin))
1026 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
1027 END DO
1028 DEALLOCATE (trho_r, trho_g, v_xc)
1029 IF (gapw_xc) THEN
1030 DO ispin = 1, nspins
1031 CALL auxbas_pw_pool%give_back_pw(trho_xc_r(ispin))
1032 CALL auxbas_pw_pool%give_back_pw(trho_xc_g(ispin))
1033 END DO
1034 DEALLOCATE (trho_xc_r, trho_xc_g)
1035 END IF
1036 IF (ASSOCIATED(v_xc_tau)) THEN
1037 DO ispin = 1, nspins
1038 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
1039 END DO
1040 DEALLOCATE (v_xc_tau)
1041 END IF
1042 IF (dft_control%do_admm) THEN
1043 IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
1044 ! add ADMM xc_section_aux terms: f_x[rhoz_ADMM]
1045 CALL get_qs_env(qs_env, admm_env=admm_env)
1046 CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit, matrix_s_aux_fit=msaux, &
1047 task_list_aux_fit=task_list)
1048 basis_type = "AUX_FIT"
1049 !
1050 NULLIFY (mpe, mhz)
1051 ALLOCATE (mpe(nspins, 1))
1052 CALL dbcsr_allocate_matrix_set(mhz, nspins, 1)
1053 DO ispin = 1, nspins
1054 ALLOCATE (mhz(ispin, 1)%matrix)
1055 CALL dbcsr_create(mhz(ispin, 1)%matrix, template=msaux(1)%matrix)
1056 CALL dbcsr_copy(mhz(ispin, 1)%matrix, msaux(1)%matrix)
1057 CALL dbcsr_set(mhz(ispin, 1)%matrix, 0.0_dp)
1058 mpe(ispin, 1)%matrix => matrix_pe_admm(ispin)%matrix
1059 END DO
1060 !
1061 ! GAPW/GAPW_XC initializations
1062 NULLIFY (local_rho_set_admm)
1063 IF (admm_env%do_gapw) THEN
1064 basis_type = "AUX_FIT_SOFT"
1065 task_list => admm_env%admm_gapw_env%task_list
1066 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
1067 CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
1068 CALL local_rho_set_create(local_rho_set_admm)
1069 CALL allocate_rho_atom_internals(local_rho_set_admm%rho_atom_set, atomic_kind_set, &
1070 admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
1071 CALL calculate_rho_atom_coeff(qs_env, matrix_pe_admm, &
1072 rho_atom_set=local_rho_set_admm%rho_atom_set, &
1073 qs_kind_set=admm_env%admm_gapw_env%admm_kind_set, &
1074 oce=admm_env%admm_gapw_env%oce, sab=sab_aux_fit, para_env=para_env)
1075 CALL prepare_gapw_den(qs_env, local_rho_set=local_rho_set_admm, &
1076 do_rho0=.false., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
1077 END IF
1078 !
1079 xc_section => admm_env%xc_section_aux
1080 !
1081 NULLIFY (rho_g_aux, rho_r_aux, rhoz_g_aux, rhoz_r_aux)
1082 CALL qs_rho_get(rho_aux_fit, rho_r=rho_r_aux, rho_g=rho_g_aux)
1083 ! rhoz_aux
1084 ALLOCATE (rhoz_r_aux(nspins), rhoz_g_aux(nspins))
1085 DO ispin = 1, nspins
1086 CALL auxbas_pw_pool%create_pw(rhoz_r_aux(ispin))
1087 CALL auxbas_pw_pool%create_pw(rhoz_g_aux(ispin))
1088 END DO
1089 DO ispin = 1, nspins
1090 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpe(ispin, 1)%matrix, &
1091 rho=rhoz_r_aux(ispin), rho_gspace=rhoz_g_aux(ispin), &
1092 basis_type=basis_type, &
1093 task_list_external=task_list)
1094 END DO
1095 !
1096 NULLIFY (v_xc)
1097 deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
1098 IF (deriv2_analytic) THEN
1099 CALL get_qs_env(qs_env=qs_env, xcint_weights=weights)
1100 NULLIFY (tau_r)
1101 CALL qs_fxc_analytic(rho_aux_fit, rhoz_r_aux, tau_r, xc_section, weights, auxbas_pw_pool, &
1102 .false., v_xc, v_xc_tau)
1103 ELSE
1104 cpabort("NYA 00007")
1105 NULLIFY (rhoz_aux)
1106 ALLOCATE (rhoz_aux)
1107 CALL qs_rho_create(rhoz_aux)
1108 CALL qs_rho_set(rhoz_aux, rho_r=rhoz_r_aux, rho_g=rhoz_g_aux)
1109 CALL qs_fxc_fdiff(ks_env, rho_aux_fit, rhoz_aux, xc_section, 6, .false., v_xc, v_xc_tau)
1110 DEALLOCATE (rhoz_aux)
1111 END IF
1112 !
1113 DO ispin = 1, nspins
1114 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
1115 CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
1116 hmat=mhz(ispin, 1), basis_type=basis_type, &
1117 calculate_forces=.false., &
1118 task_list_external=task_list)
1119 END DO
1120 DO ispin = 1, nspins
1121 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
1122 CALL auxbas_pw_pool%give_back_pw(rhoz_r_aux(ispin))
1123 CALL auxbas_pw_pool%give_back_pw(rhoz_g_aux(ispin))
1124 END DO
1125 DEALLOCATE (v_xc, rhoz_r_aux, rhoz_g_aux)
1126 !
1127 IF (admm_env%do_gapw) THEN
1128 rho_atom_set => admm_env%admm_gapw_env%local_rho_set%rho_atom_set
1129 rho1_atom_set => local_rho_set_admm%rho_atom_set
1130 CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, &
1131 para_env, kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
1132 CALL update_ks_atom(qs_env, mhz(:, 1), matrix_pe_admm, forces=.false., tddft=.false., &
1133 rho_atom_external=rho1_atom_set, &
1134 kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
1135 oce_external=admm_env%admm_gapw_env%oce, &
1136 sab_external=sab_aux_fit)
1137 END IF
1138 !
1139 nao = admm_env%nao_orb
1140 nao_aux = admm_env%nao_aux_fit
1141 ALLOCATE (dbwork)
1142 CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
1143 DO ispin = 1, nspins
1144 CALL cp_dbcsr_sm_fm_multiply(mhz(ispin, 1)%matrix, admm_env%A, &
1145 admm_env%work_aux_orb, nao)
1146 CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
1147 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
1148 admm_env%work_orb_orb)
1149 CALL dbcsr_copy(dbwork, matrix_hz(1)%matrix)
1150 CALL dbcsr_set(dbwork, 0.0_dp)
1151 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.true.)
1152 CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
1153 END DO
1154 CALL dbcsr_release(dbwork)
1155 DEALLOCATE (dbwork)
1156 CALL dbcsr_deallocate_matrix_set(mhz)
1157 DEALLOCATE (mpe)
1158 IF (admm_env%do_gapw) THEN
1159 IF (ASSOCIATED(local_rho_set_admm)) CALL local_rho_set_release(local_rho_set_admm)
1160 END IF
1161 END IF
1162 END IF
1163 IF (gapw .OR. gapw_xc) THEN
1164 IF (ASSOCIATED(local_rho_set)) CALL local_rho_set_release(local_rho_set)
1165 IF (ASSOCIATED(hartree_local)) CALL hartree_local_release(hartree_local)
1166 END IF
1167
1168 ! HFX
1169 hfx_section => section_vals_get_subs_vals(xc_section, "HF")
1170 CALL section_vals_get(hfx_section, explicit=do_hfx)
1171 IF (do_hfx) THEN
1172 CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
1173 cpassert(n_rep_hf == 1)
1174 CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
1175 i_rep_section=1)
1176 mspin = 1
1177 IF (hfx_treat_lsd_in_core) mspin = nspins
1178 !
1179 CALL get_qs_env(qs_env=qs_env, rho=rho, x_data=x_data, para_env=para_env, &
1180 s_mstruct_changed=s_mstruct_changed)
1181 distribute_fock_matrix = .true.
1182 IF (dft_control%do_admm) THEN
1183 CALL get_qs_env(qs_env, admm_env=admm_env)
1184 CALL get_admm_env(admm_env, matrix_s_aux_fit=msaux)
1185 NULLIFY (mpe, mhz)
1186 ALLOCATE (mpe(nspins, 1))
1187 CALL dbcsr_allocate_matrix_set(mhz, nspins, 1)
1188 DO ispin = 1, nspins
1189 ALLOCATE (mhz(ispin, 1)%matrix)
1190 CALL dbcsr_create(mhz(ispin, 1)%matrix, template=msaux(1)%matrix)
1191 CALL dbcsr_copy(mhz(ispin, 1)%matrix, msaux(1)%matrix)
1192 CALL dbcsr_set(mhz(ispin, 1)%matrix, 0.0_dp)
1193 mpe(ispin, 1)%matrix => matrix_pe_admm(ispin)%matrix
1194 END DO
1195 IF (x_data(1, 1)%do_hfx_ri) THEN
1196 eh1 = 0.0_dp
1197 CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpe, &
1198 geometry_did_change=s_mstruct_changed, nspins=nspins, &
1199 hf_fraction=x_data(1, 1)%general_parameter%fraction)
1200 ELSE
1201 DO ispin = 1, mspin
1202 eh1 = 0.0
1203 CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpe, hfx_section, &
1204 para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
1205 ispin=ispin)
1206 END DO
1207 END IF
1208 !
1209 cpassert(ASSOCIATED(admm_env%work_aux_orb))
1210 cpassert(ASSOCIATED(admm_env%work_orb_orb))
1211 nao = admm_env%nao_orb
1212 nao_aux = admm_env%nao_aux_fit
1213 ALLOCATE (dbwork)
1214 CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
1215 DO ispin = 1, nspins
1216 CALL cp_dbcsr_sm_fm_multiply(mhz(ispin, 1)%matrix, admm_env%A, &
1217 admm_env%work_aux_orb, nao)
1218 CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
1219 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
1220 admm_env%work_orb_orb)
1221 CALL dbcsr_copy(dbwork, matrix_hz(ispin)%matrix)
1222 CALL dbcsr_set(dbwork, 0.0_dp)
1223 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.true.)
1224 CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
1225 END DO
1226 CALL dbcsr_release(dbwork)
1227 DEALLOCATE (dbwork)
1228 CALL dbcsr_deallocate_matrix_set(mhz)
1229 DEALLOCATE (mpe)
1230 ELSE
1231 NULLIFY (mpe, mhz)
1232 ALLOCATE (mpe(nspins, 1), mhz(nspins, 1))
1233 DO ispin = 1, nspins
1234 mhz(ispin, 1)%matrix => matrix_hz(ispin)%matrix
1235 mpe(ispin, 1)%matrix => matrix_pe(ispin)%matrix
1236 END DO
1237 IF (x_data(1, 1)%do_hfx_ri) THEN
1238 eh1 = 0.0_dp
1239 CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpe, &
1240 geometry_did_change=s_mstruct_changed, nspins=nspins, &
1241 hf_fraction=x_data(1, 1)%general_parameter%fraction)
1242 ELSE
1243 DO ispin = 1, mspin
1244 eh1 = 0.0
1245 CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpe, hfx_section, &
1246 para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
1247 ispin=ispin)
1248 END DO
1249 END IF
1250 DEALLOCATE (mpe, mhz)
1251 END IF
1252 END IF
1253
1254 focc = 4.0_dp
1255 IF (nspins == 2) focc = 2.0_dp
1256 DO ispin = 1, nspins
1257 mos => gs_mos(ispin)%mos_occ
1258 CALL cp_fm_get_info(mos, ncol_global=norb)
1259 CALL cp_dbcsr_sm_fm_multiply(matrix_hz(ispin)%matrix, mos, cpmos(ispin), &
1260 norb, alpha=focc, beta=0.0_dp)
1261 END DO
1262
1263 CALL timestop(handle)
1264
1265 END SUBROUTINE tddfpt_resvec2
1266
1267! **************************************************************************************************
1268!> \brief ...
1269!> \param qs_env ...
1270!> \param matrix_pe ...
1271!> \param gs_mos ...
1272!> \param matrix_hz ...
1273!> \param cpmos ...
1274! **************************************************************************************************
1275 SUBROUTINE tddfpt_resvec2_xtb(qs_env, matrix_pe, gs_mos, matrix_hz, cpmos)
1276
1277 TYPE(qs_environment_type), POINTER :: qs_env
1278 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pe
1279 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1280 POINTER :: gs_mos
1281 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz
1282 TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: cpmos
1283
1284 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_resvec2_xtb'
1285
1286 INTEGER :: atom_a, handle, iatom, ikind, is, ispin, &
1287 na, natom, natorb, nkind, norb, ns, &
1288 nsgf, nspins
1289 INTEGER, DIMENSION(25) :: lao
1290 INTEGER, DIMENSION(5) :: occ
1291 REAL(dp), ALLOCATABLE, DIMENSION(:) :: mcharge, mcharge1
1292 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, charges, charges1
1293 REAL(kind=dp) :: focc
1294 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1295 TYPE(cp_fm_type), POINTER :: mos
1296 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
1297 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
1298 TYPE(dbcsr_type), POINTER :: s_matrix
1299 TYPE(dft_control_type), POINTER :: dft_control
1300 TYPE(mp_para_env_type), POINTER :: para_env
1301 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1302 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1303 TYPE(qs_rho_type), POINTER :: rho
1304 TYPE(xtb_atom_type), POINTER :: xtb_kind
1305
1306 CALL timeset(routinen, handle)
1307
1308 cpassert(ASSOCIATED(matrix_pe))
1309
1310 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
1311 nspins = dft_control%nspins
1312
1313 DO ispin = 1, nspins
1314 CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
1315 END DO
1316
1317 IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
1318 ! Mulliken charges
1319 CALL get_qs_env(qs_env, rho=rho, particle_set=particle_set, &
1320 matrix_s_kp=matrix_s, para_env=para_env)
1321 natom = SIZE(particle_set)
1322 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1323 ALLOCATE (mcharge(natom), charges(natom, 5))
1324 ALLOCATE (mcharge1(natom), charges1(natom, 5))
1325 charges = 0.0_dp
1326 charges1 = 0.0_dp
1327 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
1328 nkind = SIZE(atomic_kind_set)
1329 CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
1330 ALLOCATE (aocg(nsgf, natom))
1331 aocg = 0.0_dp
1332 ALLOCATE (aocg1(nsgf, natom))
1333 aocg1 = 0.0_dp
1334 p_matrix => matrix_p(:, 1)
1335 s_matrix => matrix_s(1, 1)%matrix
1336 CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
1337 CALL ao_charges(matrix_pe, s_matrix, aocg1, para_env)
1338 DO ikind = 1, nkind
1339 CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
1340 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1341 CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
1342 DO iatom = 1, na
1343 atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1344 charges(atom_a, :) = real(occ(:), kind=dp)
1345 DO is = 1, natorb
1346 ns = lao(is) + 1
1347 charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
1348 charges1(atom_a, ns) = charges1(atom_a, ns) - aocg1(is, atom_a)
1349 END DO
1350 END DO
1351 END DO
1352 DEALLOCATE (aocg, aocg1)
1353 DO iatom = 1, natom
1354 mcharge(iatom) = sum(charges(iatom, :))
1355 mcharge1(iatom) = sum(charges1(iatom, :))
1356 END DO
1357 ! Coulomb Kernel
1358 CALL xtb_coulomb_hessian(qs_env, matrix_hz, charges1, mcharge1, mcharge)
1359 !
1360 DEALLOCATE (charges, mcharge, charges1, mcharge1)
1361 END IF
1362
1363 focc = 2.0_dp
1364 IF (nspins == 2) focc = 1.0_dp
1365 DO ispin = 1, nspins
1366 mos => gs_mos(ispin)%mos_occ
1367 CALL cp_fm_get_info(mos, ncol_global=norb)
1368 CALL cp_dbcsr_sm_fm_multiply(matrix_hz(ispin)%matrix, mos, cpmos(ispin), &
1369 norb, alpha=focc, beta=0.0_dp)
1370 END DO
1371
1372 CALL timestop(handle)
1373
1374 END SUBROUTINE tddfpt_resvec2_xtb
1375
1376! **************************************************************************************************
1377!> \brief ...
1378!> \param qs_env ...
1379!> \param cpmos ...
1380!> \param work ...
1381! **************************************************************************************************
1382 SUBROUTINE tddfpt_resvec3(qs_env, cpmos, work)
1383
1384 TYPE(qs_environment_type), POINTER :: qs_env
1385 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: cpmos
1386 TYPE(tddfpt_work_matrices) :: work
1387
1388 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_resvec3'
1389
1390 INTEGER :: handle, ispin, nao, norb, nspins
1391 TYPE(cp_fm_struct_type), POINTER :: fmstruct
1392 TYPE(cp_fm_type) :: cvec, umat
1393 TYPE(cp_fm_type), POINTER :: omos
1394 TYPE(dft_control_type), POINTER :: dft_control
1395 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1396
1397 CALL timeset(routinen, handle)
1398
1399 CALL get_qs_env(qs_env, mos=mos, dft_control=dft_control)
1400 nspins = dft_control%nspins
1401
1402 DO ispin = 1, nspins
1403 CALL get_mo_set(mos(ispin), mo_coeff=omos)
1404 associate(rvecs => cpmos(ispin))
1405 CALL cp_fm_get_info(rvecs, nrow_global=nao, ncol_global=norb)
1406 CALL cp_fm_create(cvec, rvecs%matrix_struct, "cvec")
1407 CALL cp_fm_struct_create(fmstruct, context=rvecs%matrix_struct%context, nrow_global=norb, &
1408 ncol_global=norb, para_env=rvecs%matrix_struct%para_env)
1409 CALL cp_fm_create(umat, fmstruct, "umat")
1410 CALL cp_fm_struct_release(fmstruct)
1411 !
1412 CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, omos, work%S_C0(ispin), 0.0_dp, umat)
1413 CALL cp_fm_copy_general(rvecs, cvec, rvecs%matrix_struct%para_env)
1414 CALL parallel_gemm("N", "T", nao, norb, norb, 1.0_dp, cvec, umat, 0.0_dp, rvecs)
1415 END associate
1416 CALL cp_fm_release(cvec)
1417 CALL cp_fm_release(umat)
1418 END DO
1419
1420 CALL timestop(handle)
1421
1422 END SUBROUTINE tddfpt_resvec3
1423
1424! **************************************************************************************************
1425!> \brief Calculate direct tddft forces
1426!> \param qs_env ...
1427!> \param ex_env ...
1428!> \param gs_mos ...
1429!> \param kernel_env ...
1430!> \param sub_env ...
1431!> \param work_matrices ...
1432!> \param debug_forces ...
1433!> \par History
1434!> * 01.2020 screated [JGH]
1435! **************************************************************************************************
1436 SUBROUTINE tddfpt_kernel_force(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, debug_forces)
1437
1438 TYPE(qs_environment_type), POINTER :: qs_env
1439 TYPE(excited_energy_type), POINTER :: ex_env
1440 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1441 POINTER :: gs_mos
1442 TYPE(kernel_env_type), INTENT(IN) :: kernel_env
1443 TYPE(tddfpt_subgroup_env_type) :: sub_env
1444 TYPE(tddfpt_work_matrices) :: work_matrices
1445 LOGICAL, INTENT(IN) :: debug_forces
1446
1447 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_kernel_force'
1448
1449 INTEGER :: handle
1450 TYPE(dft_control_type), POINTER :: dft_control
1451 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
1452
1453 CALL timeset(routinen, handle)
1454
1455 CALL get_qs_env(qs_env, dft_control=dft_control)
1456 tddfpt_control => dft_control%tddfpt2_control
1457
1458 IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
1459 ! full Kernel
1460 CALL fhxc_force(qs_env, ex_env, gs_mos, kernel_env%full_kernel, debug_forces)
1461 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
1462 ! sTDA Kernel
1463 CALL stda_force(qs_env, ex_env, gs_mos, kernel_env%stda_kernel, sub_env, work_matrices, debug_forces)
1464 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
1465 ! nothing to be done here
1466 ex_env%matrix_wx1 => null()
1467 ELSE
1468 cpabort('Unknown kernel type')
1469 END IF
1470
1471 CALL timestop(handle)
1472
1473 END SUBROUTINE tddfpt_kernel_force
1474
1475END MODULE qs_tddfpt2_forces
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
mulliken::ao_charges
Definition mulliken.F:42
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
pw_poisson_methods::pw_poisson_solve
Definition pw_poisson_methods.F:78
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
admm_types::get_admm_env
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition admm_types.F:593
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::hehn2024
integer, save, public hehn2024
Definition bibliography.F:36
bibliography::sertcan2024
integer, save, public sertcan2024
Definition bibliography.F:36
bibliography::hehn2022
integer, save, public hehn2022
Definition bibliography.F:36
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1178
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1194
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1132
cp_dbcsr_api::dbcsr_complete_redistribute
subroutine, public dbcsr_complete_redistribute(matrix, redist)
...
Definition cp_dbcsr_api.F:319
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:253
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:552
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:214
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
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,...
Definition cp_dbcsr_operations.F:695
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_copy_general
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition cp_fm_types.F:1615
cp_fm_types::cp_fm_get_info
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:1087
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:169
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:556
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
exstates_types
Types for excited states potential energies.
Definition exstates_types.F:14
exstates_types::exstate_potential_release
subroutine, public exstate_potential_release(ex_env)
...
Definition exstates_types.F:150
hartree_local_methods
Definition hartree_local_methods.F:8
hartree_local_methods::init_coulomb_local
subroutine, public init_coulomb_local(hartree_local, natom)
...
Definition hartree_local_methods.F:72
hartree_local_methods::vh_1c_gg_integrals
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...
Definition hartree_local_methods.F:216
hartree_local_types
Definition hartree_local_types.F:8
hartree_local_types::hartree_local_release
subroutine, public hartree_local_release(hartree_local)
...
Definition hartree_local_types.F:142
hartree_local_types::hartree_local_create
subroutine, public hartree_local_create(hartree_local)
...
Definition hartree_local_types.F:128
hfx_energy_potential
Routines to calculate HFX energy and potential.
Definition hfx_energy_potential.F:14
hfx_energy_potential::integrate_four_center
subroutine, public integrate_four_center(qs_env, x_data, ks_matrix, ehfx, rho_ao, hfx_section, para_env, geometry_did_change, irep, distribute_fock_matrix, ispin, nspins)
computes four center integrals for a full basis set and updates the Kohn-Sham-Matrix and energy....
Definition hfx_energy_potential.F:137
hfx_ri
RI-methods for HFX.
Definition hfx_ri.F:12
hfx_ri::hfx_ri_update_ks
subroutine, public hfx_ri_update_ks(qs_env, ri_data, ks_matrix, ehfx, mos, rho_ao, geometry_did_change, nspins, hf_fraction)
...
Definition hfx_ri.F:1041
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::tddfpt_kernel_none
integer, parameter, public tddfpt_kernel_none
Definition input_constants.F:585
input_constants::oe_shift
integer, parameter, public oe_shift
Definition input_constants.F:591
input_constants::do_admm_aux_exch_func_none
integer, parameter, public do_admm_aux_exch_func_none
Definition input_constants.F:875
input_constants::tddfpt_kernel_full
integer, parameter, public tddfpt_kernel_full
Definition input_constants.F:585
input_constants::tddfpt_kernel_stda
integer, parameter, public tddfpt_kernel_stda
Definition input_constants.F:585
input_constants::no_sf_tddfpt
integer, parameter, public no_sf_tddfpt
Definition input_constants.F:588
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
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
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
input_section_types::section_get_lval
logical function, public section_get_lval(section_vals, keyword_name)
...
Definition input_section_types.F:1012
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
mulliken
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition mulliken.F:13
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
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
pw_methods
Definition pw_methods.F:25
pw_poisson_methods
Definition pw_poisson_methods.F:13
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_pool_types
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
pw_types
Definition pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_collocate_density::calculate_rho_elec
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
Definition qs_collocate_density.F:1715
qs_density_matrices
collects routines that calculate density matrices
Definition qs_density_matrices.F:14
qs_density_matrices::calculate_xwx_matrix
subroutine, public calculate_xwx_matrix(mos_active, xvec, s_matrix, ks_matrix, w_matrix, eval)
Calculate the excited state W matrix from the MO eigenvectors, KS matrix.
Definition qs_density_matrices.F:500
qs_density_matrices::calculate_wx_matrix
subroutine, public calculate_wx_matrix(mos_active, xvec, ks_matrix, w_matrix)
Calculate the excited state W matrix from the MO eigenvectors, KS matrix.
Definition qs_density_matrices.F:457
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
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, mimic, 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, sab_cneo, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, mimic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, rhoz_cneo_set, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, harris_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1106
qs_force_types
Definition qs_force_types.F:13
qs_force_types::sum_qs_force
subroutine, public sum_qs_force(qs_force_out, qs_force_in)
Sum up two qs_force entities qs_force_out = qs_force_out + qs_force_in.
Definition qs_force_types.F:306
qs_force_types::deallocate_qs_force
subroutine, public deallocate_qs_force(qs_force)
Deallocate a Quickstep force data structure.
Definition qs_force_types.F:135
qs_force_types::zero_qs_force
subroutine, public zero_qs_force(qs_force)
Initialize a Quickstep force data structure.
Definition qs_force_types.F:262
qs_force_types::allocate_qs_force
subroutine, public allocate_qs_force(qs_force, natom_of_kind)
Allocate a Quickstep force data structure.
Definition qs_force_types.F:81
qs_force_types::total_qs_force
subroutine, public total_qs_force(force, qs_force, atomic_kind_set)
Get current total force.
Definition qs_force_types.F:562
qs_fxc
https://en.wikipedia.org/wiki/Finite_difference_coefficient
Definition qs_fxc.F:27
qs_fxc::qs_fxc_analytic
subroutine, public qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, weights, auxbas_pw_pool, is_triplet, v_xc, v_xc_tau, spinflip)
...
Definition qs_fxc.F:96
qs_fxc::qs_fxc_fdiff
subroutine, public qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, fxc_rho, fxc_tau)
...
Definition qs_fxc.F:165
qs_gapw_densities
Definition qs_gapw_densities.F:9
qs_gapw_densities::prepare_gapw_den
subroutine, public prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external, pw_env_sub)
...
Definition qs_gapw_densities.F:57
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition qs_integrate_potential.F:22
qs_kernel_types
Definition qs_kernel_types.F:8
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:468
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:913
qs_ks_atom
routines that build the Kohn-Sham matrix contributions coming from local atomic densities
Definition qs_ks_atom.F:12
qs_ks_atom::update_ks_atom
subroutine, public update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, kind_set_external, oce_external, sab_external, kscale, kintegral, kforce, fscale)
The correction to the KS matrix due to the GAPW local terms to the hartree and XC contributions is he...
Definition qs_ks_atom.F:110
qs_ks_reference
Calculate the KS reference potentials.
Definition qs_ks_reference.F:14
qs_ks_reference::ks_ref_potential
subroutine, public ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress)
calculate the Kohn-Sham reference potential
Definition qs_ks_reference.F:95
qs_ks_reference::ks_ref_potential_atom
subroutine, public ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm, v_hartree_rspace)
calculate the Kohn-Sham GAPW reference potentials
Definition qs_ks_reference.F:297
qs_ks_types
Definition qs_ks_types.F:15
qs_local_rho_types
Definition qs_local_rho_types.F:8
qs_local_rho_types::local_rho_set_create
subroutine, public local_rho_set_create(local_rho_set)
...
Definition qs_local_rho_types.F:206
qs_local_rho_types::local_rho_set_release
subroutine, public local_rho_set_release(local_rho_set)
...
Definition qs_local_rho_types.F:227
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:399
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_oce_types
Definition qs_oce_types.F:8
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition qs_overlap.F:19
qs_overlap::build_overlap_matrix
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition qs_overlap.F:120
qs_rho0_ggrid
Definition qs_rho0_ggrid.F:8
qs_rho0_ggrid::rho0_s_grid_create
subroutine, public rho0_s_grid_create(pw_env, rho0_mpole)
...
Definition qs_rho0_ggrid.F:280
qs_rho0_ggrid::integrate_vhg0_rspace
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)
...
Definition qs_rho0_ggrid.F:342
qs_rho0_methods
Definition qs_rho0_methods.F:9
qs_rho0_methods::init_rho0
subroutine, public init_rho0(local_rho_set, qs_env, gapw_control, zcore)
...
Definition qs_rho0_methods.F:402
qs_rho0_types
Definition qs_rho0_types.F:8
qs_rho0_types::get_rho0_mpole
subroutine, public get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, mp_gau_ikind, mp_rho, norm_g0l_h, qlm_gg, qlm_car, qlm_tot, zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, max_rpgf0_s, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs)
...
Definition qs_rho0_types.F:449
qs_rho_atom_methods
Definition qs_rho_atom_methods.F:7
qs_rho_atom_methods::allocate_rho_atom_internals
subroutine, public allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
...
Definition qs_rho_atom_methods.F:916
qs_rho_atom_methods::calculate_rho_atom_coeff
subroutine, public calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
...
Definition qs_rho_atom_methods.F:420
qs_rho_atom_types
Definition qs_rho_atom_types.F:8
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_set
subroutine, public qs_rho_set(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)
...
Definition qs_rho_types.F:308
qs_rho_types::qs_rho_get
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
qs_rho_types::qs_rho_create
subroutine, public qs_rho_create(rho)
Allocates a new instance of rho.
Definition qs_rho_types.F:101
qs_tddfpt2_fhxc_forces
Definition qs_tddfpt2_fhxc_forces.F:8
qs_tddfpt2_fhxc_forces::stda_force
subroutine, public stda_force(qs_env, ex_env, gs_mos, stda_env, sub_env, work, debug_forces)
Simplified Tamm Dancoff approach (sTDA). Kernel contribution to forces.
Definition qs_tddfpt2_fhxc_forces.F:1443
qs_tddfpt2_fhxc_forces::fhxc_force
subroutine, public fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
Calculate direct tddft forces. Calculate the three last terms of the response vector in equation 49 a...
Definition qs_tddfpt2_fhxc_forces.F:169
qs_tddfpt2_forces
Definition qs_tddfpt2_forces.F:8
qs_tddfpt2_forces::tddfpt_forces_main
subroutine, public tddfpt_forces_main(qs_env, gs_mos, ex_env, kernel_env, sub_env, work_matrices)
Perform TDDFPT gradient calculation. This routine calculates the response vector R of Eq....
Definition qs_tddfpt2_forces.F:164
qs_tddfpt2_subgroups
Definition qs_tddfpt2_subgroups.F:8
qs_tddfpt2_types
Definition qs_tddfpt2_types.F:8
qs_vxc_atom
routines that build the integrals of the Vxc potential calculated for the atomic density in the basis...
Definition qs_vxc_atom.F:12
qs_vxc_atom::calculate_xc_2nd_deriv_atom
subroutine, public calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, do_tddfpt2, do_triplet, do_sf, kind_set_external)
...
Definition qs_vxc_atom.F:762
task_list_types
types for task lists
Definition task_list_types.F:14
xtb_ehess
Calculation of Coulomb Hessian contributions in xTB.
Definition xtb_ehess.F:12
xtb_ehess::xtb_coulomb_hessian
subroutine, public xtb_coulomb_hessian(qs_env, ks_matrix, charges1, mcharge1, mcharge)
...
Definition xtb_ehess.F:77
xtb_types
Definition of the xTB parameter types.
Definition xtb_types.F:20
xtb_types::get_xtb_atom_param
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, xi, kappa0, alpg, occupation, electronegativity, chmax, en, kqat2, kcn, kq)
...
Definition xtb_types.F:199
admm_types::admm_type
stores some data used in wavefunction fitting
Definition admm_types.F:120
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cp_control_types::dft_control_type
Definition cp_control_types.F:631
cp_control_types::tddfpt2_control_type
Definition cp_control_types.F:510
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:115
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
exstates_types::excited_energy_type
Contains information on the excited states energy.
Definition exstates_types.F:54
hartree_local_types::hartree_local_type
Definition hartree_local_types.F:34
hfx_types::hfx_type
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F:511
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:123
pw_pool_types::pw_pool_type
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:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kernel_types::kernel_env_type
Type to hold environments for the different kernels.
Definition qs_kernel_types.F:83
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_local_rho_types::local_rho_type
Definition qs_local_rho_types.F:45
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_oce_types::oce_matrix_type
Definition qs_oce_types.F:37
qs_rho_atom_types::rho_atom_type
Definition qs_rho_atom_types.F:23
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_tddfpt2_subgroups::tddfpt_subgroup_env_type
Parallel (sub)group environment.
Definition qs_tddfpt2_subgroups.F:103
qs_tddfpt2_types::tddfpt_ground_state_mos
Ground state molecular orbitals.
Definition qs_tddfpt2_types.F:85
qs_tddfpt2_types::tddfpt_work_matrices
Set of temporary ("work") matrices.
Definition qs_tddfpt2_types.F:117
task_list_types::task_list_type
Definition task_list_types.F:58
xtb_types::xtb_atom_type
Definition xtb_types.F:42