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