(git:34ef472)
response_solver.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 ! **************************************************************************************************
9 !> \brief Calculate the CPKS equation and the resulting forces
10 !> \par History
11 !> 03.2014 created
12 !> 09.2019 Moved from KG to Kohn-Sham
13 !> 11.2019 Moved from energy_correction
14 !> 08.2020 AO linear response solver [fbelle]
15 !> \author JGH
16 ! **************************************************************************************************
17 MODULE response_solver
18  USE admm_methods, ONLY: admm_projection_derivative
19  USE admm_types, ONLY: admm_type,&
20  get_admm_env
21  USE atomic_kind_types, ONLY: atomic_kind_type,&
22  get_atomic_kind
23  USE cell_types, ONLY: cell_type
24  USE core_ae, ONLY: build_core_ae
25  USE core_ppl, ONLY: build_core_ppl
26  USE core_ppnl, ONLY: build_core_ppnl
27  USE cp_blacs_env, ONLY: cp_blacs_env_type
28  USE cp_control_types, ONLY: dft_control_type
29  USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
30  USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
31  copy_fm_to_dbcsr,&
32  cp_dbcsr_sm_fm_multiply,&
33  dbcsr_allocate_matrix_set,&
34  dbcsr_deallocate_matrix_set
35  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
36  cp_fm_struct_release,&
37  cp_fm_struct_type
38  USE cp_fm_types, ONLY: cp_fm_create,&
39  cp_fm_init_random,&
40  cp_fm_release,&
41  cp_fm_set_all,&
42  cp_fm_to_fm,&
43  cp_fm_type
44  USE cp_log_handling, ONLY: cp_get_default_logger,&
45  cp_logger_get_default_unit_nr,&
46  cp_logger_type
47  USE dbcsr_api, ONLY: &
48  dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_multiply, &
49  dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
50  USE ec_env_types, ONLY: energy_correction_type
51  USE ec_methods, ONLY: create_kernel,&
52  ec_mos_init
53  USE ec_orth_solver, ONLY: ec_response_ao
54  USE exstates_types, ONLY: excited_energy_type
55  USE hartree_local_methods, ONLY: vh_1c_gg_integrals,&
56  init_coulomb_local
57  USE hartree_local_types, ONLY: hartree_local_create,&
58  hartree_local_release,&
59  hartree_local_type
60  USE hfx_derivatives, ONLY: derivatives_four_center
61  USE hfx_energy_potential, ONLY: integrate_four_center
62  USE hfx_ri, ONLY: hfx_ri_update_forces,&
63  hfx_ri_update_ks
64  USE hfx_types, ONLY: hfx_type
65  USE input_constants, ONLY: &
66  do_admm_aux_exch_func_none, ec_ls_solver, ec_mo_solver, kg_tnadd_atomic, kg_tnadd_embed, &
67  kg_tnadd_embed_ri, ls_s_sqrt_ns, ls_s_sqrt_proot, ot_precond_full_all, &
68  ot_precond_full_kinetic, ot_precond_full_single, ot_precond_full_single_inverse, &
69  ot_precond_none, ot_precond_s_inverse, precond_mlp, 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 kg_correction, ONLY: kg_ekin_subset
75  USE kg_environment_types, ONLY: kg_environment_type
76  USE kg_tnadd_mat, ONLY: build_tnadd_mat
77  USE kinds, ONLY: default_string_length,&
78  dp
79  USE machine, ONLY: m_flush
80  USE mathlib, ONLY: det_3x3
81  USE message_passing, ONLY: mp_para_env_type
82  USE mulliken, ONLY: ao_charges
83  USE parallel_gemm_api, ONLY: parallel_gemm
84  USE particle_types, ONLY: particle_type
85  USE physcon, ONLY: pascal
86  USE pw_env_types, ONLY: pw_env_get,&
87  pw_env_type
88  USE pw_methods, ONLY: pw_axpy,&
89  pw_copy,&
90  pw_integral_ab,&
91  pw_scale,&
92  pw_transfer,&
93  pw_zero
94  USE pw_poisson_methods, ONLY: pw_poisson_solve
95  USE pw_poisson_types, ONLY: pw_poisson_type
96  USE pw_pool_types, ONLY: pw_pool_type
97  USE pw_types, ONLY: pw_c1d_gs_type,&
98  pw_r3d_rs_type
99  USE qs_2nd_kernel_ao, ONLY: build_dm_response
100  USE qs_collocate_density, ONLY: calculate_rho_elec
101  USE qs_density_matrices, ONLY: calculate_whz_matrix,&
102  calculate_wz_matrix
103  USE qs_energy_types, ONLY: qs_energy_type
104  USE qs_environment_types, ONLY: get_qs_env,&
105  qs_environment_type,&
106  set_qs_env
107  USE qs_force_types, ONLY: qs_force_type,&
108  total_qs_force
109  USE qs_gapw_densities, ONLY: prepare_gapw_den
110  USE qs_integrate_potential, ONLY: integrate_v_core_rspace,&
111  integrate_v_rspace
112  USE qs_kind_types, ONLY: get_qs_kind,&
113  get_qs_kind_set,&
114  qs_kind_type
115  USE qs_kinetic, ONLY: build_kinetic_matrix
116  USE qs_ks_atom, ONLY: update_ks_atom
117  USE qs_ks_methods, ONLY: calc_rho_tot_gspace
118  USE qs_ks_types, ONLY: qs_ks_env_type
119  USE qs_linres_methods, ONLY: linres_solver
120  USE qs_linres_types, ONLY: linres_control_type
121  USE qs_local_rho_types, ONLY: local_rho_set_create,&
122  local_rho_set_release,&
123  local_rho_type
124  USE qs_matrix_pools, ONLY: mpools_rebuild_fm_pools
125  USE qs_mo_methods, ONLY: make_basis_sm
126  USE qs_mo_types, ONLY: deallocate_mo_set,&
127  get_mo_set,&
128  mo_set_type
129  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
130  USE qs_oce_types, ONLY: oce_matrix_type
131  USE qs_overlap, ONLY: build_overlap_matrix
132  USE qs_p_env_methods, ONLY: p_env_create,&
133  p_env_psi0_changed
134  USE qs_p_env_types, ONLY: p_env_release,&
135  qs_p_env_type
136  USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace,&
137  rho0_s_grid_create
138  USE qs_rho0_methods, ONLY: init_rho0
139  USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals,&
140  calculate_rho_atom_coeff
141  USE qs_rho_types, ONLY: qs_rho_get,&
142  qs_rho_type
143  USE qs_vxc_atom, ONLY: calculate_vxc_atom,&
144  calculate_xc_2nd_deriv_atom
145  USE task_list_types, ONLY: task_list_type
146  USE virial_methods, ONLY: one_third_sum_diag
147  USE virial_types, ONLY: virial_type
148  USE xtb_ehess, ONLY: xtb_coulomb_hessian
149  USE xtb_ehess_force, ONLY: calc_xtb_ehess_force
150  USE xtb_matrices, ONLY: xtb_hab_force
151  USE xtb_types, ONLY: get_xtb_atom_param,&
152  xtb_atom_type
153 #include "./base/base_uses.f90"
154 
155  IMPLICIT NONE
156 
157  PRIVATE
158 
159  ! Global parameters
160 
161  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'response_solver'
162 
163  PUBLIC :: response_calculation, response_equation, response_force, response_force_xtb, &
164  response_equation_new
165 
166 ! **************************************************************************************************
167 
168 CONTAINS
169 
170 ! **************************************************************************************************
171 !> \brief Initializes solver of linear response equation for energy correction
172 !> \brief Call AO or MO based linear response solver for energy correction
173 !>
174 !> \param qs_env The quickstep environment
175 !> \param ec_env The energy correction environment
176 !>
177 !> \date 01.2020
178 !> \author Fabian Belleflamme
179 ! **************************************************************************************************
180  SUBROUTINE response_calculation(qs_env, ec_env)
181  TYPE(qs_environment_type), POINTER :: qs_env
182  TYPE(energy_correction_type), POINTER :: ec_env
183 
184  CHARACTER(LEN=*), PARAMETER :: routinen = 'response_calculation'
185 
186  INTEGER :: handle, homo, ispin, nao, nao_aux, nmo, &
187  nocc, nspins, solver_method, unit_nr
188  LOGICAL :: should_stop
189  REAL(kind=dp) :: focc
190  TYPE(admm_type), POINTER :: admm_env
191  TYPE(cp_blacs_env_type), POINTER :: blacs_env
192  TYPE(cp_fm_struct_type), POINTER :: fm_struct
193  TYPE(cp_fm_type) :: sv
194  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: cpmos, mo_occ
195  TYPE(cp_fm_type), POINTER :: mo_coeff
196  TYPE(cp_logger_type), POINTER :: logger
197  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_s_aux, rho_ao
198  TYPE(dft_control_type), POINTER :: dft_control
199  TYPE(linres_control_type), POINTER :: linres_control
200  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
201  TYPE(mp_para_env_type), POINTER :: para_env
202  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
203  POINTER :: sab_orb
204  TYPE(qs_energy_type), POINTER :: energy
205  TYPE(qs_p_env_type), POINTER :: p_env
206  TYPE(qs_rho_type), POINTER :: rho
207  TYPE(section_vals_type), POINTER :: input, solver_section
208 
209  CALL timeset(routinen, handle)
210 
211  NULLIFY (admm_env, dft_control, energy, logger, matrix_s, matrix_s_aux, mo_coeff, mos, para_env, &
212  rho_ao, sab_orb, solver_section)
213 
214  ! Get useful output unit
215  logger => cp_get_default_logger()
216  IF (logger%para_env%is_source()) THEN
217  unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
218  ELSE
219  unit_nr = -1
220  END IF
221 
222  CALL get_qs_env(qs_env, &
223  dft_control=dft_control, &
224  input=input, &
225  matrix_s=matrix_s, &
226  para_env=para_env, &
227  sab_orb=sab_orb)
228  nspins = dft_control%nspins
229 
230  ! initialize linres_control
231  NULLIFY (linres_control)
232  ALLOCATE (linres_control)
233  linres_control%do_kernel = .true.
234  linres_control%lr_triplet = .false.
235  linres_control%converged = .false.
236  linres_control%energy_gap = 0.02_dp
237 
238  ! Read input
239  solver_section => section_vals_get_subs_vals(input, "DFT%ENERGY_CORRECTION%RESPONSE_SOLVER")
240  CALL section_vals_val_get(solver_section, "EPS", r_val=linres_control%eps)
241  CALL section_vals_val_get(solver_section, "EPS_FILTER", r_val=linres_control%eps_filter)
242  CALL section_vals_val_get(solver_section, "MAX_ITER", i_val=linres_control%max_iter)
243  CALL section_vals_val_get(solver_section, "METHOD", i_val=solver_method)
244  CALL section_vals_val_get(solver_section, "PRECONDITIONER", i_val=linres_control%preconditioner_type)
245  CALL section_vals_val_get(solver_section, "RESTART", l_val=linres_control%linres_restart)
246  CALL section_vals_val_get(solver_section, "RESTART_EVERY", i_val=linres_control%restart_every)
247  CALL set_qs_env(qs_env, linres_control=linres_control)
248 
249  ! Write input section of response solver
250  CALL response_solver_write_input(solver_section, linres_control, unit_nr)
251 
252  ! Allocate and initialize response density matrix Z,
253  ! and the energy weighted response density matrix
254  ! Template is the ground-state overlap matrix
255  CALL dbcsr_allocate_matrix_set(ec_env%matrix_wz, nspins)
256  CALL dbcsr_allocate_matrix_set(ec_env%matrix_z, nspins)
257  DO ispin = 1, nspins
258  ALLOCATE (ec_env%matrix_wz(ispin)%matrix)
259  ALLOCATE (ec_env%matrix_z(ispin)%matrix)
260  CALL dbcsr_create(ec_env%matrix_wz(ispin)%matrix, name="Wz MATRIX", &
261  template=matrix_s(1)%matrix)
262  CALL dbcsr_create(ec_env%matrix_z(ispin)%matrix, name="Z MATRIX", &
263  template=matrix_s(1)%matrix)
264  CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_wz(ispin)%matrix, sab_orb)
265  CALL cp_dbcsr_alloc_block_from_nbl(ec_env%matrix_z(ispin)%matrix, sab_orb)
266  CALL dbcsr_set(ec_env%matrix_wz(ispin)%matrix, 0.0_dp)
267  CALL dbcsr_set(ec_env%matrix_z(ispin)%matrix, 0.0_dp)
268  END DO
269 
270  ! MO solver requires MO's of the ground-state calculation,
271  ! The MOs environment is not allocated if LS-DFT has been used.
272  ! Introduce MOs here
273  ! Remark: MOS environment also required for creation of p_env
274  IF (dft_control%qs_control%do_ls_scf) THEN
275 
276  ! Allocate and initialize MO environment
277  CALL ec_mos_init(qs_env, matrix_s(1)%matrix)
278  CALL get_qs_env(qs_env, mos=mos, rho=rho)
279 
280  ! Get ground-state density matrix
281  CALL qs_rho_get(rho, rho_ao=rho_ao)
282 
283  DO ispin = 1, nspins
284  CALL get_mo_set(mo_set=mos(ispin), &
285  mo_coeff=mo_coeff, &
286  nmo=nmo, nao=nao, homo=homo)
287 
288  CALL cp_fm_set_all(mo_coeff, 0.0_dp)
289  CALL cp_fm_init_random(mo_coeff, nmo)
290 
291  CALL cp_fm_create(sv, mo_coeff%matrix_struct, "SV")
292  ! multiply times PS
293  ! PS*C(:,1:nomo)+C(:,nomo+1:nmo) (nomo=NINT(nelectron/maxocc))
294  CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, mo_coeff, sv, nmo)
295  CALL cp_dbcsr_sm_fm_multiply(rho_ao(ispin)%matrix, sv, mo_coeff, homo)
296  CALL cp_fm_release(sv)
297  ! and ortho the result
298  CALL make_basis_sm(mo_coeff, nmo, matrix_s(1)%matrix)
299 
300  ! rebuilds fm_pools
301  ! originally done in qs_env_setup, only when mos associated
302  NULLIFY (blacs_env)
303  CALL get_qs_env(qs_env, blacs_env=blacs_env)
304  CALL mpools_rebuild_fm_pools(qs_env%mpools, mos=mos, &
305  blacs_env=blacs_env, para_env=para_env)
306  END DO
307  END IF
308 
309  ! initialize p_env
310  ! Remark: mos environment is needed for this
311  IF (ASSOCIATED(ec_env%p_env)) THEN
312  CALL p_env_release(ec_env%p_env)
313  DEALLOCATE (ec_env%p_env)
314  NULLIFY (ec_env%p_env)
315  END IF
316  ALLOCATE (ec_env%p_env)
317  CALL p_env_create(ec_env%p_env, qs_env, orthogonal_orbitals=.true., &
318  linres_control=linres_control)
319  CALL p_env_psi0_changed(ec_env%p_env, qs_env)
320  ! Total energy overwritten, replace with Etot from energy correction
321  CALL get_qs_env(qs_env, energy=energy)
322  energy%total = ec_env%etotal
323  !
324  p_env => ec_env%p_env
325  !
326  CALL dbcsr_allocate_matrix_set(p_env%p1, nspins)
327  CALL dbcsr_allocate_matrix_set(p_env%w1, nspins)
328  DO ispin = 1, nspins
329  ALLOCATE (p_env%p1(ispin)%matrix, p_env%w1(ispin)%matrix)
330  CALL dbcsr_create(matrix=p_env%p1(ispin)%matrix, template=matrix_s(1)%matrix)
331  CALL dbcsr_create(matrix=p_env%w1(ispin)%matrix, template=matrix_s(1)%matrix)
332  CALL cp_dbcsr_alloc_block_from_nbl(p_env%p1(ispin)%matrix, sab_orb)
333  CALL cp_dbcsr_alloc_block_from_nbl(p_env%w1(ispin)%matrix, sab_orb)
334  END DO
335  IF (dft_control%do_admm) THEN
336  CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux)
337  CALL dbcsr_allocate_matrix_set(p_env%p1_admm, nspins)
338  DO ispin = 1, nspins
339  ALLOCATE (p_env%p1_admm(ispin)%matrix)
340  CALL dbcsr_create(p_env%p1_admm(ispin)%matrix, &
341  template=matrix_s_aux(1)%matrix)
342  CALL dbcsr_copy(p_env%p1_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
343  CALL dbcsr_set(p_env%p1_admm(ispin)%matrix, 0.0_dp)
344  END DO
345  END IF
346 
347  ! Choose between MO-solver and AO-solver
348  SELECT CASE (solver_method)
349  CASE (ec_mo_solver)
350 
351  ! CPKS vector cpmos - RHS of response equation as Ax + b = 0 (sign of b)
352  ! Sign is changed in linres_solver!
353  ! Projector Q applied in linres_solver!
354  CALL get_qs_env(qs_env, mos=mos)
355  ALLOCATE (cpmos(nspins), mo_occ(nspins))
356  DO ispin = 1, nspins
357  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
358  NULLIFY (fm_struct)
359  CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, &
360  template_fmstruct=mo_coeff%matrix_struct)
361  CALL cp_fm_create(cpmos(ispin), fm_struct)
362  CALL cp_fm_set_all(cpmos(ispin), 0.0_dp)
363  CALL cp_fm_create(mo_occ(ispin), fm_struct)
364  CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc)
365  CALL cp_fm_struct_release(fm_struct)
366  END DO
367 
368  focc = 2.0_dp
369  IF (nspins == 1) focc = 4.0_dp
370  DO ispin = 1, nspins
371  CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
372  CALL cp_dbcsr_sm_fm_multiply(ec_env%matrix_hz(ispin)%matrix, mo_occ(ispin), &
373  cpmos(ispin), nocc, &
374  alpha=focc, beta=0.0_dp)
375  END DO
376  CALL cp_fm_release(mo_occ)
377 
378  CALL response_equation_new(qs_env, p_env, cpmos, unit_nr)
379 
380  ! Get the response density matrix,
381  ! and energy-weighted response density matrix
382  DO ispin = 1, nspins
383  CALL dbcsr_copy(ec_env%matrix_z(ispin)%matrix, p_env%p1(ispin)%matrix)
384  CALL dbcsr_copy(ec_env%matrix_wz(ispin)%matrix, p_env%w1(ispin)%matrix)
385  END DO
386  CALL cp_fm_release(cpmos)
387 
388  CASE (ec_ls_solver)
389 
390  ! AO ortho solver
391  CALL ec_response_ao(qs_env=qs_env, &
392  p_env=p_env, &
393  matrix_hz=ec_env%matrix_hz, &
394  matrix_pz=ec_env%matrix_z, &
395  matrix_wz=ec_env%matrix_wz, &
396  iounit=unit_nr, &
397  should_stop=should_stop)
398 
399  IF (dft_control%do_admm) THEN
400  CALL get_qs_env(qs_env, admm_env=admm_env)
401  cpassert(ASSOCIATED(admm_env%work_orb_orb))
402  cpassert(ASSOCIATED(admm_env%work_aux_orb))
403  cpassert(ASSOCIATED(admm_env%work_aux_aux))
404  nao = admm_env%nao_orb
405  nao_aux = admm_env%nao_aux_fit
406  DO ispin = 1, nspins
407  CALL copy_dbcsr_to_fm(ec_env%matrix_z(ispin)%matrix, admm_env%work_orb_orb)
408  CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
409  1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
410  admm_env%work_aux_orb)
411  CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
412  1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
413  admm_env%work_aux_aux)
414  CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, p_env%p1_admm(ispin)%matrix, &
415  keep_sparsity=.true.)
416  END DO
417  END IF
418 
419  CASE DEFAULT
420  cpabort("Unknown solver for response equation requested")
421  END SELECT
422 
423  IF (dft_control%do_admm) THEN
424  CALL dbcsr_allocate_matrix_set(ec_env%z_admm, nspins)
425  DO ispin = 1, nspins
426  ALLOCATE (ec_env%z_admm(ispin)%matrix)
427  CALL dbcsr_create(matrix=ec_env%z_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
428  CALL get_qs_env(qs_env, admm_env=admm_env)
429  CALL dbcsr_copy(ec_env%z_admm(ispin)%matrix, p_env%p1_admm(ispin)%matrix)
430  END DO
431  END IF
432 
433  ! Get rid of MO environment again
434  IF (dft_control%qs_control%do_ls_scf) THEN
435  DO ispin = 1, nspins
436  CALL deallocate_mo_set(mos(ispin))
437  END DO
438  IF (ASSOCIATED(qs_env%mos)) THEN
439  DO ispin = 1, SIZE(qs_env%mos)
440  CALL deallocate_mo_set(qs_env%mos(ispin))
441  END DO
442  DEALLOCATE (qs_env%mos)
443  END IF
444  END IF
445 
446  CALL timestop(handle)
447 
448  END SUBROUTINE response_calculation
449 
450 ! **************************************************************************************************
451 !> \brief Parse the input section of the response solver
452 !> \param input Input section which controls response solver parameters
453 !> \param linres_control Environment for general setting of linear response calculation
454 !> \param unit_nr ...
455 !> \par History
456 !> 2020.05 created [Fabian Belleflamme]
457 !> \author Fabian Belleflamme
458 ! **************************************************************************************************
459  SUBROUTINE response_solver_write_input(input, linres_control, unit_nr)
460  TYPE(section_vals_type), POINTER :: input
461  TYPE(linres_control_type), POINTER :: linres_control
462  INTEGER, INTENT(IN) :: unit_nr
463 
464  CHARACTER(len=*), PARAMETER :: routinen = 'response_solver_write_input'
465 
466  INTEGER :: handle, max_iter_lanczos, s_sqrt_method, &
467  s_sqrt_order, solver_method
468  REAL(kind=dp) :: eps_lanczos
469 
470  CALL timeset(routinen, handle)
471 
472  IF (unit_nr > 0) THEN
473 
474  ! linres_control
475  WRITE (unit_nr, '(/,T2,A)') &
476  repeat("-", 30)//" Linear Response Solver "//repeat("-", 25)
477 
478  ! Which type of solver is used
479  CALL section_vals_val_get(input, "METHOD", i_val=solver_method)
480 
481  SELECT CASE (solver_method)
482  CASE (ec_ls_solver)
483  WRITE (unit_nr, '(T2,A,T61,A20)') "Solver: ", "AO-based CG-solver"
484  CASE (ec_mo_solver)
485  WRITE (unit_nr, '(T2,A,T61,A20)') "Solver: ", "MO-based CG-solver"
486  END SELECT
487 
488  WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps:", linres_control%eps
489  WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_filter:", linres_control%eps_filter
490  WRITE (unit_nr, '(T2,A,T61,I20)') "Max iter:", linres_control%max_iter
491 
492  SELECT CASE (linres_control%preconditioner_type)
493  CASE (ot_precond_full_all)
494  WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "FULL_ALL"
495  CASE (ot_precond_full_single_inverse)
496  WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "FULL_SINGLE_INVERSE"
497  CASE (ot_precond_full_single)
498  WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "FULL_SINGLE"
499  CASE (ot_precond_full_kinetic)
500  WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "FULL_KINETIC"
501  CASE (ot_precond_s_inverse)
502  WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "FULL_S_INVERSE"
503  CASE (precond_mlp)
504  WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "MULTI_LEVEL"
505  CASE (ot_precond_none)
506  WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "NONE"
507  END SELECT
508 
509  SELECT CASE (solver_method)
510  CASE (ec_ls_solver)
511 
512  CALL section_vals_val_get(input, "S_SQRT_METHOD", i_val=s_sqrt_method)
513  CALL section_vals_val_get(input, "S_SQRT_ORDER", i_val=s_sqrt_order)
514  CALL section_vals_val_get(input, "EPS_LANCZOS", r_val=eps_lanczos)
515  CALL section_vals_val_get(input, "MAX_ITER_LANCZOS", i_val=max_iter_lanczos)
516 
517  ! Response solver transforms P and KS into orthonormal basis,
518  ! reuires matrx S sqrt and its inverse
519  SELECT CASE (s_sqrt_method)
520  CASE (ls_s_sqrt_ns)
521  WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "NEWTONSCHULZ"
522  CASE (ls_s_sqrt_proot)
523  WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "PROOT"
524  CASE DEFAULT
525  cpabort("Unknown sqrt method.")
526  END SELECT
527  WRITE (unit_nr, '(T2,A,T61,I20)') "S sqrt order:", s_sqrt_order
528 
529  CASE (ec_mo_solver)
530  END SELECT
531 
532  WRITE (unit_nr, '(T2,A)') repeat("-", 79)
533 
534  CALL m_flush(unit_nr)
535  END IF
536 
537  CALL timestop(handle)
538 
539  END SUBROUTINE response_solver_write_input
540 
541 ! **************************************************************************************************
542 !> \brief Initializes vectors for MO-coefficient based linear response solver
543 !> and calculates response density, and energy-weighted response density matrix
544 !>
545 !> \param qs_env ...
546 !> \param p_env ...
547 !> \param cpmos ...
548 !> \param iounit ...
549 ! **************************************************************************************************
550  SUBROUTINE response_equation_new(qs_env, p_env, cpmos, iounit)
551  TYPE(qs_environment_type), POINTER :: qs_env
552  TYPE(qs_p_env_type) :: p_env
553  TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: cpmos
554  INTEGER, INTENT(IN) :: iounit
555 
556  CHARACTER(LEN=*), PARAMETER :: routinen = 'response_equation_new'
557 
558  INTEGER :: handle, ispin, nao, nao_aux, nocc, nspins
559  LOGICAL :: should_stop
560  TYPE(admm_type), POINTER :: admm_env
561  TYPE(cp_fm_struct_type), POINTER :: fm_struct
562  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: psi0, psi1
563  TYPE(cp_fm_type), POINTER :: mo_coeff
564  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
565  TYPE(dft_control_type), POINTER :: dft_control
566  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
567 
568  CALL timeset(routinen, handle)
569 
570  NULLIFY (dft_control, matrix_ks, mo_coeff, mos)
571 
572  CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks, &
573  matrix_s=matrix_s, mos=mos)
574  nspins = dft_control%nspins
575 
576  ! Initialize vectors:
577  ! psi0 : The ground-state MO-coefficients
578  ! psi1 : The "perturbed" linear response orbitals
579  ALLOCATE (psi0(nspins), psi1(nspins))
580  DO ispin = 1, nspins
581  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, homo=nocc)
582  NULLIFY (fm_struct)
583  CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, &
584  template_fmstruct=mo_coeff%matrix_struct)
585  CALL cp_fm_create(psi0(ispin), fm_struct)
586  CALL cp_fm_to_fm(mo_coeff, psi0(ispin), nocc)
587  CALL cp_fm_create(psi1(ispin), fm_struct)
588  CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
589  CALL cp_fm_struct_release(fm_struct)
590  END DO
591 
592  should_stop = .false.
593  ! The response solver
594  CALL linres_solver(p_env, qs_env, psi1, cpmos, psi0, iounit, should_stop)
595 
596  ! Building the response density matrix
597  DO ispin = 1, nspins
598  CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_s(1)%matrix)
599  END DO
600  CALL build_dm_response(psi0, psi1, p_env%p1)
601  DO ispin = 1, nspins
602  CALL dbcsr_scale(p_env%p1(ispin)%matrix, 0.5_dp)
603  END DO
604 
605  IF (dft_control%do_admm) THEN
606  CALL get_qs_env(qs_env, admm_env=admm_env)
607  cpassert(ASSOCIATED(admm_env%work_orb_orb))
608  cpassert(ASSOCIATED(admm_env%work_aux_orb))
609  cpassert(ASSOCIATED(admm_env%work_aux_aux))
610  nao = admm_env%nao_orb
611  nao_aux = admm_env%nao_aux_fit
612  DO ispin = 1, nspins
613  CALL copy_dbcsr_to_fm(p_env%p1(ispin)%matrix, admm_env%work_orb_orb)
614  CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
615  1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
616  admm_env%work_aux_orb)
617  CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
618  1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
619  admm_env%work_aux_aux)
620  CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, p_env%p1_admm(ispin)%matrix, &
621  keep_sparsity=.true.)
622  END DO
623  END IF
624 
625  ! Calculate Wz = 0.5*(psi1*eps*psi0^T + psi0*eps*psi1^T)
626  DO ispin = 1, nspins
627  CALL calculate_wz_matrix(mos(ispin), psi1(ispin), matrix_ks(ispin)%matrix, &
628  p_env%w1(ispin)%matrix)
629  END DO
630  DO ispin = 1, nspins
631  CALL cp_fm_release(cpmos(ispin))
632  END DO
633  CALL cp_fm_release(psi1)
634  CALL cp_fm_release(psi0)
635 
636  CALL timestop(handle)
637 
638  END SUBROUTINE response_equation_new
639 
640 ! **************************************************************************************************
641 !> \brief Initializes vectors for MO-coefficient based linear response solver
642 !> and calculates response density, and energy-weighted response density matrix
643 !>
644 !> \param qs_env ...
645 !> \param p_env ...
646 !> \param cpmos RHS of equation as Ax + b = 0 (sign of b)
647 !> \param iounit ...
648 !> \param lr_section ...
649 ! **************************************************************************************************
650  SUBROUTINE response_equation(qs_env, p_env, cpmos, iounit, lr_section)
651  TYPE(qs_environment_type), POINTER :: qs_env
652  TYPE(qs_p_env_type) :: p_env
653  TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos
654  INTEGER, INTENT(IN) :: iounit
655  TYPE(section_vals_type), OPTIONAL, POINTER :: lr_section
656 
657  CHARACTER(LEN=*), PARAMETER :: routinen = 'response_equation'
658 
659  INTEGER :: handle, ispin, nao, nao_aux, nocc, nspins
660  LOGICAL :: should_stop
661  TYPE(admm_type), POINTER :: admm_env
662  TYPE(cp_fm_struct_type), POINTER :: fm_struct
663  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: psi0, psi1
664  TYPE(cp_fm_type), POINTER :: mo_coeff
665  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_s_aux
666  TYPE(dft_control_type), POINTER :: dft_control
667  TYPE(linres_control_type), POINTER :: linres_control
668  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
669  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
670  POINTER :: sab_orb
671 
672  CALL timeset(routinen, handle)
673 
674  ! initialized linres_control
675  NULLIFY (linres_control)
676  ALLOCATE (linres_control)
677  linres_control%do_kernel = .true.
678  linres_control%lr_triplet = .false.
679  IF (PRESENT(lr_section)) THEN
680  CALL section_vals_val_get(lr_section, "RESTART", l_val=linres_control%linres_restart)
681  CALL section_vals_val_get(lr_section, "MAX_ITER", i_val=linres_control%max_iter)
682  CALL section_vals_val_get(lr_section, "EPS", r_val=linres_control%eps)
683  CALL section_vals_val_get(lr_section, "EPS_FILTER", r_val=linres_control%eps_filter)
684  CALL section_vals_val_get(lr_section, "RESTART_EVERY", i_val=linres_control%restart_every)
685  CALL section_vals_val_get(lr_section, "PRECONDITIONER", i_val=linres_control%preconditioner_type)
686  CALL section_vals_val_get(lr_section, "ENERGY_GAP", r_val=linres_control%energy_gap)
687  ELSE
688  linres_control%linres_restart = .false.
689  linres_control%max_iter = 100
690  linres_control%eps = 1.0e-10_dp
691  linres_control%eps_filter = 1.0e-15_dp
692  linres_control%restart_every = 50
693  linres_control%preconditioner_type = ot_precond_full_single_inverse
694  linres_control%energy_gap = 0.02_dp
695  END IF
696 
697  ! initialized p_env
698  CALL p_env_create(p_env, qs_env, orthogonal_orbitals=.true., &
699  linres_control=linres_control)
700  CALL set_qs_env(qs_env, linres_control=linres_control)
701  CALL p_env_psi0_changed(p_env, qs_env)
702  p_env%new_preconditioner = .true.
703 
704  CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos)
705  !
706  nspins = dft_control%nspins
707 
708  ! Initialize vectors:
709  ! psi0 : The ground-state MO-coefficients
710  ! psi1 : The "perturbed" linear response orbitals
711  ALLOCATE (psi0(nspins), psi1(nspins))
712  DO ispin = 1, nspins
713  CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, homo=nocc)
714  NULLIFY (fm_struct)
715  CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, &
716  template_fmstruct=mo_coeff%matrix_struct)
717  CALL cp_fm_create(psi0(ispin), fm_struct)
718  CALL cp_fm_to_fm(mo_coeff, psi0(ispin), nocc)
719  CALL cp_fm_create(psi1(ispin), fm_struct)
720  CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
721  CALL cp_fm_struct_release(fm_struct)
722  END DO
723 
724  should_stop = .false.
725  ! The response solver
726  CALL get_qs_env(qs_env, matrix_s=matrix_s, sab_orb=sab_orb)
727  CALL dbcsr_allocate_matrix_set(p_env%p1, nspins)
728  CALL dbcsr_allocate_matrix_set(p_env%w1, nspins)
729  DO ispin = 1, nspins
730  ALLOCATE (p_env%p1(ispin)%matrix, p_env%w1(ispin)%matrix)
731  CALL dbcsr_create(matrix=p_env%p1(ispin)%matrix, template=matrix_s(1)%matrix)
732  CALL dbcsr_create(matrix=p_env%w1(ispin)%matrix, template=matrix_s(1)%matrix)
733  CALL cp_dbcsr_alloc_block_from_nbl(p_env%p1(ispin)%matrix, sab_orb)
734  CALL cp_dbcsr_alloc_block_from_nbl(p_env%w1(ispin)%matrix, sab_orb)
735  END DO
736  IF (dft_control%do_admm) THEN
737  CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux)
738  CALL dbcsr_allocate_matrix_set(p_env%p1_admm, nspins)
739  DO ispin = 1, nspins
740  ALLOCATE (p_env%p1_admm(ispin)%matrix)
741  CALL dbcsr_create(p_env%p1_admm(ispin)%matrix, &
742  template=matrix_s_aux(1)%matrix)
743  CALL dbcsr_copy(p_env%p1_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
744  CALL dbcsr_set(p_env%p1_admm(ispin)%matrix, 0.0_dp)
745  END DO
746  END IF
747 
748  CALL linres_solver(p_env, qs_env, psi1, cpmos, psi0, iounit, should_stop)
749 
750  ! Building the response density matrix
751  DO ispin = 1, nspins
752  CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_s(1)%matrix)
753  END DO
754  CALL build_dm_response(psi0, psi1, p_env%p1)
755  DO ispin = 1, nspins
756  CALL dbcsr_scale(p_env%p1(ispin)%matrix, 0.5_dp)
757  END DO
758  IF (dft_control%do_admm) THEN
759  CALL get_qs_env(qs_env, admm_env=admm_env)
760  cpassert(ASSOCIATED(admm_env%work_orb_orb))
761  cpassert(ASSOCIATED(admm_env%work_aux_orb))
762  cpassert(ASSOCIATED(admm_env%work_aux_aux))
763  nao = admm_env%nao_orb
764  nao_aux = admm_env%nao_aux_fit
765  DO ispin = 1, nspins
766  CALL copy_dbcsr_to_fm(p_env%p1(ispin)%matrix, admm_env%work_orb_orb)
767  CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
768  1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
769  admm_env%work_aux_orb)
770  CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
771  1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
772  admm_env%work_aux_aux)
773  CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, p_env%p1_admm(ispin)%matrix, &
774  keep_sparsity=.true.)
775  END DO
776  END IF
777 
778  ! Calculate Wz = 0.5*(psi1*eps*psi0^T + psi0*eps*psi1^T)
779  CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
780  DO ispin = 1, nspins
781  CALL calculate_wz_matrix(mos(ispin), psi1(ispin), matrix_ks(ispin)%matrix, &
782  p_env%w1(ispin)%matrix)
783  END DO
784  CALL cp_fm_release(psi0)
785  CALL cp_fm_release(psi1)
786 
787  CALL timestop(handle)
788 
789  END SUBROUTINE response_equation
790 
791 ! **************************************************************************************************
792 !> \brief ...
793 !> \param qs_env ...
794 !> \param vh_rspace ...
795 !> \param vxc_rspace ...
796 !> \param vtau_rspace ...
797 !> \param vadmm_rspace ...
798 !> \param matrix_hz Right-hand-side of linear response equation
799 !> \param matrix_pz Linear response density matrix
800 !> \param matrix_pz_admm Linear response density matrix in ADMM basis
801 !> \param matrix_wz Energy-weighted linear response density
802 !> \param zehartree Hartree volume response contribution to stress tensor
803 !> \param zexc XC volume response contribution to stress tensor
804 !> \param zexc_aux_fit ADMM XC volume response contribution to stress tensor
805 !> \param rhopz_r Response density on real space grid
806 !> \param p_env ...
807 !> \param ex_env ...
808 !> \param debug ...
809 ! **************************************************************************************************
810  SUBROUTINE response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, &
811  matrix_hz, matrix_pz, matrix_pz_admm, matrix_wz, &
812  zehartree, zexc, zexc_aux_fit, rhopz_r, p_env, ex_env, debug)
813  TYPE(qs_environment_type), POINTER :: qs_env
814  TYPE(pw_r3d_rs_type), INTENT(IN) :: vh_rspace
815  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rspace, vtau_rspace, vadmm_rspace
816  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz, matrix_pz, matrix_pz_admm, &
817  matrix_wz
818  REAL(kind=dp), OPTIONAL :: zehartree, zexc, zexc_aux_fit
819  TYPE(pw_r3d_rs_type), DIMENSION(:), &
820  INTENT(INOUT), OPTIONAL :: rhopz_r
821  TYPE(qs_p_env_type), OPTIONAL :: p_env
822  TYPE(excited_energy_type), OPTIONAL, POINTER :: ex_env
823  LOGICAL, INTENT(IN), OPTIONAL :: debug
824 
825  CHARACTER(LEN=*), PARAMETER :: routinen = 'response_force'
826 
827  CHARACTER(LEN=default_string_length) :: basis_type
828  INTEGER :: handle, iounit, ispin, mspin, myfun, &
829  n_rep_hf, nao, nao_aux, natom, nder, &
830  nimages, nocc, nspins
831  INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
832  LOGICAL :: debug_forces, debug_stress, distribute_fock_matrix, do_ex, do_hfx, gapw, gapw_xc, &
833  hfx_treat_lsd_in_core, resp_only, s_mstruct_changed, use_virial
834  REAL(kind=dp) :: eh1, ehartree, ekin_mol, eps_filter, &
835  eps_ppnl, exc, exc_aux_fit, fconv, &
836  focc, hartree_gs, hartree_t
837  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ftot1, ftot2, ftot3
838  REAL(kind=dp), DIMENSION(2) :: total_rho_gs, total_rho_t
839  REAL(kind=dp), DIMENSION(3) :: fodeb
840  REAL(kind=dp), DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot, sttot2
841  TYPE(admm_type), POINTER :: admm_env
842  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
843  TYPE(cell_type), POINTER :: cell
844  TYPE(cp_logger_type), POINTER :: logger
845  TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
846  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ht, matrix_pd, matrix_pza, &
847  matrix_s, mpa, scrm
848  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, mhd, mhx, mhy, mhz, &
849  mpd, mpz
850  TYPE(dbcsr_type), POINTER :: dbwork
851  TYPE(dft_control_type), POINTER :: dft_control
852  TYPE(hartree_local_type), POINTER :: hartree_local_gs, hartree_local_t
853  TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
854  TYPE(kg_environment_type), POINTER :: kg_env
855  TYPE(local_rho_type), POINTER :: local_rho_set_f, local_rho_set_gs, &
856  local_rho_set_t, local_rho_set_vxc, &
857  local_rhoz_set_admm
858  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
859  TYPE(mp_para_env_type), POINTER :: para_env
860  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
861  POINTER :: sab_aux_fit, sab_orb, sac_ae, sac_ppl, &
862  sap_ppnl
863  TYPE(oce_matrix_type), POINTER :: oce
864  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
865  TYPE(pw_c1d_gs_type) :: rho_tot_gspace, rho_tot_gspace_gs, rho_tot_gspace_t, &
866  rhoz_tot_gspace, v_hartree_gspace_gs, v_hartree_gspace_t, zv_hartree_gspace
867  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux, rho_g_gs, rho_g_t, rhoz_g, &
868  rhoz_g_aux, rhoz_g_xc
869  TYPE(pw_c1d_gs_type), POINTER :: rho_core
870  TYPE(pw_env_type), POINTER :: pw_env
871  TYPE(pw_poisson_type), POINTER :: poisson_env
872  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
873  TYPE(pw_r3d_rs_type) :: v_hartree_rspace_gs, v_hartree_rspace_t, &
874  vhxc_rspace, zv_hartree_rspace
875  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_aux, rho_r_gs, rho_r_t, rhoz_r, &
876  rhoz_r_aux, rhoz_r_xc, tau_r_aux, &
877  tauz_r, tauz_r_xc, v_xc, v_xc_tau
878  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
879  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
880  TYPE(qs_ks_env_type), POINTER :: ks_env
881  TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit, rho_xc
882  TYPE(section_vals_type), POINTER :: hfx_section, xc_fun_section, xc_section
883  TYPE(task_list_type), POINTER :: task_list, task_list_aux_fit
884  TYPE(virial_type), POINTER :: virial
885 
886  CALL timeset(routinen, handle)
887 
888  IF (PRESENT(debug)) THEN
889  debug_forces = debug
890  debug_stress = debug
891  ELSE
892  debug_forces = .false.
893  debug_stress = .false.
894  END IF
895 
896  logger => cp_get_default_logger()
897  logger => cp_get_default_logger()
898  IF (logger%para_env%is_source()) THEN
899  iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
900  ELSE
901  iounit = -1
902  END IF
903 
904  do_ex = .false.
905  IF (PRESENT(ex_env)) do_ex = .true.
906  IF (do_ex) THEN
907  cpassert(PRESENT(p_env))
908  END IF
909 
910  NULLIFY (ks_env, sab_orb, sac_ae, sac_ppl, sap_ppnl, virial)
911  CALL get_qs_env(qs_env=qs_env, &
912  cell=cell, &
913  force=force, &
914  ks_env=ks_env, &
915  dft_control=dft_control, &
916  para_env=para_env, &
917  sab_orb=sab_orb, &
918  sac_ae=sac_ae, &
919  sac_ppl=sac_ppl, &
920  sap_ppnl=sap_ppnl, &
921  virial=virial)
922  nspins = dft_control%nspins
923  gapw = dft_control%qs_control%gapw
924  gapw_xc = dft_control%qs_control%gapw_xc
925 
926  IF (debug_forces) THEN
927  CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
928  ALLOCATE (ftot1(3, natom))
929  CALL total_qs_force(ftot1, force, atomic_kind_set)
930  END IF
931 
932  ! check for virial
933  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
934 
935  fconv = 1.0e-9_dp*pascal/cell%deth
936  IF (debug_stress .AND. use_virial) THEN
937  sttot = virial%pv_virial
938  END IF
939 
940  ! *** If LSD, then combine alpha density and beta density to
941  ! *** total density: alpha <- alpha + beta and
942  NULLIFY (mpa)
943  NULLIFY (matrix_ht)
944  IF (do_ex) THEN
945  CALL dbcsr_allocate_matrix_set(mpa, nspins)
946  DO ispin = 1, nspins
947  ALLOCATE (mpa(ispin)%matrix)
948  CALL dbcsr_create(mpa(ispin)%matrix, template=p_env%p1(ispin)%matrix)
949  CALL dbcsr_copy(mpa(ispin)%matrix, p_env%p1(ispin)%matrix)
950  CALL dbcsr_add(mpa(ispin)%matrix, ex_env%matrix_pe(ispin)%matrix, 1.0_dp, 1.0_dp)
951  CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
952  END DO
953 
954  CALL dbcsr_allocate_matrix_set(matrix_ht, nspins)
955  DO ispin = 1, nspins
956  ALLOCATE (matrix_ht(ispin)%matrix)
957  CALL dbcsr_create(matrix_ht(ispin)%matrix, template=matrix_hz(ispin)%matrix)
958  CALL dbcsr_copy(matrix_ht(ispin)%matrix, matrix_hz(ispin)%matrix)
959  CALL dbcsr_set(matrix_ht(ispin)%matrix, 0.0_dp)
960  END DO
961  ELSE
962  mpa => matrix_pz
963  END IF
964  !
965  ! START OF Tr(P+Z)Hcore
966  !
967  IF (nspins == 2) THEN
968  CALL dbcsr_add(mpa(1)%matrix, mpa(2)%matrix, 1.0_dp, 1.0_dp)
969  END IF
970  nimages = 1
971 
972  ! Kinetic energy matrix
973  NULLIFY (scrm)
974  IF (debug_forces) fodeb(1:3) = force(1)%kinetic(1:3, 1)
975  IF (debug_stress .AND. use_virial) stdeb = virial%pv_ekinetic
976  CALL build_kinetic_matrix(ks_env, matrix_t=scrm, &
977  matrix_name="KINETIC ENERGY MATRIX", &
978  basis_type="ORB", &
979  sab_nl=sab_orb, calculate_forces=.true., &
980  matrix_p=mpa(1)%matrix)
981  IF (debug_forces) THEN
982  fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
983  CALL para_env%sum(fodeb)
984  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dT ", fodeb
985  END IF
986  IF (debug_stress .AND. use_virial) THEN
987  stdeb = fconv*(virial%pv_ekinetic - stdeb)
988  CALL para_env%sum(stdeb)
989  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
990  'STRESS| Kinetic energy', one_third_sum_diag(stdeb), det_3x3(stdeb)
991  END IF
992 
993  IF (nspins == 2) THEN
994  CALL dbcsr_add(mpa(1)%matrix, mpa(2)%matrix, 1.0_dp, -1.0_dp)
995  END IF
996  CALL dbcsr_deallocate_matrix_set(scrm)
997 
998  ! Initialize a matrix scrm, later used for scratch purposes
999  CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
1000  CALL dbcsr_allocate_matrix_set(scrm, nspins)
1001  DO ispin = 1, nspins
1002  ALLOCATE (scrm(ispin)%matrix)
1003  CALL dbcsr_create(scrm(ispin)%matrix, template=matrix_s(1)%matrix)
1004  CALL dbcsr_copy(scrm(ispin)%matrix, matrix_s(1)%matrix)
1005  CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
1006  END DO
1007 
1008  CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
1009  atomic_kind_set=atomic_kind_set)
1010 
1011  NULLIFY (cell_to_index)
1012  ALLOCATE (matrix_p(nspins, 1), matrix_h(nspins, 1))
1013  DO ispin = 1, nspins
1014  matrix_p(ispin, 1)%matrix => mpa(ispin)%matrix
1015  matrix_h(ispin, 1)%matrix => scrm(ispin)%matrix
1016  END DO
1017  matrix_h(1, 1)%matrix => scrm(1)%matrix
1018 
1019  IF (ASSOCIATED(sac_ae)) THEN
1020  nder = 1
1021  IF (debug_forces) fodeb(1:3) = force(1)%all_potential(1:3, 1)
1022  IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl
1023  CALL build_core_ae(matrix_h, matrix_p, force, virial, .true., use_virial, nder, &
1024  qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
1025  nimages, cell_to_index)
1026  IF (debug_forces) THEN
1027  fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3)
1028  CALL para_env%sum(fodeb)
1029  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHae ", fodeb
1030  END IF
1031  IF (debug_stress .AND. use_virial) THEN
1032  stdeb = fconv*(virial%pv_ppl - stdeb)
1033  CALL para_env%sum(stdeb)
1034  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1035  'STRESS| Pz*dHae ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1036  END IF
1037  END IF
1038 
1039  IF (ASSOCIATED(sac_ppl)) THEN
1040  nder = 1
1041  IF (debug_forces) fodeb(1:3) = force(1)%gth_ppl(1:3, 1)
1042  IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppl
1043  CALL build_core_ppl(matrix_h, matrix_p, force, &
1044  virial, .true., use_virial, nder, &
1045  qs_kind_set, atomic_kind_set, particle_set, &
1046  sab_orb, sac_ppl, nimages, cell_to_index, "ORB")
1047  IF (debug_forces) THEN
1048  fodeb(1:3) = force(1)%gth_ppl(1:3, 1) - fodeb(1:3)
1049  CALL para_env%sum(fodeb)
1050  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHppl ", fodeb
1051  END IF
1052  IF (debug_stress .AND. use_virial) THEN
1053  stdeb = fconv*(virial%pv_ppl - stdeb)
1054  CALL para_env%sum(stdeb)
1055  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1056  'STRESS| Pz*dHppl ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1057  END IF
1058  END IF
1059  eps_ppnl = dft_control%qs_control%eps_ppnl
1060  IF (ASSOCIATED(sap_ppnl)) THEN
1061  nder = 1
1062  IF (debug_forces) fodeb(1:3) = force(1)%gth_ppnl(1:3, 1)
1063  IF (debug_stress .AND. use_virial) stdeb = virial%pv_ppnl
1064  CALL build_core_ppnl(matrix_h, matrix_p, force, &
1065  virial, .true., use_virial, nder, &
1066  qs_kind_set, atomic_kind_set, particle_set, &
1067  sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, "ORB")
1068  IF (debug_forces) THEN
1069  fodeb(1:3) = force(1)%gth_ppnl(1:3, 1) - fodeb(1:3)
1070  CALL para_env%sum(fodeb)
1071  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHppnl ", fodeb
1072  END IF
1073  IF (debug_stress .AND. use_virial) THEN
1074  stdeb = fconv*(virial%pv_ppnl - stdeb)
1075  CALL para_env%sum(stdeb)
1076  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1077  'STRESS| Pz*dHppnl ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1078  END IF
1079 
1080  END IF
1081  ! Kim-Gordon subsystem DFT
1082  ! Atomic potential for nonadditive kinetic energy contribution
1083  IF (dft_control%qs_control%do_kg) THEN
1084  IF (qs_env%kg_env%tnadd_method == kg_tnadd_atomic) THEN
1085  CALL get_qs_env(qs_env=qs_env, kg_env=kg_env, dbcsr_dist=dbcsr_dist)
1086 
1087  IF (use_virial) THEN
1088  pv_loc = virial%pv_virial
1089  END IF
1090 
1091  IF (debug_forces) fodeb(1:3) = force(1)%kinetic(1:3, 1)
1092  IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1093  CALL build_tnadd_mat(kg_env=kg_env, matrix_p=matrix_p, force=force, virial=virial, &
1094  calculate_forces=.true., use_virial=use_virial, &
1095  qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
1096  particle_set=particle_set, sab_orb=sab_orb, dbcsr_dist=dbcsr_dist)
1097  IF (debug_forces) THEN
1098  fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
1099  CALL para_env%sum(fodeb)
1100  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dTnadd ", fodeb
1101  END IF
1102  IF (debug_stress .AND. use_virial) THEN
1103  stdeb = fconv*(virial%pv_virial - stdeb)
1104  CALL para_env%sum(stdeb)
1105  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1106  'STRESS| Pz*dTnadd ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1107  END IF
1108 
1109  ! Stress-tensor update components
1110  IF (use_virial) THEN
1111  virial%pv_ekinetic = virial%pv_ekinetic + (virial%pv_virial - pv_loc)
1112  END IF
1113 
1114  END IF
1115  END IF
1116 
1117  DEALLOCATE (matrix_h)
1118  DEALLOCATE (matrix_p)
1119  CALL dbcsr_deallocate_matrix_set(scrm)
1120 
1121  ! initialize src matrix
1122  ! Necessary as build_kinetic_matrix will only allocate scrm(1)
1123  ! and not scrm(2) in open-shell case
1124  NULLIFY (scrm)
1125  CALL dbcsr_allocate_matrix_set(scrm, nspins)
1126  DO ispin = 1, nspins
1127  ALLOCATE (scrm(ispin)%matrix)
1128  CALL dbcsr_create(scrm(ispin)%matrix, template=matrix_pz(1)%matrix)
1129  CALL dbcsr_copy(scrm(ispin)%matrix, matrix_pz(ispin)%matrix)
1130  CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
1131  END DO
1132 
1133  IF (debug_forces) THEN
1134  ALLOCATE (ftot2(3, natom))
1135  CALL total_qs_force(ftot2, force, atomic_kind_set)
1136  fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
1137  CALL para_env%sum(fodeb)
1138  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHcore", fodeb
1139  END IF
1140  IF (debug_stress .AND. use_virial) THEN
1141  stdeb = fconv*(virial%pv_virial - sttot)
1142  CALL para_env%sum(stdeb)
1143  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1144  'STRESS| Stress Pz*dHcore ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1145  ! save current total viral, does not contain volume terms yet
1146  sttot2 = virial%pv_virial
1147  END IF
1148  !
1149  ! END OF Tr(P+Z)Hcore
1150  !
1151  !
1152  ! Vhxc (KS potentials calculated externally)
1153  CALL get_qs_env(qs_env, pw_env=pw_env)
1154  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
1155  !
1156  IF (dft_control%do_admm) THEN
1157  CALL get_qs_env(qs_env, admm_env=admm_env)
1158  xc_section => admm_env%xc_section_primary
1159  ELSE
1160  xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
1161  END IF
1162  xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1163  CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
1164  !
1165  IF (gapw .OR. gapw_xc) THEN
1166  NULLIFY (oce, sab_orb)
1167  CALL get_qs_env(qs_env=qs_env, oce=oce, sab_orb=sab_orb)
1168  ! set up local_rho_set for GS density
1169  NULLIFY (local_rho_set_gs)
1170  CALL get_qs_env(qs_env=qs_env, rho=rho)
1171  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1172  CALL local_rho_set_create(local_rho_set_gs)
1173  CALL allocate_rho_atom_internals(local_rho_set_gs%rho_atom_set, atomic_kind_set, &
1174  qs_kind_set, dft_control, para_env)
1175  CALL init_rho0(local_rho_set_gs, qs_env, dft_control%qs_control%gapw_control)
1176  CALL rho0_s_grid_create(pw_env, local_rho_set_gs%rho0_mpole)
1177  CALL calculate_rho_atom_coeff(qs_env, matrix_p(:, 1), local_rho_set_gs%rho_atom_set, &
1178  qs_kind_set, oce, sab_orb, para_env)
1179  CALL prepare_gapw_den(qs_env, local_rho_set_gs, do_rho0=gapw)
1180  ! set up local_rho_set for response density
1181  NULLIFY (local_rho_set_t)
1182  CALL local_rho_set_create(local_rho_set_t)
1183  CALL allocate_rho_atom_internals(local_rho_set_t%rho_atom_set, atomic_kind_set, &
1184  qs_kind_set, dft_control, para_env)
1185  CALL init_rho0(local_rho_set_t, qs_env, dft_control%qs_control%gapw_control, &
1186  zcore=0.0_dp)
1187  CALL rho0_s_grid_create(pw_env, local_rho_set_t%rho0_mpole)
1188  CALL calculate_rho_atom_coeff(qs_env, mpa(:), local_rho_set_t%rho_atom_set, &
1189  qs_kind_set, oce, sab_orb, para_env)
1190  CALL prepare_gapw_den(qs_env, local_rho_set_t, do_rho0=gapw)
1191 
1192  ! compute soft GS potential
1193  ALLOCATE (rho_r_gs(nspins), rho_g_gs(nspins))
1194  DO ispin = 1, nspins
1195  CALL auxbas_pw_pool%create_pw(rho_r_gs(ispin))
1196  CALL auxbas_pw_pool%create_pw(rho_g_gs(ispin))
1197  END DO
1198  CALL auxbas_pw_pool%create_pw(rho_tot_gspace_gs)
1199  ! compute soft GS density
1200  total_rho_gs = 0.0_dp
1201  CALL pw_zero(rho_tot_gspace_gs)
1202  DO ispin = 1, nspins
1203  CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_p(ispin, 1)%matrix, &
1204  rho=rho_r_gs(ispin), &
1205  rho_gspace=rho_g_gs(ispin), &
1206  soft_valid=(gapw .OR. gapw_xc), &
1207  total_rho=total_rho_gs(ispin))
1208  CALL pw_axpy(rho_g_gs(ispin), rho_tot_gspace_gs)
1209  END DO
1210  IF (gapw) THEN
1211  CALL get_qs_env(qs_env, natom=natom)
1212  ! add rho0 contributions to GS density (only for Coulomb) only for gapw
1213  CALL pw_axpy(local_rho_set_gs%rho0_mpole%rho0_s_gs, rho_tot_gspace_gs)
1214  IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
1215  CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
1216  CALL pw_axpy(rho_core, rho_tot_gspace_gs)
1217  END IF
1218  ! compute GS potential
1219  CALL auxbas_pw_pool%create_pw(v_hartree_gspace_gs)
1220  CALL auxbas_pw_pool%create_pw(v_hartree_rspace_gs)
1221  NULLIFY (hartree_local_gs)
1222  CALL hartree_local_create(hartree_local_gs)
1223  CALL init_coulomb_local(hartree_local_gs, natom)
1224  CALL pw_poisson_solve(poisson_env, rho_tot_gspace_gs, hartree_gs, v_hartree_gspace_gs)
1225  CALL pw_transfer(v_hartree_gspace_gs, v_hartree_rspace_gs)
1226  CALL pw_scale(v_hartree_rspace_gs, v_hartree_rspace_gs%pw_grid%dvol)
1227  END IF
1228  END IF
1229 
1230  IF (gapw) THEN
1231  ! Hartree grid PAW term
1232  cpassert(do_ex)
1233  cpassert(.NOT. use_virial)
1234  IF (debug_forces) fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1)
1235  CALL vh_1c_gg_integrals(qs_env, hartree_gs, hartree_local_gs%ecoul_1c, local_rho_set_t, para_env, tddft=.true., &
1236  local_rho_set_2nd=local_rho_set_gs, core_2nd=.false.) ! n^core for GS potential
1237  ! 1st to define integral space, 2nd for potential, integral contributions stored on local_rho_set_gs
1238  CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace_gs, para_env, calculate_forces=.true., &
1239  local_rho_set=local_rho_set_t, local_rho_set_2nd=local_rho_set_gs)
1240  IF (debug_forces) THEN
1241  fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1) - fodeb(1:3)
1242  CALL para_env%sum(fodeb)
1243  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (T+Dz)*dVh[D^GS]PAWg0", fodeb
1244  END IF
1245  END IF
1246  IF (gapw .OR. gapw_xc) THEN
1247  IF (myfun /= xc_none) THEN
1248  ! add 1c hard and soft XC contributions
1249  NULLIFY (local_rho_set_vxc)
1250  CALL local_rho_set_create(local_rho_set_vxc)
1251  CALL allocate_rho_atom_internals(local_rho_set_vxc%rho_atom_set, atomic_kind_set, &
1252  qs_kind_set, dft_control, para_env)
1253  CALL calculate_rho_atom_coeff(qs_env, matrix_p(:, 1), local_rho_set_vxc%rho_atom_set, &
1254  qs_kind_set, oce, sab_orb, para_env)
1255  CALL prepare_gapw_den(qs_env, local_rho_set_vxc, do_rho0=.false.)
1256  ! compute hard and soft atomic contributions
1257  CALL calculate_vxc_atom(qs_env, .false., exc1=hartree_gs, xc_section_external=xc_section, &
1258  rho_atom_set_external=local_rho_set_vxc%rho_atom_set)
1259  END IF ! myfun
1260  END IF ! gapw
1261 
1262  CALL auxbas_pw_pool%create_pw(vhxc_rspace)
1263  !
1264  ! Stress-tensor: integration contribution direct term
1265  ! int v_Hxc[n^in]*n^z
1266  IF (use_virial) THEN
1267  pv_loc = virial%pv_virial
1268  END IF
1269 
1270  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1271  IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1272  IF (gapw .OR. gapw_xc) THEN
1273  ! vtot = v_xc + v_hartree
1274  DO ispin = 1, nspins
1275  CALL pw_zero(vhxc_rspace)
1276  IF (gapw) THEN
1277  CALL pw_transfer(v_hartree_rspace_gs, vhxc_rspace)
1278  ELSEIF (gapw_xc) THEN
1279  CALL pw_transfer(vh_rspace, vhxc_rspace)
1280  END IF
1281  CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1282  hmat=scrm(ispin), pmat=mpa(ispin), &
1283  qs_env=qs_env, gapw=gapw, &
1284  calculate_forces=.true.)
1285  END DO
1286  IF (myfun /= xc_none) THEN
1287  DO ispin = 1, nspins
1288  CALL pw_zero(vhxc_rspace)
1289  CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
1290  CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1291  hmat=scrm(ispin), pmat=mpa(ispin), &
1292  qs_env=qs_env, gapw=(gapw .OR. gapw_xc), &
1293  calculate_forces=.true.)
1294  END DO
1295  END IF
1296  ELSE ! original GPW with Standard Hartree as Potential
1297  DO ispin = 1, nspins
1298  CALL pw_transfer(vh_rspace, vhxc_rspace)
1299  CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
1300  CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1301  hmat=scrm(ispin), pmat=mpa(ispin), &
1302  qs_env=qs_env, gapw=gapw, calculate_forces=.true.)
1303  END DO
1304  END IF
1305 
1306  IF (debug_forces) THEN
1307  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1308  CALL para_env%sum(fodeb)
1309  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (T+Dz)*dVhxc[D^GS] ", fodeb
1310  END IF
1311  IF (debug_stress .AND. use_virial) THEN
1312  stdeb = fconv*(virial%pv_virial - pv_loc)
1313  CALL para_env%sum(stdeb)
1314  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1315  'STRESS| INT Pz*dVhxc ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1316  END IF
1317 
1318  IF (gapw .OR. gapw_xc) THEN
1319  cpassert(do_ex)
1320  ! HXC term
1321  IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
1322  IF (gapw) CALL update_ks_atom(qs_env, scrm, mpa, forces=.true., tddft=.false., &
1323  rho_atom_external=local_rho_set_gs%rho_atom_set)
1324  IF (myfun /= xc_none) CALL update_ks_atom(qs_env, scrm, mpa, forces=.true., tddft=.false., &
1325  rho_atom_external=local_rho_set_vxc%rho_atom_set)
1326  IF (debug_forces) THEN
1327  fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
1328  CALL para_env%sum(fodeb)
1329  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (T+Dz)*dVhxc[D^GS]PAW ", fodeb
1330  END IF
1331  ! release local environments for GAPW
1332  IF (myfun /= xc_none) THEN
1333  IF (ASSOCIATED(local_rho_set_vxc)) CALL local_rho_set_release(local_rho_set_vxc)
1334  END IF
1335  IF (ASSOCIATED(local_rho_set_gs)) CALL local_rho_set_release(local_rho_set_gs)
1336  IF (gapw) THEN
1337  IF (ASSOCIATED(hartree_local_gs)) CALL hartree_local_release(hartree_local_gs)
1338  CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace_gs)
1339  CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace_gs)
1340  END IF
1341  CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace_gs)
1342  IF (ASSOCIATED(rho_r_gs)) THEN
1343  DO ispin = 1, nspins
1344  CALL auxbas_pw_pool%give_back_pw(rho_r_gs(ispin))
1345  END DO
1346  DEALLOCATE (rho_r_gs)
1347  END IF
1348  IF (ASSOCIATED(rho_g_gs)) THEN
1349  DO ispin = 1, nspins
1350  CALL auxbas_pw_pool%give_back_pw(rho_g_gs(ispin))
1351  END DO
1352  DEALLOCATE (rho_g_gs)
1353  END IF
1354  END IF !gapw
1355 
1356  IF (ASSOCIATED(vtau_rspace)) THEN
1357  cpassert(.NOT. (gapw .OR. gapw_xc))
1358  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1359  IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1360  DO ispin = 1, nspins
1361  CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
1362  hmat=scrm(ispin), pmat=mpa(ispin), &
1363  qs_env=qs_env, gapw=(gapw .OR. gapw_xc), &
1364  calculate_forces=.true., compute_tau=.true.)
1365  END DO
1366  IF (debug_forces) THEN
1367  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1368  CALL para_env%sum(fodeb)
1369  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dVxc_tau ", fodeb
1370  END IF
1371  IF (debug_stress .AND. use_virial) THEN
1372  stdeb = fconv*(virial%pv_virial - pv_loc)
1373  CALL para_env%sum(stdeb)
1374  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1375  'STRESS| INT Pz*dVxc_tau ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1376  END IF
1377  END IF
1378  CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
1379 
1380  ! Stress-tensor Pz*v_Hxc[Pin]
1381  IF (use_virial) THEN
1382  virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1383  END IF
1384 
1385  ! KG Embedding
1386  ! calculate kinetic energy potential and integrate with response density
1387  IF (dft_control%qs_control%do_kg) THEN
1388  IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed .OR. &
1389  qs_env%kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
1390 
1391  ekin_mol = 0.0_dp
1392  IF (use_virial) THEN
1393  pv_loc = virial%pv_virial
1394  END IF
1395 
1396  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1397  CALL kg_ekin_subset(qs_env=qs_env, &
1398  ks_matrix=scrm, &
1399  ekin_mol=ekin_mol, &
1400  calc_force=.true., &
1401  do_kernel=.false., &
1402  pmat_ext=mpa)
1403  IF (debug_forces) THEN
1404  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1405  CALL para_env%sum(fodeb)
1406  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dVkg ", fodeb
1407  END IF
1408  IF (debug_stress .AND. use_virial) THEN
1409  !IF (iounit > 0) WRITE(iounit, *) &
1410  ! "response_force | VOL 1st KG - v_KG[n_in]*n_z: ", ekin_mol
1411  stdeb = 1.0_dp*fconv*ekin_mol
1412  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1413  'STRESS| VOL KG Pz*dVKG ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1414 
1415  stdeb = fconv*(virial%pv_virial - pv_loc)
1416  CALL para_env%sum(stdeb)
1417  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1418  'STRESS| INT KG Pz*dVKG ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1419 
1420  stdeb = fconv*virial%pv_xc
1421  CALL para_env%sum(stdeb)
1422  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1423  'STRESS| GGA KG Pz*dVKG ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1424  END IF
1425  IF (use_virial) THEN
1426  ! Direct integral contribution
1427  virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1428  END IF
1429 
1430  END IF ! tnadd_method
1431  END IF ! do_kg
1432 
1433  CALL dbcsr_deallocate_matrix_set(scrm)
1434 
1435  !
1436  ! Hartree potential of response density
1437  !
1438  ALLOCATE (rhoz_r(nspins), rhoz_g(nspins))
1439  DO ispin = 1, nspins
1440  CALL auxbas_pw_pool%create_pw(rhoz_r(ispin))
1441  CALL auxbas_pw_pool%create_pw(rhoz_g(ispin))
1442  END DO
1443  CALL auxbas_pw_pool%create_pw(rhoz_tot_gspace)
1444  CALL auxbas_pw_pool%create_pw(zv_hartree_rspace)
1445  CALL auxbas_pw_pool%create_pw(zv_hartree_gspace)
1446 
1447  CALL pw_zero(rhoz_tot_gspace)
1448  DO ispin = 1, nspins
1449  CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1450  rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin), &
1451  soft_valid=gapw)
1452  CALL pw_axpy(rhoz_g(ispin), rhoz_tot_gspace)
1453  END DO
1454  IF (gapw_xc) THEN
1455  NULLIFY (tauz_r_xc)
1456  ALLOCATE (rhoz_r_xc(nspins), rhoz_g_xc(nspins))
1457  DO ispin = 1, nspins
1458  CALL auxbas_pw_pool%create_pw(rhoz_r_xc(ispin))
1459  CALL auxbas_pw_pool%create_pw(rhoz_g_xc(ispin))
1460  END DO
1461  DO ispin = 1, nspins
1462  CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1463  rho=rhoz_r_xc(ispin), rho_gspace=rhoz_g_xc(ispin), &
1464  soft_valid=gapw_xc)
1465  END DO
1466  END IF
1467 
1468  IF (ASSOCIATED(vtau_rspace)) THEN
1469  cpassert(.NOT. (gapw .OR. gapw_xc))
1470  block
1471  TYPE(pw_c1d_gs_type) :: work_g
1472  ALLOCATE (tauz_r(nspins))
1473  CALL auxbas_pw_pool%create_pw(work_g)
1474  DO ispin = 1, nspins
1475  CALL auxbas_pw_pool%create_pw(tauz_r(ispin))
1476  CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1477  rho=tauz_r(ispin), rho_gspace=work_g, &
1478  compute_tau=.true.)
1479  END DO
1480  CALL auxbas_pw_pool%give_back_pw(work_g)
1481  END block
1482  END IF
1483 
1484  !
1485  IF (PRESENT(rhopz_r)) THEN
1486  DO ispin = 1, nspins
1487  CALL pw_copy(rhoz_r(ispin), rhopz_r(ispin))
1488  END DO
1489  END IF
1490 
1491  ! Stress-tensor contribution second derivative
1492  ! Volume : int v_H[n^z]*n_in
1493  ! Volume : int epsilon_xc*n_z
1494  IF (use_virial) THEN
1495 
1496  CALL get_qs_env(qs_env, rho=rho)
1497  CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
1498 
1499  ! Get the total input density in g-space [ions + electrons]
1500  CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
1501 
1502  h_stress(:, :) = 0.0_dp
1503  ! calculate associated hartree potential
1504  ! This term appears twice in the derivation of the equations
1505  ! v_H[n_in]*n_z and v_H[n_z]*n_in
1506  ! due to symmetry we only need to call this routine once,
1507  ! and count the Volume and Green function contribution
1508  ! which is stored in h_stress twice
1509  CALL pw_poisson_solve(poisson_env, &
1510  density=rhoz_tot_gspace, & ! n_z
1511  ehartree=ehartree, &
1512  vhartree=zv_hartree_gspace, & ! v_H[n_z]
1513  h_stress=h_stress, &
1514  aux_density=rho_tot_gspace) ! n_in
1515 
1516  CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1517 
1518  ! Stress tensor Green function contribution
1519  virial%pv_ehartree = virial%pv_ehartree + 2.0_dp*h_stress/real(para_env%num_pe, dp)
1520  virial%pv_virial = virial%pv_virial + 2.0_dp*h_stress/real(para_env%num_pe, dp)
1521 
1522  IF (debug_stress) THEN
1523  stdeb = -1.0_dp*fconv*ehartree
1524  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1525  'STRESS| VOL 1st v_H[n_z]*n_in ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1526  stdeb = -1.0_dp*fconv*ehartree
1527  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1528  'STRESS| VOL 2nd v_H[n_in]*n_z ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1529  stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
1530  CALL para_env%sum(stdeb)
1531  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1532  'STRESS| GREEN 1st v_H[n_z]*n_in ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1533  stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
1534  CALL para_env%sum(stdeb)
1535  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1536  'STRESS| GREEN 2nd v_H[n_in]*n_z ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1537  END IF
1538 
1539  ! Stress tensor volume term: \int v_xc[n_in]*n_z
1540  ! vxc_rspace already scaled, we need to unscale it!
1541  exc = 0.0_dp
1542  DO ispin = 1, nspins
1543  exc = exc + pw_integral_ab(rhoz_r(ispin), vxc_rspace(ispin))/ &
1544  vxc_rspace(ispin)%pw_grid%dvol
1545  END DO
1546  IF (ASSOCIATED(vtau_rspace)) THEN
1547  DO ispin = 1, nspins
1548  exc = exc + pw_integral_ab(tauz_r(ispin), vtau_rspace(ispin))/ &
1549  vtau_rspace(ispin)%pw_grid%dvol
1550  END DO
1551  END IF
1552 
1553  ! Add KG embedding correction
1554  IF (dft_control%qs_control%do_kg) THEN
1555  IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed .OR. &
1556  qs_env%kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
1557  exc = exc - ekin_mol
1558  END IF
1559  END IF
1560 
1561  IF (debug_stress) THEN
1562  stdeb = -1.0_dp*fconv*exc
1563  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1564  'STRESS| VOL 1st eps_XC[n_in]*n_z', one_third_sum_diag(stdeb), det_3x3(stdeb)
1565  END IF
1566 
1567  ELSE ! use_virial
1568 
1569  ! calculate associated hartree potential
1570  ! contribution for both T and D^Z
1571  IF (gapw) THEN
1572  CALL pw_axpy(local_rho_set_t%rho0_mpole%rho0_s_gs, rhoz_tot_gspace)
1573  END IF
1574  CALL pw_poisson_solve(poisson_env, rhoz_tot_gspace, ehartree, zv_hartree_gspace)
1575 
1576  END IF ! use virial
1577  IF (gapw .OR. gapw_xc) THEN
1578  IF (ASSOCIATED(local_rho_set_t)) CALL local_rho_set_release(local_rho_set_t)
1579  END IF
1580 
1581  IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
1582  IF (debug_stress .AND. use_virial) stdeb = virial%pv_ehartree
1583  CALL pw_transfer(zv_hartree_gspace, zv_hartree_rspace)
1584  CALL pw_scale(zv_hartree_rspace, zv_hartree_rspace%pw_grid%dvol)
1585  ! Getting nuclear force contribution from the core charge density (not for GAPW)
1586  CALL integrate_v_core_rspace(zv_hartree_rspace, qs_env)
1587  IF (debug_forces) THEN
1588  fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
1589  CALL para_env%sum(fodeb)
1590  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vh(rhoz)*dncore ", fodeb
1591  END IF
1592  IF (debug_stress .AND. use_virial) THEN
1593  stdeb = fconv*(virial%pv_ehartree - stdeb)
1594  CALL para_env%sum(stdeb)
1595  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1596  'STRESS| INT Vh(rhoz)*dncore ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1597  END IF
1598 
1599  !
1600  IF (gapw_xc) THEN
1601  CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
1602  ELSE
1603  CALL get_qs_env(qs_env=qs_env, rho=rho)
1604  END IF
1605  IF (dft_control%do_admm) THEN
1606  CALL get_qs_env(qs_env, admm_env=admm_env)
1607  xc_section => admm_env%xc_section_primary
1608  ELSE
1609  xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
1610  END IF
1611 
1612  IF (use_virial) THEN
1613  virial%pv_xc = 0.0_dp
1614  END IF
1615  !
1616  NULLIFY (v_xc, v_xc_tau)
1617  IF (gapw_xc) THEN
1618  CALL create_kernel(qs_env, vxc=v_xc, vxc_tau=v_xc_tau, &
1619  rho=rho_xc, rho1_r=rhoz_r_xc, rho1_g=rhoz_g_xc, tau1_r=tauz_r_xc, &
1620  xc_section=xc_section, compute_virial=use_virial, virial_xc=virial%pv_xc)
1621  ELSE
1622  CALL create_kernel(qs_env, vxc=v_xc, vxc_tau=v_xc_tau, &
1623  rho=rho, rho1_r=rhoz_r, rho1_g=rhoz_g, tau1_r=tauz_r, &
1624  xc_section=xc_section, compute_virial=use_virial, virial_xc=virial%pv_xc)
1625  END IF
1626 
1627  IF (gapw .OR. gapw_xc) THEN
1628  !get local_rho_set for GS density and response potential / density
1629  NULLIFY (local_rho_set_t)
1630  CALL local_rho_set_create(local_rho_set_t)
1631  CALL allocate_rho_atom_internals(local_rho_set_t%rho_atom_set, atomic_kind_set, &
1632  qs_kind_set, dft_control, para_env)
1633  CALL init_rho0(local_rho_set_t, qs_env, dft_control%qs_control%gapw_control, &
1634  zcore=0.0_dp)
1635  CALL rho0_s_grid_create(pw_env, local_rho_set_t%rho0_mpole)
1636  CALL calculate_rho_atom_coeff(qs_env, mpa(:), local_rho_set_t%rho_atom_set, &
1637  qs_kind_set, oce, sab_orb, para_env)
1638  CALL prepare_gapw_den(qs_env, local_rho_set_t, do_rho0=gapw)
1639  NULLIFY (local_rho_set_gs)
1640  CALL local_rho_set_create(local_rho_set_gs)
1641  CALL allocate_rho_atom_internals(local_rho_set_gs%rho_atom_set, atomic_kind_set, &
1642  qs_kind_set, dft_control, para_env)
1643  CALL init_rho0(local_rho_set_gs, qs_env, dft_control%qs_control%gapw_control)
1644  CALL rho0_s_grid_create(pw_env, local_rho_set_gs%rho0_mpole)
1645  CALL calculate_rho_atom_coeff(qs_env, matrix_p(:, 1), local_rho_set_gs%rho_atom_set, &
1646  qs_kind_set, oce, sab_orb, para_env)
1647  CALL prepare_gapw_den(qs_env, local_rho_set_gs, do_rho0=gapw)
1648  ! compute response potential
1649  ALLOCATE (rho_r_t(nspins), rho_g_t(nspins))
1650  DO ispin = 1, nspins
1651  CALL auxbas_pw_pool%create_pw(rho_r_t(ispin))
1652  CALL auxbas_pw_pool%create_pw(rho_g_t(ispin))
1653  END DO
1654  CALL auxbas_pw_pool%create_pw(rho_tot_gspace_t)
1655  total_rho_t = 0.0_dp
1656  CALL pw_zero(rho_tot_gspace_t)
1657  DO ispin = 1, nspins
1658  CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1659  rho=rho_r_t(ispin), &
1660  rho_gspace=rho_g_t(ispin), &
1661  soft_valid=gapw, &
1662  total_rho=total_rho_t(ispin))
1663  CALL pw_axpy(rho_g_t(ispin), rho_tot_gspace_t)
1664  END DO
1665  ! add rho0 contributions to response density (only for Coulomb) only for gapw
1666  IF (gapw) THEN
1667  CALL pw_axpy(local_rho_set_t%rho0_mpole%rho0_s_gs, rho_tot_gspace_t)
1668  ! compute response Coulomb potential
1669  CALL auxbas_pw_pool%create_pw(v_hartree_gspace_t)
1670  CALL auxbas_pw_pool%create_pw(v_hartree_rspace_t)
1671  NULLIFY (hartree_local_t)
1672  CALL hartree_local_create(hartree_local_t)
1673  CALL init_coulomb_local(hartree_local_t, natom)
1674  CALL pw_poisson_solve(poisson_env, rho_tot_gspace_t, hartree_t, v_hartree_gspace_t)
1675  CALL pw_transfer(v_hartree_gspace_t, v_hartree_rspace_t)
1676  CALL pw_scale(v_hartree_rspace_t, v_hartree_rspace_t%pw_grid%dvol)
1677  !
1678  IF (debug_forces) fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1)
1679  CALL vh_1c_gg_integrals(qs_env, hartree_t, hartree_local_t%ecoul_1c, local_rho_set_gs, para_env, tddft=.false., &
1680  local_rho_set_2nd=local_rho_set_t, core_2nd=.true.) ! n^core for GS potential
1681  CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace_t, para_env, calculate_forces=.true., &
1682  local_rho_set=local_rho_set_gs, local_rho_set_2nd=local_rho_set_t)
1683  IF (debug_forces) THEN
1684  fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1) - fodeb(1:3)
1685  CALL para_env%sum(fodeb)
1686  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vh(T)*dncore PAWg0", fodeb
1687  END IF
1688  END IF !gapw
1689  END IF !gapw
1690 
1691  IF (gapw .OR. gapw_xc) THEN
1692  !GAPW compute atomic fxc contributions
1693  IF (myfun /= xc_none) THEN
1694  ! local_rho_set_f
1695  NULLIFY (local_rho_set_f)
1696  CALL local_rho_set_create(local_rho_set_f)
1697  CALL allocate_rho_atom_internals(local_rho_set_f%rho_atom_set, atomic_kind_set, &
1698  qs_kind_set, dft_control, para_env)
1699  CALL calculate_rho_atom_coeff(qs_env, mpa, local_rho_set_f%rho_atom_set, &
1700  qs_kind_set, oce, sab_orb, para_env)
1701  CALL prepare_gapw_den(qs_env, local_rho_set_f, do_rho0=.false.)
1702  ! add hard and soft atomic contributions
1703  CALL calculate_xc_2nd_deriv_atom(local_rho_set_gs%rho_atom_set, &
1704  local_rho_set_f%rho_atom_set, &
1705  qs_env, xc_section, para_env, &
1706  do_tddft=.false., do_triplet=.false.)
1707  END IF ! myfun
1708  END IF
1709 
1710  ! Stress-tensor XC-kernel GGA contribution
1711  IF (use_virial) THEN
1712  virial%pv_exc = virial%pv_exc + virial%pv_xc
1713  virial%pv_virial = virial%pv_virial + virial%pv_xc
1714  END IF
1715 
1716  IF (debug_stress .AND. use_virial) THEN
1717  stdeb = 1.0_dp*fconv*virial%pv_xc
1718  CALL para_env%sum(stdeb)
1719  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1720  'STRESS| GGA 2nd Pin*dK*rhoz', one_third_sum_diag(stdeb), det_3x3(stdeb)
1721  END IF
1722 
1723  ! Stress-tensor integral contribution of 2nd derivative terms
1724  IF (use_virial) THEN
1725  pv_loc = virial%pv_virial
1726  END IF
1727 
1728  CALL get_qs_env(qs_env=qs_env, rho=rho)
1729  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1730  IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1731 
1732  DO ispin = 1, nspins
1733  CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
1734  END DO
1735  IF ((.NOT. (gapw)) .AND. (.NOT. gapw_xc)) THEN
1736  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1737  DO ispin = 1, nspins
1738  CALL pw_axpy(zv_hartree_rspace, v_xc(ispin)) ! Hartree potential of response density
1739  CALL integrate_v_rspace(qs_env=qs_env, &
1740  v_rspace=v_xc(ispin), &
1741  hmat=matrix_hz(ispin), &
1742  pmat=matrix_p(ispin, 1), &
1743  gapw=.false., &
1744  calculate_forces=.true.)
1745  END DO
1746  IF (debug_forces) THEN
1747  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1748  CALL para_env%sum(fodeb)
1749  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKhxc*rhoz ", fodeb
1750  END IF
1751  ELSE
1752  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1753  IF (myfun /= xc_none) THEN
1754  DO ispin = 1, nspins
1755  CALL integrate_v_rspace(qs_env=qs_env, &
1756  v_rspace=v_xc(ispin), &
1757  hmat=matrix_hz(ispin), &
1758  pmat=matrix_p(ispin, 1), &
1759  gapw=.true., &
1760  calculate_forces=.true.)
1761  END DO
1762  END IF ! my_fun
1763  ! Coulomb T+Dz
1764  DO ispin = 1, nspins
1765  CALL pw_zero(v_xc(ispin))
1766  IF (gapw) THEN ! Hartree potential of response density
1767  CALL pw_axpy(v_hartree_rspace_t, v_xc(ispin))
1768  ELSEIF (gapw_xc) THEN
1769  CALL pw_axpy(zv_hartree_rspace, v_xc(ispin))
1770  END IF
1771  CALL integrate_v_rspace(qs_env=qs_env, &
1772  v_rspace=v_xc(ispin), &
1773  hmat=matrix_ht(ispin), &
1774  pmat=matrix_p(ispin, 1), &
1775  gapw=gapw, &
1776  calculate_forces=.true.)
1777  END DO
1778  IF (debug_forces) THEN
1779  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1780  CALL para_env%sum(fodeb)
1781  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKhxc*rhoz ", fodeb
1782  END IF
1783  END IF
1784 
1785  IF (gapw .OR. gapw_xc) THEN
1786  ! compute hard and soft atomic contributions
1787  IF (myfun /= xc_none) THEN
1788  IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
1789  CALL update_ks_atom(qs_env, matrix_hz, matrix_p, forces=.true., tddft=.false., &
1790  rho_atom_external=local_rho_set_f%rho_atom_set)
1791  IF (debug_forces) THEN
1792  fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
1793  CALL para_env%sum(fodeb)
1794  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P^GS*dKxc*(Dz+T) PAW", fodeb
1795  END IF
1796  END IF !myfun
1797  ! Coulomb contributions
1798  IF (gapw) THEN
1799  IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
1800  CALL update_ks_atom(qs_env, matrix_ht, matrix_p, forces=.true., tddft=.false., &
1801  rho_atom_external=local_rho_set_t%rho_atom_set)
1802  IF (debug_forces) THEN
1803  fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
1804  CALL para_env%sum(fodeb)
1805  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P^GS*dKh*(Dz+T) PAW", fodeb
1806  END IF
1807  END IF
1808  ! add Coulomb and XC
1809  DO ispin = 1, nspins
1810  CALL dbcsr_add(matrix_hz(ispin)%matrix, matrix_ht(ispin)%matrix, 1.0_dp, 1.0_dp)
1811  END DO
1812 
1813  ! release
1814  IF (myfun /= xc_none) THEN
1815  IF (ASSOCIATED(local_rho_set_f)) CALL local_rho_set_release(local_rho_set_f)
1816  END IF
1817  IF (ASSOCIATED(local_rho_set_t)) CALL local_rho_set_release(local_rho_set_t)
1818  IF (ASSOCIATED(local_rho_set_gs)) CALL local_rho_set_release(local_rho_set_gs)
1819  IF (gapw) THEN
1820  IF (ASSOCIATED(hartree_local_t)) CALL hartree_local_release(hartree_local_t)
1821  CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace_t)
1822  CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace_t)
1823  END IF
1824  CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace_t)
1825  DO ispin = 1, nspins
1826  CALL auxbas_pw_pool%give_back_pw(rho_r_t(ispin))
1827  CALL auxbas_pw_pool%give_back_pw(rho_g_t(ispin))
1828  END DO
1829  DEALLOCATE (rho_r_t, rho_g_t)
1830  END IF ! gapw
1831 
1832  IF (debug_stress .AND. use_virial) THEN
1833  stdeb = fconv*(virial%pv_virial - stdeb)
1834  CALL para_env%sum(stdeb)
1835  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1836  'STRESS| INT 2nd f_Hxc[Pz]*Pin', one_third_sum_diag(stdeb), det_3x3(stdeb)
1837  END IF
1838  !
1839  IF (ASSOCIATED(v_xc_tau)) THEN
1840  cpassert(.NOT. (gapw .OR. gapw_xc))
1841  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1842  IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1843  DO ispin = 1, nspins
1844  CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
1845  CALL integrate_v_rspace(qs_env=qs_env, &
1846  v_rspace=v_xc_tau(ispin), &
1847  hmat=matrix_hz(ispin), &
1848  pmat=matrix_p(ispin, 1), &
1849  compute_tau=.true., &
1850  gapw=(gapw .OR. gapw_xc), &
1851  calculate_forces=.true.)
1852  END DO
1853  IF (debug_forces) THEN
1854  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1855  CALL para_env%sum(fodeb)
1856  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKtau*tauz ", fodeb
1857  END IF
1858  END IF
1859  IF (debug_stress .AND. use_virial) THEN
1860  stdeb = fconv*(virial%pv_virial - stdeb)
1861  CALL para_env%sum(stdeb)
1862  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1863  'STRESS| INT 2nd f_xctau[Pz]*Pin', one_third_sum_diag(stdeb), det_3x3(stdeb)
1864  END IF
1865  ! Stress-tensor integral contribution of 2nd derivative terms
1866  IF (use_virial) THEN
1867  virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1868  END IF
1869 
1870  ! KG Embedding
1871  ! calculate kinetic energy kernel, folded with response density for partial integration
1872  IF (dft_control%qs_control%do_kg) THEN
1873  IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed) THEN
1874  ekin_mol = 0.0_dp
1875  IF (use_virial) THEN
1876  pv_loc = virial%pv_virial
1877  END IF
1878 
1879  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1880  IF (use_virial) virial%pv_xc = 0.0_dp
1881  CALL kg_ekin_subset(qs_env=qs_env, &
1882  ks_matrix=matrix_hz, &
1883  ekin_mol=ekin_mol, &
1884  calc_force=.true., &
1885  do_kernel=.true., &
1886  pmat_ext=matrix_pz)
1887 
1888  IF (debug_forces) THEN
1889  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1890  CALL para_env%sum(fodeb)
1891  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*d(Kkg)*rhoz ", fodeb
1892  END IF
1893  IF (debug_stress .AND. use_virial) THEN
1894  stdeb = fconv*(virial%pv_virial - pv_loc)
1895  CALL para_env%sum(stdeb)
1896  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1897  'STRESS| INT KG Pin*d(KKG)*rhoz ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1898 
1899  stdeb = fconv*(virial%pv_xc)
1900  CALL para_env%sum(stdeb)
1901  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1902  'STRESS| GGA KG Pin*d(KKG)*rhoz ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1903  END IF
1904 
1905  ! Stress tensor
1906  IF (use_virial) THEN
1907  ! XC-kernel Integral contribution
1908  virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1909 
1910  ! XC-kernel GGA contribution
1911  virial%pv_exc = virial%pv_exc - virial%pv_xc
1912  virial%pv_virial = virial%pv_virial - virial%pv_xc
1913  virial%pv_xc = 0.0_dp
1914  END IF
1915  END IF
1916  END IF
1917  CALL auxbas_pw_pool%give_back_pw(rhoz_tot_gspace)
1918  CALL auxbas_pw_pool%give_back_pw(zv_hartree_gspace)
1919  CALL auxbas_pw_pool%give_back_pw(zv_hartree_rspace)
1920  DO ispin = 1, nspins
1921  CALL auxbas_pw_pool%give_back_pw(rhoz_r(ispin))
1922  CALL auxbas_pw_pool%give_back_pw(rhoz_g(ispin))
1923  CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
1924  END DO
1925  DEALLOCATE (rhoz_r, rhoz_g, v_xc)
1926  IF (gapw_xc) THEN
1927  DO ispin = 1, nspins
1928  CALL auxbas_pw_pool%give_back_pw(rhoz_r_xc(ispin))
1929  CALL auxbas_pw_pool%give_back_pw(rhoz_g_xc(ispin))
1930  END DO
1931  DEALLOCATE (rhoz_r_xc, rhoz_g_xc)
1932  END IF
1933  IF (ASSOCIATED(v_xc_tau)) THEN
1934  DO ispin = 1, nspins
1935  CALL auxbas_pw_pool%give_back_pw(tauz_r(ispin))
1936  CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
1937  END DO
1938  DEALLOCATE (tauz_r, v_xc_tau)
1939  END IF
1940  IF (debug_forces) THEN
1941  ALLOCATE (ftot3(3, natom))
1942  CALL total_qs_force(ftot3, force, atomic_kind_set)
1943  fodeb(1:3) = ftot3(1:3, 1) - ftot2(1:3, 1)
1944  CALL para_env%sum(fodeb)
1945  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*V(rhoz)", fodeb
1946  END IF
1947  CALL dbcsr_deallocate_matrix_set(scrm)
1948  CALL dbcsr_deallocate_matrix_set(matrix_ht)
1949 
1950  ! -----------------------------------------
1951  ! Apply ADMM exchange correction
1952  ! -----------------------------------------
1953 
1954  IF (dft_control%do_admm) THEN
1955  ! volume term
1956  exc_aux_fit = 0.0_dp
1957 
1958  IF (qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
1959  ! nothing to do
1960  NULLIFY (mpz, mhz, mhx, mhy)
1961  ELSE
1962  ! add ADMM xc_section_aux terms: Pz*Vxc + P0*K0[rhoz]
1963  CALL get_qs_env(qs_env, admm_env=admm_env)
1964  CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit, matrix_s_aux_fit=scrm, &
1965  task_list_aux_fit=task_list_aux_fit)
1966  !
1967  NULLIFY (mpz, mhz, mhx, mhy)
1968  CALL dbcsr_allocate_matrix_set(mhx, nspins, 1)
1969  CALL dbcsr_allocate_matrix_set(mhy, nspins, 1)
1970  CALL dbcsr_allocate_matrix_set(mpz, nspins, 1)
1971  DO ispin = 1, nspins
1972  ALLOCATE (mhx(ispin, 1)%matrix)
1973  CALL dbcsr_create(mhx(ispin, 1)%matrix, template=scrm(1)%matrix)
1974  CALL dbcsr_copy(mhx(ispin, 1)%matrix, scrm(1)%matrix)
1975  CALL dbcsr_set(mhx(ispin, 1)%matrix, 0.0_dp)
1976  ALLOCATE (mhy(ispin, 1)%matrix)
1977  CALL dbcsr_create(mhy(ispin, 1)%matrix, template=scrm(1)%matrix)
1978  CALL dbcsr_copy(mhy(ispin, 1)%matrix, scrm(1)%matrix)
1979  CALL dbcsr_set(mhy(ispin, 1)%matrix, 0.0_dp)
1980  ALLOCATE (mpz(ispin, 1)%matrix)
1981  IF (do_ex) THEN
1982  CALL dbcsr_create(mpz(ispin, 1)%matrix, template=p_env%p1_admm(ispin)%matrix)
1983  CALL dbcsr_copy(mpz(ispin, 1)%matrix, p_env%p1_admm(ispin)%matrix)
1984  CALL dbcsr_add(mpz(ispin, 1)%matrix, ex_env%matrix_pe_admm(ispin)%matrix, &
1985  1.0_dp, 1.0_dp)
1986  ELSE
1987  CALL dbcsr_create(mpz(ispin, 1)%matrix, template=matrix_pz_admm(ispin)%matrix)
1988  CALL dbcsr_copy(mpz(ispin, 1)%matrix, matrix_pz_admm(ispin)%matrix)
1989  END IF
1990  END DO
1991  !
1992  xc_section => admm_env%xc_section_aux
1993  ! Stress-tensor: integration contribution direct term
1994  ! int Pz*v_xc[rho_admm]
1995  IF (use_virial) THEN
1996  pv_loc = virial%pv_virial
1997  END IF
1998 
1999  basis_type = "AUX_FIT"
2000  task_list => task_list_aux_fit
2001  IF (admm_env%do_gapw) THEN
2002  basis_type = "AUX_FIT_SOFT"
2003  task_list => admm_env%admm_gapw_env%task_list
2004  END IF
2005  !
2006  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2007  IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2008  DO ispin = 1, nspins
2009  CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), &
2010  hmat=mhx(ispin, 1), pmat=mpz(ispin, 1), &
2011  qs_env=qs_env, calculate_forces=.true., &
2012  basis_type=basis_type, task_list_external=task_list)
2013  END DO
2014  IF (debug_forces) THEN
2015  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2016  CALL para_env%sum(fodeb)
2017  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*Vxc(rho_admm)", fodeb
2018  END IF
2019  IF (debug_stress .AND. use_virial) THEN
2020  stdeb = fconv*(virial%pv_virial - pv_loc)
2021  CALL para_env%sum(stdeb)
2022  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2023  'STRESS| INT 1st Pz*dVxc(rho_admm) ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2024  END IF
2025  ! Stress-tensor Pz_admm*v_xc[rho_admm]
2026  IF (use_virial) THEN
2027  virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2028  END IF
2029  !
2030  IF (admm_env%do_gapw) THEN
2031  CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
2032  IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
2033  CALL update_ks_atom(qs_env, mhx(:, 1), mpz(:, 1), forces=.true., tddft=.false., &
2034  rho_atom_external=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
2035  kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
2036  oce_external=admm_env%admm_gapw_env%oce, &
2037  sab_external=sab_aux_fit)
2038  IF (debug_forces) THEN
2039  fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
2040  CALL para_env%sum(fodeb)
2041  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*Vxc(rho_admm)PAW", fodeb
2042  END IF
2043  END IF
2044  !
2045  NULLIFY (rho_g_aux, rho_r_aux, tau_r_aux)
2046  CALL qs_rho_get(rho_aux_fit, rho_r=rho_r_aux, rho_g=rho_g_aux, tau_r=tau_r_aux)
2047  ! rhoz_aux
2048  NULLIFY (rhoz_g_aux, rhoz_r_aux)
2049  ALLOCATE (rhoz_r_aux(nspins), rhoz_g_aux(nspins))
2050  DO ispin = 1, nspins
2051  CALL auxbas_pw_pool%create_pw(rhoz_r_aux(ispin))
2052  CALL auxbas_pw_pool%create_pw(rhoz_g_aux(ispin))
2053  END DO
2054  DO ispin = 1, nspins
2055  CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpz(ispin, 1)%matrix, &
2056  rho=rhoz_r_aux(ispin), rho_gspace=rhoz_g_aux(ispin), &
2057  basis_type=basis_type, task_list_external=task_list)
2058  END DO
2059  !
2060  ! Add ADMM volume contribution to stress tensor
2061  IF (use_virial) THEN
2062 
2063  ! Stress tensor volume term: \int v_xc[n_in_admm]*n_z_admm
2064  ! vadmm_rspace already scaled, we need to unscale it!
2065  DO ispin = 1, nspins
2066  exc_aux_fit = exc_aux_fit + pw_integral_ab(rhoz_r_aux(ispin), vadmm_rspace(ispin))/ &
2067  vadmm_rspace(ispin)%pw_grid%dvol
2068  END DO
2069 
2070  IF (debug_stress) THEN
2071  stdeb = -1.0_dp*fconv*exc_aux_fit
2072  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T43,2(1X,ES19.11))") &
2073  'STRESS| VOL 1st eps_XC[n_in_admm]*n_z_admm', one_third_sum_diag(stdeb), det_3x3(stdeb)
2074  END IF
2075 
2076  END IF
2077  !
2078  NULLIFY (v_xc)
2079 
2080  IF (use_virial) virial%pv_xc = 0.0_dp
2081 
2082  CALL create_kernel(qs_env=qs_env, &
2083  vxc=v_xc, &
2084  vxc_tau=v_xc_tau, &
2085  rho=rho_aux_fit, &
2086  rho1_r=rhoz_r_aux, &
2087  rho1_g=rhoz_g_aux, &
2088  tau1_r=tau_r_aux, &
2089  xc_section=xc_section, &
2090  compute_virial=use_virial, &
2091  virial_xc=virial%pv_xc)
2092 
2093  ! Stress-tensor ADMM-kernel GGA contribution
2094  IF (use_virial) THEN
2095  virial%pv_exc = virial%pv_exc + virial%pv_xc
2096  virial%pv_virial = virial%pv_virial + virial%pv_xc
2097  END IF
2098 
2099  IF (debug_stress .AND. use_virial) THEN
2100  stdeb = 1.0_dp*fconv*virial%pv_xc
2101  CALL para_env%sum(stdeb)
2102  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2103  'STRESS| GGA 2nd Pin_admm*dK*rhoz_admm', one_third_sum_diag(stdeb), det_3x3(stdeb)
2104  END IF
2105  !
2106  CALL qs_rho_get(rho_aux_fit, rho_ao_kp=matrix_p)
2107  ! Stress-tensor Pin*dK*rhoz_admm
2108  IF (use_virial) THEN
2109  virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2110  END IF
2111  IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2112  IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2113  DO ispin = 1, nspins
2114  CALL dbcsr_set(mhy(ispin, 1)%matrix, 0.0_dp)
2115  CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2116  CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
2117  hmat=mhy(ispin, 1), pmat=matrix_p(ispin, 1), &
2118  calculate_forces=.true., &
2119  basis_type=basis_type, task_list_external=task_list)
2120  END DO
2121  IF (debug_forces) THEN
2122  fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2123  CALL para_env%sum(fodeb)
2124  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dK*rhoz_admm ", fodeb
2125  END IF
2126  IF (debug_stress .AND. use_virial) THEN
2127  stdeb = fconv*(virial%pv_virial - pv_loc)
2128  CALL para_env%sum(stdeb)
2129  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2130  'STRESS| INT 2nd Pin*dK*rhoz_admm ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2131  END IF
2132  ! Stress-tensor Pin*dK*rhoz_admm
2133  IF (use_virial) THEN
2134  virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2135  END IF
2136  DO ispin = 1, nspins
2137  CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2138  CALL auxbas_pw_pool%give_back_pw(rhoz_r_aux(ispin))
2139  CALL auxbas_pw_pool%give_back_pw(rhoz_g_aux(ispin))
2140  END DO
2141  DEALLOCATE (v_xc, rhoz_r_aux, rhoz_g_aux)
2142  !
2143  IF (admm_env%do_gapw) THEN
2144  CALL local_rho_set_create(local_rhoz_set_admm)
2145  CALL allocate_rho_atom_internals(local_rhoz_set_admm%rho_atom_set, atomic_kind_set, &
2146  admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
2147  CALL calculate_rho_atom_coeff(qs_env, mpz(:, 1), local_rhoz_set_admm%rho_atom_set, &
2148  admm_env%admm_gapw_env%admm_kind_set, &
2149  admm_env%admm_gapw_env%oce, sab_aux_fit, para_env)
2150  CALL prepare_gapw_den(qs_env, local_rho_set=local_rhoz_set_admm, &
2151  do_rho0=.false., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
2152  !compute the potential due to atomic densities
2153  CALL calculate_xc_2nd_deriv_atom(admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
2154  local_rhoz_set_admm%rho_atom_set, &
2155  qs_env, xc_section, para_env, do_tddft=.false., do_triplet=.false., &
2156  kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
2157  !
2158  IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
2159  CALL update_ks_atom(qs_env, mhy(:, 1), matrix_p(:, 1), forces=.true., tddft=.false., &
2160  rho_atom_external=local_rhoz_set_admm%rho_atom_set, &
2161  kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
2162  oce_external=admm_env%admm_gapw_env%oce, &
2163  sab_external=sab_aux_fit)
2164  IF (debug_forces) THEN
2165  fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
2166  CALL para_env%sum(fodeb)
2167  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dK*rhoz_admm[PAW] ", fodeb
2168  END IF
2169  CALL local_rho_set_release(local_rhoz_set_admm)
2170  END IF
2171  !
2172  nao = admm_env%nao_orb
2173  nao_aux = admm_env%nao_aux_fit
2174  ALLOCATE (dbwork)
2175  CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
2176  DO ispin = 1, nspins
2177  CALL cp_dbcsr_sm_fm_multiply(mhy(ispin, 1)%matrix, admm_env%A, &
2178  admm_env%work_aux_orb, nao)
2179  CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
2180  1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
2181  admm_env%work_orb_orb)
2182  CALL dbcsr_copy(dbwork, matrix_hz(ispin)%matrix)
2183  CALL dbcsr_set(dbwork, 0.0_dp)
2184  CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.true.)
2185  CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
2186  END DO
2187  CALL dbcsr_release(dbwork)
2188  DEALLOCATE (dbwork)
2189  CALL dbcsr_deallocate_matrix_set(mpz)
2190  END IF ! qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none
2191  END IF ! do_admm
2192 
2193  ! -----------------------------------------
2194  ! HFX
2195  ! -----------------------------------------
2196 
2197  ! HFX
2198  hfx_section => section_vals_get_subs_vals(xc_section, "HF")
2199  CALL section_vals_get(hfx_section, explicit=do_hfx)
2200  IF (do_hfx) THEN
2201  CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
2202  cpassert(n_rep_hf == 1)
2203  CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
2204  i_rep_section=1)
2205  mspin = 1
2206  IF (hfx_treat_lsd_in_core) mspin = nspins
2207  IF (use_virial) virial%pv_fock_4c = 0.0_dp
2208  !
2209  CALL get_qs_env(qs_env=qs_env, rho=rho, x_data=x_data, &
2210  s_mstruct_changed=s_mstruct_changed)
2211  distribute_fock_matrix = .true.
2212 
2213  ! -----------------------------------------
2214  ! HFX-ADMM
2215  ! -----------------------------------------
2216  IF (dft_control%do_admm) THEN
2217  CALL get_qs_env(qs_env=qs_env, admm_env=admm_env)
2218  CALL get_admm_env(admm_env, matrix_s_aux_fit=scrm, rho_aux_fit=rho_aux_fit)
2219  CALL qs_rho_get(rho_aux_fit, rho_ao_kp=matrix_p)
2220  NULLIFY (mpz, mhz, mpd, mhd)
2221  CALL dbcsr_allocate_matrix_set(mpz, nspins, 1)
2222  CALL dbcsr_allocate_matrix_set(mhz, nspins, 1)
2223  CALL dbcsr_allocate_matrix_set(mpd, nspins, 1)
2224  CALL dbcsr_allocate_matrix_set(mhd, nspins, 1)
2225  DO ispin = 1, nspins
2226  ALLOCATE (mhz(ispin, 1)%matrix, mhd(ispin, 1)%matrix)
2227  CALL dbcsr_create(mhz(ispin, 1)%matrix, template=scrm(1)%matrix)
2228  CALL dbcsr_create(mhd(ispin, 1)%matrix, template=scrm(1)%matrix)
2229  CALL dbcsr_copy(mhz(ispin, 1)%matrix, scrm(1)%matrix)
2230  CALL dbcsr_copy(mhd(ispin, 1)%matrix, scrm(1)%matrix)
2231  CALL dbcsr_set(mhz(ispin, 1)%matrix, 0.0_dp)
2232  CALL dbcsr_set(mhd(ispin, 1)%matrix, 0.0_dp)
2233  ALLOCATE (mpz(ispin, 1)%matrix)
2234  IF (do_ex) THEN
2235  CALL dbcsr_create(mpz(ispin, 1)%matrix, template=scrm(1)%matrix)
2236  CALL dbcsr_copy(mpz(ispin, 1)%matrix, p_env%p1_admm(ispin)%matrix)
2237  CALL dbcsr_add(mpz(ispin, 1)%matrix, ex_env%matrix_pe_admm(ispin)%matrix, &
2238  1.0_dp, 1.0_dp)
2239  ELSE
2240  CALL dbcsr_create(mpz(ispin, 1)%matrix, template=scrm(1)%matrix)
2241  CALL dbcsr_copy(mpz(ispin, 1)%matrix, matrix_pz_admm(ispin)%matrix)
2242  END IF
2243  mpd(ispin, 1)%matrix => matrix_p(ispin, 1)%matrix
2244  END DO
2245  !
2246  IF (x_data(1, 1)%do_hfx_ri) THEN
2247 
2248  eh1 = 0.0_dp
2249  CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpz, &
2250  geometry_did_change=s_mstruct_changed, nspins=nspins, &
2251  hf_fraction=x_data(1, 1)%general_parameter%fraction)
2252 
2253  eh1 = 0.0_dp
2254  CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhd, eh1, rho_ao=mpd, &
2255  geometry_did_change=s_mstruct_changed, nspins=nspins, &
2256  hf_fraction=x_data(1, 1)%general_parameter%fraction)
2257 
2258  ELSE
2259  DO ispin = 1, mspin
2260  eh1 = 0.0
2261  CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpz, hfx_section, &
2262  para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
2263  ispin=ispin)
2264  END DO
2265  DO ispin = 1, mspin
2266  eh1 = 0.0
2267  CALL integrate_four_center(qs_env, x_data, mhd, eh1, mpd, hfx_section, &
2268  para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
2269  ispin=ispin)
2270  END DO
2271  END IF
2272  !
2273  CALL get_qs_env(qs_env, admm_env=admm_env)
2274  cpassert(ASSOCIATED(admm_env%work_aux_orb))
2275  cpassert(ASSOCIATED(admm_env%work_orb_orb))
2276  nao = admm_env%nao_orb
2277  nao_aux = admm_env%nao_aux_fit
2278  ALLOCATE (dbwork)
2279  CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
2280  DO ispin = 1, nspins
2281  CALL cp_dbcsr_sm_fm_multiply(mhz(ispin, 1)%matrix, admm_env%A, &
2282  admm_env%work_aux_orb, nao)
2283  CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
2284  1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
2285  admm_env%work_orb_orb)
2286  CALL dbcsr_copy(dbwork, matrix_hz(ispin)%matrix)
2287  CALL dbcsr_set(dbwork, 0.0_dp)
2288  CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.true.)
2289  CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
2290  END DO
2291  CALL dbcsr_release(dbwork)
2292  DEALLOCATE (dbwork)
2293  ! derivatives Tr (Pz [A(T)H dA/dR])
2294  IF (debug_forces) fodeb(1:3) = force(1)%overlap_admm(1:3, 1)
2295  IF (ASSOCIATED(mhx) .AND. ASSOCIATED(mhy)) THEN
2296  DO ispin = 1, nspins
2297  CALL dbcsr_add(mhd(ispin, 1)%matrix, mhx(ispin, 1)%matrix, 1.0_dp, 1.0_dp)
2298  CALL dbcsr_add(mhz(ispin, 1)%matrix, mhy(ispin, 1)%matrix, 1.0_dp, 1.0_dp)
2299  END DO
2300  END IF
2301  CALL qs_rho_get(rho, rho_ao=matrix_pd)
2302  CALL admm_projection_derivative(qs_env, mhd(:, 1), mpa)
2303  CALL admm_projection_derivative(qs_env, mhz(:, 1), matrix_pd)
2304  IF (debug_forces) THEN
2305  fodeb(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb(1:3)
2306  CALL para_env%sum(fodeb)
2307  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*hfx*S' ", fodeb
2308  END IF
2309  CALL dbcsr_deallocate_matrix_set(mpz)
2310  CALL dbcsr_deallocate_matrix_set(mhz)
2311  CALL dbcsr_deallocate_matrix_set(mhd)
2312  IF (ASSOCIATED(mhx) .AND. ASSOCIATED(mhy)) THEN
2313  CALL dbcsr_deallocate_matrix_set(mhx)
2314  CALL dbcsr_deallocate_matrix_set(mhy)
2315  END IF
2316  DEALLOCATE (mpd)
2317  ELSE
2318  ! -----------------------------------------
2319  ! conventional HFX
2320  ! -----------------------------------------
2321  ALLOCATE (mpz(nspins, 1), mhz(nspins, 1))
2322  DO ispin = 1, nspins
2323  mhz(ispin, 1)%matrix => matrix_hz(ispin)%matrix
2324  mpz(ispin, 1)%matrix => mpa(ispin)%matrix
2325  END DO
2326 
2327  IF (x_data(1, 1)%do_hfx_ri) THEN
2328 
2329  eh1 = 0.0_dp
2330  CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpz, &
2331  geometry_did_change=s_mstruct_changed, nspins=nspins, &
2332  hf_fraction=x_data(1, 1)%general_parameter%fraction)
2333  ELSE
2334  DO ispin = 1, mspin
2335  eh1 = 0.0
2336  CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpz, hfx_section, &
2337  para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
2338  ispin=ispin)
2339  END DO
2340  END IF
2341  DEALLOCATE (mhz, mpz)
2342  END IF
2343 
2344  ! -----------------------------------------
2345  ! HFX FORCES
2346  ! -----------------------------------------
2347 
2348  resp_only = .true.
2349  IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
2350  IF (dft_control%do_admm) THEN
2351  ! -----------------------------------------
2352  ! HFX-ADMM FORCES
2353  ! -----------------------------------------
2354  CALL qs_rho_get(rho_aux_fit, rho_ao_kp=matrix_p)
2355  NULLIFY (matrix_pza)
2356  CALL dbcsr_allocate_matrix_set(matrix_pza, nspins)
2357  DO ispin = 1, nspins
2358  ALLOCATE (matrix_pza(ispin)%matrix)
2359  IF (do_ex) THEN
2360  CALL dbcsr_create(matrix_pza(ispin)%matrix, template=p_env%p1_admm(ispin)%matrix)
2361  CALL dbcsr_copy(matrix_pza(ispin)%matrix, p_env%p1_admm(ispin)%matrix)
2362  CALL dbcsr_add(matrix_pza(ispin)%matrix, ex_env%matrix_pe_admm(ispin)%matrix, &
2363  1.0_dp, 1.0_dp)
2364  ELSE
2365  CALL dbcsr_create(matrix_pza(ispin)%matrix, template=matrix_pz_admm(ispin)%matrix)
2366  CALL dbcsr_copy(matrix_pza(ispin)%matrix, matrix_pz_admm(ispin)%matrix)
2367  END IF
2368  END DO
2369  IF (x_data(1, 1)%do_hfx_ri) THEN
2370 
2371  CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
2372  x_data(1, 1)%general_parameter%fraction, &
2373  rho_ao=matrix_p, rho_ao_resp=matrix_pza, &
2374  use_virial=use_virial, resp_only=resp_only)
2375  ELSE
2376  CALL derivatives_four_center(qs_env, matrix_p, matrix_pza, hfx_section, para_env, &
2377  1, use_virial, resp_only=resp_only)
2378  END IF
2379  CALL dbcsr_deallocate_matrix_set(matrix_pza)
2380  ELSE
2381  ! -----------------------------------------
2382  ! conventional HFX FORCES
2383  ! -----------------------------------------
2384  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2385  IF (x_data(1, 1)%do_hfx_ri) THEN
2386 
2387  CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
2388  x_data(1, 1)%general_parameter%fraction, &
2389  rho_ao=matrix_p, rho_ao_resp=mpa, &
2390  use_virial=use_virial, resp_only=resp_only)
2391  ELSE
2392  CALL derivatives_four_center(qs_env, matrix_p, mpa, hfx_section, para_env, &
2393  1, use_virial, resp_only=resp_only)
2394  END IF
2395  END IF ! do_admm
2396 
2397  IF (use_virial) THEN
2398  virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
2399  virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
2400  virial%pv_calculate = .false.
2401  END IF
2402 
2403  IF (debug_forces) THEN
2404  fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
2405  CALL para_env%sum(fodeb)
2406  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*hfx ", fodeb
2407  END IF
2408  IF (debug_stress .AND. use_virial) THEN
2409  stdeb = -1.0_dp*fconv*virial%pv_fock_4c
2410  CALL para_env%sum(stdeb)
2411  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2412  'STRESS| Pz*hfx ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2413  END IF
2414  END IF ! do_hfx
2415 
2416  ! Stress-tensor volume contributions
2417  ! These need to be applied at the end of qs_force
2418  IF (use_virial) THEN
2419  ! Adding mixed Hartree energy twice, due to symmetry
2420  zehartree = zehartree + 2.0_dp*ehartree
2421  zexc = zexc + exc
2422  ! ADMM contribution handled differently in qs_force
2423  IF (dft_control%do_admm) THEN
2424  zexc_aux_fit = zexc_aux_fit + exc_aux_fit
2425  END IF
2426  END IF
2427 
2428  ! Overlap matrix
2429  ! H(drho+dz) + Wz
2430  ! If ground-state density matrix solved by diagonalization, then use this
2431  IF (dft_control%qs_control%do_ls_scf) THEN
2432  ! Ground-state density has been calculated by LS
2433  eps_filter = dft_control%qs_control%eps_filter_matrix
2434  CALL calculate_whz_ao_matrix(qs_env, matrix_hz, matrix_wz, eps_filter)
2435  ELSE
2436  IF (do_ex) THEN
2437  matrix_wz => p_env%w1
2438  END IF
2439  focc = 1.0_dp
2440  IF (nspins == 1) focc = 2.0_dp
2441  CALL get_qs_env(qs_env, mos=mos)
2442  DO ispin = 1, nspins
2443  CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
2444  CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
2445  matrix_wz(ispin)%matrix, focc, nocc)
2446  END DO
2447  END IF
2448  IF (nspins == 2) THEN
2449  CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
2450  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2451  END IF
2452 
2453  IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
2454  IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
2455  NULLIFY (scrm)
2456  CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
2457  matrix_name="OVERLAP MATRIX", &
2458  basis_type_a="ORB", basis_type_b="ORB", &
2459  sab_nl=sab_orb, calculate_forces=.true., &
2460  matrix_p=matrix_wz(1)%matrix)
2461 
2462  IF (SIZE(matrix_wz, 1) == 2) THEN
2463  CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
2464  alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2465  END IF
2466 
2467  IF (debug_forces) THEN
2468  fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
2469  CALL para_env%sum(fodeb)
2470  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wz*dS ", fodeb
2471  END IF
2472  IF (debug_stress .AND. use_virial) THEN
2473  stdeb = fconv*(virial%pv_overlap - stdeb)
2474  CALL para_env%sum(stdeb)
2475  IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2476  'STRESS| WHz ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2477  END IF
2478  CALL dbcsr_deallocate_matrix_set(scrm)
2479 
2480  IF (debug_forces) THEN
2481  CALL total_qs_force(ftot2, force, atomic_kind_set)
2482  fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
2483  CALL para_env%sum(fodeb)
2484  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Response Force", fodeb
2485  fodeb(1:3) = ftot2(1:3, 1)
2486  CALL para_env%sum(fodeb)
2487  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Total Force ", fodeb
2488  DEALLOCATE (ftot1, ftot2, ftot3)
2489  END IF
2490 
2491  IF (do_ex) THEN
2492  CALL dbcsr_deallocate_matrix_set(mpa)
2493  CALL dbcsr_deallocate_matrix_set(matrix_hz)
2494  END IF
2495 
2496  CALL timestop(handle)
2497 
2498  END SUBROUTINE response_force
2499 
2500 ! **************************************************************************************************
2501 !> \brief ...
2502 !> \param qs_env ...
2503 !> \param p_env ...
2504 !> \param matrix_hz ...
2505 !> \param ex_env ...
2506 !> \param debug ...
2507 ! **************************************************************************************************
2508  SUBROUTINE response_force_xtb(qs_env, p_env, matrix_hz, ex_env, debug)
2509  TYPE(qs_environment_type), POINTER :: qs_env
2510  TYPE(qs_p_env_type) :: p_env
2511  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz
2512  TYPE(excited_energy_type), OPTIONAL, POINTER :: ex_env
2513  LOGICAL, INTENT(IN), OPTIONAL :: debug
2514 
2515  CHARACTER(LEN=*), PARAMETER :: routinen = 'response_force_xtb'
2516 
2517  INTEGER :: atom_a, handle, iatom, ikind, iounit, &
2518  is, ispin, na, natom, natorb, nimages, &
2519  nkind, nocc, ns, nsgf, nspins
2520  INTEGER, DIMENSION(25) :: lao
2521  INTEGER, DIMENSION(5) :: occ
2522  LOGICAL :: debug_forces, do_ex, use_virial
2523  REAL(kind=dp) :: focc
2524  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: mcharge, mcharge1
2525  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, charges, charges1, ftot1, &
2526  ftot2
2527  REAL(kind=dp), DIMENSION(3) :: fodeb
2528  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2529  TYPE(cp_logger_type), POINTER :: logger
2530  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pz, matrix_wz, mpa, p_matrix, scrm
2531  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
2532  TYPE(dbcsr_type), POINTER :: s_matrix
2533  TYPE(dft_control_type), POINTER :: dft_control
2534  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2535  TYPE(mp_para_env_type), POINTER :: para_env
2536  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2537  POINTER :: sab_orb
2538  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2539  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2540  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2541  TYPE(qs_ks_env_type), POINTER :: ks_env
2542  TYPE(qs_rho_type), POINTER :: rho
2543  TYPE(xtb_atom_type), POINTER :: xtb_kind
2544 
2545  CALL timeset(routinen, handle)
2546 
2547  IF (PRESENT(debug)) THEN
2548  debug_forces = debug
2549  ELSE
2550  debug_forces = .false.
2551  END IF
2552 
2553  logger => cp_get_default_logger()
2554  IF (logger%para_env%is_source()) THEN
2555  iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
2556  ELSE
2557  iounit = -1
2558  END IF
2559 
2560  do_ex = .false.
2561  IF (PRESENT(ex_env)) do_ex = .true.
2562 
2563  NULLIFY (ks_env, sab_orb)
2564  CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control, &
2565  sab_orb=sab_orb)
2566  CALL get_qs_env(qs_env=qs_env, para_env=para_env, force=force)
2567  nspins = dft_control%nspins
2568 
2569  IF (debug_forces) THEN
2570  CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
2571  ALLOCATE (ftot1(3, natom))
2572  ALLOCATE (ftot2(3, natom))
2573  CALL total_qs_force(ftot1, force, atomic_kind_set)
2574  END IF
2575 
2576  matrix_pz => p_env%p1
2577  NULLIFY (mpa)
2578  IF (do_ex) THEN
2579  CALL dbcsr_allocate_matrix_set(mpa, nspins)
2580  DO ispin = 1, nspins
2581  ALLOCATE (mpa(ispin)%matrix)
2582  CALL dbcsr_create(mpa(ispin)%matrix, template=matrix_pz(ispin)%matrix)
2583  CALL dbcsr_copy(mpa(ispin)%matrix, matrix_pz(ispin)%matrix)
2584  CALL dbcsr_add(mpa(ispin)%matrix, ex_env%matrix_pe(ispin)%matrix, 1.0_dp, 1.0_dp)
2585  CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
2586  END DO
2587  ELSE
2588  mpa => p_env%p1
2589  END IF
2590  !
2591  ! START OF Tr(P+Z)Hcore
2592  !
2593  IF (nspins == 2) THEN
2594  CALL dbcsr_add(mpa(1)%matrix, mpa(2)%matrix, 1.0_dp, 1.0_dp)
2595  END IF
2596  ! Hcore matrix
2597  IF (debug_forces) fodeb(1:3) = force(1)%all_potential(1:3, 1)
2598  CALL xtb_hab_force(qs_env, mpa(1)%matrix)
2599  IF (debug_forces) THEN
2600  fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3)
2601  CALL para_env%sum(fodeb)
2602  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHcore ", fodeb
2603  END IF
2604  IF (nspins == 2) THEN
2605  CALL dbcsr_add(mpa(1)%matrix, mpa(2)%matrix, 1.0_dp, -1.0_dp)
2606  END IF
2607  !
2608  ! END OF Tr(P+Z)Hcore
2609  !
2610  use_virial = .false.
2611  nimages = 1
2612  !
2613  ! Hartree potential of response density
2614  !
2615  IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
2616  ! Mulliken charges
2617  CALL get_qs_env(qs_env, rho=rho, particle_set=particle_set, matrix_s_kp=matrix_s)
2618  natom = SIZE(particle_set)
2619  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2620  ALLOCATE (mcharge(natom), charges(natom, 5))
2621  ALLOCATE (mcharge1(natom), charges1(natom, 5))
2622  charges = 0.0_dp
2623  charges1 = 0.0_dp
2624  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
2625  nkind = SIZE(atomic_kind_set)
2626  CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
2627  ALLOCATE (aocg(nsgf, natom))
2628  aocg = 0.0_dp
2629  ALLOCATE (aocg1(nsgf, natom))
2630  aocg1 = 0.0_dp
2631  p_matrix => matrix_p(:, 1)
2632  s_matrix => matrix_s(1, 1)%matrix
2633  CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
2634  CALL ao_charges(mpa, s_matrix, aocg1, para_env)
2635  DO ikind = 1, nkind
2636  CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
2637  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
2638  CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
2639  DO iatom = 1, na
2640  atom_a = atomic_kind_set(ikind)%atom_list(iatom)
2641  charges(atom_a, :) = real(occ(:), kind=dp)
2642  DO is = 1, natorb
2643  ns = lao(is) + 1
2644  charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
2645  charges1(atom_a, ns) = charges1(atom_a, ns) - aocg1(is, atom_a)
2646  END DO
2647  END DO
2648  END DO
2649  DEALLOCATE (aocg, aocg1)
2650  DO iatom = 1, natom
2651  mcharge(iatom) = sum(charges(iatom, :))
2652  mcharge1(iatom) = sum(charges1(iatom, :))
2653  END DO
2654  ! Coulomb Kernel
2655  CALL xtb_coulomb_hessian(qs_env, matrix_hz, charges1, mcharge1, mcharge)
2656  CALL calc_xtb_ehess_force(qs_env, p_matrix, mpa, charges, mcharge, charges1, &
2657  mcharge1, debug_forces)
2658  !
2659  DEALLOCATE (charges, mcharge, charges1, mcharge1)
2660  END IF
2661  ! Overlap matrix
2662  ! H(drho+dz) + Wz
2663  matrix_wz => p_env%w1
2664  focc = 0.5_dp
2665  IF (nspins == 1) focc = 1.0_dp
2666  CALL get_qs_env(qs_env, mos=mos)
2667  DO ispin = 1, nspins
2668  CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
2669  CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
2670  matrix_wz(ispin)%matrix, focc, nocc)
2671  END DO
2672  IF (nspins == 2) THEN
2673  CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
2674  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2675  END IF
2676  IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
2677  NULLIFY (scrm)
2678  CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
2679  matrix_name="OVERLAP MATRIX", &
2680  basis_type_a="ORB", basis_type_b="ORB", &
2681  sab_nl=sab_orb, calculate_forces=.true., &
2682  matrix_p=matrix_wz(1)%matrix)
2683  IF (debug_forces) THEN
2684  fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
2685  CALL para_env%sum(fodeb)
2686  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wz*dS ", fodeb
2687  END IF
2688  CALL dbcsr_deallocate_matrix_set(scrm)
2689 
2690  IF (debug_forces) THEN
2691  CALL total_qs_force(ftot2, force, atomic_kind_set)
2692  fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
2693  CALL para_env%sum(fodeb)
2694  IF (iounit > 0) WRITE (iounit, "(T3,A,T30,3F16.8)") "DEBUG:: Response Force", fodeb
2695  DEALLOCATE (ftot1, ftot2)
2696  END IF
2697 
2698  IF (do_ex) THEN
2699  CALL dbcsr_deallocate_matrix_set(mpa)
2700  END IF
2701 
2702  CALL timestop(handle)
2703 
2704  END SUBROUTINE response_force_xtb
2705 
2706 ! **************************************************************************************************
2707 !> \brief Win = focc*(P*(H[P_out - P_in] + H[Z] )*P)
2708 !> Langrange multiplier matrix with response and perturbation (Harris) kernel matrices
2709 !>
2710 !> \param qs_env ...
2711 !> \param matrix_hz ...
2712 !> \param matrix_whz ...
2713 !> \param eps_filter ...
2714 !> \param
2715 !> \par History
2716 !> 2020.2 created [Fabian Belleflamme]
2717 !> \author Fabian Belleflamme
2718 ! **************************************************************************************************
2719  SUBROUTINE calculate_whz_ao_matrix(qs_env, matrix_hz, matrix_whz, eps_filter)
2720 
2721  TYPE(qs_environment_type), POINTER :: qs_env
2722  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
2723  POINTER :: matrix_hz
2724  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
2725  POINTER :: matrix_whz
2726  REAL(kind=dp), INTENT(IN) :: eps_filter
2727 
2728  CHARACTER(len=*), PARAMETER :: routinen = 'calculate_whz_ao_matrix'
2729 
2730  INTEGER :: handle, ispin, nspins
2731  REAL(kind=dp) :: scaling
2732  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
2733  TYPE(dbcsr_type) :: matrix_tmp
2734  TYPE(dft_control_type), POINTER :: dft_control
2735  TYPE(mp_para_env_type), POINTER :: para_env
2736  TYPE(qs_rho_type), POINTER :: rho
2737 
2738  CALL timeset(routinen, handle)
2739 
2740  cpassert(ASSOCIATED(qs_env))
2741  cpassert(ASSOCIATED(matrix_hz))
2742  cpassert(ASSOCIATED(matrix_whz))
2743 
2744  CALL get_qs_env(qs_env=qs_env, &
2745  dft_control=dft_control, &
2746  rho=rho, &
2747  para_env=para_env)
2748  nspins = dft_control%nspins
2749  CALL qs_rho_get(rho, rho_ao=rho_ao)
2750 
2751  ! init temp matrix
2752  CALL dbcsr_create(matrix_tmp, template=matrix_hz(1)%matrix, &
2753  matrix_type=dbcsr_type_no_symmetry)
2754 
2755  !Spin factors simplify to
2756  scaling = 1.0_dp
2757  IF (nspins == 1) scaling = 0.5_dp
2758 
2759  ! Operation in MO-solver :
2760  ! Whz = focc*(CC^T*Hz*CC^T)
2761  ! focc = 2.0_dp Closed-shell
2762  ! focc = 1.0_dp Open-shell
2763 
2764  ! Operation in AO-solver :
2765  ! Whz = (scaling*P)*(focc*Hz)*(scaling*P)
2766  ! focc see above
2767  ! scaling = 0.5_dp Closed-shell (P = 2*CC^T), WHz = (0.5*P)*(2*Hz)*(0.5*P)
2768  ! scaling = 1.0_dp Open-shell, WHz = P*Hz*P
2769 
2770  ! Spin factors from Hz and P simplify to
2771  scaling = 1.0_dp
2772  IF (nspins == 1) scaling = 0.5_dp
2773 
2774  DO ispin = 1, nspins
2775 
2776  ! tmp = H*CC^T
2777  CALL dbcsr_multiply("N", "N", scaling, matrix_hz(ispin)%matrix, rho_ao(ispin)%matrix, &
2778  0.0_dp, matrix_tmp, filter_eps=eps_filter)
2779  ! WHz = CC^T*tmp
2780  ! WHz = Wz + (scaling*P)*(focc*Hz)*(scaling*P)
2781  ! WHz = Wz + scaling*(P*Hz*P)
2782  CALL dbcsr_multiply("N", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_tmp, &
2783  1.0_dp, matrix_whz(ispin)%matrix, filter_eps=eps_filter, &
2784  retain_sparsity=.true.)
2785 
2786  END DO
2787 
2788  CALL dbcsr_release(matrix_tmp)
2789 
2790  CALL timestop(handle)
2791 
2792  END SUBROUTINE
2793 
2794 ! **************************************************************************************************
2795 
2796 END MODULE response_solver
det_3x3
real(kind=dp) function det_3x3(a)
...
Definition: dumpdcd.F:1091
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
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition: atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
core_ae
Calculation of the nuclear attraction contribution to the core Hamiltonian <a|erfc|b> :we only calcul...
Definition: core_ae.F:14
core_ae::build_core_ae
subroutine, public build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index)
...
Definition: core_ae.F:84
core_ppl
Calculation of the local pseudopotential contribution to the core Hamiltonian <a|V(local)|b> = <a|Sum...
Definition: core_ppl.F:18
core_ppl::build_core_ppl
subroutine, public build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, nimages, cell_to_index, basis_type, deltaR)
...
Definition: core_ppl.F:95
core_ppnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition: core_ppnl.F:15
core_ppnl::build_core_ppnl
subroutine, public build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, basis_type, deltaR, matrix_l)
...
Definition: core_ppnl.F:89
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.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::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_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_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_init_random
subroutine, public cp_fm_init_random(matrix, ncol, start_col)
fills a matrix with random numbers
Definition: cp_fm_types.F:452
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
ec_env_types
Types needed for a for a Energy Correction.
Definition: ec_env_types.F:14
ec_methods
Routines used for Harris functional Kohn-Sham calculation.
Definition: ec_methods.F:15
ec_methods::ec_mos_init
subroutine, public ec_mos_init(qs_env, matrix_s)
Allocate and initiate molecular orbitals environment.
Definition: ec_methods.F:162
ec_methods::create_kernel
subroutine, public create_kernel(qs_env, vxc, vxc_tau, rho, rho1_r, rho1_g, tau1_r, xc_section, compute_virial, virial_xc)
Creation of second derivative xc-potential.
Definition: ec_methods.F:88
ec_orth_solver
AO-based conjugate-gradient response solver routines.
Definition: ec_orth_solver.F:15
ec_orth_solver::ec_response_ao
subroutine, public ec_response_ao(qs_env, p_env, matrix_hz, matrix_pz, matrix_wz, iounit, should_stop)
AO-based conjugate gradient linear response solver. In goes the right hand side B of the equation AZ=...
Definition: ec_orth_solver.F:357
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::precond_mlp
integer, parameter, public precond_mlp
Definition: input_constants.F:513
input_constants::kg_tnadd_embed_ri
integer, parameter, public kg_tnadd_embed_ri
Definition: input_constants.F:1130
input_constants::kg_tnadd_embed
integer, parameter, public kg_tnadd_embed
Definition: input_constants.F:1130
input_constants::ec_mo_solver
integer, parameter, public ec_mo_solver
Definition: input_constants.F:1154
input_constants::ot_precond_full_kinetic
integer, parameter, public ot_precond_full_kinetic
Definition: input_constants.F:513
input_constants::kg_tnadd_atomic
integer, parameter, public kg_tnadd_atomic
Definition: input_constants.F:1130
input_constants::do_admm_aux_exch_func_none
integer, parameter, public do_admm_aux_exch_func_none
Definition: input_constants.F:863
input_constants::ls_s_sqrt_proot
integer, parameter, public ls_s_sqrt_proot
Definition: input_constants.F:957
input_constants::ot_precond_full_single
integer, parameter, public ot_precond_full_single
Definition: input_constants.F:513
input_constants::ls_s_sqrt_ns
integer, parameter, public ls_s_sqrt_ns
Definition: input_constants.F:957
input_constants::ot_precond_none
integer, parameter, public ot_precond_none
Definition: input_constants.F:513
input_constants::ot_precond_full_single_inverse
integer, parameter, public ot_precond_full_single_inverse
Definition: input_constants.F:513
input_constants::xc_none
integer, parameter, public xc_none
Definition: input_constants.F:551
input_constants::ot_precond_s_inverse
integer, parameter, public ot_precond_s_inverse
Definition: input_constants.F:513
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition: input_constants.F:513
input_constants::ec_ls_solver
integer, parameter, public ec_ls_solver
Definition: input_constants.F:1154
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
kg_correction
Routines for a Kim-Gordon-like partitioning into molecular subunits.
Definition: kg_correction.F:14
kg_correction::kg_ekin_subset
subroutine, public kg_ekin_subset(qs_env, ks_matrix, ekin_mol, calc_force, do_kernel, pmat_ext)
Calculates the subsystem Hohenberg-Kohn kinetic energy and the forces.
Definition: kg_correction.F:88
kg_environment_types
Types needed for a Kim-Gordon-like partitioning into molecular subunits.
Definition: kg_environment_types.F:15
kg_tnadd_mat
Calculation of the local potential contribution of the nonadditive kinetic energy <a|V(local)|b> = <a...
Definition: kg_tnadd_mat.F:13
kg_tnadd_mat::build_tnadd_mat
subroutine, public build_tnadd_mat(kg_env, matrix_p, force, virial, calculate_forces, use_virial, qs_kind_set, atomic_kind_set, particle_set, sab_orb, dbcsr_dist)
...
Definition: kg_tnadd_mat.F:85
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
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
mathlib
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
mulliken
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition: mulliken.F:13
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::pascal
real(kind=dp), parameter, public pascal
Definition: physcon.F:174
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_2nd_kernel_ao
Routines to calculate 2nd order kernels from a given response density in ao basis linear response scf...
Definition: qs_2nd_kernel_ao.F:15
qs_2nd_kernel_ao::build_dm_response
subroutine, public build_dm_response(c0, c1, dm)
This routine builds response density in dbcsr format.
Definition: qs_2nd_kernel_ao.F:81
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_density_matrices
collects routines that calculate density matrices
Definition: qs_density_matrices.F:14
qs_density_matrices::calculate_whz_matrix
subroutine, public calculate_whz_matrix(c0vec, hzm, w_matrix, focc, nocc)
Calculate the Wz matrix from the MO eigenvectors, MO eigenvalues, and the MO occupation numbers....
Definition: qs_density_matrices.F:425
qs_density_matrices::calculate_wz_matrix
subroutine, public calculate_wz_matrix(mo_set, psi1, ks_matrix, w_matrix)
Calculate the response W matrix from the MO eigenvectors, MO eigenvalues, and the MO occupation numbe...
Definition: qs_density_matrices.F:364
qs_energy_types
Definition: qs_energy_types.F:14
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_force_types::total_qs_force
subroutine, public total_qs_force(force, qs_force, atomic_kind_set)
Get current total force.
Definition: qs_force_types.F:540
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_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
qs_kinetic
Calculation of kinetic energy matrix and forces.
Definition: qs_kinetic.F:15
qs_kinetic::build_kinetic_matrix
subroutine, public build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, nderivative)
Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
Definition: qs_kinetic.F:101
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_methods
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition: qs_ks_methods.F:22
qs_ks_methods::calc_rho_tot_gspace
subroutine, public calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density)
...
Definition: qs_ks_methods.F:916
qs_ks_types
Definition: qs_ks_types.F:15
qs_linres_methods
localize wavefunctions linear response scf
Definition: qs_linres_methods.F:15
qs_linres_methods::linres_solver
subroutine, public linres_solver(p_env, qs_env, psi1, h1_psi0, psi0_order, iounit, should_stop)
scf loop to optimize the first order wavefunctions (psi1) given a perturbation as an operator applied...
Definition: qs_linres_methods.F:216
qs_linres_types
Type definitiona for linear response calculations.
Definition: qs_linres_types.F:12
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_matrix_pools
wrapper for the pools of matrixes
Definition: qs_matrix_pools.F:14
qs_matrix_pools::mpools_rebuild_fm_pools
subroutine, public mpools_rebuild_fm_pools(mpools, mos, blacs_env, para_env, nmosub)
rebuilds the pools of the (ao x mo, ao x ao , mo x mo) full matrixes
Definition: qs_matrix_pools.F:211
qs_mo_methods
collects routines that perform operations directly related to MOs
Definition: qs_mo_methods.F:14
qs_mo_methods::make_basis_sm
subroutine, public make_basis_sm(vmatrix, ncol, matrix_s)
returns an S-orthonormal basis v (v^T S v ==1)
Definition: qs_mo_methods.F:87
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::deallocate_mo_set
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition: qs_mo_types.F:352
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_oce_types
Definition: qs_oce_types.F:8
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_p_env_methods
Utility functions for the perturbation calculations.
Definition: qs_p_env_methods.F:16
qs_p_env_methods::p_env_psi0_changed
subroutine, public p_env_psi0_changed(p_env, qs_env)
To be called after the value of psi0 has changed. Recalculates the quantities S_psi0 and m_epsilon.
Definition: qs_p_env_methods.F:482
qs_p_env_methods::p_env_create
subroutine, public p_env_create(p_env, qs_env, p1_option, p1_admm_option, orthogonal_orbitals, linres_control)
allocates and initializes the perturbation environment (no setup)
Definition: qs_p_env_methods.F:133
qs_p_env_types
basis types for the calculation of the perturbation of density theory.
Definition: qs_p_env_types.F:14
qs_p_env_types::p_env_release
subroutine, public p_env_release(p_env)
relases the given p_env (see doc/ReferenceCounting.html)
Definition: qs_p_env_types.F:99
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_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_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_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_xc_2nd_deriv_atom
subroutine, public calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, do_tddft, do_tddfpt2, do_triplet, kind_set_external)
...
Definition: qs_vxc_atom.F:446
qs_vxc_atom::calculate_vxc_atom
subroutine, public calculate_vxc_atom(qs_env, energy_only, exc1, gradient_atom_set, adiabatic_rescale_factor, kind_set_external, rho_atom_set_external, xc_section_external)
...
Definition: qs_vxc_atom.F:87
response_solver
Calculate the CPKS equation and the resulting forces.
Definition: response_solver.F:17
response_solver::response_equation_new
subroutine, public response_equation_new(qs_env, p_env, cpmos, iounit)
Initializes vectors for MO-coefficient based linear response solver and calculates response density,...
Definition: response_solver.F:551
response_solver::response_force_xtb
subroutine, public response_force_xtb(qs_env, p_env, matrix_hz, ex_env, debug)
...
Definition: response_solver.F:2509
response_solver::response_force
subroutine, public response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, matrix_hz, matrix_pz, matrix_pz_admm, matrix_wz, zehartree, zexc, zexc_aux_fit, rhopz_r, p_env, ex_env, debug)
...
Definition: response_solver.F:813
response_solver::response_equation
subroutine, public response_equation(qs_env, p_env, cpmos, iounit, lr_section)
Initializes vectors for MO-coefficient based linear response solver and calculates response density,...
Definition: response_solver.F:651
response_solver::response_calculation
subroutine, public response_calculation(qs_env, ec_env)
Initializes solver of linear response equation for energy correction.
Definition: response_solver.F:181
task_list_types
types for task lists
Definition: task_list_types.F:14
virial_methods
Definition: virial_methods.F:12
virial_methods::one_third_sum_diag
pure real(kind=dp) function, public one_third_sum_diag(a)
...
Definition: virial_methods.F:238
virial_types
Definition: virial_types.F:13
xtb_ehess_force
Calculation of forces for Coulomb contributions in response xTB.
Definition: xtb_ehess_force.F:12
xtb_ehess_force::calc_xtb_ehess_force
subroutine, public calc_xtb_ehess_force(qs_env, matrix_p0, matrix_p1, charges0, mcharge0, charges1, mcharge1, debug_forces)
...
Definition: xtb_ehess_force.F:86
xtb_ehess
Calculation of Coulomb Hessian contributions in xTB.
Definition: xtb_ehess.F:12
xtb_ehess::xtb_coulomb_hessian
subroutine, public xtb_coulomb_hessian(qs_env, ks_matrix, charges1, mcharge1, mcharge)
...
Definition: xtb_ehess.F:77
xtb_matrices
Calculation of Overlap and Hamiltonian matrices in xTB Reference: Stefan Grimme, Christoph Bannwarth,...
Definition: xtb_matrices.F:15
xtb_matrices::xtb_hab_force
subroutine, public xtb_hab_force(qs_env, p_matrix)
...
Definition: xtb_matrices.F:987
xtb_types
Definition of the xTB parameter types.
Definition: xtb_types.F:20
xtb_types::get_xtb_atom_param
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, hen, zeta, occupation, electronegativity, chmax)
...
Definition: xtb_types.F:175