Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : MODULE 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: Sertcan2024,&
15 : cite_reference
16 : USE cp_control_types, ONLY: dft_control_type,&
17 : tddfpt2_control_type
18 : USE cp_dbcsr_api, ONLY: &
19 : dbcsr_add, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_p_type, &
20 : dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric
21 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot
22 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
23 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
24 : copy_fm_to_dbcsr,&
25 : cp_dbcsr_plus_fm_fm_t,&
26 : cp_dbcsr_sm_fm_multiply,&
27 : dbcsr_allocate_matrix_set,&
28 : dbcsr_deallocate_matrix_set
29 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
30 : cp_fm_struct_release,&
31 : cp_fm_struct_type
32 : USE cp_fm_types, ONLY: cp_fm_copy_general,&
33 : cp_fm_create,&
34 : cp_fm_get_info,&
35 : cp_fm_release,&
36 : cp_fm_set_all,&
37 : cp_fm_type
38 : USE cp_log_handling, ONLY: cp_get_default_logger,&
39 : cp_logger_get_default_unit_nr,&
40 : cp_logger_type
41 : USE exstates_types, ONLY: excited_energy_type,&
42 : exstate_potential_release
43 : USE hartree_local_methods, ONLY: Vh_1c_gg_integrals,&
44 : init_coulomb_local
45 : USE hartree_local_types, ONLY: hartree_local_create,&
46 : hartree_local_release,&
47 : hartree_local_type
48 : USE hfx_energy_potential, ONLY: integrate_four_center
49 : USE hfx_ri, ONLY: hfx_ri_update_ks
50 : USE hfx_types, ONLY: hfx_type
51 : USE input_constants, ONLY: do_admm_aux_exch_func_none,&
52 : oe_shift,&
53 : tddfpt_kernel_full,&
54 : tddfpt_kernel_none,&
55 : tddfpt_kernel_stda
56 : USE input_section_types, ONLY: section_get_lval,&
57 : section_vals_get,&
58 : section_vals_get_subs_vals,&
59 : section_vals_type,&
60 : section_vals_val_get
61 : USE kinds, ONLY: default_string_length,&
62 : dp
63 : USE message_passing, ONLY: mp_para_env_type
64 : USE mulliken, ONLY: ao_charges
65 : USE parallel_gemm_api, ONLY: parallel_gemm
66 : USE particle_types, ONLY: particle_type
67 : USE pw_env_types, ONLY: pw_env_get,&
68 : pw_env_type
69 : USE pw_methods, ONLY: pw_axpy,&
70 : pw_scale,&
71 : pw_transfer,&
72 : pw_zero
73 : USE pw_poisson_methods, ONLY: pw_poisson_solve
74 : USE pw_poisson_types, ONLY: pw_poisson_type
75 : USE pw_pool_types, ONLY: pw_pool_type
76 : USE pw_types, ONLY: pw_c1d_gs_type,&
77 : pw_r3d_rs_type
78 : USE qs_collocate_density, ONLY: calculate_rho_elec
79 : USE qs_density_matrices, ONLY: calculate_wx_matrix,&
80 : calculate_xwx_matrix
81 : USE qs_environment_types, ONLY: get_qs_env,&
82 : qs_environment_type,&
83 : set_qs_env
84 : USE qs_force_types, ONLY: allocate_qs_force,&
85 : deallocate_qs_force,&
86 : qs_force_type,&
87 : sum_qs_force,&
88 : total_qs_force,&
89 : zero_qs_force
90 : USE qs_fxc, ONLY: qs_fxc_analytic,&
91 : qs_fxc_fdiff
92 : USE qs_gapw_densities, ONLY: prepare_gapw_den
93 : USE qs_integrate_potential, ONLY: integrate_v_rspace
94 : USE qs_kernel_types, ONLY: kernel_env_type
95 : USE qs_kind_types, ONLY: get_qs_kind,&
96 : get_qs_kind_set,&
97 : qs_kind_type
98 : USE qs_ks_atom, ONLY: update_ks_atom
99 : USE qs_ks_reference, ONLY: ks_ref_potential,&
100 : ks_ref_potential_atom
101 : USE qs_ks_types, ONLY: qs_ks_env_type
102 : USE qs_local_rho_types, ONLY: local_rho_set_create,&
103 : local_rho_set_release,&
104 : local_rho_type
105 : USE qs_mo_types, ONLY: get_mo_set,&
106 : mo_set_type
107 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
108 : USE qs_oce_types, ONLY: oce_matrix_type
109 : USE qs_overlap, ONLY: build_overlap_matrix
110 : USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace,&
111 : rho0_s_grid_create
112 : USE qs_rho0_methods, ONLY: init_rho0
113 : USE qs_rho0_types, ONLY: get_rho0_mpole
114 : USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals,&
115 : calculate_rho_atom_coeff
116 : USE qs_rho_atom_types, ONLY: rho_atom_type
117 : USE qs_rho_types, ONLY: qs_rho_create,&
118 : qs_rho_get,&
119 : qs_rho_set,&
120 : qs_rho_type
121 : USE qs_tddfpt2_fhxc_forces, ONLY: fhxc_force,&
122 : stda_force
123 : USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
124 : USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos,&
125 : tddfpt_work_matrices
126 : USE qs_vxc_atom, ONLY: calculate_xc_2nd_deriv_atom
127 : USE task_list_types, ONLY: task_list_type
128 : USE xtb_ehess, ONLY: xtb_coulomb_hessian
129 : USE xtb_types, ONLY: get_xtb_atom_param,&
130 : xtb_atom_type
131 : #include "./base/base_uses.f90"
132 :
133 : IMPLICIT NONE
134 :
135 : PRIVATE
136 :
137 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_forces'
138 :
139 : PUBLIC :: tddfpt_forces_main
140 :
141 : ! **************************************************************************************************
142 :
143 : CONTAINS
144 :
145 : ! **************************************************************************************************
146 : !> \brief Perform TDDFPT gradient calculation.
147 : !> \param qs_env Quickstep environment
148 : !> \param gs_mos ...
149 : !> \param ex_env ...
150 : !> \param kernel_env ...
151 : !> \param sub_env ...
152 : !> \param work_matrices ...
153 : !> \par History
154 : !> * 10.2022 created JHU
155 : ! **************************************************************************************************
156 568 : SUBROUTINE tddfpt_forces_main(qs_env, gs_mos, ex_env, kernel_env, sub_env, work_matrices)
157 : TYPE(qs_environment_type), POINTER :: qs_env
158 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
159 : POINTER :: gs_mos
160 : TYPE(excited_energy_type), POINTER :: ex_env
161 : TYPE(kernel_env_type) :: kernel_env
162 : TYPE(tddfpt_subgroup_env_type) :: sub_env
163 : TYPE(tddfpt_work_matrices) :: work_matrices
164 :
165 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_forces_main'
166 :
167 : INTEGER :: handle, ispin, nspins
168 : TYPE(admm_type), POINTER :: admm_env
169 : TYPE(cp_fm_struct_type), POINTER :: matrix_struct
170 568 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pe_asymm, matrix_pe_symm, &
171 568 : matrix_s, matrix_s_aux_fit
172 : TYPE(dft_control_type), POINTER :: dft_control
173 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
174 :
175 568 : CALL timeset(routineN, handle)
176 :
177 568 : CALL get_qs_env(qs_env, dft_control=dft_control)
178 :
179 568 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) CALL cite_reference(Sertcan2024)
180 :
181 568 : nspins = dft_control%nspins
182 568 : tddfpt_control => dft_control%tddfpt2_control
183 : ! rhs of linres equation
184 568 : IF (ASSOCIATED(ex_env%cpmos)) THEN
185 504 : DO ispin = 1, SIZE(ex_env%cpmos)
186 504 : CALL cp_fm_release(ex_env%cpmos(ispin))
187 : END DO
188 230 : DEALLOCATE (ex_env%cpmos)
189 : END IF
190 2368 : ALLOCATE (ex_env%cpmos(nspins))
191 1232 : DO ispin = 1, nspins
192 664 : CALL cp_fm_get_info(matrix=ex_env%evect(ispin), matrix_struct=matrix_struct)
193 664 : CALL cp_fm_create(ex_env%cpmos(ispin), matrix_struct)
194 1232 : CALL cp_fm_set_all(ex_env%cpmos(ispin), 0.0_dp)
195 : END DO
196 568 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
197 568 : NULLIFY (matrix_pe_asymm, matrix_pe_symm)
198 568 : CALL dbcsr_allocate_matrix_set(ex_env%matrix_pe, nspins)
199 568 : CALL dbcsr_allocate_matrix_set(matrix_pe_symm, nspins)
200 568 : CALL dbcsr_allocate_matrix_set(matrix_pe_asymm, nspins)
201 1232 : DO ispin = 1, nspins
202 664 : ALLOCATE (ex_env%matrix_pe(ispin)%matrix)
203 664 : CALL dbcsr_create(ex_env%matrix_pe(ispin)%matrix, template=matrix_s(1)%matrix)
204 664 : CALL dbcsr_copy(ex_env%matrix_pe(ispin)%matrix, matrix_s(1)%matrix)
205 664 : CALL dbcsr_set(ex_env%matrix_pe(ispin)%matrix, 0.0_dp)
206 :
207 664 : ALLOCATE (matrix_pe_symm(ispin)%matrix)
208 664 : CALL dbcsr_create(matrix_pe_symm(ispin)%matrix, template=matrix_s(1)%matrix)
209 664 : CALL dbcsr_copy(matrix_pe_symm(ispin)%matrix, ex_env%matrix_pe(ispin)%matrix)
210 :
211 664 : ALLOCATE (matrix_pe_asymm(ispin)%matrix)
212 : CALL dbcsr_create(matrix_pe_asymm(ispin)%matrix, template=matrix_s(1)%matrix, &
213 664 : matrix_type=dbcsr_type_antisymmetric)
214 664 : CALL dbcsr_complete_redistribute(ex_env%matrix_pe(ispin)%matrix, matrix_pe_asymm(ispin)%matrix)
215 :
216 : CALL tddfpt_resvec1(ex_env%evect(ispin), gs_mos(ispin)%mos_occ, &
217 1232 : matrix_s(1)%matrix, ex_env%matrix_pe(ispin)%matrix)
218 : END DO
219 : !
220 : ! ground state ADMM!
221 568 : IF (dft_control%do_admm) THEN
222 128 : CALL get_qs_env(qs_env, admm_env=admm_env)
223 128 : CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
224 128 : CALL dbcsr_allocate_matrix_set(ex_env%matrix_pe_admm, nspins)
225 276 : DO ispin = 1, nspins
226 148 : ALLOCATE (ex_env%matrix_pe_admm(ispin)%matrix)
227 148 : CALL dbcsr_create(ex_env%matrix_pe_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
228 148 : CALL dbcsr_copy(ex_env%matrix_pe_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
229 148 : CALL dbcsr_set(ex_env%matrix_pe_admm(ispin)%matrix, 0.0_dp)
230 : CALL tddfpt_resvec1_admm(ex_env%matrix_pe(ispin)%matrix, &
231 276 : admm_env, ex_env%matrix_pe_admm(ispin)%matrix)
232 : END DO
233 : END IF
234 : !
235 568 : CALL dbcsr_allocate_matrix_set(ex_env%matrix_hz, nspins)
236 1232 : DO ispin = 1, nspins
237 664 : ALLOCATE (ex_env%matrix_hz(ispin)%matrix)
238 664 : CALL dbcsr_create(ex_env%matrix_hz(ispin)%matrix, template=matrix_s(1)%matrix)
239 664 : CALL dbcsr_copy(ex_env%matrix_hz(ispin)%matrix, matrix_s(1)%matrix)
240 1232 : CALL dbcsr_set(ex_env%matrix_hz(ispin)%matrix, 0.0_dp)
241 : END DO
242 568 : IF (dft_control%qs_control%xtb) THEN
243 16 : CALL tddfpt_resvec2_xtb(qs_env, ex_env%matrix_pe, gs_mos, ex_env%matrix_hz, ex_env%cpmos)
244 : ELSE
245 : CALL tddfpt_resvec2(qs_env, ex_env%matrix_pe, ex_env%matrix_pe_admm, &
246 552 : gs_mos, ex_env%matrix_hz, ex_env%cpmos)
247 : END IF
248 : !
249 568 : CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1, nspins)
250 568 : CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_asymm, nspins)
251 1232 : DO ispin = 1, nspins
252 664 : ALLOCATE (ex_env%matrix_px1(ispin)%matrix)
253 664 : CALL dbcsr_create(ex_env%matrix_px1(ispin)%matrix, template=matrix_s(1)%matrix)
254 664 : CALL dbcsr_copy(ex_env%matrix_px1(ispin)%matrix, matrix_s(1)%matrix)
255 664 : CALL dbcsr_set(ex_env%matrix_px1(ispin)%matrix, 0.0_dp)
256 :
257 664 : ALLOCATE (ex_env%matrix_px1_asymm(ispin)%matrix)
258 : CALL dbcsr_create(ex_env%matrix_px1_asymm(ispin)%matrix, template=matrix_s(1)%matrix, &
259 664 : matrix_type=dbcsr_type_antisymmetric)
260 1232 : CALL dbcsr_complete_redistribute(ex_env%matrix_px1(ispin)%matrix, ex_env%matrix_px1_asymm(ispin)%matrix)
261 : END DO
262 : ! Kernel ADMM
263 568 : IF (tddfpt_control%do_admm) THEN
264 64 : CALL get_qs_env(qs_env, admm_env=admm_env)
265 64 : CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
266 64 : CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_admm, nspins)
267 64 : CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_admm_asymm, nspins)
268 132 : DO ispin = 1, nspins
269 68 : ALLOCATE (ex_env%matrix_px1_admm(ispin)%matrix)
270 68 : CALL dbcsr_create(ex_env%matrix_px1_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
271 68 : CALL dbcsr_copy(ex_env%matrix_px1_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
272 68 : CALL dbcsr_set(ex_env%matrix_px1_admm(ispin)%matrix, 0.0_dp)
273 :
274 68 : ALLOCATE (ex_env%matrix_px1_admm_asymm(ispin)%matrix)
275 : CALL dbcsr_create(ex_env%matrix_px1_admm_asymm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix, &
276 68 : matrix_type=dbcsr_type_antisymmetric)
277 : CALL dbcsr_complete_redistribute(ex_env%matrix_px1_admm(ispin)%matrix, &
278 132 : ex_env%matrix_px1_admm_asymm(ispin)%matrix)
279 : END DO
280 : END IF
281 : ! TDA forces
282 568 : CALL tddfpt_forces(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices)
283 : ! Rotate res vector cpmos into original frame of occupied orbitals
284 568 : CALL tddfpt_resvec3(qs_env, ex_env%cpmos, work_matrices)
285 :
286 568 : CALL dbcsr_deallocate_matrix_set(matrix_pe_symm)
287 568 : CALL dbcsr_deallocate_matrix_set(matrix_pe_asymm)
288 :
289 568 : CALL timestop(handle)
290 :
291 568 : END SUBROUTINE tddfpt_forces_main
292 :
293 : ! **************************************************************************************************
294 : !> \brief Calculate direct tddft forces
295 : !> \param qs_env ...
296 : !> \param ex_env ...
297 : !> \param gs_mos ...
298 : !> \param kernel_env ...
299 : !> \param sub_env ...
300 : !> \param work_matrices ...
301 : !> \par History
302 : !> * 01.2020 screated [JGH]
303 : ! **************************************************************************************************
304 568 : SUBROUTINE tddfpt_forces(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices)
305 :
306 : TYPE(qs_environment_type), POINTER :: qs_env
307 : TYPE(excited_energy_type), POINTER :: ex_env
308 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
309 : POINTER :: gs_mos
310 : TYPE(kernel_env_type), INTENT(IN) :: kernel_env
311 : TYPE(tddfpt_subgroup_env_type) :: sub_env
312 : TYPE(tddfpt_work_matrices) :: work_matrices
313 :
314 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_forces'
315 :
316 : INTEGER :: handle
317 568 : INTEGER, ALLOCATABLE, DIMENSION(:) :: natom_of_kind
318 : LOGICAL :: debug_forces
319 : REAL(KIND=dp) :: ehartree, exc
320 568 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
321 : TYPE(dft_control_type), POINTER :: dft_control
322 568 : TYPE(qs_force_type), DIMENSION(:), POINTER :: ks_force, td_force
323 :
324 568 : CALL timeset(routineN, handle)
325 :
326 : ! for extended debug output
327 568 : debug_forces = ex_env%debug_forces
328 : ! prepare force array
329 : CALL get_qs_env(qs_env, dft_control=dft_control, force=ks_force, &
330 568 : atomic_kind_set=atomic_kind_set)
331 568 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
332 568 : NULLIFY (td_force)
333 568 : CALL allocate_qs_force(td_force, natom_of_kind)
334 568 : DEALLOCATE (natom_of_kind)
335 568 : CALL zero_qs_force(td_force)
336 568 : CALL set_qs_env(qs_env, force=td_force)
337 : !
338 568 : IF (dft_control%qs_control%xtb) THEN
339 : CALL tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, &
340 16 : work_matrices, debug_forces)
341 : ELSE
342 : !
343 552 : CALL exstate_potential_release(ex_env)
344 : CALL ks_ref_potential(qs_env, ex_env%vh_rspace, ex_env%vxc_rspace, &
345 552 : ex_env%vtau_rspace, ex_env%vadmm_rspace, ehartree, exc)
346 : CALL ks_ref_potential_atom(qs_env, ex_env%local_rho_set, ex_env%local_rho_set_admm, &
347 552 : ex_env%vh_rspace)
348 : CALL tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, &
349 552 : work_matrices, debug_forces)
350 : END IF
351 : !
352 : ! add TD and KS forces
353 568 : CALL get_qs_env(qs_env, force=td_force)
354 568 : CALL sum_qs_force(ks_force, td_force)
355 568 : CALL set_qs_env(qs_env, force=ks_force)
356 568 : CALL deallocate_qs_force(td_force)
357 : !
358 568 : CALL timestop(handle)
359 :
360 568 : END SUBROUTINE tddfpt_forces
361 :
362 : ! **************************************************************************************************
363 : !> \brief Calculate direct tddft forces
364 : !> \param qs_env ...
365 : !> \param ex_env ...
366 : !> \param gs_mos ...
367 : !> \param kernel_env ...
368 : !> \param sub_env ...
369 : !> \param work_matrices ...
370 : !> \param debug_forces ...
371 : !> \par History
372 : !> * 01.2020 screated [JGH]
373 : ! **************************************************************************************************
374 568 : SUBROUTINE tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, &
375 : debug_forces)
376 :
377 : TYPE(qs_environment_type), POINTER :: qs_env
378 : TYPE(excited_energy_type), POINTER :: ex_env
379 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
380 : POINTER :: gs_mos
381 : TYPE(kernel_env_type), INTENT(IN) :: kernel_env
382 : TYPE(tddfpt_subgroup_env_type) :: sub_env
383 : TYPE(tddfpt_work_matrices) :: work_matrices
384 : LOGICAL :: debug_forces
385 :
386 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_force_direct'
387 :
388 : INTEGER :: handle, iounit, ispin, natom, norb, &
389 : nspins
390 : REAL(KIND=dp) :: evalue
391 568 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ftot1, ftot2
392 : REAL(KIND=dp), DIMENSION(3) :: fodeb
393 568 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
394 568 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: evect
395 : TYPE(cp_logger_type), POINTER :: logger
396 568 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_wx1, &
397 568 : matrix_wz, scrm
398 : TYPE(dft_control_type), POINTER :: dft_control
399 : TYPE(mp_para_env_type), POINTER :: para_env
400 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
401 568 : POINTER :: sab_orb
402 568 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
403 : TYPE(qs_ks_env_type), POINTER :: ks_env
404 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
405 :
406 568 : CALL timeset(routineN, handle)
407 :
408 568 : logger => cp_get_default_logger()
409 568 : IF (logger%para_env%is_source()) THEN
410 284 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
411 : ELSE
412 : iounit = -1
413 : END IF
414 :
415 568 : evect => ex_env%evect
416 :
417 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, para_env=para_env, &
418 568 : sab_orb=sab_orb, dft_control=dft_control, force=force)
419 568 : NULLIFY (tddfpt_control)
420 568 : tddfpt_control => dft_control%tddfpt2_control
421 568 : nspins = dft_control%nspins
422 :
423 568 : IF (debug_forces) THEN
424 58 : CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
425 174 : ALLOCATE (ftot1(3, natom))
426 58 : CALL total_qs_force(ftot1, force, atomic_kind_set)
427 : END IF
428 :
429 568 : CALL tddfpt_kernel_force(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, debug_forces)
430 :
431 : ! Overlap matrix
432 568 : matrix_wx1 => ex_env%matrix_wx1
433 568 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, matrix_ks=matrix_ks)
434 568 : NULLIFY (matrix_wz)
435 568 : CALL dbcsr_allocate_matrix_set(matrix_wz, nspins)
436 1232 : DO ispin = 1, nspins
437 664 : ALLOCATE (matrix_wz(ispin)%matrix)
438 664 : CALL dbcsr_create(matrix=matrix_wz(ispin)%matrix, template=matrix_s(1)%matrix)
439 664 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_wz(ispin)%matrix, sab_orb)
440 664 : CALL dbcsr_set(matrix_wz(ispin)%matrix, 0.0_dp)
441 664 : CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=norb)
442 664 : CALL cp_dbcsr_plus_fm_fm_t(matrix_wz(ispin)%matrix, matrix_v=evect(ispin), ncol=norb)
443 664 : evalue = ex_env%evalue
444 664 : IF (tddfpt_control%oe_corr == oe_shift) THEN
445 4 : evalue = ex_env%evalue - tddfpt_control%ev_shift
446 : END IF
447 664 : CALL dbcsr_scale(matrix_wz(ispin)%matrix, evalue)
448 : CALL calculate_wx_matrix(gs_mos(ispin)%mos_occ, evect(ispin), matrix_ks(ispin)%matrix, &
449 1896 : matrix_wz(ispin)%matrix)
450 : END DO
451 568 : IF (nspins == 2) THEN
452 : CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
453 96 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
454 : END IF
455 568 : NULLIFY (scrm)
456 742 : IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
457 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
458 : matrix_name="OVERLAP MATRIX", &
459 : basis_type_a="ORB", basis_type_b="ORB", &
460 : sab_nl=sab_orb, calculate_forces=.TRUE., &
461 568 : matrix_p=matrix_wz(1)%matrix)
462 568 : CALL dbcsr_deallocate_matrix_set(scrm)
463 568 : CALL dbcsr_deallocate_matrix_set(matrix_wz)
464 568 : IF (debug_forces) THEN
465 232 : fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
466 58 : CALL para_env%sum(fodeb)
467 58 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wx*dS ", fodeb
468 : END IF
469 :
470 : ! Overlap matrix
471 568 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, matrix_ks=matrix_ks)
472 568 : NULLIFY (matrix_wz)
473 568 : CALL dbcsr_allocate_matrix_set(matrix_wz, nspins)
474 1232 : DO ispin = 1, nspins
475 664 : ALLOCATE (matrix_wz(ispin)%matrix)
476 664 : CALL dbcsr_create(matrix=matrix_wz(ispin)%matrix, template=matrix_s(1)%matrix)
477 664 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_wz(ispin)%matrix, sab_orb)
478 664 : CALL dbcsr_set(matrix_wz(ispin)%matrix, 0.0_dp)
479 664 : CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=norb)
480 664 : evalue = ex_env%evalue
481 664 : IF (tddfpt_control%oe_corr == oe_shift) THEN
482 4 : evalue = ex_env%evalue - tddfpt_control%ev_shift
483 : END IF
484 : CALL calculate_xwx_matrix(gs_mos(ispin)%mos_occ, evect(ispin), matrix_s(1)%matrix, &
485 1896 : matrix_ks(ispin)%matrix, matrix_wz(ispin)%matrix, evalue)
486 : END DO
487 568 : IF (nspins == 2) THEN
488 : CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
489 96 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
490 : END IF
491 568 : NULLIFY (scrm)
492 742 : IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
493 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
494 : matrix_name="OVERLAP MATRIX", &
495 : basis_type_a="ORB", basis_type_b="ORB", &
496 : sab_nl=sab_orb, calculate_forces=.TRUE., &
497 568 : matrix_p=matrix_wz(1)%matrix)
498 568 : CALL dbcsr_deallocate_matrix_set(scrm)
499 568 : CALL dbcsr_deallocate_matrix_set(matrix_wz)
500 568 : IF (debug_forces) THEN
501 232 : fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
502 58 : CALL para_env%sum(fodeb)
503 58 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: xWx*dS ", fodeb
504 : END IF
505 :
506 : ! Overlap matrix
507 568 : IF (ASSOCIATED(matrix_wx1)) THEN
508 500 : IF (nspins == 2) THEN
509 : CALL dbcsr_add(matrix_wx1(1)%matrix, matrix_wx1(2)%matrix, &
510 96 : alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
511 : END IF
512 500 : NULLIFY (scrm)
513 656 : IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
514 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
515 : matrix_name="OVERLAP MATRIX", &
516 : basis_type_a="ORB", basis_type_b="ORB", &
517 : sab_nl=sab_orb, calculate_forces=.TRUE., &
518 500 : matrix_p=matrix_wx1(1)%matrix)
519 500 : CALL dbcsr_deallocate_matrix_set(scrm)
520 500 : IF (debug_forces) THEN
521 208 : fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
522 52 : CALL para_env%sum(fodeb)
523 52 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: WK*dS ", fodeb
524 : END IF
525 : END IF
526 :
527 568 : IF (debug_forces) THEN
528 174 : ALLOCATE (ftot2(3, natom))
529 58 : CALL total_qs_force(ftot2, force, atomic_kind_set)
530 232 : fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
531 58 : CALL para_env%sum(fodeb)
532 58 : IF (iounit > 0) WRITE (iounit, "(T3,A,T30,3F16.8)") "DEBUG:: Excitation Force", fodeb
533 58 : DEALLOCATE (ftot1, ftot2)
534 : END IF
535 :
536 568 : CALL timestop(handle)
537 :
538 1136 : END SUBROUTINE tddfpt_force_direct
539 :
540 : ! **************************************************************************************************
541 : !> \brief ...
542 : !> \param evect ...
543 : !> \param mos_occ ...
544 : !> \param matrix_s ...
545 : !> \param matrix_pe ...
546 : ! **************************************************************************************************
547 3320 : SUBROUTINE tddfpt_resvec1(evect, mos_occ, matrix_s, matrix_pe)
548 :
549 : TYPE(cp_fm_type), INTENT(IN) :: evect, mos_occ
550 : TYPE(dbcsr_type), POINTER :: matrix_s, matrix_pe
551 :
552 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_resvec1'
553 :
554 : INTEGER :: handle, iounit, nao, norb
555 : REAL(KIND=dp) :: tmp
556 : TYPE(cp_fm_struct_type), POINTER :: fmstruct, fmstruct2
557 : TYPE(cp_fm_type) :: cxmat, xxmat
558 : TYPE(cp_logger_type), POINTER :: logger
559 :
560 664 : CALL timeset(routineN, handle)
561 : ! X*X^T
562 664 : CALL cp_fm_get_info(mos_occ, nrow_global=nao, ncol_global=norb)
563 664 : CALL cp_dbcsr_plus_fm_fm_t(matrix_pe, matrix_v=evect, ncol=norb)
564 : ! X^T*S*X
565 664 : CALL cp_fm_get_info(evect, matrix_struct=fmstruct)
566 664 : NULLIFY (fmstruct2)
567 : CALL cp_fm_struct_create(fmstruct=fmstruct2, template_fmstruct=fmstruct, &
568 664 : nrow_global=norb, ncol_global=norb)
569 664 : CALL cp_fm_create(xxmat, matrix_struct=fmstruct2)
570 664 : CALL cp_fm_struct_release(fmstruct2)
571 664 : CALL cp_fm_create(cxmat, matrix_struct=fmstruct)
572 664 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, evect, cxmat, norb, alpha=1.0_dp, beta=0.0_dp)
573 664 : CALL parallel_gemm('T', 'N', norb, norb, nao, 1.0_dp, cxmat, evect, 0.0_dp, xxmat)
574 664 : CALL parallel_gemm('N', 'N', nao, norb, norb, 1.0_dp, mos_occ, xxmat, 0.0_dp, cxmat)
575 664 : CALL cp_fm_release(xxmat)
576 : ! C*C^T*XX
577 : CALL cp_dbcsr_plus_fm_fm_t(matrix_pe, matrix_v=mos_occ, matrix_g=cxmat, &
578 664 : ncol=norb, alpha=-1.0_dp, symmetry_mode=1)
579 664 : CALL cp_fm_release(cxmat)
580 : !
581 : ! Test for Tr(Pe*S)=0
582 664 : CALL dbcsr_dot(matrix_pe, matrix_s, tmp)
583 664 : IF (ABS(tmp) > 1.e-08_dp) THEN
584 0 : logger => cp_get_default_logger()
585 0 : IF (logger%para_env%is_source()) THEN
586 0 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
587 : ELSE
588 : iounit = -1
589 : END IF
590 0 : CPWARN("Electron count of excitation density matrix is non-zero.")
591 0 : IF (iounit > 0) THEN
592 0 : WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", tmp
593 0 : WRITE (iounit, "(T2,A,/)") REPEAT("*", 79)
594 : END IF
595 : END IF
596 : !
597 :
598 664 : CALL timestop(handle)
599 :
600 664 : END SUBROUTINE tddfpt_resvec1
601 :
602 : ! **************************************************************************************************
603 : !> \brief PA = A * P * A(T)
604 : !> \param matrix_pe ...
605 : !> \param admm_env ...
606 : !> \param matrix_pe_admm ...
607 : ! **************************************************************************************************
608 148 : SUBROUTINE tddfpt_resvec1_admm(matrix_pe, admm_env, matrix_pe_admm)
609 :
610 : TYPE(dbcsr_type), POINTER :: matrix_pe
611 : TYPE(admm_type), POINTER :: admm_env
612 : TYPE(dbcsr_type), POINTER :: matrix_pe_admm
613 :
614 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_resvec1_admm'
615 :
616 : INTEGER :: handle, nao, nao_aux
617 :
618 148 : CALL timeset(routineN, handle)
619 : !
620 148 : nao_aux = admm_env%nao_aux_fit
621 148 : nao = admm_env%nao_orb
622 : !
623 148 : CALL copy_dbcsr_to_fm(matrix_pe, admm_env%work_orb_orb)
624 : CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
625 : 1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
626 148 : admm_env%work_aux_orb)
627 : CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
628 : 1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
629 148 : admm_env%work_aux_aux)
630 148 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, matrix_pe_admm, keep_sparsity=.TRUE.)
631 : !
632 148 : CALL timestop(handle)
633 :
634 148 : END SUBROUTINE tddfpt_resvec1_admm
635 :
636 : ! **************************************************************************************************
637 : !> \brief ...
638 : !> \param qs_env ...
639 : !> \param matrix_pe ...
640 : !> \param matrix_pe_admm ...
641 : !> \param gs_mos ...
642 : !> \param matrix_hz ...
643 : !> \param cpmos ...
644 : ! **************************************************************************************************
645 552 : SUBROUTINE tddfpt_resvec2(qs_env, matrix_pe, matrix_pe_admm, gs_mos, matrix_hz, cpmos)
646 :
647 : TYPE(qs_environment_type), POINTER :: qs_env
648 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pe, matrix_pe_admm
649 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
650 : POINTER :: gs_mos
651 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz
652 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: cpmos
653 :
654 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_resvec2'
655 :
656 : CHARACTER(LEN=default_string_length) :: basis_type
657 : INTEGER :: handle, iounit, ispin, mspin, n_rep_hf, &
658 : nao, nao_aux, natom, norb, nspins
659 : LOGICAL :: deriv2_analytic, distribute_fock_matrix, &
660 : do_hfx, gapw, gapw_xc, &
661 : hfx_treat_lsd_in_core, &
662 : s_mstruct_changed
663 : REAL(KIND=dp) :: eh1, focc, rhotot, thartree
664 : REAL(KIND=dp), DIMENSION(2) :: total_rho
665 552 : REAL(KIND=dp), DIMENSION(:), POINTER :: Qlm_tot
666 : TYPE(admm_type), POINTER :: admm_env
667 552 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
668 : TYPE(cp_fm_type), POINTER :: mos
669 : TYPE(cp_logger_type), POINTER :: logger
670 552 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: msaux
671 552 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mhz, mpe
672 : TYPE(dbcsr_type), POINTER :: dbwork
673 : TYPE(dft_control_type), POINTER :: dft_control
674 : TYPE(hartree_local_type), POINTER :: hartree_local
675 552 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
676 : TYPE(local_rho_type), POINTER :: local_rho_set, local_rho_set_admm
677 : TYPE(mp_para_env_type), POINTER :: para_env
678 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
679 552 : POINTER :: sab, sab_aux_fit
680 : TYPE(oce_matrix_type), POINTER :: oce
681 : TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
682 552 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_g_aux, rhoz_g_aux, trho_g, &
683 552 : trho_xc_g
684 : TYPE(pw_env_type), POINTER :: pw_env
685 : TYPE(pw_poisson_type), POINTER :: poisson_env
686 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
687 : TYPE(pw_r3d_rs_type) :: v_hartree_rspace
688 552 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_r_aux, rhoz_r_aux, tau_r, &
689 552 : trho_r, trho_xc_r, v_xc, v_xc_tau
690 552 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
691 : TYPE(qs_ks_env_type), POINTER :: ks_env
692 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit, rho_xc, rhoz_aux, trho
693 552 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set, rho_atom_set
694 : TYPE(section_vals_type), POINTER :: hfx_section, input, xc_section
695 : TYPE(task_list_type), POINTER :: task_list
696 :
697 552 : CALL timeset(routineN, handle)
698 :
699 552 : NULLIFY (pw_env)
700 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, ks_env=ks_env, &
701 552 : dft_control=dft_control, para_env=para_env)
702 552 : CPASSERT(ASSOCIATED(pw_env))
703 552 : nspins = dft_control%nspins
704 552 : gapw = dft_control%qs_control%gapw
705 552 : gapw_xc = dft_control%qs_control%gapw_xc
706 :
707 552 : CPASSERT(.NOT. dft_control%tddfpt2_control%do_exck)
708 552 : CPASSERT(.NOT. dft_control%tddfpt2_control%do_hfxsr)
709 552 : CPASSERT(.NOT. dft_control%tddfpt2_control%do_hfxlr)
710 :
711 552 : NULLIFY (auxbas_pw_pool, poisson_env)
712 : ! gets the tmp grids
713 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
714 552 : poisson_env=poisson_env)
715 :
716 552 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
717 552 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
718 552 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
719 :
720 4056 : ALLOCATE (trho_r(nspins), trho_g(nspins))
721 1200 : DO ispin = 1, nspins
722 648 : CALL auxbas_pw_pool%create_pw(trho_r(ispin))
723 1200 : CALL auxbas_pw_pool%create_pw(trho_g(ispin))
724 : END DO
725 552 : IF (gapw_xc) THEN
726 70 : ALLOCATE (trho_xc_r(nspins), trho_xc_g(nspins))
727 28 : DO ispin = 1, nspins
728 14 : CALL auxbas_pw_pool%create_pw(trho_xc_r(ispin))
729 28 : CALL auxbas_pw_pool%create_pw(trho_xc_g(ispin))
730 : END DO
731 : END IF
732 :
733 : ! GAPW/GAPW_XC initializations
734 552 : NULLIFY (hartree_local, local_rho_set)
735 552 : IF (gapw) THEN
736 : CALL get_qs_env(qs_env, &
737 : atomic_kind_set=atomic_kind_set, &
738 : natom=natom, &
739 64 : qs_kind_set=qs_kind_set)
740 64 : CALL local_rho_set_create(local_rho_set)
741 : CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
742 64 : qs_kind_set, dft_control, para_env)
743 : CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control, &
744 64 : zcore=0.0_dp)
745 64 : CALL rho0_s_grid_create(pw_env, local_rho_set%rho0_mpole)
746 64 : CALL hartree_local_create(hartree_local)
747 64 : CALL init_coulomb_local(hartree_local, natom)
748 488 : ELSEIF (gapw_xc) THEN
749 : CALL get_qs_env(qs_env, &
750 : atomic_kind_set=atomic_kind_set, &
751 14 : qs_kind_set=qs_kind_set)
752 14 : CALL local_rho_set_create(local_rho_set)
753 : CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
754 14 : qs_kind_set, dft_control, para_env)
755 : END IF
756 :
757 552 : total_rho = 0.0_dp
758 552 : CALL pw_zero(rho_tot_gspace)
759 1200 : DO ispin = 1, nspins
760 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_pe(ispin)%matrix, &
761 : rho=trho_r(ispin), &
762 : rho_gspace=trho_g(ispin), &
763 : soft_valid=gapw, &
764 648 : total_rho=total_rho(ispin))
765 648 : CALL pw_axpy(trho_g(ispin), rho_tot_gspace)
766 1200 : IF (gapw_xc) THEN
767 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_pe(ispin)%matrix, &
768 : rho=trho_xc_r(ispin), &
769 : rho_gspace=trho_xc_g(ispin), &
770 : soft_valid=gapw_xc, &
771 14 : total_rho=rhotot)
772 : END IF
773 : END DO
774 :
775 : ! GAPW o GAPW_XC require the calculation of hard and soft local densities
776 552 : IF (gapw .OR. gapw_xc) THEN
777 78 : CALL get_qs_env(qs_env=qs_env, oce=oce, sab_orb=sab)
778 : CALL calculate_rho_atom_coeff(qs_env, matrix_pe, local_rho_set%rho_atom_set, &
779 78 : qs_kind_set, oce, sab, para_env)
780 78 : CALL prepare_gapw_den(qs_env, local_rho_set, do_rho0=gapw)
781 : END IF
782 1656 : rhotot = SUM(total_rho)
783 552 : IF (gapw) THEN
784 64 : CALL get_rho0_mpole(local_rho_set%rho0_mpole, Qlm_tot=Qlm_tot)
785 64 : rhotot = rhotot + local_rho_set%rho0_mpole%total_rho0_h
786 64 : CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rho_tot_gspace)
787 : END IF
788 :
789 552 : IF (ABS(rhotot) > 1.e-05_dp) THEN
790 20 : logger => cp_get_default_logger()
791 20 : IF (logger%para_env%is_source()) THEN
792 10 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
793 : ELSE
794 : iounit = -1
795 : END IF
796 20 : CPWARN("Real space electron count of excitation density is non-zero.")
797 20 : IF (iounit > 0) THEN
798 10 : WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", rhotot
799 10 : WRITE (iounit, "(T2,A,/)") REPEAT("*", 79)
800 : END IF
801 : END IF
802 :
803 : ! calculate associated hartree potential
804 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace, thartree, &
805 552 : v_hartree_gspace)
806 552 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
807 552 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
808 552 : IF (gapw) THEN
809 : CALL Vh_1c_gg_integrals(qs_env, thartree, hartree_local%ecoul_1c, &
810 64 : local_rho_set, para_env, tddft=.TRUE.)
811 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, &
812 : calculate_forces=.FALSE., &
813 64 : local_rho_set=local_rho_set)
814 : END IF
815 :
816 : ! Fxc*drho term
817 552 : CALL get_qs_env(qs_env, rho=rho)
818 552 : CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
819 : !
820 552 : CALL get_qs_env(qs_env, input=input)
821 552 : IF (dft_control%do_admm) THEN
822 128 : CALL get_qs_env(qs_env, admm_env=admm_env)
823 128 : xc_section => admm_env%xc_section_primary
824 : ELSE
825 424 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
826 : END IF
827 : !
828 552 : deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
829 552 : IF (deriv2_analytic) THEN
830 552 : NULLIFY (v_xc, v_xc_tau, tau_r)
831 552 : IF (gapw_xc) THEN
832 14 : CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
833 14 : CALL qs_fxc_analytic(rho_xc, trho_xc_r, tau_r, xc_section, auxbas_pw_pool, .FALSE., v_xc, v_xc_tau)
834 : ELSE
835 538 : CALL qs_fxc_analytic(rho, trho_r, tau_r, xc_section, auxbas_pw_pool, .FALSE., v_xc, v_xc_tau)
836 : END IF
837 552 : IF (gapw .OR. gapw_xc) THEN
838 78 : CALL get_qs_env(qs_env, rho_atom_set=rho_atom_set)
839 78 : rho1_atom_set => local_rho_set%rho_atom_set
840 : CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
841 78 : do_tddft=.TRUE., do_triplet=.FALSE.)
842 : END IF
843 : ELSE
844 0 : CPABORT("NYA 00006")
845 0 : NULLIFY (v_xc, trho)
846 0 : ALLOCATE (trho)
847 0 : CALL qs_rho_create(trho)
848 0 : CALL qs_rho_set(trho, rho_r=trho_r, rho_g=trho_g)
849 0 : CALL qs_fxc_fdiff(ks_env, rho, trho, xc_section, 6, .FALSE., v_xc, v_xc_tau)
850 0 : DEALLOCATE (trho)
851 : END IF
852 :
853 1200 : DO ispin = 1, nspins
854 648 : CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
855 1200 : CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
856 : END DO
857 552 : IF (gapw_xc) THEN
858 28 : DO ispin = 1, nspins
859 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_hartree_rspace, &
860 : hmat=matrix_hz(ispin), &
861 14 : calculate_forces=.FALSE.)
862 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
863 : hmat=matrix_hz(ispin), &
864 28 : gapw=gapw_xc, calculate_forces=.FALSE.)
865 : END DO
866 : ELSE
867 : ! vtot = v_xc(ispin) + v_hartree
868 1172 : DO ispin = 1, nspins
869 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
870 : hmat=matrix_hz(ispin), &
871 634 : gapw=gapw, calculate_forces=.FALSE.)
872 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_hartree_rspace, &
873 : hmat=matrix_hz(ispin), &
874 1172 : gapw=gapw, calculate_forces=.FALSE.)
875 : END DO
876 : END IF
877 552 : IF (gapw .OR. gapw_xc) THEN
878 78 : mhz(1:nspins, 1:1) => matrix_hz(1:nspins)
879 78 : mpe(1:nspins, 1:1) => matrix_pe(1:nspins)
880 : CALL update_ks_atom(qs_env, mhz, mpe, forces=.FALSE., &
881 78 : rho_atom_external=local_rho_set%rho_atom_set)
882 : END IF
883 :
884 552 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
885 552 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
886 552 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
887 1200 : DO ispin = 1, nspins
888 648 : CALL auxbas_pw_pool%give_back_pw(trho_r(ispin))
889 648 : CALL auxbas_pw_pool%give_back_pw(trho_g(ispin))
890 1200 : CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
891 : END DO
892 552 : DEALLOCATE (trho_r, trho_g, v_xc)
893 552 : IF (gapw_xc) THEN
894 28 : DO ispin = 1, nspins
895 14 : CALL auxbas_pw_pool%give_back_pw(trho_xc_r(ispin))
896 28 : CALL auxbas_pw_pool%give_back_pw(trho_xc_g(ispin))
897 : END DO
898 14 : DEALLOCATE (trho_xc_r, trho_xc_g)
899 : END IF
900 552 : IF (ASSOCIATED(v_xc_tau)) THEN
901 0 : DO ispin = 1, nspins
902 0 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
903 : END DO
904 0 : DEALLOCATE (v_xc_tau)
905 : END IF
906 552 : IF (dft_control%do_admm) THEN
907 128 : IF (qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
908 : ! nothing to do
909 : ELSE
910 : ! add ADMM xc_section_aux terms: f_x[rhoz_ADMM]
911 78 : CALL get_qs_env(qs_env, admm_env=admm_env)
912 : CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit, matrix_s_aux_fit=msaux, &
913 78 : task_list_aux_fit=task_list)
914 78 : basis_type = "AUX_FIT"
915 : !
916 78 : NULLIFY (mpe, mhz)
917 398 : ALLOCATE (mpe(nspins, 1))
918 78 : CALL dbcsr_allocate_matrix_set(mhz, nspins, 1)
919 164 : DO ispin = 1, nspins
920 86 : ALLOCATE (mhz(ispin, 1)%matrix)
921 86 : CALL dbcsr_create(mhz(ispin, 1)%matrix, template=msaux(1)%matrix)
922 86 : CALL dbcsr_copy(mhz(ispin, 1)%matrix, msaux(1)%matrix)
923 86 : CALL dbcsr_set(mhz(ispin, 1)%matrix, 0.0_dp)
924 164 : mpe(ispin, 1)%matrix => matrix_pe_admm(ispin)%matrix
925 : END DO
926 : !
927 : ! GAPW/GAPW_XC initializations
928 78 : NULLIFY (local_rho_set_admm)
929 78 : IF (admm_env%do_gapw) THEN
930 4 : basis_type = "AUX_FIT_SOFT"
931 4 : task_list => admm_env%admm_gapw_env%task_list
932 4 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
933 4 : CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
934 4 : CALL local_rho_set_create(local_rho_set_admm)
935 : CALL allocate_rho_atom_internals(local_rho_set_admm%rho_atom_set, atomic_kind_set, &
936 4 : admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
937 : CALL calculate_rho_atom_coeff(qs_env, matrix_pe_admm, &
938 : rho_atom_set=local_rho_set_admm%rho_atom_set, &
939 : qs_kind_set=admm_env%admm_gapw_env%admm_kind_set, &
940 4 : oce=admm_env%admm_gapw_env%oce, sab=sab_aux_fit, para_env=para_env)
941 : CALL prepare_gapw_den(qs_env, local_rho_set=local_rho_set_admm, &
942 4 : do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
943 : END IF
944 : !
945 78 : xc_section => admm_env%xc_section_aux
946 : !
947 78 : NULLIFY (rho_g_aux, rho_r_aux, rhoz_g_aux, rhoz_r_aux)
948 78 : CALL qs_rho_get(rho_aux_fit, rho_r=rho_r_aux, rho_g=rho_g_aux)
949 : ! rhoz_aux
950 406 : ALLOCATE (rhoz_r_aux(nspins), rhoz_g_aux(nspins))
951 164 : DO ispin = 1, nspins
952 86 : CALL auxbas_pw_pool%create_pw(rhoz_r_aux(ispin))
953 164 : CALL auxbas_pw_pool%create_pw(rhoz_g_aux(ispin))
954 : END DO
955 164 : DO ispin = 1, nspins
956 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpe(ispin, 1)%matrix, &
957 : rho=rhoz_r_aux(ispin), rho_gspace=rhoz_g_aux(ispin), &
958 : basis_type=basis_type, &
959 164 : task_list_external=task_list)
960 : END DO
961 : !
962 78 : NULLIFY (v_xc)
963 78 : deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
964 78 : IF (deriv2_analytic) THEN
965 78 : NULLIFY (tau_r)
966 78 : CALL qs_fxc_analytic(rho_aux_fit, rhoz_r_aux, tau_r, xc_section, auxbas_pw_pool, .FALSE., v_xc, v_xc_tau)
967 : ELSE
968 0 : CPABORT("NYA 00007")
969 : NULLIFY (rhoz_aux)
970 0 : ALLOCATE (rhoz_aux)
971 0 : CALL qs_rho_create(rhoz_aux)
972 0 : CALL qs_rho_set(rhoz_aux, rho_r=rhoz_r_aux, rho_g=rhoz_g_aux)
973 0 : CALL qs_fxc_fdiff(ks_env, rho_aux_fit, rhoz_aux, xc_section, 6, .FALSE., v_xc, v_xc_tau)
974 0 : DEALLOCATE (rhoz_aux)
975 : END IF
976 : !
977 164 : DO ispin = 1, nspins
978 86 : CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
979 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
980 : hmat=mhz(ispin, 1), basis_type=basis_type, &
981 : calculate_forces=.FALSE., &
982 164 : task_list_external=task_list)
983 : END DO
984 164 : DO ispin = 1, nspins
985 86 : CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
986 86 : CALL auxbas_pw_pool%give_back_pw(rhoz_r_aux(ispin))
987 164 : CALL auxbas_pw_pool%give_back_pw(rhoz_g_aux(ispin))
988 : END DO
989 78 : DEALLOCATE (v_xc, rhoz_r_aux, rhoz_g_aux)
990 : !
991 78 : IF (admm_env%do_gapw) THEN
992 4 : rho_atom_set => admm_env%admm_gapw_env%local_rho_set%rho_atom_set
993 4 : rho1_atom_set => local_rho_set_admm%rho_atom_set
994 : CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, &
995 4 : para_env, kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
996 : CALL update_ks_atom(qs_env, mhz(:, 1), matrix_pe_admm, forces=.FALSE., tddft=.FALSE., &
997 : rho_atom_external=rho1_atom_set, &
998 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
999 : oce_external=admm_env%admm_gapw_env%oce, &
1000 4 : sab_external=sab_aux_fit)
1001 : END IF
1002 : !
1003 78 : nao = admm_env%nao_orb
1004 78 : nao_aux = admm_env%nao_aux_fit
1005 78 : ALLOCATE (dbwork)
1006 78 : CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
1007 164 : DO ispin = 1, nspins
1008 : CALL cp_dbcsr_sm_fm_multiply(mhz(ispin, 1)%matrix, admm_env%A, &
1009 86 : admm_env%work_aux_orb, nao)
1010 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
1011 : 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
1012 86 : admm_env%work_orb_orb)
1013 86 : CALL dbcsr_copy(dbwork, matrix_hz(1)%matrix)
1014 86 : CALL dbcsr_set(dbwork, 0.0_dp)
1015 86 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.TRUE.)
1016 164 : CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
1017 : END DO
1018 78 : CALL dbcsr_release(dbwork)
1019 78 : DEALLOCATE (dbwork)
1020 78 : CALL dbcsr_deallocate_matrix_set(mhz)
1021 78 : DEALLOCATE (mpe)
1022 78 : IF (admm_env%do_gapw) THEN
1023 4 : IF (ASSOCIATED(local_rho_set_admm)) CALL local_rho_set_release(local_rho_set_admm)
1024 : END IF
1025 : END IF
1026 : END IF
1027 552 : IF (gapw .OR. gapw_xc) THEN
1028 78 : IF (ASSOCIATED(local_rho_set)) CALL local_rho_set_release(local_rho_set)
1029 78 : IF (ASSOCIATED(hartree_local)) CALL hartree_local_release(hartree_local)
1030 : END IF
1031 :
1032 : ! HFX
1033 552 : hfx_section => section_vals_get_subs_vals(xc_section, "HF")
1034 552 : CALL section_vals_get(hfx_section, explicit=do_hfx)
1035 552 : IF (do_hfx) THEN
1036 248 : CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
1037 248 : CPASSERT(n_rep_hf == 1)
1038 : CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
1039 248 : i_rep_section=1)
1040 248 : mspin = 1
1041 248 : IF (hfx_treat_lsd_in_core) mspin = nspins
1042 : !
1043 : CALL get_qs_env(qs_env=qs_env, rho=rho, x_data=x_data, para_env=para_env, &
1044 248 : s_mstruct_changed=s_mstruct_changed)
1045 248 : distribute_fock_matrix = .TRUE.
1046 248 : IF (dft_control%do_admm) THEN
1047 128 : CALL get_qs_env(qs_env, admm_env=admm_env)
1048 128 : CALL get_admm_env(admm_env, matrix_s_aux_fit=msaux)
1049 128 : NULLIFY (mpe, mhz)
1050 660 : ALLOCATE (mpe(nspins, 1))
1051 128 : CALL dbcsr_allocate_matrix_set(mhz, nspins, 1)
1052 276 : DO ispin = 1, nspins
1053 148 : ALLOCATE (mhz(ispin, 1)%matrix)
1054 148 : CALL dbcsr_create(mhz(ispin, 1)%matrix, template=msaux(1)%matrix)
1055 148 : CALL dbcsr_copy(mhz(ispin, 1)%matrix, msaux(1)%matrix)
1056 148 : CALL dbcsr_set(mhz(ispin, 1)%matrix, 0.0_dp)
1057 276 : mpe(ispin, 1)%matrix => matrix_pe_admm(ispin)%matrix
1058 : END DO
1059 128 : IF (x_data(1, 1)%do_hfx_ri) THEN
1060 : eh1 = 0.0_dp
1061 : CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpe, &
1062 : geometry_did_change=s_mstruct_changed, nspins=nspins, &
1063 6 : hf_fraction=x_data(1, 1)%general_parameter%fraction)
1064 : ELSE
1065 244 : DO ispin = 1, mspin
1066 : eh1 = 0.0
1067 : CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpe, hfx_section, &
1068 : para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
1069 244 : ispin=ispin)
1070 : END DO
1071 : END IF
1072 : !
1073 128 : CPASSERT(ASSOCIATED(admm_env%work_aux_orb))
1074 128 : CPASSERT(ASSOCIATED(admm_env%work_orb_orb))
1075 128 : nao = admm_env%nao_orb
1076 128 : nao_aux = admm_env%nao_aux_fit
1077 128 : ALLOCATE (dbwork)
1078 128 : CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
1079 276 : DO ispin = 1, nspins
1080 : CALL cp_dbcsr_sm_fm_multiply(mhz(ispin, 1)%matrix, admm_env%A, &
1081 148 : admm_env%work_aux_orb, nao)
1082 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
1083 : 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
1084 148 : admm_env%work_orb_orb)
1085 148 : CALL dbcsr_copy(dbwork, matrix_hz(ispin)%matrix)
1086 148 : CALL dbcsr_set(dbwork, 0.0_dp)
1087 148 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.TRUE.)
1088 276 : CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
1089 : END DO
1090 128 : CALL dbcsr_release(dbwork)
1091 128 : DEALLOCATE (dbwork)
1092 128 : CALL dbcsr_deallocate_matrix_set(mhz)
1093 128 : DEALLOCATE (mpe)
1094 : ELSE
1095 120 : NULLIFY (mpe, mhz)
1096 1000 : ALLOCATE (mpe(nspins, 1), mhz(nspins, 1))
1097 260 : DO ispin = 1, nspins
1098 140 : mhz(ispin, 1)%matrix => matrix_hz(ispin)%matrix
1099 260 : mpe(ispin, 1)%matrix => matrix_pe(ispin)%matrix
1100 : END DO
1101 120 : IF (x_data(1, 1)%do_hfx_ri) THEN
1102 : eh1 = 0.0_dp
1103 : CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpe, &
1104 : geometry_did_change=s_mstruct_changed, nspins=nspins, &
1105 18 : hf_fraction=x_data(1, 1)%general_parameter%fraction)
1106 : ELSE
1107 204 : DO ispin = 1, mspin
1108 : eh1 = 0.0
1109 : CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpe, hfx_section, &
1110 : para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
1111 204 : ispin=ispin)
1112 : END DO
1113 : END IF
1114 120 : DEALLOCATE (mpe, mhz)
1115 : END IF
1116 : END IF
1117 :
1118 552 : focc = 4.0_dp
1119 552 : IF (nspins == 2) focc = 2.0_dp
1120 1200 : DO ispin = 1, nspins
1121 648 : mos => gs_mos(ispin)%mos_occ
1122 648 : CALL cp_fm_get_info(mos, ncol_global=norb)
1123 : CALL cp_dbcsr_sm_fm_multiply(matrix_hz(ispin)%matrix, mos, cpmos(ispin), &
1124 1200 : norb, alpha=focc, beta=0.0_dp)
1125 : END DO
1126 :
1127 552 : CALL timestop(handle)
1128 :
1129 1656 : END SUBROUTINE tddfpt_resvec2
1130 :
1131 : ! **************************************************************************************************
1132 : !> \brief ...
1133 : !> \param qs_env ...
1134 : !> \param matrix_pe ...
1135 : !> \param gs_mos ...
1136 : !> \param matrix_hz ...
1137 : !> \param cpmos ...
1138 : ! **************************************************************************************************
1139 16 : SUBROUTINE tddfpt_resvec2_xtb(qs_env, matrix_pe, gs_mos, matrix_hz, cpmos)
1140 :
1141 : TYPE(qs_environment_type), POINTER :: qs_env
1142 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pe
1143 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1144 : POINTER :: gs_mos
1145 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz
1146 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: cpmos
1147 :
1148 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_resvec2_xtb'
1149 :
1150 : INTEGER :: atom_a, handle, iatom, ikind, is, ispin, &
1151 : na, natom, natorb, nkind, norb, ns, &
1152 : nsgf, nspins
1153 : INTEGER, DIMENSION(25) :: lao
1154 : INTEGER, DIMENSION(5) :: occ
1155 16 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: mcharge, mcharge1
1156 16 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, charges, charges1
1157 : REAL(KIND=dp) :: focc
1158 16 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1159 : TYPE(cp_fm_type), POINTER :: mos
1160 16 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
1161 16 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
1162 : TYPE(dbcsr_type), POINTER :: s_matrix
1163 : TYPE(dft_control_type), POINTER :: dft_control
1164 : TYPE(mp_para_env_type), POINTER :: para_env
1165 16 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1166 16 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1167 : TYPE(qs_rho_type), POINTER :: rho
1168 : TYPE(xtb_atom_type), POINTER :: xtb_kind
1169 :
1170 16 : CALL timeset(routineN, handle)
1171 :
1172 16 : CPASSERT(ASSOCIATED(matrix_pe))
1173 :
1174 16 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
1175 16 : nspins = dft_control%nspins
1176 :
1177 32 : DO ispin = 1, nspins
1178 32 : CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
1179 : END DO
1180 :
1181 16 : IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
1182 : ! Mulliken charges
1183 : CALL get_qs_env(qs_env, rho=rho, particle_set=particle_set, &
1184 14 : matrix_s_kp=matrix_s, para_env=para_env)
1185 14 : natom = SIZE(particle_set)
1186 14 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1187 70 : ALLOCATE (mcharge(natom), charges(natom, 5))
1188 42 : ALLOCATE (mcharge1(natom), charges1(natom, 5))
1189 1254 : charges = 0.0_dp
1190 1254 : charges1 = 0.0_dp
1191 14 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
1192 14 : nkind = SIZE(atomic_kind_set)
1193 14 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
1194 56 : ALLOCATE (aocg(nsgf, natom))
1195 1184 : aocg = 0.0_dp
1196 42 : ALLOCATE (aocg1(nsgf, natom))
1197 1184 : aocg1 = 0.0_dp
1198 14 : p_matrix => matrix_p(:, 1)
1199 14 : s_matrix => matrix_s(1, 1)%matrix
1200 14 : CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
1201 14 : CALL ao_charges(matrix_pe, s_matrix, aocg1, para_env)
1202 48 : DO ikind = 1, nkind
1203 34 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
1204 34 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1205 34 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
1206 316 : DO iatom = 1, na
1207 234 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1208 1404 : charges(atom_a, :) = REAL(occ(:), KIND=dp)
1209 900 : DO is = 1, natorb
1210 632 : ns = lao(is) + 1
1211 632 : charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
1212 866 : charges1(atom_a, ns) = charges1(atom_a, ns) - aocg1(is, atom_a)
1213 : END DO
1214 : END DO
1215 : END DO
1216 14 : DEALLOCATE (aocg, aocg1)
1217 248 : DO iatom = 1, natom
1218 1404 : mcharge(iatom) = SUM(charges(iatom, :))
1219 1418 : mcharge1(iatom) = SUM(charges1(iatom, :))
1220 : END DO
1221 : ! Coulomb Kernel
1222 14 : CALL xtb_coulomb_hessian(qs_env, matrix_hz, charges1, mcharge1, mcharge)
1223 : !
1224 28 : DEALLOCATE (charges, mcharge, charges1, mcharge1)
1225 : END IF
1226 :
1227 16 : focc = 2.0_dp
1228 16 : IF (nspins == 2) focc = 1.0_dp
1229 32 : DO ispin = 1, nspins
1230 16 : mos => gs_mos(ispin)%mos_occ
1231 16 : CALL cp_fm_get_info(mos, ncol_global=norb)
1232 : CALL cp_dbcsr_sm_fm_multiply(matrix_hz(ispin)%matrix, mos, cpmos(ispin), &
1233 32 : norb, alpha=focc, beta=0.0_dp)
1234 : END DO
1235 :
1236 16 : CALL timestop(handle)
1237 :
1238 32 : END SUBROUTINE tddfpt_resvec2_xtb
1239 :
1240 : ! **************************************************************************************************
1241 : !> \brief ...
1242 : !> \param qs_env ...
1243 : !> \param cpmos ...
1244 : !> \param work ...
1245 : ! **************************************************************************************************
1246 568 : SUBROUTINE tddfpt_resvec3(qs_env, cpmos, work)
1247 :
1248 : TYPE(qs_environment_type), POINTER :: qs_env
1249 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: cpmos
1250 : TYPE(tddfpt_work_matrices) :: work
1251 :
1252 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_resvec3'
1253 :
1254 : INTEGER :: handle, ispin, nao, norb, nspins
1255 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
1256 : TYPE(cp_fm_type) :: cvec, umat
1257 : TYPE(cp_fm_type), POINTER :: omos
1258 : TYPE(dft_control_type), POINTER :: dft_control
1259 568 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1260 :
1261 568 : CALL timeset(routineN, handle)
1262 :
1263 568 : CALL get_qs_env(qs_env, mos=mos, dft_control=dft_control)
1264 568 : nspins = dft_control%nspins
1265 :
1266 1232 : DO ispin = 1, nspins
1267 664 : CALL get_mo_set(mos(ispin), mo_coeff=omos)
1268 : ASSOCIATE (rvecs => cpmos(ispin))
1269 664 : CALL cp_fm_get_info(rvecs, nrow_global=nao, ncol_global=norb)
1270 664 : CALL cp_fm_create(cvec, rvecs%matrix_struct, "cvec")
1271 : CALL cp_fm_struct_create(fmstruct, context=rvecs%matrix_struct%context, nrow_global=norb, &
1272 664 : ncol_global=norb, para_env=rvecs%matrix_struct%para_env)
1273 664 : CALL cp_fm_create(umat, fmstruct, "umat")
1274 664 : CALL cp_fm_struct_release(fmstruct)
1275 : !
1276 664 : CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, omos, work%S_C0(ispin), 0.0_dp, umat)
1277 664 : CALL cp_fm_copy_general(rvecs, cvec, rvecs%matrix_struct%para_env)
1278 664 : CALL parallel_gemm("N", "T", nao, norb, norb, 1.0_dp, cvec, umat, 0.0_dp, rvecs)
1279 : END ASSOCIATE
1280 664 : CALL cp_fm_release(cvec)
1281 2560 : CALL cp_fm_release(umat)
1282 : END DO
1283 :
1284 568 : CALL timestop(handle)
1285 :
1286 568 : END SUBROUTINE tddfpt_resvec3
1287 :
1288 : ! **************************************************************************************************
1289 : !> \brief Calculate direct tddft forces
1290 : !> \param qs_env ...
1291 : !> \param ex_env ...
1292 : !> \param gs_mos ...
1293 : !> \param kernel_env ...
1294 : !> \param sub_env ...
1295 : !> \param work_matrices ...
1296 : !> \param debug_forces ...
1297 : !> \par History
1298 : !> * 01.2020 screated [JGH]
1299 : ! **************************************************************************************************
1300 568 : SUBROUTINE tddfpt_kernel_force(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, debug_forces)
1301 :
1302 : TYPE(qs_environment_type), POINTER :: qs_env
1303 : TYPE(excited_energy_type), POINTER :: ex_env
1304 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1305 : POINTER :: gs_mos
1306 : TYPE(kernel_env_type), INTENT(IN) :: kernel_env
1307 : TYPE(tddfpt_subgroup_env_type) :: sub_env
1308 : TYPE(tddfpt_work_matrices) :: work_matrices
1309 : LOGICAL, INTENT(IN) :: debug_forces
1310 :
1311 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_kernel_force'
1312 :
1313 : INTEGER :: handle
1314 : TYPE(dft_control_type), POINTER :: dft_control
1315 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
1316 :
1317 : MARK_USED(work_matrices)
1318 :
1319 568 : CALL timeset(routineN, handle)
1320 :
1321 568 : CALL get_qs_env(qs_env, dft_control=dft_control)
1322 568 : tddfpt_control => dft_control%tddfpt2_control
1323 :
1324 568 : IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
1325 : ! full Kernel
1326 342 : CALL fhxc_force(qs_env, ex_env, gs_mos, kernel_env%full_kernel, debug_forces)
1327 226 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
1328 : ! sTDA Kernel
1329 158 : CALL stda_force(qs_env, ex_env, gs_mos, kernel_env%stda_kernel, sub_env, work_matrices, debug_forces)
1330 68 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
1331 : ! nothing to be done here
1332 68 : ex_env%matrix_wx1 => NULL()
1333 : ELSE
1334 0 : CPABORT('Unknown kernel type')
1335 : END IF
1336 :
1337 568 : CALL timestop(handle)
1338 :
1339 568 : END SUBROUTINE tddfpt_kernel_force
1340 :
1341 : END MODULE qs_tddfpt2_forces
|
