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