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