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