(git:0de0cc2)
qs_tddfpt2_eigensolver.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 cp_blacs_env, ONLY: cp_blacs_env_type
10  USE cp_control_types, ONLY: tddfpt2_control_type
12  USE cp_fm_basic_linalg, ONLY: cp_fm_contracted_trace,&
13  cp_fm_scale,&
15  cp_fm_trace
21  cp_fm_struct_type
22  USE cp_fm_types, ONLY: &
24  cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type
25  USE cp_log_handling, ONLY: cp_logger_type
27  USE dbcsr_api, ONLY: dbcsr_get_info,&
28  dbcsr_p_type,&
29  dbcsr_type
33  USE input_section_types, ONLY: section_vals_type
34  USE kinds, ONLY: dp,&
35  int_8
36  USE machine, ONLY: m_flush,&
38  USE memory_utilities, ONLY: reallocate
39  USE message_passing, ONLY: mp_para_env_type
40  USE parallel_gemm_api, ONLY: parallel_gemm
41  USE physcon, ONLY: evolt
42  USE qs_environment_types, ONLY: get_qs_env,&
43  qs_environment_type
44  USE qs_kernel_types, ONLY: full_kernel_env_type,&
45  kernel_env_type
46  USE qs_scf_methods, ONLY: eigensolver
47  USE qs_tddfpt2_fhxc, ONLY: fhxc_kernel,&
54  USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
55  USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos,&
56  tddfpt_work_matrices
58 #include "./base/base_uses.f90"
59 
60  IMPLICIT NONE
61 
62  PRIVATE
63 
64  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_eigensolver'
65 
66  LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
67  ! number of first derivative components (3: d/dx, d/dy, d/dz)
68  INTEGER, PARAMETER, PRIVATE :: nderivs = 3
69  INTEGER, PARAMETER, PRIVATE :: maxspins = 2
70 
73 
74 ! **************************************************************************************************
75 
76 CONTAINS
77 
78 ! **************************************************************************************************
79 !> \brief Make TDDFPT trial vectors orthogonal to all occupied molecular orbitals.
80 !> \param evects trial vectors distributed across all processors (modified on exit)
81 !> \param S_C0_C0T matrix product S * C_0 * C_0^T, where C_0 is the ground state
82 !> wave function for each spin expressed in atomic basis set,
83 !> and S is the corresponding overlap matrix
84 !> \par History
85 !> * 05.2016 created [Sergey Chulkov]
86 !> * 05.2019 use a temporary work matrix [JHU]
87 !> \note Based on the subroutine p_preortho() which was created by Thomas Chassaing on 09.2002.
88 !> Should be useless when ground state MOs are computed with extremely high accuracy,
89 !> as all virtual orbitals are already orthogonal to the occupied ones by design.
90 !> However, when the norm of residual vectors is relatively small (e.g. less then SCF_EPS),
91 !> new Krylov's vectors seem to be random and should be orthogonalised even with respect to
92 !> the occupied MOs.
93 ! **************************************************************************************************
94  SUBROUTINE tddfpt_orthogonalize_psi1_psi0(evects, S_C0_C0T)
95  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
96  TYPE(cp_fm_type), DIMENSION(:), INTENT(in) :: s_c0_c0t
97 
98  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_orthogonalize_psi1_psi0'
99 
100  INTEGER :: handle, ispin, ivect, nactive, nao, &
101  nspins, nvects
102  TYPE(cp_fm_struct_type), POINTER :: matrix_struct
103  TYPE(cp_fm_type) :: evortho
104 
105  CALL timeset(routinen, handle)
106 
107  nspins = SIZE(evects, 1)
108  nvects = SIZE(evects, 2)
109 
110  IF (nvects > 0) THEN
111  DO ispin = 1, nspins
112  CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
113  nrow_global=nao, ncol_global=nactive)
114  CALL cp_fm_create(evortho, matrix_struct)
115  DO ivect = 1, nvects
116  ! evortho: C0 * C0^T * S * C1 == (S * C0 * C0^T)^T * C1
117  CALL parallel_gemm('T', 'N', nao, nactive, nao, 1.0_dp, s_c0_c0t(ispin), &
118  evects(ispin, ivect), 0.0_dp, evortho)
119  CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, ivect), -1.0_dp, evortho)
120  END DO
121  CALL cp_fm_release(evortho)
122  END DO
123  END IF
124 
125  CALL timestop(handle)
126 
127  END SUBROUTINE tddfpt_orthogonalize_psi1_psi0
128 
129 ! **************************************************************************************************
130 !> \brief Check that orthogonalised TDDFPT trial vectors remain orthogonal to
131 !> occupied molecular orbitals.
132 !> \param evects trial vectors
133 !> \param S_C0 matrix product S * C_0, where C_0 is the ground state wave function
134 !> for each spin in atomic basis set, and S is the corresponding overlap matrix
135 !> \param max_norm the largest possible overlap between the ground state and
136 !> excited state wave functions
137 !> \return true if trial vectors are non-orthogonal to occupied molecular orbitals
138 !> \par History
139 !> * 07.2016 created [Sergey Chulkov]
140 !> * 05.2019 use temporary work matrices [JHU]
141 ! **************************************************************************************************
142  FUNCTION tddfpt_is_nonorthogonal_psi1_psi0(evects, S_C0, max_norm) &
143  result(is_nonortho)
144  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
145  TYPE(cp_fm_type), DIMENSION(:), INTENT(in) :: s_c0
146  REAL(kind=dp), INTENT(in) :: max_norm
147  LOGICAL :: is_nonortho
148 
149  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_is_nonorthogonal_psi1_psi0'
150 
151  INTEGER :: handle, ispin, ivect, nactive, nao, &
152  nocc, nspins, nvects
153  REAL(kind=dp) :: maxabs_val
154  TYPE(cp_fm_struct_type), POINTER :: matrix_struct, matrix_struct_tmp
155  TYPE(cp_fm_type) :: aortho
156 
157  CALL timeset(routinen, handle)
158 
159  nspins = SIZE(evects, 1)
160  nvects = SIZE(evects, 2)
161 
162  is_nonortho = .false.
163 
164  loop: DO ispin = 1, nspins
165  CALL cp_fm_get_info(matrix=s_c0(ispin), ncol_global=nocc)
166  CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
167  nrow_global=nao, ncol_global=nactive)
168  CALL cp_fm_struct_create(matrix_struct_tmp, nrow_global=nocc, &
169  ncol_global=nactive, template_fmstruct=matrix_struct)
170  CALL cp_fm_create(aortho, matrix_struct_tmp)
171  CALL cp_fm_struct_release(matrix_struct_tmp)
172  DO ivect = 1, nvects
173  ! aortho = S_0^T * S * C_1
174  CALL parallel_gemm('T', 'N', nocc, nactive, nao, 1.0_dp, s_c0(ispin), &
175  evects(ispin, ivect), 0.0_dp, aortho)
176  CALL cp_fm_maxabsval(aortho, maxabs_val)
177  is_nonortho = maxabs_val > max_norm
178  IF (is_nonortho) THEN
179  CALL cp_fm_release(aortho)
180  EXIT loop
181  END IF
182  END DO
183  CALL cp_fm_release(aortho)
184  END DO loop
185 
186  CALL timestop(handle)
187 
188  END FUNCTION tddfpt_is_nonorthogonal_psi1_psi0
189 
190 ! **************************************************************************************************
191 !> \brief Make new TDDFPT trial vectors orthonormal to all previous trial vectors.
192 !> \param evects trial vectors (modified on exit)
193 !> \param nvects_new number of new trial vectors to orthogonalise
194 !> \param S_evects set of matrices to store matrix product S * evects (modified on exit)
195 !> \param matrix_s overlap matrix
196 !> \par History
197 !> * 05.2016 created [Sergey Chulkov]
198 !> * 02.2017 caching the matrix product S * evects [Sergey Chulkov]
199 !> \note \parblock
200 !> Based on the subroutines reorthogonalize() and normalize() which were originally created
201 !> by Thomas Chassaing on 03.2003.
202 !>
203 !> In order to orthogonalise a trial vector C3 = evects(:,3) with respect to previously
204 !> orthogonalised vectors C1 = evects(:,1) and C2 = evects(:,2) we need to compute the
205 !> quantity C3'' using the following formulae:
206 !> C3' = C3 - Tr(C3^T * S * C1) * C1,
207 !> C3'' = C3' - Tr(C3'^T * S * C2) * C2,
208 !> which can be expanded as:
209 !> C3'' = C3 - Tr(C3^T * S * C1) * C1 - Tr(C3^T * S * C2) * C2 +
210 !> Tr(C3^T * S * C1) * Tr(C2^T * S * C1) * C2 .
211 !> In case of unlimited float-point precision, the last term in above expression is exactly 0,
212 !> due to orthogonality condition between C1 and C2. In this case the expression could be
213 !> simplified as (taking into account the identity: Tr(A * S * B) = Tr(B * S * A)):
214 !> C3'' = C3 - Tr(C1^T * S * C3) * C1 - Tr(C2^T * S * C3) * C2 ,
215 !> which means we do not need the variable S_evects to keep the matrix products S * Ci .
216 !>
217 !> In reality, however, we deal with limited float-point precision arithmetic meaning that
218 !> the trace Tr(C2^T * S * C1) is close to 0 but does not equal to 0 exactly. The term
219 !> Tr(C3^T * S * C1) * Tr(C2^T * S * C1) * C2
220 !> can not be ignored anymore. Ignorance of this term will lead to numerical instability
221 !> when the trace Tr(C3^T * S * C1) is large enough.
222 !> \endparblock
223 ! **************************************************************************************************
224  SUBROUTINE tddfpt_orthonormalize_psi1_psi1(evects, nvects_new, S_evects, matrix_s)
225  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
226  INTEGER, INTENT(in) :: nvects_new
227  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: s_evects
228  TYPE(dbcsr_type), POINTER :: matrix_s
229 
230  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_orthonormalize_psi1_psi1'
231 
232  INTEGER :: handle, ispin, ivect, jvect, nspins, &
233  nvects_old, nvects_total
234  INTEGER, DIMENSION(maxspins) :: nactive
235  REAL(kind=dp) :: norm
236  REAL(kind=dp), DIMENSION(maxspins) :: weights
237 
238  CALL timeset(routinen, handle)
239 
240  nspins = SIZE(evects, 1)
241  nvects_total = SIZE(evects, 2)
242  nvects_old = nvects_total - nvects_new
243 
244  IF (debug_this_module) THEN
245  cpassert(SIZE(s_evects, 1) == nspins)
246  cpassert(SIZE(s_evects, 2) == nvects_total)
247  cpassert(nvects_old >= 0)
248  END IF
249 
250  DO ispin = 1, nspins
251  CALL cp_fm_get_info(matrix=evects(ispin, 1), ncol_global=nactive(ispin))
252  END DO
253 
254  DO jvect = nvects_old + 1, nvects_total
255  ! <psi1_i | psi1_j>
256  DO ivect = 1, jvect - 1
257  CALL cp_fm_trace(evects(:, jvect), s_evects(:, ivect), weights(1:nspins), accurate=.false.)
258  norm = sum(weights(1:nspins))
259 
260  DO ispin = 1, nspins
261  CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, jvect), -norm, evects(ispin, ivect))
262  END DO
263  END DO
264 
265  ! <psi1_j | psi1_j>
266  DO ispin = 1, nspins
267  CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, jvect), s_evects(ispin, jvect), &
268  ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
269  END DO
270 
271  CALL cp_fm_trace(evects(:, jvect), s_evects(:, jvect), weights(1:nspins), accurate=.false.)
272 
273  norm = sum(weights(1:nspins))
274  norm = 1.0_dp/sqrt(norm)
275 
276  DO ispin = 1, nspins
277  CALL cp_fm_scale(norm, evects(ispin, jvect))
278  CALL cp_fm_scale(norm, s_evects(ispin, jvect))
279  END DO
280  END DO
281 
282  CALL timestop(handle)
283 
284  END SUBROUTINE tddfpt_orthonormalize_psi1_psi1
285 
286 ! **************************************************************************************************
287 !> \brief Compute action matrix-vector products.
288 !> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
289 !> \param evects TDDFPT trial vectors
290 !> \param S_evects cached matrix product S * evects where S is the overlap matrix
291 !> in primary basis set
292 !> \param gs_mos molecular orbitals optimised for the ground state
293 !> \param tddfpt_control control section for tddfpt
294 !> \param matrix_ks Kohn-Sham matrix
295 !> \param qs_env Quickstep environment
296 !> \param kernel_env kernel environment
297 !> \param sub_env parallel (sub)group environment
298 !> \param work_matrices collection of work matrices (modified on exit)
299 !> \par History
300 !> * 06.2016 created [Sergey Chulkov]
301 !> * 03.2017 refactored [Sergey Chulkov]
302 ! **************************************************************************************************
303  SUBROUTINE tddfpt_compute_aop_evects(Aop_evects, evects, S_evects, gs_mos, tddfpt_control, &
304  matrix_ks, qs_env, kernel_env, &
305  sub_env, work_matrices)
306  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: aop_evects, evects, s_evects
307  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
308  INTENT(in) :: gs_mos
309  TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
310  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_ks
311  TYPE(qs_environment_type), POINTER :: qs_env
312  TYPE(kernel_env_type), INTENT(in) :: kernel_env
313  TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
314  TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices
315 
316  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_compute_Aop_evects'
317 
318  INTEGER :: handle, ispin, ivect, nspins, nvects
319  INTEGER, DIMENSION(maxspins) :: nmo_occ
320  LOGICAL :: do_admm, do_hfx, do_lri_response, &
321  is_rks_triplets, re_int
322  REAL(kind=dp) :: rcut, scale
323  TYPE(cp_fm_type) :: fm_dummy
324  TYPE(full_kernel_env_type), POINTER :: kernel_env_admm_aux
325  TYPE(mp_para_env_type), POINTER :: para_env
326 
327  CALL timeset(routinen, handle)
328 
329  nspins = SIZE(evects, 1)
330  nvects = SIZE(evects, 2)
331  do_hfx = tddfpt_control%do_hfx
332  do_admm = tddfpt_control%do_admm
333  IF (do_admm) THEN
334  kernel_env_admm_aux => kernel_env%admm_kernel
335  ELSE
336  NULLIFY (kernel_env_admm_aux)
337  END IF
338  is_rks_triplets = tddfpt_control%rks_triplets
339  do_lri_response = tddfpt_control%do_lrigpw
340 
341  IF (debug_this_module) THEN
342  cpassert(nspins > 0)
343  cpassert(SIZE(aop_evects, 1) == nspins)
344  cpassert(SIZE(aop_evects, 2) == nvects)
345  cpassert(SIZE(s_evects, 1) == nspins)
346  cpassert(SIZE(s_evects, 2) == nvects)
347  cpassert(SIZE(gs_mos) == nspins)
348  END IF
349 
350  DO ispin = 1, nspins
351  nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
352  END DO
353 
354  IF (nvects > 0) THEN
355  CALL cp_fm_get_info(evects(1, 1), para_env=para_env)
356  IF (ALLOCATED(work_matrices%evects_sub)) THEN
357  DO ivect = 1, nvects
358  DO ispin = 1, nspins
359  associate(evect => evects(ispin, ivect), work_matrix => work_matrices%evects_sub(ispin, ivect))
360  IF (ASSOCIATED(evect%matrix_struct)) THEN
361  IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
362  CALL cp_fm_copy_general(evect, work_matrix, para_env)
363  ELSE
364  CALL cp_fm_copy_general(evect, fm_dummy, para_env)
365  END IF
366  ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
367  CALL cp_fm_copy_general(fm_dummy, work_matrix, para_env)
368  ELSE
369  CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env)
370  END IF
371  END associate
372  END DO
373  END DO
374  END IF
375 
376  IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
377  ! full TDDFPT kernel
378  CALL fhxc_kernel(aop_evects, evects, is_rks_triplets, do_hfx, do_admm, qs_env, &
379  kernel_env%full_kernel, kernel_env_admm_aux, sub_env, work_matrices, &
380  tddfpt_control%admm_symm, tddfpt_control%admm_xc_correction, &
381  do_lri_response)
382  ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
383  ! sTDA kernel
384  CALL stda_kernel(aop_evects, evects, is_rks_triplets, qs_env, tddfpt_control%stda_control, &
385  kernel_env%stda_kernel, sub_env, work_matrices)
386  ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
387  ! No kernel
388  DO ivect = 1, nvects
389  DO ispin = 1, nspins
390  CALL cp_fm_set_all(aop_evects(ispin, ivect), 0.0_dp)
391  END DO
392  END DO
393  ELSE
394  cpabort("Kernel type not implemented")
395  END IF
396 
397  IF (ALLOCATED(work_matrices%evects_sub)) THEN
398  DO ivect = 1, nvects
399  DO ispin = 1, nspins
400  associate(aop_evect => aop_evects(ispin, ivect), &
401  work_matrix => work_matrices%Aop_evects_sub(ispin, ivect))
402  IF (ASSOCIATED(aop_evect%matrix_struct)) THEN
403  IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
404  CALL cp_fm_copy_general(work_matrix, aop_evect, para_env)
405  ELSE
406  CALL cp_fm_copy_general(fm_dummy, aop_evect, para_env)
407  END IF
408  ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
409  CALL cp_fm_copy_general(work_matrix, fm_dummy, para_env)
410  ELSE
411  CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env)
412  END IF
413  END associate
414  END DO
415  END DO
416  END IF
417 
418  ! orbital energy difference term
419  CALL tddfpt_apply_energy_diff(aop_evects=aop_evects, evects=evects, s_evects=s_evects, &
420  gs_mos=gs_mos, matrix_ks=matrix_ks)
421 
422  IF (do_hfx) THEN
423  IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
424  ! full TDDFPT kernel
425  CALL tddfpt_apply_hfx(aop_evects=aop_evects, evects=evects, gs_mos=gs_mos, do_admm=do_admm, &
426  qs_env=qs_env, wfm_rho_orb=work_matrices%hfx_fm_ao_ao, &
427  work_hmat_symm=work_matrices%hfx_hmat_symm, &
428  work_hmat_asymm=work_matrices%hfx_hmat_asymm, &
429  work_rho_ia_ao_symm=work_matrices%hfx_rho_ao_symm, &
430  work_rho_ia_ao_asymm=work_matrices%hfx_rho_ao_asymm)
431  ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
432  ! sTDA kernel
433  ! special treatment of HFX term
434  ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
435  ! No kernel
436  ! drop kernel contribution of HFX term
437  ELSE
438  cpabort("Kernel type not implemented")
439  END IF
440  END IF
441  ! short/long range HFX
442  IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
443  IF (tddfpt_control%do_hfxsr) THEN
444  re_int = tddfpt_control%hfxsr_re_int
445  ! symmetric dmat
446  CALL tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, &
447  kernel_env%full_kernel%admm_env, &
448  kernel_env%full_kernel%hfxsr_section, &
449  kernel_env%full_kernel%x_data, 1, re_int, &
450  work_rho_ia_ao=work_matrices%hfxsr_rho_ao_symm, &
451  work_hmat=work_matrices%hfxsr_hmat_symm, &
452  wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
453  ! antisymmetric dmat
454  CALL tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, &
455  kernel_env%full_kernel%admm_env, &
456  kernel_env%full_kernel%hfxsr_section, &
457  kernel_env%full_kernel%x_data, -1, .false., &
458  work_rho_ia_ao=work_matrices%hfxsr_rho_ao_asymm, &
459  work_hmat=work_matrices%hfxsr_hmat_asymm, &
460  wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
461  tddfpt_control%hfxsr_re_int = .false.
462  END IF
463  IF (tddfpt_control%do_hfxlr) THEN
464  rcut = tddfpt_control%hfxlr_rcut
465  scale = tddfpt_control%hfxlr_scale
466  DO ivect = 1, nvects
467  IF (ALLOCATED(work_matrices%evects_sub)) THEN
468  IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct)) THEN
469  CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
470  work_matrices%evects_sub(:, ivect), &
471  work_matrices%Aop_evects_sub(:, ivect))
472  ELSE
473  ! skip trial vectors which are assigned to different parallel groups
474  cycle
475  END IF
476  ELSE
477  CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
478  evects(:, ivect), aop_evects(:, ivect))
479  END IF
480  END DO
481  END IF
482  END IF
483 
484  END IF
485 
486  CALL timestop(handle)
487 
488  END SUBROUTINE tddfpt_compute_aop_evects
489 
490 ! **************************************************************************************************
491 !> \brief Solve eigenproblem for the reduced action matrix and find new Ritz eigenvectors and
492 !> eigenvalues.
493 !> \param ritz_vects Ritz eigenvectors (initialised on exit)
494 !> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors (initialised on exit)
495 !> \param evals Ritz eigenvalues (initialised on exit)
496 !> \param krylov_vects Krylov's vectors
497 !> \param Aop_krylov action of TDDFPT operator on Krylov's vectors
498 !> \param Atilde ...
499 !> \param nkvo ...
500 !> \param nkvn ...
501 !> \par History
502 !> * 06.2016 created [Sergey Chulkov]
503 !> * 03.2017 altered prototype, OpenMP parallelisation [Sergey Chulkov]
504 ! **************************************************************************************************
505  SUBROUTINE tddfpt_compute_ritz_vects(ritz_vects, Aop_ritz, evals, krylov_vects, Aop_krylov, &
506  Atilde, nkvo, nkvn)
507  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: ritz_vects, aop_ritz
508  REAL(kind=dp), DIMENSION(:), INTENT(out) :: evals
509  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: krylov_vects, aop_krylov
510  REAL(kind=dp), DIMENSION(:, :), POINTER :: atilde
511  INTEGER, INTENT(IN) :: nkvo, nkvn
512 
513  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_compute_ritz_vects'
514 
515  INTEGER :: handle, ikv, irv, ispin, nkvs, nrvs, &
516  nspins
517  REAL(kind=dp) :: act
518  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: at12, at21, at22, evects_atilde
519  TYPE(cp_blacs_env_type), POINTER :: blacs_env_global
520  TYPE(cp_fm_struct_type), POINTER :: fm_struct
521  TYPE(cp_fm_type) :: atilde_fm, evects_atilde_fm
522 
523  CALL timeset(routinen, handle)
524 
525  nspins = SIZE(krylov_vects, 1)
526  nkvs = SIZE(krylov_vects, 2)
527  nrvs = SIZE(ritz_vects, 2)
528  cpassert(nkvs == nkvo + nkvn)
529 
530  CALL cp_fm_get_info(krylov_vects(1, 1), context=blacs_env_global)
531 
532  CALL cp_fm_struct_create(fm_struct, nrow_global=nkvs, ncol_global=nkvs, context=blacs_env_global)
533  CALL cp_fm_create(atilde_fm, fm_struct)
534  CALL cp_fm_create(evects_atilde_fm, fm_struct)
535  CALL cp_fm_struct_release(fm_struct)
536 
537  ! *** compute upper-diagonal reduced action matrix ***
538  CALL reallocate(atilde, 1, nkvs, 1, nkvs)
539  ! TO DO: the subroutine 'cp_fm_contracted_trace' will compute all elements of
540  ! the matrix 'Atilde', however only upper-triangular elements are actually needed
541  !
542  IF (nkvo == 0) THEN
543  CALL cp_fm_contracted_trace(aop_krylov(:, 1:nkvs), krylov_vects(:, 1:nkvs), &
544  atilde(1:nkvs, 1:nkvs), accurate=.false.)
545  ELSE
546  ALLOCATE (at12(nkvn, nkvo), at21(nkvo, nkvn), at22(nkvn, nkvn))
547  CALL cp_fm_contracted_trace(aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, 1:nkvo), &
548  at12, accurate=.false.)
549  atilde(nkvo + 1:nkvs, 1:nkvo) = at12(1:nkvn, 1:nkvo)
550  CALL cp_fm_contracted_trace(aop_krylov(:, 1:nkvo), krylov_vects(:, nkvo + 1:nkvs), &
551  at21, accurate=.false.)
552  atilde(1:nkvo, nkvo + 1:nkvs) = at21(1:nkvo, 1:nkvn)
553  CALL cp_fm_contracted_trace(aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, nkvo + 1:nkvs), &
554  at22, accurate=.false.)
555  atilde(nkvo + 1:nkvs, nkvo + 1:nkvs) = at22(1:nkvn, 1:nkvn)
556  DEALLOCATE (at12, at21, at22)
557  END IF
558  atilde = 0.5_dp*(atilde + transpose(atilde))
559  CALL cp_fm_set_submatrix(atilde_fm, atilde)
560 
561  ! *** solve an eigenproblem for the reduced matrix ***
562  CALL choose_eigv_solver(atilde_fm, evects_atilde_fm, evals(1:nkvs))
563 
564  ALLOCATE (evects_atilde(nkvs, nrvs))
565  CALL cp_fm_get_submatrix(evects_atilde_fm, evects_atilde, start_row=1, start_col=1, n_rows=nkvs, n_cols=nrvs)
566  CALL cp_fm_release(evects_atilde_fm)
567  CALL cp_fm_release(atilde_fm)
568 
569 !$OMP PARALLEL DO DEFAULT(NONE), &
570 !$OMP PRIVATE(act, ikv, irv, ispin), &
571 !$OMP SHARED(Aop_krylov, Aop_ritz, krylov_vects, evects_Atilde, nkvs, nrvs, nspins, ritz_vects)
572  DO irv = 1, nrvs
573  DO ispin = 1, nspins
574  CALL cp_fm_set_all(ritz_vects(ispin, irv), 0.0_dp)
575  CALL cp_fm_set_all(aop_ritz(ispin, irv), 0.0_dp)
576  END DO
577 
578  DO ikv = 1, nkvs
579  act = evects_atilde(ikv, irv)
580  DO ispin = 1, nspins
581  CALL cp_fm_scale_and_add(1.0_dp, ritz_vects(ispin, irv), &
582  act, krylov_vects(ispin, ikv))
583  CALL cp_fm_scale_and_add(1.0_dp, aop_ritz(ispin, irv), &
584  act, aop_krylov(ispin, ikv))
585  END DO
586  END DO
587  END DO
588 !$OMP END PARALLEL DO
589 
590  DEALLOCATE (evects_atilde)
591 
592  CALL timestop(handle)
593 
594  END SUBROUTINE tddfpt_compute_ritz_vects
595 
596 ! **************************************************************************************************
597 !> \brief Expand Krylov space by computing residual vectors.
598 !> \param residual_vects residual vectors (modified on exit)
599 !> \param evals Ritz eigenvalues (modified on exit)
600 !> \param ritz_vects Ritz eigenvectors
601 !> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors
602 !> \param gs_mos molecular orbitals optimised for the ground state
603 !> \param matrix_s overlap matrix
604 !> \par History
605 !> * 06.2016 created [Sergey Chulkov]
606 !> * 03.2017 refactored to achieve significant performance gain [Sergey Chulkov]
607 ! **************************************************************************************************
608  SUBROUTINE tddfpt_compute_residual_vects(residual_vects, evals, ritz_vects, Aop_ritz, gs_mos, &
609  matrix_s)
610  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: residual_vects
611  REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
612  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: ritz_vects, aop_ritz
613  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
614  INTENT(in) :: gs_mos
615  TYPE(dbcsr_type), POINTER :: matrix_s
616 
617  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_compute_residual_vects'
618  REAL(kind=dp), PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*epsilon(1.0_dp)
619 
620  INTEGER :: handle, icol_local, irow_local, irv, &
621  ispin, nao, ncols_local, nrows_local, &
622  nrvs, nspins
623  INTEGER, DIMENSION(:), POINTER :: col_indices_local, row_indices_local
624  INTEGER, DIMENSION(maxspins) :: nactive, nmo_virt
625  REAL(kind=dp) :: e_occ_plus_lambda, eref, lambda
626  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
627  POINTER :: weights_ldata
628  TYPE(cp_fm_struct_type), POINTER :: ao_mo_struct, virt_mo_struct
629  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: awork, vomat
630 
631  CALL timeset(routinen, handle)
632 
633  nspins = SIZE(residual_vects, 1)
634  nrvs = SIZE(residual_vects, 2)
635 
636  IF (nrvs > 0) THEN
637  CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
638  ALLOCATE (awork(nspins), vomat(nspins))
639  DO ispin = 1, nspins
640  nmo_virt(ispin) = SIZE(gs_mos(ispin)%evals_virt)
641  !
642  CALL cp_fm_get_info(matrix=ritz_vects(ispin, 1), matrix_struct=ao_mo_struct, &
643  ncol_global=nactive(ispin))
644  CALL cp_fm_create(awork(ispin), ao_mo_struct)
645  !
646  CALL cp_fm_struct_create(virt_mo_struct, nrow_global=nmo_virt(ispin), &
647  ncol_global=nactive(ispin), template_fmstruct=ao_mo_struct)
648  CALL cp_fm_create(vomat(ispin), virt_mo_struct)
649  CALL cp_fm_struct_release(virt_mo_struct)
650  END DO
651 
652  ! *** actually compute residual vectors ***
653  DO irv = 1, nrvs
654  lambda = evals(irv)
655 
656  DO ispin = 1, nspins
657  CALL cp_fm_get_info(vomat(ispin), nrow_local=nrows_local, &
658  ncol_local=ncols_local, row_indices=row_indices_local, &
659  col_indices=col_indices_local, local_data=weights_ldata)
660 
661  ! awork := Ab(ispin, irv) - evals(irv) b(ispin, irv), where 'b' is a Ritz vector
662  CALL cp_dbcsr_sm_fm_multiply(matrix_s, ritz_vects(ispin, irv), awork(ispin), &
663  ncol=nactive(ispin), alpha=-lambda, beta=0.0_dp)
664  CALL cp_fm_scale_and_add(1.0_dp, awork(ispin), 1.0_dp, aop_ritz(ispin, irv))
665  !
666  CALL parallel_gemm('T', 'N', nmo_virt(ispin), nactive(ispin), nao, 1.0_dp, gs_mos(ispin)%mos_virt, &
667  awork(ispin), 0.0_dp, vomat(ispin))
668 
669  DO icol_local = 1, ncols_local
670  e_occ_plus_lambda = gs_mos(ispin)%evals_occ(col_indices_local(icol_local)) + lambda
671 
672  DO irow_local = 1, nrows_local
673  eref = gs_mos(ispin)%evals_virt(row_indices_local(irow_local)) - e_occ_plus_lambda
674 
675  ! eref = e_virt - e_occ - lambda = e_virt - e_occ - (eref_scale*lambda + (1-eref_scale)*lambda);
676  ! eref_new = e_virt - e_occ - eref_scale*lambda = eref + (1 - eref_scale)*lambda
677  IF (abs(eref) < threshold) &
678  eref = eref + (1.0_dp - eref_scale)*lambda
679 
680  weights_ldata(irow_local, icol_local) = weights_ldata(irow_local, icol_local)/eref
681  END DO
682  END DO
683 
684  CALL parallel_gemm('N', 'N', nao, nactive(ispin), nmo_virt(ispin), 1.0_dp, gs_mos(ispin)%mos_virt, &
685  vomat(ispin), 0.0_dp, residual_vects(ispin, irv))
686  END DO
687  END DO
688  !
689  CALL cp_fm_release(awork)
690  CALL cp_fm_release(vomat)
691  END IF
692 
693  CALL timestop(handle)
694 
695  END SUBROUTINE tddfpt_compute_residual_vects
696 
697 ! **************************************************************************************************
698 !> \brief Perform Davidson iterations.
699 !> \param evects TDDFPT trial vectors (modified on exit)
700 !> \param evals TDDFPT eigenvalues (modified on exit)
701 !> \param S_evects cached matrix product S * evects (modified on exit)
702 !> \param gs_mos molecular orbitals optimised for the ground state
703 !> \param tddfpt_control TDDFPT control parameters
704 !> \param matrix_ks Kohn-Sham matrix
705 !> \param qs_env Quickstep environment
706 !> \param kernel_env kernel environment
707 !> \param sub_env parallel (sub)group environment
708 !> \param logger CP2K logger
709 !> \param iter_unit I/O unit to write basic iteration information
710 !> \param energy_unit I/O unit to write detailed energy information
711 !> \param tddfpt_print_section TDDFPT print input section (need to write TDDFPT restart files)
712 !> \param work_matrices collection of work matrices (modified on exit)
713 !> \return energy convergence achieved (in Hartree)
714 !> \par History
715 !> * 03.2017 code related to Davidson eigensolver has been moved here from the main subroutine
716 !> tddfpt() [Sergey Chulkov]
717 !> \note Based on the subroutines apply_op() and iterative_solver() originally created by
718 !> Thomas Chassaing in 2002.
719 ! **************************************************************************************************
720  FUNCTION tddfpt_davidson_solver(evects, evals, S_evects, gs_mos, tddfpt_control, &
721  matrix_ks, qs_env, kernel_env, &
722  sub_env, logger, iter_unit, energy_unit, &
723  tddfpt_print_section, work_matrices) RESULT(conv)
724  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout) :: evects
725  REAL(kind=dp), DIMENSION(:), INTENT(inout) :: evals
726  TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout) :: s_evects
727  TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
728  INTENT(in) :: gs_mos
729  TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
730  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
731  TYPE(qs_environment_type), POINTER :: qs_env
732  TYPE(kernel_env_type), INTENT(in) :: kernel_env
733  TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
734  TYPE(cp_logger_type), POINTER :: logger
735  INTEGER, INTENT(in) :: iter_unit, energy_unit
736  TYPE(section_vals_type), POINTER :: tddfpt_print_section
737  TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices
738  REAL(kind=dp) :: conv
739 
740  CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_davidson_solver'
741 
742  INTEGER :: handle, ispin, istate, iter, &
743  max_krylov_vects, nspins, nstates, &
744  nstates_conv, nvects_exist, nvects_new
745  INTEGER(kind=int_8) :: nstates_total
746  LOGICAL :: is_nonortho
747  REAL(kind=dp) :: t1, t2
748  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_last
749  REAL(kind=dp), DIMENSION(:, :), POINTER :: atilde
750  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: aop_krylov, aop_ritz, krylov_vects, &
751  s_krylov
752  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
753 
754  CALL timeset(routinen, handle)
755 
756  nspins = SIZE(gs_mos)
757  nstates = tddfpt_control%nstates
758  nstates_total = tddfpt_total_number_of_states(gs_mos)
759 
760  IF (debug_this_module) THEN
761  cpassert(SIZE(evects, 1) == nspins)
762  cpassert(SIZE(evects, 2) == nstates)
763  cpassert(SIZE(evals) == nstates)
764  END IF
765 
766  CALL get_qs_env(qs_env, matrix_s=matrix_s)
767 
768  ! adjust the number of Krylov vectors
769  max_krylov_vects = tddfpt_control%nkvs
770  IF (max_krylov_vects < nstates) max_krylov_vects = nstates
771  IF (int(max_krylov_vects, kind=int_8) > nstates_total) max_krylov_vects = int(nstates_total)
772 
773  ALLOCATE (aop_ritz(nspins, nstates))
774  DO istate = 1, nstates
775  DO ispin = 1, nspins
776  CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, aop_ritz(ispin, istate))
777  END DO
778  END DO
779 
780  ALLOCATE (evals_last(max_krylov_vects))
781  ALLOCATE (aop_krylov(nspins, max_krylov_vects), krylov_vects(nspins, max_krylov_vects), &
782  s_krylov(nspins, max_krylov_vects))
783 
784  DO istate = 1, nstates
785  DO ispin = 1, nspins
786  CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, krylov_vects(ispin, istate))
787  CALL cp_fm_to_fm(evects(ispin, istate), krylov_vects(ispin, istate))
788 
789  CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, s_krylov(ispin, istate))
790  CALL cp_fm_to_fm(s_evects(ispin, istate), s_krylov(ispin, istate))
791 
792  CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, aop_krylov(ispin, istate))
793  END DO
794  END DO
795 
796  nvects_exist = 0
797  nvects_new = nstates
798 
799  t1 = m_walltime()
800 
801  ALLOCATE (atilde(1, 1))
802 
803  DO
804  ! davidson iteration
805  CALL cp_iterate(logger%iter_info, iter_nr_out=iter)
806 
807  CALL tddfpt_compute_aop_evects(aop_evects=aop_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), &
808  evects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
809  s_evects=s_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), &
810  gs_mos=gs_mos, tddfpt_control=tddfpt_control, &
811  matrix_ks=matrix_ks, &
812  qs_env=qs_env, kernel_env=kernel_env, &
813  sub_env=sub_env, &
814  work_matrices=work_matrices)
815 
816  CALL tddfpt_compute_ritz_vects(ritz_vects=evects, aop_ritz=aop_ritz, &
817  evals=evals_last(1:nvects_exist + nvects_new), &
818  krylov_vects=krylov_vects(:, 1:nvects_exist + nvects_new), &
819  aop_krylov=aop_krylov(:, 1:nvects_exist + nvects_new), &
820  atilde=atilde, nkvo=nvects_exist, nkvn=nvects_new)
821 
822  CALL tddfpt_write_restart(evects=evects, evals=evals_last(1:nstates), gs_mos=gs_mos, &
823  logger=logger, tddfpt_print_section=tddfpt_print_section)
824 
825  conv = maxval(abs(evals_last(1:nstates) - evals(1:nstates)))
826 
827  nvects_exist = nvects_exist + nvects_new
828  IF (nvects_exist + nvects_new > max_krylov_vects) &
829  nvects_new = max_krylov_vects - nvects_exist
830  IF (iter >= tddfpt_control%niters) nvects_new = 0
831 
832  IF (conv > tddfpt_control%conv .AND. nvects_new > 0) THEN
833  ! compute residual vectors for the next iteration
834  DO istate = 1, nvects_new
835  DO ispin = 1, nspins
836  CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, &
837  krylov_vects(ispin, nvects_exist + istate))
838  CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, &
839  s_krylov(ispin, nvects_exist + istate))
840  CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, &
841  aop_krylov(ispin, nvects_exist + istate))
842  END DO
843  END DO
844 
845  CALL tddfpt_compute_residual_vects(residual_vects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
846  evals=evals_last(1:nvects_new), &
847  ritz_vects=evects(:, 1:nvects_new), aop_ritz=aop_ritz(:, 1:nvects_new), &
848  gs_mos=gs_mos, matrix_s=matrix_s(1)%matrix)
849 
850  CALL tddfpt_orthogonalize_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
851  work_matrices%S_C0_C0T)
852 
853  CALL tddfpt_orthonormalize_psi1_psi1(krylov_vects(:, 1:nvects_exist + nvects_new), nvects_new, &
854  s_krylov(:, 1:nvects_exist + nvects_new), matrix_s(1)%matrix)
855 
856  is_nonortho = tddfpt_is_nonorthogonal_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
857  work_matrices%S_C0, tddfpt_control%orthogonal_eps)
858  ELSE
859  ! convergence or the maximum number of Krylov vectors have been achieved
860  nvects_new = 0
861  is_nonortho = .false.
862  END IF
863 
864  t2 = m_walltime()
865  IF (energy_unit > 0) THEN
866  WRITE (energy_unit, '(/,4X,A,T14,A,T36,A)') "State", "Exc. energy (eV)", "Convergence (eV)"
867  DO istate = 1, nstates
868  WRITE (energy_unit, '(1X,I8,T12,F14.7,T38,ES11.4)') istate, &
869  evals_last(istate)*evolt, (evals_last(istate) - evals(istate))*evolt
870  END DO
871  WRITE (energy_unit, *)
872  CALL m_flush(energy_unit)
873  END IF
874 
875  IF (iter_unit > 0) THEN
876  nstates_conv = 0
877  DO istate = 1, nstates
878  IF (abs(evals_last(istate) - evals(istate)) <= tddfpt_control%conv) &
879  nstates_conv = nstates_conv + 1
880  END DO
881 
882  WRITE (iter_unit, '(T7,I8,T24,F7.1,T40,ES11.4,T66,I8)') iter, t2 - t1, conv, nstates_conv
883  CALL m_flush(iter_unit)
884  END IF
885 
886  t1 = t2
887  evals(1:nstates) = evals_last(1:nstates)
888 
889  ! nvects_new == 0 if iter >= tddfpt_control%niters
890  IF (nvects_new == 0 .OR. is_nonortho) THEN
891  ! restart Davidson iterations
892  CALL tddfpt_orthogonalize_psi1_psi0(evects, work_matrices%S_C0_C0T)
893  CALL tddfpt_orthonormalize_psi1_psi1(evects, nstates, s_evects, matrix_s(1)%matrix)
894 
895  EXIT
896  END IF
897  END DO
898 
899  DEALLOCATE (atilde)
900 
901  DO istate = nvects_exist + nvects_new, 1, -1
902  DO ispin = nspins, 1, -1
903  CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, aop_krylov(ispin, istate))
904  CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, s_krylov(ispin, istate))
905  CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, krylov_vects(ispin, istate))
906  END DO
907  END DO
908  DEALLOCATE (aop_krylov, krylov_vects, s_krylov)
909  DEALLOCATE (evals_last)
910 
911  DO istate = nstates, 1, -1
912  DO ispin = nspins, 1, -1
913  CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_occ(ispin)%pool, aop_ritz(ispin, istate))
914  END DO
915  END DO
916  DEALLOCATE (aop_ritz)
917 
918  CALL timestop(handle)
919 
920  END FUNCTION tddfpt_davidson_solver
921 
922 END MODULE qs_tddfpt2_eigensolver
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
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
basic linear algebra operations for full matrices
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition: cp_fm_diag.F:17
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition: cp_fm_diag.F:208
pool for for elements that are retained and released
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
subroutine, public fm_pool_give_back_fm(pool, element)
returns the element to the pool
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_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition: cp_fm_types.F:768
subroutine, public cp_fm_maxabsval(matrix, a_max, ir_max, ic_max)
find the maximum absolute value of the matrix element maxval(abs(matrix))
Definition: cp_fm_types.F:1064
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_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition: cp_fm_types.F:901
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 ...
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
subroutine, public cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
adds one to the actual iteration
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public tddfpt_kernel_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
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public int_8
Definition: kinds.F:54
integer, parameter, public dp
Definition: kinds.F:34
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition: machine.F:123
Utility routines for the memory handling.
Interface to the message passing library MPI.
basic linear algebra operations for full matrixes
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public evolt
Definition: physcon.F:183
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.
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
subroutine, public eigensolver(matrix_ks_fm, mo_set, ortho, work, cholesky_method, do_level_shift, level_shift, matrix_u_fm, use_jacobi)
Diagonalise the Kohn-Sham matrix to get a new set of MO eigen- vectors and MO eigenvalues....
real(kind=dp) function, public tddfpt_davidson_solver(evects, evals, S_evects, gs_mos, tddfpt_control, matrix_ks, qs_env, kernel_env, sub_env, logger, iter_unit, energy_unit, tddfpt_print_section, work_matrices)
Perform Davidson iterations.
subroutine, public tddfpt_orthogonalize_psi1_psi0(evects, S_C0_C0T)
Make TDDFPT trial vectors orthogonal to all occupied molecular orbitals.
subroutine, public tddfpt_orthonormalize_psi1_psi1(evects, nvects_new, S_evects, matrix_s)
Make new TDDFPT trial vectors orthonormal to all previous trial vectors.
subroutine, public stda_kernel(Aop_evects, evects, is_rks_triplets, qs_env, stda_control, stda_env, sub_env, work_matrices)
Compute action matrix-vector products with the sTDA Kernel.
subroutine, public fhxc_kernel(Aop_evects, evects, is_rks_triplets, do_hfx, do_admm, qs_env, kernel_env, kernel_env_admm_aux, sub_env, work_matrices, admm_symm, admm_xc_correction, do_lrigpw)
Compute action matrix-vector products with the FHxc Kernel.
subroutine, public tddfpt_apply_energy_diff(Aop_evects, evects, S_evects, gs_mos, matrix_ks)
Apply orbital energy difference term: Aop_evects(spin,state) += KS(spin) * evects(spin,...
subroutine, public tddfpt_apply_hfx(Aop_evects, evects, gs_mos, do_admm, qs_env, work_rho_ia_ao_symm, work_hmat_symm, work_rho_ia_ao_asymm, work_hmat_asymm, wfm_rho_orb)
Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
subroutine, public tddfpt_apply_hfxsr_kernel(Aop_evects, evects, gs_mos, qs_env, admm_env, hfx_section, x_data, symmetry, recalc_integrals, work_rho_ia_ao, work_hmat, wfm_rho_orb)
Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
subroutine, public tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, hfx_scale, work, X, res)
...Calculate the HFXLR kernel contribution by contracting the Lowdin MO coefficients – transition cha...
subroutine, public tddfpt_write_restart(evects, evals, gs_mos, logger, tddfpt_print_section)
Write Ritz vectors to a binary restart file.
pure integer(kind=int_8) function, public tddfpt_total_number_of_states(gs_mos)
Compute the number of possible singly excited states (occ -> virt)