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