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