(git:ccc2433)
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-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 MODULE qs_tddfpt2_fhxc_forces
9  USE admm_methods, ONLY: admm_projection_derivative
10  USE admm_types, ONLY: admm_type,&
11  get_admm_env
12  USE atomic_kind_types, ONLY: atomic_kind_type,&
13  get_atomic_kind_set
14  USE cell_types, ONLY: cell_type,&
15  pbc
16  USE cp_control_types, ONLY: dft_control_type,&
17  stda_control_type,&
18  tddfpt2_control_type
19  USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
20  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
21  copy_fm_to_dbcsr,&
22  cp_dbcsr_plus_fm_fm_t,&
23  cp_dbcsr_sm_fm_multiply,&
24  dbcsr_allocate_matrix_set,&
25  dbcsr_deallocate_matrix_set
26  USE cp_fm_basic_linalg, ONLY: cp_fm_row_scale,&
27  cp_fm_schur_product
28  USE cp_fm_pool_types, ONLY: fm_pool_create_fm,&
29  fm_pool_give_back_fm
30  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
31  cp_fm_struct_release,&
32  cp_fm_struct_type
33  USE cp_fm_types, ONLY: cp_fm_create,&
34  cp_fm_get_info,&
35  cp_fm_release,&
36  cp_fm_to_fm,&
37  cp_fm_type,&
38  cp_fm_vectorssum
39  USE cp_log_handling, ONLY: cp_get_default_logger,&
40  cp_logger_get_default_unit_nr,&
41  cp_logger_type
42  USE dbcsr_api, ONLY: &
43  dbcsr_add, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_filter, &
44  dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
45  dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, &
46  dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type, &
47  dbcsr_type_antisymmetric, dbcsr_type_no_symmetry
48  USE ewald_environment_types, ONLY: ewald_env_get,&
49  ewald_environment_type
50  USE ewald_methods_tb, ONLY: tb_ewald_overlap,&
51  tb_spme_evaluate
52  USE ewald_pw_types, ONLY: ewald_pw_type
53  USE exstates_types, ONLY: excited_energy_type
54  USE hartree_local_methods, ONLY: vh_1c_gg_integrals,&
55  init_coulomb_local
56  USE hartree_local_types, ONLY: hartree_local_create,&
57  hartree_local_release,&
58  hartree_local_type
59  USE hfx_derivatives, ONLY: derivatives_four_center
60  USE hfx_energy_potential, ONLY: integrate_four_center
61  USE hfx_ri, ONLY: hfx_ri_update_forces,&
62  hfx_ri_update_ks
63  USE hfx_types, ONLY: hfx_type
64  USE input_constants, ONLY: do_admm_aux_exch_func_none,&
65  tddfpt_kernel_full,&
66  xc_kernel_method_analytic,&
67  xc_kernel_method_best,&
68  xc_kernel_method_numeric,&
69  xc_none
70  USE input_section_types, ONLY: section_vals_get,&
71  section_vals_get_subs_vals,&
72  section_vals_type,&
73  section_vals_val_get
74  USE kinds, ONLY: default_string_length,&
75  dp
76  USE mathconstants, ONLY: oorootpi
77  USE message_passing, ONLY: mp_para_env_type
78  USE parallel_gemm_api, ONLY: parallel_gemm
79  USE particle_methods, ONLY: get_particle_set
80  USE particle_types, ONLY: particle_type
81  USE pw_env_types, ONLY: pw_env_get,&
82  pw_env_type
83  USE pw_methods, ONLY: pw_axpy,&
84  pw_scale,&
85  pw_transfer,&
86  pw_zero
87  USE pw_poisson_methods, ONLY: pw_poisson_solve
88  USE pw_poisson_types, ONLY: pw_poisson_type
89  USE pw_pool_types, ONLY: pw_pool_type
90  USE pw_types, ONLY: pw_c1d_gs_type,&
91  pw_r3d_rs_type
92  USE qs_collocate_density, ONLY: calculate_rho_elec
93  USE qs_environment_types, ONLY: get_qs_env,&
94  qs_environment_type,&
95  set_qs_env
96  USE qs_force_types, ONLY: qs_force_type
97  USE qs_fxc, ONLY: qs_fgxc_create,&
98  qs_fgxc_gdiff,&
99  qs_fgxc_release
100  USE qs_gapw_densities, ONLY: prepare_gapw_den
101  USE qs_integrate_potential, ONLY: integrate_v_rspace
102  USE qs_kernel_types, ONLY: full_kernel_env_type
103  USE qs_kind_types, ONLY: qs_kind_type
104  USE qs_ks_atom, ONLY: update_ks_atom
105  USE qs_ks_types, ONLY: qs_ks_env_type
106  USE qs_local_rho_types, ONLY: local_rho_set_create,&
107  local_rho_set_release,&
108  local_rho_type
109  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
110  USE qs_oce_methods, ONLY: build_oce_matrices
111  USE qs_oce_types, ONLY: allocate_oce_set,&
112  create_oce_set,&
113  oce_matrix_type
114  USE qs_overlap, ONLY: build_overlap_matrix
115  USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace,&
116  rho0_s_grid_create
117  USE qs_rho0_methods, ONLY: init_rho0
118  USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals,&
119  calculate_rho_atom_coeff
120  USE qs_rho_atom_types, ONLY: rho_atom_type
121  USE qs_rho_types, ONLY: qs_rho_create,&
122  qs_rho_get,&
123  qs_rho_set,&
124  qs_rho_type
125  USE qs_tddfpt2_stda_types, ONLY: stda_env_type
126  USE qs_tddfpt2_stda_utils, ONLY: get_lowdin_x,&
127  setup_gamma
128  USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
129  USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos,&
130  tddfpt_work_matrices
131  USE qs_vxc_atom, ONLY: calculate_gfxc_atom,&
132  gfxc_atom_diff
133  USE task_list_types, ONLY: task_list_type
134  USE util, ONLY: get_limit
135  USE virial_types, ONLY: virial_type
136 #include "./base/base_uses.f90"
137 
138  IMPLICIT NONE
139 
140  PRIVATE
141 
142  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_fhxc_forces'
143 
144  PUBLIC :: fhxc_force, stda_force
145 
146 ! **************************************************************************************************
147 
148 CONTAINS
149 
150 ! **************************************************************************************************
151 !> \brief Calculate direct tddft forces
152 !> \param qs_env ...
153 !> \param ex_env ...
154 !> \param gs_mos ...
155 !> \param full_kernel ...
156 !> \param debug_forces ...
157 !> \par History
158 !> * 01.2020 screated [JGH]
159 ! **************************************************************************************************
160  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, blk, ewald_type, handle, i, ia, iatom, idimk, ikind, iounit, is, &
1289  ispin, jatom, jkind, jspin, nao, natom, norb, nsgf, nspins
1290  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, first_sgf, kind_of, &
1291  last_sgf
1292  INTEGER, DIMENSION(2) :: nactive, nlim
1293  LOGICAL :: calculate_forces, do_coulomb, do_ewald, &
1294  found, is_rks_triplets, use_virial
1295  REAL(kind=dp) :: alpha, bp, dgabr, dr, eta, fdim, gabr, &
1296  hfx, rbeta, spinfac, xfac
1297  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tcharge, tv
1298  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: gtcharge
1299  REAL(kind=dp), DIMENSION(3) :: fij, fodeb, rij
1300  REAL(kind=dp), DIMENSION(:, :), POINTER :: gab, pblock
1301  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1302  TYPE(cell_type), POINTER :: cell
1303  TYPE(cp_fm_struct_type), POINTER :: fmstruct, fmstruct_mat, fmstructjspin
1304  TYPE(cp_fm_type) :: cvcmat, cvec, cvecjspin, t0matrix, &
1305  t1matrix, vcvec, xvec
1306  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: xtransformed
1307  TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos, x
1308  TYPE(cp_fm_type), POINTER :: ct, ctjspin, ucmatrix, uxmatrix
1309  TYPE(cp_logger_type), POINTER :: logger
1310  TYPE(dbcsr_iterator_type) :: iter
1311  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: gamma_matrix, matrix_plo, matrix_s, &
1312  matrix_wx1, scrm
1313  TYPE(dbcsr_type) :: pdens, ptrans
1314  TYPE(dbcsr_type), POINTER :: tempmat
1315  TYPE(dft_control_type), POINTER :: dft_control
1316  TYPE(ewald_environment_type), POINTER :: ewald_env
1317  TYPE(ewald_pw_type), POINTER :: ewald_pw
1318  TYPE(mp_para_env_type), POINTER :: para_env
1319  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1320  POINTER :: n_list, sab_orb
1321  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1322  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1323  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1324  TYPE(qs_ks_env_type), POINTER :: ks_env
1325  TYPE(stda_control_type), POINTER :: stda_control
1326  TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
1327  TYPE(virial_type), POINTER :: virial
1328 
1329  CALL timeset(routinen, handle)
1330 
1331  cpassert(ASSOCIATED(ex_env))
1332  cpassert(ASSOCIATED(gs_mos))
1333 
1334  logger => cp_get_default_logger()
1335  IF (logger%para_env%is_source()) THEN
1336  iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
1337  ELSE
1338  iounit = -1
1339  END IF
1340 
1341  CALL get_qs_env(qs_env, dft_control=dft_control)
1342  tddfpt_control => dft_control%tddfpt2_control
1343  stda_control => tddfpt_control%stda_control
1344  nspins = dft_control%nspins
1345  is_rks_triplets = tddfpt_control%rks_triplets .AND. (nspins == 1)
1346 
1347  x => ex_env%evect
1348 
1349  nactive(:) = stda_env%nactive(:)
1350  xfac = 2.0_dp
1351  spinfac = 2.0_dp
1352  IF (nspins == 2) spinfac = 1.0_dp
1353  NULLIFY (matrix_wx1)
1354  CALL dbcsr_allocate_matrix_set(matrix_wx1, nspins)
1355  NULLIFY (matrix_plo)
1356  CALL dbcsr_allocate_matrix_set(matrix_plo, nspins)
1357 
1358  IF (nspins == 1 .AND. is_rks_triplets) THEN
1359  do_coulomb = .false.
1360  ELSE
1361  do_coulomb = .true.
1362  END IF
1363  do_ewald = stda_control%do_ewald
1364 
1365  CALL get_qs_env(qs_env, para_env=para_env, force=force)
1366 
1367  CALL get_qs_env(qs_env, natom=natom, cell=cell, &
1368  particle_set=particle_set, qs_kind_set=qs_kind_set)
1369  ALLOCATE (first_sgf(natom))
1370  ALLOCATE (last_sgf(natom))
1371  CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf)
1372 
1373  CALL get_qs_env(qs_env, ks_env=ks_env, matrix_s=matrix_s, sab_orb=sab_orb, atomic_kind_set=atomic_kind_set)
1374  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
1375 
1376  ! calculate Loewdin transformed Davidson trial vector tilde(X)=S^1/2*X
1377  ! and tilde(tilde(X))=S^1/2_A*tilde(X)_A
1378  ALLOCATE (xtransformed(nspins))
1379  DO ispin = 1, nspins
1380  NULLIFY (fmstruct)
1381  ct => work%ctransformed(ispin)
1382  CALL cp_fm_get_info(ct, matrix_struct=fmstruct)
1383  CALL cp_fm_create(matrix=xtransformed(ispin), matrix_struct=fmstruct, name="XTRANSFORMED")
1384  END DO
1385  CALL get_lowdin_x(work%shalf, x, xtransformed)
1386 
1387  ALLOCATE (tcharge(natom), gtcharge(natom, 4))
1388 
1389  cpmos => ex_env%cpmos
1390 
1391  DO ispin = 1, nspins
1392  ct => work%ctransformed(ispin)
1393  CALL cp_fm_get_info(ct, matrix_struct=fmstruct, nrow_global=nsgf)
1394  ALLOCATE (tv(nsgf))
1395  CALL cp_fm_create(cvec, fmstruct)
1396  CALL cp_fm_create(xvec, fmstruct)
1397  !
1398  ALLOCATE (matrix_wx1(ispin)%matrix)
1399  CALL dbcsr_create(matrix=matrix_wx1(ispin)%matrix, template=matrix_s(1)%matrix)
1400  CALL cp_dbcsr_alloc_block_from_nbl(matrix_wx1(ispin)%matrix, sab_orb)
1401  CALL dbcsr_set(matrix_wx1(ispin)%matrix, 0.0_dp)
1402  ALLOCATE (matrix_plo(ispin)%matrix)
1403  CALL dbcsr_create(matrix=matrix_plo(ispin)%matrix, template=matrix_s(1)%matrix)
1404  CALL cp_dbcsr_alloc_block_from_nbl(matrix_plo(ispin)%matrix, sab_orb)
1405  CALL dbcsr_set(matrix_plo(ispin)%matrix, 0.0_dp)
1406  !
1407  ! *** Coulomb contribution
1408  !
1409  IF (do_coulomb) THEN
1410  !
1411  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1412  !
1413  tcharge(:) = 0.0_dp
1414  DO jspin = 1, nspins
1415  ctjspin => work%ctransformed(jspin)
1416  CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin)
1417  CALL cp_fm_get_info(ctjspin, matrix_struct=fmstructjspin, nrow_global=nsgf)
1418  CALL cp_fm_create(cvecjspin, fmstructjspin)
1419  ! CV(mu,j) = CT(mu,j)*XT(mu,j)
1420  CALL cp_fm_schur_product(ctjspin, xtransformed(jspin), cvecjspin)
1421  ! TV(mu) = SUM_j CV(mu,j)
1422  CALL cp_fm_vectorssum(cvecjspin, tv, "R")
1423  ! contract charges
1424  ! TC(a) = SUM_(mu of a) TV(mu)
1425  DO ia = 1, natom
1426  DO is = first_sgf(ia), last_sgf(ia)
1427  tcharge(ia) = tcharge(ia) + tv(is)
1428  END DO
1429  END DO
1430  CALL cp_fm_release(cvecjspin)
1431  END DO !jspin
1432  ! Apply tcharge*gab -> gtcharge
1433  ! gT(b) = SUM_a g(a,b)*TC(a)
1434  ! gab = work%gamma_exchange(1)%matrix
1435  gtcharge = 0.0_dp
1436  ! short range contribution
1437  NULLIFY (gamma_matrix)
1438  CALL setup_gamma(qs_env, stda_env, sub_env, gamma_matrix, ndim=4)
1439  tempmat => gamma_matrix(1)%matrix
1440  CALL dbcsr_iterator_start(iter, tempmat)
1441  DO WHILE (dbcsr_iterator_blocks_left(iter))
1442  CALL dbcsr_iterator_next_block(iter, iatom, jatom, gab, blk)
1443  gtcharge(iatom, 1) = gtcharge(iatom, 1) + gab(1, 1)*tcharge(jatom)
1444  IF (iatom /= jatom) THEN
1445  gtcharge(jatom, 1) = gtcharge(jatom, 1) + gab(1, 1)*tcharge(iatom)
1446  END IF
1447  DO idimk = 2, 4
1448  fdim = -1.0_dp
1449  CALL dbcsr_get_block_p(matrix=gamma_matrix(idimk)%matrix, &
1450  row=iatom, col=jatom, block=gab, found=found)
1451  IF (found) THEN
1452  gtcharge(iatom, idimk) = gtcharge(iatom, idimk) + gab(1, 1)*tcharge(jatom)
1453  IF (iatom /= jatom) THEN
1454  gtcharge(jatom, idimk) = gtcharge(jatom, idimk) + fdim*gab(1, 1)*tcharge(iatom)
1455  END IF
1456  END IF
1457  END DO
1458  END DO
1459  CALL dbcsr_iterator_stop(iter)
1460  CALL dbcsr_deallocate_matrix_set(gamma_matrix)
1461  ! Ewald long range contribution
1462  IF (do_ewald) THEN
1463  ewald_env => work%ewald_env
1464  ewald_pw => work%ewald_pw
1465  CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
1466  CALL get_qs_env(qs_env=qs_env, sab_orb=n_list, virial=virial)
1467  use_virial = .false.
1468  calculate_forces = .false.
1469  CALL tb_ewald_overlap(gtcharge, tcharge, alpha, n_list, virial, use_virial)
1470  CALL tb_spme_evaluate(ewald_env, ewald_pw, particle_set, cell, &
1471  gtcharge, tcharge, calculate_forces, virial, use_virial)
1472  ! add self charge interaction contribution
1473  IF (para_env%is_source()) THEN
1474  gtcharge(:, 1) = gtcharge(:, 1) - 2._dp*alpha*oorootpi*tcharge(:)
1475  END IF
1476  ELSE
1477  nlim = get_limit(natom, para_env%num_pe, para_env%mepos)
1478  DO iatom = nlim(1), nlim(2)
1479  DO jatom = 1, iatom - 1
1480  rij = particle_set(iatom)%r - particle_set(jatom)%r
1481  rij = pbc(rij, cell)
1482  dr = sqrt(sum(rij(:)**2))
1483  IF (dr > 1.e-6_dp) THEN
1484  gtcharge(iatom, 1) = gtcharge(iatom, 1) + tcharge(jatom)/dr
1485  gtcharge(jatom, 1) = gtcharge(jatom, 1) + tcharge(iatom)/dr
1486  DO idimk = 2, 4
1487  gtcharge(iatom, idimk) = gtcharge(iatom, idimk) + rij(idimk - 1)*tcharge(jatom)/dr**3
1488  gtcharge(jatom, idimk) = gtcharge(jatom, idimk) - rij(idimk - 1)*tcharge(iatom)/dr**3
1489  END DO
1490  END IF
1491  END DO
1492  END DO
1493  END IF
1494  CALL para_env%sum(gtcharge(:, 1))
1495  ! expand charges
1496  ! TV(mu) = TC(a of mu)
1497  tv(1:nsgf) = 0.0_dp
1498  DO ia = 1, natom
1499  DO is = first_sgf(ia), last_sgf(ia)
1500  tv(is) = gtcharge(ia, 1)
1501  END DO
1502  END DO
1503  !
1504  DO iatom = 1, natom
1505  ikind = kind_of(iatom)
1506  atom_i = atom_of_kind(iatom)
1507  DO i = 1, 3
1508  fij(i) = spinfac*spinfac*gtcharge(iatom, i + 1)*tcharge(iatom)
1509  END DO
1510  force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
1511  force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
1512  force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
1513  END DO
1514  !
1515  IF (debug_forces) THEN
1516  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1517  CALL para_env%sum(fodeb)
1518  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Coul[X] ", fodeb
1519  END IF
1520  norb = nactive(ispin)
1521  ! forces from Lowdin charge derivative
1522  CALL cp_fm_get_info(work%S_C0_C0T(ispin), matrix_struct=fmstruct)
1523  CALL cp_fm_create(t0matrix, matrix_struct=fmstruct, name="T0 SCRATCH")
1524  CALL cp_fm_create(t1matrix, matrix_struct=fmstruct, name="T1 SCRATCH")
1525  ALLOCATE (ucmatrix)
1526  CALL fm_pool_create_fm(work%fm_pool_ao_mo_occ(ispin)%pool, ucmatrix)
1527  ALLOCATE (uxmatrix)
1528  CALL fm_pool_create_fm(work%fm_pool_ao_mo_occ(ispin)%pool, uxmatrix)
1529  ct => work%ctransformed(ispin)
1530  CALL cp_fm_to_fm(ct, cvec)
1531  CALL cp_fm_row_scale(cvec, tv)
1532  CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1533  cvec, 0.0_dp, ucmatrix)
1534  CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1535  x(ispin), 0.0_dp, uxmatrix)
1536  CALL parallel_gemm('N', 'T', nsgf, nsgf, norb, 1.0_dp, uxmatrix, ucmatrix, 0.0_dp, t0matrix)
1537  CALL cp_fm_to_fm(xtransformed(ispin), cvec)
1538  CALL cp_fm_row_scale(cvec, tv)
1539  CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1540  cvec, 0.0_dp, uxmatrix)
1541  CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1542  gs_mos(ispin)%mos_occ, 0.0_dp, ucmatrix)
1543  CALL parallel_gemm('N', 'T', nsgf, nsgf, norb, 1.0_dp, ucmatrix, uxmatrix, 1.0_dp, t0matrix)
1544  CALL cp_fm_schur_product(work%slambda, t0matrix, t1matrix)
1545  !
1546  CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, spinfac, work%S_eigenvectors, t1matrix, &
1547  0.0_dp, t0matrix)
1548  CALL cp_dbcsr_plus_fm_fm_t(matrix_plo(ispin)%matrix, matrix_v=t0matrix, &
1549  matrix_g=work%S_eigenvectors, ncol=nsgf, alpha=2.0_dp, symmetry_mode=1)
1550  CALL fm_pool_give_back_fm(work%fm_pool_ao_mo_occ(ispin)%pool, ucmatrix)
1551  DEALLOCATE (ucmatrix)
1552  CALL fm_pool_give_back_fm(work%fm_pool_ao_mo_occ(ispin)%pool, uxmatrix)
1553  DEALLOCATE (uxmatrix)
1554  CALL cp_fm_release(t0matrix)
1555  CALL cp_fm_release(t1matrix)
1556  !
1557  ! CV(mu,i) = TV(mu)*XT(mu,i)
1558  CALL cp_fm_to_fm(xtransformed(ispin), cvec)
1559  CALL cp_fm_row_scale(cvec, tv)
1560  CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, cpmos(ispin), norb, 2.0_dp*spinfac, 1.0_dp)
1561  ! CV(mu,i) = TV(mu)*CT(mu,i)
1562  ct => work%ctransformed(ispin)
1563  CALL cp_fm_to_fm(ct, cvec)
1564  CALL cp_fm_row_scale(cvec, tv)
1565  ! Shalf(nu,mu)*CV(mu,i)
1566  CALL cp_fm_get_info(cvec, matrix_struct=fmstruct, nrow_global=nao)
1567  CALL cp_fm_create(vcvec, fmstruct)
1568  CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, vcvec, norb, 1.0_dp, 0.0_dp)
1569  CALL cp_fm_struct_create(fmstruct_mat, context=fmstruct%context, nrow_global=norb, &
1570  ncol_global=norb, para_env=fmstruct%para_env)
1571  CALL cp_fm_create(cvcmat, fmstruct_mat)
1572  CALL cp_fm_struct_release(fmstruct_mat)
1573  CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, gs_mos(ispin)%mos_occ, vcvec, 0.0_dp, cvcmat)
1574  CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, x(ispin), cvcmat, 0.0_dp, vcvec)
1575  CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
1576  nactive(ispin), alpha=-2.0_dp*spinfac, beta=1.0_dp)
1577  ! wx1
1578  alpha = 2.0_dp
1579  CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=gs_mos(ispin)%mos_occ, &
1580  matrix_g=vcvec, ncol=norb, alpha=2.0_dp*alpha, symmetry_mode=1)
1581  CALL cp_fm_release(vcvec)
1582  CALL cp_fm_release(cvcmat)
1583  END IF
1584  !
1585  ! *** Exchange contribution
1586  !
1587  IF (stda_env%do_exchange) THEN
1588  !
1589  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1590  !
1591  norb = nactive(ispin)
1592  !
1593  tempmat => work%shalf
1594  CALL dbcsr_create(pdens, template=tempmat, matrix_type=dbcsr_type_no_symmetry)
1595  ! P(nu,mu) = SUM_j XT(nu,j)*CT(mu,j)
1596  ct => work%ctransformed(ispin)
1597  CALL dbcsr_set(pdens, 0.0_dp)
1598  CALL cp_dbcsr_plus_fm_fm_t(pdens, xtransformed(ispin), ct, nactive(ispin), &
1599  1.0_dp, keep_sparsity=.false.)
1600  CALL dbcsr_filter(pdens, stda_env%eps_td_filter)
1601  ! Apply PP*gab -> PP; gab = gamma_coulomb
1602  ! P(nu,mu) = P(nu,mu)*g(a of nu,b of mu)
1603  bp = stda_env%beta_param
1604  hfx = stda_env%hfx_fraction
1605  CALL dbcsr_iterator_start(iter, pdens)
1606  DO WHILE (dbcsr_iterator_blocks_left(iter))
1607  CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock, blk)
1608  rij = particle_set(iatom)%r - particle_set(jatom)%r
1609  rij = pbc(rij, cell)
1610  dr = sqrt(sum(rij(:)**2))
1611  ikind = kind_of(iatom)
1612  jkind = kind_of(jatom)
1613  eta = (stda_env%kind_param_set(ikind)%kind_param%hardness_param + &
1614  stda_env%kind_param_set(jkind)%kind_param%hardness_param)/2.0_dp
1615  rbeta = dr**bp
1616  IF (hfx > 0.0_dp) THEN
1617  gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
1618  ELSE
1619  IF (dr < 1.0e-6_dp) THEN
1620  gabr = 0.0_dp
1621  ELSE
1622  gabr = 1._dp/dr
1623  END IF
1624  END IF
1625  ! gabr = (1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp)
1626  ! forces
1627  IF (dr > 1.0e-6_dp) THEN
1628  IF (hfx > 0.0_dp) THEN
1629  dgabr = -(1._dp/bp)*(1._dp/(rbeta + (hfx*eta)**(-bp)))**(1._dp/bp + 1._dp)
1630  dgabr = bp*rbeta/dr**2*dgabr
1631  dgabr = sum(pblock**2)*dgabr
1632  ELSE
1633  dgabr = -1.0_dp/dr**3
1634  dgabr = sum(pblock**2)*dgabr
1635  END IF
1636  atom_i = atom_of_kind(iatom)
1637  atom_j = atom_of_kind(jatom)
1638  DO i = 1, 3
1639  fij(i) = dgabr*rij(i)
1640  END DO
1641  force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
1642  force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
1643  force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
1644  force(jkind)%rho_elec(1, atom_j) = force(jkind)%rho_elec(1, atom_j) + fij(1)
1645  force(jkind)%rho_elec(2, atom_j) = force(jkind)%rho_elec(2, atom_j) + fij(2)
1646  force(jkind)%rho_elec(3, atom_j) = force(jkind)%rho_elec(3, atom_j) + fij(3)
1647  END IF
1648  !
1649  pblock = gabr*pblock
1650  END DO
1651  CALL dbcsr_iterator_stop(iter)
1652  !
1653  ! Transpose pdens matrix
1654  CALL dbcsr_create(ptrans, template=pdens)
1655  CALL dbcsr_transposed(ptrans, pdens)
1656  !
1657  ! forces from Lowdin charge derivative
1658  CALL cp_fm_get_info(work%S_C0_C0T(ispin), matrix_struct=fmstruct)
1659  CALL cp_fm_create(t0matrix, matrix_struct=fmstruct, name="T0 SCRATCH")
1660  CALL cp_fm_create(t1matrix, matrix_struct=fmstruct, name="T1 SCRATCH")
1661  ALLOCATE (ucmatrix)
1662  CALL fm_pool_create_fm(work%fm_pool_ao_mo_occ(ispin)%pool, ucmatrix)
1663  ALLOCATE (uxmatrix)
1664  CALL fm_pool_create_fm(work%fm_pool_ao_mo_occ(ispin)%pool, uxmatrix)
1665  ct => work%ctransformed(ispin)
1666  CALL cp_dbcsr_sm_fm_multiply(pdens, ct, cvec, norb, 1.0_dp, 0.0_dp)
1667  CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1668  cvec, 0.0_dp, ucmatrix)
1669  CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1670  x(ispin), 0.0_dp, uxmatrix)
1671  CALL parallel_gemm('N', 'T', nsgf, nsgf, norb, 1.0_dp, uxmatrix, ucmatrix, 0.0_dp, t0matrix)
1672  CALL cp_dbcsr_sm_fm_multiply(ptrans, xtransformed(ispin), cvec, norb, 1.0_dp, 0.0_dp)
1673  CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1674  cvec, 0.0_dp, uxmatrix)
1675  CALL parallel_gemm('T', 'N', nsgf, norb, nsgf, 1.0_dp, work%S_eigenvectors, &
1676  gs_mos(ispin)%mos_occ, 0.0_dp, ucmatrix)
1677  CALL parallel_gemm('N', 'T', nsgf, nsgf, norb, 1.0_dp, ucmatrix, uxmatrix, 1.0_dp, t0matrix)
1678  CALL cp_fm_schur_product(work%slambda, t0matrix, t1matrix)
1679  CALL parallel_gemm('N', 'N', nsgf, nsgf, nsgf, -1.0_dp, work%S_eigenvectors, t1matrix, &
1680  0.0_dp, t0matrix)
1681  CALL cp_dbcsr_plus_fm_fm_t(matrix_plo(ispin)%matrix, matrix_v=t0matrix, &
1682  matrix_g=work%S_eigenvectors, ncol=nsgf, alpha=2.0_dp, symmetry_mode=1)
1683  CALL fm_pool_give_back_fm(work%fm_pool_ao_mo_occ(ispin)%pool, ucmatrix)
1684  DEALLOCATE (ucmatrix)
1685  CALL fm_pool_give_back_fm(work%fm_pool_ao_mo_occ(ispin)%pool, uxmatrix)
1686  DEALLOCATE (uxmatrix)
1687  CALL cp_fm_release(t0matrix)
1688  CALL cp_fm_release(t1matrix)
1689 
1690  ! RHS contribution to response matrix
1691  ! CV(nu,i) = P(nu,mu)*XT(mu,i)
1692  CALL cp_dbcsr_sm_fm_multiply(ptrans, xtransformed(ispin), cvec, norb, 1.0_dp, 0.0_dp)
1693  CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, cpmos(ispin), norb, &
1694  alpha=-xfac, beta=1.0_dp)
1695  !
1696  CALL cp_fm_get_info(cvec, matrix_struct=fmstruct, nrow_global=nao)
1697  CALL cp_fm_create(vcvec, fmstruct)
1698  ! CV(nu,i) = P(nu,mu)*CT(mu,i)
1699  CALL cp_dbcsr_sm_fm_multiply(ptrans, ct, cvec, norb, 1.0_dp, 0.0_dp)
1700  CALL cp_dbcsr_sm_fm_multiply(work%shalf, cvec, vcvec, norb, 1.0_dp, 0.0_dp)
1701  CALL cp_fm_struct_create(fmstruct_mat, context=fmstruct%context, nrow_global=norb, &
1702  ncol_global=norb, para_env=fmstruct%para_env)
1703  CALL cp_fm_create(cvcmat, fmstruct_mat)
1704  CALL cp_fm_struct_release(fmstruct_mat)
1705  CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, gs_mos(ispin)%mos_occ, vcvec, 0.0_dp, cvcmat)
1706  CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, x(ispin), cvcmat, 0.0_dp, vcvec)
1707  CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, vcvec, cpmos(ispin), &
1708  norb, alpha=xfac, beta=1.0_dp)
1709  ! wx1
1710  IF (nspins == 2) THEN
1711  alpha = -2.0_dp
1712  ELSE
1713  alpha = -1.0_dp
1714  END IF
1715  CALL cp_dbcsr_plus_fm_fm_t(matrix_wx1(ispin)%matrix, matrix_v=gs_mos(ispin)%mos_occ, &
1716  matrix_g=vcvec, &
1717  ncol=norb, alpha=2.0_dp*alpha, symmetry_mode=1)
1718  CALL cp_fm_release(vcvec)
1719  CALL cp_fm_release(cvcmat)
1720  !
1721  CALL dbcsr_release(pdens)
1722  CALL dbcsr_release(ptrans)
1723  !
1724  IF (debug_forces) THEN
1725  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1726  CALL para_env%sum(fodeb)
1727  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Exch[X] ", fodeb
1728  END IF
1729  END IF
1730  !
1731  CALL cp_fm_release(cvec)
1732  CALL cp_fm_release(xvec)
1733  DEALLOCATE (tv)
1734  END DO
1735 
1736  CALL cp_fm_release(xtransformed)
1737  DEALLOCATE (tcharge, gtcharge)
1738  DEALLOCATE (first_sgf, last_sgf)
1739 
1740  ! Lowdin forces
1741  IF (nspins == 2) THEN
1742  CALL dbcsr_add(matrix_plo(1)%matrix, matrix_plo(2)%matrix, &
1743  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
1744  END IF
1745  CALL dbcsr_scale(matrix_plo(1)%matrix, -1.0_dp)
1746  NULLIFY (scrm)
1747  IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
1748  CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
1749  matrix_name="OVERLAP MATRIX", &
1750  basis_type_a="ORB", basis_type_b="ORB", &
1751  sab_nl=sab_orb, calculate_forces=.true., &
1752  matrix_p=matrix_plo(1)%matrix)
1753  CALL dbcsr_deallocate_matrix_set(scrm)
1754  CALL dbcsr_deallocate_matrix_set(matrix_plo)
1755  IF (debug_forces) THEN
1756  fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
1757  CALL para_env%sum(fodeb)
1758  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Lowdin ", fodeb
1759  END IF
1760 
1761  IF (ASSOCIATED(ex_env%matrix_wx1)) CALL dbcsr_deallocate_matrix_set(ex_env%matrix_wx1)
1762  ex_env%matrix_wx1 => matrix_wx1
1763 
1764  CALL timestop(handle)
1765 
1766  END SUBROUTINE stda_force
1767 
1768 ! **************************************************************************************************
1769 
1770 END MODULE qs_tddfpt2_fhxc_forces
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
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_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:744
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:208
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:887
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:118
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:1945
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:728
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:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
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:1016
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:1252
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:167
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:123
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, atprop)
...
Definition: ewald_methods_tb.F:364
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, atprop)
...
Definition: ewald_methods_tb.F:90
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:1033
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:3036
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::xc_kernel_method_best
integer, parameter, public xc_kernel_method_best
Definition: input_constants.F:699
input_constants::xc_kernel_method_analytic
integer, parameter, public xc_kernel_method_analytic
Definition: input_constants.F:699
input_constants::do_admm_aux_exch_func_none
integer, parameter, public do_admm_aux_exch_func_none
Definition: input_constants.F:863
input_constants::tddfpt_kernel_full
integer, parameter, public tddfpt_kernel_full
Definition: input_constants.F:581
input_constants::xc_kernel_method_numeric
integer, parameter, public xc_kernel_method_numeric
Definition: input_constants.F:699
input_constants::xc_none
integer, parameter, public xc_none
Definition: input_constants.F:551
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:724
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:697
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:1040
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:1710
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_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, 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, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
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, 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, rhs)
Set the QUICKSTEP environment.
Definition: qs_environment_types.F:1035
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)
...
Definition: qs_gapw_densities.F:49
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::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)
...
Definition: qs_rho0_ggrid.F:316
qs_rho0_ggrid::rho0_s_grid_create
subroutine, public rho0_s_grid_create(pw_env, rho0_mpole)
...
Definition: qs_rho0_ggrid.F:260
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:99
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:755
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:1155
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