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