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