Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculate the CPKS equation and the resulting forces
10 : !> \par History
11 : !> 03.2014 created
12 : !> 09.2019 Moved from KG to Kohn-Sham
13 : !> 11.2019 Moved from energy_correction
14 : !> 08.2020 AO linear response solver [fbelle]
15 : !> \author JGH
16 : ! **************************************************************************************************
17 : MODULE response_solver
18 : USE accint_weights_forces, ONLY: accint_weight_force
19 : USE admm_methods, ONLY: admm_projection_derivative
20 : USE admm_types, ONLY: admm_type,&
21 : get_admm_env
22 : USE atomic_kind_types, ONLY: atomic_kind_type,&
23 : get_atomic_kind
24 : USE cell_types, ONLY: cell_type
25 : USE cp_blacs_env, ONLY: cp_blacs_env_type
26 : USE cp_control_types, ONLY: dft_control_type
27 : USE cp_dbcsr_api, ONLY: &
28 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_multiply, &
29 : dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
30 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
31 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
32 : copy_fm_to_dbcsr,&
33 : cp_dbcsr_sm_fm_multiply,&
34 : dbcsr_allocate_matrix_set,&
35 : dbcsr_deallocate_matrix_set
36 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
37 : cp_fm_struct_release,&
38 : cp_fm_struct_type
39 : USE cp_fm_types, ONLY: cp_fm_create,&
40 : cp_fm_init_random,&
41 : cp_fm_release,&
42 : cp_fm_set_all,&
43 : cp_fm_to_fm,&
44 : cp_fm_type
45 : USE cp_log_handling, ONLY: cp_get_default_logger,&
46 : cp_logger_get_default_unit_nr,&
47 : cp_logger_type
48 : USE ec_env_types, ONLY: energy_correction_type
49 : USE ec_methods, ONLY: ec_mos_init
50 : USE ec_orth_solver, ONLY: ec_response_ao
51 : USE exstates_types, ONLY: excited_energy_type
52 : USE hartree_local_methods, ONLY: Vh_1c_gg_integrals,&
53 : init_coulomb_local
54 : USE hartree_local_types, ONLY: hartree_local_create,&
55 : hartree_local_release,&
56 : hartree_local_type
57 : USE hfx_derivatives, ONLY: derivatives_four_center
58 : USE hfx_energy_potential, ONLY: integrate_four_center
59 : USE hfx_ri, ONLY: hfx_ri_update_forces,&
60 : hfx_ri_update_ks
61 : USE hfx_types, ONLY: hfx_type
62 : USE input_constants, ONLY: &
63 : do_admm_aux_exch_func_none, ec_functional_ext, ec_ls_solver, ec_mo_solver, &
64 : kg_tnadd_atomic, kg_tnadd_embed, kg_tnadd_embed_ri, ls_s_sqrt_ns, ls_s_sqrt_proot, &
65 : ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, &
66 : ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, precond_mlp, xc_none
67 : USE input_section_types, ONLY: section_vals_get,&
68 : section_vals_get_subs_vals,&
69 : section_vals_type,&
70 : section_vals_val_get
71 : USE kg_correction, ONLY: kg_ekin_subset
72 : USE kg_environment_types, ONLY: kg_environment_type
73 : USE kg_tnadd_mat, ONLY: build_tnadd_mat
74 : USE kinds, ONLY: default_string_length,&
75 : dp
76 : USE machine, ONLY: m_flush
77 : USE mathlib, ONLY: det_3x3
78 : USE message_passing, ONLY: mp_para_env_type
79 : USE mulliken, ONLY: ao_charges
80 : USE parallel_gemm_api, ONLY: parallel_gemm
81 : USE particle_types, ONLY: particle_type
82 : USE physcon, ONLY: pascal
83 : USE pw_env_types, ONLY: pw_env_get,&
84 : pw_env_type
85 : USE pw_methods, ONLY: pw_axpy,&
86 : pw_copy,&
87 : pw_integral_ab,&
88 : pw_scale,&
89 : pw_transfer,&
90 : pw_zero
91 : USE pw_poisson_methods, ONLY: pw_poisson_solve
92 : USE pw_poisson_types, ONLY: pw_poisson_type
93 : USE pw_pool_types, ONLY: pw_pool_type
94 : USE pw_types, ONLY: pw_c1d_gs_type,&
95 : pw_r3d_rs_type
96 : USE qs_2nd_kernel_ao, ONLY: build_dm_response
97 : USE qs_collocate_density, ONLY: calculate_rho_elec
98 : USE qs_core_matrices, ONLY: core_matrices,&
99 : kinetic_energy_matrix
100 : USE qs_density_matrices, ONLY: calculate_whz_matrix,&
101 : calculate_wz_matrix
102 : USE qs_energy_types, ONLY: qs_energy_type
103 : USE qs_environment_types, ONLY: get_qs_env,&
104 : qs_environment_type,&
105 : set_qs_env
106 : USE qs_force_types, ONLY: qs_force_type,&
107 : total_qs_force
108 : USE qs_fxc, ONLY: qs_fxc_create
109 : USE qs_gapw_densities, ONLY: prepare_gapw_den
110 : USE qs_integrate_potential, ONLY: integrate_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_atom_types, ONLY: rho_atom_type
141 : USE qs_rho_types, ONLY: qs_rho_create,&
142 : qs_rho_get,&
143 : qs_rho_set,&
144 : qs_rho_type
145 : USE qs_vxc_atom, ONLY: calculate_vxc_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 2560 : 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 668 : 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 668 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: psi0, psi1
702 : TYPE(cp_fm_type), POINTER :: mo_coeff
703 668 : 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 668 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
707 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
708 668 : POINTER :: sab_orb
709 :
710 668 : CALL timeset(routineN, handle)
711 :
712 : ! initialized linres_control
713 : NULLIFY (linres_control)
714 668 : ALLOCATE (linres_control)
715 668 : linres_control%do_kernel = .TRUE.
716 : linres_control%lr_triplet = .FALSE.
717 668 : IF (PRESENT(lr_section)) THEN
718 668 : CALL section_vals_val_get(lr_section, "RESTART", l_val=linres_control%linres_restart)
719 668 : CALL section_vals_val_get(lr_section, "MAX_ITER", i_val=linres_control%max_iter)
720 668 : CALL section_vals_val_get(lr_section, "EPS", r_val=linres_control%eps)
721 668 : CALL section_vals_val_get(lr_section, "EPS_FILTER", r_val=linres_control%eps_filter)
722 668 : CALL section_vals_val_get(lr_section, "RESTART_EVERY", i_val=linres_control%restart_every)
723 668 : CALL section_vals_val_get(lr_section, "PRECONDITIONER", i_val=linres_control%preconditioner_type)
724 668 : 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 668 : linres_control=linres_control)
738 668 : CALL set_qs_env(qs_env, linres_control=linres_control)
739 668 : CALL p_env_psi0_changed(p_env, qs_env)
740 668 : p_env%new_preconditioner = .TRUE.
741 :
742 668 : CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos)
743 : !
744 668 : nspins = dft_control%nspins
745 :
746 : ! Initialize vectors:
747 : ! psi0 : The ground-state MO-coefficients
748 : ! psi1 : The "perturbed" linear response orbitals
749 4244 : ALLOCATE (psi0(nspins), psi1(nspins))
750 1454 : DO ispin = 1, nspins
751 786 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, homo=nocc)
752 786 : NULLIFY (fm_struct)
753 : CALL cp_fm_struct_create(fm_struct, ncol_global=nocc, &
754 786 : template_fmstruct=mo_coeff%matrix_struct)
755 786 : CALL cp_fm_create(psi0(ispin), fm_struct)
756 786 : CALL cp_fm_to_fm(mo_coeff, psi0(ispin), nocc)
757 786 : CALL cp_fm_create(psi1(ispin), fm_struct)
758 786 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
759 2240 : CALL cp_fm_struct_release(fm_struct)
760 : END DO
761 :
762 668 : should_stop = .FALSE.
763 : ! The response solver
764 668 : CALL get_qs_env(qs_env, matrix_s=matrix_s, sab_orb=sab_orb)
765 668 : CALL dbcsr_allocate_matrix_set(p_env%p1, nspins)
766 668 : CALL dbcsr_allocate_matrix_set(p_env%w1, nspins)
767 1454 : DO ispin = 1, nspins
768 786 : ALLOCATE (p_env%p1(ispin)%matrix, p_env%w1(ispin)%matrix)
769 786 : CALL dbcsr_create(matrix=p_env%p1(ispin)%matrix, template=matrix_s(1)%matrix)
770 786 : CALL dbcsr_create(matrix=p_env%w1(ispin)%matrix, template=matrix_s(1)%matrix)
771 786 : CALL cp_dbcsr_alloc_block_from_nbl(p_env%p1(ispin)%matrix, sab_orb)
772 1454 : CALL cp_dbcsr_alloc_block_from_nbl(p_env%w1(ispin)%matrix, sab_orb)
773 : END DO
774 668 : 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 668 : should_stop, silent=silent)
788 :
789 : ! Building the response density matrix
790 1454 : DO ispin = 1, nspins
791 1454 : CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_s(1)%matrix)
792 : END DO
793 668 : CALL build_dm_response(psi0, psi1, p_env%p1)
794 1454 : DO ispin = 1, nspins
795 1454 : CALL dbcsr_scale(p_env%p1(ispin)%matrix, 0.5_dp)
796 : END DO
797 668 : 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 668 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
819 1454 : DO ispin = 1, nspins
820 : CALL calculate_wz_matrix(mos(ispin), psi1(ispin), matrix_ks(ispin)%matrix, &
821 1454 : p_env%w1(ispin)%matrix)
822 : END DO
823 668 : CALL cp_fm_release(psi0)
824 668 : CALL cp_fm_release(psi1)
825 :
826 668 : CALL timestop(handle)
827 :
828 2004 : 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, do_onecenter, &
873 : gapw, gapw_xc, hfx_treat_lsd_in_core, needs_tau_response, needs_tau_response_aux, &
874 : resp_only, 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 :: kind_set, 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 1146 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set
923 : TYPE(section_vals_type), POINTER :: hfx_section, xc_fun_section, xc_section
924 : TYPE(task_list_type), POINTER :: task_list, task_list_aux_fit
925 : TYPE(virial_type), POINTER :: virial
926 : TYPE(xc_rho_cflags_type) :: needs
927 :
928 1146 : CALL timeset(routineN, handle)
929 :
930 1146 : IF (PRESENT(debug)) THEN
931 1146 : debug_forces = debug
932 1146 : debug_stress = debug
933 : ELSE
934 0 : debug_forces = .FALSE.
935 0 : debug_stress = .FALSE.
936 : END IF
937 :
938 1146 : logger => cp_get_default_logger()
939 1146 : IF (logger%para_env%is_source()) THEN
940 573 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
941 : ELSE
942 : iounit = -1
943 : END IF
944 :
945 1146 : do_ex = .FALSE.
946 1146 : IF (PRESENT(ex_env)) do_ex = .TRUE.
947 : IF (do_ex) THEN
948 642 : CPASSERT(PRESENT(p_env))
949 : END IF
950 :
951 1146 : NULLIFY (ks_env, sab_orb, virial)
952 : CALL get_qs_env(qs_env=qs_env, &
953 : cell=cell, &
954 : force=force, &
955 : ks_env=ks_env, &
956 : dft_control=dft_control, &
957 : para_env=para_env, &
958 : sab_orb=sab_orb, &
959 1146 : virial=virial)
960 1146 : nspins = dft_control%nspins
961 1146 : gapw = dft_control%qs_control%gapw
962 1146 : gapw_xc = dft_control%qs_control%gapw_xc
963 :
964 1146 : IF (debug_forces) THEN
965 166 : CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
966 498 : ALLOCATE (ftot1(3, natom))
967 166 : CALL total_qs_force(ftot1, force, atomic_kind_set)
968 : END IF
969 :
970 : ! check for virial
971 1146 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
972 :
973 1146 : IF (use_virial .AND. do_ex) THEN
974 0 : CALL cp_abort(__LOCATION__, "Stress Tensor not available for TDDFT calculations.")
975 : END IF
976 :
977 1146 : fconv = 1.0E-9_dp*pascal/cell%deth
978 1146 : IF (debug_stress .AND. use_virial) THEN
979 0 : sttot = virial%pv_virial
980 : END IF
981 :
982 : ! *** If LSD, then combine alpha density and beta density to
983 : ! *** total density: alpha <- alpha + beta and
984 1146 : NULLIFY (mpa)
985 1146 : NULLIFY (matrix_ht)
986 1146 : IF (do_ex) THEN
987 642 : CALL dbcsr_allocate_matrix_set(mpa, nspins)
988 1392 : DO ispin = 1, nspins
989 750 : ALLOCATE (mpa(ispin)%matrix)
990 750 : CALL dbcsr_create(mpa(ispin)%matrix, template=p_env%p1(ispin)%matrix)
991 750 : CALL dbcsr_copy(mpa(ispin)%matrix, p_env%p1(ispin)%matrix)
992 750 : CALL dbcsr_add(mpa(ispin)%matrix, ex_env%matrix_pe(ispin)%matrix, 1.0_dp, 1.0_dp)
993 1392 : CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
994 : END DO
995 : ELSE
996 504 : mpa => matrix_pz
997 : END IF
998 : !
999 1146 : IF (do_ex .OR. (gapw .OR. gapw_xc)) THEN
1000 702 : CALL dbcsr_allocate_matrix_set(matrix_ht, nspins)
1001 1514 : DO ispin = 1, nspins
1002 812 : ALLOCATE (matrix_ht(ispin)%matrix)
1003 812 : CALL dbcsr_create(matrix_ht(ispin)%matrix, template=matrix_hz(ispin)%matrix)
1004 812 : CALL dbcsr_copy(matrix_ht(ispin)%matrix, matrix_hz(ispin)%matrix)
1005 1958 : CALL dbcsr_set(matrix_ht(ispin)%matrix, 0.0_dp)
1006 : END DO
1007 : END IF
1008 : !
1009 : ! START OF Tr[(P+Z)Hcore]
1010 : !
1011 :
1012 : ! Kinetic energy matrix
1013 1146 : NULLIFY (scrm2)
1014 1146 : mpa2(1:nspins, 1:1) => mpa(1:nspins)
1015 : CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm2, matrix_p=mpa2, &
1016 : matrix_name="KINETIC ENERGY MATRIX", &
1017 : basis_type="ORB", &
1018 : sab_orb=sab_orb, calculate_forces=.TRUE., &
1019 1146 : debug_forces=debug_forces, debug_stress=debug_stress)
1020 1146 : CALL dbcsr_deallocate_matrix_set(scrm2)
1021 :
1022 : ! Initialize a matrix scrm, later used for scratch purposes
1023 1146 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
1024 1146 : NULLIFY (scrm)
1025 1146 : CALL dbcsr_allocate_matrix_set(scrm, nspins)
1026 2402 : DO ispin = 1, nspins
1027 1256 : ALLOCATE (scrm(ispin)%matrix)
1028 1256 : CALL dbcsr_create(scrm(ispin)%matrix, template=matrix_s(1)%matrix)
1029 1256 : CALL dbcsr_copy(scrm(ispin)%matrix, matrix_s(1)%matrix)
1030 2402 : CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
1031 : END DO
1032 :
1033 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
1034 1146 : atomic_kind_set=atomic_kind_set)
1035 :
1036 9388 : ALLOCATE (matrix_p(nspins, 1), matrix_h(nspins, 1))
1037 2402 : DO ispin = 1, nspins
1038 1256 : matrix_p(ispin, 1)%matrix => mpa(ispin)%matrix
1039 2402 : matrix_h(ispin, 1)%matrix => scrm(ispin)%matrix
1040 : END DO
1041 1146 : matrix_h(1, 1)%matrix => scrm(1)%matrix
1042 :
1043 1146 : nder = 1
1044 : CALL core_matrices(qs_env, matrix_h, matrix_p, .TRUE., nder, &
1045 1146 : debug_forces=debug_forces, debug_stress=debug_stress)
1046 :
1047 : ! Kim-Gordon subsystem DFT
1048 : ! Atomic potential for nonadditive kinetic energy contribution
1049 1146 : IF (dft_control%qs_control%do_kg) THEN
1050 24 : IF (qs_env%kg_env%tnadd_method == kg_tnadd_atomic) THEN
1051 12 : CALL get_qs_env(qs_env=qs_env, kg_env=kg_env, dbcsr_dist=dbcsr_dist)
1052 :
1053 12 : IF (use_virial) THEN
1054 130 : pv_loc = virial%pv_virial
1055 : END IF
1056 :
1057 12 : IF (debug_forces) fodeb(1:3) = force(1)%kinetic(1:3, 1)
1058 12 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1059 : CALL build_tnadd_mat(kg_env=kg_env, matrix_p=matrix_p, force=force, virial=virial, &
1060 : calculate_forces=.TRUE., use_virial=use_virial, &
1061 : qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
1062 12 : particle_set=particle_set, sab_orb=sab_orb, dbcsr_dist=dbcsr_dist)
1063 12 : IF (debug_forces) THEN
1064 0 : fodeb(1:3) = force(1)%kinetic(1:3, 1) - fodeb(1:3)
1065 0 : CALL para_env%sum(fodeb)
1066 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dTnadd ", fodeb
1067 : END IF
1068 12 : IF (debug_stress .AND. use_virial) THEN
1069 0 : stdeb = fconv*(virial%pv_virial - stdeb)
1070 0 : CALL para_env%sum(stdeb)
1071 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1072 0 : 'STRESS| Pz*dTnadd ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1073 : END IF
1074 :
1075 : ! Stress-tensor update components
1076 12 : IF (use_virial) THEN
1077 130 : virial%pv_ekinetic = virial%pv_ekinetic + (virial%pv_virial - pv_loc)
1078 : END IF
1079 :
1080 : END IF
1081 : END IF
1082 :
1083 1146 : DEALLOCATE (matrix_h)
1084 1146 : DEALLOCATE (matrix_p)
1085 1146 : CALL dbcsr_deallocate_matrix_set(scrm)
1086 :
1087 : ! initialize src matrix
1088 : ! Necessary as build_kinetic_matrix will only allocate scrm(1)
1089 : ! and not scrm(2) in open-shell case
1090 1146 : NULLIFY (scrm)
1091 1146 : CALL dbcsr_allocate_matrix_set(scrm, nspins)
1092 2402 : DO ispin = 1, nspins
1093 1256 : ALLOCATE (scrm(ispin)%matrix)
1094 1256 : CALL dbcsr_create(scrm(ispin)%matrix, template=matrix_pz(1)%matrix)
1095 1256 : CALL dbcsr_copy(scrm(ispin)%matrix, matrix_pz(ispin)%matrix)
1096 2402 : CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
1097 : END DO
1098 :
1099 1146 : IF (debug_forces) THEN
1100 498 : ALLOCATE (ftot2(3, natom))
1101 166 : CALL total_qs_force(ftot2, force, atomic_kind_set)
1102 664 : fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
1103 166 : CALL para_env%sum(fodeb)
1104 166 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (T+Dz)*dHcore", fodeb
1105 : END IF
1106 1146 : IF (debug_stress .AND. use_virial) THEN
1107 0 : stdeb = fconv*(virial%pv_virial - sttot)
1108 0 : CALL para_env%sum(stdeb)
1109 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1110 0 : 'STRESS| Stress Pz*dHcore ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1111 : ! save current total viral, does not contain volume terms yet
1112 0 : sttot2 = virial%pv_virial
1113 : END IF
1114 : !
1115 : ! END OF Tr(P+Z)Hcore
1116 : !
1117 : !
1118 : ! Vhxc (KS potentials calculated externally)
1119 1146 : CALL get_qs_env(qs_env, pw_env=pw_env)
1120 1146 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
1121 : !
1122 1146 : IF (dft_control%do_admm) THEN
1123 256 : CALL get_qs_env(qs_env, admm_env=admm_env)
1124 256 : xc_section => admm_env%xc_section_primary
1125 : ELSE
1126 890 : xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
1127 : END IF
1128 1146 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1129 1146 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
1130 1146 : needs = xc_functionals_get_needs(xc_fun_section, (nspins == 2), .TRUE.)
1131 1146 : needs_tau_response = needs%tau .OR. needs%tau_spin
1132 : !
1133 1146 : IF (gapw .OR. gapw_xc) THEN
1134 216 : NULLIFY (oce, sab_orb)
1135 216 : CALL get_qs_env(qs_env=qs_env, oce=oce, sab_orb=sab_orb)
1136 : ! set up local_rho_set for GS density
1137 216 : NULLIFY (local_rho_set_gs)
1138 216 : CALL get_qs_env(qs_env=qs_env, rho=rho)
1139 216 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1140 216 : CALL local_rho_set_create(local_rho_set_gs)
1141 : CALL allocate_rho_atom_internals(local_rho_set_gs%rho_atom_set, atomic_kind_set, &
1142 216 : qs_kind_set, dft_control, para_env)
1143 216 : CALL init_rho0(local_rho_set_gs, qs_env, dft_control%qs_control%gapw_control)
1144 216 : CALL rho0_s_grid_create(pw_env, local_rho_set_gs%rho0_mpole)
1145 : CALL calculate_rho_atom_coeff(qs_env, matrix_p(:, 1), local_rho_set_gs%rho_atom_set, &
1146 216 : qs_kind_set, oce, sab_orb, para_env)
1147 216 : CALL prepare_gapw_den(qs_env, local_rho_set_gs, do_rho0=gapw)
1148 : ! set up local_rho_set for response density
1149 216 : NULLIFY (local_rho_set_t)
1150 216 : CALL local_rho_set_create(local_rho_set_t)
1151 : CALL allocate_rho_atom_internals(local_rho_set_t%rho_atom_set, atomic_kind_set, &
1152 216 : qs_kind_set, dft_control, para_env)
1153 : CALL init_rho0(local_rho_set_t, qs_env, dft_control%qs_control%gapw_control, &
1154 216 : zcore=0.0_dp)
1155 216 : CALL rho0_s_grid_create(pw_env, local_rho_set_t%rho0_mpole)
1156 : CALL calculate_rho_atom_coeff(qs_env, mpa(:), local_rho_set_t%rho_atom_set, &
1157 216 : qs_kind_set, oce, sab_orb, para_env)
1158 216 : CALL prepare_gapw_den(qs_env, local_rho_set_t, do_rho0=gapw)
1159 :
1160 : ! compute soft GS potential
1161 1516 : ALLOCATE (rho_r_gs(nspins), rho_g_gs(nspins))
1162 434 : DO ispin = 1, nspins
1163 218 : CALL auxbas_pw_pool%create_pw(rho_r_gs(ispin))
1164 434 : CALL auxbas_pw_pool%create_pw(rho_g_gs(ispin))
1165 : END DO
1166 216 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace_gs)
1167 : ! compute soft GS density
1168 216 : total_rho_gs = 0.0_dp
1169 216 : CALL pw_zero(rho_tot_gspace_gs)
1170 434 : DO ispin = 1, nspins
1171 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_p(ispin, 1)%matrix, &
1172 : rho=rho_r_gs(ispin), &
1173 : rho_gspace=rho_g_gs(ispin), &
1174 : soft_valid=(gapw .OR. gapw_xc), &
1175 218 : total_rho=total_rho_gs(ispin))
1176 434 : CALL pw_axpy(rho_g_gs(ispin), rho_tot_gspace_gs)
1177 : END DO
1178 216 : IF (gapw) THEN
1179 176 : CALL get_qs_env(qs_env, natom=natom)
1180 : ! add rho0 contributions to GS density (only for Coulomb) only for gapw
1181 176 : CALL pw_axpy(local_rho_set_gs%rho0_mpole%rho0_s_gs, rho_tot_gspace_gs)
1182 176 : IF (ASSOCIATED(local_rho_set_gs%rho0_mpole%rhoz_cneo_s_gs)) THEN
1183 0 : CALL pw_axpy(local_rho_set_gs%rho0_mpole%rhoz_cneo_s_gs, rho_tot_gspace_gs)
1184 : END IF
1185 176 : IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
1186 8 : CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
1187 8 : CALL pw_axpy(rho_core, rho_tot_gspace_gs)
1188 : END IF
1189 : ! compute GS potential
1190 176 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace_gs)
1191 176 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace_gs)
1192 176 : NULLIFY (hartree_local_gs)
1193 176 : CALL hartree_local_create(hartree_local_gs)
1194 176 : CALL init_coulomb_local(hartree_local_gs, natom)
1195 176 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace_gs, hartree_gs, v_hartree_gspace_gs)
1196 176 : CALL pw_transfer(v_hartree_gspace_gs, v_hartree_rspace_gs)
1197 176 : CALL pw_scale(v_hartree_rspace_gs, v_hartree_rspace_gs%pw_grid%dvol)
1198 : END IF
1199 : END IF
1200 :
1201 1146 : IF (gapw) THEN
1202 : ! Hartree grid PAW term
1203 176 : CPASSERT(.NOT. use_virial)
1204 584 : IF (debug_forces) fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1)
1205 : CALL Vh_1c_gg_integrals(qs_env, hartree_gs, hartree_local_gs%ecoul_1c, local_rho_set_t, para_env, tddft=.TRUE., &
1206 176 : local_rho_set_2nd=local_rho_set_gs, core_2nd=.FALSE.) ! n^core for GS potential
1207 : ! 1st to define integral space, 2nd for potential, integral contributions stored on local_rho_set_gs
1208 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace_gs, para_env, calculate_forces=.TRUE., &
1209 176 : local_rho_set=local_rho_set_t, local_rho_set_2nd=local_rho_set_gs)
1210 176 : IF (debug_forces) THEN
1211 544 : fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1) - fodeb(1:3)
1212 136 : CALL para_env%sum(fodeb)
1213 136 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (T+Dz)*dVh[D^GS]PAWg0", fodeb
1214 : END IF
1215 : END IF
1216 1146 : IF (gapw .OR. gapw_xc) THEN
1217 216 : IF (myfun /= xc_none) THEN
1218 : ! add 1c hard and soft XC contributions
1219 192 : NULLIFY (local_rho_set_vxc)
1220 192 : CALL local_rho_set_create(local_rho_set_vxc)
1221 : CALL allocate_rho_atom_internals(local_rho_set_vxc%rho_atom_set, atomic_kind_set, &
1222 192 : qs_kind_set, dft_control, para_env)
1223 : CALL calculate_rho_atom_coeff(qs_env, matrix_p(:, 1), local_rho_set_vxc%rho_atom_set, &
1224 192 : qs_kind_set, oce, sab_orb, para_env)
1225 192 : CALL prepare_gapw_den(qs_env, local_rho_set_vxc, do_rho0=.FALSE.)
1226 : ! compute hard and soft atomic contributions
1227 : CALL calculate_vxc_atom(qs_env, .FALSE., exc1=hartree_gs, xc_section_external=xc_section, &
1228 192 : rho_atom_set_external=local_rho_set_vxc%rho_atom_set)
1229 : END IF ! myfun
1230 : END IF ! gapw
1231 :
1232 1146 : CALL auxbas_pw_pool%create_pw(vhxc_rspace)
1233 : !
1234 : ! Stress-tensor: integration contribution direct term
1235 : ! int v_Hxc[n^in]*n^z
1236 1146 : IF (use_virial) THEN
1237 2236 : pv_loc = virial%pv_virial
1238 : END IF
1239 :
1240 1644 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1241 1146 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1242 1146 : IF (gapw .OR. gapw_xc) THEN
1243 : ! vtot = v_xc + v_hartree
1244 434 : DO ispin = 1, nspins
1245 218 : CALL pw_zero(vhxc_rspace)
1246 218 : IF (gapw) THEN
1247 178 : CALL pw_transfer(v_hartree_rspace_gs, vhxc_rspace)
1248 40 : ELSE IF (gapw_xc) THEN
1249 40 : CALL pw_transfer(vh_rspace, vhxc_rspace)
1250 : END IF
1251 : CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1252 : hmat=scrm(ispin), pmat=mpa(ispin), &
1253 : qs_env=qs_env, gapw=gapw, &
1254 434 : calculate_forces=.TRUE.)
1255 : END DO
1256 216 : IF (myfun /= xc_none) THEN
1257 386 : DO ispin = 1, nspins
1258 194 : CALL pw_zero(vhxc_rspace)
1259 194 : CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
1260 : CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1261 : hmat=scrm(ispin), pmat=mpa(ispin), &
1262 : qs_env=qs_env, gapw=(gapw .OR. gapw_xc), &
1263 386 : calculate_forces=.TRUE.)
1264 : END DO
1265 : END IF
1266 : ELSE ! original GPW with Standard Hartree as Potential
1267 1968 : DO ispin = 1, nspins
1268 1038 : CALL pw_transfer(vh_rspace, vhxc_rspace)
1269 1038 : CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
1270 : CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1271 : hmat=scrm(ispin), pmat=mpa(ispin), &
1272 1968 : qs_env=qs_env, gapw=gapw, calculate_forces=.TRUE.)
1273 : END DO
1274 : END IF
1275 :
1276 1146 : IF (debug_forces) THEN
1277 664 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1278 166 : CALL para_env%sum(fodeb)
1279 166 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (T+Dz)*dVhxc[D^GS] ", fodeb
1280 : END IF
1281 1146 : IF (debug_stress .AND. use_virial) THEN
1282 0 : stdeb = fconv*(virial%pv_virial - pv_loc)
1283 0 : CALL para_env%sum(stdeb)
1284 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1285 0 : 'STRESS| INT Pz*dVhxc ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1286 : END IF
1287 :
1288 1146 : IF (gapw .OR. gapw_xc) THEN
1289 : ! HXC term
1290 702 : IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
1291 216 : IF (gapw) CALL update_ks_atom(qs_env, scrm, mpa, forces=.TRUE., tddft=.FALSE., &
1292 176 : rho_atom_external=local_rho_set_gs%rho_atom_set)
1293 216 : IF (myfun /= xc_none) CALL update_ks_atom(qs_env, scrm, mpa, forces=.TRUE., tddft=.FALSE., &
1294 192 : rho_atom_external=local_rho_set_vxc%rho_atom_set)
1295 216 : IF (debug_forces) THEN
1296 648 : fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
1297 162 : CALL para_env%sum(fodeb)
1298 162 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (T+Dz)*dVhxc[D^GS]PAW ", fodeb
1299 : END IF
1300 : ! release local environments for GAPW
1301 216 : IF (myfun /= xc_none) THEN
1302 192 : IF (ASSOCIATED(local_rho_set_vxc)) CALL local_rho_set_release(local_rho_set_vxc)
1303 : END IF
1304 216 : IF (ASSOCIATED(local_rho_set_gs)) CALL local_rho_set_release(local_rho_set_gs)
1305 216 : IF (gapw) THEN
1306 176 : IF (ASSOCIATED(hartree_local_gs)) CALL hartree_local_release(hartree_local_gs)
1307 176 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace_gs)
1308 176 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace_gs)
1309 : END IF
1310 216 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace_gs)
1311 216 : IF (ASSOCIATED(rho_r_gs)) THEN
1312 434 : DO ispin = 1, nspins
1313 434 : CALL auxbas_pw_pool%give_back_pw(rho_r_gs(ispin))
1314 : END DO
1315 216 : DEALLOCATE (rho_r_gs)
1316 : END IF
1317 216 : IF (ASSOCIATED(rho_g_gs)) THEN
1318 434 : DO ispin = 1, nspins
1319 434 : CALL auxbas_pw_pool%give_back_pw(rho_g_gs(ispin))
1320 : END DO
1321 216 : DEALLOCATE (rho_g_gs)
1322 : END IF
1323 : END IF !gapw
1324 :
1325 1146 : IF (ASSOCIATED(vtau_rspace)) THEN
1326 32 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1327 32 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1328 64 : DO ispin = 1, nspins
1329 : CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
1330 : hmat=scrm(ispin), pmat=mpa(ispin), &
1331 : qs_env=qs_env, gapw=(gapw .OR. gapw_xc), &
1332 96 : calculate_forces=.TRUE., compute_tau=.TRUE.)
1333 : END DO
1334 32 : IF (debug_forces) THEN
1335 0 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1336 0 : CALL para_env%sum(fodeb)
1337 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dVxc_tau ", fodeb
1338 : END IF
1339 32 : IF (debug_stress .AND. use_virial) THEN
1340 0 : stdeb = fconv*(virial%pv_virial - pv_loc)
1341 0 : CALL para_env%sum(stdeb)
1342 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1343 0 : 'STRESS| INT Pz*dVxc_tau ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1344 : END IF
1345 : END IF
1346 1146 : CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
1347 :
1348 : ! Stress-tensor Pz*v_Hxc[Pin]
1349 1146 : IF (use_virial) THEN
1350 2236 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1351 : END IF
1352 :
1353 : ! KG Embedding
1354 : ! calculate kinetic energy potential and integrate with response density
1355 1146 : IF (dft_control%qs_control%do_kg) THEN
1356 24 : IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed .OR. &
1357 : qs_env%kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
1358 :
1359 : ekin_mol = 0.0_dp
1360 12 : IF (use_virial) THEN
1361 104 : pv_loc = virial%pv_virial
1362 : END IF
1363 :
1364 12 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1365 : CALL kg_ekin_subset(qs_env=qs_env, &
1366 : ks_matrix=scrm, &
1367 : ekin_mol=ekin_mol, &
1368 : calc_force=.TRUE., &
1369 : do_kernel=.FALSE., &
1370 12 : pmat_ext=mpa)
1371 12 : IF (debug_forces) THEN
1372 0 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1373 0 : CALL para_env%sum(fodeb)
1374 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dVkg ", fodeb
1375 : END IF
1376 12 : IF (debug_stress .AND. use_virial) THEN
1377 : !IF (iounit > 0) WRITE(iounit, *) &
1378 : ! "response_force | VOL 1st KG - v_KG[n_in]*n_z: ", ekin_mol
1379 0 : stdeb = 1.0_dp*fconv*ekin_mol
1380 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1381 0 : 'STRESS| VOL KG Pz*dVKG ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1382 :
1383 0 : stdeb = fconv*(virial%pv_virial - pv_loc)
1384 0 : CALL para_env%sum(stdeb)
1385 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1386 0 : 'STRESS| INT KG Pz*dVKG ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1387 :
1388 0 : stdeb = fconv*virial%pv_xc
1389 0 : CALL para_env%sum(stdeb)
1390 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1391 0 : 'STRESS| GGA KG Pz*dVKG ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1392 : END IF
1393 12 : IF (use_virial) THEN
1394 : ! Direct integral contribution
1395 104 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1396 : END IF
1397 :
1398 : END IF ! tnadd_method
1399 : END IF ! do_kg
1400 :
1401 1146 : CALL dbcsr_deallocate_matrix_set(scrm)
1402 :
1403 : !
1404 : ! Hartree potential of response density
1405 : !
1406 8242 : ALLOCATE (rhoz_r(nspins), rhoz_g(nspins))
1407 2402 : DO ispin = 1, nspins
1408 1256 : CALL auxbas_pw_pool%create_pw(rhoz_r(ispin))
1409 2402 : CALL auxbas_pw_pool%create_pw(rhoz_g(ispin))
1410 : END DO
1411 1146 : CALL auxbas_pw_pool%create_pw(rhoz_tot_gspace)
1412 1146 : CALL auxbas_pw_pool%create_pw(zv_hartree_rspace)
1413 1146 : CALL auxbas_pw_pool%create_pw(zv_hartree_gspace)
1414 :
1415 1146 : CALL pw_zero(rhoz_tot_gspace)
1416 2402 : DO ispin = 1, nspins
1417 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1418 : rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin), &
1419 1256 : soft_valid=gapw)
1420 2402 : CALL pw_axpy(rhoz_g(ispin), rhoz_tot_gspace)
1421 : END DO
1422 1146 : NULLIFY (tauz_r, tauz_r_xc)
1423 1146 : IF (gapw_xc) THEN
1424 200 : ALLOCATE (rhoz_r_xc(nspins), rhoz_g_xc(nspins))
1425 80 : DO ispin = 1, nspins
1426 40 : CALL auxbas_pw_pool%create_pw(rhoz_r_xc(ispin))
1427 80 : CALL auxbas_pw_pool%create_pw(rhoz_g_xc(ispin))
1428 : END DO
1429 80 : DO ispin = 1, nspins
1430 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1431 : rho=rhoz_r_xc(ispin), rho_gspace=rhoz_g_xc(ispin), &
1432 80 : soft_valid=gapw_xc)
1433 : END DO
1434 : END IF
1435 :
1436 1146 : IF (needs_tau_response) THEN
1437 : BLOCK
1438 : TYPE(pw_c1d_gs_type) :: work_g
1439 96 : ALLOCATE (tauz_r(nspins))
1440 32 : CALL auxbas_pw_pool%create_pw(work_g)
1441 64 : DO ispin = 1, nspins
1442 32 : CALL auxbas_pw_pool%create_pw(tauz_r(ispin))
1443 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1444 : rho=tauz_r(ispin), rho_gspace=work_g, &
1445 64 : soft_valid=gapw, compute_tau=.TRUE.)
1446 : END DO
1447 64 : CALL auxbas_pw_pool%give_back_pw(work_g)
1448 : END BLOCK
1449 32 : IF (gapw_xc) THEN
1450 : BLOCK
1451 : TYPE(pw_c1d_gs_type) :: work_g
1452 0 : ALLOCATE (tauz_r_xc(nspins))
1453 0 : CALL auxbas_pw_pool%create_pw(work_g)
1454 0 : DO ispin = 1, nspins
1455 0 : CALL auxbas_pw_pool%create_pw(tauz_r_xc(ispin))
1456 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1457 : rho=tauz_r_xc(ispin), rho_gspace=work_g, &
1458 0 : soft_valid=gapw_xc, compute_tau=.TRUE.)
1459 : END DO
1460 0 : CALL auxbas_pw_pool%give_back_pw(work_g)
1461 : END BLOCK
1462 : END IF
1463 : END IF
1464 :
1465 : !
1466 1146 : IF (PRESENT(rhopz_r)) THEN
1467 1010 : DO ispin = 1, nspins
1468 1010 : CALL pw_copy(rhoz_r(ispin), rhopz_r(ispin))
1469 : END DO
1470 : END IF
1471 :
1472 1146 : ALLOCATE (rho1)
1473 1146 : CALL qs_rho_create(rho1)
1474 1146 : IF (gapw_xc) THEN
1475 40 : CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
1476 40 : rho0 => rho_xc
1477 40 : IF (ASSOCIATED(tauz_r_xc)) THEN
1478 : CALL qs_rho_set(rho1, rho_r=rhoz_r_xc, rho_g=rhoz_g_xc, tau_r=tauz_r_xc, &
1479 0 : rho_r_valid=.TRUE., rho_g_valid=.TRUE., tau_r_valid=.TRUE.)
1480 : ELSE
1481 : CALL qs_rho_set(rho1, rho_r=rhoz_r_xc, rho_g=rhoz_g_xc, &
1482 40 : rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
1483 : END IF
1484 : ELSE
1485 1106 : CALL get_qs_env(qs_env=qs_env, rho=rho)
1486 1106 : rho0 => rho
1487 1106 : IF (ASSOCIATED(tauz_r)) THEN
1488 : CALL qs_rho_set(rho1, rho_r=rhoz_r, rho_g=rhoz_g, tau_r=tauz_r, &
1489 32 : rho_r_valid=.TRUE., rho_g_valid=.TRUE., tau_r_valid=.TRUE.)
1490 : ELSE
1491 : CALL qs_rho_set(rho1, rho_r=rhoz_r, rho_g=rhoz_g, &
1492 1074 : rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
1493 : END IF
1494 : END IF
1495 :
1496 1146 : IF (dft_control%qs_control%gapw_control%accurate_xcint) THEN
1497 : ! GAPW Accurate integration
1498 310 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1499 94 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1500 : !
1501 94 : CALL accint_weight_force(qs_env, rho0, rho1, 1, xc_section)
1502 : !
1503 94 : IF (debug_forces) THEN
1504 288 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1505 72 : CALL para_env%sum(fodeb)
1506 72 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*Vxc*dw ", fodeb
1507 : END IF
1508 94 : IF (debug_stress .AND. use_virial) THEN
1509 0 : stdeb = fconv*(virial%pv_virial - stdeb)
1510 0 : CALL para_env%sum(stdeb)
1511 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1512 0 : 'STRESS| INT Pz*dVxc*dw ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1513 : END IF
1514 : END IF
1515 :
1516 : ! Stress-tensor contribution second derivative
1517 : ! Volume : int v_H[n^z]*n_in
1518 : ! Volume : int epsilon_xc*n_z
1519 1146 : IF (use_virial) THEN
1520 :
1521 172 : CALL get_qs_env(qs_env, rho=rho)
1522 172 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
1523 :
1524 : ! Get the total input density in g-space [ions + electrons]
1525 172 : CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
1526 :
1527 172 : h_stress(:, :) = 0.0_dp
1528 : ! calculate associated hartree potential
1529 : ! This term appears twice in the derivation of the equations
1530 : ! v_H[n_in]*n_z and v_H[n_z]*n_in
1531 : ! due to symmetry we only need to call this routine once,
1532 : ! and count the Volume and Green function contribution
1533 : ! which is stored in h_stress twice
1534 : CALL pw_poisson_solve(poisson_env, &
1535 : density=rhoz_tot_gspace, & ! n_z
1536 : ehartree=ehartree, &
1537 : vhartree=zv_hartree_gspace, & ! v_H[n_z]
1538 : h_stress=h_stress, &
1539 172 : aux_density=rho_tot_gspace) ! n_in
1540 :
1541 172 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1542 :
1543 : ! Stress tensor Green function contribution
1544 2236 : virial%pv_ehartree = virial%pv_ehartree + 2.0_dp*h_stress/REAL(para_env%num_pe, dp)
1545 2236 : virial%pv_virial = virial%pv_virial + 2.0_dp*h_stress/REAL(para_env%num_pe, dp)
1546 :
1547 172 : IF (debug_stress) THEN
1548 0 : stdeb = -1.0_dp*fconv*ehartree
1549 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1550 0 : 'STRESS| VOL 1st v_H[n_z]*n_in ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1551 0 : stdeb = -1.0_dp*fconv*ehartree
1552 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1553 0 : 'STRESS| VOL 2nd v_H[n_in]*n_z ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1554 0 : stdeb = fconv*(h_stress/REAL(para_env%num_pe, dp))
1555 0 : CALL para_env%sum(stdeb)
1556 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1557 0 : 'STRESS| GREEN 1st v_H[n_z]*n_in ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1558 0 : stdeb = fconv*(h_stress/REAL(para_env%num_pe, dp))
1559 0 : CALL para_env%sum(stdeb)
1560 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1561 0 : 'STRESS| GREEN 2nd v_H[n_in]*n_z ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1562 : END IF
1563 :
1564 : ! Stress tensor volume term: \int v_xc[n_in]*n_z
1565 : ! vxc_rspace already scaled, we need to unscale it!
1566 172 : exc = 0.0_dp
1567 344 : DO ispin = 1, nspins
1568 : exc = exc + pw_integral_ab(rhoz_r(ispin), vxc_rspace(ispin))/ &
1569 344 : vxc_rspace(ispin)%pw_grid%dvol
1570 : END DO
1571 172 : IF (ASSOCIATED(vtau_rspace)) THEN
1572 32 : DO ispin = 1, nspins
1573 : exc = exc + pw_integral_ab(tauz_r(ispin), vtau_rspace(ispin))/ &
1574 32 : vtau_rspace(ispin)%pw_grid%dvol
1575 : END DO
1576 : END IF
1577 :
1578 : ! Add KG embedding correction
1579 172 : IF (dft_control%qs_control%do_kg) THEN
1580 18 : IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed .OR. &
1581 : qs_env%kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
1582 8 : exc = exc - ekin_mol
1583 : END IF
1584 : END IF
1585 :
1586 172 : IF (debug_stress) THEN
1587 0 : stdeb = -1.0_dp*fconv*exc
1588 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1589 0 : 'STRESS| VOL 1st eps_XC[n_in]*n_z', one_third_sum_diag(stdeb), det_3x3(stdeb)
1590 : END IF
1591 :
1592 : ELSE ! use_virial
1593 :
1594 : ! calculate associated hartree potential
1595 : ! contribution for both T and D^Z
1596 974 : IF (gapw) THEN
1597 176 : CALL pw_axpy(local_rho_set_t%rho0_mpole%rho0_s_gs, rhoz_tot_gspace)
1598 176 : IF (ASSOCIATED(local_rho_set_t%rho0_mpole%rhoz_cneo_s_gs)) THEN
1599 0 : CALL pw_axpy(local_rho_set_t%rho0_mpole%rhoz_cneo_s_gs, rhoz_tot_gspace)
1600 : END IF
1601 : END IF
1602 974 : CALL pw_poisson_solve(poisson_env, rhoz_tot_gspace, ehartree, zv_hartree_gspace)
1603 :
1604 : END IF ! use virial
1605 1146 : IF (gapw .OR. gapw_xc) THEN
1606 216 : IF (ASSOCIATED(local_rho_set_t)) CALL local_rho_set_release(local_rho_set_t)
1607 : END IF
1608 :
1609 1644 : IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
1610 1146 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_ehartree
1611 1146 : CALL pw_transfer(zv_hartree_gspace, zv_hartree_rspace)
1612 1146 : CALL pw_scale(zv_hartree_rspace, zv_hartree_rspace%pw_grid%dvol)
1613 : ! Getting nuclear force contribution from the core charge density (not for GAPW)
1614 1146 : CALL integrate_v_core_rspace(zv_hartree_rspace, qs_env)
1615 1146 : IF (debug_forces) THEN
1616 664 : fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
1617 166 : CALL para_env%sum(fodeb)
1618 166 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vh(rhoz)*dncore ", fodeb
1619 : END IF
1620 1146 : IF (debug_stress .AND. use_virial) THEN
1621 0 : stdeb = fconv*(virial%pv_ehartree - stdeb)
1622 0 : CALL para_env%sum(stdeb)
1623 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1624 0 : 'STRESS| INT Vh(rhoz)*dncore ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1625 : END IF
1626 :
1627 : !
1628 1146 : IF (gapw_xc) THEN
1629 40 : CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
1630 : ELSE
1631 1106 : CALL get_qs_env(qs_env=qs_env, rho=rho)
1632 : END IF
1633 1146 : IF (dft_control%do_admm) THEN
1634 256 : CALL get_qs_env(qs_env, admm_env=admm_env)
1635 256 : xc_section => admm_env%xc_section_primary
1636 : ELSE
1637 890 : xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
1638 : END IF
1639 :
1640 1146 : IF (use_virial) THEN
1641 2236 : virial%pv_xc = 0.0_dp
1642 : END IF
1643 :
1644 1146 : IF (gapw .OR. gapw_xc) THEN
1645 : !get local_rho_set for GS density and response potential / density
1646 216 : NULLIFY (local_rho_set_t)
1647 216 : CALL local_rho_set_create(local_rho_set_t)
1648 : CALL allocate_rho_atom_internals(local_rho_set_t%rho_atom_set, atomic_kind_set, &
1649 216 : qs_kind_set, dft_control, para_env)
1650 : CALL init_rho0(local_rho_set_t, qs_env, dft_control%qs_control%gapw_control, &
1651 216 : zcore=0.0_dp)
1652 216 : CALL rho0_s_grid_create(pw_env, local_rho_set_t%rho0_mpole)
1653 : CALL calculate_rho_atom_coeff(qs_env, mpa(:), local_rho_set_t%rho_atom_set, &
1654 216 : qs_kind_set, oce, sab_orb, para_env)
1655 216 : CALL prepare_gapw_den(qs_env, local_rho_set_t, do_rho0=gapw)
1656 216 : NULLIFY (local_rho_set_gs)
1657 216 : CALL local_rho_set_create(local_rho_set_gs)
1658 : CALL allocate_rho_atom_internals(local_rho_set_gs%rho_atom_set, atomic_kind_set, &
1659 216 : qs_kind_set, dft_control, para_env)
1660 216 : CALL init_rho0(local_rho_set_gs, qs_env, dft_control%qs_control%gapw_control)
1661 216 : CALL rho0_s_grid_create(pw_env, local_rho_set_gs%rho0_mpole)
1662 : CALL calculate_rho_atom_coeff(qs_env, matrix_p(:, 1), local_rho_set_gs%rho_atom_set, &
1663 216 : qs_kind_set, oce, sab_orb, para_env)
1664 216 : CALL prepare_gapw_den(qs_env, local_rho_set_gs, do_rho0=gapw)
1665 : ! compute response potential
1666 1084 : ALLOCATE (rho_r_t(nspins), rho_g_t(nspins))
1667 434 : DO ispin = 1, nspins
1668 218 : CALL auxbas_pw_pool%create_pw(rho_r_t(ispin))
1669 434 : CALL auxbas_pw_pool%create_pw(rho_g_t(ispin))
1670 : END DO
1671 216 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace_t)
1672 216 : total_rho_t = 0.0_dp
1673 216 : CALL pw_zero(rho_tot_gspace_t)
1674 434 : DO ispin = 1, nspins
1675 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpa(ispin)%matrix, &
1676 : rho=rho_r_t(ispin), &
1677 : rho_gspace=rho_g_t(ispin), &
1678 : soft_valid=gapw, &
1679 218 : total_rho=total_rho_t(ispin))
1680 434 : CALL pw_axpy(rho_g_t(ispin), rho_tot_gspace_t)
1681 : END DO
1682 : ! add rho0 contributions to response density (only for Coulomb) only for gapw
1683 216 : IF (gapw) THEN
1684 176 : CALL pw_axpy(local_rho_set_t%rho0_mpole%rho0_s_gs, rho_tot_gspace_t)
1685 176 : IF (ASSOCIATED(local_rho_set_t%rho0_mpole%rhoz_cneo_s_gs)) THEN
1686 0 : CALL pw_axpy(local_rho_set_t%rho0_mpole%rhoz_cneo_s_gs, rho_tot_gspace_t)
1687 : END IF
1688 : ! compute response Coulomb potential
1689 176 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace_t)
1690 176 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace_t)
1691 176 : NULLIFY (hartree_local_t)
1692 176 : CALL hartree_local_create(hartree_local_t)
1693 176 : CALL init_coulomb_local(hartree_local_t, natom)
1694 176 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace_t, hartree_t, v_hartree_gspace_t)
1695 176 : CALL pw_transfer(v_hartree_gspace_t, v_hartree_rspace_t)
1696 176 : CALL pw_scale(v_hartree_rspace_t, v_hartree_rspace_t%pw_grid%dvol)
1697 : !
1698 584 : IF (debug_forces) fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1)
1699 : CALL Vh_1c_gg_integrals(qs_env, hartree_t, hartree_local_t%ecoul_1c, local_rho_set_gs, para_env, tddft=.FALSE., &
1700 176 : local_rho_set_2nd=local_rho_set_t, core_2nd=.TRUE.) ! n^core for GS potential
1701 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace_t, para_env, calculate_forces=.TRUE., &
1702 176 : local_rho_set=local_rho_set_gs, local_rho_set_2nd=local_rho_set_t)
1703 176 : IF (debug_forces) THEN
1704 544 : fodeb(1:3) = force(1)%g0s_Vh_elec(1:3, 1) - fodeb(1:3)
1705 136 : CALL para_env%sum(fodeb)
1706 136 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vh(T)*dncore PAWg0", fodeb
1707 : END IF
1708 : END IF !gapw
1709 : END IF !gapw
1710 :
1711 1146 : do_onecenter = .FALSE.
1712 1146 : NULLIFY (rho0_atom_set, rho1_atom_set)
1713 1146 : IF (gapw .OR. gapw_xc) THEN
1714 : !GAPW compute atomic fxc contributions
1715 216 : IF (myfun /= xc_none) THEN
1716 : ! local_rho_set_f
1717 192 : NULLIFY (local_rho_set_f)
1718 192 : CALL local_rho_set_create(local_rho_set_f)
1719 : CALL allocate_rho_atom_internals(local_rho_set_f%rho_atom_set, atomic_kind_set, &
1720 192 : qs_kind_set, dft_control, para_env)
1721 : CALL calculate_rho_atom_coeff(qs_env, mpa, local_rho_set_f%rho_atom_set, &
1722 192 : qs_kind_set, oce, sab_orb, para_env)
1723 192 : CALL prepare_gapw_den(qs_env, local_rho_set_f, do_rho0=.FALSE.)
1724 192 : rho0_atom_set => local_rho_set_gs%rho_atom_set
1725 192 : rho1_atom_set => local_rho_set_f%rho_atom_set
1726 192 : do_onecenter = .TRUE.
1727 : END IF ! myfun
1728 : END IF
1729 :
1730 1146 : NULLIFY (v_xc, v_xc_tau)
1731 : CALL qs_fxc_create(qs_env, rho0, rho1, rho0_atom_set, xc_section, do_onecenter, &
1732 : v_xc, v_xc_tau, rho1_atom_set, &
1733 1146 : compute_virial=use_virial, virial_xc=virial%pv_xc)
1734 1146 : DEALLOCATE (rho1)
1735 :
1736 : ! Stress-tensor XC-kernel GGA contribution
1737 1146 : IF (use_virial) THEN
1738 2236 : virial%pv_exc = virial%pv_exc + virial%pv_xc
1739 2236 : virial%pv_virial = virial%pv_virial + virial%pv_xc
1740 : END IF
1741 :
1742 1146 : IF (debug_stress .AND. use_virial) THEN
1743 0 : stdeb = 1.0_dp*fconv*virial%pv_xc
1744 0 : CALL para_env%sum(stdeb)
1745 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1746 0 : 'STRESS| GGA 2nd Pin*dK*rhoz', one_third_sum_diag(stdeb), det_3x3(stdeb)
1747 : END IF
1748 :
1749 : ! Stress-tensor integral contribution of 2nd derivative terms
1750 1146 : IF (use_virial) THEN
1751 2236 : pv_loc = virial%pv_virial
1752 : END IF
1753 :
1754 1146 : CALL get_qs_env(qs_env=qs_env, rho=rho)
1755 1146 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1756 1146 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1757 :
1758 2402 : DO ispin = 1, nspins
1759 2402 : CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
1760 : END DO
1761 1146 : IF ((.NOT. (gapw)) .AND. (.NOT. gapw_xc)) THEN
1762 942 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1763 1968 : DO ispin = 1, nspins
1764 1038 : CALL pw_axpy(zv_hartree_rspace, v_xc(ispin)) ! Hartree potential of response density
1765 : CALL integrate_v_rspace(qs_env=qs_env, &
1766 : v_rspace=v_xc(ispin), &
1767 : hmat=matrix_hz(ispin), &
1768 : pmat=matrix_p(ispin, 1), &
1769 : gapw=.FALSE., &
1770 1968 : calculate_forces=.TRUE.)
1771 : END DO
1772 930 : IF (debug_forces) THEN
1773 16 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1774 4 : CALL para_env%sum(fodeb)
1775 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKhxc*rhoz ", fodeb
1776 : END IF
1777 : ELSE
1778 702 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1779 216 : IF (myfun /= xc_none) THEN
1780 386 : DO ispin = 1, nspins
1781 : CALL integrate_v_rspace(qs_env=qs_env, &
1782 : v_rspace=v_xc(ispin), &
1783 : hmat=matrix_hz(ispin), &
1784 : pmat=matrix_p(ispin, 1), &
1785 : gapw=.TRUE., &
1786 386 : calculate_forces=.TRUE.)
1787 : END DO
1788 : END IF ! my_fun
1789 : ! Coulomb T+Dz
1790 434 : DO ispin = 1, nspins
1791 218 : CALL pw_zero(v_xc(ispin))
1792 218 : IF (gapw) THEN ! Hartree potential of response density
1793 178 : CALL pw_axpy(v_hartree_rspace_t, v_xc(ispin))
1794 40 : ELSE IF (gapw_xc) THEN
1795 40 : CALL pw_axpy(zv_hartree_rspace, v_xc(ispin))
1796 : END IF
1797 : CALL integrate_v_rspace(qs_env=qs_env, &
1798 : v_rspace=v_xc(ispin), &
1799 : hmat=matrix_ht(ispin), &
1800 : pmat=matrix_p(ispin, 1), &
1801 : gapw=gapw, &
1802 434 : calculate_forces=.TRUE.)
1803 : END DO
1804 216 : IF (debug_forces) THEN
1805 648 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1806 162 : CALL para_env%sum(fodeb)
1807 162 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKhxc*rhoz ", fodeb
1808 : END IF
1809 : END IF
1810 :
1811 1146 : IF (gapw .OR. gapw_xc) THEN
1812 : ! compute hard and soft atomic contributions
1813 216 : IF (myfun /= xc_none) THEN
1814 606 : IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
1815 : CALL update_ks_atom(qs_env, matrix_hz, matrix_p, forces=.TRUE., tddft=.FALSE., &
1816 192 : rho_atom_external=local_rho_set_f%rho_atom_set)
1817 192 : IF (debug_forces) THEN
1818 552 : fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
1819 138 : CALL para_env%sum(fodeb)
1820 138 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P^GS*dKxc*(Dz+T) PAW", fodeb
1821 : END IF
1822 : END IF !myfun
1823 : ! Coulomb contributions
1824 216 : IF (gapw) THEN
1825 584 : IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
1826 : CALL update_ks_atom(qs_env, matrix_ht, matrix_p, forces=.TRUE., tddft=.FALSE., &
1827 176 : rho_atom_external=local_rho_set_t%rho_atom_set)
1828 176 : IF (debug_forces) THEN
1829 544 : fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
1830 136 : CALL para_env%sum(fodeb)
1831 136 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P^GS*dKh*(Dz+T) PAW", fodeb
1832 : END IF
1833 : END IF
1834 : ! add Coulomb and XC
1835 434 : DO ispin = 1, nspins
1836 434 : CALL dbcsr_add(matrix_hz(ispin)%matrix, matrix_ht(ispin)%matrix, 1.0_dp, 1.0_dp)
1837 : END DO
1838 :
1839 : ! release
1840 216 : IF (myfun /= xc_none) THEN
1841 192 : IF (ASSOCIATED(local_rho_set_f)) CALL local_rho_set_release(local_rho_set_f)
1842 : END IF
1843 216 : IF (ASSOCIATED(local_rho_set_t)) CALL local_rho_set_release(local_rho_set_t)
1844 216 : IF (ASSOCIATED(local_rho_set_gs)) CALL local_rho_set_release(local_rho_set_gs)
1845 216 : IF (gapw) THEN
1846 176 : IF (ASSOCIATED(hartree_local_t)) CALL hartree_local_release(hartree_local_t)
1847 176 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace_t)
1848 176 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace_t)
1849 : END IF
1850 216 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace_t)
1851 434 : DO ispin = 1, nspins
1852 218 : CALL auxbas_pw_pool%give_back_pw(rho_r_t(ispin))
1853 434 : CALL auxbas_pw_pool%give_back_pw(rho_g_t(ispin))
1854 : END DO
1855 216 : DEALLOCATE (rho_r_t, rho_g_t)
1856 : END IF ! gapw
1857 :
1858 1146 : IF (debug_stress .AND. use_virial) THEN
1859 0 : stdeb = fconv*(virial%pv_virial - stdeb)
1860 0 : CALL para_env%sum(stdeb)
1861 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1862 0 : 'STRESS| INT 2nd f_Hxc[Pz]*Pin', one_third_sum_diag(stdeb), det_3x3(stdeb)
1863 : END IF
1864 : !
1865 1146 : IF (ASSOCIATED(v_xc_tau)) THEN
1866 32 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1867 32 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
1868 64 : DO ispin = 1, nspins
1869 32 : CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
1870 : CALL integrate_v_rspace(qs_env=qs_env, &
1871 : v_rspace=v_xc_tau(ispin), &
1872 : hmat=matrix_hz(ispin), &
1873 : pmat=matrix_p(ispin, 1), &
1874 : compute_tau=.TRUE., &
1875 : gapw=(gapw .OR. gapw_xc), &
1876 96 : calculate_forces=.TRUE.)
1877 : END DO
1878 32 : IF (debug_forces) THEN
1879 0 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1880 0 : CALL para_env%sum(fodeb)
1881 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dKtau*tauz ", fodeb
1882 : END IF
1883 : END IF
1884 1146 : IF (debug_stress .AND. use_virial) THEN
1885 0 : stdeb = fconv*(virial%pv_virial - stdeb)
1886 0 : CALL para_env%sum(stdeb)
1887 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1888 0 : 'STRESS| INT 2nd f_xctau[Pz]*Pin', one_third_sum_diag(stdeb), det_3x3(stdeb)
1889 : END IF
1890 : ! Stress-tensor integral contribution of 2nd derivative terms
1891 1146 : IF (use_virial) THEN
1892 2236 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1893 : END IF
1894 :
1895 : ! KG Embedding
1896 : ! calculate kinetic energy kernel, folded with response density for partial integration
1897 1146 : IF (dft_control%qs_control%do_kg) THEN
1898 24 : IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed) THEN
1899 : ekin_mol = 0.0_dp
1900 12 : IF (use_virial) THEN
1901 104 : pv_loc = virial%pv_virial
1902 : END IF
1903 :
1904 12 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
1905 108 : IF (use_virial) virial%pv_xc = 0.0_dp
1906 : CALL kg_ekin_subset(qs_env=qs_env, &
1907 : ks_matrix=matrix_hz, &
1908 : ekin_mol=ekin_mol, &
1909 : calc_force=.TRUE., &
1910 : do_kernel=.TRUE., &
1911 12 : pmat_ext=matrix_pz)
1912 :
1913 12 : IF (debug_forces) THEN
1914 0 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
1915 0 : CALL para_env%sum(fodeb)
1916 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*d(Kkg)*rhoz ", fodeb
1917 : END IF
1918 12 : IF (debug_stress .AND. use_virial) THEN
1919 0 : stdeb = fconv*(virial%pv_virial - pv_loc)
1920 0 : CALL para_env%sum(stdeb)
1921 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1922 0 : 'STRESS| INT KG Pin*d(KKG)*rhoz ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1923 :
1924 0 : stdeb = fconv*(virial%pv_xc)
1925 0 : CALL para_env%sum(stdeb)
1926 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
1927 0 : 'STRESS| GGA KG Pin*d(KKG)*rhoz ', one_third_sum_diag(stdeb), det_3x3(stdeb)
1928 : END IF
1929 :
1930 : ! Stress tensor
1931 12 : IF (use_virial) THEN
1932 : ! XC-kernel Integral contribution
1933 104 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1934 :
1935 : ! XC-kernel GGA contribution
1936 104 : virial%pv_exc = virial%pv_exc - virial%pv_xc
1937 104 : virial%pv_virial = virial%pv_virial - virial%pv_xc
1938 104 : virial%pv_xc = 0.0_dp
1939 : END IF
1940 : END IF
1941 : END IF
1942 1146 : CALL auxbas_pw_pool%give_back_pw(rhoz_tot_gspace)
1943 1146 : CALL auxbas_pw_pool%give_back_pw(zv_hartree_gspace)
1944 1146 : CALL auxbas_pw_pool%give_back_pw(zv_hartree_rspace)
1945 2402 : DO ispin = 1, nspins
1946 1256 : CALL auxbas_pw_pool%give_back_pw(rhoz_r(ispin))
1947 1256 : CALL auxbas_pw_pool%give_back_pw(rhoz_g(ispin))
1948 2402 : CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
1949 : END DO
1950 1146 : DEALLOCATE (rhoz_r, rhoz_g, v_xc)
1951 1146 : IF (gapw_xc) THEN
1952 80 : DO ispin = 1, nspins
1953 40 : CALL auxbas_pw_pool%give_back_pw(rhoz_r_xc(ispin))
1954 80 : CALL auxbas_pw_pool%give_back_pw(rhoz_g_xc(ispin))
1955 : END DO
1956 40 : DEALLOCATE (rhoz_r_xc, rhoz_g_xc)
1957 : END IF
1958 1146 : IF (ASSOCIATED(v_xc_tau)) THEN
1959 64 : DO ispin = 1, nspins
1960 32 : CALL auxbas_pw_pool%give_back_pw(tauz_r(ispin))
1961 64 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
1962 : END DO
1963 32 : DEALLOCATE (tauz_r, v_xc_tau)
1964 32 : IF (ASSOCIATED(tauz_r_xc)) THEN
1965 0 : DO ispin = 1, nspins
1966 0 : CALL auxbas_pw_pool%give_back_pw(tauz_r_xc(ispin))
1967 : END DO
1968 0 : DEALLOCATE (tauz_r_xc)
1969 : END IF
1970 : END IF
1971 1146 : IF (debug_forces) THEN
1972 498 : ALLOCATE (ftot3(3, natom))
1973 166 : CALL total_qs_force(ftot3, force, atomic_kind_set)
1974 664 : fodeb(1:3) = ftot3(1:3, 1) - ftot2(1:3, 1)
1975 166 : CALL para_env%sum(fodeb)
1976 166 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*V(rhoz)", fodeb
1977 : END IF
1978 1146 : CALL dbcsr_deallocate_matrix_set(scrm)
1979 1146 : CALL dbcsr_deallocate_matrix_set(matrix_ht)
1980 :
1981 : ! -----------------------------------------
1982 : ! Apply ADMM exchange correction
1983 : ! -----------------------------------------
1984 :
1985 1146 : IF (dft_control%do_admm) THEN
1986 : ! volume term
1987 256 : exc_aux_fit = 0.0_dp
1988 :
1989 256 : IF (qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
1990 : ! nothing to do
1991 112 : NULLIFY (mpz, mhz, mhx, mhy)
1992 : ELSE
1993 : ! add ADMM xc_section_aux terms: Pz*Vxc + P0*K0[rhoz]
1994 144 : CALL get_qs_env(qs_env, admm_env=admm_env)
1995 : CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit, matrix_s_aux_fit=scrm, &
1996 144 : task_list_aux_fit=task_list_aux_fit)
1997 : !
1998 144 : NULLIFY (mpz, mhz, mhx, mhy)
1999 144 : CALL dbcsr_allocate_matrix_set(mhx, nspins, 1)
2000 144 : CALL dbcsr_allocate_matrix_set(mhy, nspins, 1)
2001 144 : CALL dbcsr_allocate_matrix_set(mpz, nspins, 1)
2002 296 : DO ispin = 1, nspins
2003 152 : ALLOCATE (mhx(ispin, 1)%matrix)
2004 152 : CALL dbcsr_create(mhx(ispin, 1)%matrix, template=scrm(1)%matrix)
2005 152 : CALL dbcsr_copy(mhx(ispin, 1)%matrix, scrm(1)%matrix)
2006 152 : CALL dbcsr_set(mhx(ispin, 1)%matrix, 0.0_dp)
2007 152 : ALLOCATE (mhy(ispin, 1)%matrix)
2008 152 : CALL dbcsr_create(mhy(ispin, 1)%matrix, template=scrm(1)%matrix)
2009 152 : CALL dbcsr_copy(mhy(ispin, 1)%matrix, scrm(1)%matrix)
2010 152 : CALL dbcsr_set(mhy(ispin, 1)%matrix, 0.0_dp)
2011 152 : ALLOCATE (mpz(ispin, 1)%matrix)
2012 296 : IF (do_ex) THEN
2013 94 : CALL dbcsr_create(mpz(ispin, 1)%matrix, template=p_env%p1_admm(ispin)%matrix)
2014 94 : CALL dbcsr_copy(mpz(ispin, 1)%matrix, p_env%p1_admm(ispin)%matrix)
2015 : CALL dbcsr_add(mpz(ispin, 1)%matrix, ex_env%matrix_pe_admm(ispin)%matrix, &
2016 94 : 1.0_dp, 1.0_dp)
2017 : ELSE
2018 58 : CALL dbcsr_create(mpz(ispin, 1)%matrix, template=matrix_pz_admm(ispin)%matrix)
2019 58 : CALL dbcsr_copy(mpz(ispin, 1)%matrix, matrix_pz_admm(ispin)%matrix)
2020 : END IF
2021 : END DO
2022 : !
2023 144 : xc_section => admm_env%xc_section_aux
2024 144 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
2025 144 : needs = xc_functionals_get_needs(xc_fun_section, (nspins == 2), .TRUE.)
2026 144 : needs_tau_response_aux = needs%tau .OR. needs%tau_spin
2027 : ! Stress-tensor: integration contribution direct term
2028 : ! int Pz*v_xc[rho_admm]
2029 144 : IF (use_virial) THEN
2030 260 : pv_loc = virial%pv_virial
2031 : END IF
2032 :
2033 144 : basis_type = "AUX_FIT"
2034 144 : task_list => task_list_aux_fit
2035 144 : IF (admm_env%do_gapw) THEN
2036 14 : basis_type = "AUX_FIT_SOFT"
2037 14 : task_list => admm_env%admm_gapw_env%task_list
2038 : END IF
2039 : !
2040 180 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2041 144 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2042 296 : DO ispin = 1, nspins
2043 : CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), &
2044 : hmat=mhx(ispin, 1), pmat=mpz(ispin, 1), &
2045 : qs_env=qs_env, calculate_forces=.TRUE., &
2046 152 : basis_type=basis_type, task_list_external=task_list)
2047 296 : IF (ASSOCIATED(vadmm_tau_rspace)) THEN
2048 : CALL integrate_v_rspace(v_rspace=vadmm_tau_rspace(ispin), &
2049 : hmat=mhx(ispin, 1), pmat=mpz(ispin, 1), &
2050 : qs_env=qs_env, calculate_forces=.TRUE., compute_tau=.TRUE., &
2051 0 : basis_type=basis_type, task_list_external=task_list)
2052 : END IF
2053 : END DO
2054 144 : IF (debug_forces) THEN
2055 48 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2056 12 : CALL para_env%sum(fodeb)
2057 12 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*Vxc(rho_admm)", fodeb
2058 : END IF
2059 144 : IF (debug_stress .AND. use_virial) THEN
2060 0 : stdeb = fconv*(virial%pv_virial - pv_loc)
2061 0 : CALL para_env%sum(stdeb)
2062 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
2063 0 : 'STRESS| INT 1st Pz*dVxc(rho_admm) ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2064 : END IF
2065 : ! Stress-tensor Pz_admm*v_xc[rho_admm]
2066 144 : IF (use_virial) THEN
2067 260 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2068 : END IF
2069 : !
2070 144 : IF (admm_env%do_gapw) THEN
2071 14 : CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
2072 50 : IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
2073 : CALL update_ks_atom(qs_env, mhx(:, 1), mpz(:, 1), forces=.TRUE., tddft=.FALSE., &
2074 : rho_atom_external=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
2075 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
2076 : oce_external=admm_env%admm_gapw_env%oce, &
2077 14 : sab_external=sab_aux_fit)
2078 14 : IF (debug_forces) THEN
2079 48 : fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
2080 12 : CALL para_env%sum(fodeb)
2081 12 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*Vxc(rho_admm)PAW", fodeb
2082 : END IF
2083 : END IF
2084 : !
2085 : ! rhoz_aux
2086 144 : NULLIFY (rhoz_g_aux, rhoz_r_aux, rhoz_tau_r_aux)
2087 1024 : ALLOCATE (rhoz_r_aux(nspins), rhoz_g_aux(nspins))
2088 296 : DO ispin = 1, nspins
2089 152 : CALL auxbas_pw_pool%create_pw(rhoz_r_aux(ispin))
2090 296 : CALL auxbas_pw_pool%create_pw(rhoz_g_aux(ispin))
2091 : END DO
2092 296 : DO ispin = 1, nspins
2093 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpz(ispin, 1)%matrix, &
2094 : rho=rhoz_r_aux(ispin), rho_gspace=rhoz_g_aux(ispin), &
2095 296 : basis_type=basis_type, task_list_external=task_list)
2096 : END DO
2097 144 : IF (needs_tau_response_aux .OR. ASSOCIATED(vadmm_tau_rspace)) THEN
2098 : BLOCK
2099 : TYPE(pw_c1d_gs_type) :: work_g
2100 0 : ALLOCATE (rhoz_tau_r_aux(nspins))
2101 0 : CALL auxbas_pw_pool%create_pw(work_g)
2102 0 : DO ispin = 1, nspins
2103 0 : CALL auxbas_pw_pool%create_pw(rhoz_tau_r_aux(ispin))
2104 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpz(ispin, 1)%matrix, &
2105 : rho=rhoz_tau_r_aux(ispin), rho_gspace=work_g, &
2106 : basis_type=basis_type, task_list_external=task_list, &
2107 0 : compute_tau=.TRUE.)
2108 : END DO
2109 0 : CALL auxbas_pw_pool%give_back_pw(work_g)
2110 : END BLOCK
2111 : END IF
2112 : !
2113 : ! Add ADMM volume contribution to stress tensor
2114 144 : IF (use_virial) THEN
2115 :
2116 : ! Stress tensor volume term: \int v_xc[n_in_admm]*n_z_admm
2117 : ! vadmm_rspace already scaled, we need to unscale it!
2118 40 : DO ispin = 1, nspins
2119 : exc_aux_fit = exc_aux_fit + pw_integral_ab(rhoz_r_aux(ispin), vadmm_rspace(ispin))/ &
2120 40 : vadmm_rspace(ispin)%pw_grid%dvol
2121 : END DO
2122 20 : IF (ASSOCIATED(vadmm_tau_rspace) .AND. ASSOCIATED(rhoz_tau_r_aux)) THEN
2123 0 : DO ispin = 1, nspins
2124 : exc_aux_fit = exc_aux_fit + pw_integral_ab(rhoz_tau_r_aux(ispin), vadmm_tau_rspace(ispin))/ &
2125 0 : vadmm_tau_rspace(ispin)%pw_grid%dvol
2126 : END DO
2127 : END IF
2128 :
2129 20 : IF (debug_stress) THEN
2130 0 : stdeb = -1.0_dp*fconv*exc_aux_fit
2131 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T43,2(1X,ES19.11))") &
2132 0 : 'STRESS| VOL 1st eps_XC[n_in_admm]*n_z_admm', one_third_sum_diag(stdeb), det_3x3(stdeb)
2133 : END IF
2134 :
2135 : END IF
2136 : !
2137 144 : NULLIFY (v_xc, v_xc_tau)
2138 :
2139 384 : IF (use_virial) virial%pv_xc = 0.0_dp
2140 :
2141 144 : NULLIFY (rho0_atom_set, rho1_atom_set)
2142 144 : kind_set => qs_kind_set
2143 144 : IF (admm_env%do_gapw) THEN
2144 14 : kind_set => admm_env%admm_gapw_env%admm_kind_set
2145 14 : CALL local_rho_set_create(local_rhoz_set_admm)
2146 : CALL allocate_rho_atom_internals(local_rhoz_set_admm%rho_atom_set, atomic_kind_set, &
2147 14 : kind_set, dft_control, para_env)
2148 : CALL calculate_rho_atom_coeff(qs_env, mpz(:, 1), local_rhoz_set_admm%rho_atom_set, &
2149 14 : kind_set, admm_env%admm_gapw_env%oce, sab_aux_fit, para_env)
2150 : CALL prepare_gapw_den(qs_env, local_rho_set=local_rhoz_set_admm, &
2151 14 : do_rho0=.FALSE., kind_set_external=kind_set)
2152 14 : rho0_atom_set => admm_env%admm_gapw_env%local_rho_set%rho_atom_set
2153 14 : rho1_atom_set => local_rhoz_set_admm%rho_atom_set
2154 14 : do_onecenter = .TRUE.
2155 : END IF
2156 :
2157 144 : ALLOCATE (rho1)
2158 144 : CALL qs_rho_create(rho1)
2159 144 : IF (ASSOCIATED(rhoz_tau_r_aux)) THEN
2160 : CALL qs_rho_set(rho1, rho_r=rhoz_r_aux, rho_g=rhoz_g_aux, tau_r=rhoz_tau_r_aux, &
2161 0 : rho_r_valid=.TRUE., rho_g_valid=.TRUE., tau_r_valid=.TRUE.)
2162 : ELSE
2163 : CALL qs_rho_set(rho1, rho_r=rhoz_r_aux, rho_g=rhoz_g_aux, &
2164 144 : rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
2165 : END IF
2166 : CALL qs_fxc_create(qs_env, rho_aux_fit, rho1, rho0_atom_set, xc_section, do_onecenter, &
2167 : v_xc, v_xc_tau, rho1_atom_set, &
2168 : kind_set_external=kind_set, &
2169 144 : compute_virial=use_virial, virial_xc=virial%pv_xc)
2170 :
2171 : ! Stress-tensor ADMM-kernel GGA contribution
2172 144 : IF (use_virial) THEN
2173 260 : virial%pv_exc = virial%pv_exc + virial%pv_xc
2174 260 : virial%pv_virial = virial%pv_virial + virial%pv_xc
2175 : END IF
2176 :
2177 144 : IF (debug_stress .AND. use_virial) THEN
2178 0 : stdeb = 1.0_dp*fconv*virial%pv_xc
2179 0 : CALL para_env%sum(stdeb)
2180 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
2181 0 : 'STRESS| GGA 2nd Pin_admm*dK*rhoz_admm', one_third_sum_diag(stdeb), det_3x3(stdeb)
2182 : END IF
2183 : !
2184 144 : CALL qs_rho_get(rho_aux_fit, rho_ao_kp=matrix_p)
2185 : ! Stress-tensor Pin*dK*rhoz_admm
2186 144 : IF (use_virial) THEN
2187 260 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2188 : END IF
2189 180 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2190 144 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2191 296 : DO ispin = 1, nspins
2192 152 : CALL dbcsr_set(mhy(ispin, 1)%matrix, 0.0_dp)
2193 152 : CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2194 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
2195 : hmat=mhy(ispin, 1), pmat=matrix_p(ispin, 1), &
2196 : calculate_forces=.TRUE., &
2197 296 : basis_type=basis_type, task_list_external=task_list)
2198 : END DO
2199 144 : IF (ASSOCIATED(v_xc_tau)) THEN
2200 0 : DO ispin = 1, nspins
2201 0 : CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2202 : CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc_tau(ispin), &
2203 : hmat=mhy(ispin, 1), pmat=matrix_p(ispin, 1), &
2204 : calculate_forces=.TRUE., compute_tau=.TRUE., &
2205 0 : basis_type=basis_type, task_list_external=task_list)
2206 : END DO
2207 : END IF
2208 144 : IF (debug_forces) THEN
2209 48 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2210 12 : CALL para_env%sum(fodeb)
2211 12 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dK*rhoz_admm ", fodeb
2212 : END IF
2213 144 : IF (debug_stress .AND. use_virial) THEN
2214 0 : stdeb = fconv*(virial%pv_virial - pv_loc)
2215 0 : CALL para_env%sum(stdeb)
2216 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
2217 0 : 'STRESS| INT 2nd Pin*dK*rhoz_admm ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2218 : END IF
2219 : ! Stress-tensor Pin*dK*rhoz_admm
2220 144 : IF (use_virial) THEN
2221 260 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
2222 : END IF
2223 : ! GAPW ADMM XC correction integrate weight contribution to force
2224 144 : IF (admm_env%do_gapw) THEN
2225 50 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
2226 14 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_virial
2227 : !
2228 14 : CALL accint_weight_force(qs_env, rho_aux_fit, rho1, 1, xc_section)
2229 : !
2230 14 : IF (debug_forces) THEN
2231 48 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
2232 12 : CALL para_env%sum(fodeb)
2233 12 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: dKxc*rhoz_admm*dw ", fodeb
2234 : END IF
2235 14 : IF (debug_stress .AND. use_virial) THEN
2236 0 : stdeb = fconv*(virial%pv_virial - stdeb)
2237 0 : CALL para_env%sum(stdeb)
2238 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
2239 0 : 'STRESS| dKxc*rhoz_admm*dw', one_third_sum_diag(stdeb), det_3x3(stdeb)
2240 : END IF
2241 : END IF
2242 : ! return ADMM response densities and potentials
2243 296 : DO ispin = 1, nspins
2244 152 : CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2245 152 : IF (ASSOCIATED(v_xc_tau)) CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2246 152 : CALL auxbas_pw_pool%give_back_pw(rhoz_r_aux(ispin))
2247 152 : CALL auxbas_pw_pool%give_back_pw(rhoz_g_aux(ispin))
2248 296 : IF (ASSOCIATED(rhoz_tau_r_aux)) CALL auxbas_pw_pool%give_back_pw(rhoz_tau_r_aux(ispin))
2249 : END DO
2250 144 : DEALLOCATE (v_xc, rhoz_r_aux, rhoz_g_aux)
2251 144 : IF (ASSOCIATED(v_xc_tau)) DEALLOCATE (v_xc_tau)
2252 144 : IF (ASSOCIATED(rhoz_tau_r_aux)) DEALLOCATE (rhoz_tau_r_aux)
2253 144 : DEALLOCATE (rho1)
2254 : !
2255 144 : IF (admm_env%do_gapw) THEN
2256 50 : IF (debug_forces) fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1)
2257 : CALL update_ks_atom(qs_env, mhy(:, 1), matrix_p(:, 1), forces=.TRUE., tddft=.FALSE., &
2258 : rho_atom_external=local_rhoz_set_admm%rho_atom_set, &
2259 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
2260 : oce_external=admm_env%admm_gapw_env%oce, &
2261 14 : sab_external=sab_aux_fit)
2262 14 : IF (debug_forces) THEN
2263 48 : fodeb(1:3) = force(1)%Vhxc_atom(1:3, 1) - fodeb(1:3)
2264 12 : CALL para_env%sum(fodeb)
2265 12 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pin*dK*rhoz_admm[PAW] ", fodeb
2266 : END IF
2267 14 : CALL local_rho_set_release(local_rhoz_set_admm)
2268 : END IF
2269 : !
2270 144 : nao = admm_env%nao_orb
2271 144 : nao_aux = admm_env%nao_aux_fit
2272 144 : ALLOCATE (dbwork)
2273 144 : CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
2274 296 : DO ispin = 1, nspins
2275 : CALL cp_dbcsr_sm_fm_multiply(mhy(ispin, 1)%matrix, admm_env%A, &
2276 152 : admm_env%work_aux_orb, nao)
2277 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
2278 : 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
2279 152 : admm_env%work_orb_orb)
2280 152 : CALL dbcsr_copy(dbwork, matrix_hz(ispin)%matrix)
2281 152 : CALL dbcsr_set(dbwork, 0.0_dp)
2282 152 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.TRUE.)
2283 296 : CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
2284 : END DO
2285 144 : CALL dbcsr_release(dbwork)
2286 144 : DEALLOCATE (dbwork)
2287 288 : CALL dbcsr_deallocate_matrix_set(mpz)
2288 : END IF ! qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none
2289 : END IF ! do_admm
2290 :
2291 : ! -----------------------------------------
2292 : ! HFX
2293 : ! -----------------------------------------
2294 :
2295 : ! HFX
2296 1146 : hfx_section => section_vals_get_subs_vals(xc_section, "HF")
2297 1146 : CALL section_vals_get(hfx_section, explicit=do_hfx)
2298 1146 : IF (do_hfx) THEN
2299 490 : CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
2300 490 : CPASSERT(n_rep_hf == 1)
2301 : CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
2302 490 : i_rep_section=1)
2303 490 : mspin = 1
2304 490 : IF (hfx_treat_lsd_in_core) mspin = nspins
2305 1306 : IF (use_virial) virial%pv_fock_4c = 0.0_dp
2306 : !
2307 : CALL get_qs_env(qs_env=qs_env, rho=rho, x_data=x_data, &
2308 490 : s_mstruct_changed=s_mstruct_changed)
2309 490 : distribute_fock_matrix = .TRUE.
2310 :
2311 : ! -----------------------------------------
2312 : ! HFX-ADMM
2313 : ! -----------------------------------------
2314 490 : IF (dft_control%do_admm) THEN
2315 256 : CALL get_qs_env(qs_env=qs_env, admm_env=admm_env)
2316 256 : CALL get_admm_env(admm_env, matrix_s_aux_fit=scrm, rho_aux_fit=rho_aux_fit)
2317 256 : CALL qs_rho_get(rho_aux_fit, rho_ao_kp=matrix_p)
2318 256 : NULLIFY (mpz, mhz, mpd, mhd)
2319 256 : CALL dbcsr_allocate_matrix_set(mpz, nspins, 1)
2320 256 : CALL dbcsr_allocate_matrix_set(mhz, nspins, 1)
2321 256 : CALL dbcsr_allocate_matrix_set(mpd, nspins, 1)
2322 256 : CALL dbcsr_allocate_matrix_set(mhd, nspins, 1)
2323 532 : DO ispin = 1, nspins
2324 276 : ALLOCATE (mhz(ispin, 1)%matrix, mhd(ispin, 1)%matrix)
2325 276 : CALL dbcsr_create(mhz(ispin, 1)%matrix, template=scrm(1)%matrix)
2326 276 : CALL dbcsr_create(mhd(ispin, 1)%matrix, template=scrm(1)%matrix)
2327 276 : CALL dbcsr_copy(mhz(ispin, 1)%matrix, scrm(1)%matrix)
2328 276 : CALL dbcsr_copy(mhd(ispin, 1)%matrix, scrm(1)%matrix)
2329 276 : CALL dbcsr_set(mhz(ispin, 1)%matrix, 0.0_dp)
2330 276 : CALL dbcsr_set(mhd(ispin, 1)%matrix, 0.0_dp)
2331 276 : ALLOCATE (mpz(ispin, 1)%matrix)
2332 276 : IF (do_ex) THEN
2333 162 : CALL dbcsr_create(mpz(ispin, 1)%matrix, template=scrm(1)%matrix)
2334 162 : CALL dbcsr_copy(mpz(ispin, 1)%matrix, p_env%p1_admm(ispin)%matrix)
2335 : CALL dbcsr_add(mpz(ispin, 1)%matrix, ex_env%matrix_pe_admm(ispin)%matrix, &
2336 162 : 1.0_dp, 1.0_dp)
2337 : ELSE
2338 114 : CALL dbcsr_create(mpz(ispin, 1)%matrix, template=scrm(1)%matrix)
2339 114 : CALL dbcsr_copy(mpz(ispin, 1)%matrix, matrix_pz_admm(ispin)%matrix)
2340 : END IF
2341 532 : mpd(ispin, 1)%matrix => matrix_p(ispin, 1)%matrix
2342 : END DO
2343 : !
2344 256 : IF (x_data(1, 1)%do_hfx_ri) THEN
2345 :
2346 : eh1 = 0.0_dp
2347 : CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpz, &
2348 : geometry_did_change=s_mstruct_changed, nspins=nspins, &
2349 6 : hf_fraction=x_data(1, 1)%general_parameter%fraction)
2350 :
2351 : eh1 = 0.0_dp
2352 : CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhd, eh1, rho_ao=mpd, &
2353 : geometry_did_change=s_mstruct_changed, nspins=nspins, &
2354 6 : hf_fraction=x_data(1, 1)%general_parameter%fraction)
2355 :
2356 : ELSE
2357 500 : DO ispin = 1, mspin
2358 : eh1 = 0.0
2359 : CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpz, hfx_section, &
2360 : para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
2361 500 : ispin=ispin)
2362 : END DO
2363 500 : DO ispin = 1, mspin
2364 : eh1 = 0.0
2365 : CALL integrate_four_center(qs_env, x_data, mhd, eh1, mpd, hfx_section, &
2366 : para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
2367 500 : ispin=ispin)
2368 : END DO
2369 : END IF
2370 : !
2371 256 : CALL get_qs_env(qs_env, admm_env=admm_env)
2372 256 : CPASSERT(ASSOCIATED(admm_env%work_aux_orb))
2373 256 : CPASSERT(ASSOCIATED(admm_env%work_orb_orb))
2374 256 : nao = admm_env%nao_orb
2375 256 : nao_aux = admm_env%nao_aux_fit
2376 256 : ALLOCATE (dbwork)
2377 256 : CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
2378 532 : DO ispin = 1, nspins
2379 : CALL cp_dbcsr_sm_fm_multiply(mhz(ispin, 1)%matrix, admm_env%A, &
2380 276 : admm_env%work_aux_orb, nao)
2381 : CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
2382 : 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
2383 276 : admm_env%work_orb_orb)
2384 276 : CALL dbcsr_copy(dbwork, matrix_hz(ispin)%matrix)
2385 276 : CALL dbcsr_set(dbwork, 0.0_dp)
2386 276 : CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.TRUE.)
2387 532 : CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
2388 : END DO
2389 256 : CALL dbcsr_release(dbwork)
2390 256 : DEALLOCATE (dbwork)
2391 : ! derivatives Tr (Pz [A(T)H dA/dR])
2392 328 : IF (debug_forces) fodeb(1:3) = force(1)%overlap_admm(1:3, 1)
2393 256 : IF (ASSOCIATED(mhx) .AND. ASSOCIATED(mhy)) THEN
2394 296 : DO ispin = 1, nspins
2395 152 : CALL dbcsr_add(mhd(ispin, 1)%matrix, mhx(ispin, 1)%matrix, 1.0_dp, 1.0_dp)
2396 296 : CALL dbcsr_add(mhz(ispin, 1)%matrix, mhy(ispin, 1)%matrix, 1.0_dp, 1.0_dp)
2397 : END DO
2398 : END IF
2399 256 : CALL qs_rho_get(rho, rho_ao=matrix_pd)
2400 256 : CALL admm_projection_derivative(qs_env, mhd(:, 1), mpa)
2401 256 : CALL admm_projection_derivative(qs_env, mhz(:, 1), matrix_pd)
2402 256 : IF (debug_forces) THEN
2403 96 : fodeb(1:3) = force(1)%overlap_admm(1:3, 1) - fodeb(1:3)
2404 24 : CALL para_env%sum(fodeb)
2405 24 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*hfx*S' ", fodeb
2406 : END IF
2407 256 : CALL dbcsr_deallocate_matrix_set(mpz)
2408 256 : CALL dbcsr_deallocate_matrix_set(mhz)
2409 256 : CALL dbcsr_deallocate_matrix_set(mhd)
2410 256 : IF (ASSOCIATED(mhx) .AND. ASSOCIATED(mhy)) THEN
2411 144 : CALL dbcsr_deallocate_matrix_set(mhx)
2412 144 : CALL dbcsr_deallocate_matrix_set(mhy)
2413 : END IF
2414 256 : DEALLOCATE (mpd)
2415 : ELSE
2416 : ! -----------------------------------------
2417 : ! conventional HFX
2418 : ! -----------------------------------------
2419 2154 : ALLOCATE (mpz(nspins, 1), mhz(nspins, 1))
2420 492 : DO ispin = 1, nspins
2421 258 : mhz(ispin, 1)%matrix => matrix_hz(ispin)%matrix
2422 492 : mpz(ispin, 1)%matrix => mpa(ispin)%matrix
2423 : END DO
2424 :
2425 234 : IF (x_data(1, 1)%do_hfx_ri) THEN
2426 :
2427 : eh1 = 0.0_dp
2428 : CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpz, &
2429 : geometry_did_change=s_mstruct_changed, nspins=nspins, &
2430 18 : hf_fraction=x_data(1, 1)%general_parameter%fraction)
2431 : ELSE
2432 432 : DO ispin = 1, mspin
2433 : eh1 = 0.0
2434 : CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpz, hfx_section, &
2435 : para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
2436 432 : ispin=ispin)
2437 : END DO
2438 : END IF
2439 234 : DEALLOCATE (mhz, mpz)
2440 : END IF
2441 :
2442 : ! -----------------------------------------
2443 : ! HFX FORCES
2444 : ! -----------------------------------------
2445 :
2446 490 : resp_only = .TRUE.
2447 676 : IF (debug_forces) fodeb(1:3) = force(1)%fock_4c(1:3, 1)
2448 490 : IF (dft_control%do_admm) THEN
2449 : ! -----------------------------------------
2450 : ! HFX-ADMM FORCES
2451 : ! -----------------------------------------
2452 256 : CALL qs_rho_get(rho_aux_fit, rho_ao_kp=matrix_p)
2453 256 : NULLIFY (matrix_pza)
2454 256 : CALL dbcsr_allocate_matrix_set(matrix_pza, nspins)
2455 532 : DO ispin = 1, nspins
2456 276 : ALLOCATE (matrix_pza(ispin)%matrix)
2457 532 : IF (do_ex) THEN
2458 162 : CALL dbcsr_create(matrix_pza(ispin)%matrix, template=p_env%p1_admm(ispin)%matrix)
2459 162 : CALL dbcsr_copy(matrix_pza(ispin)%matrix, p_env%p1_admm(ispin)%matrix)
2460 : CALL dbcsr_add(matrix_pza(ispin)%matrix, ex_env%matrix_pe_admm(ispin)%matrix, &
2461 162 : 1.0_dp, 1.0_dp)
2462 : ELSE
2463 114 : CALL dbcsr_create(matrix_pza(ispin)%matrix, template=matrix_pz_admm(ispin)%matrix)
2464 114 : CALL dbcsr_copy(matrix_pza(ispin)%matrix, matrix_pz_admm(ispin)%matrix)
2465 : END IF
2466 : END DO
2467 256 : IF (x_data(1, 1)%do_hfx_ri) THEN
2468 :
2469 : CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
2470 : x_data(1, 1)%general_parameter%fraction, &
2471 : rho_ao=matrix_p, rho_ao_resp=matrix_pza, &
2472 6 : use_virial=use_virial, resp_only=resp_only)
2473 : ELSE
2474 : CALL derivatives_four_center(qs_env, matrix_p, matrix_pza, hfx_section, para_env, &
2475 250 : 1, use_virial, resp_only=resp_only)
2476 : END IF
2477 256 : CALL dbcsr_deallocate_matrix_set(matrix_pza)
2478 : ELSE
2479 : ! -----------------------------------------
2480 : ! conventional HFX FORCES
2481 : ! -----------------------------------------
2482 234 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2483 234 : IF (x_data(1, 1)%do_hfx_ri) THEN
2484 :
2485 : CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
2486 : x_data(1, 1)%general_parameter%fraction, &
2487 : rho_ao=matrix_p, rho_ao_resp=mpa, &
2488 18 : use_virial=use_virial, resp_only=resp_only)
2489 : ELSE
2490 : CALL derivatives_four_center(qs_env, matrix_p, mpa, hfx_section, para_env, &
2491 216 : 1, use_virial, resp_only=resp_only)
2492 : END IF
2493 : END IF ! do_admm
2494 :
2495 490 : IF (use_virial) THEN
2496 884 : virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
2497 884 : virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
2498 68 : virial%pv_calculate = .FALSE.
2499 : END IF
2500 :
2501 490 : IF (debug_forces) THEN
2502 248 : fodeb(1:3) = force(1)%fock_4c(1:3, 1) - fodeb(1:3)
2503 62 : CALL para_env%sum(fodeb)
2504 62 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*hfx ", fodeb
2505 : END IF
2506 490 : IF (debug_stress .AND. use_virial) THEN
2507 0 : stdeb = -1.0_dp*fconv*virial%pv_fock_4c
2508 0 : CALL para_env%sum(stdeb)
2509 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
2510 0 : 'STRESS| Pz*hfx ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2511 : END IF
2512 : END IF ! do_hfx
2513 :
2514 : ! Stress-tensor volume contributions
2515 : ! These need to be applied at the end of qs_force
2516 1146 : IF (use_virial) THEN
2517 : ! Adding mixed Hartree energy twice, due to symmetry
2518 172 : zehartree = zehartree + 2.0_dp*ehartree
2519 172 : zexc = zexc + exc
2520 : ! ADMM contribution handled differently in qs_force
2521 172 : IF (dft_control%do_admm) THEN
2522 38 : zexc_aux_fit = zexc_aux_fit + exc_aux_fit
2523 : END IF
2524 : END IF
2525 :
2526 : ! Overlap matrix
2527 : ! H(drho+dz) + Wz
2528 : ! If ground-state density matrix solved by diagonalization, then use this
2529 1146 : IF (dft_control%qs_control%do_ls_scf) THEN
2530 : ! Ground-state density has been calculated by LS
2531 10 : eps_filter = dft_control%qs_control%eps_filter_matrix
2532 10 : CALL calculate_whz_ao_matrix(qs_env, matrix_hz, matrix_wz, eps_filter)
2533 : ELSE
2534 1136 : IF (do_ex) THEN
2535 642 : matrix_wz => p_env%w1
2536 : END IF
2537 1136 : focc = 1.0_dp
2538 1136 : IF (nspins == 1) focc = 2.0_dp
2539 1136 : CALL get_qs_env(qs_env, mos=mos)
2540 2382 : DO ispin = 1, nspins
2541 1246 : CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
2542 : CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
2543 2382 : matrix_wz(ispin)%matrix, focc, nocc)
2544 : END DO
2545 : END IF
2546 1146 : IF (nspins == 2) THEN
2547 : CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
2548 110 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2549 : END IF
2550 :
2551 1644 : IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
2552 1146 : IF (debug_stress .AND. use_virial) stdeb = virial%pv_overlap
2553 1146 : NULLIFY (scrm)
2554 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
2555 : matrix_name="OVERLAP MATRIX", &
2556 : basis_type_a="ORB", basis_type_b="ORB", &
2557 : sab_nl=sab_orb, calculate_forces=.TRUE., &
2558 1146 : matrix_p=matrix_wz(1)%matrix)
2559 :
2560 1146 : IF (SIZE(matrix_wz, 1) == 2) THEN
2561 : CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
2562 110 : alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
2563 : END IF
2564 :
2565 1146 : IF (debug_forces) THEN
2566 664 : fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
2567 166 : CALL para_env%sum(fodeb)
2568 166 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wz*dS ", fodeb
2569 : END IF
2570 1146 : IF (debug_stress .AND. use_virial) THEN
2571 0 : stdeb = fconv*(virial%pv_overlap - stdeb)
2572 0 : CALL para_env%sum(stdeb)
2573 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
2574 0 : 'STRESS| WHz ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2575 : END IF
2576 1146 : CALL dbcsr_deallocate_matrix_set(scrm)
2577 :
2578 1146 : IF (debug_forces) THEN
2579 166 : CALL total_qs_force(ftot2, force, atomic_kind_set)
2580 664 : fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
2581 166 : CALL para_env%sum(fodeb)
2582 166 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Response Force", fodeb
2583 664 : fodeb(1:3) = ftot2(1:3, 1)
2584 166 : CALL para_env%sum(fodeb)
2585 166 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Total Force ", fodeb
2586 166 : DEALLOCATE (ftot1, ftot2, ftot3)
2587 : END IF
2588 1146 : IF (debug_stress .AND. use_virial) THEN
2589 0 : stdeb = fconv*(virial%pv_virial - sttot)
2590 0 : CALL para_env%sum(stdeb)
2591 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
2592 0 : 'STRESS| Stress Response ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2593 0 : stdeb = fconv*(virial%pv_virial)
2594 0 : CALL para_env%sum(stdeb)
2595 0 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,A,T41,2(1X,ES19.11))") &
2596 0 : 'STRESS| Total Stress ', one_third_sum_diag(stdeb), det_3x3(stdeb)
2597 : IF (iounit > 0) WRITE (UNIT=iounit, FMT="(T2,3(1X,ES19.11))") &
2598 0 : stdeb(1, 1), stdeb(2, 2), stdeb(3, 3)
2599 : unitstr = "bar"
2600 : END IF
2601 :
2602 1146 : IF (do_ex) THEN
2603 642 : CALL dbcsr_deallocate_matrix_set(mpa)
2604 642 : CALL dbcsr_deallocate_matrix_set(matrix_hz)
2605 : END IF
2606 :
2607 1146 : CALL timestop(handle)
2608 :
2609 5730 : END SUBROUTINE response_force
2610 :
2611 : ! **************************************************************************************************
2612 : !> \brief ...
2613 : !> \param qs_env ...
2614 : !> \param p_env ...
2615 : !> \param matrix_hz ...
2616 : !> \param ex_env ...
2617 : !> \param debug ...
2618 : ! **************************************************************************************************
2619 26 : SUBROUTINE response_force_xtb(qs_env, p_env, matrix_hz, ex_env, debug)
2620 : TYPE(qs_environment_type), POINTER :: qs_env
2621 : TYPE(qs_p_env_type) :: p_env
2622 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz
2623 : TYPE(excited_energy_type), OPTIONAL, POINTER :: ex_env
2624 : LOGICAL, INTENT(IN), OPTIONAL :: debug
2625 :
2626 : CHARACTER(LEN=*), PARAMETER :: routineN = 'response_force_xtb'
2627 :
2628 : INTEGER :: atom_a, handle, iatom, ikind, iounit, &
2629 : is, ispin, na, natom, natorb, nimages, &
2630 : nkind, nocc, ns, nsgf, nspins
2631 : INTEGER, DIMENSION(25) :: lao
2632 : INTEGER, DIMENSION(5) :: occ
2633 : LOGICAL :: debug_forces, do_ex, use_virial
2634 : REAL(KIND=dp) :: focc
2635 26 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mcharge, mcharge1
2636 26 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, charges, charges1, ftot1, &
2637 26 : ftot2
2638 : REAL(KIND=dp), DIMENSION(3) :: fodeb
2639 26 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2640 : TYPE(cp_logger_type), POINTER :: logger
2641 26 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pz, matrix_wz, mpa, p_matrix, scrm
2642 26 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
2643 : TYPE(dbcsr_type), POINTER :: s_matrix
2644 : TYPE(dft_control_type), POINTER :: dft_control
2645 26 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2646 : TYPE(mp_para_env_type), POINTER :: para_env
2647 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2648 26 : POINTER :: sab_orb
2649 26 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2650 26 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2651 26 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2652 : TYPE(qs_ks_env_type), POINTER :: ks_env
2653 : TYPE(qs_rho_type), POINTER :: rho
2654 : TYPE(xtb_atom_type), POINTER :: xtb_kind
2655 :
2656 26 : CALL timeset(routineN, handle)
2657 :
2658 26 : IF (PRESENT(debug)) THEN
2659 26 : debug_forces = debug
2660 : ELSE
2661 0 : debug_forces = .FALSE.
2662 : END IF
2663 :
2664 26 : logger => cp_get_default_logger()
2665 26 : IF (logger%para_env%is_source()) THEN
2666 13 : iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
2667 : ELSE
2668 : iounit = -1
2669 : END IF
2670 :
2671 26 : do_ex = .FALSE.
2672 26 : IF (PRESENT(ex_env)) do_ex = .TRUE.
2673 :
2674 26 : NULLIFY (ks_env, sab_orb)
2675 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control, &
2676 26 : sab_orb=sab_orb)
2677 26 : CALL get_qs_env(qs_env=qs_env, para_env=para_env, force=force)
2678 26 : nspins = dft_control%nspins
2679 :
2680 26 : IF (debug_forces) THEN
2681 0 : CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
2682 0 : ALLOCATE (ftot1(3, natom))
2683 0 : ALLOCATE (ftot2(3, natom))
2684 0 : CALL total_qs_force(ftot1, force, atomic_kind_set)
2685 : END IF
2686 :
2687 26 : matrix_pz => p_env%p1
2688 26 : NULLIFY (mpa)
2689 26 : IF (do_ex) THEN
2690 26 : CALL dbcsr_allocate_matrix_set(mpa, nspins)
2691 62 : DO ispin = 1, nspins
2692 36 : ALLOCATE (mpa(ispin)%matrix)
2693 36 : CALL dbcsr_create(mpa(ispin)%matrix, template=matrix_pz(ispin)%matrix)
2694 36 : CALL dbcsr_copy(mpa(ispin)%matrix, matrix_pz(ispin)%matrix)
2695 36 : CALL dbcsr_add(mpa(ispin)%matrix, ex_env%matrix_pe(ispin)%matrix, 1.0_dp, 1.0_dp)
2696 62 : CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
2697 : END DO
2698 : ELSE
2699 0 : mpa => p_env%p1
2700 : END IF
2701 : !
2702 : ! START OF Tr(P+Z)Hcore
2703 : !
2704 26 : IF (nspins == 2) THEN
2705 10 : CALL dbcsr_add(mpa(1)%matrix, mpa(2)%matrix, 1.0_dp, 1.0_dp)
2706 : END IF
2707 : ! Hcore matrix
2708 26 : IF (debug_forces) fodeb(1:3) = force(1)%all_potential(1:3, 1)
2709 26 : CALL build_xtb_hab_force(qs_env, mpa(1)%matrix)
2710 26 : IF (debug_forces) THEN
2711 0 : fodeb(1:3) = force(1)%all_potential(1:3, 1) - fodeb(1:3)
2712 0 : CALL para_env%sum(fodeb)
2713 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*dHcore ", fodeb
2714 : END IF
2715 26 : IF (nspins == 2) THEN
2716 10 : CALL dbcsr_add(mpa(1)%matrix, mpa(2)%matrix, 1.0_dp, -1.0_dp)
2717 : END IF
2718 : !
2719 : ! END OF Tr(P+Z)Hcore
2720 : !
2721 26 : use_virial = .FALSE.
2722 26 : nimages = 1
2723 : !
2724 : ! Hartree potential of response density
2725 : !
2726 26 : IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
2727 : ! Mulliken charges
2728 24 : CALL get_qs_env(qs_env, rho=rho, particle_set=particle_set, matrix_s_kp=matrix_s)
2729 24 : natom = SIZE(particle_set)
2730 24 : CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
2731 120 : ALLOCATE (mcharge(natom), charges(natom, 5))
2732 72 : ALLOCATE (mcharge1(natom), charges1(natom, 5))
2733 24 : charges = 0.0_dp
2734 24 : charges1 = 0.0_dp
2735 24 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
2736 24 : nkind = SIZE(atomic_kind_set)
2737 24 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
2738 96 : ALLOCATE (aocg(nsgf, natom))
2739 24 : aocg = 0.0_dp
2740 72 : ALLOCATE (aocg1(nsgf, natom))
2741 24 : aocg1 = 0.0_dp
2742 24 : p_matrix => matrix_p(:, 1)
2743 24 : s_matrix => matrix_s(1, 1)%matrix
2744 24 : CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
2745 24 : CALL ao_charges(mpa, s_matrix, aocg1, para_env)
2746 78 : DO ikind = 1, nkind
2747 54 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
2748 54 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
2749 54 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
2750 396 : DO iatom = 1, na
2751 264 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
2752 1584 : charges(atom_a, :) = REAL(occ(:), KIND=dp)
2753 1030 : DO is = 1, natorb
2754 712 : ns = lao(is) + 1
2755 712 : charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
2756 976 : charges1(atom_a, ns) = charges1(atom_a, ns) - aocg1(is, atom_a)
2757 : END DO
2758 : END DO
2759 : END DO
2760 24 : DEALLOCATE (aocg, aocg1)
2761 288 : DO iatom = 1, natom
2762 1584 : mcharge(iatom) = SUM(charges(iatom, :))
2763 1608 : mcharge1(iatom) = SUM(charges1(iatom, :))
2764 : END DO
2765 : ! Coulomb Kernel
2766 24 : CALL xtb_coulomb_hessian(qs_env, matrix_hz, charges1, mcharge1, mcharge, mpa)
2767 : CALL calc_xtb_ehess_force(qs_env, p_matrix, mpa, charges, mcharge, charges1, &
2768 24 : mcharge1, debug_forces)
2769 : !
2770 48 : DEALLOCATE (charges, mcharge, charges1, mcharge1)
2771 : END IF
2772 : ! Overlap matrix
2773 : ! H(drho+dz) + Wz
2774 26 : matrix_wz => p_env%w1
2775 26 : focc = 0.5_dp
2776 26 : IF (nspins == 1) focc = 1.0_dp
2777 26 : CALL get_qs_env(qs_env, mos=mos)
2778 62 : DO ispin = 1, nspins
2779 36 : CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
2780 : CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
2781 62 : matrix_wz(ispin)%matrix, focc, nocc)
2782 : END DO
2783 26 : IF (nspins == 2) THEN
2784 : CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
2785 10 : alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
2786 : END IF
2787 26 : IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
2788 26 : NULLIFY (scrm)
2789 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
2790 : matrix_name="OVERLAP MATRIX", &
2791 : basis_type_a="ORB", basis_type_b="ORB", &
2792 : sab_nl=sab_orb, calculate_forces=.TRUE., &
2793 26 : matrix_p=matrix_wz(1)%matrix)
2794 26 : IF (debug_forces) THEN
2795 0 : fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
2796 0 : CALL para_env%sum(fodeb)
2797 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wz*dS ", fodeb
2798 : END IF
2799 26 : CALL dbcsr_deallocate_matrix_set(scrm)
2800 :
2801 26 : IF (debug_forces) THEN
2802 0 : CALL total_qs_force(ftot2, force, atomic_kind_set)
2803 0 : fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
2804 0 : CALL para_env%sum(fodeb)
2805 0 : IF (iounit > 0) WRITE (iounit, "(T3,A,T30,3F16.8)") "DEBUG:: Response Force", fodeb
2806 0 : DEALLOCATE (ftot1, ftot2)
2807 : END IF
2808 :
2809 26 : IF (do_ex) THEN
2810 26 : CALL dbcsr_deallocate_matrix_set(mpa)
2811 : END IF
2812 :
2813 26 : CALL timestop(handle)
2814 :
2815 52 : END SUBROUTINE response_force_xtb
2816 :
2817 : ! **************************************************************************************************
2818 : !> \brief Win = focc*(P*(H[P_out - P_in] + H[Z] )*P)
2819 : !> Langrange multiplier matrix with response and perturbation (Harris) kernel matrices
2820 : !>
2821 : !> \param qs_env ...
2822 : !> \param matrix_hz ...
2823 : !> \param matrix_whz ...
2824 : !> \param eps_filter ...
2825 : !> \param
2826 : !> \par History
2827 : !> 2020.2 created [Fabian Belleflamme]
2828 : !> \author Fabian Belleflamme
2829 : ! **************************************************************************************************
2830 10 : SUBROUTINE calculate_whz_ao_matrix(qs_env, matrix_hz, matrix_whz, eps_filter)
2831 :
2832 : TYPE(qs_environment_type), POINTER :: qs_env
2833 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
2834 : POINTER :: matrix_hz
2835 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
2836 : POINTER :: matrix_whz
2837 : REAL(KIND=dp), INTENT(IN) :: eps_filter
2838 :
2839 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_whz_ao_matrix'
2840 :
2841 : INTEGER :: handle, ispin, nspins
2842 : REAL(KIND=dp) :: scaling
2843 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
2844 : TYPE(dbcsr_type) :: matrix_tmp
2845 : TYPE(dft_control_type), POINTER :: dft_control
2846 : TYPE(mp_para_env_type), POINTER :: para_env
2847 : TYPE(qs_rho_type), POINTER :: rho
2848 :
2849 10 : CALL timeset(routineN, handle)
2850 :
2851 10 : CPASSERT(ASSOCIATED(qs_env))
2852 10 : CPASSERT(ASSOCIATED(matrix_hz))
2853 10 : CPASSERT(ASSOCIATED(matrix_whz))
2854 :
2855 : CALL get_qs_env(qs_env=qs_env, &
2856 : dft_control=dft_control, &
2857 : rho=rho, &
2858 10 : para_env=para_env)
2859 10 : nspins = dft_control%nspins
2860 10 : CALL qs_rho_get(rho, rho_ao=rho_ao)
2861 :
2862 : ! init temp matrix
2863 : CALL dbcsr_create(matrix_tmp, template=matrix_hz(1)%matrix, &
2864 10 : matrix_type=dbcsr_type_no_symmetry)
2865 :
2866 : !Spin factors simplify to
2867 10 : scaling = 1.0_dp
2868 10 : IF (nspins == 1) scaling = 0.5_dp
2869 :
2870 : ! Operation in MO-solver :
2871 : ! Whz = focc*(CC^T*Hz*CC^T)
2872 : ! focc = 2.0_dp Closed-shell
2873 : ! focc = 1.0_dp Open-shell
2874 :
2875 : ! Operation in AO-solver :
2876 : ! Whz = (scaling*P)*(focc*Hz)*(scaling*P)
2877 : ! focc see above
2878 : ! scaling = 0.5_dp Closed-shell (P = 2*CC^T), WHz = (0.5*P)*(2*Hz)*(0.5*P)
2879 : ! scaling = 1.0_dp Open-shell, WHz = P*Hz*P
2880 :
2881 : ! Spin factors from Hz and P simplify to
2882 : scaling = 1.0_dp
2883 10 : IF (nspins == 1) scaling = 0.5_dp
2884 :
2885 20 : DO ispin = 1, nspins
2886 :
2887 : ! tmp = H*CC^T
2888 : CALL dbcsr_multiply("N", "N", scaling, matrix_hz(ispin)%matrix, rho_ao(ispin)%matrix, &
2889 10 : 0.0_dp, matrix_tmp, filter_eps=eps_filter)
2890 : ! WHz = CC^T*tmp
2891 : ! WHz = Wz + (scaling*P)*(focc*Hz)*(scaling*P)
2892 : ! WHz = Wz + scaling*(P*Hz*P)
2893 : CALL dbcsr_multiply("N", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_tmp, &
2894 : 1.0_dp, matrix_whz(ispin)%matrix, filter_eps=eps_filter, &
2895 20 : retain_sparsity=.TRUE.)
2896 :
2897 : END DO
2898 :
2899 10 : CALL dbcsr_release(matrix_tmp)
2900 :
2901 10 : CALL timestop(handle)
2902 :
2903 10 : END SUBROUTINE calculate_whz_ao_matrix
2904 :
2905 : ! **************************************************************************************************
2906 :
2907 : END MODULE response_solver
|
