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