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