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