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