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