Line data Source code
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 :
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: no_sf_tddfpt,&
31 : tddfpt_kernel_full,&
32 : tddfpt_kernel_none,&
33 : tddfpt_kernel_stda
34 : USE input_section_types, ONLY: section_vals_type
35 : USE kinds, ONLY: dp,&
36 : int_8
37 : USE machine, ONLY: m_flush,&
38 : m_walltime
39 : USE memory_utilities, ONLY: reallocate
40 : USE message_passing, ONLY: mp_para_env_type
41 : USE parallel_gemm_api, ONLY: parallel_gemm
42 : USE physcon, ONLY: evolt
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_kernel_types, ONLY: full_kernel_env_type,&
46 : kernel_env_type
47 : USE qs_scf_methods, ONLY: eigensolver
48 : USE qs_tddfpt2_bse_utils, ONLY: tddfpt_apply_bse,&
49 : zeroth_order_gw
50 : USE qs_tddfpt2_fhxc, ONLY: fhxc_kernel,&
51 : stda_kernel
52 : USE qs_tddfpt2_operators, ONLY: tddfpt_apply_energy_diff,&
53 : tddfpt_apply_hfx,&
54 : tddfpt_apply_hfxlr_kernel,&
55 : tddfpt_apply_hfxsr_kernel
56 : USE qs_tddfpt2_restart, ONLY: tddfpt_write_restart
57 : USE qs_tddfpt2_smearing_methods, ONLY: add_smearing_aterm,&
58 : compute_fermib,&
59 : orthogonalize_smeared_occupation
60 : USE qs_tddfpt2_subgroups, ONLY: tddfpt_subgroup_env_type
61 : USE qs_tddfpt2_types, ONLY: tddfpt_ground_state_mos,&
62 : tddfpt_work_matrices
63 : USE qs_tddfpt2_utils, ONLY: tddfpt_total_number_of_states
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 :
77 : PUBLIC :: tddfpt_davidson_solver, tddfpt_orthogonalize_psi1_psi0, &
78 : tddfpt_orthonormalize_psi1_psi1
79 :
80 : ! **************************************************************************************************
81 :
82 : CONTAINS
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 8340 : SUBROUTINE tddfpt_orthogonalize_psi1_psi0(evects, S_C0_C0T, qs_env, gs_mos, evals, &
106 8340 : 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 8340 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos
123 : TYPE(mp_para_env_type), POINTER :: para_env
124 :
125 8340 : CALL timeset(routineN, handle)
126 :
127 8340 : nspins = SIZE(evects, 1)
128 8340 : nvects = SIZE(evects, 2)
129 :
130 8340 : IF (nvects > 0) THEN
131 8340 : IF (.NOT. tddfpt_control%do_smearing) THEN
132 17890 : DO ispin = 1, nspins
133 9564 : IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
134 : spin = ispin
135 : ELSE
136 166 : spin = 2
137 : END IF
138 : CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
139 9564 : nrow_global=nao, ncol_global=nactive)
140 9564 : CALL cp_fm_create(evortho, matrix_struct)
141 34536 : 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 24972 : evects(ispin, ivect), 0.0_dp, evortho)
145 34536 : CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, ivect), -1.0_dp, evortho)
146 : END DO
147 27454 : CALL cp_fm_release(evortho)
148 : END DO
149 : ELSE
150 14 : NULLIFY (para_env)
151 14 : CALL get_qs_env(qs_env, para_env=para_env)
152 14 : NULLIFY (mos)
153 56 : ALLOCATE (mos(nspins))
154 28 : DO ispin = 1, nspins
155 : CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
156 14 : nrow_global=nao, ncol_global=nactive)
157 14 : CALL cp_fm_create(mos(ispin), matrix_struct)
158 42 : CALL cp_fm_copy_general(gs_mos(ispin)%mos_occ, mos(ispin), para_env)
159 : END DO
160 28 : DO ivect = 1, nvects
161 14 : CALL compute_fermib(qs_env, gs_mos, evals(ivect))
162 28 : CALL orthogonalize_smeared_occupation(evects(:, ivect), qs_env, mos, S_C0)
163 : END DO
164 28 : DO ispin = 1, nspins
165 28 : CALL cp_fm_release(mos(ispin))
166 : END DO
167 14 : DEALLOCATE (mos)
168 : END IF
169 : END IF
170 :
171 8340 : CALL timestop(handle)
172 :
173 8340 : 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 5352 : 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 5352 : CALL timeset(routineN, handle)
206 :
207 5352 : nspins = SIZE(evects, 1)
208 5352 : nvects = SIZE(evects, 2)
209 :
210 5352 : is_nonortho = .FALSE.
211 :
212 11560 : loop: DO ispin = 1, nspins
213 :
214 6218 : IF (spinflip == no_sf_tddfpt) THEN
215 : spin = ispin
216 : ELSE
217 106 : spin = 2
218 : END IF
219 :
220 6218 : 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 6218 : nrow_global=nao, ncol_global=nactive)
223 : CALL cp_fm_struct_create(matrix_struct_tmp, nrow_global=nocc, &
224 6218 : ncol_global=nactive, template_fmstruct=matrix_struct)
225 6218 : CALL cp_fm_create(aortho, matrix_struct_tmp)
226 6218 : CALL cp_fm_struct_release(matrix_struct_tmp)
227 22070 : 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 15862 : evects(ispin, ivect), 0.0_dp, aortho)
231 15862 : CALL cp_fm_maxabsval(aortho, maxabs_val)
232 15862 : is_nonortho = maxabs_val > max_norm
233 22070 : IF (is_nonortho) THEN
234 10 : CALL cp_fm_release(aortho)
235 10 : EXIT loop
236 : END IF
237 : END DO
238 23986 : CALL cp_fm_release(aortho)
239 : END DO loop
240 :
241 5352 : CALL timestop(handle)
242 :
243 5352 : 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 8340 : 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 8340 : CALL timeset(routineN, handle)
294 :
295 8340 : nspins = SIZE(evects, 1)
296 8340 : nvects_total = SIZE(evects, 2)
297 8340 : 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 17918 : DO ispin = 1, nspins
306 17918 : CALL cp_fm_get_info(matrix=evects(ispin, 1), ncol_global=nactive(ispin))
307 : END DO
308 :
309 29246 : DO jvect = nvects_old + 1, nvects_total
310 : ! Orthogonalization <psi1_i | psi1_j>
311 187254 : DO ivect = 1, jvect - 1
312 166348 : CALL cp_fm_trace(evects(:, jvect), S_evects(:, ivect), weights(1:nspins), accurate=.FALSE.)
313 373064 : norm = SUM(weights(1:nspins))
314 :
315 393970 : DO ispin = 1, nspins
316 373064 : 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 45892 : DO ispin = 1, nspins
322 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, jvect), S_evects(ispin, jvect), &
323 45892 : ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
324 : END DO
325 :
326 20906 : CALL cp_fm_trace(evects(:, jvect), S_evects(:, jvect), weights(1:nspins), accurate=.FALSE.)
327 :
328 45892 : norm = SUM(weights(1:nspins))
329 20906 : norm = 1.0_dp/SQRT(norm)
330 :
331 54232 : DO ispin = 1, nspins
332 24986 : CALL cp_fm_scale(norm, evects(ispin, jvect))
333 45892 : CALL cp_fm_scale(norm, S_evects(ispin, jvect))
334 : END DO
335 : END DO
336 :
337 8340 : CALL timestop(handle)
338 :
339 8340 : END SUBROUTINE tddfpt_orthonormalize_psi1_psi1
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 6884 : SUBROUTINE tddfpt_compute_Aop_evects(Aop_evects, evects, S_evects, gs_mos, tddfpt_control, &
360 6884 : 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 6884 : CALL timeset(routineN, handle)
388 :
389 6884 : nspins = SIZE(gs_mos, 1)
390 6884 : nvects = SIZE(evects, 2)
391 6884 : do_hfx = tddfpt_control%do_hfx
392 6884 : do_admm = tddfpt_control%do_admm
393 6884 : IF (do_admm) THEN
394 872 : kernel_env_admm_aux => kernel_env%admm_kernel
395 : ELSE
396 6012 : NULLIFY (kernel_env_admm_aux)
397 : END IF
398 6884 : is_rks_triplets = tddfpt_control%rks_triplets
399 6884 : do_lri_response = tddfpt_control%do_lrigpw
400 6884 : do_bse = tddfpt_control%do_bse
401 6884 : do_bse_w_only = tddfpt_control%do_bse_w_only
402 6884 : do_bse_gw_only = tddfpt_control%do_bse_gw_only
403 6884 : IF (do_bse) do_hfx = .FALSE.
404 6884 : IF (do_bse_w_only) do_hfx = .FALSE.
405 6884 : 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 14954 : DO ispin = 1, nspins
417 6884 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
418 : END DO
419 :
420 6884 : IF (nvects > 0) THEN
421 6884 : CALL cp_fm_get_info(evects(1, 1), para_env=para_env)
422 6884 : IF (ALLOCATED(work_matrices%evects_sub)) THEN
423 24 : DO ivect = 1, nvects
424 40 : DO ispin = 1, SIZE(evects, 1)
425 16 : ASSOCIATE (evect => evects(ispin, ivect), work_matrix => work_matrices%evects_sub(ispin, ivect))
426 16 : IF (ASSOCIATED(evect%matrix_struct)) THEN
427 16 : IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
428 8 : CALL cp_fm_copy_general(evect, work_matrix, para_env)
429 : ELSE
430 8 : CALL cp_fm_copy_general(evect, fm_dummy, para_env)
431 : END IF
432 0 : ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
433 0 : CALL cp_fm_copy_general(fm_dummy, work_matrix, para_env)
434 : ELSE
435 0 : 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 6884 : 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 4344 : do_lri_response, tddfpt_mgrid=tddfpt_control%mgrid_is_explicit)
448 2540 : 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 2416 : kernel_env%stda_kernel, sub_env, work_matrices)
452 124 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
453 : ! No kernel
454 400 : DO ivect = 1, nvects
455 676 : DO ispin = 1, SIZE(evects, 1)
456 552 : CALL cp_fm_set_all(Aop_evects(ispin, ivect), 0.0_dp)
457 : END DO
458 : END DO
459 : ELSE
460 0 : CPABORT("Kernel type not implemented")
461 : END IF
462 :
463 : ! Add HFXLR before returning the subgroup kernel to the global matrices.
464 6884 : IF (tddfpt_control%kernel == tddfpt_kernel_full .AND. tddfpt_control%do_hfxlr) THEN
465 36 : rcut = tddfpt_control%hfxlr_rcut
466 36 : scale = tddfpt_control%hfxlr_scale
467 108 : DO ivect = 1, nvects
468 108 : IF (ALLOCATED(work_matrices%evects_sub)) THEN
469 0 : IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct)) THEN
470 : CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
471 : work_matrices%evects_sub(:, ivect), &
472 0 : work_matrices%Aop_evects_sub(:, ivect))
473 : END IF
474 : ELSE
475 : CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
476 72 : evects(:, ivect), Aop_evects(:, ivect))
477 : END IF
478 : END DO
479 : END IF
480 :
481 6884 : IF (ALLOCATED(work_matrices%evects_sub)) THEN
482 24 : DO ivect = 1, nvects
483 40 : DO ispin = 1, SIZE(evects, 1)
484 : ASSOCIATE (Aop_evect => Aop_evects(ispin, ivect), &
485 16 : work_matrix => work_matrices%Aop_evects_sub(ispin, ivect))
486 16 : IF (ASSOCIATED(Aop_evect%matrix_struct)) THEN
487 16 : IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
488 8 : CALL cp_fm_copy_general(work_matrix, Aop_evect, para_env)
489 : ELSE
490 8 : CALL cp_fm_copy_general(fm_dummy, Aop_evect, para_env)
491 : END IF
492 0 : ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
493 0 : CALL cp_fm_copy_general(work_matrix, fm_dummy, para_env)
494 : ELSE
495 0 : CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env)
496 : END IF
497 : END ASSOCIATE
498 : END DO
499 : END DO
500 : END IF
501 :
502 : ! orbital energy difference term
503 6884 : IF ((.NOT. do_bse) .AND. (.NOT. do_bse_gw_only)) THEN
504 : CALL tddfpt_apply_energy_diff(Aop_evects=Aop_evects, evects=evects, S_evects=S_evects, &
505 6864 : gs_mos=gs_mos, matrix_ks=matrix_ks, tddfpt_control=tddfpt_control)
506 : ELSE
507 : CALL zeroth_order_gw(qs_env=qs_env, Aop_evects=Aop_evects, evects=evects, S_evects=S_evects, &
508 20 : gs_mos=gs_mos, matrix_s=matrix_s, matrix_ks=matrix_ks)
509 : END IF
510 :
511 : ! if smeared occupation, then add aCCSX here
512 6884 : IF (tddfpt_control%do_smearing) THEN
513 24 : DO ivect = 1, nvects
514 36 : DO ispin = 1, nspins
515 : CALL add_smearing_aterm(qs_env, Aop_evects(ispin, ivect), evects(ispin, ivect), &
516 : S_evects(ispin, ivect), gs_mos(ispin)%mos_occ, &
517 24 : tddfpt_control%smeared_occup(ispin)%fermia, matrix_s)
518 : END DO
519 : END DO
520 : END IF
521 :
522 6884 : IF (do_hfx) THEN
523 1556 : IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
524 : ! full TDDFPT kernel
525 : CALL tddfpt_apply_hfx(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, do_admm=do_admm, &
526 : qs_env=qs_env, wfm_rho_orb=work_matrices%hfx_fm_ao_ao, &
527 : work_hmat_symm=work_matrices%hfx_hmat_symm, &
528 : work_hmat_asymm=work_matrices%hfx_hmat_asymm, &
529 : work_rho_ia_ao_symm=work_matrices%hfx_rho_ao_symm, &
530 1556 : work_rho_ia_ao_asymm=work_matrices%hfx_rho_ao_asymm)
531 0 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
532 : ! sTDA kernel
533 : ! special treatment of HFX term
534 0 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
535 : ! No kernel
536 : ! drop kernel contribution of HFX term
537 : ELSE
538 0 : CPABORT("Kernel type not implemented")
539 : END IF
540 : END IF
541 : ! short range HFX
542 6884 : IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
543 4344 : IF (tddfpt_control%do_hfxsr) THEN
544 22 : re_int = tddfpt_control%hfxsr_re_int
545 : ! symmetric dmat
546 : CALL tddfpt_apply_hfxsr_kernel(Aop_evects, evects, gs_mos, qs_env, &
547 : kernel_env%full_kernel%admm_env, &
548 : kernel_env%full_kernel%hfxsr_section, &
549 : kernel_env%full_kernel%x_data, 1, re_int, &
550 : work_rho_ia_ao=work_matrices%hfxsr_rho_ao_symm, &
551 : work_hmat=work_matrices%hfxsr_hmat_symm, &
552 22 : wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
553 : ! antisymmetric dmat
554 : CALL tddfpt_apply_hfxsr_kernel(Aop_evects, evects, gs_mos, qs_env, &
555 : kernel_env%full_kernel%admm_env, &
556 : kernel_env%full_kernel%hfxsr_section, &
557 : kernel_env%full_kernel%x_data, -1, .FALSE., &
558 : work_rho_ia_ao=work_matrices%hfxsr_rho_ao_asymm, &
559 : work_hmat=work_matrices%hfxsr_hmat_asymm, &
560 22 : wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
561 22 : tddfpt_control%hfxsr_re_int = .FALSE.
562 : END IF
563 : END IF
564 6884 : IF ((do_bse) .OR. (do_bse_w_only)) THEN
565 : ! add dynamical screening
566 : ! CALL tddfpt_apply_bse_debug(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, qs_env=qs_env)
567 : CALL tddfpt_apply_bse(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, qs_env=qs_env, &
568 36 : S_evects=S_evects)
569 : END IF
570 :
571 : END IF
572 :
573 6884 : CALL timestop(handle)
574 :
575 6884 : END SUBROUTINE tddfpt_compute_Aop_evects
576 :
577 : ! **************************************************************************************************
578 : !> \brief Solve eigenproblem for the reduced action matrix and find new Ritz eigenvectors and
579 : !> eigenvalues.
580 : !> \param ritz_vects Ritz eigenvectors (initialised on exit)
581 : !> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors (initialised on exit)
582 : !> \param evals Ritz eigenvalues (initialised on exit)
583 : !> \param krylov_vects Krylov's vectors
584 : !> \param Aop_krylov action of TDDFPT operator on Krylov's vectors
585 : !> \param Atilde TDDFPT matrix projected into the Krylov's vectors subspace
586 : !> \param nkvo ...
587 : !> \param nkvn ...
588 : !> \par History
589 : !> * 06.2016 created [Sergey Chulkov]
590 : !> * 03.2017 altered prototype, OpenMP parallelisation [Sergey Chulkov]
591 : ! **************************************************************************************************
592 6884 : SUBROUTINE tddfpt_compute_ritz_vects(ritz_vects, Aop_ritz, evals, krylov_vects, Aop_krylov, &
593 : Atilde, nkvo, nkvn)
594 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: ritz_vects, Aop_ritz
595 : REAL(kind=dp), DIMENSION(:), INTENT(out) :: evals
596 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: krylov_vects, Aop_krylov
597 : REAL(kind=dp), DIMENSION(:, :), POINTER :: Atilde
598 : INTEGER, INTENT(IN) :: nkvo, nkvn
599 :
600 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_ritz_vects'
601 :
602 : INTEGER :: handle, ikv, irv, ispin, nkvs, nrvs, &
603 : nspins
604 : REAL(kind=dp) :: act
605 6884 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: at12, at21, at22, evects_Atilde
606 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_global
607 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
608 : TYPE(cp_fm_type) :: Atilde_fm, evects_Atilde_fm
609 :
610 6884 : CALL timeset(routineN, handle)
611 :
612 6884 : nspins = SIZE(krylov_vects, 1)
613 6884 : nkvs = SIZE(krylov_vects, 2)
614 6884 : nrvs = SIZE(ritz_vects, 2)
615 6884 : CPASSERT(nkvs == nkvo + nkvn)
616 :
617 6884 : CALL cp_fm_get_info(krylov_vects(1, 1), context=blacs_env_global)
618 :
619 6884 : CALL cp_fm_struct_create(fm_struct, nrow_global=nkvs, ncol_global=nkvs, context=blacs_env_global)
620 6884 : CALL cp_fm_create(Atilde_fm, fm_struct, set_zero=.TRUE.)
621 6884 : CALL cp_fm_create(evects_Atilde_fm, fm_struct, set_zero=.TRUE.)
622 6884 : CALL cp_fm_struct_release(fm_struct)
623 :
624 : ! *** compute upper-diagonal reduced action matrix ***
625 6884 : CALL reallocate(Atilde, 1, nkvs, 1, nkvs)
626 : ! TO DO: the subroutine 'cp_fm_contracted_trace' will compute all elements of
627 : ! the matrix 'Atilde', however only upper-triangular elements are actually needed
628 : !
629 6884 : IF (nkvo == 0) THEN
630 : CALL cp_fm_contracted_trace(Aop_krylov(:, 1:nkvs), krylov_vects(:, 1:nkvs), &
631 1542 : Atilde(1:nkvs, 1:nkvs), accurate=.FALSE.)
632 : ELSE
633 48078 : ALLOCATE (at12(nkvn, nkvo), at21(nkvo, nkvn), at22(nkvn, nkvn))
634 : CALL cp_fm_contracted_trace(Aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, 1:nkvo), &
635 5342 : at12, accurate=.FALSE.)
636 181862 : Atilde(nkvo + 1:nkvs, 1:nkvo) = at12(1:nkvn, 1:nkvo)
637 : CALL cp_fm_contracted_trace(Aop_krylov(:, 1:nkvo), krylov_vects(:, nkvo + 1:nkvs), &
638 5342 : at21, accurate=.FALSE.)
639 153324 : Atilde(1:nkvo, nkvo + 1:nkvs) = at21(1:nkvo, 1:nkvn)
640 : CALL cp_fm_contracted_trace(Aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, nkvo + 1:nkvs), &
641 5342 : at22, accurate=.FALSE.)
642 67854 : Atilde(nkvo + 1:nkvs, nkvo + 1:nkvs) = at22(1:nkvn, 1:nkvn)
643 5342 : DEALLOCATE (at12, at21, at22)
644 : END IF
645 2105392 : Atilde = 0.5_dp*(Atilde + TRANSPOSE(Atilde))
646 6884 : CALL cp_fm_set_submatrix(Atilde_fm, Atilde)
647 :
648 : ! *** solve an eigenproblem for the reduced matrix ***
649 6884 : CALL choose_eigv_solver(Atilde_fm, evects_Atilde_fm, evals(1:nkvs))
650 :
651 27536 : ALLOCATE (evects_Atilde(nkvs, nrvs))
652 6884 : CALL cp_fm_get_submatrix(evects_Atilde_fm, evects_Atilde, start_row=1, start_col=1, n_rows=nkvs, n_cols=nrvs)
653 6884 : CALL cp_fm_release(evects_Atilde_fm)
654 6884 : CALL cp_fm_release(Atilde_fm)
655 :
656 : !$OMP PARALLEL DO DEFAULT(NONE), &
657 : !$OMP PRIVATE(act, ikv, irv, ispin), &
658 6884 : !$OMP SHARED(Aop_krylov, Aop_ritz, krylov_vects, evects_Atilde, nkvs, nrvs, nspins, ritz_vects)
659 : DO irv = 1, nrvs
660 : DO ispin = 1, nspins
661 : CALL cp_fm_set_all(ritz_vects(ispin, irv), 0.0_dp)
662 : CALL cp_fm_set_all(Aop_ritz(ispin, irv), 0.0_dp)
663 : END DO
664 :
665 : DO ikv = 1, nkvs
666 : act = evects_Atilde(ikv, irv)
667 : DO ispin = 1, nspins
668 : CALL cp_fm_scale_and_add(1.0_dp, ritz_vects(ispin, irv), &
669 : act, krylov_vects(ispin, ikv))
670 : CALL cp_fm_scale_and_add(1.0_dp, Aop_ritz(ispin, irv), &
671 : act, Aop_krylov(ispin, ikv))
672 : END DO
673 : END DO
674 : END DO
675 : !$OMP END PARALLEL DO
676 :
677 6884 : DEALLOCATE (evects_Atilde)
678 :
679 6884 : CALL timestop(handle)
680 :
681 13768 : END SUBROUTINE tddfpt_compute_ritz_vects
682 :
683 : ! **************************************************************************************************
684 : !> \brief Expand Krylov space by computing residual vectors.
685 : !> \param residual_vects residual vectors (modified on exit)
686 : !> \param evals Ritz eigenvalues (modified on exit)
687 : !> \param ritz_vects Ritz eigenvectors
688 : !> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors
689 : !> \param gs_mos molecular orbitals optimised for the ground state
690 : !> \param matrix_s overlap matrix
691 : !> \param tddfpt_control ...
692 : !> \par History
693 : !> * 06.2016 created [Sergey Chulkov]
694 : !> * 03.2017 refactored to achieve significant performance gain [Sergey Chulkov]
695 : ! **************************************************************************************************
696 5352 : SUBROUTINE tddfpt_compute_residual_vects(residual_vects, evals, ritz_vects, Aop_ritz, gs_mos, &
697 : matrix_s, tddfpt_control)
698 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: residual_vects
699 : REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
700 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: ritz_vects, Aop_ritz
701 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
702 : INTENT(in) :: gs_mos
703 : TYPE(dbcsr_type), POINTER :: matrix_s
704 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
705 :
706 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_residual_vects'
707 : REAL(kind=dp), PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*EPSILON(1.0_dp)
708 :
709 : INTEGER :: handle, ica, icb, icol_local, &
710 : irow_local, irv, ispin, nao, &
711 : ncols_local, nrows_local, nrvs, &
712 : nspins, spin2, spinflip
713 5352 : INTEGER, DIMENSION(:), POINTER :: col_indices_local, row_indices_local
714 : INTEGER, DIMENSION(maxspins) :: nactive, nmo_virt
715 : REAL(kind=dp) :: e_occ_plus_lambda, eref, lambda
716 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
717 5352 : POINTER :: weights_ldata
718 : TYPE(cp_fm_struct_type), POINTER :: ao_mo_struct, virt_mo_struct
719 5352 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: awork, vomat
720 :
721 5352 : CALL timeset(routineN, handle)
722 :
723 5352 : nspins = SIZE(residual_vects, 1)
724 5352 : nrvs = SIZE(residual_vects, 2)
725 5352 : spinflip = tddfpt_control%spinflip
726 :
727 5352 : IF (nrvs > 0) THEN
728 5352 : CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
729 33844 : ALLOCATE (awork(nspins), vomat(nspins))
730 11570 : DO ispin = 1, nspins
731 6218 : IF (spinflip == no_sf_tddfpt) THEN
732 : spin2 = ispin
733 : ELSE
734 106 : spin2 = 2
735 : END IF
736 6218 : nmo_virt(spin2) = SIZE(gs_mos(spin2)%evals_virt)
737 : !
738 : CALL cp_fm_get_info(matrix=ritz_vects(ispin, 1), matrix_struct=ao_mo_struct, &
739 6218 : ncol_global=nactive(ispin))
740 6218 : CALL cp_fm_create(awork(ispin), ao_mo_struct)
741 : !
742 : CALL cp_fm_struct_create(virt_mo_struct, nrow_global=nmo_virt(spin2), &
743 6218 : ncol_global=nactive(ispin), template_fmstruct=ao_mo_struct)
744 6218 : CALL cp_fm_create(vomat(ispin), virt_mo_struct)
745 11570 : CALL cp_fm_struct_release(virt_mo_struct)
746 : END DO
747 :
748 : ! *** actually compute residual vectors ***
749 18406 : DO irv = 1, nrvs
750 13054 : lambda = evals(irv)
751 :
752 34288 : DO ispin = 1, nspins
753 15882 : IF (spinflip == no_sf_tddfpt) THEN
754 : spin2 = ispin
755 : ELSE
756 498 : spin2 = 2
757 : END IF
758 : CALL cp_fm_get_info(vomat(ispin), nrow_local=nrows_local, &
759 : ncol_local=ncols_local, row_indices=row_indices_local, &
760 15882 : col_indices=col_indices_local, local_data=weights_ldata)
761 :
762 : ! awork := Ab(ispin, irv) - evals(irv) b(ispin, irv), where 'b' is a Ritz vector
763 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, ritz_vects(ispin, irv), awork(ispin), &
764 15882 : ncol=nactive(ispin), alpha=-lambda, beta=0.0_dp)
765 15882 : CALL cp_fm_scale_and_add(1.0_dp, awork(ispin), 1.0_dp, Aop_ritz(ispin, irv))
766 : !
767 : CALL parallel_gemm('T', 'N', nmo_virt(spin2), nactive(ispin), nao, 1.0_dp, gs_mos(spin2)%mos_virt, &
768 15882 : awork(ispin), 0.0_dp, vomat(ispin))
769 :
770 : ! Apply Davidson preconditioner to the residue vectors vomat to obtain new directions
771 127452 : DO icol_local = 1, ncols_local
772 111570 : ica = col_indices_local(icol_local)
773 111570 : icb = gs_mos(ispin)%index_active(ica)
774 111570 : e_occ_plus_lambda = gs_mos(ispin)%evals_occ(icb) + lambda
775 :
776 3562507 : DO irow_local = 1, nrows_local
777 3435055 : eref = gs_mos(spin2)%evals_virt(row_indices_local(irow_local)) - e_occ_plus_lambda
778 :
779 : ! eref = e_virt - e_occ - lambda = e_virt - e_occ - (eref_scale*lambda + (1-eref_scale)*lambda);
780 : ! eref_new = e_virt - e_occ - eref_scale*lambda = eref + (1 - eref_scale)*lambda
781 3435055 : IF (ABS(eref) < threshold) THEN
782 74 : eref = eref + (1.0_dp - eref_scale)*lambda
783 : END IF
784 :
785 3546625 : 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 44818 : vomat(ispin), 0.0_dp, residual_vects(ispin, irv))
791 : END DO
792 : END DO
793 : !
794 5352 : CALL cp_fm_release(awork)
795 10704 : CALL cp_fm_release(vomat)
796 : END IF
797 :
798 5352 : CALL timestop(handle)
799 :
800 10704 : 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 1542 : 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 1542 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_last
854 1542 : REAL(kind=dp), DIMENSION(:, :), POINTER :: Atilde
855 1542 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: Aop_krylov, Aop_ritz, krylov_vects, &
856 1542 : S_krylov
857 1542 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
858 :
859 1542 : CALL timeset(routineN, handle)
860 :
861 1542 : nspins = SIZE(evects, 1)
862 1542 : nstates = tddfpt_control%nstates
863 1542 : 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 1542 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
872 :
873 : ! adjust the number of Krylov vectors
874 1542 : max_krylov_vects = MIN(MAX(tddfpt_control%nkvs, nstates), INT(nstates_total))
875 :
876 15136 : ALLOCATE (Aop_ritz(nspins, nstates))
877 5708 : DO istate = 1, nstates
878 10510 : DO ispin = 1, nspins
879 8968 : 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 4626 : ALLOCATE (evals_last(max_krylov_vects))
884 : ALLOCATE (Aop_krylov(nspins, max_krylov_vects), krylov_vects(nspins, max_krylov_vects), &
885 1023666 : S_krylov(nspins, max_krylov_vects))
886 :
887 5708 : DO istate = 1, nstates
888 10510 : DO ispin = 1, nspins
889 4802 : CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, krylov_vects(ispin, istate))
890 4802 : CALL cp_fm_to_fm(evects(ispin, istate), krylov_vects(ispin, istate))
891 :
892 4802 : CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, S_krylov(ispin, istate))
893 4802 : CALL cp_fm_to_fm(S_evects(ispin, istate), S_krylov(ispin, istate))
894 :
895 8968 : 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 1542 : nvects_exist = 0
900 1542 : nvects_new = nstates
901 :
902 1542 : t1 = m_walltime()
903 :
904 1542 : ALLOCATE (Atilde(1, 1))
905 :
906 6884 : DO
907 : ! davidson iteration
908 6884 : 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 6884 : 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 6884 : 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 6884 : logger=logger, tddfpt_print_section=tddfpt_print_section)
929 :
930 24898 : conv = MAXVAL(ABS(evals_last(1:nstates) - evals(1:nstates)))
931 :
932 6884 : nvects_exist = nvects_exist + nvects_new
933 6884 : IF (nvects_exist + nvects_new > max_krylov_vects) THEN
934 506 : nvects_new = max_krylov_vects - nvects_exist
935 : END IF
936 6884 : IF (iter >= tddfpt_control%niters) nvects_new = 0
937 :
938 6884 : IF (conv > tddfpt_control%conv .AND. nvects_new > 0) THEN
939 : ! compute residual vectors for the next iteration
940 18406 : DO istate = 1, nvects_new
941 34288 : DO ispin = 1, nspins
942 : CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
943 15882 : krylov_vects(ispin, nvects_exist + istate))
944 : CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
945 15882 : S_krylov(ispin, nvects_exist + istate))
946 : CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
947 28936 : Aop_krylov(ispin, nvects_exist + istate))
948 : END DO
949 : END DO
950 :
951 : CALL tddfpt_compute_residual_vects(residual_vects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
952 : evals=evals_last(1:nvects_new), &
953 : ritz_vects=evects(:, 1:nvects_new), Aop_ritz=Aop_ritz(:, 1:nvects_new), &
954 5352 : gs_mos=gs_mos, matrix_s=matrix_s(1)%matrix, tddfpt_control=tddfpt_control)
955 :
956 : CALL tddfpt_orthogonalize_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
957 : work_matrices%S_C0_C0T, qs_env, &
958 5352 : gs_mos, evals(1:nstates), tddfpt_control, work_matrices%S_C0)
959 :
960 : CALL tddfpt_orthonormalize_psi1_psi1(krylov_vects(:, 1:nvects_exist + nvects_new), nvects_new, &
961 5352 : S_krylov(:, 1:nvects_exist + nvects_new), matrix_s(1)%matrix)
962 :
963 : is_nonortho = tddfpt_is_nonorthogonal_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
964 : work_matrices%S_C0, tddfpt_control%orthogonal_eps, &
965 5352 : tddfpt_control%spinflip)
966 : ELSE
967 : ! convergence or the maximum number of Krylov vectors have been achieved
968 1532 : nvects_new = 0
969 1532 : is_nonortho = .FALSE.
970 : END IF
971 :
972 6884 : t2 = m_walltime()
973 6884 : IF (energy_unit > 0) THEN
974 443 : WRITE (energy_unit, '(/,4X,A,T14,A,T36,A)') "State", "Exc. energy (eV)", "Convergence (eV)"
975 1159 : DO istate = 1, nstates
976 716 : WRITE (energy_unit, '(1X,I8,T12,F14.7,T38,ES11.4)') istate, &
977 1875 : evals_last(istate)*evolt, (evals_last(istate) - evals(istate))*evolt
978 : END DO
979 443 : WRITE (energy_unit, *)
980 443 : CALL m_flush(energy_unit)
981 : END IF
982 :
983 6884 : IF (iter_unit > 0) THEN
984 3442 : nstates_conv = 0
985 12449 : DO istate = 1, nstates
986 12449 : IF (ABS(evals_last(istate) - evals(istate)) <= tddfpt_control%conv) THEN
987 3099 : nstates_conv = nstates_conv + 1
988 : END IF
989 : END DO
990 :
991 3442 : WRITE (iter_unit, '(T7,I8,T24,F7.1,T40,ES11.4,T66,I8)') iter, t2 - t1, conv, nstates_conv
992 3442 : CALL m_flush(iter_unit)
993 : END IF
994 :
995 24898 : t1 = t2
996 24898 : evals(1:nstates) = evals_last(1:nstates)
997 :
998 : ! nvects_new == 0 if iter >= tddfpt_control%niters
999 6884 : IF (nvects_new == 0 .OR. is_nonortho) THEN
1000 : ! restart Davidson iterations
1001 : CALL tddfpt_orthogonalize_psi1_psi0(evects, work_matrices%S_C0_C0T, qs_env, &
1002 : gs_mos, &
1003 1542 : evals(1:nstates), tddfpt_control, work_matrices%S_C0)
1004 1542 : CALL tddfpt_orthonormalize_psi1_psi1(evects, nstates, S_evects, matrix_s(1)%matrix)
1005 :
1006 : EXIT
1007 : END IF
1008 : END DO
1009 :
1010 1542 : DEALLOCATE (Atilde)
1011 :
1012 18762 : DO istate = nvects_exist + nvects_new, 1, -1
1013 39446 : DO ispin = nspins, 1, -1
1014 20684 : CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, Aop_krylov(ispin, istate))
1015 20684 : CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, S_krylov(ispin, istate))
1016 37904 : CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, krylov_vects(ispin, istate))
1017 : END DO
1018 : END DO
1019 1542 : DEALLOCATE (Aop_krylov, krylov_vects, S_krylov)
1020 1542 : DEALLOCATE (evals_last)
1021 :
1022 5708 : DO istate = nstates, 1, -1
1023 10510 : DO ispin = nspins, 1, -1
1024 8968 : CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, Aop_ritz(ispin, istate))
1025 : END DO
1026 : END DO
1027 1542 : DEALLOCATE (Aop_ritz)
1028 :
1029 1542 : CALL timestop(handle)
1030 :
1031 1542 : END FUNCTION tddfpt_davidson_solver
1032 :
1033 : END MODULE qs_tddfpt2_eigensolver
|
