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