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