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