(git:e7e05ae)
qs_tddfpt2_forces.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
9  USE admm_types, ONLY: admm_type,&
11  USE atomic_kind_types, ONLY: atomic_kind_type,&
14  USE cp_control_types, ONLY: dft_control_type,&
15  tddfpt2_control_type
25  cp_fm_struct_type
26  USE cp_fm_types, ONLY: cp_fm_copy_general,&
27  cp_fm_create,&
29  cp_fm_release,&
31  cp_fm_type
34  cp_logger_type
35  USE dbcsr_api, ONLY: &
36  dbcsr_add, dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_dot, dbcsr_p_type, &
37  dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric
38  USE exstates_types, ONLY: excited_energy_type,&
44  hartree_local_type
46  USE hfx_ri, ONLY: hfx_ri_update_ks
47  USE hfx_types, ONLY: hfx_type
49  oe_shift,&
56  section_vals_type,&
58  USE kinds, ONLY: default_string_length,&
59  dp
60  USE message_passing, ONLY: mp_para_env_type
61  USE mulliken, ONLY: ao_charges
62  USE parallel_gemm_api, ONLY: parallel_gemm
63  USE particle_types, ONLY: particle_type
64  USE pw_env_types, ONLY: pw_env_get,&
65  pw_env_type
66  USE pw_methods, ONLY: pw_axpy,&
67  pw_scale,&
68  pw_transfer,&
69  pw_zero
70  USE pw_poisson_methods, ONLY: pw_poisson_solve
71  USE pw_poisson_types, ONLY: pw_poisson_type
72  USE pw_pool_types, ONLY: pw_pool_type
73  USE pw_types, ONLY: pw_c1d_gs_type,&
74  pw_r3d_rs_type
78  USE qs_environment_types, ONLY: get_qs_env,&
79  qs_environment_type,&
83  qs_force_type,&
84  sum_qs_force,&
87  USE qs_fxc, ONLY: qs_fxc_analytic,&
90  USE qs_integrate_potential, ONLY: integrate_v_rspace
91  USE qs_kernel_types, ONLY: kernel_env_type
92  USE qs_kind_types, ONLY: get_qs_kind,&
94  qs_kind_type
95  USE qs_ks_atom, ONLY: update_ks_atom
98  USE qs_ks_types, ONLY: qs_ks_env_type
101  local_rho_type
102  USE qs_mo_types, ONLY: get_mo_set,&
103  mo_set_type
104  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
105  USE qs_oce_types, ONLY: oce_matrix_type
109  USE qs_rho0_methods, ONLY: init_rho0
110  USE qs_rho0_types, ONLY: get_rho0_mpole
113  USE qs_rho_atom_types, ONLY: rho_atom_type
114  USE qs_rho_types, ONLY: qs_rho_create,&
115  qs_rho_get,&
116  qs_rho_set,&
117  qs_rho_type
119  stda_force
120  USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
121  USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos,&
122  tddfpt_work_matrices
124  USE task_list_types, ONLY: task_list_type
125  USE xtb_ehess, ONLY: xtb_coulomb_hessian
126  USE xtb_types, ONLY: get_xtb_atom_param,&
127  xtb_atom_type
128 #include "./base/base_uses.f90"
129 
130  IMPLICIT NONE
131 
132  PRIVATE
133 
134  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_forces'
135 
136  PUBLIC :: tddfpt_forces_main
137 
138 ! **************************************************************************************************
139 
140 CONTAINS
141 
142 ! **************************************************************************************************
143 !> \brief Perform TDDFPT gradient calculation.
144 !> \param qs_env Quickstep environment
145 !> \param gs_mos ...
146 !> \param ex_env ...
147 !> \param kernel_env ...
148 !> \param sub_env ...
149 !> \param work_matrices ...
150 !> \par History
151 !> * 10.2022 created JHU
152 ! **************************************************************************************************
153  SUBROUTINE tddfpt_forces_main(qs_env, gs_mos, ex_env, kernel_env, sub_env, work_matrices)
154  TYPE(qs_environment_type), POINTER :: qs_env
155  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
156  POINTER :: gs_mos
157  TYPE(excited_energy_type), POINTER :: ex_env
158  TYPE(kernel_env_type) :: kernel_env
159  TYPE(tddfpt_subgroup_env_type) :: sub_env
160  TYPE(tddfpt_work_matrices) :: work_matrices
161 
162  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_forces_main'
163 
164  INTEGER :: handle, ispin, nspins
165  TYPE(admm_type), POINTER :: admm_env
166  TYPE(cp_fm_struct_type), POINTER :: matrix_struct
167  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pe_asymm, matrix_pe_symm, &
168  matrix_s, matrix_s_aux_fit
169  TYPE(dft_control_type), POINTER :: dft_control
170  TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
171 
172  CALL timeset(routinen, handle)
173 
174  CALL get_qs_env(qs_env, dft_control=dft_control)
175  nspins = dft_control%nspins
176  tddfpt_control => dft_control%tddfpt2_control
177  ! rhs of linres equation
178  IF (ASSOCIATED(ex_env%cpmos)) THEN
179  DO ispin = 1, SIZE(ex_env%cpmos)
180  CALL cp_fm_release(ex_env%cpmos(ispin))
181  END DO
182  DEALLOCATE (ex_env%cpmos)
183  END IF
184  ALLOCATE (ex_env%cpmos(nspins))
185  DO ispin = 1, nspins
186  CALL cp_fm_get_info(matrix=ex_env%evect(ispin), matrix_struct=matrix_struct)
187  CALL cp_fm_create(ex_env%cpmos(ispin), matrix_struct)
188  CALL cp_fm_set_all(ex_env%cpmos(ispin), 0.0_dp)
189  END DO
190  CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
191  NULLIFY (matrix_pe_asymm, matrix_pe_symm)
192  CALL dbcsr_allocate_matrix_set(ex_env%matrix_pe, nspins)
193  CALL dbcsr_allocate_matrix_set(matrix_pe_symm, nspins)
194  CALL dbcsr_allocate_matrix_set(matrix_pe_asymm, nspins)
195  DO ispin = 1, nspins
196  ALLOCATE (ex_env%matrix_pe(ispin)%matrix)
197  CALL dbcsr_create(ex_env%matrix_pe(ispin)%matrix, template=matrix_s(1)%matrix)
198  CALL dbcsr_copy(ex_env%matrix_pe(ispin)%matrix, matrix_s(1)%matrix)
199  CALL dbcsr_set(ex_env%matrix_pe(ispin)%matrix, 0.0_dp)
200 
201  ALLOCATE (matrix_pe_symm(ispin)%matrix)
202  CALL dbcsr_create(matrix_pe_symm(ispin)%matrix, template=matrix_s(1)%matrix)
203  CALL dbcsr_copy(matrix_pe_symm(ispin)%matrix, ex_env%matrix_pe(ispin)%matrix)
204 
205  ALLOCATE (matrix_pe_asymm(ispin)%matrix)
206  CALL dbcsr_create(matrix_pe_asymm(ispin)%matrix, template=matrix_s(1)%matrix, &
207  matrix_type=dbcsr_type_antisymmetric)
208  CALL dbcsr_complete_redistribute(ex_env%matrix_pe(ispin)%matrix, matrix_pe_asymm(ispin)%matrix)
209 
210  CALL tddfpt_resvec1(ex_env%evect(ispin), gs_mos(ispin)%mos_occ, &
211  matrix_s(1)%matrix, ex_env%matrix_pe(ispin)%matrix)
212  END DO
213  !
214  ! ground state ADMM!
215  IF (dft_control%do_admm) THEN
216  CALL get_qs_env(qs_env, admm_env=admm_env)
217  CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
218  CALL dbcsr_allocate_matrix_set(ex_env%matrix_pe_admm, nspins)
219  DO ispin = 1, nspins
220  ALLOCATE (ex_env%matrix_pe_admm(ispin)%matrix)
221  CALL dbcsr_create(ex_env%matrix_pe_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
222  CALL dbcsr_copy(ex_env%matrix_pe_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
223  CALL dbcsr_set(ex_env%matrix_pe_admm(ispin)%matrix, 0.0_dp)
224  CALL tddfpt_resvec1_admm(ex_env%matrix_pe(ispin)%matrix, &
225  admm_env, ex_env%matrix_pe_admm(ispin)%matrix)
226  END DO
227  END IF
228  !
229  CALL dbcsr_allocate_matrix_set(ex_env%matrix_hz, nspins)
230  DO ispin = 1, nspins
231  ALLOCATE (ex_env%matrix_hz(ispin)%matrix)
232  CALL dbcsr_create(ex_env%matrix_hz(ispin)%matrix, template=matrix_s(1)%matrix)
233  CALL dbcsr_copy(ex_env%matrix_hz(ispin)%matrix, matrix_s(1)%matrix)
234  CALL dbcsr_set(ex_env%matrix_hz(ispin)%matrix, 0.0_dp)
235  END DO
236  IF (dft_control%qs_control%xtb) THEN
237  CALL tddfpt_resvec2_xtb(qs_env, ex_env%matrix_pe, gs_mos, ex_env%matrix_hz, ex_env%cpmos)
238  ELSE
239  CALL tddfpt_resvec2(qs_env, ex_env%matrix_pe, ex_env%matrix_pe_admm, &
240  gs_mos, ex_env%matrix_hz, ex_env%cpmos)
241  END IF
242  !
243  CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1, nspins)
244  CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_asymm, nspins)
245  DO ispin = 1, nspins
246  ALLOCATE (ex_env%matrix_px1(ispin)%matrix)
247  CALL dbcsr_create(ex_env%matrix_px1(ispin)%matrix, template=matrix_s(1)%matrix)
248  CALL dbcsr_copy(ex_env%matrix_px1(ispin)%matrix, matrix_s(1)%matrix)
249  CALL dbcsr_set(ex_env%matrix_px1(ispin)%matrix, 0.0_dp)
250 
251  ALLOCATE (ex_env%matrix_px1_asymm(ispin)%matrix)
252  CALL dbcsr_create(ex_env%matrix_px1_asymm(ispin)%matrix, template=matrix_s(1)%matrix, &
253  matrix_type=dbcsr_type_antisymmetric)
254  CALL dbcsr_complete_redistribute(ex_env%matrix_px1(ispin)%matrix, ex_env%matrix_px1_asymm(ispin)%matrix)
255  END DO
256  ! Kernel ADMM
257  IF (tddfpt_control%do_admm) THEN
258  CALL get_qs_env(qs_env, admm_env=admm_env)
259  CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit)
260  CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_admm, nspins)
261  CALL dbcsr_allocate_matrix_set(ex_env%matrix_px1_admm_asymm, nspins)
262  DO ispin = 1, nspins
263  ALLOCATE (ex_env%matrix_px1_admm(ispin)%matrix)
264  CALL dbcsr_create(ex_env%matrix_px1_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
265  CALL dbcsr_copy(ex_env%matrix_px1_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
266  CALL dbcsr_set(ex_env%matrix_px1_admm(ispin)%matrix, 0.0_dp)
267 
268  ALLOCATE (ex_env%matrix_px1_admm_asymm(ispin)%matrix)
269  CALL dbcsr_create(ex_env%matrix_px1_admm_asymm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix, &
270  matrix_type=dbcsr_type_antisymmetric)
271  CALL dbcsr_complete_redistribute(ex_env%matrix_px1_admm(ispin)%matrix, &
272  ex_env%matrix_px1_admm_asymm(ispin)%matrix)
273  END DO
274  END IF
275  ! TDA forces
276  CALL tddfpt_forces(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices)
277  ! Rotate res vector cpmos into original frame of occupied orbitals
278  CALL tddfpt_resvec3(qs_env, ex_env%cpmos, work_matrices)
279 
280  CALL dbcsr_deallocate_matrix_set(matrix_pe_symm)
281  CALL dbcsr_deallocate_matrix_set(matrix_pe_asymm)
282 
283  CALL timestop(handle)
284 
285  END SUBROUTINE tddfpt_forces_main
286 
287 ! **************************************************************************************************
288 !> \brief Calculate direct tddft forces
289 !> \param qs_env ...
290 !> \param ex_env ...
291 !> \param gs_mos ...
292 !> \param kernel_env ...
293 !> \param sub_env ...
294 !> \param work_matrices ...
295 !> \par History
296 !> * 01.2020 screated [JGH]
297 ! **************************************************************************************************
298  SUBROUTINE tddfpt_forces(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices)
299 
300  TYPE(qs_environment_type), POINTER :: qs_env
301  TYPE(excited_energy_type), POINTER :: ex_env
302  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
303  POINTER :: gs_mos
304  TYPE(kernel_env_type), INTENT(IN) :: kernel_env
305  TYPE(tddfpt_subgroup_env_type) :: sub_env
306  TYPE(tddfpt_work_matrices) :: work_matrices
307 
308  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_forces'
309 
310  INTEGER :: handle
311  INTEGER, ALLOCATABLE, DIMENSION(:) :: natom_of_kind
312  LOGICAL :: debug_forces
313  REAL(kind=dp) :: ehartree, exc
314  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
315  TYPE(dft_control_type), POINTER :: dft_control
316  TYPE(qs_force_type), DIMENSION(:), POINTER :: ks_force, td_force
317 
318  CALL timeset(routinen, handle)
319 
320  ! for extended debug output
321  debug_forces = ex_env%debug_forces
322  ! prepare force array
323  CALL get_qs_env(qs_env, dft_control=dft_control, force=ks_force, &
324  atomic_kind_set=atomic_kind_set)
325  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
326  NULLIFY (td_force)
327  CALL allocate_qs_force(td_force, natom_of_kind)
328  DEALLOCATE (natom_of_kind)
329  CALL zero_qs_force(td_force)
330  CALL set_qs_env(qs_env, force=td_force)
331  !
332  IF (dft_control%qs_control%xtb) THEN
333  CALL tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, &
334  work_matrices, debug_forces)
335  ELSE
336  !
337  CALL exstate_potential_release(ex_env)
338  CALL ks_ref_potential(qs_env, ex_env%vh_rspace, ex_env%vxc_rspace, &
339  ex_env%vtau_rspace, ex_env%vadmm_rspace, ehartree, exc)
340  CALL ks_ref_potential_atom(qs_env, ex_env%local_rho_set, ex_env%local_rho_set_admm, &
341  ex_env%vh_rspace)
342  CALL tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, &
343  work_matrices, debug_forces)
344  END IF
345  !
346  ! add TD and KS forces
347  CALL get_qs_env(qs_env, force=td_force)
348  CALL sum_qs_force(ks_force, td_force)
349  CALL set_qs_env(qs_env, force=ks_force)
350  CALL deallocate_qs_force(td_force)
351  !
352  CALL timestop(handle)
353 
354  END SUBROUTINE tddfpt_forces
355 
356 ! **************************************************************************************************
357 !> \brief Calculate direct tddft forces
358 !> \param qs_env ...
359 !> \param ex_env ...
360 !> \param gs_mos ...
361 !> \param kernel_env ...
362 !> \param sub_env ...
363 !> \param work_matrices ...
364 !> \param debug_forces ...
365 !> \par History
366 !> * 01.2020 screated [JGH]
367 ! **************************************************************************************************
368  SUBROUTINE tddfpt_force_direct(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, &
369  debug_forces)
370 
371  TYPE(qs_environment_type), POINTER :: qs_env
372  TYPE(excited_energy_type), POINTER :: ex_env
373  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
374  POINTER :: gs_mos
375  TYPE(kernel_env_type), INTENT(IN) :: kernel_env
376  TYPE(tddfpt_subgroup_env_type) :: sub_env
377  TYPE(tddfpt_work_matrices) :: work_matrices
378  LOGICAL :: debug_forces
379 
380  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_force_direct'
381 
382  INTEGER :: handle, iounit, ispin, natom, norb, &
383  nspins
384  REAL(kind=dp) :: evalue
385  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ftot1, ftot2
386  REAL(kind=dp), DIMENSION(3) :: fodeb
387  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
388  TYPE(cp_fm_type), DIMENSION(:), POINTER :: evect
389  TYPE(cp_logger_type), POINTER :: logger
390  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_wx1, &
391  matrix_wz, scrm
392  TYPE(dft_control_type), POINTER :: dft_control
393  TYPE(mp_para_env_type), POINTER :: para_env
394  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
395  POINTER :: sab_orb
396  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
397  TYPE(qs_ks_env_type), POINTER :: ks_env
398  TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
399 
400  CALL timeset(routinen, handle)
401 
402  logger => cp_get_default_logger()
403  IF (logger%para_env%is_source()) THEN
404  iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
405  ELSE
406  iounit = -1
407  END IF
408 
409  evect => ex_env%evect
410 
411  CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, para_env=para_env, &
412  sab_orb=sab_orb, dft_control=dft_control, force=force)
413  NULLIFY (tddfpt_control)
414  tddfpt_control => dft_control%tddfpt2_control
415  nspins = dft_control%nspins
416 
417  IF (debug_forces) THEN
418  CALL get_qs_env(qs_env, natom=natom, atomic_kind_set=atomic_kind_set)
419  ALLOCATE (ftot1(3, natom))
420  CALL total_qs_force(ftot1, force, atomic_kind_set)
421  END IF
422 
423  CALL tddfpt_kernel_force(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, debug_forces)
424 
425  ! Overlap matrix
426  matrix_wx1 => ex_env%matrix_wx1
427  CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, matrix_ks=matrix_ks)
428  NULLIFY (matrix_wz)
429  CALL dbcsr_allocate_matrix_set(matrix_wz, nspins)
430  DO ispin = 1, nspins
431  ALLOCATE (matrix_wz(ispin)%matrix)
432  CALL dbcsr_create(matrix=matrix_wz(ispin)%matrix, template=matrix_s(1)%matrix)
433  CALL cp_dbcsr_alloc_block_from_nbl(matrix_wz(ispin)%matrix, sab_orb)
434  CALL dbcsr_set(matrix_wz(ispin)%matrix, 0.0_dp)
435  CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=norb)
436  CALL cp_dbcsr_plus_fm_fm_t(matrix_wz(ispin)%matrix, matrix_v=evect(ispin), ncol=norb)
437  evalue = ex_env%evalue
438  IF (tddfpt_control%oe_corr == oe_shift) THEN
439  evalue = ex_env%evalue - tddfpt_control%ev_shift
440  END IF
441  CALL dbcsr_scale(matrix_wz(ispin)%matrix, evalue)
442  CALL calculate_wx_matrix(gs_mos(ispin)%mos_occ, evect(ispin), matrix_ks(ispin)%matrix, &
443  matrix_wz(ispin)%matrix)
444  END DO
445  IF (nspins == 2) THEN
446  CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
447  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
448  END IF
449  NULLIFY (scrm)
450  IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
451  CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
452  matrix_name="OVERLAP MATRIX", &
453  basis_type_a="ORB", basis_type_b="ORB", &
454  sab_nl=sab_orb, calculate_forces=.true., &
455  matrix_p=matrix_wz(1)%matrix)
456  CALL dbcsr_deallocate_matrix_set(scrm)
457  CALL dbcsr_deallocate_matrix_set(matrix_wz)
458  IF (debug_forces) THEN
459  fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
460  CALL para_env%sum(fodeb)
461  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Wx*dS ", fodeb
462  END IF
463 
464  ! Overlap matrix
465  CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, matrix_ks=matrix_ks)
466  NULLIFY (matrix_wz)
467  CALL dbcsr_allocate_matrix_set(matrix_wz, nspins)
468  DO ispin = 1, nspins
469  ALLOCATE (matrix_wz(ispin)%matrix)
470  CALL dbcsr_create(matrix=matrix_wz(ispin)%matrix, template=matrix_s(1)%matrix)
471  CALL cp_dbcsr_alloc_block_from_nbl(matrix_wz(ispin)%matrix, sab_orb)
472  CALL dbcsr_set(matrix_wz(ispin)%matrix, 0.0_dp)
473  CALL cp_fm_get_info(gs_mos(ispin)%mos_occ, ncol_global=norb)
474  evalue = ex_env%evalue
475  IF (tddfpt_control%oe_corr == oe_shift) THEN
476  evalue = ex_env%evalue - tddfpt_control%ev_shift
477  END IF
478  CALL calculate_xwx_matrix(gs_mos(ispin)%mos_occ, evect(ispin), matrix_s(1)%matrix, &
479  matrix_ks(ispin)%matrix, matrix_wz(ispin)%matrix, evalue)
480  END DO
481  IF (nspins == 2) THEN
482  CALL dbcsr_add(matrix_wz(1)%matrix, matrix_wz(2)%matrix, &
483  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
484  END IF
485  NULLIFY (scrm)
486  IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
487  CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
488  matrix_name="OVERLAP MATRIX", &
489  basis_type_a="ORB", basis_type_b="ORB", &
490  sab_nl=sab_orb, calculate_forces=.true., &
491  matrix_p=matrix_wz(1)%matrix)
492  CALL dbcsr_deallocate_matrix_set(scrm)
493  CALL dbcsr_deallocate_matrix_set(matrix_wz)
494  IF (debug_forces) THEN
495  fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
496  CALL para_env%sum(fodeb)
497  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: xWx*dS ", fodeb
498  END IF
499 
500  ! Overlap matrix
501  IF (ASSOCIATED(matrix_wx1)) THEN
502  IF (nspins == 2) THEN
503  CALL dbcsr_add(matrix_wx1(1)%matrix, matrix_wx1(2)%matrix, &
504  alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
505  END IF
506  NULLIFY (scrm)
507  IF (debug_forces) fodeb(1:3) = force(1)%overlap(1:3, 1)
508  CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
509  matrix_name="OVERLAP MATRIX", &
510  basis_type_a="ORB", basis_type_b="ORB", &
511  sab_nl=sab_orb, calculate_forces=.true., &
512  matrix_p=matrix_wx1(1)%matrix)
513  CALL dbcsr_deallocate_matrix_set(scrm)
514  IF (debug_forces) THEN
515  fodeb(1:3) = force(1)%overlap(1:3, 1) - fodeb(1:3)
516  CALL para_env%sum(fodeb)
517  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: WK*dS ", fodeb
518  END IF
519  END IF
520 
521  IF (debug_forces) THEN
522  ALLOCATE (ftot2(3, natom))
523  CALL total_qs_force(ftot2, force, atomic_kind_set)
524  fodeb(1:3) = ftot2(1:3, 1) - ftot1(1:3, 1)
525  CALL para_env%sum(fodeb)
526  IF (iounit > 0) WRITE (iounit, "(T3,A,T30,3F16.8)") "DEBUG:: Excitation Force", fodeb
527  DEALLOCATE (ftot1, ftot2)
528  END IF
529 
530  CALL timestop(handle)
531 
532  END SUBROUTINE tddfpt_force_direct
533 
534 ! **************************************************************************************************
535 !> \brief ...
536 !> \param evect ...
537 !> \param mos_occ ...
538 !> \param matrix_s ...
539 !> \param matrix_pe ...
540 ! **************************************************************************************************
541  SUBROUTINE tddfpt_resvec1(evect, mos_occ, matrix_s, matrix_pe)
542 
543  TYPE(cp_fm_type), INTENT(IN) :: evect, mos_occ
544  TYPE(dbcsr_type), POINTER :: matrix_s, matrix_pe
545 
546  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_resvec1'
547 
548  INTEGER :: handle, iounit, nao, norb
549  REAL(kind=dp) :: tmp
550  TYPE(cp_fm_struct_type), POINTER :: fmstruct, fmstruct2
551  TYPE(cp_fm_type) :: cxmat, xxmat
552  TYPE(cp_logger_type), POINTER :: logger
553 
554  CALL timeset(routinen, handle)
555  ! X*X^T
556  CALL cp_fm_get_info(mos_occ, nrow_global=nao, ncol_global=norb)
557  CALL cp_dbcsr_plus_fm_fm_t(matrix_pe, matrix_v=evect, ncol=norb)
558  ! X^T*S*X
559  CALL cp_fm_get_info(evect, matrix_struct=fmstruct)
560  NULLIFY (fmstruct2)
561  CALL cp_fm_struct_create(fmstruct=fmstruct2, template_fmstruct=fmstruct, &
562  nrow_global=norb, ncol_global=norb)
563  CALL cp_fm_create(xxmat, matrix_struct=fmstruct2)
564  CALL cp_fm_struct_release(fmstruct2)
565  CALL cp_fm_create(cxmat, matrix_struct=fmstruct)
566  CALL cp_dbcsr_sm_fm_multiply(matrix_s, evect, cxmat, norb, alpha=1.0_dp, beta=0.0_dp)
567  CALL parallel_gemm('T', 'N', norb, norb, nao, 1.0_dp, cxmat, evect, 0.0_dp, xxmat)
568  CALL parallel_gemm('N', 'N', nao, norb, norb, 1.0_dp, mos_occ, xxmat, 0.0_dp, cxmat)
569  CALL cp_fm_release(xxmat)
570  ! C*C^T*XX
571  CALL cp_dbcsr_plus_fm_fm_t(matrix_pe, matrix_v=mos_occ, matrix_g=cxmat, &
572  ncol=norb, alpha=-1.0_dp, symmetry_mode=1)
573  CALL cp_fm_release(cxmat)
574  !
575  ! Test for Tr(Pe*S)=0
576  CALL dbcsr_dot(matrix_pe, matrix_s, tmp)
577  IF (abs(tmp) > 1.e-08_dp) THEN
578  logger => cp_get_default_logger()
579  IF (logger%para_env%is_source()) THEN
580  iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
581  ELSE
582  iounit = -1
583  END IF
584  cpwarn("Electron count of excitation density matrix is non-zero.")
585  IF (iounit > 0) THEN
586  WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", tmp
587  WRITE (iounit, "(T2,A,/)") repeat("*", 79)
588  END IF
589  END IF
590  !
591 
592  CALL timestop(handle)
593 
594  END SUBROUTINE tddfpt_resvec1
595 
596 ! **************************************************************************************************
597 !> \brief PA = A * P * A(T)
598 !> \param matrix_pe ...
599 !> \param admm_env ...
600 !> \param matrix_pe_admm ...
601 ! **************************************************************************************************
602  SUBROUTINE tddfpt_resvec1_admm(matrix_pe, admm_env, matrix_pe_admm)
603 
604  TYPE(dbcsr_type), POINTER :: matrix_pe
605  TYPE(admm_type), POINTER :: admm_env
606  TYPE(dbcsr_type), POINTER :: matrix_pe_admm
607 
608  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_resvec1_admm'
609 
610  INTEGER :: handle, nao, nao_aux
611 
612  CALL timeset(routinen, handle)
613  !
614  nao_aux = admm_env%nao_aux_fit
615  nao = admm_env%nao_orb
616  !
617  CALL copy_dbcsr_to_fm(matrix_pe, admm_env%work_orb_orb)
618  CALL parallel_gemm('N', 'N', nao_aux, nao, nao, &
619  1.0_dp, admm_env%A, admm_env%work_orb_orb, 0.0_dp, &
620  admm_env%work_aux_orb)
621  CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, &
622  1.0_dp, admm_env%work_aux_orb, admm_env%A, 0.0_dp, &
623  admm_env%work_aux_aux)
624  CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, matrix_pe_admm, keep_sparsity=.true.)
625  !
626  CALL timestop(handle)
627 
628  END SUBROUTINE tddfpt_resvec1_admm
629 
630 ! **************************************************************************************************
631 !> \brief ...
632 !> \param qs_env ...
633 !> \param matrix_pe ...
634 !> \param matrix_pe_admm ...
635 !> \param gs_mos ...
636 !> \param matrix_hz ...
637 !> \param cpmos ...
638 ! **************************************************************************************************
639  SUBROUTINE tddfpt_resvec2(qs_env, matrix_pe, matrix_pe_admm, gs_mos, matrix_hz, cpmos)
640 
641  TYPE(qs_environment_type), POINTER :: qs_env
642  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pe, matrix_pe_admm
643  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
644  POINTER :: gs_mos
645  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz
646  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: cpmos
647 
648  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_resvec2'
649 
650  CHARACTER(LEN=default_string_length) :: basis_type
651  INTEGER :: handle, iounit, ispin, mspin, n_rep_hf, &
652  nao, nao_aux, natom, norb, nspins
653  LOGICAL :: deriv2_analytic, distribute_fock_matrix, &
654  do_hfx, gapw, gapw_xc, &
655  hfx_treat_lsd_in_core, &
656  s_mstruct_changed
657  REAL(kind=dp) :: eh1, focc, rhotot, thartree
658  REAL(kind=dp), DIMENSION(2) :: total_rho
659  REAL(kind=dp), DIMENSION(:), POINTER :: qlm_tot
660  TYPE(admm_type), POINTER :: admm_env
661  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
662  TYPE(cp_fm_type), POINTER :: mos
663  TYPE(cp_logger_type), POINTER :: logger
664  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: msaux
665  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mhz, mpe
666  TYPE(dbcsr_type), POINTER :: dbwork
667  TYPE(dft_control_type), POINTER :: dft_control
668  TYPE(hartree_local_type), POINTER :: hartree_local
669  TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
670  TYPE(local_rho_type), POINTER :: local_rho_set, local_rho_set_admm
671  TYPE(mp_para_env_type), POINTER :: para_env
672  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
673  POINTER :: sab, sab_aux_fit
674  TYPE(oce_matrix_type), POINTER :: oce
675  TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
676  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_g_aux, rhoz_g_aux, trho_g, &
677  trho_xc_g
678  TYPE(pw_env_type), POINTER :: pw_env
679  TYPE(pw_poisson_type), POINTER :: poisson_env
680  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
681  TYPE(pw_r3d_rs_type) :: v_hartree_rspace
682  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_r_aux, rhoz_r_aux, tau_r, &
683  trho_r, trho_xc_r, v_xc, v_xc_tau
684  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
685  TYPE(qs_ks_env_type), POINTER :: ks_env
686  TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit, rho_xc, rhoz_aux, trho
687  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho1_atom_set, rho_atom_set
688  TYPE(section_vals_type), POINTER :: hfx_section, input, xc_section
689  TYPE(task_list_type), POINTER :: task_list
690 
691  CALL timeset(routinen, handle)
692 
693  NULLIFY (pw_env)
694  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, ks_env=ks_env, &
695  dft_control=dft_control, para_env=para_env)
696  cpassert(ASSOCIATED(pw_env))
697  nspins = dft_control%nspins
698  gapw = dft_control%qs_control%gapw
699  gapw_xc = dft_control%qs_control%gapw_xc
700 
701  cpassert(.NOT. dft_control%tddfpt2_control%do_exck)
702  cpassert(.NOT. dft_control%tddfpt2_control%do_hfxsr)
703  cpassert(.NOT. dft_control%tddfpt2_control%do_hfxlr)
704 
705  NULLIFY (auxbas_pw_pool, poisson_env)
706  ! gets the tmp grids
707  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
708  poisson_env=poisson_env)
709 
710  CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
711  CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
712  CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
713 
714  ALLOCATE (trho_r(nspins), trho_g(nspins))
715  DO ispin = 1, nspins
716  CALL auxbas_pw_pool%create_pw(trho_r(ispin))
717  CALL auxbas_pw_pool%create_pw(trho_g(ispin))
718  END DO
719  IF (gapw_xc) THEN
720  ALLOCATE (trho_xc_r(nspins), trho_xc_g(nspins))
721  DO ispin = 1, nspins
722  CALL auxbas_pw_pool%create_pw(trho_xc_r(ispin))
723  CALL auxbas_pw_pool%create_pw(trho_xc_g(ispin))
724  END DO
725  END IF
726 
727  ! GAPW/GAPW_XC initializations
728  NULLIFY (hartree_local, local_rho_set)
729  IF (gapw) THEN
730  CALL get_qs_env(qs_env, &
731  atomic_kind_set=atomic_kind_set, &
732  natom=natom, &
733  qs_kind_set=qs_kind_set)
734  CALL local_rho_set_create(local_rho_set)
735  CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
736  qs_kind_set, dft_control, para_env)
737  CALL init_rho0(local_rho_set, qs_env, dft_control%qs_control%gapw_control, &
738  zcore=0.0_dp)
739  CALL rho0_s_grid_create(pw_env, local_rho_set%rho0_mpole)
740  CALL hartree_local_create(hartree_local)
741  CALL init_coulomb_local(hartree_local, natom)
742  ELSEIF (gapw_xc) THEN
743  CALL get_qs_env(qs_env, &
744  atomic_kind_set=atomic_kind_set, &
745  qs_kind_set=qs_kind_set)
746  CALL local_rho_set_create(local_rho_set)
747  CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
748  qs_kind_set, dft_control, para_env)
749  END IF
750 
751  total_rho = 0.0_dp
752  CALL pw_zero(rho_tot_gspace)
753  DO ispin = 1, nspins
754  CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_pe(ispin)%matrix, &
755  rho=trho_r(ispin), &
756  rho_gspace=trho_g(ispin), &
757  soft_valid=gapw, &
758  total_rho=total_rho(ispin))
759  CALL pw_axpy(trho_g(ispin), rho_tot_gspace)
760  IF (gapw_xc) THEN
761  CALL calculate_rho_elec(ks_env=ks_env, matrix_p=matrix_pe(ispin)%matrix, &
762  rho=trho_xc_r(ispin), &
763  rho_gspace=trho_xc_g(ispin), &
764  soft_valid=gapw_xc, &
765  total_rho=rhotot)
766  END IF
767  END DO
768 
769  ! GAPW o GAPW_XC require the calculation of hard and soft local densities
770  IF (gapw .OR. gapw_xc) THEN
771  CALL get_qs_env(qs_env=qs_env, oce=oce, sab_orb=sab)
772  CALL calculate_rho_atom_coeff(qs_env, matrix_pe, local_rho_set%rho_atom_set, &
773  qs_kind_set, oce, sab, para_env)
774  CALL prepare_gapw_den(qs_env, local_rho_set, do_rho0=gapw)
775  END IF
776  rhotot = sum(total_rho)
777  IF (gapw) THEN
778  CALL get_rho0_mpole(local_rho_set%rho0_mpole, qlm_tot=qlm_tot)
779  rhotot = rhotot + local_rho_set%rho0_mpole%total_rho0_h
780  CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rho_tot_gspace)
781  END IF
782 
783  IF (abs(rhotot) > 1.e-05_dp) THEN
784  logger => cp_get_default_logger()
785  IF (logger%para_env%is_source()) THEN
786  iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
787  ELSE
788  iounit = -1
789  END IF
790  cpwarn("Real space electron count of excitation density is non-zero.")
791  IF (iounit > 0) THEN
792  WRITE (iounit, "(T2,A,T61,G20.10)") "Measured electron count is ", rhotot
793  WRITE (iounit, "(T2,A,/)") repeat("*", 79)
794  END IF
795  END IF
796 
797  ! calculate associated hartree potential
798  CALL pw_poisson_solve(poisson_env, rho_tot_gspace, thartree, &
799  v_hartree_gspace)
800  CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
801  CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
802  IF (gapw) THEN
803  CALL vh_1c_gg_integrals(qs_env, thartree, hartree_local%ecoul_1c, &
804  local_rho_set, para_env, tddft=.true.)
805  CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, &
806  calculate_forces=.false., &
807  local_rho_set=local_rho_set)
808  END IF
809 
810  ! Fxc*drho term
811  CALL get_qs_env(qs_env, rho=rho)
812  CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
813  !
814  CALL get_qs_env(qs_env, input=input)
815  IF (dft_control%do_admm) THEN
816  CALL get_qs_env(qs_env, admm_env=admm_env)
817  xc_section => admm_env%xc_section_primary
818  ELSE
819  xc_section => section_vals_get_subs_vals(input, "DFT%XC")
820  END IF
821  !
822  deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
823  IF (deriv2_analytic) THEN
824  NULLIFY (v_xc, v_xc_tau, tau_r)
825  IF (gapw_xc) THEN
826  CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc)
827  CALL qs_fxc_analytic(rho_xc, trho_xc_r, tau_r, xc_section, auxbas_pw_pool, .false., v_xc, v_xc_tau)
828  ELSE
829  CALL qs_fxc_analytic(rho, trho_r, tau_r, xc_section, auxbas_pw_pool, .false., v_xc, v_xc_tau)
830  END IF
831  IF (gapw .OR. gapw_xc) THEN
832  CALL get_qs_env(qs_env, rho_atom_set=rho_atom_set)
833  rho1_atom_set => local_rho_set%rho_atom_set
834  CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
835  do_tddft=.true., do_triplet=.false.)
836  END IF
837  ELSE
838  cpabort("NYA 00006")
839  NULLIFY (v_xc, trho)
840  ALLOCATE (trho)
841  CALL qs_rho_create(trho)
842  CALL qs_rho_set(trho, rho_r=trho_r, rho_g=trho_g)
843  CALL qs_fxc_fdiff(ks_env, rho, trho, xc_section, 6, .false., v_xc, v_xc_tau)
844  DEALLOCATE (trho)
845  END IF
846 
847  DO ispin = 1, nspins
848  CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
849  CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
850  END DO
851  IF (gapw_xc) THEN
852  DO ispin = 1, nspins
853  CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_hartree_rspace, &
854  hmat=matrix_hz(ispin), &
855  calculate_forces=.false.)
856  CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
857  hmat=matrix_hz(ispin), &
858  gapw=gapw_xc, calculate_forces=.false.)
859  END DO
860  ELSE
861  ! vtot = v_xc(ispin) + v_hartree
862  DO ispin = 1, nspins
863  CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
864  hmat=matrix_hz(ispin), &
865  gapw=gapw, calculate_forces=.false.)
866  CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_hartree_rspace, &
867  hmat=matrix_hz(ispin), &
868  gapw=gapw, calculate_forces=.false.)
869  END DO
870  END IF
871  IF (gapw .OR. gapw_xc) THEN
872  mhz(1:nspins, 1:1) => matrix_hz(1:nspins)
873  mpe(1:nspins, 1:1) => matrix_pe(1:nspins)
874  CALL update_ks_atom(qs_env, mhz, mpe, forces=.false., &
875  rho_atom_external=local_rho_set%rho_atom_set)
876  END IF
877 
878  CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
879  CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
880  CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
881  DO ispin = 1, nspins
882  CALL auxbas_pw_pool%give_back_pw(trho_r(ispin))
883  CALL auxbas_pw_pool%give_back_pw(trho_g(ispin))
884  CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
885  END DO
886  DEALLOCATE (trho_r, trho_g, v_xc)
887  IF (gapw_xc) THEN
888  DO ispin = 1, nspins
889  CALL auxbas_pw_pool%give_back_pw(trho_xc_r(ispin))
890  CALL auxbas_pw_pool%give_back_pw(trho_xc_g(ispin))
891  END DO
892  DEALLOCATE (trho_xc_r, trho_xc_g)
893  END IF
894  IF (ASSOCIATED(v_xc_tau)) THEN
895  DO ispin = 1, nspins
896  CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
897  END DO
898  DEALLOCATE (v_xc_tau)
899  END IF
900  IF (dft_control%do_admm) THEN
901  IF (qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none) THEN
902  ! nothing to do
903  ELSE
904  ! add ADMM xc_section_aux terms: f_x[rhoz_ADMM]
905  CALL get_qs_env(qs_env, admm_env=admm_env)
906  CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit, matrix_s_aux_fit=msaux, &
907  task_list_aux_fit=task_list)
908  basis_type = "AUX_FIT"
909  !
910  NULLIFY (mpe, mhz)
911  ALLOCATE (mpe(nspins, 1))
912  CALL dbcsr_allocate_matrix_set(mhz, nspins, 1)
913  DO ispin = 1, nspins
914  ALLOCATE (mhz(ispin, 1)%matrix)
915  CALL dbcsr_create(mhz(ispin, 1)%matrix, template=msaux(1)%matrix)
916  CALL dbcsr_copy(mhz(ispin, 1)%matrix, msaux(1)%matrix)
917  CALL dbcsr_set(mhz(ispin, 1)%matrix, 0.0_dp)
918  mpe(ispin, 1)%matrix => matrix_pe_admm(ispin)%matrix
919  END DO
920  !
921  ! GAPW/GAPW_XC initializations
922  NULLIFY (local_rho_set_admm)
923  IF (admm_env%do_gapw) THEN
924  basis_type = "AUX_FIT_SOFT"
925  task_list => admm_env%admm_gapw_env%task_list
926  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
927  CALL get_admm_env(admm_env, sab_aux_fit=sab_aux_fit)
928  CALL local_rho_set_create(local_rho_set_admm)
929  CALL allocate_rho_atom_internals(local_rho_set_admm%rho_atom_set, atomic_kind_set, &
930  admm_env%admm_gapw_env%admm_kind_set, dft_control, para_env)
931  CALL calculate_rho_atom_coeff(qs_env, matrix_pe_admm, &
932  rho_atom_set=local_rho_set_admm%rho_atom_set, &
933  qs_kind_set=admm_env%admm_gapw_env%admm_kind_set, &
934  oce=admm_env%admm_gapw_env%oce, sab=sab_aux_fit, para_env=para_env)
935  CALL prepare_gapw_den(qs_env, local_rho_set=local_rho_set_admm, &
936  do_rho0=.false., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
937  END IF
938  !
939  xc_section => admm_env%xc_section_aux
940  !
941  NULLIFY (rho_g_aux, rho_r_aux, rhoz_g_aux, rhoz_r_aux)
942  CALL qs_rho_get(rho_aux_fit, rho_r=rho_r_aux, rho_g=rho_g_aux)
943  ! rhoz_aux
944  ALLOCATE (rhoz_r_aux(nspins), rhoz_g_aux(nspins))
945  DO ispin = 1, nspins
946  CALL auxbas_pw_pool%create_pw(rhoz_r_aux(ispin))
947  CALL auxbas_pw_pool%create_pw(rhoz_g_aux(ispin))
948  END DO
949  DO ispin = 1, nspins
950  CALL calculate_rho_elec(ks_env=ks_env, matrix_p=mpe(ispin, 1)%matrix, &
951  rho=rhoz_r_aux(ispin), rho_gspace=rhoz_g_aux(ispin), &
952  basis_type=basis_type, &
953  task_list_external=task_list)
954  END DO
955  !
956  NULLIFY (v_xc)
957  deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
958  IF (deriv2_analytic) THEN
959  NULLIFY (tau_r)
960  CALL qs_fxc_analytic(rho_aux_fit, rhoz_r_aux, tau_r, xc_section, auxbas_pw_pool, .false., v_xc, v_xc_tau)
961  ELSE
962  cpabort("NYA 00007")
963  NULLIFY (rhoz_aux)
964  ALLOCATE (rhoz_aux)
965  CALL qs_rho_create(rhoz_aux)
966  CALL qs_rho_set(rhoz_aux, rho_r=rhoz_r_aux, rho_g=rhoz_g_aux)
967  CALL qs_fxc_fdiff(ks_env, rho_aux_fit, rhoz_aux, xc_section, 6, .false., v_xc, v_xc_tau)
968  DEALLOCATE (rhoz_aux)
969  END IF
970  !
971  DO ispin = 1, nspins
972  CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
973  CALL integrate_v_rspace(qs_env=qs_env, v_rspace=v_xc(ispin), &
974  hmat=mhz(ispin, 1), basis_type=basis_type, &
975  calculate_forces=.false., &
976  task_list_external=task_list)
977  END DO
978  DO ispin = 1, nspins
979  CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
980  CALL auxbas_pw_pool%give_back_pw(rhoz_r_aux(ispin))
981  CALL auxbas_pw_pool%give_back_pw(rhoz_g_aux(ispin))
982  END DO
983  DEALLOCATE (v_xc, rhoz_r_aux, rhoz_g_aux)
984  !
985  IF (admm_env%do_gapw) THEN
986  rho_atom_set => admm_env%admm_gapw_env%local_rho_set%rho_atom_set
987  rho1_atom_set => local_rho_set_admm%rho_atom_set
988  CALL calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, &
989  para_env, kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
990  CALL update_ks_atom(qs_env, mhz(:, 1), matrix_pe_admm, forces=.false., tddft=.false., &
991  rho_atom_external=rho1_atom_set, &
992  kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
993  oce_external=admm_env%admm_gapw_env%oce, &
994  sab_external=sab_aux_fit)
995  END IF
996  !
997  nao = admm_env%nao_orb
998  nao_aux = admm_env%nao_aux_fit
999  ALLOCATE (dbwork)
1000  CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
1001  DO ispin = 1, nspins
1002  CALL cp_dbcsr_sm_fm_multiply(mhz(ispin, 1)%matrix, admm_env%A, &
1003  admm_env%work_aux_orb, nao)
1004  CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
1005  1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
1006  admm_env%work_orb_orb)
1007  CALL dbcsr_copy(dbwork, matrix_hz(1)%matrix)
1008  CALL dbcsr_set(dbwork, 0.0_dp)
1009  CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.true.)
1010  CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
1011  END DO
1012  CALL dbcsr_release(dbwork)
1013  DEALLOCATE (dbwork)
1014  CALL dbcsr_deallocate_matrix_set(mhz)
1015  DEALLOCATE (mpe)
1016  IF (admm_env%do_gapw) THEN
1017  IF (ASSOCIATED(local_rho_set_admm)) CALL local_rho_set_release(local_rho_set_admm)
1018  END IF
1019  END IF
1020  END IF
1021  IF (gapw .OR. gapw_xc) THEN
1022  IF (ASSOCIATED(local_rho_set)) CALL local_rho_set_release(local_rho_set)
1023  IF (ASSOCIATED(hartree_local)) CALL hartree_local_release(hartree_local)
1024  END IF
1025 
1026  ! HFX
1027  hfx_section => section_vals_get_subs_vals(xc_section, "HF")
1028  CALL section_vals_get(hfx_section, explicit=do_hfx)
1029  IF (do_hfx) THEN
1030  CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
1031  cpassert(n_rep_hf == 1)
1032  CALL section_vals_val_get(hfx_section, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
1033  i_rep_section=1)
1034  mspin = 1
1035  IF (hfx_treat_lsd_in_core) mspin = nspins
1036  !
1037  CALL get_qs_env(qs_env=qs_env, rho=rho, x_data=x_data, para_env=para_env, &
1038  s_mstruct_changed=s_mstruct_changed)
1039  distribute_fock_matrix = .true.
1040  IF (dft_control%do_admm) THEN
1041  CALL get_qs_env(qs_env, admm_env=admm_env)
1042  CALL get_admm_env(admm_env, matrix_s_aux_fit=msaux)
1043  NULLIFY (mpe, mhz)
1044  ALLOCATE (mpe(nspins, 1))
1045  CALL dbcsr_allocate_matrix_set(mhz, nspins, 1)
1046  DO ispin = 1, nspins
1047  ALLOCATE (mhz(ispin, 1)%matrix)
1048  CALL dbcsr_create(mhz(ispin, 1)%matrix, template=msaux(1)%matrix)
1049  CALL dbcsr_copy(mhz(ispin, 1)%matrix, msaux(1)%matrix)
1050  CALL dbcsr_set(mhz(ispin, 1)%matrix, 0.0_dp)
1051  mpe(ispin, 1)%matrix => matrix_pe_admm(ispin)%matrix
1052  END DO
1053  IF (x_data(1, 1)%do_hfx_ri) THEN
1054  eh1 = 0.0_dp
1055  CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpe, &
1056  geometry_did_change=s_mstruct_changed, nspins=nspins, &
1057  hf_fraction=x_data(1, 1)%general_parameter%fraction)
1058  ELSE
1059  DO ispin = 1, mspin
1060  eh1 = 0.0
1061  CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpe, hfx_section, &
1062  para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
1063  ispin=ispin)
1064  END DO
1065  END IF
1066  !
1067  cpassert(ASSOCIATED(admm_env%work_aux_orb))
1068  cpassert(ASSOCIATED(admm_env%work_orb_orb))
1069  nao = admm_env%nao_orb
1070  nao_aux = admm_env%nao_aux_fit
1071  ALLOCATE (dbwork)
1072  CALL dbcsr_create(dbwork, template=matrix_hz(1)%matrix)
1073  DO ispin = 1, nspins
1074  CALL cp_dbcsr_sm_fm_multiply(mhz(ispin, 1)%matrix, admm_env%A, &
1075  admm_env%work_aux_orb, nao)
1076  CALL parallel_gemm('T', 'N', nao, nao, nao_aux, &
1077  1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
1078  admm_env%work_orb_orb)
1079  CALL dbcsr_copy(dbwork, matrix_hz(ispin)%matrix)
1080  CALL dbcsr_set(dbwork, 0.0_dp)
1081  CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbwork, keep_sparsity=.true.)
1082  CALL dbcsr_add(matrix_hz(ispin)%matrix, dbwork, 1.0_dp, 1.0_dp)
1083  END DO
1084  CALL dbcsr_release(dbwork)
1085  DEALLOCATE (dbwork)
1086  CALL dbcsr_deallocate_matrix_set(mhz)
1087  DEALLOCATE (mpe)
1088  ELSE
1089  NULLIFY (mpe, mhz)
1090  ALLOCATE (mpe(nspins, 1), mhz(nspins, 1))
1091  DO ispin = 1, nspins
1092  mhz(ispin, 1)%matrix => matrix_hz(ispin)%matrix
1093  mpe(ispin, 1)%matrix => matrix_pe(ispin)%matrix
1094  END DO
1095  IF (x_data(1, 1)%do_hfx_ri) THEN
1096  eh1 = 0.0_dp
1097  CALL hfx_ri_update_ks(qs_env, x_data(1, 1)%ri_data, mhz, eh1, rho_ao=mpe, &
1098  geometry_did_change=s_mstruct_changed, nspins=nspins, &
1099  hf_fraction=x_data(1, 1)%general_parameter%fraction)
1100  ELSE
1101  DO ispin = 1, mspin
1102  eh1 = 0.0
1103  CALL integrate_four_center(qs_env, x_data, mhz, eh1, mpe, hfx_section, &
1104  para_env, s_mstruct_changed, 1, distribute_fock_matrix, &
1105  ispin=ispin)
1106  END DO
1107  END IF
1108  DEALLOCATE (mpe, mhz)
1109  END IF
1110  END IF
1111 
1112  focc = 4.0_dp
1113  IF (nspins == 2) focc = 2.0_dp
1114  DO ispin = 1, nspins
1115  mos => gs_mos(ispin)%mos_occ
1116  CALL cp_fm_get_info(mos, ncol_global=norb)
1117  CALL cp_dbcsr_sm_fm_multiply(matrix_hz(ispin)%matrix, mos, cpmos(ispin), &
1118  norb, alpha=focc, beta=0.0_dp)
1119  END DO
1120 
1121  CALL timestop(handle)
1122 
1123  END SUBROUTINE tddfpt_resvec2
1124 
1125 ! **************************************************************************************************
1126 !> \brief ...
1127 !> \param qs_env ...
1128 !> \param matrix_pe ...
1129 !> \param gs_mos ...
1130 !> \param matrix_hz ...
1131 !> \param cpmos ...
1132 ! **************************************************************************************************
1133  SUBROUTINE tddfpt_resvec2_xtb(qs_env, matrix_pe, gs_mos, matrix_hz, cpmos)
1134 
1135  TYPE(qs_environment_type), POINTER :: qs_env
1136  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_pe
1137  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1138  POINTER :: gs_mos
1139  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_hz
1140  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: cpmos
1141 
1142  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_resvec2_xtb'
1143 
1144  INTEGER :: atom_a, handle, iatom, ikind, is, ispin, &
1145  na, natom, natorb, nkind, norb, ns, &
1146  nsgf, nspins
1147  INTEGER, DIMENSION(25) :: lao
1148  INTEGER, DIMENSION(5) :: occ
1149  REAL(dp), ALLOCATABLE, DIMENSION(:) :: mcharge, mcharge1
1150  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, charges, charges1
1151  REAL(kind=dp) :: focc
1152  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1153  TYPE(cp_fm_type), POINTER :: mos
1154  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
1155  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_s
1156  TYPE(dbcsr_type), POINTER :: s_matrix
1157  TYPE(dft_control_type), POINTER :: dft_control
1158  TYPE(mp_para_env_type), POINTER :: para_env
1159  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1160  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1161  TYPE(qs_rho_type), POINTER :: rho
1162  TYPE(xtb_atom_type), POINTER :: xtb_kind
1163 
1164  CALL timeset(routinen, handle)
1165 
1166  cpassert(ASSOCIATED(matrix_pe))
1167 
1168  CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
1169  nspins = dft_control%nspins
1170 
1171  DO ispin = 1, nspins
1172  CALL dbcsr_set(matrix_hz(ispin)%matrix, 0.0_dp)
1173  END DO
1174 
1175  IF (dft_control%qs_control%xtb_control%coulomb_interaction) THEN
1176  ! Mulliken charges
1177  CALL get_qs_env(qs_env, rho=rho, particle_set=particle_set, &
1178  matrix_s_kp=matrix_s, para_env=para_env)
1179  natom = SIZE(particle_set)
1180  CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
1181  ALLOCATE (mcharge(natom), charges(natom, 5))
1182  ALLOCATE (mcharge1(natom), charges1(natom, 5))
1183  charges = 0.0_dp
1184  charges1 = 0.0_dp
1185  CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
1186  nkind = SIZE(atomic_kind_set)
1187  CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
1188  ALLOCATE (aocg(nsgf, natom))
1189  aocg = 0.0_dp
1190  ALLOCATE (aocg1(nsgf, natom))
1191  aocg1 = 0.0_dp
1192  p_matrix => matrix_p(:, 1)
1193  s_matrix => matrix_s(1, 1)%matrix
1194  CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
1195  CALL ao_charges(matrix_pe, s_matrix, aocg1, para_env)
1196  DO ikind = 1, nkind
1197  CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
1198  CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
1199  CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, occupation=occ)
1200  DO iatom = 1, na
1201  atom_a = atomic_kind_set(ikind)%atom_list(iatom)
1202  charges(atom_a, :) = real(occ(:), kind=dp)
1203  DO is = 1, natorb
1204  ns = lao(is) + 1
1205  charges(atom_a, ns) = charges(atom_a, ns) - aocg(is, atom_a)
1206  charges1(atom_a, ns) = charges1(atom_a, ns) - aocg1(is, atom_a)
1207  END DO
1208  END DO
1209  END DO
1210  DEALLOCATE (aocg, aocg1)
1211  DO iatom = 1, natom
1212  mcharge(iatom) = sum(charges(iatom, :))
1213  mcharge1(iatom) = sum(charges1(iatom, :))
1214  END DO
1215  ! Coulomb Kernel
1216  CALL xtb_coulomb_hessian(qs_env, matrix_hz, charges1, mcharge1, mcharge)
1217  !
1218  DEALLOCATE (charges, mcharge, charges1, mcharge1)
1219  END IF
1220 
1221  focc = 2.0_dp
1222  IF (nspins == 2) focc = 1.0_dp
1223  DO ispin = 1, nspins
1224  mos => gs_mos(ispin)%mos_occ
1225  CALL cp_fm_get_info(mos, ncol_global=norb)
1226  CALL cp_dbcsr_sm_fm_multiply(matrix_hz(ispin)%matrix, mos, cpmos(ispin), &
1227  norb, alpha=focc, beta=0.0_dp)
1228  END DO
1229 
1230  CALL timestop(handle)
1231 
1232  END SUBROUTINE tddfpt_resvec2_xtb
1233 
1234 ! **************************************************************************************************
1235 !> \brief ...
1236 !> \param qs_env ...
1237 !> \param cpmos ...
1238 !> \param work ...
1239 ! **************************************************************************************************
1240  SUBROUTINE tddfpt_resvec3(qs_env, cpmos, work)
1241 
1242  TYPE(qs_environment_type), POINTER :: qs_env
1243  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: cpmos
1244  TYPE(tddfpt_work_matrices) :: work
1245 
1246  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_resvec3'
1247 
1248  INTEGER :: handle, ispin, nao, norb, nspins
1249  TYPE(cp_fm_struct_type), POINTER :: fmstruct
1250  TYPE(cp_fm_type) :: cvec, umat
1251  TYPE(cp_fm_type), POINTER :: omos
1252  TYPE(dft_control_type), POINTER :: dft_control
1253  TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1254 
1255  CALL timeset(routinen, handle)
1256 
1257  CALL get_qs_env(qs_env, mos=mos, dft_control=dft_control)
1258  nspins = dft_control%nspins
1259 
1260  DO ispin = 1, nspins
1261  CALL get_mo_set(mos(ispin), mo_coeff=omos)
1262  associate(rvecs => cpmos(ispin))
1263  CALL cp_fm_get_info(rvecs, nrow_global=nao, ncol_global=norb)
1264  CALL cp_fm_create(cvec, rvecs%matrix_struct, "cvec")
1265  CALL cp_fm_struct_create(fmstruct, context=rvecs%matrix_struct%context, nrow_global=norb, &
1266  ncol_global=norb, para_env=rvecs%matrix_struct%para_env)
1267  CALL cp_fm_create(umat, fmstruct, "umat")
1268  CALL cp_fm_struct_release(fmstruct)
1269  !
1270  CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, omos, work%S_C0(ispin), 0.0_dp, umat)
1271  CALL cp_fm_copy_general(rvecs, cvec, rvecs%matrix_struct%para_env)
1272  CALL parallel_gemm("N", "T", nao, norb, norb, 1.0_dp, cvec, umat, 0.0_dp, rvecs)
1273  END associate
1274  CALL cp_fm_release(cvec)
1275  CALL cp_fm_release(umat)
1276  END DO
1277 
1278  CALL timestop(handle)
1279 
1280  END SUBROUTINE tddfpt_resvec3
1281 
1282 ! **************************************************************************************************
1283 !> \brief Calculate direct tddft forces
1284 !> \param qs_env ...
1285 !> \param ex_env ...
1286 !> \param gs_mos ...
1287 !> \param kernel_env ...
1288 !> \param sub_env ...
1289 !> \param work_matrices ...
1290 !> \param debug_forces ...
1291 !> \par History
1292 !> * 01.2020 screated [JGH]
1293 ! **************************************************************************************************
1294  SUBROUTINE tddfpt_kernel_force(qs_env, ex_env, gs_mos, kernel_env, sub_env, work_matrices, debug_forces)
1295 
1296  TYPE(qs_environment_type), POINTER :: qs_env
1297  TYPE(excited_energy_type), POINTER :: ex_env
1298  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
1299  POINTER :: gs_mos
1300  TYPE(kernel_env_type), INTENT(IN) :: kernel_env
1301  TYPE(tddfpt_subgroup_env_type) :: sub_env
1302  TYPE(tddfpt_work_matrices) :: work_matrices
1303  LOGICAL, INTENT(IN) :: debug_forces
1304 
1305  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_kernel_force'
1306 
1307  INTEGER :: handle
1308  TYPE(dft_control_type), POINTER :: dft_control
1309  TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
1310 
1311  mark_used(work_matrices)
1312 
1313  CALL timeset(routinen, handle)
1314 
1315  CALL get_qs_env(qs_env, dft_control=dft_control)
1316  tddfpt_control => dft_control%tddfpt2_control
1317 
1318  IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
1319  ! full Kernel
1320  CALL fhxc_force(qs_env, ex_env, gs_mos, kernel_env%full_kernel, debug_forces)
1321  ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
1322  ! sTDA Kernel
1323  CALL stda_force(qs_env, ex_env, gs_mos, kernel_env%stda_kernel, sub_env, work_matrices, debug_forces)
1324  ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
1325  ! nothing to be done here
1326  ex_env%matrix_wx1 => null()
1327  ELSE
1328  cpabort('Unknown kernel type')
1329  END IF
1330 
1331  CALL timestop(handle)
1332 
1333  END SUBROUTINE tddfpt_kernel_force
1334 
1335 END MODULE qs_tddfpt2_forces
Types and set/get functions for auxiliary density matrix methods.
Definition: admm_types.F:15
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
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
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.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
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:125
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
Definition: cp_fm_types.F:1538
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
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:535
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
various routines to log and control the output. The idea is that decisions about where to log should ...
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
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Types for excited states potential energies.
subroutine, public exstate_potential_release(ex_env)
...
subroutine, public init_coulomb_local(hartree_local, natom)
...
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...
subroutine, public hartree_local_release(hartree_local)
...
subroutine, public hartree_local_create(hartree_local)
...
Routines to calculate HFX energy and potential.
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)
computes four center integrals for a full basis set and updates the Kohn-Sham-Matrix and energy....
RI-methods for HFX.
Definition: hfx_ri.F:12
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:1033
Types and set/get functions for HFX.
Definition: hfx_types.F:15
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public tddfpt_kernel_none
integer, parameter, public oe_shift
integer, parameter, public do_admm_aux_exch_func_none
integer, parameter, public tddfpt_kernel_full
integer, parameter, public tddfpt_kernel_stda
objects that represent the structure of input sections and the data contained in an input section
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
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
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
logical function, public section_get_lval(section_vals, keyword_name)
...
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Interface to the message passing library MPI.
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition: mulliken.F:13
basic linear algebra operations for full matrixes
Define the data structure for the particle information.
container for various plainwaves related things
Definition: pw_env_types.F:14
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
functions related to the poisson solver on regular grids
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
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
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
collects routines that calculate density matrices
subroutine, public calculate_wx_matrix(mos_occ, xvec, ks_matrix, w_matrix)
Calculate the excited state W matrix from the MO eigenvectors, KS matrix.
subroutine, public calculate_xwx_matrix(mos_occ, xvec, s_matrix, ks_matrix, w_matrix, eval)
Calculate the excited state W matrix from the MO eigenvectors, KS matrix.
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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, 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, 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, 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, rhs)
Set the QUICKSTEP environment.
subroutine, public sum_qs_force(qs_force_out, qs_force_in)
Sum up two qs_force entities qs_force_out = qs_force_out + qs_force_in.
subroutine, public deallocate_qs_force(qs_force)
Deallocate a Quickstep force data structure.
subroutine, public zero_qs_force(qs_force)
Initialize a Quickstep force data structure.
subroutine, public allocate_qs_force(qs_force, natom_of_kind)
Allocate a Quickstep force data structure.
subroutine, public total_qs_force(force, qs_force, atomic_kind_set)
Get current total force.
https://en.wikipedia.org/wiki/Finite_difference_coefficient
Definition: qs_fxc.F:27
subroutine, public qs_fxc_analytic(rho0, rho1_r, tau1_r, xc_section, auxbas_pw_pool, is_triplet, v_xc, v_xc_tau)
...
Definition: qs_fxc.F:85
subroutine, public qs_fxc_fdiff(ks_env, rho0_struct, rho1_struct, xc_section, accuracy, is_triplet, fxc_rho, fxc_tau)
...
Definition: qs_fxc.F:146
subroutine, public prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external)
...
Integrate single or product functions over a potential on a RS grid.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, 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_r3d_rs_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, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
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)
Get attributes of an atomic kind set.
routines that build the Kohn-Sham matrix contributions coming from local atomic densities
Definition: qs_ks_atom.F:12
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
Calculate the KS reference potentials.
subroutine, public ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress)
calculate the Kohn-Sham reference potential
subroutine, public ks_ref_potential_atom(qs_env, local_rho_set, local_rho_set_admm, v_hartree_rspace)
calculate the Kohn-Sham GAPW reference potentials
subroutine, public local_rho_set_create(local_rho_set)
...
subroutine, public local_rho_set_release(local_rho_set)
...
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
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:397
Define the neighbor list data types and the corresponding functionality.
Calculation of overlap matrix, its derivatives and forces.
Definition: qs_overlap.F:19
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
subroutine, public integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, local_rho_set_2nd, atener, kforce)
...
subroutine, public rho0_s_grid_create(pw_env, rho0_mpole)
...
subroutine, public init_rho0(local_rho_set, qs_env, gapw_control, zcore)
...
subroutine, public get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, mp_gau_ikind, mp_rho, norm_g0l_h, Qlm_gg, Qlm_car, Qlm_tot, zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, max_rpgf0_s, rho0_s_rs, rho0_s_gs)
...
subroutine, public allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
...
subroutine, public calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
subroutine, public qs_rho_set(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
...
Definition: qs_rho_types.F:308
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
subroutine, public qs_rho_create(rho)
Allocates a new instance of rho.
Definition: qs_rho_types.F:99
subroutine, public stda_force(qs_env, ex_env, gs_mos, stda_env, sub_env, work, debug_forces)
Simplified Tamm Dancoff approach (sTDA). Kernel contribution to forces.
subroutine, public fhxc_force(qs_env, ex_env, gs_mos, full_kernel, debug_forces)
Calculate direct tddft forces.
subroutine, public tddfpt_forces_main(qs_env, gs_mos, ex_env, kernel_env, sub_env, work_matrices)
Perform TDDFPT gradient calculation.
routines that build the integrals of the Vxc potential calculated for the atomic density in the basis...
Definition: qs_vxc_atom.F:12
subroutine, public calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, do_tddft, do_tddfpt2, do_triplet, kind_set_external)
...
Definition: qs_vxc_atom.F:446
types for task lists
Calculation of Coulomb Hessian contributions in xTB.
Definition: xtb_ehess.F:12
subroutine, public xtb_coulomb_hessian(qs_env, ks_matrix, charges1, mcharge1, mcharge)
...
Definition: xtb_ehess.F:77
Definition of the xTB parameter types.
Definition: xtb_types.F:20
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, occupation, electronegativity, chmax)
...
Definition: xtb_types.F:175