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