(git:374b731)
Loading...
Searching...
No Matches
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
8MODULE qs_tddfpt2_eigensolver
9 USE cp_blacs_env, ONLY: cp_blacs_env_type
10 USE cp_control_types, ONLY: tddfpt2_control_type
11 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
12 USE cp_fm_basic_linalg, ONLY: cp_fm_contracted_trace,&
13 cp_fm_scale,&
14 cp_fm_scale_and_add,&
15 cp_fm_trace
16 USE cp_fm_diag, ONLY: choose_eigv_solver
17 USE cp_fm_pool_types, ONLY: fm_pool_create_fm,&
18 fm_pool_give_back_fm
19 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
20 cp_fm_struct_release,&
21 cp_fm_struct_type
22 USE cp_fm_types, ONLY: &
23 cp_fm_copy_general, cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_maxabsval, &
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
26 USE cp_output_handling, ONLY: cp_iterate
27 USE dbcsr_api, ONLY: dbcsr_get_info,&
28 dbcsr_p_type,&
29 dbcsr_type
30 USE input_constants, ONLY: tddfpt_kernel_full,&
31 tddfpt_kernel_none,&
32 tddfpt_kernel_stda
33 USE input_section_types, ONLY: section_vals_type
34 USE kinds, ONLY: dp,&
35 int_8
36 USE machine, ONLY: m_flush,&
37 m_walltime
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,&
48 stda_kernel
49 USE qs_tddfpt2_operators, ONLY: tddfpt_apply_energy_diff,&
50 tddfpt_apply_hfx,&
51 tddfpt_apply_hfxlr_kernel,&
52 tddfpt_apply_hfxsr_kernel
53 USE qs_tddfpt2_restart, ONLY: tddfpt_write_restart
54 USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
55 USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos,&
56 tddfpt_work_matrices
57 USE qs_tddfpt2_utils, ONLY: tddfpt_total_number_of_states
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
71 PUBLIC :: tddfpt_davidson_solver, tddfpt_orthogonalize_psi1_psi0, &
72 tddfpt_orthonormalize_psi1_psi1
73
74! **************************************************************************************************
75
76CONTAINS
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
922END MODULE qs_tddfpt2_eigensolver
cp_fm_basic_linalg::cp_fm_contracted_trace
Definition cp_fm_basic_linalg.F:84
cp_fm_basic_linalg::cp_fm_trace
Definition cp_fm_basic_linalg.F:74
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
memory_utilities::reallocate
Definition memory_utilities.F:27
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:744
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
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....
Definition cp_fm_basic_linalg.F:167
cp_fm_basic_linalg::cp_fm_scale
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
Definition cp_fm_basic_linalg.F:1680
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::choose_eigv_solver
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
cp_fm_pool_types
pool for for elements that are retained and released
Definition cp_fm_pool_types.F:14
cp_fm_pool_types::fm_pool_create_fm
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
Definition cp_fm_pool_types.F:185
cp_fm_pool_types::fm_pool_give_back_fm
subroutine, public fm_pool_give_back_fm(pool, element)
returns the element to the pool
Definition cp_fm_pool_types.F:224
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_copy_general
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
cp_fm_types::cp_fm_get_info
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
cp_fm_types::cp_fm_set_submatrix
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
cp_fm_types::cp_fm_maxabsval
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
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:535
cp_fm_types::cp_fm_get_submatrix
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
cp_fm_types::cp_fm_create
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
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_iterate
subroutine, public cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
adds one to the actual iteration
Definition cp_output_handling.F:574
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::tddfpt_kernel_none
integer, parameter, public tddfpt_kernel_none
Definition input_constants.F:581
input_constants::tddfpt_kernel_full
integer, parameter, public tddfpt_kernel_full
Definition input_constants.F:581
input_constants::tddfpt_kernel_stda
integer, parameter, public tddfpt_kernel_stda
Definition input_constants.F:581
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:106
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:123
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition physcon.F:183
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, 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.
Definition qs_environment_types.F:502
qs_kernel_types
Definition qs_kernel_types.F:8
qs_scf_methods
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
Definition qs_scf_methods.F:19
qs_scf_methods::eigensolver
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....
Definition qs_scf_methods.F:168
qs_tddfpt2_eigensolver
Definition qs_tddfpt2_eigensolver.F:8
qs_tddfpt2_eigensolver::tddfpt_orthonormalize_psi1_psi1
subroutine, public tddfpt_orthonormalize_psi1_psi1(evects, nvects_new, s_evects, matrix_s)
Make new TDDFPT trial vectors orthonormal to all previous trial vectors.
Definition qs_tddfpt2_eigensolver.F:225
qs_tddfpt2_eigensolver::tddfpt_orthogonalize_psi1_psi0
subroutine, public tddfpt_orthogonalize_psi1_psi0(evects, s_c0_c0t)
Make TDDFPT trial vectors orthogonal to all occupied molecular orbitals.
Definition qs_tddfpt2_eigensolver.F:95
qs_tddfpt2_eigensolver::tddfpt_davidson_solver
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.
Definition qs_tddfpt2_eigensolver.F:724
qs_tddfpt2_fhxc
Definition qs_tddfpt2_fhxc.F:8
qs_tddfpt2_fhxc::fhxc_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.
Definition qs_tddfpt2_fhxc.F:98
qs_tddfpt2_fhxc::stda_kernel
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.
Definition qs_tddfpt2_fhxc.F:492
qs_tddfpt2_operators
Definition qs_tddfpt2_operators.F:8
qs_tddfpt2_operators::tddfpt_apply_hfx
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.
Definition qs_tddfpt2_operators.F:497
qs_tddfpt2_operators::tddfpt_apply_hfxlr_kernel
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...
Definition qs_tddfpt2_operators.F:745
qs_tddfpt2_operators::tddfpt_apply_energy_diff
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,...
Definition qs_tddfpt2_operators.F:102
qs_tddfpt2_operators::tddfpt_apply_hfxsr_kernel
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.
Definition qs_tddfpt2_operators.F:663
qs_tddfpt2_restart
Definition qs_tddfpt2_restart.F:8
qs_tddfpt2_restart::tddfpt_write_restart
subroutine, public tddfpt_write_restart(evects, evals, gs_mos, logger, tddfpt_print_section)
Write Ritz vectors to a binary restart file.
Definition qs_tddfpt2_restart.F:77
qs_tddfpt2_subgroups
Definition qs_tddfpt2_subgroups.F:8
qs_tddfpt2_types
Definition qs_tddfpt2_types.F:8
qs_tddfpt2_utils
Definition qs_tddfpt2_utils.F:8
qs_tddfpt2_utils::tddfpt_total_number_of_states
pure integer(kind=int_8) function, public tddfpt_total_number_of_states(gs_mos)
Compute the number of possible singly excited states (occ -> virt)
Definition qs_tddfpt2_utils.F:735
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::tddfpt2_control_type
Definition cp_control_types.F:465
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_kernel_types::full_kernel_env_type
Collection of variables required to evaluate adiabatic TDDFPT kernel.
Definition qs_kernel_types.F:46
qs_kernel_types::kernel_env_type
Type to hold environments for the different kernels.
Definition qs_kernel_types.F:83
qs_tddfpt2_subgroups::tddfpt_subgroup_env_type
Parallel (sub)group environment.
Definition qs_tddfpt2_subgroups.F:102
qs_tddfpt2_types::tddfpt_ground_state_mos
Ground state molecular orbitals.
Definition qs_tddfpt2_types.F:86
qs_tddfpt2_types::tddfpt_work_matrices
Set of temporary ("work") matrices.
Definition qs_tddfpt2_types.F:110