(git:7a2a232)
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-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
11 USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
25 USE cp_fm_types, ONLY: &
30 USE input_constants, ONLY: no_sf_tddfpt,&
35 USE kinds, ONLY: dp,&
36 int_8
37 USE machine, ONLY: m_flush,&
42 USE physcon, ONLY: evolt
50 USE qs_tddfpt2_fhxc, ONLY: fhxc_kernel,&
64#include "./base/base_uses.f90"
65
66 IMPLICIT NONE
67
68 PRIVATE
69
70 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_eigensolver'
71
72 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
73 ! number of first derivative components (3: d/dx, d/dy, d/dz)
74 INTEGER, PARAMETER, PRIVATE :: nderivs = 3
75 INTEGER, PARAMETER, PRIVATE :: maxspins = 2
76
79
80! **************************************************************************************************
81
82CONTAINS
83
84! **************************************************************************************************
85!> \brief Make TDDFPT trial vectors orthogonal to all occupied molecular orbitals.
86!> \param evects trial vectors distributed across all processors (modified on exit)
87!> \param S_C0_C0T matrix product S * C_0 * C_0^T, where C_0 is the ground state
88!> wave function for each spin expressed in atomic basis set,
89!> and S is the corresponding overlap matrix
90!> \param qs_env ...
91!> \param gs_mos ...
92!> \param evals ...
93!> \param tddfpt_control ...
94!> \param S_C0 ...
95!> \par History
96!> * 05.2016 created [Sergey Chulkov]
97!> * 05.2019 use a temporary work matrix [JHU]
98!> \note Based on the subroutine p_preortho() which was created by Thomas Chassaing on 09.2002.
99!> Should be useless when ground state MOs are computed with extremely high accuracy,
100!> as all virtual orbitals are already orthogonal to the occupied ones by design.
101!> However, when the norm of residual vectors is relatively small (e.g. less then SCF_EPS),
102!> new Krylov's vectors seem to be random and should be orthogonalised even with respect to
103!> the occupied MOs.
104! **************************************************************************************************
105 SUBROUTINE tddfpt_orthogonalize_psi1_psi0(evects, S_C0_C0T, qs_env, gs_mos, evals, &
106 tddfpt_control, S_C0)
107 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
108 TYPE(cp_fm_type), DIMENSION(:), INTENT(in) :: s_c0_c0t
109 TYPE(qs_environment_type), POINTER :: qs_env
110 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
111 INTENT(in) :: gs_mos
112 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
113 TYPE(tddfpt2_control_type), INTENT(in), POINTER :: tddfpt_control
114 TYPE(cp_fm_type), DIMENSION(:), INTENT(in) :: s_c0
115
116 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_orthogonalize_psi1_psi0'
117
118 INTEGER :: handle, ispin, ivect, nactive, nao, &
119 nspins, nvects, spin
120 TYPE(cp_fm_struct_type), POINTER :: matrix_struct
121 TYPE(cp_fm_type) :: evortho
122 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos
123 TYPE(mp_para_env_type), POINTER :: para_env
124
125 CALL timeset(routinen, handle)
126
127 nspins = SIZE(evects, 1)
128 nvects = SIZE(evects, 2)
129
130 IF (nvects > 0) THEN
131 IF (.NOT. tddfpt_control%do_smearing) THEN
132 DO ispin = 1, nspins
133 IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
134 spin = ispin
135 ELSE
136 spin = 2
137 END IF
138 CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
139 nrow_global=nao, ncol_global=nactive)
140 CALL cp_fm_create(evortho, matrix_struct)
141 DO ivect = 1, nvects
142 ! evortho: C0 * C0^T * S * C1 == (S * C0 * C0^T)^T * C1
143 CALL parallel_gemm('T', 'N', nao, nactive, nao, 1.0_dp, s_c0_c0t(spin), &
144 evects(ispin, ivect), 0.0_dp, evortho)
145 CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, ivect), -1.0_dp, evortho)
146 END DO
147 CALL cp_fm_release(evortho)
148 END DO
149 ELSE
150 NULLIFY (para_env)
151 CALL get_qs_env(qs_env, para_env=para_env)
152 NULLIFY (mos)
153 ALLOCATE (mos(nspins))
154 DO ispin = 1, nspins
155 CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
156 nrow_global=nao, ncol_global=nactive)
157 CALL cp_fm_create(mos(ispin), matrix_struct)
158 CALL cp_fm_copy_general(gs_mos(ispin)%mos_occ, mos(ispin), para_env)
159 END DO
160 DO ivect = 1, nvects
161 CALL compute_fermib(qs_env, gs_mos, evals(ivect))
162 CALL orthogonalize_smeared_occupation(evects(:, ivect), qs_env, mos, s_c0)
163 END DO
164 DO ispin = 1, nspins
165 CALL cp_fm_release(mos(ispin))
166 END DO
167 DEALLOCATE (mos)
168 END IF
169 END IF
170
171 CALL timestop(handle)
172
173 END SUBROUTINE tddfpt_orthogonalize_psi1_psi0
174
175! **************************************************************************************************
176!> \brief Check that orthogonalised TDDFPT trial vectors remain orthogonal to
177!> occupied molecular orbitals.
178!> \param evects trial vectors
179!> \param S_C0 matrix product S * C_0, where C_0 is the ground state wave function
180!> for each spin in atomic basis set, and S is the corresponding overlap matrix
181!> \param max_norm the largest possible overlap between the ground state and
182!> excited state wave functions
183!> \param spinflip ...
184!> \return true if trial vectors are non-orthogonal to occupied molecular orbitals
185!> \par History
186!> * 07.2016 created [Sergey Chulkov]
187!> * 05.2019 use temporary work matrices [JHU]
188! **************************************************************************************************
189 FUNCTION tddfpt_is_nonorthogonal_psi1_psi0(evects, S_C0, max_norm, spinflip) &
190 result(is_nonortho)
191 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
192 TYPE(cp_fm_type), DIMENSION(:), INTENT(in) :: s_c0
193 REAL(kind=dp), INTENT(in) :: max_norm
194 INTEGER, INTENT(in) :: spinflip
195 LOGICAL :: is_nonortho
196
197 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_is_nonorthogonal_psi1_psi0'
198
199 INTEGER :: handle, ispin, ivect, nactive, nao, &
200 nocc, nspins, nvects, spin
201 REAL(kind=dp) :: maxabs_val
202 TYPE(cp_fm_struct_type), POINTER :: matrix_struct, matrix_struct_tmp
203 TYPE(cp_fm_type) :: aortho
204
205 CALL timeset(routinen, handle)
206
207 nspins = SIZE(evects, 1)
208 nvects = SIZE(evects, 2)
209
210 is_nonortho = .false.
211
212 loop: DO ispin = 1, nspins
213
214 IF (spinflip == no_sf_tddfpt) THEN
215 spin = ispin
216 ELSE
217 spin = 2
218 END IF
219
220 CALL cp_fm_get_info(matrix=s_c0(spin), ncol_global=nocc)
221 CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
222 nrow_global=nao, ncol_global=nactive)
223 CALL cp_fm_struct_create(matrix_struct_tmp, nrow_global=nocc, &
224 ncol_global=nactive, template_fmstruct=matrix_struct)
225 CALL cp_fm_create(aortho, matrix_struct_tmp)
226 CALL cp_fm_struct_release(matrix_struct_tmp)
227 DO ivect = 1, nvects
228 ! aortho = S_C0^T * S * C_1
229 CALL parallel_gemm('T', 'N', nocc, nactive, nao, 1.0_dp, s_c0(spin), &
230 evects(ispin, ivect), 0.0_dp, aortho)
231 CALL cp_fm_maxabsval(aortho, maxabs_val)
232 is_nonortho = maxabs_val > max_norm
233 IF (is_nonortho) THEN
234 CALL cp_fm_release(aortho)
235 EXIT loop
236 END IF
237 END DO
238 CALL cp_fm_release(aortho)
239 END DO loop
240
241 CALL timestop(handle)
242
243 END FUNCTION tddfpt_is_nonorthogonal_psi1_psi0
244
245! **************************************************************************************************
246!> \brief Make new TDDFPT trial vectors orthonormal to all previous trial vectors.
247!> \param evects trial vectors (modified on exit)
248!> \param nvects_new number of new trial vectors to orthogonalise
249!> \param S_evects set of matrices to store matrix product S * evects (modified on exit)
250!> \param matrix_s overlap matrix
251!> \par History
252!> * 05.2016 created [Sergey Chulkov]
253!> * 02.2017 caching the matrix product S * evects [Sergey Chulkov]
254!> \note \parblock
255!> Based on the subroutines reorthogonalize() and normalize() which were originally created
256!> by Thomas Chassaing on 03.2003.
257!>
258!> In order to orthogonalise a trial vector C3 = evects(:,3) with respect to previously
259!> orthogonalised vectors C1 = evects(:,1) and C2 = evects(:,2) we need to compute the
260!> quantity C3'' using the following formulae:
261!> C3' = C3 - Tr(C3^T * S * C1) * C1,
262!> C3'' = C3' - Tr(C3'^T * S * C2) * C2,
263!> which can be expanded as:
264!> C3'' = C3 - Tr(C3^T * S * C1) * C1 - Tr(C3^T * S * C2) * C2 +
265!> Tr(C3^T * S * C1) * Tr(C2^T * S * C1) * C2 .
266!> In case of unlimited float-point precision, the last term in above expression is exactly 0,
267!> due to orthogonality condition between C1 and C2. In this case the expression could be
268!> simplified as (taking into account the identity: Tr(A * S * B) = Tr(B * S * A)):
269!> C3'' = C3 - Tr(C1^T * S * C3) * C1 - Tr(C2^T * S * C3) * C2 ,
270!> which means we do not need the variable S_evects to keep the matrix products S * Ci .
271!>
272!> In reality, however, we deal with limited float-point precision arithmetic meaning that
273!> the trace Tr(C2^T * S * C1) is close to 0 but does not equal to 0 exactly. The term
274!> Tr(C3^T * S * C1) * Tr(C2^T * S * C1) * C2
275!> can not be ignored anymore. Ignorance of this term will lead to numerical instability
276!> when the trace Tr(C3^T * S * C1) is large enough.
277!> \endparblock
278! **************************************************************************************************
279 SUBROUTINE tddfpt_orthonormalize_psi1_psi1(evects, nvects_new, S_evects, matrix_s)
280 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: evects
281 INTEGER, INTENT(in) :: nvects_new
282 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: s_evects
283 TYPE(dbcsr_type), POINTER :: matrix_s
284
285 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_orthonormalize_psi1_psi1'
286
287 INTEGER :: handle, ispin, ivect, jvect, nspins, &
288 nvects_old, nvects_total
289 INTEGER, DIMENSION(maxspins) :: nactive
290 REAL(kind=dp) :: norm
291 REAL(kind=dp), DIMENSION(maxspins) :: weights
292
293 CALL timeset(routinen, handle)
294
295 nspins = SIZE(evects, 1)
296 nvects_total = SIZE(evects, 2)
297 nvects_old = nvects_total - nvects_new
298
299 IF (debug_this_module) THEN
300 cpassert(SIZE(s_evects, 1) == nspins)
301 cpassert(SIZE(s_evects, 2) == nvects_total)
302 cpassert(nvects_old >= 0)
303 END IF
304
305 DO ispin = 1, nspins
306 CALL cp_fm_get_info(matrix=evects(ispin, 1), ncol_global=nactive(ispin))
307 END DO
308
309 DO jvect = nvects_old + 1, nvects_total
310 ! Orthogonalization <psi1_i | psi1_j>
311 DO ivect = 1, jvect - 1
312 CALL cp_fm_trace(evects(:, jvect), s_evects(:, ivect), weights(1:nspins), accurate=.false.)
313 norm = sum(weights(1:nspins))
314
315 DO ispin = 1, nspins
316 CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, jvect), -norm, evects(ispin, ivect))
317 END DO
318 END DO
319
320 ! Normalization <psi1_j | psi1_j> = 1
321 DO ispin = 1, nspins
322 CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, jvect), s_evects(ispin, jvect), &
323 ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
324 END DO
325
326 CALL cp_fm_trace(evects(:, jvect), s_evects(:, jvect), weights(1:nspins), accurate=.false.)
327
328 norm = sum(weights(1:nspins))
329 norm = 1.0_dp/sqrt(norm)
330
331 DO ispin = 1, nspins
332 CALL cp_fm_scale(norm, evects(ispin, jvect))
333 CALL cp_fm_scale(norm, s_evects(ispin, jvect))
334 END DO
335 END DO
336
337 CALL timestop(handle)
338
340
341! **************************************************************************************************
342!> \brief Compute action matrix-vector products.
343!> \param Aop_evects action of TDDFPT operator on trial vectors (modified on exit)
344!> \param evects TDDFPT trial vectors
345!> \param S_evects cached matrix product S * evects where S is the overlap matrix
346!> in primary basis set
347!> \param gs_mos molecular orbitals optimised for the ground state
348!> \param tddfpt_control control section for tddfpt
349!> \param matrix_ks Kohn-Sham matrix
350!> \param qs_env Quickstep environment
351!> \param kernel_env kernel environment
352!> \param sub_env parallel (sub)group environment
353!> \param work_matrices collection of work matrices (modified on exit)
354!> \param matrix_s ...
355!> \par History
356!> * 06.2016 created [Sergey Chulkov]
357!> * 03.2017 refactored [Sergey Chulkov]
358! **************************************************************************************************
359 SUBROUTINE tddfpt_compute_aop_evects(Aop_evects, evects, S_evects, gs_mos, tddfpt_control, &
360 matrix_ks, qs_env, kernel_env, &
361 sub_env, work_matrices, matrix_s)
362 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(INOUT) :: aop_evects
363 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: evects, s_evects
364 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
365 INTENT(in) :: gs_mos
366 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
367 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(in) :: matrix_ks
368 TYPE(qs_environment_type), POINTER :: qs_env
369 TYPE(kernel_env_type), INTENT(in) :: kernel_env
370 TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
371 TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices
372 TYPE(dbcsr_type), POINTER :: matrix_s
373
374 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_compute_Aop_evects'
375
376 INTEGER :: handle, ispin, ivect, nspins, nvects
377 INTEGER, DIMENSION(maxspins) :: nmo_occ
378 LOGICAL :: do_admm, do_bse, do_bse_gw_only, &
379 do_bse_w_only, do_hfx, &
380 do_lri_response, is_rks_triplets, &
381 re_int
382 REAL(kind=dp) :: rcut, scale
383 TYPE(cp_fm_type) :: fm_dummy
384 TYPE(full_kernel_env_type), POINTER :: kernel_env_admm_aux
385 TYPE(mp_para_env_type), POINTER :: para_env
386
387 CALL timeset(routinen, handle)
388
389 nspins = SIZE(gs_mos, 1)
390 nvects = SIZE(evects, 2)
391 do_hfx = tddfpt_control%do_hfx
392 do_admm = tddfpt_control%do_admm
393 IF (do_admm) THEN
394 kernel_env_admm_aux => kernel_env%admm_kernel
395 ELSE
396 NULLIFY (kernel_env_admm_aux)
397 END IF
398 is_rks_triplets = tddfpt_control%rks_triplets
399 do_lri_response = tddfpt_control%do_lrigpw
400 do_bse = tddfpt_control%do_bse
401 do_bse_w_only = tddfpt_control%do_bse_w_only
402 do_bse_gw_only = tddfpt_control%do_bse_gw_only
403 IF (do_bse) do_hfx = .false.
404 IF (do_bse_w_only) do_hfx = .false.
405 IF (do_bse_gw_only) do_hfx = .false.
406
407 IF (debug_this_module) THEN
408 cpassert(nspins > 0)
409 cpassert(SIZE(aop_evects, 1) == SIZE(evects, 1))
410 cpassert(SIZE(s_evects, 1) == SIZE(evects, 1))
411 cpassert(SIZE(aop_evects, 2) == nvects)
412 cpassert(SIZE(s_evects, 2) == nvects)
413 cpassert(SIZE(gs_mos) == nspins)
414 END IF
415
416 DO ispin = 1, nspins
417 nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
418 END DO
419
420 IF (nvects > 0) THEN
421 CALL cp_fm_get_info(evects(1, 1), para_env=para_env)
422 IF (ALLOCATED(work_matrices%evects_sub)) THEN
423 DO ivect = 1, nvects
424 DO ispin = 1, SIZE(evects, 1)
425 associate(evect => evects(ispin, ivect), work_matrix => work_matrices%evects_sub(ispin, ivect))
426 IF (ASSOCIATED(evect%matrix_struct)) THEN
427 IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
428 CALL cp_fm_copy_general(evect, work_matrix, para_env)
429 ELSE
430 CALL cp_fm_copy_general(evect, fm_dummy, para_env)
431 END IF
432 ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
433 CALL cp_fm_copy_general(fm_dummy, work_matrix, para_env)
434 ELSE
435 CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env)
436 END IF
437 END associate
438 END DO
439 END DO
440 END IF
441
442 IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
443 ! full TDDFPT kernel
444 CALL fhxc_kernel(aop_evects, evects, is_rks_triplets, do_hfx, do_admm, qs_env, &
445 kernel_env%full_kernel, kernel_env_admm_aux, sub_env, work_matrices, &
446 tddfpt_control%admm_symm, tddfpt_control%admm_xc_correction, &
447 do_lri_response, tddfpt_mgrid=tddfpt_control%mgrid_is_explicit)
448 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
449 ! sTDA kernel
450 CALL stda_kernel(aop_evects, evects, is_rks_triplets, qs_env, tddfpt_control%stda_control, &
451 kernel_env%stda_kernel, sub_env, work_matrices)
452 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
453 ! No kernel
454 DO ivect = 1, nvects
455 DO ispin = 1, SIZE(evects, 1)
456 CALL cp_fm_set_all(aop_evects(ispin, ivect), 0.0_dp)
457 END DO
458 END DO
459 ELSE
460 cpabort("Kernel type not implemented")
461 END IF
462
463 IF (ALLOCATED(work_matrices%evects_sub)) THEN
464 DO ivect = 1, nvects
465 DO ispin = 1, SIZE(evects, 1)
466 associate(aop_evect => aop_evects(ispin, ivect), &
467 work_matrix => work_matrices%Aop_evects_sub(ispin, ivect))
468 IF (ASSOCIATED(aop_evect%matrix_struct)) THEN
469 IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
470 CALL cp_fm_copy_general(work_matrix, aop_evect, para_env)
471 ELSE
472 CALL cp_fm_copy_general(fm_dummy, aop_evect, para_env)
473 END IF
474 ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
475 CALL cp_fm_copy_general(work_matrix, fm_dummy, para_env)
476 ELSE
477 CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env)
478 END IF
479 END associate
480 END DO
481 END DO
482 END IF
483
484 ! orbital energy difference term
485 IF ((.NOT. do_bse) .AND. (.NOT. do_bse_gw_only)) THEN
486 CALL tddfpt_apply_energy_diff(aop_evects=aop_evects, evects=evects, s_evects=s_evects, &
487 gs_mos=gs_mos, matrix_ks=matrix_ks, tddfpt_control=tddfpt_control)
488 ELSE
489 CALL zeroth_order_gw(qs_env=qs_env, aop_evects=aop_evects, evects=evects, s_evects=s_evects, &
490 gs_mos=gs_mos, matrix_s=matrix_s, matrix_ks=matrix_ks)
491 END IF
492
493 ! if smeared occupation, then add aCCSX here
494 IF (tddfpt_control%do_smearing) THEN
495 DO ivect = 1, nvects
496 DO ispin = 1, nspins
497 CALL add_smearing_aterm(qs_env, aop_evects(ispin, ivect), evects(ispin, ivect), &
498 s_evects(ispin, ivect), gs_mos(ispin)%mos_occ, &
499 tddfpt_control%smeared_occup(ispin)%fermia, matrix_s)
500 END DO
501 END DO
502 END IF
503
504 IF (do_hfx) THEN
505 IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
506 ! full TDDFPT kernel
507 CALL tddfpt_apply_hfx(aop_evects=aop_evects, evects=evects, gs_mos=gs_mos, do_admm=do_admm, &
508 qs_env=qs_env, wfm_rho_orb=work_matrices%hfx_fm_ao_ao, &
509 work_hmat_symm=work_matrices%hfx_hmat_symm, &
510 work_hmat_asymm=work_matrices%hfx_hmat_asymm, &
511 work_rho_ia_ao_symm=work_matrices%hfx_rho_ao_symm, &
512 work_rho_ia_ao_asymm=work_matrices%hfx_rho_ao_asymm)
513 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
514 ! sTDA kernel
515 ! special treatment of HFX term
516 ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
517 ! No kernel
518 ! drop kernel contribution of HFX term
519 ELSE
520 cpabort("Kernel type not implemented")
521 END IF
522 END IF
523 ! short/long range HFX
524 IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
525 IF (tddfpt_control%do_hfxsr) THEN
526 re_int = tddfpt_control%hfxsr_re_int
527 ! symmetric dmat
528 CALL tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, &
529 kernel_env%full_kernel%admm_env, &
530 kernel_env%full_kernel%hfxsr_section, &
531 kernel_env%full_kernel%x_data, 1, re_int, &
532 work_rho_ia_ao=work_matrices%hfxsr_rho_ao_symm, &
533 work_hmat=work_matrices%hfxsr_hmat_symm, &
534 wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
535 ! antisymmetric dmat
536 CALL tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, &
537 kernel_env%full_kernel%admm_env, &
538 kernel_env%full_kernel%hfxsr_section, &
539 kernel_env%full_kernel%x_data, -1, .false., &
540 work_rho_ia_ao=work_matrices%hfxsr_rho_ao_asymm, &
541 work_hmat=work_matrices%hfxsr_hmat_asymm, &
542 wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
543 tddfpt_control%hfxsr_re_int = .false.
544 END IF
545 IF (tddfpt_control%do_hfxlr) THEN
546 rcut = tddfpt_control%hfxlr_rcut
547 scale = tddfpt_control%hfxlr_scale
548 DO ivect = 1, nvects
549 IF (ALLOCATED(work_matrices%evects_sub)) THEN
550 IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct)) THEN
551 CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
552 work_matrices%evects_sub(:, ivect), &
553 work_matrices%Aop_evects_sub(:, ivect))
554 ELSE
555 ! skip trial vectors which are assigned to different parallel groups
556 cycle
557 END IF
558 ELSE
559 CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
560 evects(:, ivect), aop_evects(:, ivect))
561 END IF
562 END DO
563 END IF
564 END IF
565 IF ((do_bse) .OR. (do_bse_w_only)) THEN
566 ! add dynamical screening
567 ! CALL tddfpt_apply_bse_debug(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, qs_env=qs_env)
568 CALL tddfpt_apply_bse(aop_evects=aop_evects, evects=evects, gs_mos=gs_mos, qs_env=qs_env, &
569 s_evects=s_evects)
570 END IF
571
572 END IF
573
574 CALL timestop(handle)
575
576 END SUBROUTINE tddfpt_compute_aop_evects
577
578! **************************************************************************************************
579!> \brief Solve eigenproblem for the reduced action matrix and find new Ritz eigenvectors and
580!> eigenvalues.
581!> \param ritz_vects Ritz eigenvectors (initialised on exit)
582!> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors (initialised on exit)
583!> \param evals Ritz eigenvalues (initialised on exit)
584!> \param krylov_vects Krylov's vectors
585!> \param Aop_krylov action of TDDFPT operator on Krylov's vectors
586!> \param Atilde TDDFPT matrix projected into the Krylov's vectors subspace
587!> \param nkvo ...
588!> \param nkvn ...
589!> \par History
590!> * 06.2016 created [Sergey Chulkov]
591!> * 03.2017 altered prototype, OpenMP parallelisation [Sergey Chulkov]
592! **************************************************************************************************
593 SUBROUTINE tddfpt_compute_ritz_vects(ritz_vects, Aop_ritz, evals, krylov_vects, Aop_krylov, &
594 Atilde, nkvo, nkvn)
595 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: ritz_vects, aop_ritz
596 REAL(kind=dp), DIMENSION(:), INTENT(out) :: evals
597 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: krylov_vects, aop_krylov
598 REAL(kind=dp), DIMENSION(:, :), POINTER :: atilde
599 INTEGER, INTENT(IN) :: nkvo, nkvn
600
601 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_compute_ritz_vects'
602
603 INTEGER :: handle, ikv, irv, ispin, nkvs, nrvs, &
604 nspins
605 REAL(kind=dp) :: act
606 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: at12, at21, at22, evects_atilde
607 TYPE(cp_blacs_env_type), POINTER :: blacs_env_global
608 TYPE(cp_fm_struct_type), POINTER :: fm_struct
609 TYPE(cp_fm_type) :: atilde_fm, evects_atilde_fm
610
611 CALL timeset(routinen, handle)
612
613 nspins = SIZE(krylov_vects, 1)
614 nkvs = SIZE(krylov_vects, 2)
615 nrvs = SIZE(ritz_vects, 2)
616 cpassert(nkvs == nkvo + nkvn)
617
618 CALL cp_fm_get_info(krylov_vects(1, 1), context=blacs_env_global)
619
620 CALL cp_fm_struct_create(fm_struct, nrow_global=nkvs, ncol_global=nkvs, context=blacs_env_global)
621 CALL cp_fm_create(atilde_fm, fm_struct, set_zero=.true.)
622 CALL cp_fm_create(evects_atilde_fm, fm_struct, set_zero=.true.)
623 CALL cp_fm_struct_release(fm_struct)
624
625 ! *** compute upper-diagonal reduced action matrix ***
626 CALL reallocate(atilde, 1, nkvs, 1, nkvs)
627 ! TO DO: the subroutine 'cp_fm_contracted_trace' will compute all elements of
628 ! the matrix 'Atilde', however only upper-triangular elements are actually needed
629 !
630 IF (nkvo == 0) THEN
631 CALL cp_fm_contracted_trace(aop_krylov(:, 1:nkvs), krylov_vects(:, 1:nkvs), &
632 atilde(1:nkvs, 1:nkvs), accurate=.false.)
633 ELSE
634 ALLOCATE (at12(nkvn, nkvo), at21(nkvo, nkvn), at22(nkvn, nkvn))
635 CALL cp_fm_contracted_trace(aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, 1:nkvo), &
636 at12, accurate=.false.)
637 atilde(nkvo + 1:nkvs, 1:nkvo) = at12(1:nkvn, 1:nkvo)
638 CALL cp_fm_contracted_trace(aop_krylov(:, 1:nkvo), krylov_vects(:, nkvo + 1:nkvs), &
639 at21, accurate=.false.)
640 atilde(1:nkvo, nkvo + 1:nkvs) = at21(1:nkvo, 1:nkvn)
641 CALL cp_fm_contracted_trace(aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, nkvo + 1:nkvs), &
642 at22, accurate=.false.)
643 atilde(nkvo + 1:nkvs, nkvo + 1:nkvs) = at22(1:nkvn, 1:nkvn)
644 DEALLOCATE (at12, at21, at22)
645 END IF
646 atilde = 0.5_dp*(atilde + transpose(atilde))
647 CALL cp_fm_set_submatrix(atilde_fm, atilde)
648
649 ! *** solve an eigenproblem for the reduced matrix ***
650 CALL choose_eigv_solver(atilde_fm, evects_atilde_fm, evals(1:nkvs))
651
652 ALLOCATE (evects_atilde(nkvs, nrvs))
653 CALL cp_fm_get_submatrix(evects_atilde_fm, evects_atilde, start_row=1, start_col=1, n_rows=nkvs, n_cols=nrvs)
654 CALL cp_fm_release(evects_atilde_fm)
655 CALL cp_fm_release(atilde_fm)
656
657!$OMP PARALLEL DO DEFAULT(NONE), &
658!$OMP PRIVATE(act, ikv, irv, ispin), &
659!$OMP SHARED(Aop_krylov, Aop_ritz, krylov_vects, evects_Atilde, nkvs, nrvs, nspins, ritz_vects)
660 DO irv = 1, nrvs
661 DO ispin = 1, nspins
662 CALL cp_fm_set_all(ritz_vects(ispin, irv), 0.0_dp)
663 CALL cp_fm_set_all(aop_ritz(ispin, irv), 0.0_dp)
664 END DO
665
666 DO ikv = 1, nkvs
667 act = evects_atilde(ikv, irv)
668 DO ispin = 1, nspins
669 CALL cp_fm_scale_and_add(1.0_dp, ritz_vects(ispin, irv), &
670 act, krylov_vects(ispin, ikv))
671 CALL cp_fm_scale_and_add(1.0_dp, aop_ritz(ispin, irv), &
672 act, aop_krylov(ispin, ikv))
673 END DO
674 END DO
675 END DO
676!$OMP END PARALLEL DO
677
678 DEALLOCATE (evects_atilde)
679
680 CALL timestop(handle)
681
682 END SUBROUTINE tddfpt_compute_ritz_vects
683
684! **************************************************************************************************
685!> \brief Expand Krylov space by computing residual vectors.
686!> \param residual_vects residual vectors (modified on exit)
687!> \param evals Ritz eigenvalues (modified on exit)
688!> \param ritz_vects Ritz eigenvectors
689!> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors
690!> \param gs_mos molecular orbitals optimised for the ground state
691!> \param matrix_s overlap matrix
692!> \param tddfpt_control ...
693!> \par History
694!> * 06.2016 created [Sergey Chulkov]
695!> * 03.2017 refactored to achieve significant performance gain [Sergey Chulkov]
696! **************************************************************************************************
697 SUBROUTINE tddfpt_compute_residual_vects(residual_vects, evals, ritz_vects, Aop_ritz, gs_mos, &
698 matrix_s, tddfpt_control)
699 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: residual_vects
700 REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
701 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: ritz_vects, aop_ritz
702 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
703 INTENT(in) :: gs_mos
704 TYPE(dbcsr_type), POINTER :: matrix_s
705 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
706
707 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_compute_residual_vects'
708 REAL(kind=dp), PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*epsilon(1.0_dp)
709
710 INTEGER :: handle, ica, icb, icol_local, &
711 irow_local, irv, ispin, nao, &
712 ncols_local, nrows_local, nrvs, &
713 nspins, spin2, spinflip
714 INTEGER, DIMENSION(:), POINTER :: col_indices_local, row_indices_local
715 INTEGER, DIMENSION(maxspins) :: nactive, nmo_virt
716 REAL(kind=dp) :: e_occ_plus_lambda, eref, lambda
717 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
718 POINTER :: weights_ldata
719 TYPE(cp_fm_struct_type), POINTER :: ao_mo_struct, virt_mo_struct
720 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: awork, vomat
721
722 CALL timeset(routinen, handle)
723
724 nspins = SIZE(residual_vects, 1)
725 nrvs = SIZE(residual_vects, 2)
726 spinflip = tddfpt_control%spinflip
727
728 IF (nrvs > 0) THEN
729 CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
730 ALLOCATE (awork(nspins), vomat(nspins))
731 DO ispin = 1, nspins
732 IF (spinflip == no_sf_tddfpt) THEN
733 spin2 = ispin
734 ELSE
735 spin2 = 2
736 END IF
737 nmo_virt(spin2) = SIZE(gs_mos(spin2)%evals_virt)
738 !
739 CALL cp_fm_get_info(matrix=ritz_vects(ispin, 1), matrix_struct=ao_mo_struct, &
740 ncol_global=nactive(ispin))
741 CALL cp_fm_create(awork(ispin), ao_mo_struct)
742 !
743 CALL cp_fm_struct_create(virt_mo_struct, nrow_global=nmo_virt(spin2), &
744 ncol_global=nactive(ispin), template_fmstruct=ao_mo_struct)
745 CALL cp_fm_create(vomat(ispin), virt_mo_struct)
746 CALL cp_fm_struct_release(virt_mo_struct)
747 END DO
748
749 ! *** actually compute residual vectors ***
750 DO irv = 1, nrvs
751 lambda = evals(irv)
752
753 DO ispin = 1, nspins
754 IF (spinflip == no_sf_tddfpt) THEN
755 spin2 = ispin
756 ELSE
757 spin2 = 2
758 END IF
759 CALL cp_fm_get_info(vomat(ispin), nrow_local=nrows_local, &
760 ncol_local=ncols_local, row_indices=row_indices_local, &
761 col_indices=col_indices_local, local_data=weights_ldata)
762
763 ! awork := Ab(ispin, irv) - evals(irv) b(ispin, irv), where 'b' is a Ritz vector
764 CALL cp_dbcsr_sm_fm_multiply(matrix_s, ritz_vects(ispin, irv), awork(ispin), &
765 ncol=nactive(ispin), alpha=-lambda, beta=0.0_dp)
766 CALL cp_fm_scale_and_add(1.0_dp, awork(ispin), 1.0_dp, aop_ritz(ispin, irv))
767 !
768 CALL parallel_gemm('T', 'N', nmo_virt(spin2), nactive(ispin), nao, 1.0_dp, gs_mos(spin2)%mos_virt, &
769 awork(ispin), 0.0_dp, vomat(ispin))
770
771 ! Apply Davidson preconditioner to the residue vectors vomat to obtain new directions
772 DO icol_local = 1, ncols_local
773 ica = col_indices_local(icol_local)
774 icb = gs_mos(ispin)%index_active(ica)
775 e_occ_plus_lambda = gs_mos(ispin)%evals_occ(icb) + lambda
776
777 DO irow_local = 1, nrows_local
778 eref = gs_mos(spin2)%evals_virt(row_indices_local(irow_local)) - e_occ_plus_lambda
779
780 ! eref = e_virt - e_occ - lambda = e_virt - e_occ - (eref_scale*lambda + (1-eref_scale)*lambda);
781 ! eref_new = e_virt - e_occ - eref_scale*lambda = eref + (1 - eref_scale)*lambda
782 IF (abs(eref) < threshold) &
783 eref = eref + (1.0_dp - eref_scale)*lambda
784
785 weights_ldata(irow_local, icol_local) = weights_ldata(irow_local, icol_local)/eref
786 END DO
787 END DO
788
789 CALL parallel_gemm('N', 'N', nao, nactive(ispin), nmo_virt(spin2), 1.0_dp, gs_mos(spin2)%mos_virt, &
790 vomat(ispin), 0.0_dp, residual_vects(ispin, irv))
791 END DO
792 END DO
793 !
794 CALL cp_fm_release(awork)
795 CALL cp_fm_release(vomat)
796 END IF
797
798 CALL timestop(handle)
799
800 END SUBROUTINE tddfpt_compute_residual_vects
801
802! **************************************************************************************************
803!> \brief Perform Davidson iterations.
804!> \param evects TDDFPT trial vectors (modified on exit)
805!> \param evals TDDFPT eigenvalues (modified on exit)
806!> \param S_evects cached matrix product S * evects (modified on exit)
807!> \param gs_mos molecular orbitals optimised for the ground state
808!> \param tddfpt_control TDDFPT control parameters
809!> \param matrix_ks Kohn-Sham matrix
810!> \param qs_env Quickstep environment
811!> \param kernel_env kernel environment
812!> \param sub_env parallel (sub)group environment
813!> \param logger CP2K logger
814!> \param iter_unit I/O unit to write basic iteration information
815!> \param energy_unit I/O unit to write detailed energy information
816!> \param tddfpt_print_section TDDFPT print input section (need to write TDDFPT restart files)
817!> \param work_matrices collection of work matrices (modified on exit)
818!> \return energy convergence achieved (in Hartree)
819!> \par History
820!> * 03.2017 code related to Davidson eigensolver has been moved here from the main subroutine
821!> tddfpt() [Sergey Chulkov]
822!> \note Based on the subroutines apply_op() and iterative_solver() originally created by
823!> Thomas Chassaing in 2002.
824! **************************************************************************************************
825 FUNCTION tddfpt_davidson_solver(evects, evals, S_evects, gs_mos, tddfpt_control, &
826 matrix_ks, qs_env, kernel_env, &
827 sub_env, logger, iter_unit, energy_unit, &
828 tddfpt_print_section, work_matrices) RESULT(conv)
829 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout) :: evects
830 REAL(kind=dp), DIMENSION(:), INTENT(inout) :: evals
831 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(inout) :: s_evects
832 TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
833 INTENT(in) :: gs_mos
834 TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
835 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
836 TYPE(qs_environment_type), POINTER :: qs_env
837 TYPE(kernel_env_type), INTENT(in) :: kernel_env
838 TYPE(tddfpt_subgroup_env_type), INTENT(in) :: sub_env
839 TYPE(cp_logger_type), POINTER :: logger
840 INTEGER, INTENT(in) :: iter_unit, energy_unit
841 TYPE(section_vals_type), POINTER :: tddfpt_print_section
842 TYPE(tddfpt_work_matrices), INTENT(inout) :: work_matrices
843 REAL(kind=dp) :: conv
844
845 CHARACTER(LEN=*), PARAMETER :: routinen = 'tddfpt_davidson_solver'
846
847 INTEGER :: handle, ispin, istate, iter, &
848 max_krylov_vects, nspins, nstates, &
849 nstates_conv, nvects_exist, nvects_new
850 INTEGER(kind=int_8) :: nstates_total
851 LOGICAL :: is_nonortho
852 REAL(kind=dp) :: t1, t2
853 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_last
854 REAL(kind=dp), DIMENSION(:, :), POINTER :: atilde
855 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: aop_krylov, aop_ritz, krylov_vects, &
856 s_krylov
857 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
858
859 CALL timeset(routinen, handle)
860
861 nspins = SIZE(evects, 1)
862 nstates = tddfpt_control%nstates
863 nstates_total = tddfpt_total_number_of_states(tddfpt_control, gs_mos)
864
865 IF (debug_this_module) THEN
866 cpassert(SIZE(evects, 1) == nspins)
867 cpassert(SIZE(evects, 2) == nstates)
868 cpassert(SIZE(evals) == nstates)
869 END IF
870
871 CALL get_qs_env(qs_env, matrix_s=matrix_s)
872
873 ! adjust the number of Krylov vectors
874 max_krylov_vects = min(max(tddfpt_control%nkvs, nstates), int(nstates_total))
875
876 ALLOCATE (aop_ritz(nspins, nstates))
877 DO istate = 1, nstates
878 DO ispin = 1, nspins
879 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, aop_ritz(ispin, istate))
880 END DO
881 END DO
882
883 ALLOCATE (evals_last(max_krylov_vects))
884 ALLOCATE (aop_krylov(nspins, max_krylov_vects), krylov_vects(nspins, max_krylov_vects), &
885 s_krylov(nspins, max_krylov_vects))
886
887 DO istate = 1, nstates
888 DO ispin = 1, nspins
889 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, krylov_vects(ispin, istate))
890 CALL cp_fm_to_fm(evects(ispin, istate), krylov_vects(ispin, istate))
891
892 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, s_krylov(ispin, istate))
893 CALL cp_fm_to_fm(s_evects(ispin, istate), s_krylov(ispin, istate))
894
895 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, aop_krylov(ispin, istate))
896 END DO
897 END DO
898
899 nvects_exist = 0
900 nvects_new = nstates
901
902 t1 = m_walltime()
903
904 ALLOCATE (atilde(1, 1))
905
906 DO
907 ! davidson iteration
908 CALL cp_iterate(logger%iter_info, iter_nr_out=iter)
909
910 ! Matrix-vector operations
911 CALL tddfpt_compute_aop_evects(aop_evects=aop_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), &
912 evects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
913 s_evects=s_krylov(:, nvects_exist + 1:nvects_exist + nvects_new), &
914 gs_mos=gs_mos, tddfpt_control=tddfpt_control, &
915 matrix_ks=matrix_ks, &
916 qs_env=qs_env, kernel_env=kernel_env, &
917 sub_env=sub_env, &
918 work_matrices=work_matrices, &
919 matrix_s=matrix_s(1)%matrix)
920
921 CALL tddfpt_compute_ritz_vects(ritz_vects=evects, aop_ritz=aop_ritz, &
922 evals=evals_last(1:nvects_exist + nvects_new), &
923 krylov_vects=krylov_vects(:, 1:nvects_exist + nvects_new), &
924 aop_krylov=aop_krylov(:, 1:nvects_exist + nvects_new), &
925 atilde=atilde, nkvo=nvects_exist, nkvn=nvects_new)
926
927 CALL tddfpt_write_restart(evects=evects, evals=evals_last(1:nstates), gs_mos=gs_mos, &
928 logger=logger, tddfpt_print_section=tddfpt_print_section)
929
930 conv = maxval(abs(evals_last(1:nstates) - evals(1:nstates)))
931
932 nvects_exist = nvects_exist + nvects_new
933 IF (nvects_exist + nvects_new > max_krylov_vects) &
934 nvects_new = max_krylov_vects - nvects_exist
935 IF (iter >= tddfpt_control%niters) nvects_new = 0
936
937 IF (conv > tddfpt_control%conv .AND. nvects_new > 0) THEN
938 ! compute residual vectors for the next iteration
939 DO istate = 1, nvects_new
940 DO ispin = 1, nspins
941 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
942 krylov_vects(ispin, nvects_exist + istate))
943 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
944 s_krylov(ispin, nvects_exist + istate))
945 CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
946 aop_krylov(ispin, nvects_exist + istate))
947 END DO
948 END DO
949
950 CALL tddfpt_compute_residual_vects(residual_vects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
951 evals=evals_last(1:nvects_new), &
952 ritz_vects=evects(:, 1:nvects_new), aop_ritz=aop_ritz(:, 1:nvects_new), &
953 gs_mos=gs_mos, matrix_s=matrix_s(1)%matrix, tddfpt_control=tddfpt_control)
954
955 CALL tddfpt_orthogonalize_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
956 work_matrices%S_C0_C0T, qs_env, &
957 gs_mos, evals(1:nstates), tddfpt_control, work_matrices%S_C0)
958
959 CALL tddfpt_orthonormalize_psi1_psi1(krylov_vects(:, 1:nvects_exist + nvects_new), nvects_new, &
960 s_krylov(:, 1:nvects_exist + nvects_new), matrix_s(1)%matrix)
961
962 is_nonortho = tddfpt_is_nonorthogonal_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
963 work_matrices%S_C0, tddfpt_control%orthogonal_eps, &
964 tddfpt_control%spinflip)
965 ELSE
966 ! convergence or the maximum number of Krylov vectors have been achieved
967 nvects_new = 0
968 is_nonortho = .false.
969 END IF
970
971 t2 = m_walltime()
972 IF (energy_unit > 0) THEN
973 WRITE (energy_unit, '(/,4X,A,T14,A,T36,A)') "State", "Exc. energy (eV)", "Convergence (eV)"
974 DO istate = 1, nstates
975 WRITE (energy_unit, '(1X,I8,T12,F14.7,T38,ES11.4)') istate, &
976 evals_last(istate)*evolt, (evals_last(istate) - evals(istate))*evolt
977 END DO
978 WRITE (energy_unit, *)
979 CALL m_flush(energy_unit)
980 END IF
981
982 IF (iter_unit > 0) THEN
983 nstates_conv = 0
984 DO istate = 1, nstates
985 IF (abs(evals_last(istate) - evals(istate)) <= tddfpt_control%conv) &
986 nstates_conv = nstates_conv + 1
987 END DO
988
989 WRITE (iter_unit, '(T7,I8,T24,F7.1,T40,ES11.4,T66,I8)') iter, t2 - t1, conv, nstates_conv
990 CALL m_flush(iter_unit)
991 END IF
992
993 t1 = t2
994 evals(1:nstates) = evals_last(1:nstates)
995
996 ! nvects_new == 0 if iter >= tddfpt_control%niters
997 IF (nvects_new == 0 .OR. is_nonortho) THEN
998 ! restart Davidson iterations
999 CALL tddfpt_orthogonalize_psi1_psi0(evects, work_matrices%S_C0_C0T, qs_env, &
1000 gs_mos, &
1001 evals(1:nstates), tddfpt_control, work_matrices%S_C0)
1002 CALL tddfpt_orthonormalize_psi1_psi1(evects, nstates, s_evects, matrix_s(1)%matrix)
1003
1004 EXIT
1005 END IF
1006 END DO
1007
1008 DEALLOCATE (atilde)
1009
1010 DO istate = nvects_exist + nvects_new, 1, -1
1011 DO ispin = nspins, 1, -1
1012 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, aop_krylov(ispin, istate))
1013 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, s_krylov(ispin, istate))
1014 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, krylov_vects(ispin, istate))
1015 END DO
1016 END DO
1017 DEALLOCATE (aop_krylov, krylov_vects, s_krylov)
1018 DEALLOCATE (evals_last)
1019
1020 DO istate = nstates, 1, -1
1021 DO ispin = nspins, 1, -1
1022 CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, aop_ritz(ispin, istate))
1023 END DO
1024 END DO
1025 DEALLOCATE (aop_ritz)
1026
1027 CALL timestop(handle)
1028
1029 END FUNCTION tddfpt_davidson_solver
1030
1031END MODULE qs_tddfpt2_eigensolver
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
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)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:239
pool for for elements that are retained and released
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
subroutine, public fm_pool_give_back_fm(pool, element)
returns the element to the pool
represent the structure of a full matrix
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
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_copy_general(source, destination, para_env)
General copy of a fm matrix to another fm matrix. Uses non-blocking MPI rather than ScaLAPACK.
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
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
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...
subroutine, public cp_fm_maxabsval(matrix, a_max, ir_max, ic_max)
find the maximum absolute value of the matrix element maxval(abs(matrix))
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
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,...
various routines to log and control the output. The idea is that decisions about where to log should ...
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
subroutine, public cp_iterate(iteration_info, last, iter_nr, increment, iter_nr_out)
adds one to the actual iteration
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public tddfpt_kernel_none
integer, parameter, public tddfpt_kernel_full
integer, parameter, public tddfpt_kernel_stda
integer, parameter, public no_sf_tddfpt
objects that represent the structure of input sections and the data contained in an input section
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public int_8
Definition kinds.F:54
integer, parameter, public dp
Definition kinds.F:34
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:136
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:153
Utility routines for the memory handling.
Interface to the message passing library MPI.
basic linear algebra operations for full matrixes
Definition of physical constants:
Definition physcon.F:68
real(kind=dp), parameter, public evolt
Definition physcon.F:183
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
subroutine, public eigensolver(matrix_ks_fm, mo_set, ortho, work, cholesky_method, do_level_shift, level_shift, matrix_u_fm, use_jacobi)
Diagonalise the Kohn-Sham matrix to get a new set of MO eigen- vectors and MO eigenvalues....
subroutine, public zeroth_order_gw(qs_env, aop_evects, evects, s_evects, gs_mos, matrix_s, matrix_ks)
...
subroutine, public tddfpt_apply_bse(aop_evects, evects, gs_mos, qs_env, s_evects)
...
subroutine, public tddfpt_orthonormalize_psi1_psi1(evects, nvects_new, s_evects, matrix_s)
Make new TDDFPT trial vectors orthonormal to all previous trial vectors.
subroutine, public 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.
real(kind=dp) function, public tddfpt_davidson_solver(evects, evals, s_evects, gs_mos, tddfpt_control, matrix_ks, qs_env, kernel_env, sub_env, logger, iter_unit, energy_unit, tddfpt_print_section, work_matrices)
Perform Davidson iterations.
subroutine, public 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.
subroutine, public stda_kernel(aop_evects, evects, is_rks_triplets, qs_env, stda_control, stda_env, sub_env, work_matrices)
Compute action matrix-vector products with the sTDA Kernel.
subroutine, public tddfpt_apply_hfx(aop_evects, evects, gs_mos, do_admm, qs_env, work_rho_ia_ao_symm, work_hmat_symm, work_rho_ia_ao_asymm, work_hmat_asymm, wfm_rho_orb)
Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
subroutine, public tddfpt_apply_energy_diff(aop_evects, evects, s_evects, gs_mos, matrix_ks, tddfpt_control)
Apply orbital energy difference term: Aop_evects(spin,state) += KS(spin) * evects(spin,...
subroutine, public tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, hfx_scale, work, x, res)
...Calculate the HFXLR kernel contribution by contracting the Lowdin MO coefficients – transition cha...
subroutine, public tddfpt_apply_hfxsr_kernel(aop_evects, evects, gs_mos, qs_env, admm_env, hfx_section, x_data, symmetry, recalc_integrals, work_rho_ia_ao, work_hmat, wfm_rho_orb)
Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
subroutine, public tddfpt_write_restart(evects, evals, gs_mos, logger, tddfpt_print_section)
Write Ritz vectors to a binary restart file.
subroutine, public add_smearing_aterm(qs_env, aop_evects, evects, s_evects, mos_occ, fermia, matrix_s)
...
subroutine, public compute_fermib(qs_env, gs_mos, evals)
...
subroutine, public orthogonalize_smeared_occupation(evects, qs_env, mos_occ, s_c0)
...
pure integer(kind=int_8) function, public tddfpt_total_number_of_states(tddfpt_control, gs_mos)
Compute the number of possible singly excited states (occ -> virt)
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment
Collection of variables required to evaluate adiabatic TDDFPT kernel.
Type to hold environments for the different kernels.
Ground state molecular orbitals.
Set of temporary ("work") matrices.