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