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