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