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