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_fhxc_forces
9 : USE accint_weights_forces, ONLY: accint_weight_force
10 : USE admm_methods, ONLY: admm_projection_derivative
11 : USE admm_types, ONLY: admm_type,&
12 : get_admm_env
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind_set
15 : USE cell_types, ONLY: cell_type,&
16 : pbc
17 : USE cp_control_types, ONLY: dft_control_type,&
18 : stda_control_type,&
19 : tddfpt2_control_type
20 : USE cp_dbcsr_api, ONLY: &
21 : dbcsr_add, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_filter, &
22 : dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
23 : dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
24 : dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type, &
25 : dbcsr_type_antisymmetric, dbcsr_type_no_symmetry
26 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
27 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 : copy_fm_to_dbcsr,&
29 : cp_dbcsr_plus_fm_fm_t,&
30 : cp_dbcsr_sm_fm_multiply,&
31 : dbcsr_allocate_matrix_set,&
32 : dbcsr_deallocate_matrix_set
33 : USE cp_fm_basic_linalg, ONLY: cp_fm_add_columns,&
34 : cp_fm_geadd,&
35 : cp_fm_row_scale,&
36 : cp_fm_schur_product
37 : USE cp_fm_pool_types, ONLY: fm_pool_create_fm,&
38 : fm_pool_give_back_fm
39 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
40 : cp_fm_struct_release,&
41 : cp_fm_struct_type
42 : USE cp_fm_types, ONLY: cp_fm_create,&
43 : cp_fm_get_info,&
44 : cp_fm_release,&
45 : cp_fm_to_fm,&
46 : cp_fm_type,&
47 : cp_fm_vectorssum
48 : USE cp_log_handling, ONLY: cp_get_default_logger,&
49 : cp_logger_get_default_unit_nr,&
50 : cp_logger_type
51 : USE ewald_environment_types, ONLY: ewald_env_get,&
52 : ewald_environment_type
53 : USE ewald_methods_tb, ONLY: tb_ewald_overlap,&
54 : tb_spme_evaluate
55 : USE ewald_pw_types, ONLY: ewald_pw_type
56 : USE exstates_types, ONLY: excited_energy_type
57 : USE hartree_local_methods, ONLY: Vh_1c_gg_integrals,&
58 : init_coulomb_local
59 : USE hartree_local_types, ONLY: hartree_local_create,&
60 : hartree_local_release,&
61 : hartree_local_type
62 : USE hfx_derivatives, ONLY: derivatives_four_center
63 : USE hfx_energy_potential, ONLY: integrate_four_center
64 : USE hfx_ri, ONLY: hfx_ri_update_forces,&
65 : hfx_ri_update_ks
66 : USE hfx_types, ONLY: hfx_type
67 : USE input_constants, ONLY: do_admm_aux_exch_func_none,&
68 : no_sf_tddfpt,&
69 : tddfpt_kernel_full,&
70 : tddfpt_sf_col,&
71 : xc_kernel_method_analytic,&
72 : xc_kernel_method_best,&
73 : xc_kernel_method_numeric,&
74 : xc_none
75 : USE input_section_types, ONLY: section_vals_get,&
76 : section_vals_get_subs_vals,&
77 : section_vals_type,&
78 : section_vals_val_get
79 : USE kinds, ONLY: default_string_length,&
80 : dp
81 : USE mathconstants, ONLY: oorootpi
82 : USE message_passing, ONLY: mp_para_env_type
83 : USE parallel_gemm_api, ONLY: parallel_gemm
84 : USE particle_methods, ONLY: get_particle_set
85 : USE particle_types, ONLY: particle_type
86 : USE pw_env_types, ONLY: pw_env_get,&
87 : pw_env_type
88 : USE pw_methods, ONLY: pw_axpy,&
89 : pw_scale,&
90 : pw_transfer,&
91 : pw_zero
92 : USE pw_poisson_methods, ONLY: pw_poisson_solve
93 : USE pw_poisson_types, ONLY: pw_poisson_type
94 : USE pw_pool_types, ONLY: pw_pool_type
95 : USE pw_types, ONLY: pw_c1d_gs_type,&
96 : pw_r3d_rs_type
97 : USE qs_collocate_density, ONLY: calculate_rho_elec
98 : USE qs_environment_types, ONLY: get_qs_env,&
99 : qs_environment_type,&
100 : set_qs_env
101 : USE qs_force_types, ONLY: qs_force_type
102 : USE qs_fxc, ONLY: qs_fgxc_analytic,&
103 : qs_fgxc_create,&
104 : qs_fgxc_gdiff,&
105 : qs_fgxc_release
106 : USE qs_gapw_densities, ONLY: prepare_gapw_den
107 : USE qs_integrate_potential, ONLY: integrate_v_rspace
108 : USE qs_kernel_types, ONLY: full_kernel_env_type
109 : USE qs_kind_types, ONLY: qs_kind_type
110 : USE qs_ks_atom, ONLY: update_ks_atom
111 : USE qs_ks_types, ONLY: qs_ks_env_type
112 : USE qs_local_rho_types, ONLY: local_rho_set_create,&
113 : local_rho_set_release,&
114 : local_rho_type
115 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
116 : USE qs_oce_methods, ONLY: build_oce_matrices
117 : USE qs_oce_types, ONLY: allocate_oce_set,&
118 : create_oce_set,&
119 : oce_matrix_type
120 : USE qs_overlap, ONLY: build_overlap_matrix
121 : USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace,&
122 : rho0_s_grid_create
123 : USE qs_rho0_methods, ONLY: init_rho0
124 : USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals,&
125 : calculate_rho_atom_coeff
126 : USE qs_rho_atom_types, ONLY: rho_atom_type
127 : USE qs_rho_types, ONLY: qs_rho_create,&
128 : qs_rho_get,&
129 : qs_rho_set,&
130 : qs_rho_type
131 : USE qs_tddfpt2_stda_types, ONLY: stda_env_type
132 : USE qs_tddfpt2_stda_utils, ONLY: get_lowdin_x,&
133 : setup_gamma
134 : USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
135 : USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos,&
136 : tddfpt_work_matrices
137 : USE qs_vxc_atom, ONLY: calculate_gfxc_atom,&
138 : gfxc_atom_diff
139 : USE task_list_types, ONLY: task_list_type
140 : USE util, ONLY: get_limit
141 : USE virial_types, ONLY: virial_type
142 : USE xc_derivatives, ONLY: xc_functionals_get_needs
143 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
144 : #include "./base/base_uses.f90"
145 :
146 : IMPLICIT NONE
147 :
148 : PRIVATE
149 :
150 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_fhxc_forces'
151 :
152 : PUBLIC :: fhxc_force, stda_force
153 :
154 : ! **************************************************************************************************
155 :
156 : CONTAINS
157 :
158 : ! **************************************************************************************************
159 : !> \brief Calculate direct tddft forces. Calculate the three last terms of the response vector
160 : !> in equation 49 and the first term of \Lambda_munu in equation 51 in
161 : !> J. Chem. Theory Comput. 2022, 18, 7, 4186–4202 (https://doi.org/10.1021/acs.jctc.2c00144)
162 : !> \param qs_env Holds all system information relevant for the calculation.
163 : !> \param ex_env Holds the response vector ex_env%cpmos and Lambda ex_env%matrix_wx1.
164 : !> \param gs_mos MO coefficients of the ground state.
165 : !> \param full_kernel ...
166 : !> \param debug_forces ...
167 : !> \par History
168 : !> * 01.2020 screated [JGH]
169 : ! **************************************************************************************************
170 420 : SUBROUTINE fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
171 :
172 : TYPE(qs_environment_type), POINTER :: qs_env
173 : TYPE(excited_energy_type), POINTER :: ex_env
174 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
175 : POINTER :: gs_mos
176 : TYPE(full_kernel_env_type), INTENT(IN) :: full_kernel
177 : LOGICAL, INTENT(IN) :: debug_forces
178 :
179 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fhxc_force'
180 :
181 : CHARACTER(LEN=default_string_length) :: basis_type
182 : INTEGER :: handle, ia, ib, iounit, ispin, mspin, &
183 : myfun, n_rep_hf, nactive(2), nao, &
184 : nao_aux, natom, nkind, norb(2), nsev, &
185 : nspins, order, spin
186 : LOGICAL :: distribute_fock_matrix, do_admm, do_analytic, do_hfx, do_hfxlr, do_hfxsr, &
187 : do_numeric, do_res, do_sf, gapw, gapw_xc, hfx_treat_lsd_in_core, is_rks_triplets, &
188 : needs_tau_response, needs_tau_response_aux, s_mstruct_changed, use_virial
189 : REAL(KIND=dp) :: eh1, eh1c, eps_delta, eps_fit, focc, &
190 : fscal, fval, kval, xehartree
191 : REAL(KIND=dp), DIMENSION(3) :: fodeb
192 : TYPE(admm_type), POINTER :: admm_env
193 420 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
194 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
195 : TYPE(cp_fm_type) :: avamat, avcmat, cpscr, cvcmat, vavec, &
196 : vcvec
197 420 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos, evect
198 : TYPE(cp_fm_type), POINTER :: mos, mos2, mosa, mosa2
199 : TYPE(cp_logger_type), POINTER :: logger
200 420 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_fx, matrix_gx, matrix_hfx, &
201 420 : matrix_hfx_admm, matrix_hfx_admm_asymm, matrix_hfx_asymm, matrix_hx, matrix_p, &
202 420 : matrix_p_admm, matrix_px1, matrix_px1_admm, matrix_px1_admm_asymm, matrix_px1_asymm, &
203 420 : matrix_s, matrix_s_aux_fit, matrix_wx1, mdum, mfx, mgx
204 420 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mhe, mpe, mpga
205 : TYPE(dbcsr_type), POINTER :: dbwork, dbwork_asymm
206 : TYPE(dft_control_type), POINTER :: dft_control
207 : TYPE(hartree_local_type), POINTER :: hartree_local
208 420 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
209 : TYPE(local_rho_type), POINTER :: local_rho_set, local_rho_set_admm, local_rho_set_f, &
210 : local_rho_set_f_admm, local_rho_set_g, local_rho_set_g_admm
211 : TYPE(mp_para_env_type), POINTER :: para_env
212 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
213 420 : POINTER :: sab, sab_aux_fit, sab_orb, sap_oce
214 : TYPE(oce_matrix_type), POINTER :: oce
215 420 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
216 : TYPE(pw_c1d_gs_type) :: rhox_tot_gspace, xv_hartree_gspace
217 420 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux, rhox_g, rhox_g_aux, &
218 420 : rhox_tau_g, rhox_tau_g_aux, rhoxx_g, &
219 420 : rhoxx_tau_g
220 : TYPE(pw_env_type), POINTER :: pw_env
221 : TYPE(pw_poisson_type), POINTER :: poisson_env
222 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
223 : TYPE(pw_r3d_rs_type) :: xv_hartree_rspace
224 420 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rho, fxc_tau, gxc_rho, gxc_tau, &
225 420 : rho_r_aux, rhox_r, rhox_r_aux, &
226 420 : rhox_tau_r, rhox_tau_r_aux, rhoxx_r, &
227 420 : rhoxx_tau_r
228 : TYPE(pw_r3d_rs_type), POINTER :: weights
229 420 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
230 420 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
231 : TYPE(qs_ks_env_type), POINTER :: ks_env
232 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit, rhox, rhox_aux
233 420 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho_atom_set_f, &
234 420 : rho_atom_set_g
235 : TYPE(section_vals_type), POINTER :: hfx_section, xc_fun_section, xc_section
236 : TYPE(task_list_type), POINTER :: task_list
237 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
238 : TYPE(xc_rho_cflags_type) :: needs
239 :
240 420 : CALL timeset(routineN, handle)
241 :
242 420 : logger => cp_get_default_logger()
243 420 : IF (logger%para_env%is_source()) THEN
244 210 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
245 : ELSE
246 : iounit = -1
247 : END IF
248 :
249 420 : CALL get_qs_env(qs_env, dft_control=dft_control)
250 420 : tddfpt_control => dft_control%tddfpt2_control
251 420 : nspins = dft_control%nspins
252 420 : is_rks_triplets = tddfpt_control%rks_triplets .AND. (nspins == 1)
253 420 : IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
254 408 : do_sf = .FALSE.
255 : ELSE
256 12 : do_sf = .TRUE.
257 : END IF
258 420 : CPASSERT(tddfpt_control%kernel == tddfpt_kernel_full)
259 420 : do_hfx = tddfpt_control%do_hfx
260 420 : do_hfxsr = tddfpt_control%do_hfxsr
261 420 : do_hfxlr = tddfpt_control%do_hfxlr
262 420 : do_admm = tddfpt_control%do_admm
263 420 : gapw = dft_control%qs_control%gapw
264 420 : gapw_xc = dft_control%qs_control%gapw_xc
265 420 : xc_section => full_kernel%xc_section
266 420 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
267 420 : needs = xc_functionals_get_needs(xc_fun_section, (nspins == 2), .TRUE.)
268 420 : needs_tau_response = needs%tau .OR. needs%tau_spin
269 420 : needs_tau_response_aux = .FALSE.
270 :
271 420 : evect => ex_env%evect
272 420 : matrix_px1 => ex_env%matrix_px1
273 420 : matrix_px1_admm => ex_env%matrix_px1_admm
274 420 : matrix_px1_asymm => ex_env%matrix_px1_asymm
275 420 : matrix_px1_admm_asymm => ex_env%matrix_px1_admm_asymm
276 :
277 420 : focc = 1.0_dp
278 420 : IF (nspins == 2) focc = 0.5_dp
279 420 : nsev = SIZE(evect, 1)
280 910 : DO ispin = 1, nsev
281 490 : CALL cp_fm_get_info(evect(ispin), ncol_global=nactive(ispin))
282 : ! Calculate (C*X^T + X*C^T)/2
283 490 : CALL dbcsr_set(matrix_px1(ispin)%matrix, 0.0_dp)
284 : CALL cp_dbcsr_plus_fm_fm_t(matrix_px1(ispin)%matrix, &
285 : matrix_v=evect(ispin), &
286 : matrix_g=gs_mos(ispin)%mos_active, &
287 490 : ncol=nactive(ispin), alpha=2.0_dp*focc, symmetry_mode=1)
288 :
289 : ! Calculate (C*X^T - X*C^T)/2
290 490 : CALL dbcsr_set(matrix_px1_asymm(ispin)%matrix, 0.0_dp)
291 : CALL cp_dbcsr_plus_fm_fm_t(matrix_px1_asymm(ispin)%matrix, &
292 : matrix_v=gs_mos(ispin)%mos_active, &
293 : matrix_g=evect(ispin), &
294 : ncol=nactive(ispin), alpha=2.0_dp*focc, &
295 910 : symmetry_mode=-1)
296 : END DO
297 : !
298 420 : CALL get_qs_env(qs_env, ks_env=ks_env, pw_env=pw_env, para_env=para_env)
299 420 : CALL get_qs_env(qs_env, xcint_weights=weights)
300 : !
301 420 : NULLIFY (hartree_local, local_rho_set, local_rho_set_admm)
302 420 : IF (gapw .OR. gapw_xc) THEN
303 124 : IF (nspins == 2) THEN
304 0 : DO ispin = 1, nsev
305 0 : CALL dbcsr_scale(matrix_px1(ispin)%matrix, 2.0_dp)
306 : END DO
307 : END IF
308 : CALL get_qs_env(qs_env, &
309 : atomic_kind_set=atomic_kind_set, &
310 124 : qs_kind_set=qs_kind_set)
311 124 : CALL local_rho_set_create(local_rho_set)
312 : CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
313 124 : qs_kind_set, dft_control, para_env)
314 124 : IF (gapw) THEN
315 104 : CALL get_qs_env(qs_env, natom=natom)
316 : CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control, &
317 104 : zcore=0.0_dp)
318 104 : CALL rho0_s_grid_create(pw_env, local_rho_set%rho0_mpole)
319 104 : CALL hartree_local_create(hartree_local)
320 104 : CALL init_coulomb_local(hartree_local, natom)
321 : END IF
322 :
323 124 : CALL get_qs_env(qs_env=qs_env, oce=oce, sap_oce=sap_oce, sab_orb=sab)
324 124 : CALL create_oce_set(oce)
325 124 : CALL get_qs_env(qs_env=qs_env, nkind=nkind, particle_set=particle_set)
326 124 : CALL allocate_oce_set(oce, nkind)
327 124 : eps_fit = dft_control%qs_control%gapw_control%eps_fit
328 : CALL build_oce_matrices(oce%intac, .TRUE., 1, qs_kind_set, particle_set, &
329 124 : sap_oce, eps_fit)
330 124 : CALL set_qs_env(qs_env, oce=oce)
331 :
332 124 : mpga(1:nsev, 1:1) => matrix_px1(1:nsev)
333 : CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set%rho_atom_set, &
334 124 : qs_kind_set, oce, sab, para_env)
335 124 : CALL prepare_gapw_den(qs_env, local_rho_set, do_rho0=gapw)
336 : !
337 124 : CALL local_rho_set_create(local_rho_set_f)
338 : CALL allocate_rho_atom_internals(local_rho_set_f%rho_atom_set, atomic_kind_set, &
339 124 : qs_kind_set, dft_control, para_env)
340 : CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set_f%rho_atom_set, &
341 124 : qs_kind_set, oce, sab, para_env)
342 124 : CALL prepare_gapw_den(qs_env, local_rho_set_f, do_rho0=.FALSE.)
343 : !
344 124 : CALL local_rho_set_create(local_rho_set_g)
345 : CALL allocate_rho_atom_internals(local_rho_set_g%rho_atom_set, atomic_kind_set, &
346 124 : qs_kind_set, dft_control, para_env)
347 : CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set_g%rho_atom_set, &
348 124 : qs_kind_set, oce, sab, para_env)
349 124 : CALL prepare_gapw_den(qs_env, local_rho_set_g, do_rho0=.FALSE.)
350 124 : IF (nspins == 2) THEN
351 0 : DO ispin = 1, nsev
352 0 : CALL dbcsr_scale(matrix_px1(ispin)%matrix, 0.5_dp)
353 : END DO
354 : END IF
355 : END IF
356 : !
357 420 : IF (do_admm) THEN
358 78 : CALL get_qs_env(qs_env, admm_env=admm_env)
359 78 : nao_aux = admm_env%nao_aux_fit
360 78 : nao = admm_env%nao_orb
361 : ! Fit the symmetrized and antisymmetrized matrices
362 160 : DO ispin = 1, nsev
363 :
364 82 : CALL copy_dbcsr_to_fm(matrix_px1(ispin)%matrix, admm_env%work_orb_orb)
365 : CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
366 : 1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
367 82 : admm_env%work_aux_orb)
368 : CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
369 : 1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
370 82 : admm_env%work_aux_aux)
371 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, matrix_px1_admm(ispin)%matrix, &
372 82 : keep_sparsity=.TRUE.)
373 :
374 82 : CALL copy_dbcsr_to_fm(matrix_px1_asymm(ispin)%matrix, admm_env%work_orb_orb)
375 : CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
376 : 1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
377 82 : admm_env%work_aux_orb)
378 : CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
379 : 1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
380 82 : admm_env%work_aux_aux)
381 : CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, matrix_px1_admm_asymm(ispin)%matrix, &
382 160 : keep_sparsity=.TRUE.)
383 : END DO
384 : !
385 78 : IF (admm_env%do_gapw) THEN
386 24 : IF (do_admm .AND. tddfpt_control%admm_xc_correction) THEN
387 18 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
388 : ! nothing to do
389 : ELSE
390 6 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
391 6 : CALL local_rho_set_create(local_rho_set_admm)
392 : CALL allocate_rho_atom_internals(local_rho_set_admm%rho_atom_set, atomic_kind_set, &
393 6 : admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
394 6 : mpga(1:nsev, 1:1) => matrix_px1_admm(1:nsev)
395 6 : CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
396 : CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set_admm%rho_atom_set, &
397 : admm_env%admm_gapw_env%admm_kind_set, &
398 6 : admm_env%admm_gapw_env%oce, sab_aux_fit, para_env)
399 : CALL prepare_gapw_den(qs_env, local_rho_set_admm, do_rho0=.FALSE., &
400 6 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
401 : !
402 6 : CALL local_rho_set_create(local_rho_set_f_admm)
403 : CALL allocate_rho_atom_internals(local_rho_set_f_admm%rho_atom_set, atomic_kind_set, &
404 6 : admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
405 : CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set_f_admm%rho_atom_set, &
406 : admm_env%admm_gapw_env%admm_kind_set, &
407 6 : admm_env%admm_gapw_env%oce, sab_aux_fit, para_env)
408 : CALL prepare_gapw_den(qs_env, local_rho_set_f_admm, do_rho0=.FALSE., &
409 6 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
410 : !
411 6 : CALL local_rho_set_create(local_rho_set_g_admm)
412 : CALL allocate_rho_atom_internals(local_rho_set_g_admm%rho_atom_set, atomic_kind_set, &
413 6 : admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
414 : CALL calculate_rho_atom_coeff(qs_env, mpga, local_rho_set_g_admm%rho_atom_set, &
415 : admm_env%admm_gapw_env%admm_kind_set, &
416 6 : admm_env%admm_gapw_env%oce, sab_aux_fit, para_env)
417 : CALL prepare_gapw_den(qs_env, local_rho_set_g_admm, do_rho0=.FALSE., &
418 6 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
419 : END IF
420 : END IF
421 : END IF
422 : END IF
423 : !
424 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
425 420 : poisson_env=poisson_env)
426 :
427 420 : NULLIFY (rhox_tau_g, rhox_tau_r, rhoxx_tau_g, rhoxx_tau_r)
428 3080 : ALLOCATE (rhox_r(nsev), rhox_g(nsev))
429 910 : DO ispin = 1, SIZE(evect, 1)
430 490 : CALL auxbas_pw_pool%create_pw(rhox_r(ispin))
431 910 : CALL auxbas_pw_pool%create_pw(rhox_g(ispin))
432 : END DO
433 420 : IF (needs_tau_response) THEN
434 0 : ALLOCATE (rhox_tau_r(nsev), rhox_tau_g(nsev))
435 0 : DO ispin = 1, SIZE(evect, 1)
436 0 : CALL auxbas_pw_pool%create_pw(rhox_tau_r(ispin))
437 0 : CALL auxbas_pw_pool%create_pw(rhox_tau_g(ispin))
438 : END DO
439 : END IF
440 420 : CALL auxbas_pw_pool%create_pw(rhox_tot_gspace)
441 :
442 420 : CALL pw_zero(rhox_tot_gspace)
443 910 : DO ispin = 1, nsev
444 : ! Calculate gridpoint values of the density associated to 2*matrix_px1 = C*X^T + X*C^T
445 490 : IF (nspins == 2) CALL dbcsr_scale(matrix_px1(ispin)%matrix, 2.0_dp)
446 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_px1(ispin)%matrix, &
447 : rho=rhox_r(ispin), rho_gspace=rhox_g(ispin), &
448 490 : soft_valid=gapw)
449 490 : IF (needs_tau_response) THEN
450 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_px1(ispin)%matrix, &
451 : rho=rhox_tau_r(ispin), rho_gspace=rhox_tau_g(ispin), &
452 0 : soft_valid=gapw, compute_tau=.TRUE.)
453 : END IF
454 : ! rhox_tot_gspace contains the values on the grid points of rhox = sum_munu 4D^X_munu*mu(r)*nu(r)
455 490 : CALL pw_axpy(rhox_g(ispin), rhox_tot_gspace)
456 : ! Recover matrix_px1 = (C*X^T + X*C^T)/2
457 910 : IF (nspins == 2) CALL dbcsr_scale(matrix_px1(ispin)%matrix, 0.5_dp)
458 : END DO
459 :
460 420 : IF (gapw_xc) THEN
461 100 : ALLOCATE (rhoxx_r(nsev), rhoxx_g(nsev))
462 40 : DO ispin = 1, nsev
463 20 : CALL auxbas_pw_pool%create_pw(rhoxx_r(ispin))
464 40 : CALL auxbas_pw_pool%create_pw(rhoxx_g(ispin))
465 : END DO
466 20 : IF (needs_tau_response) THEN
467 0 : ALLOCATE (rhoxx_tau_r(nsev), rhoxx_tau_g(nsev))
468 0 : DO ispin = 1, nsev
469 0 : CALL auxbas_pw_pool%create_pw(rhoxx_tau_r(ispin))
470 0 : CALL auxbas_pw_pool%create_pw(rhoxx_tau_g(ispin))
471 : END DO
472 : END IF
473 40 : DO ispin = 1, nsev
474 20 : IF (nspins == 2) CALL dbcsr_scale(matrix_px1(ispin)%matrix, 2.0_dp)
475 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_px1(ispin)%matrix, &
476 : rho=rhoxx_r(ispin), rho_gspace=rhoxx_g(ispin), &
477 20 : soft_valid=gapw_xc)
478 20 : IF (needs_tau_response) THEN
479 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_px1(ispin)%matrix, &
480 : rho=rhoxx_tau_r(ispin), rho_gspace=rhoxx_tau_g(ispin), &
481 0 : soft_valid=gapw_xc, compute_tau=.TRUE.)
482 : END IF
483 40 : IF (nspins == 2) CALL dbcsr_scale(matrix_px1(ispin)%matrix, 0.5_dp)
484 : END DO
485 : END IF
486 :
487 420 : CALL get_qs_env(qs_env, matrix_s=matrix_s, force=force)
488 :
489 420 : IF (.NOT. (is_rks_triplets .OR. do_sf)) THEN
490 362 : CALL auxbas_pw_pool%create_pw(xv_hartree_rspace)
491 362 : CALL auxbas_pw_pool%create_pw(xv_hartree_gspace)
492 : ! calculate associated hartree potential
493 362 : IF (gapw) THEN
494 88 : CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rhox_tot_gspace)
495 88 : IF (ASSOCIATED(local_rho_set%rho0_mpole%rhoz_cneo_s_gs)) THEN
496 0 : CALL pw_axpy(local_rho_set%rho0_mpole%rhoz_cneo_s_gs, rhox_tot_gspace)
497 : END IF
498 : END IF
499 : CALL pw_poisson_solve(poisson_env, rhox_tot_gspace, xehartree, &
500 362 : xv_hartree_gspace)
501 362 : CALL pw_transfer(xv_hartree_gspace, xv_hartree_rspace)
502 362 : CALL pw_scale(xv_hartree_rspace, xv_hartree_rspace%pw_grid%dvol)
503 : !
504 632 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
505 362 : NULLIFY (matrix_hx)
506 362 : CALL dbcsr_allocate_matrix_set(matrix_hx, nspins)
507 794 : DO ispin = 1, nspins
508 432 : ALLOCATE (matrix_hx(ispin)%matrix)
509 432 : CALL dbcsr_create(matrix_hx(ispin)%matrix, template=matrix_s(1)%matrix)
510 432 : CALL dbcsr_copy(matrix_hx(ispin)%matrix, matrix_s(1)%matrix)
511 432 : CALL dbcsr_set(matrix_hx(ispin)%matrix, 0.0_dp)
512 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=xv_hartree_rspace, &
513 : pmat=matrix_px1(ispin), hmat=matrix_hx(ispin), &
514 794 : gapw=gapw, calculate_forces=.TRUE.)
515 : END DO
516 362 : IF (debug_forces) THEN
517 360 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
518 90 : CALL para_env%sum(fodeb)
519 90 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dKh[Dx] ", fodeb
520 : END IF
521 362 : IF (gapw) THEN
522 322 : IF (debug_forces) fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1)
523 : CALL Vh_1c_gg_integrals(qs_env, eh1c, hartree_local%ecoul_1c, local_rho_set, para_env, tddft=.TRUE., &
524 88 : core_2nd=.TRUE.)
525 88 : IF (nspins == 1) THEN
526 88 : kval = 1.0_dp
527 : ELSE
528 0 : kval = 0.5_dp
529 : END IF
530 : CALL integrate_vhg0_rspace(qs_env, xv_hartree_rspace, para_env, calculate_forces=.TRUE., &
531 88 : local_rho_set=local_rho_set, kforce=kval)
532 88 : IF (debug_forces) THEN
533 312 : fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1) - fodeb(1:3)
534 78 : CALL para_env%sum(fodeb)
535 78 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dKh[Dx]PAWg0", fodeb
536 : END IF
537 322 : IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
538 : CALL update_ks_atom(qs_env, matrix_hx, matrix_px1, forces=.TRUE., &
539 88 : rho_atom_external=local_rho_set%rho_atom_set)
540 88 : IF (debug_forces) THEN
541 312 : fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
542 78 : CALL para_env%sum(fodeb)
543 78 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dKh[Dx]PAW", fodeb
544 : END IF
545 : END IF
546 : END IF
547 :
548 : ! XC
549 420 : IF (full_kernel%do_exck) THEN
550 0 : CPABORT("NYA")
551 : END IF
552 420 : NULLIFY (fxc_rho, fxc_tau, gxc_rho, gxc_tau)
553 420 : xc_section => full_kernel%xc_section
554 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", &
555 420 : i_val=myfun)
556 420 : IF (.NOT. ((myfun == xc_none) .OR. (tddfpt_control%spinflip == tddfpt_sf_col))) THEN
557 298 : SELECT CASE (ex_env%xc_kernel_method)
558 : CASE (xc_kernel_method_best)
559 : do_analytic = .TRUE.
560 0 : do_numeric = .TRUE.
561 : CASE (xc_kernel_method_analytic)
562 0 : do_analytic = .TRUE.
563 0 : do_numeric = .FALSE.
564 : CASE (xc_kernel_method_numeric)
565 0 : do_analytic = .FALSE.
566 0 : do_numeric = .TRUE.
567 : CASE DEFAULT
568 298 : CPABORT("invalid xc_kernel_method")
569 : END SELECT
570 298 : order = ex_env%diff_order
571 298 : eps_delta = ex_env%eps_delta_rho
572 :
573 298 : IF (gapw_xc) THEN
574 20 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, rho_xc=rho)
575 : ELSE
576 278 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, rho=rho)
577 : END IF
578 298 : CALL qs_rho_get(rho, rho_ao=matrix_p)
579 : NULLIFY (rhox)
580 298 : ALLOCATE (rhox)
581 : ! Create rhox object to collect all information on matrix_px1, including its values on the
582 : ! grid points
583 298 : CALL qs_rho_create(rhox)
584 298 : IF (gapw_xc) THEN
585 20 : IF (needs_tau_response) THEN
586 : CALL qs_rho_set(rho_struct=rhox, rho_ao=matrix_px1, rho_r=rhoxx_r, rho_g=rhoxx_g, &
587 : tau_r=rhoxx_tau_r, tau_g=rhoxx_tau_g, &
588 : rho_r_valid=.TRUE., rho_g_valid=.TRUE., &
589 0 : tau_r_valid=.TRUE., tau_g_valid=.TRUE.)
590 : ELSE
591 : CALL qs_rho_set(rho_struct=rhox, rho_ao=matrix_px1, rho_r=rhoxx_r, rho_g=rhoxx_g, &
592 20 : rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
593 : END IF
594 : ELSE
595 278 : IF (needs_tau_response) THEN
596 : CALL qs_rho_set(rho_struct=rhox, rho_ao=matrix_px1, rho_r=rhox_r, rho_g=rhox_g, &
597 : tau_r=rhox_tau_r, tau_g=rhox_tau_g, &
598 : rho_r_valid=.TRUE., rho_g_valid=.TRUE., &
599 0 : tau_r_valid=.TRUE., tau_g_valid=.TRUE.)
600 : ELSE
601 : CALL qs_rho_set(rho_struct=rhox, rho_ao=matrix_px1, rho_r=rhox_r, rho_g=rhox_g, &
602 278 : rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
603 : END IF
604 : END IF
605 : ! Calculate the exchange-correlation kernel derivative contribution, notice that for open-shell
606 : ! rhox_r contains a factor of 2!
607 298 : IF (do_analytic .AND. .NOT. do_numeric) THEN
608 0 : IF (.NOT. do_sf) THEN
609 0 : CPABORT("Analytic 3rd EXC derivatives not available")
610 : ELSE !TODO
611 0 : CALL get_qs_env(qs_env, xcint_weights=weights)
612 : CALL qs_fgxc_analytic(rho, rhox, xc_section, weights, auxbas_pw_pool, &
613 0 : fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip=do_sf)
614 : END IF
615 298 : ELSEIF (do_numeric) THEN
616 298 : IF (do_analytic) THEN
617 : CALL qs_fgxc_gdiff(ks_env, rho, rhox, xc_section, order, eps_delta, is_rks_triplets, &
618 298 : weights, fxc_rho, fxc_tau, gxc_rho, gxc_tau, spinflip=do_sf)
619 : ELSE
620 : CALL qs_fgxc_create(ks_env, rho, rhox, xc_section, order, is_rks_triplets, &
621 0 : fxc_rho, fxc_tau, gxc_rho, gxc_tau)
622 : END IF
623 : ELSE
624 0 : CPABORT("FHXC forces analytic/numeric")
625 : END IF
626 :
627 298 : IF (nspins == 2) THEN
628 132 : DO ispin = 1, nspins
629 88 : CALL pw_scale(gxc_rho(ispin), 0.5_dp)
630 132 : IF (ASSOCIATED(gxc_tau)) CALL pw_scale(gxc_tau(ispin), 0.5_dp)
631 : END DO
632 : END IF
633 298 : IF (gapw .OR. gapw_xc) THEN
634 108 : IF (do_analytic .AND. .NOT. do_numeric) THEN
635 0 : CPABORT("Analytic 3rd EXC derivatives not available")
636 108 : ELSEIF (do_numeric) THEN
637 108 : IF (do_analytic) THEN
638 : CALL gfxc_atom_diff(qs_env, ex_env%local_rho_set%rho_atom_set, &
639 : local_rho_set_f%rho_atom_set, local_rho_set_g%rho_atom_set, &
640 108 : qs_kind_set, xc_section, is_rks_triplets, order, eps_delta)
641 : ELSE
642 : CALL calculate_gfxc_atom(qs_env, ex_env%local_rho_set%rho_atom_set, &
643 : local_rho_set_f%rho_atom_set, local_rho_set_g%rho_atom_set, &
644 0 : qs_kind_set, xc_section, is_rks_triplets, order)
645 : END IF
646 : ELSE
647 0 : CPABORT("FHXC forces analytic/numeric")
648 : END IF
649 : END IF
650 :
651 568 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
652 298 : NULLIFY (matrix_fx)
653 298 : CALL dbcsr_allocate_matrix_set(matrix_fx, SIZE(fxc_rho))
654 632 : DO ispin = 1, SIZE(fxc_rho, 1)
655 334 : ALLOCATE (matrix_fx(ispin)%matrix)
656 334 : CALL dbcsr_create(matrix_fx(ispin)%matrix, template=matrix_s(1)%matrix)
657 334 : CALL dbcsr_copy(matrix_fx(ispin)%matrix, matrix_s(1)%matrix)
658 334 : CALL dbcsr_set(matrix_fx(ispin)%matrix, 0.0_dp)
659 334 : CALL pw_scale(fxc_rho(ispin), fxc_rho(ispin)%pw_grid%dvol)
660 : ! Calculate 2sum_sigmatau<munu|fxc|sigmatau>D^X_sigmatau
661 : ! fxc_rho here containes a factor of 2
662 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=fxc_rho(ispin), &
663 : pmat=matrix_px1(ispin), hmat=matrix_fx(ispin), &
664 858 : gapw=(gapw .OR. gapw_xc), calculate_forces=.TRUE.)
665 : END DO
666 298 : IF (debug_forces) THEN
667 360 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
668 90 : CALL para_env%sum(fodeb)
669 90 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dfxc[Dx] ", fodeb
670 : END IF
671 :
672 568 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
673 298 : NULLIFY (matrix_gx)
674 298 : CALL dbcsr_allocate_matrix_set(matrix_gx, nspins)
675 : ! Calculate exchange-correlation kernel derivative 2<D^X D^X|gxc|mu nu>
676 : ! gxc comes with a factor of 4, so a factor of 1/2 is introduced
677 640 : DO ispin = 1, nspins
678 342 : ALLOCATE (matrix_gx(ispin)%matrix)
679 342 : CALL dbcsr_create(matrix_gx(ispin)%matrix, template=matrix_s(1)%matrix)
680 342 : CALL dbcsr_copy(matrix_gx(ispin)%matrix, matrix_s(1)%matrix)
681 342 : CALL dbcsr_set(matrix_gx(ispin)%matrix, 0.0_dp)
682 342 : CALL pw_scale(gxc_rho(ispin), gxc_rho(ispin)%pw_grid%dvol)
683 342 : CALL pw_scale(gxc_rho(ispin), 0.5_dp)
684 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=gxc_rho(ispin), &
685 : pmat=matrix_p(ispin), hmat=matrix_gx(ispin), &
686 576 : gapw=(gapw .OR. gapw_xc), calculate_forces=.TRUE.)
687 640 : CALL dbcsr_scale(matrix_gx(ispin)%matrix, 2.0_dp)
688 : END DO
689 298 : IF (debug_forces) THEN
690 360 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
691 90 : CALL para_env%sum(fodeb)
692 90 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dgxc[Dx]", fodeb
693 : END IF
694 298 : CALL qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
695 :
696 : ! grid weight contribution to forces
697 568 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
698 298 : CALL accint_weight_force(qs_env, rho, rhox, 2, xc_section, is_rks_triplets)
699 298 : IF (debug_forces) THEN
700 360 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
701 90 : CALL para_env%sum(fodeb)
702 90 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*fxc[Dx]*dw", fodeb
703 : END IF
704 :
705 : ! Well, this is a hack :-(
706 : ! When qs_rho_set() was called on rhox it assumed ownership of the passed arrays.
707 : ! However, these arrays actually belong to ex_env. Hence, we can not call qs_rho_release()
708 : ! because this would release the arrays. Instead we're simply going to deallocate rhox.
709 298 : DEALLOCATE (rhox)
710 :
711 298 : IF (gapw .OR. gapw_xc) THEN
712 378 : IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
713 : CALL update_ks_atom(qs_env, matrix_fx, matrix_px1, forces=.TRUE., tddft=.TRUE., &
714 : rho_atom_external=local_rho_set_f%rho_atom_set, &
715 108 : kintegral=1.0_dp, kforce=1.0_dp)
716 108 : IF (debug_forces) THEN
717 360 : fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
718 90 : CALL para_env%sum(fodeb)
719 90 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dfxc[Dx]PAW ", fodeb
720 : END IF
721 378 : IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
722 108 : IF (nspins == 1) THEN
723 : CALL update_ks_atom(qs_env, matrix_gx, matrix_p, forces=.TRUE., tddft=.TRUE., &
724 : rho_atom_external=local_rho_set_g%rho_atom_set, &
725 108 : kscale=0.5_dp)
726 : ELSE
727 : CALL update_ks_atom(qs_env, matrix_gx, matrix_p, forces=.TRUE., &
728 : rho_atom_external=local_rho_set_g%rho_atom_set, &
729 0 : kintegral=0.5_dp, kforce=0.25_dp)
730 : END IF
731 108 : IF (debug_forces) THEN
732 360 : fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
733 90 : CALL para_env%sum(fodeb)
734 90 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dgxc[Dx]PAW ", fodeb
735 : END IF
736 : END IF
737 : END IF
738 :
739 : ! ADMM XC correction Exc[rho_admm]
740 420 : IF (do_admm .AND. tddfpt_control%admm_xc_correction .AND. (tddfpt_control%spinflip /= tddfpt_sf_col)) THEN
741 62 : IF (admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
742 : ! nothing to do
743 : ELSE
744 36 : IF (.NOT. tddfpt_control%admm_symm) THEN
745 0 : CALL cp_warn(__LOCATION__, "Forces need symmetric ADMM kernel corrections")
746 0 : CPABORT("ADMM KERNEL CORRECTION")
747 : END IF
748 36 : xc_section => admm_env%xc_section_aux
749 36 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
750 36 : needs = xc_functionals_get_needs(xc_fun_section, (nspins == 2), .TRUE.)
751 36 : needs_tau_response_aux = needs%tau .OR. needs%tau_spin
752 : CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit, matrix_s_aux_fit=matrix_s_aux_fit, &
753 36 : task_list_aux_fit=task_list)
754 36 : basis_type = "AUX_FIT"
755 36 : IF (admm_env%do_gapw) THEN
756 6 : basis_type = "AUX_FIT_SOFT"
757 6 : task_list => admm_env%admm_gapw_env%task_list
758 : END IF
759 : !
760 36 : NULLIFY (mfx, mgx)
761 36 : CALL dbcsr_allocate_matrix_set(mfx, nsev)
762 36 : CALL dbcsr_allocate_matrix_set(mgx, nspins)
763 72 : DO ispin = 1, nsev
764 36 : ALLOCATE (mfx(ispin)%matrix)
765 36 : CALL dbcsr_create(mfx(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
766 36 : CALL dbcsr_copy(mfx(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
767 72 : CALL dbcsr_set(mfx(ispin)%matrix, 0.0_dp)
768 : END DO
769 72 : DO ispin = 1, nspins
770 36 : ALLOCATE (mgx(ispin)%matrix)
771 36 : CALL dbcsr_create(mgx(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
772 36 : CALL dbcsr_copy(mgx(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
773 72 : CALL dbcsr_set(mgx(ispin)%matrix, 0.0_dp)
774 : END DO
775 :
776 : ! ADMM density and response density
777 36 : NULLIFY (rho_g_aux, rho_r_aux, rhox_g_aux, rhox_r_aux, rhox_tau_g_aux, rhox_tau_r_aux)
778 36 : CALL qs_rho_get(rho_aux_fit, rho_r=rho_r_aux, rho_g=rho_g_aux)
779 36 : CALL qs_rho_get(rho_aux_fit, rho_ao=matrix_p_admm)
780 : ! rhox_aux
781 252 : ALLOCATE (rhox_r_aux(nsev), rhox_g_aux(nsev))
782 72 : DO ispin = 1, nsev
783 36 : CALL auxbas_pw_pool%create_pw(rhox_r_aux(ispin))
784 72 : CALL auxbas_pw_pool%create_pw(rhox_g_aux(ispin))
785 : END DO
786 36 : IF (needs_tau_response_aux) THEN
787 0 : ALLOCATE (rhox_tau_r_aux(nsev), rhox_tau_g_aux(nsev))
788 0 : DO ispin = 1, nsev
789 0 : CALL auxbas_pw_pool%create_pw(rhox_tau_r_aux(ispin))
790 0 : CALL auxbas_pw_pool%create_pw(rhox_tau_g_aux(ispin))
791 : END DO
792 : END IF
793 72 : DO ispin = 1, nsev
794 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_px1_admm(ispin)%matrix, &
795 : rho=rhox_r_aux(ispin), rho_gspace=rhox_g_aux(ispin), &
796 : basis_type=basis_type, &
797 36 : task_list_external=task_list)
798 72 : IF (needs_tau_response_aux) THEN
799 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_px1_admm(ispin)%matrix, &
800 : rho=rhox_tau_r_aux(ispin), &
801 : rho_gspace=rhox_tau_g_aux(ispin), &
802 : basis_type=basis_type, task_list_external=task_list, &
803 0 : compute_tau=.TRUE.)
804 : END IF
805 : END DO
806 : !
807 : NULLIFY (rhox_aux)
808 36 : ALLOCATE (rhox_aux)
809 36 : CALL qs_rho_create(rhox_aux)
810 36 : IF (needs_tau_response_aux) THEN
811 : CALL qs_rho_set(rho_struct=rhox_aux, rho_ao=matrix_px1_admm, &
812 : rho_r=rhox_r_aux, rho_g=rhox_g_aux, &
813 : tau_r=rhox_tau_r_aux, tau_g=rhox_tau_g_aux, &
814 : rho_r_valid=.TRUE., rho_g_valid=.TRUE., &
815 0 : tau_r_valid=.TRUE., tau_g_valid=.TRUE.)
816 : ELSE
817 : CALL qs_rho_set(rho_struct=rhox_aux, rho_ao=matrix_px1_admm, &
818 : rho_r=rhox_r_aux, rho_g=rhox_g_aux, &
819 36 : rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
820 : END IF
821 36 : IF (do_analytic .AND. .NOT. do_numeric) THEN
822 0 : CPABORT("Analytic 3rd derivatives of EXC not available")
823 36 : ELSEIF (do_numeric) THEN
824 36 : IF (do_analytic) THEN
825 : CALL qs_fgxc_gdiff(ks_env, rho_aux_fit, rhox_aux, xc_section, order, eps_delta, &
826 36 : is_rks_triplets, weights, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
827 : ELSE
828 : CALL qs_fgxc_create(ks_env, rho_aux_fit, rhox_aux, xc_section, &
829 0 : order, is_rks_triplets, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
830 : END IF
831 : ELSE
832 0 : CPABORT("FHXC forces analytic/numeric")
833 : END IF
834 :
835 : ! grid weight contribution to forces ADMM XC correction term
836 54 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
837 : !
838 36 : CALL accint_weight_force(qs_env, rho_aux_fit, rhox_aux, 2, xc_section, is_rks_triplets)
839 : !
840 36 : IF (debug_forces) THEN
841 24 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
842 6 : CALL para_env%sum(fodeb)
843 6 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx_admm*fxc[Dx_admm]*dw", fodeb
844 : END IF
845 :
846 : ! Well, this is a hack :-(
847 : ! When qs_rho_set() was called on rhox_aux it assumed ownership of the passed arrays.
848 : ! However, these arrays actually belong to ex_env. Hence, we can not call qs_rho_release()
849 : ! because this would release the arrays. Instead we're simply going to deallocate rhox_aux.
850 36 : DEALLOCATE (rhox_aux)
851 :
852 72 : DO ispin = 1, nsev
853 36 : CALL auxbas_pw_pool%give_back_pw(rhox_r_aux(ispin))
854 72 : CALL auxbas_pw_pool%give_back_pw(rhox_g_aux(ispin))
855 : END DO
856 36 : DEALLOCATE (rhox_r_aux, rhox_g_aux)
857 36 : IF (needs_tau_response_aux) THEN
858 0 : DO ispin = 1, nsev
859 0 : CALL auxbas_pw_pool%give_back_pw(rhox_tau_r_aux(ispin))
860 0 : CALL auxbas_pw_pool%give_back_pw(rhox_tau_g_aux(ispin))
861 : END DO
862 0 : DEALLOCATE (rhox_tau_r_aux, rhox_tau_g_aux)
863 : END IF
864 36 : fscal = 1.0_dp
865 36 : IF (nspins == 2) fscal = 2.0_dp
866 : !
867 54 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
868 72 : DO ispin = 1, nsev
869 36 : CALL pw_scale(fxc_rho(ispin), fxc_rho(ispin)%pw_grid%dvol)
870 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=fxc_rho(ispin), &
871 : hmat=mfx(ispin), &
872 : pmat=matrix_px1_admm(ispin), &
873 : basis_type=basis_type, &
874 : calculate_forces=.TRUE., &
875 : force_adm=fscal, &
876 72 : task_list_external=task_list)
877 : END DO
878 36 : IF (debug_forces) THEN
879 24 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
880 6 : CALL para_env%sum(fodeb)
881 6 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dfxc[Dx]ADMM", fodeb
882 : END IF
883 :
884 54 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
885 72 : DO ispin = 1, nsev
886 36 : CALL pw_scale(gxc_rho(ispin), gxc_rho(ispin)%pw_grid%dvol)
887 36 : CALL pw_scale(gxc_rho(ispin), 0.5_dp)
888 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=gxc_rho(ispin), &
889 : hmat=mgx(ispin), &
890 : pmat=matrix_p_admm(ispin), &
891 : basis_type=basis_type, &
892 : calculate_forces=.TRUE., &
893 : force_adm=fscal, &
894 36 : task_list_external=task_list)
895 72 : CALL dbcsr_scale(mgx(ispin)%matrix, 2.0_dp)
896 : END DO
897 36 : IF (debug_forces) THEN
898 24 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
899 6 : CALL para_env%sum(fodeb)
900 6 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dgxc[Dx]ADMM", fodeb
901 : END IF
902 36 : CALL qs_fgxc_release(ks_env, fxc_rho, fxc_tau, gxc_rho, gxc_tau)
903 : !
904 36 : IF (admm_env%do_gapw) THEN
905 6 : CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
906 6 : rho_atom_set => admm_env%admm_gapw_env%local_rho_set%rho_atom_set
907 6 : rho_atom_set_f => local_rho_set_f_admm%rho_atom_set
908 6 : rho_atom_set_g => local_rho_set_g_admm%rho_atom_set
909 :
910 6 : IF (do_analytic .AND. .NOT. do_numeric) THEN
911 0 : CPABORT("Analytic 3rd EXC derivatives not available")
912 6 : ELSEIF (do_numeric) THEN
913 6 : IF (do_analytic) THEN
914 : CALL gfxc_atom_diff(qs_env, rho_atom_set, &
915 : rho_atom_set_f, rho_atom_set_g, &
916 : admm_env%admm_gapw_env%admm_kind_set, xc_section, &
917 6 : is_rks_triplets, order, eps_delta)
918 : ELSE
919 : CALL calculate_gfxc_atom(qs_env, rho_atom_set, &
920 : rho_atom_set_f, rho_atom_set_g, &
921 : admm_env%admm_gapw_env%admm_kind_set, xc_section, &
922 0 : is_rks_triplets, order)
923 : END IF
924 : ELSE
925 0 : CPABORT("FHXC forces analytic/numeric")
926 : END IF
927 :
928 24 : IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
929 6 : IF (nspins == 1) THEN
930 : CALL update_ks_atom(qs_env, mfx, matrix_px1_admm, forces=.TRUE., &
931 : rho_atom_external=rho_atom_set_f, &
932 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
933 : oce_external=admm_env%admm_gapw_env%oce, sab_external=sab_aux_fit, &
934 6 : kintegral=2.0_dp, kforce=0.5_dp)
935 : ELSE
936 : CALL update_ks_atom(qs_env, mfx, matrix_px1_admm, forces=.TRUE., &
937 : rho_atom_external=rho_atom_set_f, &
938 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
939 : oce_external=admm_env%admm_gapw_env%oce, sab_external=sab_aux_fit, &
940 0 : kintegral=2.0_dp, kforce=1.0_dp)
941 : END IF
942 6 : IF (debug_forces) THEN
943 24 : fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
944 6 : CALL para_env%sum(fodeb)
945 6 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dfxc[Dx]ADMM-PAW ", fodeb
946 : END IF
947 24 : IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
948 6 : IF (nspins == 1) THEN
949 : CALL update_ks_atom(qs_env, mgx, matrix_p, forces=.TRUE., &
950 : rho_atom_external=rho_atom_set_g, &
951 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
952 : oce_external=admm_env%admm_gapw_env%oce, sab_external=sab_aux_fit, &
953 6 : kintegral=1.0_dp, kforce=0.5_dp)
954 : ELSE
955 : CALL update_ks_atom(qs_env, mgx, matrix_p, forces=.TRUE., &
956 : rho_atom_external=rho_atom_set_g, &
957 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
958 : oce_external=admm_env%admm_gapw_env%oce, sab_external=sab_aux_fit, &
959 0 : kintegral=1.0_dp, kforce=1.0_dp)
960 : END IF
961 6 : IF (debug_forces) THEN
962 24 : fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
963 6 : CALL para_env%sum(fodeb)
964 6 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dgxc[Dx]ADMM-PAW ", fodeb
965 : END IF
966 : END IF
967 : !
968 : ! A' fx A - Forces
969 : !
970 54 : IF (debug_forces) fodeb(1:3) = force(1)%overlap_admm(1:3, 1)
971 36 : fval = 2.0_dp*REAL(nspins, KIND=dp)
972 36 : CALL admm_projection_derivative(qs_env, mfx, matrix_px1, fval)
973 36 : IF (debug_forces) THEN
974 24 : fodeb(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb(1:3)
975 6 : CALL para_env%sum(fodeb)
976 6 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dfXC(P)*S' ", fodeb
977 : END IF
978 54 : IF (debug_forces) fodeb(1:3) = force(1)%overlap_admm(1:3, 1)
979 : fval = 2.0_dp*REAL(nspins, KIND=dp)
980 36 : CALL admm_projection_derivative(qs_env, mgx, matrix_p, fval)
981 36 : IF (debug_forces) THEN
982 24 : fodeb(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb(1:3)
983 6 : CALL para_env%sum(fodeb)
984 6 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dgXC(P)*S' ", fodeb
985 : END IF
986 : !
987 : ! Add ADMM fx/gx to the full basis fx/gx
988 36 : fscal = 1.0_dp
989 36 : IF (nspins == 2) fscal = 2.0_dp
990 36 : nao = admm_env%nao_orb
991 36 : nao_aux = admm_env%nao_aux_fit
992 36 : ALLOCATE (dbwork)
993 36 : CALL dbcsr_create(dbwork, template=matrix_fx(1)%matrix)
994 72 : DO ispin = 1, nsev
995 : ! fx
996 : CALL cp_dbcsr_sm_fm_multiply(mfx(ispin)%matrix, admm_env%A, &
997 36 : admm_env%work_aux_orb, nao)
998 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
999 : 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
1000 36 : admm_env%work_orb_orb)
1001 36 : CALL dbcsr_copy(dbwork, matrix_fx(1)%matrix)
1002 36 : CALL dbcsr_set(dbwork, 0.0_dp)
1003 36 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.TRUE.)
1004 36 : CALL dbcsr_add(matrix_fx(ispin)%matrix, dbwork, 1.0_dp, fscal)
1005 : ! gx
1006 : CALL cp_dbcsr_sm_fm_multiply(mgx(ispin)%matrix, admm_env%A, &
1007 36 : admm_env%work_aux_orb, nao)
1008 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
1009 : 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
1010 36 : admm_env%work_orb_orb)
1011 36 : CALL dbcsr_set(dbwork, 0.0_dp)
1012 36 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.TRUE.)
1013 72 : CALL dbcsr_add(matrix_gx(ispin)%matrix, dbwork, 1.0_dp, fscal)
1014 : END DO
1015 36 : CALL dbcsr_release(dbwork)
1016 36 : DEALLOCATE (dbwork)
1017 36 : CALL dbcsr_deallocate_matrix_set(mfx)
1018 72 : CALL dbcsr_deallocate_matrix_set(mgx)
1019 :
1020 : END IF
1021 : END IF
1022 :
1023 910 : DO ispin = 1, nsev
1024 490 : CALL auxbas_pw_pool%give_back_pw(rhox_r(ispin))
1025 910 : CALL auxbas_pw_pool%give_back_pw(rhox_g(ispin))
1026 : END DO
1027 420 : DEALLOCATE (rhox_r, rhox_g)
1028 420 : IF (needs_tau_response) THEN
1029 0 : DO ispin = 1, nsev
1030 0 : CALL auxbas_pw_pool%give_back_pw(rhox_tau_r(ispin))
1031 0 : CALL auxbas_pw_pool%give_back_pw(rhox_tau_g(ispin))
1032 : END DO
1033 0 : DEALLOCATE (rhox_tau_r, rhox_tau_g)
1034 : END IF
1035 420 : CALL auxbas_pw_pool%give_back_pw(rhox_tot_gspace)
1036 420 : IF (gapw_xc) THEN
1037 40 : DO ispin = 1, nsev
1038 20 : CALL auxbas_pw_pool%give_back_pw(rhoxx_r(ispin))
1039 40 : CALL auxbas_pw_pool%give_back_pw(rhoxx_g(ispin))
1040 : END DO
1041 20 : DEALLOCATE (rhoxx_r, rhoxx_g)
1042 20 : IF (needs_tau_response) THEN
1043 0 : DO ispin = 1, nsev
1044 0 : CALL auxbas_pw_pool%give_back_pw(rhoxx_tau_r(ispin))
1045 0 : CALL auxbas_pw_pool%give_back_pw(rhoxx_tau_g(ispin))
1046 : END DO
1047 0 : DEALLOCATE (rhoxx_tau_r, rhoxx_tau_g)
1048 : END IF
1049 : END IF
1050 420 : IF (.NOT. (is_rks_triplets .OR. do_sf)) THEN
1051 362 : CALL auxbas_pw_pool%give_back_pw(xv_hartree_rspace)
1052 362 : CALL auxbas_pw_pool%give_back_pw(xv_hartree_gspace)
1053 : END IF
1054 :
1055 : ! HFX
1056 420 : IF (do_hfx) THEN
1057 150 : NULLIFY (matrix_hfx, matrix_hfx_asymm)
1058 150 : CALL dbcsr_allocate_matrix_set(matrix_hfx, nsev)
1059 150 : CALL dbcsr_allocate_matrix_set(matrix_hfx_asymm, nsev)
1060 312 : DO ispin = 1, nsev
1061 162 : ALLOCATE (matrix_hfx(ispin)%matrix)
1062 162 : CALL dbcsr_create(matrix_hfx(ispin)%matrix, template=matrix_s(1)%matrix)
1063 162 : CALL dbcsr_copy(matrix_hfx(ispin)%matrix, matrix_s(1)%matrix)
1064 162 : CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
1065 :
1066 162 : ALLOCATE (matrix_hfx_asymm(ispin)%matrix)
1067 : CALL dbcsr_create(matrix_hfx_asymm(ispin)%matrix, template=matrix_s(1)%matrix, &
1068 162 : matrix_type=dbcsr_type_antisymmetric)
1069 312 : CALL dbcsr_complete_redistribute(matrix_hfx(ispin)%matrix, matrix_hfx_asymm(ispin)%matrix)
1070 : END DO
1071 : !
1072 150 : xc_section => full_kernel%xc_section
1073 150 : hfx_section => section_vals_get_subs_vals(xc_section, "HF")
1074 150 : CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
1075 150 : CPASSERT(n_rep_hf == 1)
1076 : CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
1077 150 : i_rep_section=1)
1078 150 : mspin = 1
1079 150 : IF (hfx_treat_lsd_in_core) mspin = nsev
1080 : !
1081 150 : CALL get_qs_env(qs_env=qs_env, x_data=x_data, s_mstruct_changed=s_mstruct_changed)
1082 150 : distribute_fock_matrix = .TRUE.
1083 : !
1084 150 : IF (do_admm) THEN
1085 78 : CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
1086 78 : NULLIFY (matrix_hfx_admm, matrix_hfx_admm_asymm)
1087 78 : CALL dbcsr_allocate_matrix_set(matrix_hfx_admm, nsev)
1088 78 : CALL dbcsr_allocate_matrix_set(matrix_hfx_admm_asymm, nsev)
1089 160 : DO ispin = 1, nsev
1090 82 : ALLOCATE (matrix_hfx_admm(ispin)%matrix)
1091 82 : CALL dbcsr_create(matrix_hfx_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
1092 82 : CALL dbcsr_copy(matrix_hfx_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
1093 82 : CALL dbcsr_set(matrix_hfx_admm(ispin)%matrix, 0.0_dp)
1094 :
1095 82 : ALLOCATE (matrix_hfx_admm_asymm(ispin)%matrix)
1096 : CALL dbcsr_create(matrix_hfx_admm_asymm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix, &
1097 82 : matrix_type=dbcsr_type_antisymmetric)
1098 160 : CALL dbcsr_complete_redistribute(matrix_hfx_admm(ispin)%matrix, matrix_hfx_admm_asymm(ispin)%matrix)
1099 : END DO
1100 : !
1101 78 : NULLIFY (mpe, mhe)
1102 632 : ALLOCATE (mpe(nsev, 1), mhe(nsev, 1))
1103 160 : DO ispin = 1, nsev
1104 82 : mhe(ispin, 1)%matrix => matrix_hfx_admm(ispin)%matrix
1105 160 : mpe(ispin, 1)%matrix => matrix_px1_admm(ispin)%matrix
1106 : END DO
1107 78 : IF (x_data(1, 1)%do_hfx_ri) THEN
1108 : eh1 = 0.0_dp
1109 : CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhe, eh1, rho_ao=mpe, &
1110 : geometry_did_change=s_mstruct_changed, nspins=nspins, &
1111 6 : hf_fraction=x_data(1, 1)%general_parameter%fraction)
1112 : ELSE
1113 144 : DO ispin = 1, mspin
1114 : eh1 = 0.0
1115 : CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
1116 : para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
1117 144 : ispin=ispin, nspins=nsev)
1118 : END DO
1119 : END IF
1120 : !anti-symmetric density
1121 160 : DO ispin = 1, nsev
1122 82 : mhe(ispin, 1)%matrix => matrix_hfx_admm_asymm(ispin)%matrix
1123 160 : mpe(ispin, 1)%matrix => matrix_px1_admm_asymm(ispin)%matrix
1124 : END DO
1125 78 : IF (x_data(1, 1)%do_hfx_ri) THEN
1126 : eh1 = 0.0_dp
1127 : CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhe, eh1, rho_ao=mpe, &
1128 : geometry_did_change=s_mstruct_changed, nspins=nspins, &
1129 6 : hf_fraction=x_data(1, 1)%general_parameter%fraction)
1130 : ELSE
1131 144 : DO ispin = 1, mspin
1132 : eh1 = 0.0
1133 : CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
1134 : para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
1135 144 : ispin=ispin, nspins=nsev)
1136 : END DO
1137 : END IF
1138 : !
1139 78 : nao = admm_env%nao_orb
1140 78 : nao_aux = admm_env%nao_aux_fit
1141 78 : ALLOCATE (dbwork, dbwork_asymm)
1142 78 : CALL dbcsr_create(dbwork, template=matrix_hfx(1)%matrix)
1143 78 : CALL dbcsr_create(dbwork_asymm, template=matrix_hfx(1)%matrix, matrix_type=dbcsr_type_antisymmetric)
1144 160 : DO ispin = 1, nsev
1145 : CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_admm(ispin)%matrix, admm_env%A, &
1146 82 : admm_env%work_aux_orb, nao)
1147 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
1148 : 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
1149 82 : admm_env%work_orb_orb)
1150 82 : CALL dbcsr_copy(dbwork, matrix_hfx(1)%matrix)
1151 82 : CALL dbcsr_set(dbwork, 0.0_dp)
1152 82 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.TRUE.)
1153 82 : CALL dbcsr_add(matrix_hfx(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
1154 : !anti-symmetric case
1155 : CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_admm_asymm(ispin)%matrix, admm_env%A, &
1156 82 : admm_env%work_aux_orb, nao)
1157 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
1158 : 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
1159 82 : admm_env%work_orb_orb)
1160 82 : CALL dbcsr_copy(dbwork_asymm, matrix_hfx_asymm(1)%matrix)
1161 82 : CALL dbcsr_set(dbwork_asymm, 0.0_dp)
1162 82 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork_asymm, keep_sparsity=.TRUE.)
1163 160 : CALL dbcsr_add(matrix_hfx_asymm(ispin)%matrix, dbwork_asymm, 1.0_dp, 1.0_dp)
1164 : END DO
1165 78 : CALL dbcsr_release(dbwork)
1166 78 : CALL dbcsr_release(dbwork_asymm)
1167 78 : DEALLOCATE (dbwork, dbwork_asymm)
1168 : ! forces
1169 : ! ADMM Projection force
1170 120 : IF (debug_forces) fodeb(1:3) = force(1)%overlap_admm(1:3, 1)
1171 78 : fval = 4.0_dp*REAL(nspins, KIND=dp)*0.5_dp !0.5 for symm/anti-symm
1172 78 : CALL admm_projection_derivative(qs_env, matrix_hfx_admm, matrix_px1, fval)
1173 78 : CALL admm_projection_derivative(qs_env, matrix_hfx_admm_asymm, matrix_px1_asymm, fval)
1174 78 : IF (debug_forces) THEN
1175 56 : fodeb(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb(1:3)
1176 14 : CALL para_env%sum(fodeb)
1177 14 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*Hx(P)*S' ", fodeb
1178 : END IF
1179 : !
1180 78 : use_virial = .FALSE.
1181 78 : NULLIFY (mdum)
1182 78 : fval = 2.0_dp*REAL(nspins, KIND=dp)*0.5_dp !0.5 factor because of symemtry/anti-symmetry
1183 : ! For SF TDDFT integrate_four_center and derivatives_four_center routines introduce a factor of 1/2
1184 78 : IF (do_sf) fval = fval*2.0_dp
1185 120 : IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
1186 160 : DO ispin = 1, nsev
1187 160 : mpe(ispin, 1)%matrix => matrix_px1_admm(ispin)%matrix
1188 : END DO
1189 78 : IF (x_data(1, 1)%do_hfx_ri) THEN
1190 : CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
1191 : x_data(1, 1)%general_parameter%fraction, &
1192 : rho_ao=mpe, rho_ao_resp=mdum, &
1193 6 : use_virial=use_virial, rescale_factor=fval)
1194 : ELSE
1195 : CALL derivatives_four_center(qs_env, mpe, mdum, hfx_section, para_env, 1, use_virial, &
1196 72 : adiabatic_rescale_factor=fval, nspins=nsev)
1197 : END IF
1198 160 : DO ispin = 1, nsev
1199 160 : mpe(ispin, 1)%matrix => matrix_px1_admm_asymm(ispin)%matrix
1200 : END DO
1201 78 : IF (x_data(1, 1)%do_hfx_ri) THEN
1202 : CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
1203 : x_data(1, 1)%general_parameter%fraction, &
1204 : rho_ao=mpe, rho_ao_resp=mdum, &
1205 6 : use_virial=use_virial, rescale_factor=fval)
1206 : ELSE
1207 : CALL derivatives_four_center(qs_env, mpe, mdum, hfx_section, para_env, 1, use_virial, &
1208 72 : adiabatic_rescale_factor=fval, nspins=SIZE(mpe, 1))
1209 : END IF
1210 78 : IF (debug_forces) THEN
1211 56 : fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
1212 14 : CALL para_env%sum(fodeb)
1213 14 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dhfx'*Dx ", fodeb
1214 : END IF
1215 : !
1216 78 : DEALLOCATE (mpe, mhe)
1217 : !
1218 78 : CALL dbcsr_deallocate_matrix_set(matrix_hfx_admm)
1219 78 : CALL dbcsr_deallocate_matrix_set(matrix_hfx_admm_asymm)
1220 : ELSE
1221 72 : NULLIFY (mpe, mhe)
1222 592 : ALLOCATE (mpe(nsev, 1), mhe(nsev, 1))
1223 152 : DO ispin = 1, nsev
1224 80 : mhe(ispin, 1)%matrix => matrix_hfx(ispin)%matrix
1225 152 : mpe(ispin, 1)%matrix => matrix_px1(ispin)%matrix
1226 : END DO
1227 72 : IF (x_data(1, 1)%do_hfx_ri) THEN
1228 : eh1 = 0.0_dp
1229 : CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhe, eh1, rho_ao=mpe, &
1230 : geometry_did_change=s_mstruct_changed, nspins=nspins, &
1231 18 : hf_fraction=x_data(1, 1)%general_parameter%fraction)
1232 : ELSE
1233 108 : DO ispin = 1, mspin
1234 : eh1 = 0.0
1235 : CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
1236 : para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
1237 108 : ispin=ispin, nspins=SIZE(mpe, 1))
1238 : END DO
1239 : END IF
1240 :
1241 : !anti-symmetric density matrix
1242 152 : DO ispin = 1, nsev
1243 80 : mhe(ispin, 1)%matrix => matrix_hfx_asymm(ispin)%matrix
1244 152 : mpe(ispin, 1)%matrix => matrix_px1_asymm(ispin)%matrix
1245 : END DO
1246 72 : IF (x_data(1, 1)%do_hfx_ri) THEN
1247 : eh1 = 0.0_dp
1248 : CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhe, eh1, rho_ao=mpe, &
1249 : geometry_did_change=s_mstruct_changed, nspins=nspins, &
1250 18 : hf_fraction=x_data(1, 1)%general_parameter%fraction)
1251 : ELSE
1252 108 : DO ispin = 1, mspin
1253 : eh1 = 0.0
1254 : CALL integrate_four_center(qs_env, x_data, mhe, eh1, mpe, hfx_section, &
1255 : para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
1256 108 : ispin=ispin, nspins=SIZE(mpe, 1))
1257 : END DO
1258 : END IF
1259 : ! forces
1260 72 : use_virial = .FALSE.
1261 72 : NULLIFY (mdum)
1262 72 : fval = 2.0_dp*REAL(nspins, KIND=dp)*0.5_dp !extra 0.5 factor because of symmetry/antisymemtry
1263 : ! For SF TDDFT integrate_four_center and derivatives_four_center routines introduce a factor of 1/2
1264 72 : IF (do_sf) fval = fval*2.0_dp
1265 120 : IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
1266 152 : DO ispin = 1, nsev
1267 152 : mpe(ispin, 1)%matrix => matrix_px1(ispin)%matrix
1268 : END DO
1269 72 : IF (x_data(1, 1)%do_hfx_ri) THEN
1270 : CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
1271 : x_data(1, 1)%general_parameter%fraction, &
1272 : rho_ao=mpe, rho_ao_resp=mdum, &
1273 18 : use_virial=use_virial, rescale_factor=fval)
1274 : ELSE
1275 : CALL derivatives_four_center(qs_env, mpe, mdum, hfx_section, para_env, 1, use_virial, &
1276 54 : adiabatic_rescale_factor=fval, nspins=SIZE(mpe, 1))
1277 : END IF
1278 152 : DO ispin = 1, nsev
1279 152 : mpe(ispin, 1)%matrix => matrix_px1_asymm(ispin)%matrix
1280 : END DO
1281 72 : IF (x_data(1, 1)%do_hfx_ri) THEN
1282 : CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
1283 : x_data(1, 1)%general_parameter%fraction, &
1284 : rho_ao=mpe, rho_ao_resp=mdum, &
1285 18 : use_virial=use_virial, rescale_factor=fval)
1286 : ELSE
1287 : CALL derivatives_four_center(qs_env, mpe, mdum, hfx_section, para_env, 1, use_virial, &
1288 54 : adiabatic_rescale_factor=fval, nspins=SIZE(mpe, 1))
1289 : END IF
1290 72 : IF (debug_forces) THEN
1291 64 : fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
1292 16 : CALL para_env%sum(fodeb)
1293 16 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Dx*dhfx*Dx ", fodeb
1294 : END IF
1295 : !
1296 72 : DEALLOCATE (mpe, mhe)
1297 : END IF
1298 150 : fval = 2.0_dp*REAL(nspins, KIND=dp)*0.5_dp !extra 0.5 because of symm/antisymm
1299 : ! For SF TDDFT integrate_four_center and derivatives_four_center routines introduce a factor of 1/2
1300 150 : IF (do_sf) fval = fval*2.0_dp
1301 312 : DO ispin = 1, nsev
1302 162 : CALL dbcsr_scale(matrix_hfx(ispin)%matrix, fval)
1303 312 : CALL dbcsr_scale(matrix_hfx_asymm(ispin)%matrix, fval)
1304 : END DO
1305 : END IF
1306 :
1307 420 : IF (gapw .OR. gapw_xc) THEN
1308 124 : CALL local_rho_set_release(local_rho_set)
1309 124 : CALL local_rho_set_release(local_rho_set_f)
1310 124 : CALL local_rho_set_release(local_rho_set_g)
1311 124 : IF (gapw) THEN
1312 104 : CALL hartree_local_release(hartree_local)
1313 : END IF
1314 : END IF
1315 420 : IF (do_admm) THEN
1316 78 : IF (admm_env%do_gapw) THEN
1317 24 : IF (tddfpt_control%admm_xc_correction) THEN
1318 18 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
1319 6 : CALL local_rho_set_release(local_rho_set_admm)
1320 6 : CALL local_rho_set_release(local_rho_set_f_admm)
1321 6 : CALL local_rho_set_release(local_rho_set_g_admm)
1322 : END IF
1323 : END IF
1324 : END IF
1325 : END IF
1326 :
1327 : ! HFX short range
1328 420 : IF (do_hfxsr) THEN
1329 0 : CPABORT("HFXSR not implemented")
1330 : END IF
1331 : ! HFX long range
1332 420 : IF (do_hfxlr) THEN
1333 0 : CPABORT("HFXLR not implemented")
1334 : END IF
1335 :
1336 420 : CALL get_qs_env(qs_env, sab_orb=sab_orb)
1337 420 : NULLIFY (matrix_wx1)
1338 420 : CALL dbcsr_allocate_matrix_set(matrix_wx1, nspins)
1339 420 : cpmos => ex_env%cpmos
1340 420 : focc = 2.0_dp
1341 420 : IF (nspins == 2) focc = 1.0_dp
1342 :
1343 : ! Initialize mos and dimensions of occupied space
1344 : ! In the following comments mos is referred to as Ca and mos2 as Cb
1345 420 : spin = 1
1346 420 : mos => gs_mos(1)%mos_occ
1347 420 : mosa => gs_mos(1)%mos_active
1348 420 : norb(1) = gs_mos(1)%nmo_occ
1349 420 : nactive(1) = gs_mos(1)%nmo_active
1350 420 : IF (nspins == 2) THEN
1351 82 : mos2 => gs_mos(2)%mos_occ
1352 82 : mosa2 => gs_mos(2)%mos_active
1353 82 : norb(2) = gs_mos(2)%nmo_occ
1354 82 : nactive(2) = gs_mos(2)%nmo_active
1355 : END IF
1356 : ! Build response vector, Eq. 49, and the third term of \Lamda_munu, Eq. 51
1357 922 : DO ispin = 1, nspins
1358 :
1359 502 : IF (nactive(ispin) == norb(ispin)) THEN
1360 502 : do_res = .FALSE.
1361 2570 : DO ia = 1, nactive(ispin)
1362 2570 : CPASSERT(ia == gs_mos(ispin)%index_active(ia))
1363 : END DO
1364 : ELSE
1365 : do_res = .TRUE.
1366 : END IF
1367 :
1368 : ! Initialize mos and dimensions of occupied space
1369 502 : IF (.NOT. do_sf) THEN
1370 478 : spin = ispin
1371 478 : mos => gs_mos(ispin)%mos_occ
1372 478 : mos2 => gs_mos(ispin)%mos_occ
1373 478 : mosa => gs_mos(ispin)%mos_active
1374 478 : mosa2 => gs_mos(ispin)%mos_active
1375 : END IF
1376 :
1377 : ! Create working fields for the response vector
1378 502 : CALL cp_fm_create(cpscr, gs_mos(ispin)%mos_active%matrix_struct, "cpscr", set_zero=.TRUE.)
1379 : !
1380 502 : CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, matrix_struct=fm_struct, nrow_global=nao)
1381 : !
1382 502 : CALL cp_fm_create(avamat, fm_struct, nrow=nactive(spin), ncol=nactive(spin))
1383 502 : CALL cp_fm_create(avcmat, fm_struct, nrow=nactive(spin), ncol=norb(spin))
1384 502 : CALL cp_fm_create(cvcmat, fm_struct, nrow=norb(spin), ncol=norb(spin))
1385 : !
1386 502 : CALL cp_fm_create(vcvec, gs_mos(ispin)%mos_occ%matrix_struct, "vcvec")
1387 502 : CALL cp_fm_create(vavec, gs_mos(ispin)%mos_active%matrix_struct, "vavec")
1388 :
1389 : ! Allocate and initialize the Lambda matrix
1390 502 : ALLOCATE (matrix_wx1(ispin)%matrix)
1391 502 : CALL dbcsr_create(matrix=matrix_wx1(ispin)%matrix, template=matrix_s(1)%matrix)
1392 502 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_wx1(ispin)%matrix, sab_orb)
1393 502 : CALL dbcsr_set(matrix_wx1(ispin)%matrix, 0.0_dp)
1394 :
1395 : ! Add Hartree contributions to the perturbation vector
1396 502 : IF (.NOT. (is_rks_triplets .OR. do_sf)) THEN
1397 : CALL cp_dbcsr_sm_fm_multiply(matrix_hx(ispin)%matrix, evect(ispin), &
1398 432 : cpscr, nactive(ispin), alpha=focc, beta=1.0_dp)
1399 : CALL cp_dbcsr_sm_fm_multiply(matrix_hx(ispin)%matrix, mos, vcvec, norb(ispin), &
1400 432 : alpha=1.0_dp, beta=0.0_dp)
1401 : CALL parallel_gemm("T", "N", nactive(ispin), norb(ispin), nao, 1.0_dp, &
1402 432 : mosa, vcvec, 0.0_dp, avcmat)
1403 : CALL parallel_gemm("N", "N", nao, norb(ispin), nactive(ispin), 1.0_dp, &
1404 432 : evect(ispin), avcmat, 0.0_dp, vcvec)
1405 432 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), norb(ispin), alpha=-focc, beta=1.0_dp)
1406 : !
1407 : CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=mos, matrix_g=vcvec, &
1408 432 : ncol=norb(ispin), alpha=2.0_dp, symmetry_mode=1)
1409 : END IF
1410 : ! Add exchange-correlation kernel and exchange-correlation kernel derivative contributions to the response vector
1411 502 : IF ((myfun /= xc_none) .AND. (tddfpt_control%spinflip /= tddfpt_sf_col)) THEN
1412 :
1413 : ! XC Kernel contributions
1414 : ! For spin-flip excitations this is the only contribution to the alpha response vector
1415 342 : IF (.NOT. (do_sf .AND. (ispin == 2))) THEN
1416 : ! F*X
1417 : CALL cp_dbcsr_sm_fm_multiply(matrix_fx(spin)%matrix, evect(spin), &
1418 334 : cpscr, nactive(ispin), alpha=focc, beta=1.0_dp)
1419 : END IF
1420 : ! For spin-flip excitations this is the only contribution to the beta response vector
1421 342 : IF (.NOT. (do_sf .AND. (ispin == 1))) THEN
1422 : ! F*Cb
1423 : CALL cp_dbcsr_sm_fm_multiply(matrix_fx(spin)%matrix, mos2, vcvec, &
1424 334 : norb(ispin), alpha=1.0_dp, beta=0.0_dp)
1425 : ! Ca^T*F*Cb
1426 : CALL parallel_gemm("T", "N", nactive(spin), norb(ispin), nao, 1.0_dp, &
1427 334 : mosa, vcvec, 0.0_dp, avcmat)
1428 : ! X*Ca^T*F*Cb
1429 : CALL parallel_gemm("N", "N", nao, norb(ispin), nactive(spin), 1.0_dp, &
1430 334 : evect(spin), avcmat, 0.0_dp, vcvec)
1431 : ! -S*X*Ca^T*F*Cb
1432 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
1433 334 : norb(ispin), alpha=-focc, beta=1.0_dp)
1434 : ! Add contributions to the \Lambda_munu for the perturbed overlap matrix term, third term of Eq. 51
1435 : ! 2X*Ca^T*F*Cb*Cb^T
1436 : CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=vcvec, matrix_g=gs_mos(ispin)%mos_occ, &
1437 334 : ncol=norb(ispin), alpha=2.0_dp, symmetry_mode=1)
1438 : END IF
1439 : !
1440 :
1441 : ! XC g (third functional derivative) contributions
1442 : ! g*Ca*focc
1443 : CALL cp_dbcsr_sm_fm_multiply(matrix_gx(ispin)%matrix, gs_mos(ispin)%mos_occ, &
1444 342 : cpmos(ispin), norb(ispin), alpha=focc, beta=1.0_dp)
1445 : ! Add contributions to the \Lambda_munu for the perturbed overlap matrix term, third term of Eq. 51
1446 : ! g*Ca
1447 : CALL cp_dbcsr_sm_fm_multiply(matrix_gx(ispin)%matrix, gs_mos(ispin)%mos_occ, vcvec, norb(ispin), &
1448 342 : alpha=1.0_dp, beta=0.0_dp)
1449 : ! Ca^T*g*Ca
1450 342 : CALL parallel_gemm("T", "N", norb(ispin), norb(ispin), nao, 1.0_dp, gs_mos(ispin)%mos_occ, vcvec, 0.0_dp, cvcmat)
1451 : ! Ca*Ca^T*g*Ca
1452 342 : CALL parallel_gemm("N", "N", nao, norb(ispin), norb(ispin), 1.0_dp, gs_mos(ispin)%mos_occ, cvcmat, 0.0_dp, vcvec)
1453 : ! Ca*Ca^T*g*Ca*Ca^T
1454 : CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=vcvec, matrix_g=gs_mos(ispin)%mos_occ, &
1455 342 : ncol=norb(ispin), alpha=1.0_dp, symmetry_mode=1)
1456 : !
1457 : END IF
1458 : ! Add Fock contributions to the response vector
1459 502 : IF (do_hfx) THEN
1460 : ! For spin-flip excitations this is the only contribution to the alpha response vector
1461 166 : IF (.NOT. (do_sf .AND. (ispin == 2))) THEN
1462 : ! F^sym*X
1463 : CALL cp_dbcsr_sm_fm_multiply(matrix_hfx(spin)%matrix, evect(spin), &
1464 162 : cpscr, nactive(spin), alpha=focc, beta=1.0_dp)
1465 : ! F^asym*X
1466 : CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_asymm(spin)%matrix, evect(spin), &
1467 162 : cpscr, nactive(spin), alpha=focc, beta=1.0_dp)
1468 : END IF
1469 :
1470 : ! For spin-flip excitations this is the only contribution to the beta response vector
1471 166 : IF (.NOT. (do_sf .AND. (ispin == 1))) THEN
1472 : ! F^sym*Cb
1473 : CALL cp_dbcsr_sm_fm_multiply(matrix_hfx(spin)%matrix, mos2, vcvec, norb(ispin), &
1474 162 : alpha=1.0_dp, beta=0.0_dp)
1475 : ! -F^asym*Cb
1476 : CALL cp_dbcsr_sm_fm_multiply(matrix_hfx_asymm(spin)%matrix, mos2, vcvec, norb(ispin), &
1477 162 : alpha=1.0_dp, beta=1.0_dp)
1478 : ! Ca^T*F*Cb
1479 : CALL parallel_gemm("T", "N", nactive(spin), norb(ispin), nao, 1.0_dp, &
1480 162 : mosa, vcvec, 0.0_dp, avcmat)
1481 : ! X*Ca^T*F*Cb
1482 : CALL parallel_gemm("N", "N", nao, norb(ispin), nactive(spin), 1.0_dp, &
1483 162 : evect(spin), avcmat, 0.0_dp, vcvec)
1484 : ! -S*X*Ca^T*F*Cb
1485 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
1486 162 : norb(ispin), alpha=-focc, beta=1.0_dp)
1487 : ! Add contributions to the \Lambda_munu for the perturbed overlap matrix term, third term of Eq. 51
1488 : ! 2X*Ca^T*F*Cb*Cb^T
1489 : CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=vcvec, matrix_g=mos2, &
1490 162 : ncol=norb(ispin), alpha=2.0_dp, symmetry_mode=1)
1491 : END IF
1492 : END IF
1493 : !
1494 502 : IF (do_res) THEN
1495 0 : DO ia = 1, nactive(ispin)
1496 0 : ib = gs_mos(ispin)%index_active(ia)
1497 0 : CALL cp_fm_add_columns(cpscr, cpmos(ispin), 1, 1.0_dp, ia, ib)
1498 : END DO
1499 : ELSE
1500 502 : CALL cp_fm_geadd(1.0_dp, "N", cpscr, 1.0_dp, cpmos(ispin))
1501 : END IF
1502 : !
1503 502 : CALL cp_fm_release(cpscr)
1504 502 : CALL cp_fm_release(avamat)
1505 502 : CALL cp_fm_release(avcmat)
1506 502 : CALL cp_fm_release(cvcmat)
1507 502 : CALL cp_fm_release(vcvec)
1508 1424 : CALL cp_fm_release(vavec)
1509 : END DO
1510 :
1511 420 : IF (.NOT. (is_rks_triplets .OR. do_sf)) THEN
1512 362 : CALL dbcsr_deallocate_matrix_set(matrix_hx)
1513 : END IF
1514 420 : IF (ASSOCIATED(ex_env%matrix_wx1)) CALL dbcsr_deallocate_matrix_set(ex_env%matrix_wx1)
1515 420 : ex_env%matrix_wx1 => matrix_wx1
1516 420 : IF (.NOT. ((myfun == xc_none) .OR. (tddfpt_control%spinflip == tddfpt_sf_col))) THEN
1517 298 : CALL dbcsr_deallocate_matrix_set(matrix_fx)
1518 298 : CALL dbcsr_deallocate_matrix_set(matrix_gx)
1519 : END IF
1520 420 : IF (do_hfx) THEN
1521 150 : CALL dbcsr_deallocate_matrix_set(matrix_hfx)
1522 150 : CALL dbcsr_deallocate_matrix_set(matrix_hfx_asymm)
1523 : END IF
1524 :
1525 420 : CALL timestop(handle)
1526 :
1527 840 : END SUBROUTINE fhxc_force
1528 :
1529 : ! **************************************************************************************************
1530 : !> \brief Simplified Tamm Dancoff approach (sTDA). Kernel contribution to forces
1531 : !> \param qs_env ...
1532 : !> \param ex_env ...
1533 : !> \param gs_mos ...
1534 : !> \param stda_env ...
1535 : !> \param sub_env ...
1536 : !> \param work ...
1537 : !> \param debug_forces ...
1538 : ! **************************************************************************************************
1539 160 : SUBROUTINE stda_force(qs_env, ex_env, gs_mos, stda_env, sub_env, work, debug_forces)
1540 :
1541 : TYPE(qs_environment_type), POINTER :: qs_env
1542 : TYPE(excited_energy_type), POINTER :: ex_env
1543 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1544 : POINTER :: gs_mos
1545 : TYPE(stda_env_type), POINTER :: stda_env
1546 : TYPE(tddfpt_subgroup_env_type) :: sub_env
1547 : TYPE(tddfpt_work_matrices) :: work
1548 : LOGICAL, INTENT(IN) :: debug_forces
1549 :
1550 : CHARACTER(len=*), PARAMETER :: routineN = 'stda_force'
1551 :
1552 : INTEGER :: atom_i, atom_j, ewald_type, handle, i, &
1553 : ia, iatom, idimk, ikind, iounit, is, &
1554 : ispin, jatom, jkind, jspin, nao, &
1555 : natom, norb, nsgf, nspins
1556 160 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, first_sgf, kind_of, &
1557 160 : last_sgf
1558 : INTEGER, DIMENSION(2) :: nactive, nlim
1559 : LOGICAL :: calculate_forces, do_coulomb, do_ewald, &
1560 : found, is_rks_triplets, use_virial
1561 : REAL(KIND=dp) :: alpha, bp, dgabr, dr, eta, fdim, gabr, &
1562 : hfx, rbeta, spinfac, xfac
1563 160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tcharge, tv
1564 160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: gtcharge
1565 : REAL(KIND=dp), DIMENSION(3) :: fij, fodeb, rij
1566 160 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: gab, pblock
1567 160 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1568 : TYPE(cell_type), POINTER :: cell
1569 : TYPE(cp_fm_struct_type), POINTER :: fmstruct, fmstruct_mat, fmstructjspin
1570 : TYPE(cp_fm_type) :: cvcmat, cvec, cvecjspin, t0matrix, &
1571 : t1matrix, vcvec, xvec
1572 160 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: xtransformed
1573 160 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos, X
1574 : TYPE(cp_fm_type), POINTER :: ct, ctjspin, ucmatrix, uxmatrix
1575 : TYPE(cp_logger_type), POINTER :: logger
1576 : TYPE(dbcsr_iterator_type) :: iter
1577 160 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: gamma_matrix, matrix_plo, matrix_s, &
1578 160 : matrix_wx1, scrm
1579 : TYPE(dbcsr_type) :: pdens, ptrans
1580 : TYPE(dbcsr_type), POINTER :: tempmat
1581 : TYPE(dft_control_type), POINTER :: dft_control
1582 : TYPE(ewald_environment_type), POINTER :: ewald_env
1583 : TYPE(ewald_pw_type), POINTER :: ewald_pw
1584 : TYPE(mp_para_env_type), POINTER :: para_env
1585 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1586 160 : POINTER :: n_list, sab_orb
1587 160 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1588 160 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1589 160 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1590 : TYPE(qs_ks_env_type), POINTER :: ks_env
1591 : TYPE(stda_control_type), POINTER :: stda_control
1592 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
1593 : TYPE(virial_type), POINTER :: virial
1594 :
1595 160 : CALL timeset(routineN, handle)
1596 :
1597 160 : CPASSERT(ASSOCIATED(ex_env))
1598 160 : CPASSERT(ASSOCIATED(gs_mos))
1599 :
1600 160 : logger => cp_get_default_logger()
1601 160 : IF (logger%para_env%is_source()) THEN
1602 80 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1603 : ELSE
1604 : iounit = -1
1605 : END IF
1606 :
1607 160 : CALL get_qs_env(qs_env, dft_control=dft_control)
1608 160 : tddfpt_control => dft_control%tddfpt2_control
1609 160 : stda_control => tddfpt_control%stda_control
1610 160 : nspins = dft_control%nspins
1611 160 : is_rks_triplets = tddfpt_control%rks_triplets .AND. (nspins == 1)
1612 :
1613 160 : X => ex_env%evect
1614 :
1615 480 : nactive(:) = stda_env%nactive(:)
1616 160 : xfac = 2.0_dp
1617 160 : spinfac = 2.0_dp
1618 160 : IF (nspins == 2) spinfac = 1.0_dp
1619 160 : NULLIFY (matrix_wx1)
1620 160 : CALL dbcsr_allocate_matrix_set(matrix_wx1, nspins)
1621 160 : NULLIFY (matrix_plo)
1622 160 : CALL dbcsr_allocate_matrix_set(matrix_plo, nspins)
1623 :
1624 160 : IF (nspins == 1 .AND. is_rks_triplets) THEN
1625 : do_coulomb = .FALSE.
1626 : ELSE
1627 144 : do_coulomb = .TRUE.
1628 : END IF
1629 160 : do_ewald = stda_control%do_ewald
1630 :
1631 160 : CALL get_qs_env(qs_env, para_env=para_env, force=force)
1632 :
1633 : CALL get_qs_env(qs_env, natom=natom, cell=cell, &
1634 160 : particle_set=particle_set, qs_kind_set=qs_kind_set)
1635 480 : ALLOCATE (first_sgf(natom))
1636 320 : ALLOCATE (last_sgf(natom))
1637 160 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf)
1638 :
1639 160 : CALL get_qs_env(qs_env, ks_env=ks_env, matrix_s=matrix_s, sab_orb=sab_orb, atomic_kind_set=atomic_kind_set)
1640 160 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
1641 :
1642 : ! calculate Loewdin transformed Davidson trial vector tilde(X)=S^1/2*X
1643 : ! and tilde(tilde(X))=S^1/2_A*tilde(X)_A
1644 666 : ALLOCATE (xtransformed(nspins))
1645 346 : DO ispin = 1, nspins
1646 186 : NULLIFY (fmstruct)
1647 186 : ct => work%ctransformed(ispin)
1648 186 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct)
1649 346 : CALL cp_fm_create(matrix=xtransformed(ispin), matrix_struct=fmstruct, name="XTRANSFORMED")
1650 : END DO
1651 160 : CALL get_lowdin_x(work%shalf, X, xtransformed)
1652 :
1653 800 : ALLOCATE (tcharge(natom), gtcharge(natom, 4))
1654 :
1655 160 : cpmos => ex_env%cpmos
1656 :
1657 346 : DO ispin = 1, nspins
1658 186 : ct => work%ctransformed(ispin)
1659 186 : CALL cp_fm_get_info(ct, matrix_struct=fmstruct, nrow_global=nsgf)
1660 558 : ALLOCATE (tv(nsgf))
1661 186 : CALL cp_fm_create(cvec, fmstruct)
1662 186 : CALL cp_fm_create(xvec, fmstruct)
1663 : !
1664 186 : ALLOCATE (matrix_wx1(ispin)%matrix)
1665 186 : CALL dbcsr_create(matrix=matrix_wx1(ispin)%matrix, template=matrix_s(1)%matrix)
1666 186 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_wx1(ispin)%matrix, sab_orb)
1667 186 : CALL dbcsr_set(matrix_wx1(ispin)%matrix, 0.0_dp)
1668 186 : ALLOCATE (matrix_plo(ispin)%matrix)
1669 186 : CALL dbcsr_create(matrix=matrix_plo(ispin)%matrix, template=matrix_s(1)%matrix)
1670 186 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_plo(ispin)%matrix, sab_orb)
1671 186 : CALL dbcsr_set(matrix_plo(ispin)%matrix, 0.0_dp)
1672 : !
1673 : ! *** Coulomb contribution
1674 : !
1675 186 : IF (do_coulomb) THEN
1676 : !
1677 182 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1678 : !
1679 878 : tcharge(:) = 0.0_dp
1680 392 : DO jspin = 1, nspins
1681 222 : ctjspin => work%ctransformed(jspin)
1682 222 : CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin)
1683 222 : CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin, nrow_global=nsgf)
1684 222 : CALL cp_fm_create(cvecjspin, fmstructjspin)
1685 : ! CV(mu,j) = CT(mu,j)*XT(mu,j)
1686 222 : CALL cp_fm_schur_product(ctjspin, xtransformed(jspin), cvecjspin)
1687 : ! TV(mu) = SUM_j CV(mu,j)
1688 222 : CALL cp_fm_vectorssum(cvecjspin, tv, "R")
1689 : ! contract charges
1690 : ! TC(a) = SUM_(mu of a) TV(mu)
1691 1086 : DO ia = 1, natom
1692 5746 : DO is = first_sgf(ia), last_sgf(ia)
1693 5524 : tcharge(ia) = tcharge(ia) + tv(is)
1694 : END DO
1695 : END DO
1696 614 : CALL cp_fm_release(cvecjspin)
1697 : END DO !jspin
1698 : ! Apply tcharge*gab -> gtcharge
1699 : ! gT(b) = SUM_a g(a,b)*TC(a)
1700 : ! gab = work%gamma_exchange(1)%matrix
1701 3682 : gtcharge = 0.0_dp
1702 : ! short range contribution
1703 170 : NULLIFY (gamma_matrix)
1704 170 : CALL setup_gamma(qs_env, stda_env, sub_env, gamma_matrix, ndim=4)
1705 170 : tempmat => gamma_matrix(1)%matrix
1706 170 : CALL dbcsr_iterator_start(iter, tempmat)
1707 5331 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1708 5161 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, gab)
1709 5161 : gtcharge(iatom, 1) = gtcharge(iatom, 1) + gab(1, 1)*tcharge(jatom)
1710 5161 : IF (iatom /= jatom) THEN
1711 4807 : gtcharge(jatom, 1) = gtcharge(jatom, 1) + gab(1, 1)*tcharge(iatom)
1712 : END IF
1713 20814 : DO idimk = 2, 4
1714 15483 : fdim = -1.0_dp
1715 : CALL dbcsr_get_block_p(matrix=gamma_matrix(idimk)%matrix, &
1716 15483 : row=iatom, col=jatom, block=gab, found=found)
1717 20644 : IF (found) THEN
1718 15483 : gtcharge(iatom, idimk) = gtcharge(iatom, idimk) + gab(1, 1)*tcharge(jatom)
1719 15483 : IF (iatom /= jatom) THEN
1720 14421 : gtcharge(jatom, idimk) = gtcharge(jatom, idimk) + fdim*gab(1, 1)*tcharge(iatom)
1721 : END IF
1722 : END IF
1723 : END DO
1724 : END DO
1725 170 : CALL dbcsr_iterator_stop(iter)
1726 170 : CALL dbcsr_deallocate_matrix_set(gamma_matrix)
1727 : ! Ewald long range contribution
1728 170 : IF (do_ewald) THEN
1729 40 : ewald_env => work%ewald_env
1730 40 : ewald_pw => work%ewald_pw
1731 40 : CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
1732 40 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list, virial=virial)
1733 40 : use_virial = .FALSE.
1734 40 : calculate_forces = .FALSE.
1735 40 : CALL tb_ewald_overlap(gtcharge, tcharge, alpha, n_list, virial, use_virial)
1736 : CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
1737 40 : gtcharge, tcharge, calculate_forces, virial, use_virial)
1738 : ! add self charge interaction contribution
1739 40 : IF (para_env%is_source()) THEN
1740 173 : gtcharge(:, 1) = gtcharge(:, 1) - 2._dp*alpha*oorootpi*tcharge(:)
1741 : END IF
1742 : ELSE
1743 130 : nlim = get_limit(natom, para_env%num_pe, para_env%mepos)
1744 331 : DO iatom = nlim(1), nlim(2)
1745 544 : DO jatom = 1, iatom - 1
1746 852 : rij = particle_set(iatom)%r - particle_set(jatom)%r
1747 852 : rij = pbc(rij, cell)
1748 852 : dr = SQRT(SUM(rij(:)**2))
1749 414 : IF (dr > 1.e-6_dp) THEN
1750 213 : gtcharge(iatom, 1) = gtcharge(iatom, 1) + tcharge(jatom)/dr
1751 213 : gtcharge(jatom, 1) = gtcharge(jatom, 1) + tcharge(iatom)/dr
1752 852 : DO idimk = 2, 4
1753 639 : gtcharge(iatom, idimk) = gtcharge(iatom, idimk) + rij(idimk - 1)*tcharge(jatom)/dr**3
1754 852 : gtcharge(jatom, idimk) = gtcharge(jatom, idimk) - rij(idimk - 1)*tcharge(iatom)/dr**3
1755 : END DO
1756 : END IF
1757 : END DO
1758 : END DO
1759 : END IF
1760 170 : CALL para_env%sum(gtcharge(:, 1))
1761 : ! expand charges
1762 : ! TV(mu) = TC(a of mu)
1763 3946 : tv(1:nsgf) = 0.0_dp
1764 878 : DO ia = 1, natom
1765 4654 : DO is = first_sgf(ia), last_sgf(ia)
1766 4484 : tv(is) = gtcharge(ia, 1)
1767 : END DO
1768 : END DO
1769 : !
1770 878 : DO iatom = 1, natom
1771 708 : ikind = kind_of(iatom)
1772 708 : atom_i = atom_of_kind(iatom)
1773 2832 : DO i = 1, 3
1774 2832 : fij(i) = spinfac*spinfac*gtcharge(iatom, i + 1)*tcharge(iatom)
1775 : END DO
1776 708 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
1777 708 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
1778 878 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
1779 : END DO
1780 : !
1781 170 : IF (debug_forces) THEN
1782 16 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1783 4 : CALL para_env%sum(fodeb)
1784 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Coul[X] ", fodeb
1785 : END IF
1786 170 : norb = nactive(ispin)
1787 : ! forces from Lowdin charge derivative
1788 170 : CALL cp_fm_get_info(work%S_C0_C0T(ispin), matrix_struct=fmstruct)
1789 170 : CALL cp_fm_create(t0matrix, matrix_struct=fmstruct, name="T0 SCRATCH")
1790 170 : CALL cp_fm_create(t1matrix, matrix_struct=fmstruct, name="T1 SCRATCH")
1791 170 : ALLOCATE (ucmatrix)
1792 170 : CALL fm_pool_create_fm(work%fm_pool_ao_mo_active(ispin)%pool, ucmatrix)
1793 170 : ALLOCATE (uxmatrix)
1794 170 : CALL fm_pool_create_fm(work%fm_pool_ao_mo_active(ispin)%pool, uxmatrix)
1795 170 : ct => work%ctransformed(ispin)
1796 170 : CALL cp_fm_to_fm(ct, cvec)
1797 170 : CALL cp_fm_row_scale(cvec, tv)
1798 : CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1799 170 : cvec, 0.0_dp, ucmatrix)
1800 : CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1801 170 : X(ispin), 0.0_dp, uxmatrix)
1802 170 : CALL parallel_gemm('N', 'T', nsgf, nsgf, norb, 1.0_dp, uxmatrix, ucmatrix, 0.0_dp, t0matrix)
1803 170 : CALL cp_fm_to_fm(xtransformed(ispin), cvec)
1804 170 : CALL cp_fm_row_scale(cvec, tv)
1805 : CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1806 170 : cvec, 0.0_dp, uxmatrix)
1807 : CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1808 170 : gs_mos(ispin)%mos_occ, 0.0_dp, ucmatrix)
1809 170 : CALL parallel_gemm('N', 'T', nsgf, nsgf, norb, 1.0_dp, ucmatrix, uxmatrix, 1.0_dp, t0matrix)
1810 170 : CALL cp_fm_schur_product(work%slambda, t0matrix, t1matrix)
1811 : !
1812 : CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, spinfac, work%S_eigenvectors, t1matrix, &
1813 170 : 0.0_dp, t0matrix)
1814 : CALL cp_dbcsr_plus_fm_fm_t(matrix_plo(ispin)%matrix, matrix_v=t0matrix, &
1815 170 : matrix_g=work%S_eigenvectors, ncol=nsgf, alpha=2.0_dp, symmetry_mode=1)
1816 170 : CALL fm_pool_give_back_fm(work%fm_pool_ao_mo_active(ispin)%pool, ucmatrix)
1817 170 : DEALLOCATE (ucmatrix)
1818 170 : CALL fm_pool_give_back_fm(work%fm_pool_ao_mo_active(ispin)%pool, uxmatrix)
1819 170 : DEALLOCATE (uxmatrix)
1820 170 : CALL cp_fm_release(t0matrix)
1821 170 : CALL cp_fm_release(t1matrix)
1822 : !
1823 : ! CV(mu,i) = TV(mu)*XT(mu,i)
1824 170 : CALL cp_fm_to_fm(xtransformed(ispin), cvec)
1825 170 : CALL cp_fm_row_scale(cvec, tv)
1826 170 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, cpmos(ispin), norb, 2.0_dp*spinfac, 1.0_dp)
1827 : ! CV(mu,i) = TV(mu)*CT(mu,i)
1828 170 : ct => work%ctransformed(ispin)
1829 170 : CALL cp_fm_to_fm(ct, cvec)
1830 170 : CALL cp_fm_row_scale(cvec, tv)
1831 : ! Shalf(nu,mu)*CV(mu,i)
1832 170 : CALL cp_fm_get_info(cvec, matrix_struct=fmstruct, nrow_global=nao)
1833 170 : CALL cp_fm_create(vcvec, fmstruct)
1834 170 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, vcvec, norb, 1.0_dp, 0.0_dp)
1835 : CALL cp_fm_struct_create(fmstruct_mat, context=fmstruct%context, nrow_global=norb, &
1836 170 : ncol_global=norb, para_env=fmstruct%para_env)
1837 170 : CALL cp_fm_create(cvcmat, fmstruct_mat)
1838 170 : CALL cp_fm_struct_release(fmstruct_mat)
1839 170 : CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, gs_mos(ispin)%mos_occ, vcvec, 0.0_dp, cvcmat)
1840 170 : CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, X(ispin), cvcmat, 0.0_dp, vcvec)
1841 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
1842 170 : nactive(ispin), alpha=-2.0_dp*spinfac, beta=1.0_dp)
1843 : ! wx1
1844 170 : alpha = 2.0_dp
1845 : CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=gs_mos(ispin)%mos_occ, &
1846 170 : matrix_g=vcvec, ncol=norb, alpha=2.0_dp*alpha, symmetry_mode=1)
1847 170 : CALL cp_fm_release(vcvec)
1848 170 : CALL cp_fm_release(cvcmat)
1849 : END IF
1850 : !
1851 : ! *** Exchange contribution
1852 : !
1853 186 : IF (stda_env%do_exchange) THEN
1854 : !
1855 174 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1856 : !
1857 162 : norb = nactive(ispin)
1858 : !
1859 162 : tempmat => work%shalf
1860 162 : CALL dbcsr_create(pdens, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
1861 : ! P(nu,mu) = SUM_j XT(nu,j)*CT(mu,j)
1862 162 : ct => work%ctransformed(ispin)
1863 162 : CALL dbcsr_set(pdens, 0.0_dp)
1864 : CALL cp_dbcsr_plus_fm_fm_t(pdens, xtransformed(ispin), ct, nactive(ispin), &
1865 162 : 1.0_dp, keep_sparsity=.FALSE.)
1866 162 : CALL dbcsr_filter(pdens, stda_env%eps_td_filter)
1867 : ! Apply PP*gab -> PP; gab = gamma_coulomb
1868 : ! P(nu,mu) = P(nu,mu)*g(a of nu,b of mu)
1869 162 : bp = stda_env%beta_param
1870 162 : hfx = stda_env%hfx_fraction
1871 162 : CALL dbcsr_iterator_start(iter, pdens)
1872 10077 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1873 9915 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock)
1874 39660 : rij = particle_set(iatom)%r - particle_set(jatom)%r
1875 39660 : rij = pbc(rij, cell)
1876 39660 : dr = SQRT(SUM(rij(:)**2))
1877 9915 : ikind = kind_of(iatom)
1878 9915 : jkind = kind_of(jatom)
1879 : eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
1880 9915 : stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
1881 9915 : rbeta = dr**bp
1882 9915 : IF (hfx > 0.0_dp) THEN
1883 9847 : gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
1884 : ELSE
1885 68 : IF (dr < 1.0e-6_dp) THEN
1886 : gabr = 0.0_dp
1887 : ELSE
1888 48 : gabr = 1._dp/dr
1889 : END IF
1890 : END IF
1891 : ! gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
1892 : ! forces
1893 9895 : IF (dr > 1.0e-6_dp) THEN
1894 9592 : IF (hfx > 0.0_dp) THEN
1895 9544 : dgabr = -(1._dp/bp)*(1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp + 1._dp)
1896 9544 : dgabr = bp*rbeta/dr**2*dgabr
1897 112572 : dgabr = SUM(pblock**2)*dgabr
1898 : ELSE
1899 48 : dgabr = -1.0_dp/dr**3
1900 3504 : dgabr = SUM(pblock**2)*dgabr
1901 : END IF
1902 9592 : atom_i = atom_of_kind(iatom)
1903 9592 : atom_j = atom_of_kind(jatom)
1904 38368 : DO i = 1, 3
1905 38368 : fij(i) = dgabr*rij(i)
1906 : END DO
1907 9592 : force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
1908 9592 : force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
1909 9592 : force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
1910 9592 : force(jkind)%rho_elec(1, atom_j) = force(jkind)%rho_elec(1, atom_j) + fij(1)
1911 9592 : force(jkind)%rho_elec(2, atom_j) = force(jkind)%rho_elec(2, atom_j) + fij(2)
1912 9592 : force(jkind)%rho_elec(3, atom_j) = force(jkind)%rho_elec(3, atom_j) + fij(3)
1913 : END IF
1914 : !
1915 133487 : pblock = gabr*pblock
1916 : END DO
1917 162 : CALL dbcsr_iterator_stop(iter)
1918 : !
1919 : ! Transpose pdens matrix
1920 162 : CALL dbcsr_create(ptrans, template=pdens)
1921 162 : CALL dbcsr_transposed(ptrans, pdens)
1922 : !
1923 : ! forces from Lowdin charge derivative
1924 162 : CALL cp_fm_get_info(work%S_C0_C0T(ispin), matrix_struct=fmstruct)
1925 162 : CALL cp_fm_create(t0matrix, matrix_struct=fmstruct, name="T0 SCRATCH")
1926 162 : CALL cp_fm_create(t1matrix, matrix_struct=fmstruct, name="T1 SCRATCH")
1927 162 : ALLOCATE (ucmatrix)
1928 162 : CALL fm_pool_create_fm(work%fm_pool_ao_mo_active(ispin)%pool, ucmatrix)
1929 162 : ALLOCATE (uxmatrix)
1930 162 : CALL fm_pool_create_fm(work%fm_pool_ao_mo_active(ispin)%pool, uxmatrix)
1931 162 : ct => work%ctransformed(ispin)
1932 162 : CALL cp_dbcsr_sm_fm_multiply(pdens, ct, cvec, norb, 1.0_dp, 0.0_dp)
1933 : CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1934 162 : cvec, 0.0_dp, ucmatrix)
1935 : CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1936 162 : X(ispin), 0.0_dp, uxmatrix)
1937 162 : CALL parallel_gemm('N', 'T', nsgf, nsgf, norb, 1.0_dp, uxmatrix, ucmatrix, 0.0_dp, t0matrix)
1938 162 : CALL cp_dbcsr_sm_fm_multiply(ptrans, xtransformed(ispin), cvec, norb, 1.0_dp, 0.0_dp)
1939 : CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1940 162 : cvec, 0.0_dp, uxmatrix)
1941 : CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1942 162 : gs_mos(ispin)%mos_occ, 0.0_dp, ucmatrix)
1943 162 : CALL parallel_gemm('N', 'T', nsgf, nsgf, norb, 1.0_dp, ucmatrix, uxmatrix, 1.0_dp, t0matrix)
1944 162 : CALL cp_fm_schur_product(work%slambda, t0matrix, t1matrix)
1945 : CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, -1.0_dp, work%S_eigenvectors, t1matrix, &
1946 162 : 0.0_dp, t0matrix)
1947 : CALL cp_dbcsr_plus_fm_fm_t(matrix_plo(ispin)%matrix, matrix_v=t0matrix, &
1948 162 : matrix_g=work%S_eigenvectors, ncol=nsgf, alpha=2.0_dp, symmetry_mode=1)
1949 162 : CALL fm_pool_give_back_fm(work%fm_pool_ao_mo_active(ispin)%pool, ucmatrix)
1950 162 : DEALLOCATE (ucmatrix)
1951 162 : CALL fm_pool_give_back_fm(work%fm_pool_ao_mo_active(ispin)%pool, uxmatrix)
1952 162 : DEALLOCATE (uxmatrix)
1953 162 : CALL cp_fm_release(t0matrix)
1954 162 : CALL cp_fm_release(t1matrix)
1955 :
1956 : ! RHS contribution to response matrix
1957 : ! CV(nu,i) = P(nu,mu)*XT(mu,i)
1958 162 : CALL cp_dbcsr_sm_fm_multiply(ptrans, xtransformed(ispin), cvec, norb, 1.0_dp, 0.0_dp)
1959 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, cpmos(ispin), norb, &
1960 162 : alpha=-xfac, beta=1.0_dp)
1961 : !
1962 162 : CALL cp_fm_get_info(cvec, matrix_struct=fmstruct, nrow_global=nao)
1963 162 : CALL cp_fm_create(vcvec, fmstruct)
1964 : ! CV(nu,i) = P(nu,mu)*CT(mu,i)
1965 162 : CALL cp_dbcsr_sm_fm_multiply(ptrans, ct, cvec, norb, 1.0_dp, 0.0_dp)
1966 162 : CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, vcvec, norb, 1.0_dp, 0.0_dp)
1967 : CALL cp_fm_struct_create(fmstruct_mat, context=fmstruct%context, nrow_global=norb, &
1968 162 : ncol_global=norb, para_env=fmstruct%para_env)
1969 162 : CALL cp_fm_create(cvcmat, fmstruct_mat)
1970 162 : CALL cp_fm_struct_release(fmstruct_mat)
1971 162 : CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, gs_mos(ispin)%mos_occ, vcvec, 0.0_dp, cvcmat)
1972 162 : CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, X(ispin), cvcmat, 0.0_dp, vcvec)
1973 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
1974 162 : norb, alpha=xfac, beta=1.0_dp)
1975 : ! wx1
1976 162 : IF (nspins == 2) THEN
1977 44 : alpha = -2.0_dp
1978 : ELSE
1979 118 : alpha = -1.0_dp
1980 : END IF
1981 : CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=gs_mos(ispin)%mos_occ, &
1982 : matrix_g=vcvec, &
1983 162 : ncol=norb, alpha=2.0_dp*alpha, symmetry_mode=1)
1984 162 : CALL cp_fm_release(vcvec)
1985 162 : CALL cp_fm_release(cvcmat)
1986 : !
1987 162 : CALL dbcsr_release(pdens)
1988 162 : CALL dbcsr_release(ptrans)
1989 : !
1990 162 : IF (debug_forces) THEN
1991 16 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1992 4 : CALL para_env%sum(fodeb)
1993 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Exch[X] ", fodeb
1994 : END IF
1995 : END IF
1996 : !
1997 186 : CALL cp_fm_release(cvec)
1998 186 : CALL cp_fm_release(xvec)
1999 718 : DEALLOCATE (tv)
2000 : END DO
2001 :
2002 160 : CALL cp_fm_release(xtransformed)
2003 160 : DEALLOCATE (tcharge, gtcharge)
2004 160 : DEALLOCATE (first_sgf, last_sgf)
2005 :
2006 : ! Lowdin forces
2007 160 : IF (nspins == 2) THEN
2008 : CALL dbcsr_add(matrix_plo(1)%matrix, matrix_plo(2)%matrix, &
2009 26 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2010 : END IF
2011 160 : CALL dbcsr_scale(matrix_plo(1)%matrix, -1.0_dp)
2012 160 : NULLIFY (scrm)
2013 172 : IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
2014 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
2015 : matrix_name="OVERLAP MATRIX", &
2016 : basis_type_a="ORB", basis_type_b="ORB", &
2017 : sab_nl=sab_orb, calculate_forces=.TRUE., &
2018 160 : matrix_p=matrix_plo(1)%matrix)
2019 160 : CALL dbcsr_deallocate_matrix_set(scrm)
2020 160 : CALL dbcsr_deallocate_matrix_set(matrix_plo)
2021 160 : IF (debug_forces) THEN
2022 16 : fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
2023 4 : CALL para_env%sum(fodeb)
2024 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Lowdin ", fodeb
2025 : END IF
2026 :
2027 160 : IF (ASSOCIATED(ex_env%matrix_wx1)) CALL dbcsr_deallocate_matrix_set(ex_env%matrix_wx1)
2028 160 : ex_env%matrix_wx1 => matrix_wx1
2029 :
2030 160 : CALL timestop(handle)
2031 :
2032 320 : END SUBROUTINE stda_force
2033 :
2034 : ! **************************************************************************************************
2035 :
2036 : END MODULE qs_tddfpt2_fhxc_forces
|
