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