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