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