(git:936074a)
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_functional_ext, ec_ls_solver, ec_mo_solver, &
64 kg_tnadd_atomic, kg_tnadd_embed, kg_tnadd_embed_ri, ls_s_sqrt_ns, ls_s_sqrt_proot, &
65 ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, &
66 ot_precond_full_single_inverse, 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 IF (ec_env%energy_functional == ec_functional_ext) THEN
396 cpabort("AO Response Solver NYA for External Functional")
397 END IF
398
399 ! AO ortho solver
400 CALL ec_response_ao(qs_env=qs_env, &
401 p_env=p_env, &
402 matrix_hz=ec_env%matrix_hz, &
403 matrix_pz=ec_env%matrix_z, &
404 matrix_wz=ec_env%matrix_wz, &
405 iounit=unit_nr, &
406 should_stop=should_stop)
407
408 IF (dft_control%do_admm) THEN
409 CALL get_qs_env(qs_env, admm_env=admm_env)
410 cpassert(ASSOCIATED(admm_env%work_orb_orb))
411 cpassert(ASSOCIATED(admm_env%work_aux_orb))
412 cpassert(ASSOCIATED(admm_env%work_aux_aux))
413 nao = admm_env%nao_orb
414 nao_aux = admm_env%nao_aux_fit
415 DO ispin = 1, nspins
416 CALL copy_dbcsr_to_fm(ec_env%matrix_z(ispin)%matrix, admm_env%work_orb_orb)
417 CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
418 1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
419 admm_env%work_aux_orb)
420 CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
421 1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
422 admm_env%work_aux_aux)
423 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, p_env%p1_admm(ispin)%matrix, &
424 keep_sparsity=.true.)
425 END DO
426 END IF
427
428 CASE DEFAULT
429 cpabort("Unknown solver for response equation requested")
430 END SELECT
431
432 IF (dft_control%do_admm) THEN
433 CALL dbcsr_allocate_matrix_set(ec_env%z_admm, nspins)
434 DO ispin = 1, nspins
435 ALLOCATE (ec_env%z_admm(ispin)%matrix)
436 CALL dbcsr_create(matrix=ec_env%z_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
437 CALL get_qs_env(qs_env, admm_env=admm_env)
438 CALL dbcsr_copy(ec_env%z_admm(ispin)%matrix, p_env%p1_admm(ispin)%matrix)
439 END DO
440 END IF
441
442 ! Get rid of MO environment again
443 IF (dft_control%qs_control%do_ls_scf) THEN
444 DO ispin = 1, nspins
445 CALL deallocate_mo_set(mos(ispin))
446 END DO
447 IF (ASSOCIATED(qs_env%mos)) THEN
448 DO ispin = 1, SIZE(qs_env%mos)
449 CALL deallocate_mo_set(qs_env%mos(ispin))
450 END DO
451 DEALLOCATE (qs_env%mos)
452 END IF
453 END IF
454
455 CALL timestop(handle)
456
457 END SUBROUTINE response_calculation
458
459! **************************************************************************************************
460!> \brief Parse the input section of the response solver
461!> \param input Input section which controls response solver parameters
462!> \param linres_control Environment for general setting of linear response calculation
463!> \param unit_nr ...
464!> \par History
465!> 2020.05 created [Fabian Belleflamme]
466!> \author Fabian Belleflamme
467! **************************************************************************************************
468 SUBROUTINE response_solver_write_input(input, linres_control, unit_nr)
469 TYPE(section_vals_type), POINTER :: input
470 TYPE(linres_control_type), POINTER :: linres_control
471 INTEGER, INTENT(IN) :: unit_nr
472
473 CHARACTER(len=*), PARAMETER :: routinen = 'response_solver_write_input'
474
475 INTEGER :: handle, max_iter_lanczos, s_sqrt_method, &
476 s_sqrt_order, solver_method
477 REAL(kind=dp) :: eps_lanczos
478
479 CALL timeset(routinen, handle)
480
481 IF (unit_nr > 0) THEN
482
483 ! linres_control
484 WRITE (unit_nr, '(/,T2,A)') &
485 repeat("-", 30)//" Linear Response Solver "//repeat("-", 25)
486
487 ! Which type of solver is used
488 CALL section_vals_val_get(input, "METHOD", i_val=solver_method)
489
490 SELECT CASE (solver_method)
491 CASE (ec_ls_solver)
492 WRITE (unit_nr, '(T2,A,T61,A20)') "Solver: ", "AO-based CG-solver"
493 CASE (ec_mo_solver)
494 WRITE (unit_nr, '(T2,A,T61,A20)') "Solver: ", "MO-based CG-solver"
495 END SELECT
496
497 WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps:", linres_control%eps
498 WRITE (unit_nr, '(T2,A,T61,E20.3)') "eps_filter:", linres_control%eps_filter
499 WRITE (unit_nr, '(T2,A,T61,I20)') "Max iter:", linres_control%max_iter
500
501 SELECT CASE (linres_control%preconditioner_type)
502 CASE (ot_precond_full_all)
503 WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "FULL_ALL"
504 CASE (ot_precond_full_single_inverse)
505 WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "FULL_SINGLE_INVERSE"
506 CASE (ot_precond_full_single)
507 WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "FULL_SINGLE"
508 CASE (ot_precond_full_kinetic)
509 WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "FULL_KINETIC"
510 CASE (ot_precond_s_inverse)
511 WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "FULL_S_INVERSE"
512 CASE (precond_mlp)
513 WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "MULTI_LEVEL"
514 CASE (ot_precond_none)
515 WRITE (unit_nr, '(T2,A,T61,A20)') "Preconditioner: ", "NONE"
516 END SELECT
517
518 SELECT CASE (solver_method)
519 CASE (ec_ls_solver)
520
521 CALL section_vals_val_get(input, "S_SQRT_METHOD", i_val=s_sqrt_method)
522 CALL section_vals_val_get(input, "S_SQRT_ORDER", i_val=s_sqrt_order)
523 CALL section_vals_val_get(input, "EPS_LANCZOS", r_val=eps_lanczos)
524 CALL section_vals_val_get(input, "MAX_ITER_LANCZOS", i_val=max_iter_lanczos)
525
526 ! Response solver transforms P and KS into orthonormal basis,
527 ! reuires matrx S sqrt and its inverse
528 SELECT CASE (s_sqrt_method)
529 CASE (ls_s_sqrt_ns)
530 WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "NEWTONSCHULZ"
531 CASE (ls_s_sqrt_proot)
532 WRITE (unit_nr, '(T2,A,T61,A20)') "S sqrt method:", "PROOT"
533 CASE DEFAULT
534 cpabort("Unknown sqrt method.")
535 END SELECT
536 WRITE (unit_nr, '(T2,A,T61,I20)') "S sqrt order:", s_sqrt_order
537
538 CASE (ec_mo_solver)
539 END SELECT
540
541 WRITE (unit_nr, '(T2,A)') repeat("-", 79)
542
543 CALL m_flush(unit_nr)
544 END IF
545
546 CALL timestop(handle)
547
548 END SUBROUTINE response_solver_write_input
549
550! **************************************************************************************************
551!> \brief Initializes vectors for MO-coefficient based linear response solver
552!> and calculates response density, and energy-weighted response density matrix
553!>
554!> \param qs_env ...
555!> \param p_env ...
556!> \param cpmos ...
557!> \param iounit ...
558! **************************************************************************************************
559 SUBROUTINE response_equation_new(qs_env, p_env, cpmos, iounit)
560 TYPE(qs_environment_type), POINTER :: qs_env
561 TYPE(qs_p_env_type) :: p_env
562 TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: cpmos
563 INTEGER, INTENT(IN) :: iounit
564
565 CHARACTER(LEN=*), PARAMETER :: routinen = 'response_equation_new'
566
567 INTEGER :: handle, ispin, nao, nao_aux, nocc, nspins
568 LOGICAL :: should_stop
569 TYPE(admm_type), POINTER :: admm_env
570 TYPE(cp_fm_struct_type), POINTER :: fm_struct
571 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: psi0, psi1
572 TYPE(cp_fm_type), POINTER :: mo_coeff
573 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
574 TYPE(dft_control_type), POINTER :: dft_control
575 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
576
577 CALL timeset(routinen, handle)
578
579 NULLIFY (dft_control, matrix_ks, mo_coeff, mos)
580
581 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks=matrix_ks, &
582 matrix_s=matrix_s, mos=mos)
583 nspins = dft_control%nspins
584
585 ! Initialize vectors:
586 ! psi0 : The ground-state MO-coefficients
587 ! psi1 : The "perturbed" linear response orbitals
588 ALLOCATE (psi0(nspins), psi1(nspins))
589 DO ispin = 1, nspins
590 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, homo=nocc)
591 NULLIFY (fm_struct)
592 CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, &
593 template_fmstruct=mo_coeff%matrix_struct)
594 CALL cp_fm_create(psi0(ispin), fm_struct)
595 CALL cp_fm_to_fm(mo_coeff, psi0(ispin), nocc)
596 CALL cp_fm_create(psi1(ispin), fm_struct)
597 CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
598 CALL cp_fm_struct_release(fm_struct)
599 END DO
600
601 should_stop = .false.
602 ! The response solver
603 CALL linres_solver(p_env, qs_env, psi1, cpmos, psi0, iounit, should_stop)
604
605 ! Building the response density matrix
606 DO ispin = 1, nspins
607 CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_s(1)%matrix)
608 END DO
609 CALL build_dm_response(psi0, psi1, p_env%p1)
610 DO ispin = 1, nspins
611 CALL dbcsr_scale(p_env%p1(ispin)%matrix, 0.5_dp)
612 END DO
613
614 IF (dft_control%do_admm) THEN
615 CALL get_qs_env(qs_env, admm_env=admm_env)
616 cpassert(ASSOCIATED(admm_env%work_orb_orb))
617 cpassert(ASSOCIATED(admm_env%work_aux_orb))
618 cpassert(ASSOCIATED(admm_env%work_aux_aux))
619 nao = admm_env%nao_orb
620 nao_aux = admm_env%nao_aux_fit
621 DO ispin = 1, nspins
622 CALL copy_dbcsr_to_fm(p_env%p1(ispin)%matrix, admm_env%work_orb_orb)
623 CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
624 1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
625 admm_env%work_aux_orb)
626 CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
627 1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
628 admm_env%work_aux_aux)
629 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, p_env%p1_admm(ispin)%matrix, &
630 keep_sparsity=.true.)
631 END DO
632 END IF
633
634 ! Calculate Wz = 0.5*(psi1*eps*psi0^T + psi0*eps*psi1^T)
635 DO ispin = 1, nspins
636 CALL calculate_wz_matrix(mos(ispin), psi1(ispin), matrix_ks(ispin)%matrix, &
637 p_env%w1(ispin)%matrix)
638 END DO
639 DO ispin = 1, nspins
640 CALL cp_fm_release(cpmos(ispin))
641 END DO
642 CALL cp_fm_release(psi1)
643 CALL cp_fm_release(psi0)
644
645 CALL timestop(handle)
646
647 END SUBROUTINE response_equation_new
648
649! **************************************************************************************************
650!> \brief Initializes vectors for MO-coefficient based linear response solver
651!> and calculates response density, and energy-weighted response density matrix
652!> J. Chem. Theory Comput. 2022, 18, 4186−4202 (https://doi.org/10.1021/acs.jctc.2c00144)
653!>
654!> \param qs_env ...
655!> \param p_env Holds the two results of this routine, p_env%p1 = CZ^T + ZC^T,
656!> p_env%w1 = 0.5\sum_i(C_i*\epsilon_i*Z_i^T + Z_i*\epsilon_i*C_i^T)
657!> \param cpmos RHS of equation as Ax + b = 0 (sign of b)
658!> \param iounit ...
659!> \param lr_section ...
660! **************************************************************************************************
661 SUBROUTINE response_equation(qs_env, p_env, cpmos, iounit, lr_section)
662 TYPE(qs_environment_type), POINTER :: qs_env
663 TYPE(qs_p_env_type) :: p_env
664 TYPE(cp_fm_type), DIMENSION(:), POINTER :: cpmos
665 INTEGER, INTENT(IN) :: iounit
666 TYPE(section_vals_type), OPTIONAL, POINTER :: lr_section
667
668 CHARACTER(LEN=*), PARAMETER :: routinen = 'response_equation'
669
670 INTEGER :: handle, ispin, nao, nao_aux, nocc, nspins
671 LOGICAL :: should_stop
672 TYPE(admm_type), POINTER :: admm_env
673 TYPE(cp_fm_struct_type), POINTER :: fm_struct
674 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: psi0, psi1
675 TYPE(cp_fm_type), POINTER :: mo_coeff
676 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_s_aux
677 TYPE(dft_control_type), POINTER :: dft_control
678 TYPE(linres_control_type), POINTER :: linres_control
679 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
680 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
681 POINTER :: sab_orb
682
683 CALL timeset(routinen, handle)
684
685 ! initialized linres_control
686 NULLIFY (linres_control)
687 ALLOCATE (linres_control)
688 linres_control%do_kernel = .true.
689 linres_control%lr_triplet = .false.
690 IF (PRESENT(lr_section)) THEN
691 CALL section_vals_val_get(lr_section, "RESTART", l_val=linres_control%linres_restart)
692 CALL section_vals_val_get(lr_section, "MAX_ITER", i_val=linres_control%max_iter)
693 CALL section_vals_val_get(lr_section, "EPS", r_val=linres_control%eps)
694 CALL section_vals_val_get(lr_section, "EPS_FILTER", r_val=linres_control%eps_filter)
695 CALL section_vals_val_get(lr_section, "RESTART_EVERY", i_val=linres_control%restart_every)
696 CALL section_vals_val_get(lr_section, "PRECONDITIONER", i_val=linres_control%preconditioner_type)
697 CALL section_vals_val_get(lr_section, "ENERGY_GAP", r_val=linres_control%energy_gap)
698 ELSE
699 linres_control%linres_restart = .false.
700 linres_control%max_iter = 100
701 linres_control%eps = 1.0e-10_dp
702 linres_control%eps_filter = 1.0e-15_dp
703 linres_control%restart_every = 50
704 linres_control%preconditioner_type = ot_precond_full_single_inverse
705 linres_control%energy_gap = 0.02_dp
706 END IF
707
708 ! initialized p_env
709 CALL p_env_create(p_env, qs_env, orthogonal_orbitals=.true., &
710 linres_control=linres_control)
711 CALL set_qs_env(qs_env, linres_control=linres_control)
712 CALL p_env_psi0_changed(p_env, qs_env)
713 p_env%new_preconditioner = .true.
714
715 CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos)
716 !
717 nspins = dft_control%nspins
718
719 ! Initialize vectors:
720 ! psi0 : The ground-state MO-coefficients
721 ! psi1 : The "perturbed" linear response orbitals
722 ALLOCATE (psi0(nspins), psi1(nspins))
723 DO ispin = 1, nspins
724 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, homo=nocc)
725 NULLIFY (fm_struct)
726 CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, &
727 template_fmstruct=mo_coeff%matrix_struct)
728 CALL cp_fm_create(psi0(ispin), fm_struct)
729 CALL cp_fm_to_fm(mo_coeff, psi0(ispin), nocc)
730 CALL cp_fm_create(psi1(ispin), fm_struct)
731 CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
732 CALL cp_fm_struct_release(fm_struct)
733 END DO
734
735 should_stop = .false.
736 ! The response solver
737 CALL get_qs_env(qs_env, matrix_s=matrix_s, sab_orb=sab_orb)
738 CALL dbcsr_allocate_matrix_set(p_env%p1, nspins)
739 CALL dbcsr_allocate_matrix_set(p_env%w1, nspins)
740 DO ispin = 1, nspins
741 ALLOCATE (p_env%p1(ispin)%matrix, p_env%w1(ispin)%matrix)
742 CALL dbcsr_create(matrix=p_env%p1(ispin)%matrix, template=matrix_s(1)%matrix)
743 CALL dbcsr_create(matrix=p_env%w1(ispin)%matrix, template=matrix_s(1)%matrix)
744 CALL cp_dbcsr_alloc_block_from_nbl(p_env%p1(ispin)%matrix, sab_orb)
745 CALL cp_dbcsr_alloc_block_from_nbl(p_env%w1(ispin)%matrix, sab_orb)
746 END DO
747 IF (dft_control%do_admm) THEN
748 CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux)
749 CALL dbcsr_allocate_matrix_set(p_env%p1_admm, nspins)
750 DO ispin = 1, nspins
751 ALLOCATE (p_env%p1_admm(ispin)%matrix)
752 CALL dbcsr_create(p_env%p1_admm(ispin)%matrix, &
753 template=matrix_s_aux(1)%matrix)
754 CALL dbcsr_copy(p_env%p1_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
755 CALL dbcsr_set(p_env%p1_admm(ispin)%matrix, 0.0_dp)
756 END DO
757 END IF
758
759 CALL linres_solver(p_env, qs_env, psi1, cpmos, psi0, iounit, should_stop)
760
761 ! Building the response density matrix
762 DO ispin = 1, nspins
763 CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_s(1)%matrix)
764 END DO
765 CALL build_dm_response(psi0, psi1, p_env%p1)
766 DO ispin = 1, nspins
767 CALL dbcsr_scale(p_env%p1(ispin)%matrix, 0.5_dp)
768 END DO
769 IF (dft_control%do_admm) THEN
770 CALL get_qs_env(qs_env, admm_env=admm_env)
771 cpassert(ASSOCIATED(admm_env%work_orb_orb))
772 cpassert(ASSOCIATED(admm_env%work_aux_orb))
773 cpassert(ASSOCIATED(admm_env%work_aux_aux))
774 nao = admm_env%nao_orb
775 nao_aux = admm_env%nao_aux_fit
776 DO ispin = 1, nspins
777 CALL copy_dbcsr_to_fm(p_env%p1(ispin)%matrix, admm_env%work_orb_orb)
778 CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
779 1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
780 admm_env%work_aux_orb)
781 CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
782 1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
783 admm_env%work_aux_aux)
784 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, p_env%p1_admm(ispin)%matrix, &
785 keep_sparsity=.true.)
786 END DO
787 END IF
788
789 ! Calculate the second term of Eq. 51 Wz = 0.5*(psi1*eps*psi0^T + psi0*eps*psi1^T)
790 CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
791 DO ispin = 1, nspins
792 CALL calculate_wz_matrix(mos(ispin), psi1(ispin), matrix_ks(ispin)%matrix, &
793 p_env%w1(ispin)%matrix)
794 END DO
795 CALL cp_fm_release(psi0)
796 CALL cp_fm_release(psi1)
797
798 CALL timestop(handle)
799
800 END SUBROUTINE response_equation
801
802! **************************************************************************************************
803!> \brief ...
804!> \param qs_env ...
805!> \param vh_rspace ...
806!> \param vxc_rspace ...
807!> \param vtau_rspace ...
808!> \param vadmm_rspace ...
809!> \param matrix_hz Right-hand-side of linear response equation
810!> \param matrix_pz Linear response density matrix
811!> \param matrix_pz_admm Linear response density matrix in ADMM basis
812!> \param matrix_wz Energy-weighted linear response density
813!> \param zehartree Hartree volume response contribution to stress tensor
814!> \param zexc XC volume response contribution to stress tensor
815!> \param zexc_aux_fit ADMM XC volume response contribution to stress tensor
816!> \param rhopz_r Response density on real space grid
817!> \param p_env ...
818!> \param ex_env ...
819!> \param debug ...
820! **************************************************************************************************
821 SUBROUTINE response_force(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, &
822 matrix_hz, matrix_pz, matrix_pz_admm, matrix_wz, &
823 zehartree, zexc, zexc_aux_fit, rhopz_r, p_env, ex_env, debug)
824 TYPE(qs_environment_type), POINTER :: qs_env
825 TYPE(pw_r3d_rs_type), INTENT(IN) :: vh_rspace
826 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rspace, vtau_rspace, vadmm_rspace
827 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz, matrix_pz, matrix_pz_admm, &
828 matrix_wz
829 REAL(kind=dp), OPTIONAL :: zehartree, zexc, zexc_aux_fit
830 TYPE(pw_r3d_rs_type), DIMENSION(:), &
831 INTENT(INOUT), OPTIONAL :: rhopz_r
832 TYPE(qs_p_env_type), OPTIONAL :: p_env
833 TYPE(excited_energy_type), OPTIONAL, POINTER :: ex_env
834 LOGICAL, INTENT(IN), OPTIONAL :: debug
835
836 CHARACTER(LEN=*), PARAMETER :: routinen = 'response_force'
837
838 CHARACTER(LEN=default_string_length) :: basis_type, unitstr
839 INTEGER :: handle, iounit, ispin, mspin, myfun, &
840 n_rep_hf, nao, nao_aux, natom, nder, &
841 nocc, nspins
842 LOGICAL :: debug_forces, debug_stress, distribute_fock_matrix, do_ex, do_hfx, gapw, gapw_xc, &
843 hfx_treat_lsd_in_core, resp_only, s_mstruct_changed, use_virial
844 REAL(kind=dp) :: eh1, ehartree, ekin_mol, eps_filter, &
845 exc, exc_aux_fit, fconv, focc, &
846 hartree_gs, hartree_t
847 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ftot1, ftot2, ftot3
848 REAL(kind=dp), DIMENSION(2) :: total_rho_gs, total_rho_t
849 REAL(kind=dp), DIMENSION(3) :: fodeb
850 REAL(kind=dp), DIMENSION(3, 3) :: h_stress, pv_loc, stdeb, sttot, sttot2
851 TYPE(admm_type), POINTER :: admm_env
852 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
853 TYPE(cell_type), POINTER :: cell
854 TYPE(cp_logger_type), POINTER :: logger
855 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
856 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ht, matrix_pd, matrix_pza, &
857 matrix_s, mpa, scrm
858 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p, mhd, mhx, mhy, mhz, &
859 mpa2, mpd, mpz, scrm2
860 TYPE(dbcsr_type), POINTER :: dbwork
861 TYPE(dft_control_type), POINTER :: dft_control
862 TYPE(hartree_local_type), POINTER :: hartree_local_gs, hartree_local_t
863 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
864 TYPE(kg_environment_type), POINTER :: kg_env
865 TYPE(local_rho_type), POINTER :: local_rho_set_f, local_rho_set_gs, &
866 local_rho_set_t, local_rho_set_vxc, &
867 local_rhoz_set_admm
868 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
869 TYPE(mp_para_env_type), POINTER :: para_env
870 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
871 POINTER :: sab_aux_fit, sab_orb
872 TYPE(oce_matrix_type), POINTER :: oce
873 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
874 TYPE(pw_c1d_gs_type) :: rho_tot_gspace, rho_tot_gspace_gs, rho_tot_gspace_t, &
875 rhoz_tot_gspace, v_hartree_gspace_gs, v_hartree_gspace_t, zv_hartree_gspace
876 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_aux, rho_g_gs, rho_g_t, rhoz_g, &
877 rhoz_g_aux, rhoz_g_xc
878 TYPE(pw_c1d_gs_type), POINTER :: rho_core
879 TYPE(pw_env_type), POINTER :: pw_env
880 TYPE(pw_poisson_type), POINTER :: poisson_env
881 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
882 TYPE(pw_r3d_rs_type) :: v_hartree_rspace_gs, v_hartree_rspace_t, &
883 vhxc_rspace, zv_hartree_rspace
884 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_aux, rho_r_gs, rho_r_t, rhoz_r, &
885 rhoz_r_aux, rhoz_r_xc, tau_r_aux, &
886 tauz_r, tauz_r_xc, v_xc, v_xc_tau
887 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
888 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
889 TYPE(qs_ks_env_type), POINTER :: ks_env
890 TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit, rho_xc
891 TYPE(section_vals_type), POINTER :: hfx_section, xc_fun_section, xc_section
892 TYPE(task_list_type), POINTER :: task_list, task_list_aux_fit
893 TYPE(virial_type), POINTER :: virial
894
895 CALL timeset(routinen, handle)
896
897 IF (PRESENT(debug)) THEN
898 debug_forces = debug
899 debug_stress = debug
900 ELSE
901 debug_forces = .false.
902 debug_stress = .false.
903 END IF
904
905 logger => cp_get_default_logger()
906 IF (logger%para_env%is_source()) THEN
907 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
908 ELSE
909 iounit = -1
910 END IF
911
912 do_ex = .false.
913 IF (PRESENT(ex_env)) do_ex = .true.
914 IF (do_ex) THEN
915 cpassert(PRESENT(p_env))
916 END IF
917
918 NULLIFY (ks_env, sab_orb, virial)
919 CALL get_qs_env(qs_env=qs_env, &
920 cell=cell, &
921 force=force, &
922 ks_env=ks_env, &
923 dft_control=dft_control, &
924 para_env=para_env, &
925 sab_orb=sab_orb, &
926 virial=virial)
927 nspins = dft_control%nspins
928 gapw = dft_control%qs_control%gapw
929 gapw_xc = dft_control%qs_control%gapw_xc
930
931 IF (debug_forces) THEN
932 CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
933 ALLOCATE (ftot1(3, natom))
934 CALL total_qs_force(ftot1, force, atomic_kind_set)
935 END IF
936
937 ! check for virial
938 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
939
940 IF (use_virial .AND. do_ex) THEN
941 CALL cp_abort(__location__, "Stress Tensor not available for TDDFT calculations.")
942 END IF
943
944 fconv = 1.0e-9_dp*pascal/cell%deth
945 IF (debug_stress .AND. use_virial) THEN
946 sttot = virial%pv_virial
947 END IF
948
949 ! *** If LSD, then combine alpha density and beta density to
950 ! *** total density: alpha <- alpha + beta and
951 NULLIFY (mpa)
952 NULLIFY (matrix_ht)
953 IF (do_ex) THEN
954 CALL dbcsr_allocate_matrix_set(mpa, nspins)
955 DO ispin = 1, nspins
956 ALLOCATE (mpa(ispin)%matrix)
957 CALL dbcsr_create(mpa(ispin)%matrix, template=p_env%p1(ispin)%matrix)
958 CALL dbcsr_copy(mpa(ispin)%matrix, p_env%p1(ispin)%matrix)
959 CALL dbcsr_add(mpa(ispin)%matrix, ex_env%matrix_pe(ispin)%matrix, 1.0_dp, 1.0_dp)
960 CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
961 END DO
962 ELSE
963 mpa => matrix_pz
964 END IF
965 !
966 IF (do_ex .OR. (gapw .OR. gapw_xc)) THEN
967 CALL dbcsr_allocate_matrix_set(matrix_ht, nspins)
968 DO ispin = 1, nspins
969 ALLOCATE (matrix_ht(ispin)%matrix)
970 CALL dbcsr_create(matrix_ht(ispin)%matrix, template=matrix_hz(ispin)%matrix)
971 CALL dbcsr_copy(matrix_ht(ispin)%matrix, matrix_hz(ispin)%matrix)
972 CALL dbcsr_set(matrix_ht(ispin)%matrix, 0.0_dp)
973 END DO
974 END IF
975 !
976 ! START OF Tr[(P+Z)Hcore]
977 !
978
979 ! Kinetic energy matrix
980 NULLIFY (scrm2)
981 mpa2(1:nspins, 1:1) => mpa(1:nspins)
982 CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm2, matrix_p=mpa2, &
983 matrix_name="KINETIC ENERGY MATRIX", &
984 basis_type="ORB", &
985 sab_orb=sab_orb, calculate_forces=.true., &
986 debug_forces=debug_forces, debug_stress=debug_stress)
987 CALL dbcsr_deallocate_matrix_set(scrm2)
988
989 ! Initialize a matrix scrm, later used for scratch purposes
990 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
991 NULLIFY (scrm)
992 CALL dbcsr_allocate_matrix_set(scrm, nspins)
993 DO ispin = 1, nspins
994 ALLOCATE (scrm(ispin)%matrix)
995 CALL dbcsr_create(scrm(ispin)%matrix, template=matrix_s(1)%matrix)
996 CALL dbcsr_copy(scrm(ispin)%matrix, matrix_s(1)%matrix)
997 CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
998 END DO
999
1000 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
1001 atomic_kind_set=atomic_kind_set)
1002
1003 ALLOCATE (matrix_p(nspins, 1), matrix_h(nspins, 1))
1004 DO ispin = 1, nspins
1005 matrix_p(ispin, 1)%matrix => mpa(ispin)%matrix
1006 matrix_h(ispin, 1)%matrix => scrm(ispin)%matrix
1007 END DO
1008 matrix_h(1, 1)%matrix => scrm(1)%matrix
1009
1010 nder = 1
1011 CALL core_matrices(qs_env, matrix_h, matrix_p, .true., nder, &
1012 debug_forces=debug_forces, debug_stress=debug_stress)
1013
1014 ! Kim-Gordon subsystem DFT
1015 ! Atomic potential for nonadditive kinetic energy contribution
1016 IF (dft_control%qs_control%do_kg) THEN
1017 IF (qs_env%kg_env%tnadd_method == kg_tnadd_atomic) THEN
1018 CALL get_qs_env(qs_env=qs_env, kg_env=kg_env, dbcsr_dist=dbcsr_dist)
1019
1020 IF (use_virial) THEN
1021 pv_loc = virial%pv_virial
1022 END IF
1023
1024 IF (debug_forces) fodeb(1:3) = force(1)%kinetic(1:3, 1)
1025 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1026 CALL build_tnadd_mat(kg_env=kg_env, matrix_p=matrix_p, force=force, virial=virial, &
1027 calculate_forces=.true., use_virial=use_virial, &
1028 qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
1029 particle_set=particle_set, sab_orb=sab_orb, dbcsr_dist=dbcsr_dist)
1030 IF (debug_forces) THEN
1031 fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
1032 CALL para_env%sum(fodeb)
1033 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dTnadd ", fodeb
1034 END IF
1035 IF (debug_stress .AND. use_virial) THEN
1036 stdeb = fconv*(virial%pv_virial - stdeb)
1037 CALL para_env%sum(stdeb)
1038 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1039 'STRESS| Pz*dTnadd ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1040 END IF
1041
1042 ! Stress-tensor update components
1043 IF (use_virial) THEN
1044 virial%pv_ekinetic = virial%pv_ekinetic + (virial%pv_virial - pv_loc)
1045 END IF
1046
1047 END IF
1048 END IF
1049
1050 DEALLOCATE (matrix_h)
1051 DEALLOCATE (matrix_p)
1052 CALL dbcsr_deallocate_matrix_set(scrm)
1053
1054 ! initialize src matrix
1055 ! Necessary as build_kinetic_matrix will only allocate scrm(1)
1056 ! and not scrm(2) in open-shell case
1057 NULLIFY (scrm)
1058 CALL dbcsr_allocate_matrix_set(scrm, nspins)
1059 DO ispin = 1, nspins
1060 ALLOCATE (scrm(ispin)%matrix)
1061 CALL dbcsr_create(scrm(ispin)%matrix, template=matrix_pz(1)%matrix)
1062 CALL dbcsr_copy(scrm(ispin)%matrix, matrix_pz(ispin)%matrix)
1063 CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
1064 END DO
1065
1066 IF (debug_forces) THEN
1067 ALLOCATE (ftot2(3, natom))
1068 CALL total_qs_force(ftot2, force, atomic_kind_set)
1069 fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
1070 CALL para_env%sum(fodeb)
1071 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (T+Dz)*dHcore", fodeb
1072 END IF
1073 IF (debug_stress .AND. use_virial) THEN
1074 stdeb = fconv*(virial%pv_virial - sttot)
1075 CALL para_env%sum(stdeb)
1076 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1077 'STRESS| Stress Pz*dHcore ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1078 ! save current total viral, does not contain volume terms yet
1079 sttot2 = virial%pv_virial
1080 END IF
1081 !
1082 ! END OF Tr(P+Z)Hcore
1083 !
1084 !
1085 ! Vhxc (KS potentials calculated externally)
1086 CALL get_qs_env(qs_env, pw_env=pw_env)
1087 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
1088 !
1089 IF (dft_control%do_admm) THEN
1090 CALL get_qs_env(qs_env, admm_env=admm_env)
1091 xc_section => admm_env%xc_section_primary
1092 ELSE
1093 xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
1094 END IF
1095 xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1096 CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
1097 !
1098 IF (gapw .OR. gapw_xc) THEN
1099 NULLIFY (oce, sab_orb)
1100 CALL get_qs_env(qs_env=qs_env, oce=oce, sab_orb=sab_orb)
1101 ! set up local_rho_set for GS density
1102 NULLIFY (local_rho_set_gs)
1103 CALL get_qs_env(qs_env=qs_env, rho=rho)
1104 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1105 CALL local_rho_set_create(local_rho_set_gs)
1106 CALL allocate_rho_atom_internals(local_rho_set_gs%rho_atom_set, atomic_kind_set, &
1107 qs_kind_set, dft_control, para_env)
1108 CALL init_rho0(local_rho_set_gs, qs_env, dft_control%qs_control%gapw_control)
1109 CALL rho0_s_grid_create(pw_env, local_rho_set_gs%rho0_mpole)
1110 CALL calculate_rho_atom_coeff(qs_env, matrix_p(:, 1), local_rho_set_gs%rho_atom_set, &
1111 qs_kind_set, oce, sab_orb, para_env)
1112 CALL prepare_gapw_den(qs_env, local_rho_set_gs, do_rho0=gapw)
1113 ! set up local_rho_set for response density
1114 NULLIFY (local_rho_set_t)
1115 CALL local_rho_set_create(local_rho_set_t)
1116 CALL allocate_rho_atom_internals(local_rho_set_t%rho_atom_set, atomic_kind_set, &
1117 qs_kind_set, dft_control, para_env)
1118 CALL init_rho0(local_rho_set_t, qs_env, dft_control%qs_control%gapw_control, &
1119 zcore=0.0_dp)
1120 CALL rho0_s_grid_create(pw_env, local_rho_set_t%rho0_mpole)
1121 CALL calculate_rho_atom_coeff(qs_env, mpa(:), local_rho_set_t%rho_atom_set, &
1122 qs_kind_set, oce, sab_orb, para_env)
1123 CALL prepare_gapw_den(qs_env, local_rho_set_t, do_rho0=gapw)
1124
1125 ! compute soft GS potential
1126 ALLOCATE (rho_r_gs(nspins), rho_g_gs(nspins))
1127 DO ispin = 1, nspins
1128 CALL auxbas_pw_pool%create_pw(rho_r_gs(ispin))
1129 CALL auxbas_pw_pool%create_pw(rho_g_gs(ispin))
1130 END DO
1131 CALL auxbas_pw_pool%create_pw(rho_tot_gspace_gs)
1132 ! compute soft GS density
1133 total_rho_gs = 0.0_dp
1134 CALL pw_zero(rho_tot_gspace_gs)
1135 DO ispin = 1, nspins
1136 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_p(ispin, 1)%matrix, &
1137 rho=rho_r_gs(ispin), &
1138 rho_gspace=rho_g_gs(ispin), &
1139 soft_valid=(gapw .OR. gapw_xc), &
1140 total_rho=total_rho_gs(ispin))
1141 CALL pw_axpy(rho_g_gs(ispin), rho_tot_gspace_gs)
1142 END DO
1143 IF (gapw) THEN
1144 CALL get_qs_env(qs_env, natom=natom)
1145 ! add rho0 contributions to GS density (only for Coulomb) only for gapw
1146 CALL pw_axpy(local_rho_set_gs%rho0_mpole%rho0_s_gs, rho_tot_gspace_gs)
1147 IF (ASSOCIATED(local_rho_set_gs%rho0_mpole%rhoz_cneo_s_gs)) THEN
1148 CALL pw_axpy(local_rho_set_gs%rho0_mpole%rhoz_cneo_s_gs, rho_tot_gspace_gs)
1149 END IF
1150 IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
1151 CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
1152 CALL pw_axpy(rho_core, rho_tot_gspace_gs)
1153 END IF
1154 ! compute GS potential
1155 CALL auxbas_pw_pool%create_pw(v_hartree_gspace_gs)
1156 CALL auxbas_pw_pool%create_pw(v_hartree_rspace_gs)
1157 NULLIFY (hartree_local_gs)
1158 CALL hartree_local_create(hartree_local_gs)
1159 CALL init_coulomb_local(hartree_local_gs, natom)
1160 CALL pw_poisson_solve(poisson_env, rho_tot_gspace_gs, hartree_gs, v_hartree_gspace_gs)
1161 CALL pw_transfer(v_hartree_gspace_gs, v_hartree_rspace_gs)
1162 CALL pw_scale(v_hartree_rspace_gs, v_hartree_rspace_gs%pw_grid%dvol)
1163 END IF
1164 END IF
1165
1166 IF (gapw) THEN
1167 ! Hartree grid PAW term
1168 cpassert(.NOT. use_virial)
1169 IF (debug_forces) fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1)
1170 CALL vh_1c_gg_integrals(qs_env, hartree_gs, hartree_local_gs%ecoul_1c, local_rho_set_t, para_env, tddft=.true., &
1171 local_rho_set_2nd=local_rho_set_gs, core_2nd=.false.) ! n^core for GS potential
1172 ! 1st to define integral space, 2nd for potential, integral contributions stored on local_rho_set_gs
1173 CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace_gs, para_env, calculate_forces=.true., &
1174 local_rho_set=local_rho_set_t, local_rho_set_2nd=local_rho_set_gs)
1175 IF (debug_forces) THEN
1176 fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1) - fodeb(1:3)
1177 CALL para_env%sum(fodeb)
1178 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (T+Dz)*dVh[D^GS]PAWg0", fodeb
1179 END IF
1180 END IF
1181 IF (gapw .OR. gapw_xc) THEN
1182 IF (myfun /= xc_none) THEN
1183 ! add 1c hard and soft XC contributions
1184 NULLIFY (local_rho_set_vxc)
1185 CALL local_rho_set_create(local_rho_set_vxc)
1186 CALL allocate_rho_atom_internals(local_rho_set_vxc%rho_atom_set, atomic_kind_set, &
1187 qs_kind_set, dft_control, para_env)
1188 CALL calculate_rho_atom_coeff(qs_env, matrix_p(:, 1), local_rho_set_vxc%rho_atom_set, &
1189 qs_kind_set, oce, sab_orb, para_env)
1190 CALL prepare_gapw_den(qs_env, local_rho_set_vxc, do_rho0=.false.)
1191 ! compute hard and soft atomic contributions
1192 CALL calculate_vxc_atom(qs_env, .false., exc1=hartree_gs, xc_section_external=xc_section, &
1193 rho_atom_set_external=local_rho_set_vxc%rho_atom_set)
1194 END IF ! myfun
1195 END IF ! gapw
1196
1197 CALL auxbas_pw_pool%create_pw(vhxc_rspace)
1198 !
1199 ! Stress-tensor: integration contribution direct term
1200 ! int v_Hxc[n^in]*n^z
1201 IF (use_virial) THEN
1202 pv_loc = virial%pv_virial
1203 END IF
1204
1205 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1206 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1207 IF (gapw .OR. gapw_xc) THEN
1208 ! vtot = v_xc + v_hartree
1209 DO ispin = 1, nspins
1210 CALL pw_zero(vhxc_rspace)
1211 IF (gapw) THEN
1212 CALL pw_transfer(v_hartree_rspace_gs, vhxc_rspace)
1213 ELSEIF (gapw_xc) THEN
1214 CALL pw_transfer(vh_rspace, vhxc_rspace)
1215 END IF
1216 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1217 hmat=scrm(ispin), pmat=mpa(ispin), &
1218 qs_env=qs_env, gapw=gapw, &
1219 calculate_forces=.true.)
1220 END DO
1221 IF (myfun /= xc_none) THEN
1222 DO ispin = 1, nspins
1223 CALL pw_zero(vhxc_rspace)
1224 CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
1225 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1226 hmat=scrm(ispin), pmat=mpa(ispin), &
1227 qs_env=qs_env, gapw=(gapw .OR. gapw_xc), &
1228 calculate_forces=.true.)
1229 END DO
1230 END IF
1231 ELSE ! original GPW with Standard Hartree as Potential
1232 DO ispin = 1, nspins
1233 CALL pw_transfer(vh_rspace, vhxc_rspace)
1234 CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
1235 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1236 hmat=scrm(ispin), pmat=mpa(ispin), &
1237 qs_env=qs_env, gapw=gapw, calculate_forces=.true.)
1238 END DO
1239 END IF
1240
1241 IF (debug_forces) THEN
1242 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1243 CALL para_env%sum(fodeb)
1244 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (T+Dz)*dVhxc[D^GS] ", fodeb
1245 END IF
1246 IF (debug_stress .AND. use_virial) THEN
1247 stdeb = fconv*(virial%pv_virial - pv_loc)
1248 CALL para_env%sum(stdeb)
1249 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1250 'STRESS| INT Pz*dVhxc ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1251 END IF
1252
1253 IF (gapw .OR. gapw_xc) THEN
1254 ! HXC term
1255 IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
1256 IF (gapw) CALL update_ks_atom(qs_env, scrm, mpa, forces=.true., tddft=.false., &
1257 rho_atom_external=local_rho_set_gs%rho_atom_set)
1258 IF (myfun /= xc_none) CALL update_ks_atom(qs_env, scrm, mpa, forces=.true., tddft=.false., &
1259 rho_atom_external=local_rho_set_vxc%rho_atom_set)
1260 IF (debug_forces) THEN
1261 fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
1262 CALL para_env%sum(fodeb)
1263 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (T+Dz)*dVhxc[D^GS]PAW ", fodeb
1264 END IF
1265 ! release local environments for GAPW
1266 IF (myfun /= xc_none) THEN
1267 IF (ASSOCIATED(local_rho_set_vxc)) CALL local_rho_set_release(local_rho_set_vxc)
1268 END IF
1269 IF (ASSOCIATED(local_rho_set_gs)) CALL local_rho_set_release(local_rho_set_gs)
1270 IF (gapw) THEN
1271 IF (ASSOCIATED(hartree_local_gs)) CALL hartree_local_release(hartree_local_gs)
1272 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace_gs)
1273 CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace_gs)
1274 END IF
1275 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace_gs)
1276 IF (ASSOCIATED(rho_r_gs)) THEN
1277 DO ispin = 1, nspins
1278 CALL auxbas_pw_pool%give_back_pw(rho_r_gs(ispin))
1279 END DO
1280 DEALLOCATE (rho_r_gs)
1281 END IF
1282 IF (ASSOCIATED(rho_g_gs)) THEN
1283 DO ispin = 1, nspins
1284 CALL auxbas_pw_pool%give_back_pw(rho_g_gs(ispin))
1285 END DO
1286 DEALLOCATE (rho_g_gs)
1287 END IF
1288 END IF !gapw
1289
1290 IF (ASSOCIATED(vtau_rspace)) THEN
1291 cpassert(.NOT. (gapw .OR. gapw_xc))
1292 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1293 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1294 DO ispin = 1, nspins
1295 CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
1296 hmat=scrm(ispin), pmat=mpa(ispin), &
1297 qs_env=qs_env, gapw=(gapw .OR. gapw_xc), &
1298 calculate_forces=.true., compute_tau=.true.)
1299 END DO
1300 IF (debug_forces) THEN
1301 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1302 CALL para_env%sum(fodeb)
1303 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dVxc_tau ", fodeb
1304 END IF
1305 IF (debug_stress .AND. use_virial) THEN
1306 stdeb = fconv*(virial%pv_virial - pv_loc)
1307 CALL para_env%sum(stdeb)
1308 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1309 'STRESS| INT Pz*dVxc_tau ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1310 END IF
1311 END IF
1312 CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
1313
1314 ! Stress-tensor Pz*v_Hxc[Pin]
1315 IF (use_virial) THEN
1316 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1317 END IF
1318
1319 ! KG Embedding
1320 ! calculate kinetic energy potential and integrate with response density
1321 IF (dft_control%qs_control%do_kg) THEN
1322 IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed .OR. &
1323 qs_env%kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
1324
1325 ekin_mol = 0.0_dp
1326 IF (use_virial) THEN
1327 pv_loc = virial%pv_virial
1328 END IF
1329
1330 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1331 CALL kg_ekin_subset(qs_env=qs_env, &
1332 ks_matrix=scrm, &
1333 ekin_mol=ekin_mol, &
1334 calc_force=.true., &
1335 do_kernel=.false., &
1336 pmat_ext=mpa)
1337 IF (debug_forces) THEN
1338 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1339 CALL para_env%sum(fodeb)
1340 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dVkg ", fodeb
1341 END IF
1342 IF (debug_stress .AND. use_virial) THEN
1343 !IF (iounit > 0) WRITE(iounit, *) &
1344 ! "response_force | VOL 1st KG - v_KG[n_in]*n_z: ", ekin_mol
1345 stdeb = 1.0_dp*fconv*ekin_mol
1346 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1347 'STRESS| VOL KG Pz*dVKG ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1348
1349 stdeb = fconv*(virial%pv_virial - pv_loc)
1350 CALL para_env%sum(stdeb)
1351 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1352 'STRESS| INT KG Pz*dVKG ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1353
1354 stdeb = fconv*virial%pv_xc
1355 CALL para_env%sum(stdeb)
1356 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1357 'STRESS| GGA KG Pz*dVKG ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1358 END IF
1359 IF (use_virial) THEN
1360 ! Direct integral contribution
1361 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1362 END IF
1363
1364 END IF ! tnadd_method
1365 END IF ! do_kg
1366
1367 CALL dbcsr_deallocate_matrix_set(scrm)
1368
1369 !
1370 ! Hartree potential of response density
1371 !
1372 ALLOCATE (rhoz_r(nspins), rhoz_g(nspins))
1373 DO ispin = 1, nspins
1374 CALL auxbas_pw_pool%create_pw(rhoz_r(ispin))
1375 CALL auxbas_pw_pool%create_pw(rhoz_g(ispin))
1376 END DO
1377 CALL auxbas_pw_pool%create_pw(rhoz_tot_gspace)
1378 CALL auxbas_pw_pool%create_pw(zv_hartree_rspace)
1379 CALL auxbas_pw_pool%create_pw(zv_hartree_gspace)
1380
1381 CALL pw_zero(rhoz_tot_gspace)
1382 DO ispin = 1, nspins
1383 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1384 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin), &
1385 soft_valid=gapw)
1386 CALL pw_axpy(rhoz_g(ispin), rhoz_tot_gspace)
1387 END DO
1388 IF (gapw_xc) THEN
1389 NULLIFY (tauz_r_xc)
1390 ALLOCATE (rhoz_r_xc(nspins), rhoz_g_xc(nspins))
1391 DO ispin = 1, nspins
1392 CALL auxbas_pw_pool%create_pw(rhoz_r_xc(ispin))
1393 CALL auxbas_pw_pool%create_pw(rhoz_g_xc(ispin))
1394 END DO
1395 DO ispin = 1, nspins
1396 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1397 rho=rhoz_r_xc(ispin), rho_gspace=rhoz_g_xc(ispin), &
1398 soft_valid=gapw_xc)
1399 END DO
1400 END IF
1401
1402 IF (ASSOCIATED(vtau_rspace)) THEN
1403 cpassert(.NOT. (gapw .OR. gapw_xc))
1404 block
1405 TYPE(pw_c1d_gs_type) :: work_g
1406 ALLOCATE (tauz_r(nspins))
1407 CALL auxbas_pw_pool%create_pw(work_g)
1408 DO ispin = 1, nspins
1409 CALL auxbas_pw_pool%create_pw(tauz_r(ispin))
1410 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1411 rho=tauz_r(ispin), rho_gspace=work_g, &
1412 compute_tau=.true.)
1413 END DO
1414 CALL auxbas_pw_pool%give_back_pw(work_g)
1415 END block
1416 END IF
1417
1418 !
1419 IF (PRESENT(rhopz_r)) THEN
1420 DO ispin = 1, nspins
1421 CALL pw_copy(rhoz_r(ispin), rhopz_r(ispin))
1422 END DO
1423 END IF
1424
1425 ! Stress-tensor contribution second derivative
1426 ! Volume : int v_H[n^z]*n_in
1427 ! Volume : int epsilon_xc*n_z
1428 IF (use_virial) THEN
1429
1430 CALL get_qs_env(qs_env, rho=rho)
1431 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
1432
1433 ! Get the total input density in g-space [ions + electrons]
1434 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
1435
1436 h_stress(:, :) = 0.0_dp
1437 ! calculate associated hartree potential
1438 ! This term appears twice in the derivation of the equations
1439 ! v_H[n_in]*n_z and v_H[n_z]*n_in
1440 ! due to symmetry we only need to call this routine once,
1441 ! and count the Volume and Green function contribution
1442 ! which is stored in h_stress twice
1443 CALL pw_poisson_solve(poisson_env, &
1444 density=rhoz_tot_gspace, & ! n_z
1445 ehartree=ehartree, &
1446 vhartree=zv_hartree_gspace, & ! v_H[n_z]
1447 h_stress=h_stress, &
1448 aux_density=rho_tot_gspace) ! n_in
1449
1450 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1451
1452 ! Stress tensor Green function contribution
1453 virial%pv_ehartree = virial%pv_ehartree + 2.0_dp*h_stress/real(para_env%num_pe, dp)
1454 virial%pv_virial = virial%pv_virial + 2.0_dp*h_stress/real(para_env%num_pe, dp)
1455
1456 IF (debug_stress) THEN
1457 stdeb = -1.0_dp*fconv*ehartree
1458 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1459 'STRESS| VOL 1st v_H[n_z]*n_in ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1460 stdeb = -1.0_dp*fconv*ehartree
1461 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1462 'STRESS| VOL 2nd v_H[n_in]*n_z ', 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 1st v_H[n_z]*n_in ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1467 stdeb = fconv*(h_stress/real(para_env%num_pe, dp))
1468 CALL para_env%sum(stdeb)
1469 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1470 'STRESS| GREEN 2nd v_H[n_in]*n_z ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1471 END IF
1472
1473 ! Stress tensor volume term: \int v_xc[n_in]*n_z
1474 ! vxc_rspace already scaled, we need to unscale it!
1475 exc = 0.0_dp
1476 DO ispin = 1, nspins
1477 exc = exc + pw_integral_ab(rhoz_r(ispin), vxc_rspace(ispin))/ &
1478 vxc_rspace(ispin)%pw_grid%dvol
1479 END DO
1480 IF (ASSOCIATED(vtau_rspace)) THEN
1481 DO ispin = 1, nspins
1482 exc = exc + pw_integral_ab(tauz_r(ispin), vtau_rspace(ispin))/ &
1483 vtau_rspace(ispin)%pw_grid%dvol
1484 END DO
1485 END IF
1486
1487 ! Add KG embedding correction
1488 IF (dft_control%qs_control%do_kg) THEN
1489 IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed .OR. &
1490 qs_env%kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
1491 exc = exc - ekin_mol
1492 END IF
1493 END IF
1494
1495 IF (debug_stress) THEN
1496 stdeb = -1.0_dp*fconv*exc
1497 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1498 'STRESS| VOL 1st eps_XC[n_in]*n_z', one_third_sum_diag(stdeb), det_3x3(stdeb)
1499 END IF
1500
1501 ELSE ! use_virial
1502
1503 ! calculate associated hartree potential
1504 ! contribution for both T and D^Z
1505 IF (gapw) THEN
1506 CALL pw_axpy(local_rho_set_t%rho0_mpole%rho0_s_gs, rhoz_tot_gspace)
1507 IF (ASSOCIATED(local_rho_set_t%rho0_mpole%rhoz_cneo_s_gs)) THEN
1508 CALL pw_axpy(local_rho_set_t%rho0_mpole%rhoz_cneo_s_gs, rhoz_tot_gspace)
1509 END IF
1510 END IF
1511 CALL pw_poisson_solve(poisson_env, rhoz_tot_gspace, ehartree, zv_hartree_gspace)
1512
1513 END IF ! use virial
1514 IF (gapw .OR. gapw_xc) THEN
1515 IF (ASSOCIATED(local_rho_set_t)) CALL local_rho_set_release(local_rho_set_t)
1516 END IF
1517
1518 IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
1519 IF (debug_stress .AND. use_virial) stdeb = virial%pv_ehartree
1520 CALL pw_transfer(zv_hartree_gspace, zv_hartree_rspace)
1521 CALL pw_scale(zv_hartree_rspace, zv_hartree_rspace%pw_grid%dvol)
1522 ! Getting nuclear force contribution from the core charge density (not for GAPW)
1523 CALL integrate_v_core_rspace(zv_hartree_rspace, qs_env)
1524 IF (debug_forces) THEN
1525 fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
1526 CALL para_env%sum(fodeb)
1527 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vh(rhoz)*dncore ", fodeb
1528 END IF
1529 IF (debug_stress .AND. use_virial) THEN
1530 stdeb = fconv*(virial%pv_ehartree - stdeb)
1531 CALL para_env%sum(stdeb)
1532 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1533 'STRESS| INT Vh(rhoz)*dncore ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1534 END IF
1535
1536 !
1537 IF (gapw_xc) THEN
1538 CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
1539 ELSE
1540 CALL get_qs_env(qs_env=qs_env, rho=rho)
1541 END IF
1542 IF (dft_control%do_admm) THEN
1543 CALL get_qs_env(qs_env, admm_env=admm_env)
1544 xc_section => admm_env%xc_section_primary
1545 ELSE
1546 xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
1547 END IF
1548
1549 IF (use_virial) THEN
1550 virial%pv_xc = 0.0_dp
1551 END IF
1552 !
1553 NULLIFY (v_xc, v_xc_tau)
1554 IF (gapw_xc) THEN
1555 CALL create_kernel(qs_env, vxc=v_xc, vxc_tau=v_xc_tau, &
1556 rho=rho_xc, rho1_r=rhoz_r_xc, rho1_g=rhoz_g_xc, tau1_r=tauz_r_xc, &
1557 xc_section=xc_section, compute_virial=use_virial, virial_xc=virial%pv_xc)
1558 ELSE
1559 CALL create_kernel(qs_env, vxc=v_xc, vxc_tau=v_xc_tau, &
1560 rho=rho, rho1_r=rhoz_r, rho1_g=rhoz_g, tau1_r=tauz_r, &
1561 xc_section=xc_section, compute_virial=use_virial, virial_xc=virial%pv_xc)
1562 END IF
1563
1564 IF (gapw .OR. gapw_xc) THEN
1565 !get local_rho_set for GS density and response potential / density
1566 NULLIFY (local_rho_set_t)
1567 CALL local_rho_set_create(local_rho_set_t)
1568 CALL allocate_rho_atom_internals(local_rho_set_t%rho_atom_set, atomic_kind_set, &
1569 qs_kind_set, dft_control, para_env)
1570 CALL init_rho0(local_rho_set_t, qs_env, dft_control%qs_control%gapw_control, &
1571 zcore=0.0_dp)
1572 CALL rho0_s_grid_create(pw_env, local_rho_set_t%rho0_mpole)
1573 CALL calculate_rho_atom_coeff(qs_env, mpa(:), local_rho_set_t%rho_atom_set, &
1574 qs_kind_set, oce, sab_orb, para_env)
1575 CALL prepare_gapw_den(qs_env, local_rho_set_t, do_rho0=gapw)
1576 NULLIFY (local_rho_set_gs)
1577 CALL local_rho_set_create(local_rho_set_gs)
1578 CALL allocate_rho_atom_internals(local_rho_set_gs%rho_atom_set, atomic_kind_set, &
1579 qs_kind_set, dft_control, para_env)
1580 CALL init_rho0(local_rho_set_gs, qs_env, dft_control%qs_control%gapw_control)
1581 CALL rho0_s_grid_create(pw_env, local_rho_set_gs%rho0_mpole)
1582 CALL calculate_rho_atom_coeff(qs_env, matrix_p(:, 1), local_rho_set_gs%rho_atom_set, &
1583 qs_kind_set, oce, sab_orb, para_env)
1584 CALL prepare_gapw_den(qs_env, local_rho_set_gs, do_rho0=gapw)
1585 ! compute response potential
1586 ALLOCATE (rho_r_t(nspins), rho_g_t(nspins))
1587 DO ispin = 1, nspins
1588 CALL auxbas_pw_pool%create_pw(rho_r_t(ispin))
1589 CALL auxbas_pw_pool%create_pw(rho_g_t(ispin))
1590 END DO
1591 CALL auxbas_pw_pool%create_pw(rho_tot_gspace_t)
1592 total_rho_t = 0.0_dp
1593 CALL pw_zero(rho_tot_gspace_t)
1594 DO ispin = 1, nspins
1595 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1596 rho=rho_r_t(ispin), &
1597 rho_gspace=rho_g_t(ispin), &
1598 soft_valid=gapw, &
1599 total_rho=total_rho_t(ispin))
1600 CALL pw_axpy(rho_g_t(ispin), rho_tot_gspace_t)
1601 END DO
1602 ! add rho0 contributions to response density (only for Coulomb) only for gapw
1603 IF (gapw) THEN
1604 CALL pw_axpy(local_rho_set_t%rho0_mpole%rho0_s_gs, rho_tot_gspace_t)
1605 IF (ASSOCIATED(local_rho_set_t%rho0_mpole%rhoz_cneo_s_gs)) THEN
1606 CALL pw_axpy(local_rho_set_t%rho0_mpole%rhoz_cneo_s_gs, rho_tot_gspace_t)
1607 END IF
1608 ! compute response Coulomb potential
1609 CALL auxbas_pw_pool%create_pw(v_hartree_gspace_t)
1610 CALL auxbas_pw_pool%create_pw(v_hartree_rspace_t)
1611 NULLIFY (hartree_local_t)
1612 CALL hartree_local_create(hartree_local_t)
1613 CALL init_coulomb_local(hartree_local_t, natom)
1614 CALL pw_poisson_solve(poisson_env, rho_tot_gspace_t, hartree_t, v_hartree_gspace_t)
1615 CALL pw_transfer(v_hartree_gspace_t, v_hartree_rspace_t)
1616 CALL pw_scale(v_hartree_rspace_t, v_hartree_rspace_t%pw_grid%dvol)
1617 !
1618 IF (debug_forces) fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1)
1619 CALL vh_1c_gg_integrals(qs_env, hartree_t, hartree_local_t%ecoul_1c, local_rho_set_gs, para_env, tddft=.false., &
1620 local_rho_set_2nd=local_rho_set_t, core_2nd=.true.) ! n^core for GS potential
1621 CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace_t, para_env, calculate_forces=.true., &
1622 local_rho_set=local_rho_set_gs, local_rho_set_2nd=local_rho_set_t)
1623 IF (debug_forces) THEN
1624 fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1) - fodeb(1:3)
1625 CALL para_env%sum(fodeb)
1626 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vh(T)*dncore PAWg0", fodeb
1627 END IF
1628 END IF !gapw
1629 END IF !gapw
1630
1631 IF (gapw .OR. gapw_xc) THEN
1632 !GAPW compute atomic fxc contributions
1633 IF (myfun /= xc_none) THEN
1634 ! local_rho_set_f
1635 NULLIFY (local_rho_set_f)
1636 CALL local_rho_set_create(local_rho_set_f)
1637 CALL allocate_rho_atom_internals(local_rho_set_f%rho_atom_set, atomic_kind_set, &
1638 qs_kind_set, dft_control, para_env)
1639 CALL calculate_rho_atom_coeff(qs_env, mpa, local_rho_set_f%rho_atom_set, &
1640 qs_kind_set, oce, sab_orb, para_env)
1641 CALL prepare_gapw_den(qs_env, local_rho_set_f, do_rho0=.false.)
1642 ! add hard and soft atomic contributions
1643 CALL calculate_xc_2nd_deriv_atom(local_rho_set_gs%rho_atom_set, &
1644 local_rho_set_f%rho_atom_set, &
1645 qs_env, xc_section, para_env, &
1646 do_triplet=.false.)
1647 END IF ! myfun
1648 END IF
1649
1650 ! Stress-tensor XC-kernel GGA contribution
1651 IF (use_virial) THEN
1652 virial%pv_exc = virial%pv_exc + virial%pv_xc
1653 virial%pv_virial = virial%pv_virial + virial%pv_xc
1654 END IF
1655
1656 IF (debug_stress .AND. use_virial) THEN
1657 stdeb = 1.0_dp*fconv*virial%pv_xc
1658 CALL para_env%sum(stdeb)
1659 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1660 'STRESS| GGA 2nd Pin*dK*rhoz', one_third_sum_diag(stdeb), det_3x3(stdeb)
1661 END IF
1662
1663 ! Stress-tensor integral contribution of 2nd derivative terms
1664 IF (use_virial) THEN
1665 pv_loc = virial%pv_virial
1666 END IF
1667
1668 CALL get_qs_env(qs_env=qs_env, rho=rho)
1669 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1670 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1671
1672 DO ispin = 1, nspins
1673 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
1674 END DO
1675 IF ((.NOT. (gapw)) .AND. (.NOT. gapw_xc)) THEN
1676 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1677 DO ispin = 1, nspins
1678 CALL pw_axpy(zv_hartree_rspace, v_xc(ispin)) ! Hartree potential of response density
1679 CALL integrate_v_rspace(qs_env=qs_env, &
1680 v_rspace=v_xc(ispin), &
1681 hmat=matrix_hz(ispin), &
1682 pmat=matrix_p(ispin, 1), &
1683 gapw=.false., &
1684 calculate_forces=.true.)
1685 END DO
1686 IF (debug_forces) THEN
1687 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1688 CALL para_env%sum(fodeb)
1689 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKhxc*rhoz ", fodeb
1690 END IF
1691 ELSE
1692 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1693 IF (myfun /= xc_none) THEN
1694 DO ispin = 1, nspins
1695 CALL integrate_v_rspace(qs_env=qs_env, &
1696 v_rspace=v_xc(ispin), &
1697 hmat=matrix_hz(ispin), &
1698 pmat=matrix_p(ispin, 1), &
1699 gapw=.true., &
1700 calculate_forces=.true.)
1701 END DO
1702 END IF ! my_fun
1703 ! Coulomb T+Dz
1704 DO ispin = 1, nspins
1705 CALL pw_zero(v_xc(ispin))
1706 IF (gapw) THEN ! Hartree potential of response density
1707 CALL pw_axpy(v_hartree_rspace_t, v_xc(ispin))
1708 ELSEIF (gapw_xc) THEN
1709 CALL pw_axpy(zv_hartree_rspace, v_xc(ispin))
1710 END IF
1711 CALL integrate_v_rspace(qs_env=qs_env, &
1712 v_rspace=v_xc(ispin), &
1713 hmat=matrix_ht(ispin), &
1714 pmat=matrix_p(ispin, 1), &
1715 gapw=gapw, &
1716 calculate_forces=.true.)
1717 END DO
1718 IF (debug_forces) THEN
1719 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1720 CALL para_env%sum(fodeb)
1721 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKhxc*rhoz ", fodeb
1722 END IF
1723 END IF
1724
1725 IF (gapw .OR. gapw_xc) THEN
1726 ! compute hard and soft atomic contributions
1727 IF (myfun /= xc_none) THEN
1728 IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
1729 CALL update_ks_atom(qs_env, matrix_hz, matrix_p, forces=.true., tddft=.false., &
1730 rho_atom_external=local_rho_set_f%rho_atom_set)
1731 IF (debug_forces) THEN
1732 fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
1733 CALL para_env%sum(fodeb)
1734 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P^GS*dKxc*(Dz+T) PAW", fodeb
1735 END IF
1736 END IF !myfun
1737 ! Coulomb contributions
1738 IF (gapw) THEN
1739 IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
1740 CALL update_ks_atom(qs_env, matrix_ht, matrix_p, forces=.true., tddft=.false., &
1741 rho_atom_external=local_rho_set_t%rho_atom_set)
1742 IF (debug_forces) THEN
1743 fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
1744 CALL para_env%sum(fodeb)
1745 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P^GS*dKh*(Dz+T) PAW", fodeb
1746 END IF
1747 END IF
1748 ! add Coulomb and XC
1749 DO ispin = 1, nspins
1750 CALL dbcsr_add(matrix_hz(ispin)%matrix, matrix_ht(ispin)%matrix, 1.0_dp, 1.0_dp)
1751 END DO
1752
1753 ! release
1754 IF (myfun /= xc_none) THEN
1755 IF (ASSOCIATED(local_rho_set_f)) CALL local_rho_set_release(local_rho_set_f)
1756 END IF
1757 IF (ASSOCIATED(local_rho_set_t)) CALL local_rho_set_release(local_rho_set_t)
1758 IF (ASSOCIATED(local_rho_set_gs)) CALL local_rho_set_release(local_rho_set_gs)
1759 IF (gapw) THEN
1760 IF (ASSOCIATED(hartree_local_t)) CALL hartree_local_release(hartree_local_t)
1761 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace_t)
1762 CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace_t)
1763 END IF
1764 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace_t)
1765 DO ispin = 1, nspins
1766 CALL auxbas_pw_pool%give_back_pw(rho_r_t(ispin))
1767 CALL auxbas_pw_pool%give_back_pw(rho_g_t(ispin))
1768 END DO
1769 DEALLOCATE (rho_r_t, rho_g_t)
1770 END IF ! gapw
1771
1772 IF (debug_stress .AND. use_virial) THEN
1773 stdeb = fconv*(virial%pv_virial - stdeb)
1774 CALL para_env%sum(stdeb)
1775 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1776 'STRESS| INT 2nd f_Hxc[Pz]*Pin', one_third_sum_diag(stdeb), det_3x3(stdeb)
1777 END IF
1778 !
1779 IF (ASSOCIATED(v_xc_tau)) THEN
1780 cpassert(.NOT. (gapw .OR. gapw_xc))
1781 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1782 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1783 DO ispin = 1, nspins
1784 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
1785 CALL integrate_v_rspace(qs_env=qs_env, &
1786 v_rspace=v_xc_tau(ispin), &
1787 hmat=matrix_hz(ispin), &
1788 pmat=matrix_p(ispin, 1), &
1789 compute_tau=.true., &
1790 gapw=(gapw .OR. gapw_xc), &
1791 calculate_forces=.true.)
1792 END DO
1793 IF (debug_forces) THEN
1794 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1795 CALL para_env%sum(fodeb)
1796 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKtau*tauz ", fodeb
1797 END IF
1798 END IF
1799 IF (debug_stress .AND. use_virial) THEN
1800 stdeb = fconv*(virial%pv_virial - stdeb)
1801 CALL para_env%sum(stdeb)
1802 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1803 'STRESS| INT 2nd f_xctau[Pz]*Pin', one_third_sum_diag(stdeb), det_3x3(stdeb)
1804 END IF
1805 ! Stress-tensor integral contribution of 2nd derivative terms
1806 IF (use_virial) THEN
1807 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1808 END IF
1809
1810 ! KG Embedding
1811 ! calculate kinetic energy kernel, folded with response density for partial integration
1812 IF (dft_control%qs_control%do_kg) THEN
1813 IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed) THEN
1814 ekin_mol = 0.0_dp
1815 IF (use_virial) THEN
1816 pv_loc = virial%pv_virial
1817 END IF
1818
1819 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1820 IF (use_virial) virial%pv_xc = 0.0_dp
1821 CALL kg_ekin_subset(qs_env=qs_env, &
1822 ks_matrix=matrix_hz, &
1823 ekin_mol=ekin_mol, &
1824 calc_force=.true., &
1825 do_kernel=.true., &
1826 pmat_ext=matrix_pz)
1827
1828 IF (debug_forces) THEN
1829 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1830 CALL para_env%sum(fodeb)
1831 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*d(Kkg)*rhoz ", fodeb
1832 END IF
1833 IF (debug_stress .AND. use_virial) THEN
1834 stdeb = fconv*(virial%pv_virial - pv_loc)
1835 CALL para_env%sum(stdeb)
1836 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1837 'STRESS| INT KG Pin*d(KKG)*rhoz ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1838
1839 stdeb = fconv*(virial%pv_xc)
1840 CALL para_env%sum(stdeb)
1841 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1842 'STRESS| GGA KG Pin*d(KKG)*rhoz ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1843 END IF
1844
1845 ! Stress tensor
1846 IF (use_virial) THEN
1847 ! XC-kernel Integral contribution
1848 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1849
1850 ! XC-kernel GGA contribution
1851 virial%pv_exc = virial%pv_exc - virial%pv_xc
1852 virial%pv_virial = virial%pv_virial - virial%pv_xc
1853 virial%pv_xc = 0.0_dp
1854 END IF
1855 END IF
1856 END IF
1857 CALL auxbas_pw_pool%give_back_pw(rhoz_tot_gspace)
1858 CALL auxbas_pw_pool%give_back_pw(zv_hartree_gspace)
1859 CALL auxbas_pw_pool%give_back_pw(zv_hartree_rspace)
1860 DO ispin = 1, nspins
1861 CALL auxbas_pw_pool%give_back_pw(rhoz_r(ispin))
1862 CALL auxbas_pw_pool%give_back_pw(rhoz_g(ispin))
1863 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
1864 END DO
1865 DEALLOCATE (rhoz_r, rhoz_g, v_xc)
1866 IF (gapw_xc) THEN
1867 DO ispin = 1, nspins
1868 CALL auxbas_pw_pool%give_back_pw(rhoz_r_xc(ispin))
1869 CALL auxbas_pw_pool%give_back_pw(rhoz_g_xc(ispin))
1870 END DO
1871 DEALLOCATE (rhoz_r_xc, rhoz_g_xc)
1872 END IF
1873 IF (ASSOCIATED(v_xc_tau)) THEN
1874 DO ispin = 1, nspins
1875 CALL auxbas_pw_pool%give_back_pw(tauz_r(ispin))
1876 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
1877 END DO
1878 DEALLOCATE (tauz_r, v_xc_tau)
1879 END IF
1880 IF (debug_forces) THEN
1881 ALLOCATE (ftot3(3, natom))
1882 CALL total_qs_force(ftot3, force, atomic_kind_set)
1883 fodeb(1:3) = ftot3(1:3, 1) - ftot2(1:3, 1)
1884 CALL para_env%sum(fodeb)
1885 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*V(rhoz)", fodeb
1886 END IF
1887 CALL dbcsr_deallocate_matrix_set(scrm)
1888 CALL dbcsr_deallocate_matrix_set(matrix_ht)
1889
1890 ! -----------------------------------------
1891 ! Apply ADMM exchange correction
1892 ! -----------------------------------------
1893
1894 IF (dft_control%do_admm) THEN
1895 ! volume term
1896 exc_aux_fit = 0.0_dp
1897
1898 IF (qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
1899 ! nothing to do
1900 NULLIFY (mpz, mhz, mhx, mhy)
1901 ELSE
1902 ! add ADMM xc_section_aux terms: Pz*Vxc + P0*K0[rhoz]
1903 CALL get_qs_env(qs_env, admm_env=admm_env)
1904 CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit, matrix_s_aux_fit=scrm, &
1905 task_list_aux_fit=task_list_aux_fit)
1906 !
1907 NULLIFY (mpz, mhz, mhx, mhy)
1908 CALL dbcsr_allocate_matrix_set(mhx, nspins, 1)
1909 CALL dbcsr_allocate_matrix_set(mhy, nspins, 1)
1910 CALL dbcsr_allocate_matrix_set(mpz, nspins, 1)
1911 DO ispin = 1, nspins
1912 ALLOCATE (mhx(ispin, 1)%matrix)
1913 CALL dbcsr_create(mhx(ispin, 1)%matrix, template=scrm(1)%matrix)
1914 CALL dbcsr_copy(mhx(ispin, 1)%matrix, scrm(1)%matrix)
1915 CALL dbcsr_set(mhx(ispin, 1)%matrix, 0.0_dp)
1916 ALLOCATE (mhy(ispin, 1)%matrix)
1917 CALL dbcsr_create(mhy(ispin, 1)%matrix, template=scrm(1)%matrix)
1918 CALL dbcsr_copy(mhy(ispin, 1)%matrix, scrm(1)%matrix)
1919 CALL dbcsr_set(mhy(ispin, 1)%matrix, 0.0_dp)
1920 ALLOCATE (mpz(ispin, 1)%matrix)
1921 IF (do_ex) THEN
1922 CALL dbcsr_create(mpz(ispin, 1)%matrix, template=p_env%p1_admm(ispin)%matrix)
1923 CALL dbcsr_copy(mpz(ispin, 1)%matrix, p_env%p1_admm(ispin)%matrix)
1924 CALL dbcsr_add(mpz(ispin, 1)%matrix, ex_env%matrix_pe_admm(ispin)%matrix, &
1925 1.0_dp, 1.0_dp)
1926 ELSE
1927 CALL dbcsr_create(mpz(ispin, 1)%matrix, template=matrix_pz_admm(ispin)%matrix)
1928 CALL dbcsr_copy(mpz(ispin, 1)%matrix, matrix_pz_admm(ispin)%matrix)
1929 END IF
1930 END DO
1931 !
1932 xc_section => admm_env%xc_section_aux
1933 ! Stress-tensor: integration contribution direct term
1934 ! int Pz*v_xc[rho_admm]
1935 IF (use_virial) THEN
1936 pv_loc = virial%pv_virial
1937 END IF
1938
1939 basis_type = "AUX_FIT"
1940 task_list => task_list_aux_fit
1941 IF (admm_env%do_gapw) THEN
1942 basis_type = "AUX_FIT_SOFT"
1943 task_list => admm_env%admm_gapw_env%task_list
1944 END IF
1945 !
1946 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1947 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1948 DO ispin = 1, nspins
1949 CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), &
1950 hmat=mhx(ispin, 1), pmat=mpz(ispin, 1), &
1951 qs_env=qs_env, calculate_forces=.true., &
1952 basis_type=basis_type, task_list_external=task_list)
1953 END DO
1954 IF (debug_forces) THEN
1955 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1956 CALL para_env%sum(fodeb)
1957 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*Vxc(rho_admm)", fodeb
1958 END IF
1959 IF (debug_stress .AND. use_virial) THEN
1960 stdeb = fconv*(virial%pv_virial - pv_loc)
1961 CALL para_env%sum(stdeb)
1962 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
1963 'STRESS| INT 1st Pz*dVxc(rho_admm) ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1964 END IF
1965 ! Stress-tensor Pz_admm*v_xc[rho_admm]
1966 IF (use_virial) THEN
1967 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1968 END IF
1969 !
1970 IF (admm_env%do_gapw) THEN
1971 CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
1972 IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
1973 CALL update_ks_atom(qs_env, mhx(:, 1), mpz(:, 1), forces=.true., tddft=.false., &
1974 rho_atom_external=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
1975 kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
1976 oce_external=admm_env%admm_gapw_env%oce, &
1977 sab_external=sab_aux_fit)
1978 IF (debug_forces) THEN
1979 fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
1980 CALL para_env%sum(fodeb)
1981 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*Vxc(rho_admm)PAW", fodeb
1982 END IF
1983 END IF
1984 !
1985 NULLIFY (rho_g_aux, rho_r_aux, tau_r_aux)
1986 CALL qs_rho_get(rho_aux_fit, rho_r=rho_r_aux, rho_g=rho_g_aux, tau_r=tau_r_aux)
1987 ! rhoz_aux
1988 NULLIFY (rhoz_g_aux, rhoz_r_aux)
1989 ALLOCATE (rhoz_r_aux(nspins), rhoz_g_aux(nspins))
1990 DO ispin = 1, nspins
1991 CALL auxbas_pw_pool%create_pw(rhoz_r_aux(ispin))
1992 CALL auxbas_pw_pool%create_pw(rhoz_g_aux(ispin))
1993 END DO
1994 DO ispin = 1, nspins
1995 CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpz(ispin, 1)%matrix, &
1996 rho=rhoz_r_aux(ispin), rho_gspace=rhoz_g_aux(ispin), &
1997 basis_type=basis_type, task_list_external=task_list)
1998 END DO
1999 !
2000 ! Add ADMM volume contribution to stress tensor
2001 IF (use_virial) THEN
2002
2003 ! Stress tensor volume term: \int v_xc[n_in_admm]*n_z_admm
2004 ! vadmm_rspace already scaled, we need to unscale it!
2005 DO ispin = 1, nspins
2006 exc_aux_fit = exc_aux_fit + pw_integral_ab(rhoz_r_aux(ispin), vadmm_rspace(ispin))/ &
2007 vadmm_rspace(ispin)%pw_grid%dvol
2008 END DO
2009
2010 IF (debug_stress) THEN
2011 stdeb = -1.0_dp*fconv*exc_aux_fit
2012 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T43,2(1X,ES19.11))") &
2013 'STRESS| VOL 1st eps_XC[n_in_admm]*n_z_admm', one_third_sum_diag(stdeb), det_3x3(stdeb)
2014 END IF
2015
2016 END IF
2017 !
2018 NULLIFY (v_xc)
2019
2020 IF (use_virial) virial%pv_xc = 0.0_dp
2021
2022 CALL create_kernel(qs_env=qs_env, &
2023 vxc=v_xc, &
2024 vxc_tau=v_xc_tau, &
2025 rho=rho_aux_fit, &
2026 rho1_r=rhoz_r_aux, &
2027 rho1_g=rhoz_g_aux, &
2028 tau1_r=tau_r_aux, &
2029 xc_section=xc_section, &
2030 compute_virial=use_virial, &
2031 virial_xc=virial%pv_xc)
2032
2033 ! Stress-tensor ADMM-kernel GGA contribution
2034 IF (use_virial) THEN
2035 virial%pv_exc = virial%pv_exc + virial%pv_xc
2036 virial%pv_virial = virial%pv_virial + virial%pv_xc
2037 END IF
2038
2039 IF (debug_stress .AND. use_virial) THEN
2040 stdeb = 1.0_dp*fconv*virial%pv_xc
2041 CALL para_env%sum(stdeb)
2042 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2043 'STRESS| GGA 2nd Pin_admm*dK*rhoz_admm', one_third_sum_diag(stdeb), det_3x3(stdeb)
2044 END IF
2045 !
2046 CALL qs_rho_get(rho_aux_fit, rho_ao_kp=matrix_p)
2047 ! Stress-tensor Pin*dK*rhoz_admm
2048 IF (use_virial) THEN
2049 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2050 END IF
2051 IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2052 IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2053 DO ispin = 1, nspins
2054 CALL dbcsr_set(mhy(ispin, 1)%matrix, 0.0_dp)
2055 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2056 CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
2057 hmat=mhy(ispin, 1), pmat=matrix_p(ispin, 1), &
2058 calculate_forces=.true., &
2059 basis_type=basis_type, task_list_external=task_list)
2060 END DO
2061 IF (debug_forces) THEN
2062 fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2063 CALL para_env%sum(fodeb)
2064 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dK*rhoz_admm ", fodeb
2065 END IF
2066 IF (debug_stress .AND. use_virial) THEN
2067 stdeb = fconv*(virial%pv_virial - pv_loc)
2068 CALL para_env%sum(stdeb)
2069 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2070 'STRESS| INT 2nd Pin*dK*rhoz_admm ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2071 END IF
2072 ! Stress-tensor Pin*dK*rhoz_admm
2073 IF (use_virial) THEN
2074 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2075 END IF
2076 DO ispin = 1, nspins
2077 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2078 CALL auxbas_pw_pool%give_back_pw(rhoz_r_aux(ispin))
2079 CALL auxbas_pw_pool%give_back_pw(rhoz_g_aux(ispin))
2080 END DO
2081 DEALLOCATE (v_xc, rhoz_r_aux, rhoz_g_aux)
2082 !
2083 IF (admm_env%do_gapw) THEN
2084 CALL local_rho_set_create(local_rhoz_set_admm)
2085 CALL allocate_rho_atom_internals(local_rhoz_set_admm%rho_atom_set, atomic_kind_set, &
2086 admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
2087 CALL calculate_rho_atom_coeff(qs_env, mpz(:, 1), local_rhoz_set_admm%rho_atom_set, &
2088 admm_env%admm_gapw_env%admm_kind_set, &
2089 admm_env%admm_gapw_env%oce, sab_aux_fit, para_env)
2090 CALL prepare_gapw_den(qs_env, local_rho_set=local_rhoz_set_admm, &
2091 do_rho0=.false., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
2092 !compute the potential due to atomic densities
2093 CALL calculate_xc_2nd_deriv_atom(admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
2094 local_rhoz_set_admm%rho_atom_set, &
2095 qs_env, xc_section, para_env, do_triplet=.false., &
2096 kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
2097 !
2098 IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
2099 CALL update_ks_atom(qs_env, mhy(:, 1), matrix_p(:, 1), forces=.true., tddft=.false., &
2100 rho_atom_external=local_rhoz_set_admm%rho_atom_set, &
2101 kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
2102 oce_external=admm_env%admm_gapw_env%oce, &
2103 sab_external=sab_aux_fit)
2104 IF (debug_forces) THEN
2105 fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
2106 CALL para_env%sum(fodeb)
2107 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dK*rhoz_admm[PAW] ", fodeb
2108 END IF
2109 CALL local_rho_set_release(local_rhoz_set_admm)
2110 END IF
2111 !
2112 nao = admm_env%nao_orb
2113 nao_aux = admm_env%nao_aux_fit
2114 ALLOCATE (dbwork)
2115 CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
2116 DO ispin = 1, nspins
2117 CALL cp_dbcsr_sm_fm_multiply(mhy(ispin, 1)%matrix, admm_env%A, &
2118 admm_env%work_aux_orb, nao)
2119 CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
2120 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
2121 admm_env%work_orb_orb)
2122 CALL dbcsr_copy(dbwork, matrix_hz(ispin)%matrix)
2123 CALL dbcsr_set(dbwork, 0.0_dp)
2124 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.true.)
2125 CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
2126 END DO
2127 CALL dbcsr_release(dbwork)
2128 DEALLOCATE (dbwork)
2129 CALL dbcsr_deallocate_matrix_set(mpz)
2130 END IF ! qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none
2131 END IF ! do_admm
2132
2133 ! -----------------------------------------
2134 ! HFX
2135 ! -----------------------------------------
2136
2137 ! HFX
2138 hfx_section => section_vals_get_subs_vals(xc_section, "HF")
2139 CALL section_vals_get(hfx_section, explicit=do_hfx)
2140 IF (do_hfx) THEN
2141 CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
2142 cpassert(n_rep_hf == 1)
2143 CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
2144 i_rep_section=1)
2145 mspin = 1
2146 IF (hfx_treat_lsd_in_core) mspin = nspins
2147 IF (use_virial) virial%pv_fock_4c = 0.0_dp
2148 !
2149 CALL get_qs_env(qs_env=qs_env, rho=rho, x_data=x_data, &
2150 s_mstruct_changed=s_mstruct_changed)
2151 distribute_fock_matrix = .true.
2152
2153 ! -----------------------------------------
2154 ! HFX-ADMM
2155 ! -----------------------------------------
2156 IF (dft_control%do_admm) THEN
2157 CALL get_qs_env(qs_env=qs_env, admm_env=admm_env)
2158 CALL get_admm_env(admm_env, matrix_s_aux_fit=scrm, rho_aux_fit=rho_aux_fit)
2159 CALL qs_rho_get(rho_aux_fit, rho_ao_kp=matrix_p)
2160 NULLIFY (mpz, mhz, mpd, mhd)
2161 CALL dbcsr_allocate_matrix_set(mpz, nspins, 1)
2162 CALL dbcsr_allocate_matrix_set(mhz, nspins, 1)
2163 CALL dbcsr_allocate_matrix_set(mpd, nspins, 1)
2164 CALL dbcsr_allocate_matrix_set(mhd, nspins, 1)
2165 DO ispin = 1, nspins
2166 ALLOCATE (mhz(ispin, 1)%matrix, mhd(ispin, 1)%matrix)
2167 CALL dbcsr_create(mhz(ispin, 1)%matrix, template=scrm(1)%matrix)
2168 CALL dbcsr_create(mhd(ispin, 1)%matrix, template=scrm(1)%matrix)
2169 CALL dbcsr_copy(mhz(ispin, 1)%matrix, scrm(1)%matrix)
2170 CALL dbcsr_copy(mhd(ispin, 1)%matrix, scrm(1)%matrix)
2171 CALL dbcsr_set(mhz(ispin, 1)%matrix, 0.0_dp)
2172 CALL dbcsr_set(mhd(ispin, 1)%matrix, 0.0_dp)
2173 ALLOCATE (mpz(ispin, 1)%matrix)
2174 IF (do_ex) THEN
2175 CALL dbcsr_create(mpz(ispin, 1)%matrix, template=scrm(1)%matrix)
2176 CALL dbcsr_copy(mpz(ispin, 1)%matrix, p_env%p1_admm(ispin)%matrix)
2177 CALL dbcsr_add(mpz(ispin, 1)%matrix, ex_env%matrix_pe_admm(ispin)%matrix, &
2178 1.0_dp, 1.0_dp)
2179 ELSE
2180 CALL dbcsr_create(mpz(ispin, 1)%matrix, template=scrm(1)%matrix)
2181 CALL dbcsr_copy(mpz(ispin, 1)%matrix, matrix_pz_admm(ispin)%matrix)
2182 END IF
2183 mpd(ispin, 1)%matrix => matrix_p(ispin, 1)%matrix
2184 END DO
2185 !
2186 IF (x_data(1, 1)%do_hfx_ri) THEN
2187
2188 eh1 = 0.0_dp
2189 CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpz, &
2190 geometry_did_change=s_mstruct_changed, nspins=nspins, &
2191 hf_fraction=x_data(1, 1)%general_parameter%fraction)
2192
2193 eh1 = 0.0_dp
2194 CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhd, eh1, rho_ao=mpd, &
2195 geometry_did_change=s_mstruct_changed, nspins=nspins, &
2196 hf_fraction=x_data(1, 1)%general_parameter%fraction)
2197
2198 ELSE
2199 DO ispin = 1, mspin
2200 eh1 = 0.0
2201 CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpz, hfx_section, &
2202 para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
2203 ispin=ispin)
2204 END DO
2205 DO ispin = 1, mspin
2206 eh1 = 0.0
2207 CALL integrate_four_center(qs_env, x_data, mhd, eh1, mpd, hfx_section, &
2208 para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
2209 ispin=ispin)
2210 END DO
2211 END IF
2212 !
2213 CALL get_qs_env(qs_env, admm_env=admm_env)
2214 cpassert(ASSOCIATED(admm_env%work_aux_orb))
2215 cpassert(ASSOCIATED(admm_env%work_orb_orb))
2216 nao = admm_env%nao_orb
2217 nao_aux = admm_env%nao_aux_fit
2218 ALLOCATE (dbwork)
2219 CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
2220 DO ispin = 1, nspins
2221 CALL cp_dbcsr_sm_fm_multiply(mhz(ispin, 1)%matrix, admm_env%A, &
2222 admm_env%work_aux_orb, nao)
2223 CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
2224 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
2225 admm_env%work_orb_orb)
2226 CALL dbcsr_copy(dbwork, matrix_hz(ispin)%matrix)
2227 CALL dbcsr_set(dbwork, 0.0_dp)
2228 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.true.)
2229 CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
2230 END DO
2231 CALL dbcsr_release(dbwork)
2232 DEALLOCATE (dbwork)
2233 ! derivatives Tr (Pz [A(T)H dA/dR])
2234 IF (debug_forces) fodeb(1:3) = force(1)%overlap_admm(1:3, 1)
2235 IF (ASSOCIATED(mhx) .AND. ASSOCIATED(mhy)) THEN
2236 DO ispin = 1, nspins
2237 CALL dbcsr_add(mhd(ispin, 1)%matrix, mhx(ispin, 1)%matrix, 1.0_dp, 1.0_dp)
2238 CALL dbcsr_add(mhz(ispin, 1)%matrix, mhy(ispin, 1)%matrix, 1.0_dp, 1.0_dp)
2239 END DO
2240 END IF
2241 CALL qs_rho_get(rho, rho_ao=matrix_pd)
2242 CALL admm_projection_derivative(qs_env, mhd(:, 1), mpa)
2243 CALL admm_projection_derivative(qs_env, mhz(:, 1), matrix_pd)
2244 IF (debug_forces) THEN
2245 fodeb(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb(1:3)
2246 CALL para_env%sum(fodeb)
2247 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*hfx*S' ", fodeb
2248 END IF
2249 CALL dbcsr_deallocate_matrix_set(mpz)
2250 CALL dbcsr_deallocate_matrix_set(mhz)
2251 CALL dbcsr_deallocate_matrix_set(mhd)
2252 IF (ASSOCIATED(mhx) .AND. ASSOCIATED(mhy)) THEN
2253 CALL dbcsr_deallocate_matrix_set(mhx)
2254 CALL dbcsr_deallocate_matrix_set(mhy)
2255 END IF
2256 DEALLOCATE (mpd)
2257 ELSE
2258 ! -----------------------------------------
2259 ! conventional HFX
2260 ! -----------------------------------------
2261 ALLOCATE (mpz(nspins, 1), mhz(nspins, 1))
2262 DO ispin = 1, nspins
2263 mhz(ispin, 1)%matrix => matrix_hz(ispin)%matrix
2264 mpz(ispin, 1)%matrix => mpa(ispin)%matrix
2265 END DO
2266
2267 IF (x_data(1, 1)%do_hfx_ri) THEN
2268
2269 eh1 = 0.0_dp
2270 CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpz, &
2271 geometry_did_change=s_mstruct_changed, nspins=nspins, &
2272 hf_fraction=x_data(1, 1)%general_parameter%fraction)
2273 ELSE
2274 DO ispin = 1, mspin
2275 eh1 = 0.0
2276 CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpz, hfx_section, &
2277 para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
2278 ispin=ispin)
2279 END DO
2280 END IF
2281 DEALLOCATE (mhz, mpz)
2282 END IF
2283
2284 ! -----------------------------------------
2285 ! HFX FORCES
2286 ! -----------------------------------------
2287
2288 resp_only = .true.
2289 IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
2290 IF (dft_control%do_admm) THEN
2291 ! -----------------------------------------
2292 ! HFX-ADMM FORCES
2293 ! -----------------------------------------
2294 CALL qs_rho_get(rho_aux_fit, rho_ao_kp=matrix_p)
2295 NULLIFY (matrix_pza)
2296 CALL dbcsr_allocate_matrix_set(matrix_pza, nspins)
2297 DO ispin = 1, nspins
2298 ALLOCATE (matrix_pza(ispin)%matrix)
2299 IF (do_ex) THEN
2300 CALL dbcsr_create(matrix_pza(ispin)%matrix, template=p_env%p1_admm(ispin)%matrix)
2301 CALL dbcsr_copy(matrix_pza(ispin)%matrix, p_env%p1_admm(ispin)%matrix)
2302 CALL dbcsr_add(matrix_pza(ispin)%matrix, ex_env%matrix_pe_admm(ispin)%matrix, &
2303 1.0_dp, 1.0_dp)
2304 ELSE
2305 CALL dbcsr_create(matrix_pza(ispin)%matrix, template=matrix_pz_admm(ispin)%matrix)
2306 CALL dbcsr_copy(matrix_pza(ispin)%matrix, matrix_pz_admm(ispin)%matrix)
2307 END IF
2308 END DO
2309 IF (x_data(1, 1)%do_hfx_ri) THEN
2310
2311 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
2312 x_data(1, 1)%general_parameter%fraction, &
2313 rho_ao=matrix_p, rho_ao_resp=matrix_pza, &
2314 use_virial=use_virial, resp_only=resp_only)
2315 ELSE
2316 CALL derivatives_four_center(qs_env, matrix_p, matrix_pza, hfx_section, para_env, &
2317 1, use_virial, resp_only=resp_only)
2318 END IF
2319 CALL dbcsr_deallocate_matrix_set(matrix_pza)
2320 ELSE
2321 ! -----------------------------------------
2322 ! conventional HFX FORCES
2323 ! -----------------------------------------
2324 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2325 IF (x_data(1, 1)%do_hfx_ri) THEN
2326
2327 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
2328 x_data(1, 1)%general_parameter%fraction, &
2329 rho_ao=matrix_p, rho_ao_resp=mpa, &
2330 use_virial=use_virial, resp_only=resp_only)
2331 ELSE
2332 CALL derivatives_four_center(qs_env, matrix_p, mpa, hfx_section, para_env, &
2333 1, use_virial, resp_only=resp_only)
2334 END IF
2335 END IF ! do_admm
2336
2337 IF (use_virial) THEN
2338 virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
2339 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
2340 virial%pv_calculate = .false.
2341 END IF
2342
2343 IF (debug_forces) THEN
2344 fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
2345 CALL para_env%sum(fodeb)
2346 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*hfx ", fodeb
2347 END IF
2348 IF (debug_stress .AND. use_virial) THEN
2349 stdeb = -1.0_dp*fconv*virial%pv_fock_4c
2350 CALL para_env%sum(stdeb)
2351 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2352 'STRESS| Pz*hfx ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2353 END IF
2354 END IF ! do_hfx
2355
2356 ! Stress-tensor volume contributions
2357 ! These need to be applied at the end of qs_force
2358 IF (use_virial) THEN
2359 ! Adding mixed Hartree energy twice, due to symmetry
2360 zehartree = zehartree + 2.0_dp*ehartree
2361 zexc = zexc + exc
2362 ! ADMM contribution handled differently in qs_force
2363 IF (dft_control%do_admm) THEN
2364 zexc_aux_fit = zexc_aux_fit + exc_aux_fit
2365 END IF
2366 END IF
2367
2368 ! Overlap matrix
2369 ! H(drho+dz) + Wz
2370 ! If ground-state density matrix solved by diagonalization, then use this
2371 IF (dft_control%qs_control%do_ls_scf) THEN
2372 ! Ground-state density has been calculated by LS
2373 eps_filter = dft_control%qs_control%eps_filter_matrix
2374 CALL calculate_whz_ao_matrix(qs_env, matrix_hz, matrix_wz, eps_filter)
2375 ELSE
2376 IF (do_ex) THEN
2377 matrix_wz => p_env%w1
2378 END IF
2379 focc = 1.0_dp
2380 IF (nspins == 1) focc = 2.0_dp
2381 CALL get_qs_env(qs_env, mos=mos)
2382 DO ispin = 1, nspins
2383 CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
2384 CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
2385 matrix_wz(ispin)%matrix, focc, nocc)
2386 END DO
2387 END IF
2388 IF (nspins == 2) THEN
2389 CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
2390 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2391 END IF
2392
2393 IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
2394 IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
2395 NULLIFY (scrm)
2396 CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
2397 matrix_name="OVERLAP MATRIX", &
2398 basis_type_a="ORB", basis_type_b="ORB", &
2399 sab_nl=sab_orb, calculate_forces=.true., &
2400 matrix_p=matrix_wz(1)%matrix)
2401
2402 IF (SIZE(matrix_wz, 1) == 2) THEN
2403 CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
2404 alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2405 END IF
2406
2407 IF (debug_forces) THEN
2408 fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
2409 CALL para_env%sum(fodeb)
2410 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wz*dS ", fodeb
2411 END IF
2412 IF (debug_stress .AND. use_virial) THEN
2413 stdeb = fconv*(virial%pv_overlap - stdeb)
2414 CALL para_env%sum(stdeb)
2415 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2416 'STRESS| WHz ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2417 END IF
2418 CALL dbcsr_deallocate_matrix_set(scrm)
2419
2420 IF (debug_forces) THEN
2421 CALL total_qs_force(ftot2, force, atomic_kind_set)
2422 fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
2423 CALL para_env%sum(fodeb)
2424 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Response Force", fodeb
2425 fodeb(1:3) = ftot2(1:3, 1)
2426 CALL para_env%sum(fodeb)
2427 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Total Force ", fodeb
2428 DEALLOCATE (ftot1, ftot2, ftot3)
2429 END IF
2430 IF (debug_stress .AND. use_virial) THEN
2431 stdeb = fconv*(virial%pv_virial - sttot)
2432 CALL para_env%sum(stdeb)
2433 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2434 'STRESS| Stress Response ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2435 stdeb = fconv*(virial%pv_virial)
2436 CALL para_env%sum(stdeb)
2437 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,A,T41,2(1X,ES19.11))") &
2438 'STRESS| Total Stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2439 IF (iounit > 0) WRITE (unit=iounit, fmt="(T2,3(1X,ES19.11))") &
2440 stdeb(1, 1), stdeb(2, 2), stdeb(3, 3)
2441 unitstr = "bar"
2442 END IF
2443
2444 IF (do_ex) THEN
2445 CALL dbcsr_deallocate_matrix_set(mpa)
2446 CALL dbcsr_deallocate_matrix_set(matrix_hz)
2447 END IF
2448
2449 CALL timestop(handle)
2450
2451 END SUBROUTINE response_force
2452
2453! **************************************************************************************************
2454!> \brief ...
2455!> \param qs_env ...
2456!> \param p_env ...
2457!> \param matrix_hz ...
2458!> \param ex_env ...
2459!> \param debug ...
2460! **************************************************************************************************
2461 SUBROUTINE response_force_xtb(qs_env, p_env, matrix_hz, ex_env, debug)
2462 TYPE(qs_environment_type), POINTER :: qs_env
2463 TYPE(qs_p_env_type) :: p_env
2464 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz
2465 TYPE(excited_energy_type), OPTIONAL, POINTER :: ex_env
2466 LOGICAL, INTENT(IN), OPTIONAL :: debug
2467
2468 CHARACTER(LEN=*), PARAMETER :: routinen = 'response_force_xtb'
2469
2470 INTEGER :: atom_a, handle, iatom, ikind, iounit, &
2471 is, ispin, na, natom, natorb, nimages, &
2472 nkind, nocc, ns, nsgf, nspins
2473 INTEGER, DIMENSION(25) :: lao
2474 INTEGER, DIMENSION(5) :: occ
2475 LOGICAL :: debug_forces, do_ex, use_virial
2476 REAL(kind=dp) :: focc
2477 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: mcharge, mcharge1
2478 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, charges, charges1, ftot1, &
2479 ftot2
2480 REAL(kind=dp), DIMENSION(3) :: fodeb
2481 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2482 TYPE(cp_logger_type), POINTER :: logger
2483 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pz, matrix_wz, mpa, p_matrix, scrm
2484 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
2485 TYPE(dbcsr_type), POINTER :: s_matrix
2486 TYPE(dft_control_type), POINTER :: dft_control
2487 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2488 TYPE(mp_para_env_type), POINTER :: para_env
2489 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2490 POINTER :: sab_orb
2491 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2492 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2493 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2494 TYPE(qs_ks_env_type), POINTER :: ks_env
2495 TYPE(qs_rho_type), POINTER :: rho
2496 TYPE(xtb_atom_type), POINTER :: xtb_kind
2497
2498 CALL timeset(routinen, handle)
2499
2500 IF (PRESENT(debug)) THEN
2501 debug_forces = debug
2502 ELSE
2503 debug_forces = .false.
2504 END IF
2505
2506 logger => cp_get_default_logger()
2507 IF (logger%para_env%is_source()) THEN
2508 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
2509 ELSE
2510 iounit = -1
2511 END IF
2512
2513 do_ex = .false.
2514 IF (PRESENT(ex_env)) do_ex = .true.
2515
2516 NULLIFY (ks_env, sab_orb)
2517 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control, &
2518 sab_orb=sab_orb)
2519 CALL get_qs_env(qs_env=qs_env, para_env=para_env, force=force)
2520 nspins = dft_control%nspins
2521
2522 IF (debug_forces) THEN
2523 CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
2524 ALLOCATE (ftot1(3, natom))
2525 ALLOCATE (ftot2(3, natom))
2526 CALL total_qs_force(ftot1, force, atomic_kind_set)
2527 END IF
2528
2529 matrix_pz => p_env%p1
2530 NULLIFY (mpa)
2531 IF (do_ex) THEN
2532 CALL dbcsr_allocate_matrix_set(mpa, nspins)
2533 DO ispin = 1, nspins
2534 ALLOCATE (mpa(ispin)%matrix)
2535 CALL dbcsr_create(mpa(ispin)%matrix, template=matrix_pz(ispin)%matrix)
2536 CALL dbcsr_copy(mpa(ispin)%matrix, matrix_pz(ispin)%matrix)
2537 CALL dbcsr_add(mpa(ispin)%matrix, ex_env%matrix_pe(ispin)%matrix, 1.0_dp, 1.0_dp)
2538 CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
2539 END DO
2540 ELSE
2541 mpa => p_env%p1
2542 END IF
2543 !
2544 ! START OF Tr(P+Z)Hcore
2545 !
2546 IF (nspins == 2) THEN
2547 CALL dbcsr_add(mpa(1)%matrix, mpa(2)%matrix, 1.0_dp, 1.0_dp)
2548 END IF
2549 ! Hcore matrix
2550 IF (debug_forces) fodeb(1:3) = force(1)%all_potential(1:3, 1)
2551 CALL build_xtb_hab_force(qs_env, mpa(1)%matrix)
2552 IF (debug_forces) THEN
2553 fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3)
2554 CALL para_env%sum(fodeb)
2555 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHcore ", fodeb
2556 END IF
2557 IF (nspins == 2) THEN
2558 CALL dbcsr_add(mpa(1)%matrix, mpa(2)%matrix, 1.0_dp, -1.0_dp)
2559 END IF
2560 !
2561 ! END OF Tr(P+Z)Hcore
2562 !
2563 use_virial = .false.
2564 nimages = 1
2565 !
2566 ! Hartree potential of response density
2567 !
2568 IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
2569 ! Mulliken charges
2570 CALL get_qs_env(qs_env, rho=rho, particle_set=particle_set, matrix_s_kp=matrix_s)
2571 natom = SIZE(particle_set)
2572 CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2573 ALLOCATE (mcharge(natom), charges(natom, 5))
2574 ALLOCATE (mcharge1(natom), charges1(natom, 5))
2575 charges = 0.0_dp
2576 charges1 = 0.0_dp
2577 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
2578 nkind = SIZE(atomic_kind_set)
2579 CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
2580 ALLOCATE (aocg(nsgf, natom))
2581 aocg = 0.0_dp
2582 ALLOCATE (aocg1(nsgf, natom))
2583 aocg1 = 0.0_dp
2584 p_matrix => matrix_p(:, 1)
2585 s_matrix => matrix_s(1, 1)%matrix
2586 CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
2587 CALL ao_charges(mpa, s_matrix, aocg1, para_env)
2588 DO ikind = 1, nkind
2589 CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
2590 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
2591 CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
2592 DO iatom = 1, na
2593 atom_a = atomic_kind_set(ikind)%atom_list(iatom)
2594 charges(atom_a, :) = real(occ(:), kind=dp)
2595 DO is = 1, natorb
2596 ns = lao(is) + 1
2597 charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
2598 charges1(atom_a, ns) = charges1(atom_a, ns) - aocg1(is, atom_a)
2599 END DO
2600 END DO
2601 END DO
2602 DEALLOCATE (aocg, aocg1)
2603 DO iatom = 1, natom
2604 mcharge(iatom) = sum(charges(iatom, :))
2605 mcharge1(iatom) = sum(charges1(iatom, :))
2606 END DO
2607 ! Coulomb Kernel
2608 CALL xtb_coulomb_hessian(qs_env, matrix_hz, charges1, mcharge1, mcharge)
2609 CALL calc_xtb_ehess_force(qs_env, p_matrix, mpa, charges, mcharge, charges1, &
2610 mcharge1, debug_forces)
2611 !
2612 DEALLOCATE (charges, mcharge, charges1, mcharge1)
2613 END IF
2614 ! Overlap matrix
2615 ! H(drho+dz) + Wz
2616 matrix_wz => p_env%w1
2617 focc = 0.5_dp
2618 IF (nspins == 1) focc = 1.0_dp
2619 CALL get_qs_env(qs_env, mos=mos)
2620 DO ispin = 1, nspins
2621 CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
2622 CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
2623 matrix_wz(ispin)%matrix, focc, nocc)
2624 END DO
2625 IF (nspins == 2) THEN
2626 CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
2627 alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2628 END IF
2629 IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
2630 NULLIFY (scrm)
2631 CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
2632 matrix_name="OVERLAP MATRIX", &
2633 basis_type_a="ORB", basis_type_b="ORB", &
2634 sab_nl=sab_orb, calculate_forces=.true., &
2635 matrix_p=matrix_wz(1)%matrix)
2636 IF (debug_forces) THEN
2637 fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
2638 CALL para_env%sum(fodeb)
2639 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wz*dS ", fodeb
2640 END IF
2641 CALL dbcsr_deallocate_matrix_set(scrm)
2642
2643 IF (debug_forces) THEN
2644 CALL total_qs_force(ftot2, force, atomic_kind_set)
2645 fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
2646 CALL para_env%sum(fodeb)
2647 IF (iounit > 0) WRITE (iounit, "(T3,A,T30,3F16.8)") "DEBUG:: Response Force", fodeb
2648 DEALLOCATE (ftot1, ftot2)
2649 END IF
2650
2651 IF (do_ex) THEN
2652 CALL dbcsr_deallocate_matrix_set(mpa)
2653 END IF
2654
2655 CALL timestop(handle)
2656
2657 END SUBROUTINE response_force_xtb
2658
2659! **************************************************************************************************
2660!> \brief Win = focc*(P*(H[P_out - P_in] + H[Z] )*P)
2661!> Langrange multiplier matrix with response and perturbation (Harris) kernel matrices
2662!>
2663!> \param qs_env ...
2664!> \param matrix_hz ...
2665!> \param matrix_whz ...
2666!> \param eps_filter ...
2667!> \param
2668!> \par History
2669!> 2020.2 created [Fabian Belleflamme]
2670!> \author Fabian Belleflamme
2671! **************************************************************************************************
2672 SUBROUTINE calculate_whz_ao_matrix(qs_env, matrix_hz, matrix_whz, eps_filter)
2673
2674 TYPE(qs_environment_type), POINTER :: qs_env
2675 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
2676 POINTER :: matrix_hz
2677 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
2678 POINTER :: matrix_whz
2679 REAL(kind=dp), INTENT(IN) :: eps_filter
2680
2681 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_whz_ao_matrix'
2682
2683 INTEGER :: handle, ispin, nspins
2684 REAL(kind=dp) :: scaling
2685 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
2686 TYPE(dbcsr_type) :: matrix_tmp
2687 TYPE(dft_control_type), POINTER :: dft_control
2688 TYPE(mp_para_env_type), POINTER :: para_env
2689 TYPE(qs_rho_type), POINTER :: rho
2690
2691 CALL timeset(routinen, handle)
2692
2693 cpassert(ASSOCIATED(qs_env))
2694 cpassert(ASSOCIATED(matrix_hz))
2695 cpassert(ASSOCIATED(matrix_whz))
2696
2697 CALL get_qs_env(qs_env=qs_env, &
2698 dft_control=dft_control, &
2699 rho=rho, &
2700 para_env=para_env)
2701 nspins = dft_control%nspins
2702 CALL qs_rho_get(rho, rho_ao=rho_ao)
2703
2704 ! init temp matrix
2705 CALL dbcsr_create(matrix_tmp, template=matrix_hz(1)%matrix, &
2706 matrix_type=dbcsr_type_no_symmetry)
2707
2708 !Spin factors simplify to
2709 scaling = 1.0_dp
2710 IF (nspins == 1) scaling = 0.5_dp
2711
2712 ! Operation in MO-solver :
2713 ! Whz = focc*(CC^T*Hz*CC^T)
2714 ! focc = 2.0_dp Closed-shell
2715 ! focc = 1.0_dp Open-shell
2716
2717 ! Operation in AO-solver :
2718 ! Whz = (scaling*P)*(focc*Hz)*(scaling*P)
2719 ! focc see above
2720 ! scaling = 0.5_dp Closed-shell (P = 2*CC^T), WHz = (0.5*P)*(2*Hz)*(0.5*P)
2721 ! scaling = 1.0_dp Open-shell, WHz = P*Hz*P
2722
2723 ! Spin factors from Hz and P simplify to
2724 scaling = 1.0_dp
2725 IF (nspins == 1) scaling = 0.5_dp
2726
2727 DO ispin = 1, nspins
2728
2729 ! tmp = H*CC^T
2730 CALL dbcsr_multiply("N", "N", scaling, matrix_hz(ispin)%matrix, rho_ao(ispin)%matrix, &
2731 0.0_dp, matrix_tmp, filter_eps=eps_filter)
2732 ! WHz = CC^T*tmp
2733 ! WHz = Wz + (scaling*P)*(focc*Hz)*(scaling*P)
2734 ! WHz = Wz + scaling*(P*Hz*P)
2735 CALL dbcsr_multiply("N", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_tmp, &
2736 1.0_dp, matrix_whz(ispin)%matrix, filter_eps=eps_filter, &
2737 retain_sparsity=.true.)
2738
2739 END DO
2740
2741 CALL dbcsr_release(matrix_tmp)
2742
2743 CALL timestop(handle)
2744
2745 END SUBROUTINE calculate_whz_ao_matrix
2746
2747! **************************************************************************************************
2748
2749END 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:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
mathlib::det_3x3
Definition mathlib.F:69
mulliken::ao_charges
Definition mulliken.F:42
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
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_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:166
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:533
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:450
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::ec_functional_ext
integer, parameter, public ec_functional_ext
Definition input_constants.F:1164
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:136
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:407
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, monovalent, 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:468
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:913
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:956
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:354
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:399
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:560
response_solver::response_force_xtb
subroutine, public response_force_xtb(qs_env, p_env, matrix_hz, ex_env, debug)
...
Definition response_solver.F:2462
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:824
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:662
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:60
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:624
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:114
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:60
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:721
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:47
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