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