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