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 7774 : SUBROUTINE tddfpt_orthogonalize_psi1_psi0(evects, S_C0_C0T, qs_env, gs_mos, evals, &
106 7774 : 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 7774 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos
123 : TYPE(mp_para_env_type), POINTER :: para_env
124 :
125 7774 : CALL timeset(routineN, handle)
126 :
127 7774 : nspins = SIZE(evects, 1)
128 7774 : nvects = SIZE(evects, 2)
129 :
130 7774 : IF (nvects > 0) THEN
131 7774 : IF (.NOT. tddfpt_control%do_smearing) THEN
132 16360 : DO ispin = 1, nspins
133 8600 : IF (tddfpt_control%spinflip == no_sf_tddfpt) THEN
134 : spin = ispin
135 : ELSE
136 122 : spin = 2
137 : END IF
138 : CALL cp_fm_get_info(matrix=evects(ispin, 1), matrix_struct=matrix_struct, &
139 8600 : nrow_global=nao, ncol_global=nactive)
140 8600 : CALL cp_fm_create(evortho, matrix_struct)
141 30220 : 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 21620 : evects(ispin, ivect), 0.0_dp, evortho)
145 30220 : CALL cp_fm_scale_and_add(1.0_dp, evects(ispin, ivect), -1.0_dp, evortho)
146 : END DO
147 24960 : 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 7774 : CALL timestop(handle)
172 :
173 7774 : 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 4968 : 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 4968 : CALL timeset(routineN, handle)
206 :
207 4968 : nspins = SIZE(evects, 1)
208 4968 : nvects = SIZE(evects, 2)
209 :
210 4968 : is_nonortho = .FALSE.
211 :
212 10508 : loop: DO ispin = 1, nspins
213 :
214 5548 : IF (spinflip == no_sf_tddfpt) THEN
215 : spin = ispin
216 : ELSE
217 74 : spin = 2
218 : END IF
219 :
220 5548 : 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 5548 : nrow_global=nao, ncol_global=nactive)
223 : CALL cp_fm_struct_create(matrix_struct_tmp, nrow_global=nocc, &
224 5548 : ncol_global=nactive, template_fmstruct=matrix_struct)
225 5548 : CALL cp_fm_create(aortho, matrix_struct_tmp)
226 5548 : CALL cp_fm_struct_release(matrix_struct_tmp)
227 19072 : 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 13532 : evects(ispin, ivect), 0.0_dp, aortho)
231 13532 : CALL cp_fm_maxabsval(aortho, maxabs_val)
232 13532 : is_nonortho = maxabs_val > max_norm
233 19072 : IF (is_nonortho) THEN
234 8 : CALL cp_fm_release(aortho)
235 8 : EXIT loop
236 : END IF
237 : END DO
238 21596 : CALL cp_fm_release(aortho)
239 : END DO loop
240 :
241 4968 : CALL timestop(handle)
242 :
243 4968 : 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 7774 : 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 7774 : CALL timeset(routineN, handle)
294 :
295 7774 : nspins = SIZE(evects, 1)
296 7774 : nvects_total = SIZE(evects, 2)
297 7774 : 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 16388 : DO ispin = 1, nspins
306 16388 : CALL cp_fm_get_info(matrix=evects(ispin, 1), ncol_global=nactive(ispin))
307 : END DO
308 :
309 26832 : DO jvect = nvects_old + 1, nvects_total
310 : ! Orthogonalization <psi1_i | psi1_j>
311 166644 : DO ivect = 1, jvect - 1
312 147586 : CALL cp_fm_trace(evects(:, jvect), S_evects(:, ivect), weights(1:nspins), accurate=.FALSE.)
313 318186 : norm = SUM(weights(1:nspins))
314 :
315 337244 : DO ispin = 1, nspins
316 318186 : 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 40692 : DO ispin = 1, nspins
322 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, evects(ispin, jvect), S_evects(ispin, jvect), &
323 40692 : ncol=nactive(ispin), alpha=1.0_dp, beta=0.0_dp)
324 : END DO
325 :
326 19058 : CALL cp_fm_trace(evects(:, jvect), S_evects(:, jvect), weights(1:nspins), accurate=.FALSE.)
327 :
328 40692 : norm = SUM(weights(1:nspins))
329 19058 : norm = 1.0_dp/SQRT(norm)
330 :
331 48466 : DO ispin = 1, nspins
332 21634 : CALL cp_fm_scale(norm, evects(ispin, jvect))
333 40692 : CALL cp_fm_scale(norm, S_evects(ispin, jvect))
334 : END DO
335 : END DO
336 :
337 7774 : CALL timestop(handle)
338 :
339 7774 : 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 6408 : SUBROUTINE tddfpt_compute_Aop_evects(Aop_evects, evects, S_evects, gs_mos, tddfpt_control, &
360 6408 : 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 6408 : CALL timeset(routineN, handle)
388 :
389 6408 : nspins = SIZE(gs_mos, 1)
390 6408 : nvects = SIZE(evects, 2)
391 6408 : do_hfx = tddfpt_control%do_hfx
392 6408 : do_admm = tddfpt_control%do_admm
393 6408 : IF (do_admm) THEN
394 856 : kernel_env_admm_aux => kernel_env%admm_kernel
395 : ELSE
396 5552 : NULLIFY (kernel_env_admm_aux)
397 : END IF
398 6408 : is_rks_triplets = tddfpt_control%rks_triplets
399 6408 : do_lri_response = tddfpt_control%do_lrigpw
400 6408 : do_bse = tddfpt_control%do_bse
401 6408 : do_bse_w_only = tddfpt_control%do_bse_w_only
402 6408 : do_bse_gw_only = tddfpt_control%do_bse_gw_only
403 6408 : IF (do_bse) do_hfx = .FALSE.
404 6408 : IF (do_bse_w_only) do_hfx = .FALSE.
405 6408 : 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 13622 : DO ispin = 1, nspins
417 6408 : nmo_occ(ispin) = SIZE(gs_mos(ispin)%evals_occ)
418 : END DO
419 :
420 6408 : IF (nvects > 0) THEN
421 6408 : CALL cp_fm_get_info(evects(1, 1), para_env=para_env)
422 6408 : 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 6408 : 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 3980 : do_lri_response, tddfpt_mgrid=tddfpt_control%mgrid_is_explicit)
448 2428 : 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 2304 : 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 6408 : IF (ALLOCATED(work_matrices%evects_sub)) THEN
464 24 : DO ivect = 1, nvects
465 40 : DO ispin = 1, SIZE(evects, 1)
466 : ASSOCIATE (Aop_evect => Aop_evects(ispin, ivect), &
467 16 : work_matrix => work_matrices%Aop_evects_sub(ispin, ivect))
468 16 : IF (ASSOCIATED(Aop_evect%matrix_struct)) THEN
469 16 : IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
470 8 : CALL cp_fm_copy_general(work_matrix, Aop_evect, para_env)
471 : ELSE
472 8 : CALL cp_fm_copy_general(fm_dummy, Aop_evect, para_env)
473 : END IF
474 0 : ELSE IF (ASSOCIATED(work_matrix%matrix_struct)) THEN
475 0 : CALL cp_fm_copy_general(work_matrix, fm_dummy, para_env)
476 : ELSE
477 0 : CALL cp_fm_copy_general(fm_dummy, fm_dummy, para_env)
478 : END IF
479 : END ASSOCIATE
480 : END DO
481 : END DO
482 : END IF
483 :
484 : ! orbital energy difference term
485 6408 : IF ((.NOT. do_bse) .AND. (.NOT. do_bse_gw_only)) THEN
486 : CALL tddfpt_apply_energy_diff(Aop_evects=Aop_evects, evects=evects, S_evects=S_evects, &
487 6388 : gs_mos=gs_mos, matrix_ks=matrix_ks, tddfpt_control=tddfpt_control)
488 : ELSE
489 : CALL zeroth_order_gw(qs_env=qs_env, Aop_evects=Aop_evects, evects=evects, S_evects=S_evects, &
490 20 : gs_mos=gs_mos, matrix_s=matrix_s, matrix_ks=matrix_ks)
491 : END IF
492 :
493 : ! if smeared occupation, then add aCCSX here
494 6408 : IF (tddfpt_control%do_smearing) THEN
495 24 : DO ivect = 1, nvects
496 36 : DO ispin = 1, nspins
497 : CALL add_smearing_aterm(qs_env, Aop_evects(ispin, ivect), evects(ispin, ivect), &
498 : S_evects(ispin, ivect), gs_mos(ispin)%mos_occ, &
499 24 : tddfpt_control%smeared_occup(ispin)%fermia, matrix_s)
500 : END DO
501 : END DO
502 : END IF
503 :
504 6408 : IF (do_hfx) THEN
505 1546 : IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
506 : ! full TDDFPT kernel
507 : CALL tddfpt_apply_hfx(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, do_admm=do_admm, &
508 : qs_env=qs_env, wfm_rho_orb=work_matrices%hfx_fm_ao_ao, &
509 : work_hmat_symm=work_matrices%hfx_hmat_symm, &
510 : work_hmat_asymm=work_matrices%hfx_hmat_asymm, &
511 : work_rho_ia_ao_symm=work_matrices%hfx_rho_ao_symm, &
512 1546 : work_rho_ia_ao_asymm=work_matrices%hfx_rho_ao_asymm)
513 0 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_stda) THEN
514 : ! sTDA kernel
515 : ! special treatment of HFX term
516 0 : ELSE IF (tddfpt_control%kernel == tddfpt_kernel_none) THEN
517 : ! No kernel
518 : ! drop kernel contribution of HFX term
519 : ELSE
520 0 : CPABORT("Kernel type not implemented")
521 : END IF
522 : END IF
523 : ! short/long range HFX
524 6408 : IF (tddfpt_control%kernel == tddfpt_kernel_full) THEN
525 3980 : IF (tddfpt_control%do_hfxsr) THEN
526 22 : re_int = tddfpt_control%hfxsr_re_int
527 : ! symmetric dmat
528 : CALL tddfpt_apply_hfxsr_kernel(Aop_evects, evects, gs_mos, qs_env, &
529 : kernel_env%full_kernel%admm_env, &
530 : kernel_env%full_kernel%hfxsr_section, &
531 : kernel_env%full_kernel%x_data, 1, re_int, &
532 : work_rho_ia_ao=work_matrices%hfxsr_rho_ao_symm, &
533 : work_hmat=work_matrices%hfxsr_hmat_symm, &
534 22 : wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
535 : ! antisymmetric dmat
536 : CALL tddfpt_apply_hfxsr_kernel(Aop_evects, evects, gs_mos, qs_env, &
537 : kernel_env%full_kernel%admm_env, &
538 : kernel_env%full_kernel%hfxsr_section, &
539 : kernel_env%full_kernel%x_data, -1, .FALSE., &
540 : work_rho_ia_ao=work_matrices%hfxsr_rho_ao_asymm, &
541 : work_hmat=work_matrices%hfxsr_hmat_asymm, &
542 22 : wfm_rho_orb=work_matrices%hfxsr_fm_ao_ao)
543 22 : tddfpt_control%hfxsr_re_int = .FALSE.
544 : END IF
545 3980 : IF (tddfpt_control%do_hfxlr) THEN
546 36 : rcut = tddfpt_control%hfxlr_rcut
547 36 : scale = tddfpt_control%hfxlr_scale
548 108 : DO ivect = 1, nvects
549 108 : IF (ALLOCATED(work_matrices%evects_sub)) THEN
550 0 : IF (ASSOCIATED(work_matrices%evects_sub(1, ivect)%matrix_struct)) THEN
551 : CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
552 : work_matrices%evects_sub(:, ivect), &
553 0 : work_matrices%Aop_evects_sub(:, ivect))
554 : ELSE
555 : ! skip trial vectors which are assigned to different parallel groups
556 : CYCLE
557 : END IF
558 : ELSE
559 : CALL tddfpt_apply_hfxlr_kernel(qs_env, sub_env, rcut, scale, work_matrices, &
560 72 : evects(:, ivect), Aop_evects(:, ivect))
561 : END IF
562 : END DO
563 : END IF
564 : END IF
565 6408 : IF ((do_bse) .OR. (do_bse_w_only)) THEN
566 : ! add dynamical screening
567 : ! CALL tddfpt_apply_bse_debug(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, qs_env=qs_env)
568 : CALL tddfpt_apply_bse(Aop_evects=Aop_evects, evects=evects, gs_mos=gs_mos, qs_env=qs_env, &
569 36 : S_evects=S_evects)
570 : END IF
571 :
572 : END IF
573 :
574 6408 : CALL timestop(handle)
575 :
576 6408 : END SUBROUTINE tddfpt_compute_Aop_evects
577 :
578 : ! **************************************************************************************************
579 : !> \brief Solve eigenproblem for the reduced action matrix and find new Ritz eigenvectors and
580 : !> eigenvalues.
581 : !> \param ritz_vects Ritz eigenvectors (initialised on exit)
582 : !> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors (initialised on exit)
583 : !> \param evals Ritz eigenvalues (initialised on exit)
584 : !> \param krylov_vects Krylov's vectors
585 : !> \param Aop_krylov action of TDDFPT operator on Krylov's vectors
586 : !> \param Atilde TDDFPT matrix projected into the Krylov's vectors subspace
587 : !> \param nkvo ...
588 : !> \param nkvn ...
589 : !> \par History
590 : !> * 06.2016 created [Sergey Chulkov]
591 : !> * 03.2017 altered prototype, OpenMP parallelisation [Sergey Chulkov]
592 : ! **************************************************************************************************
593 6408 : SUBROUTINE tddfpt_compute_ritz_vects(ritz_vects, Aop_ritz, evals, krylov_vects, Aop_krylov, &
594 : Atilde, nkvo, nkvn)
595 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: ritz_vects, Aop_ritz
596 : REAL(kind=dp), DIMENSION(:), INTENT(out) :: evals
597 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: krylov_vects, Aop_krylov
598 : REAL(kind=dp), DIMENSION(:, :), POINTER :: Atilde
599 : INTEGER, INTENT(IN) :: nkvo, nkvn
600 :
601 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_ritz_vects'
602 :
603 : INTEGER :: handle, ikv, irv, ispin, nkvs, nrvs, &
604 : nspins
605 : REAL(kind=dp) :: act
606 6408 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: at12, at21, at22, evects_Atilde
607 : TYPE(cp_blacs_env_type), POINTER :: blacs_env_global
608 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
609 : TYPE(cp_fm_type) :: Atilde_fm, evects_Atilde_fm
610 :
611 6408 : CALL timeset(routineN, handle)
612 :
613 6408 : nspins = SIZE(krylov_vects, 1)
614 6408 : nkvs = SIZE(krylov_vects, 2)
615 6408 : nrvs = SIZE(ritz_vects, 2)
616 6408 : CPASSERT(nkvs == nkvo + nkvn)
617 :
618 6408 : CALL cp_fm_get_info(krylov_vects(1, 1), context=blacs_env_global)
619 :
620 6408 : CALL cp_fm_struct_create(fm_struct, nrow_global=nkvs, ncol_global=nkvs, context=blacs_env_global)
621 6408 : CALL cp_fm_create(Atilde_fm, fm_struct, set_zero=.TRUE.)
622 6408 : CALL cp_fm_create(evects_Atilde_fm, fm_struct, set_zero=.TRUE.)
623 6408 : CALL cp_fm_struct_release(fm_struct)
624 :
625 : ! *** compute upper-diagonal reduced action matrix ***
626 6408 : CALL reallocate(Atilde, 1, nkvs, 1, nkvs)
627 : ! TO DO: the subroutine 'cp_fm_contracted_trace' will compute all elements of
628 : ! the matrix 'Atilde', however only upper-triangular elements are actually needed
629 : !
630 6408 : IF (nkvo == 0) THEN
631 : CALL cp_fm_contracted_trace(Aop_krylov(:, 1:nkvs), krylov_vects(:, 1:nkvs), &
632 1448 : Atilde(1:nkvs, 1:nkvs), accurate=.FALSE.)
633 : ELSE
634 44640 : ALLOCATE (at12(nkvn, nkvo), at21(nkvo, nkvn), at22(nkvn, nkvn))
635 : CALL cp_fm_contracted_trace(Aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, 1:nkvo), &
636 4960 : at12, accurate=.FALSE.)
637 161840 : Atilde(nkvo + 1:nkvs, 1:nkvo) = at12(1:nkvn, 1:nkvo)
638 : CALL cp_fm_contracted_trace(Aop_krylov(:, 1:nkvo), krylov_vects(:, nkvo + 1:nkvs), &
639 4960 : at21, accurate=.FALSE.)
640 136000 : Atilde(1:nkvo, nkvo + 1:nkvs) = at21(1:nkvo, 1:nkvn)
641 : CALL cp_fm_contracted_trace(Aop_krylov(:, nkvo + 1:nkvs), krylov_vects(:, nkvo + 1:nkvs), &
642 4960 : at22, accurate=.FALSE.)
643 60648 : Atilde(nkvo + 1:nkvs, nkvo + 1:nkvs) = at22(1:nkvn, 1:nkvn)
644 4960 : DEALLOCATE (at12, at21, at22)
645 : END IF
646 1887104 : Atilde = 0.5_dp*(Atilde + TRANSPOSE(Atilde))
647 6408 : CALL cp_fm_set_submatrix(Atilde_fm, Atilde)
648 :
649 : ! *** solve an eigenproblem for the reduced matrix ***
650 6408 : CALL choose_eigv_solver(Atilde_fm, evects_Atilde_fm, evals(1:nkvs))
651 :
652 25632 : ALLOCATE (evects_Atilde(nkvs, nrvs))
653 6408 : CALL cp_fm_get_submatrix(evects_Atilde_fm, evects_Atilde, start_row=1, start_col=1, n_rows=nkvs, n_cols=nrvs)
654 6408 : CALL cp_fm_release(evects_Atilde_fm)
655 6408 : CALL cp_fm_release(Atilde_fm)
656 :
657 : !$OMP PARALLEL DO DEFAULT(NONE), &
658 : !$OMP PRIVATE(act, ikv, irv, ispin), &
659 6408 : !$OMP SHARED(Aop_krylov, Aop_ritz, krylov_vects, evects_Atilde, nkvs, nrvs, nspins, ritz_vects)
660 : DO irv = 1, nrvs
661 : DO ispin = 1, nspins
662 : CALL cp_fm_set_all(ritz_vects(ispin, irv), 0.0_dp)
663 : CALL cp_fm_set_all(Aop_ritz(ispin, irv), 0.0_dp)
664 : END DO
665 :
666 : DO ikv = 1, nkvs
667 : act = evects_Atilde(ikv, irv)
668 : DO ispin = 1, nspins
669 : CALL cp_fm_scale_and_add(1.0_dp, ritz_vects(ispin, irv), &
670 : act, krylov_vects(ispin, ikv))
671 : CALL cp_fm_scale_and_add(1.0_dp, Aop_ritz(ispin, irv), &
672 : act, Aop_krylov(ispin, ikv))
673 : END DO
674 : END DO
675 : END DO
676 : !$OMP END PARALLEL DO
677 :
678 6408 : DEALLOCATE (evects_Atilde)
679 :
680 6408 : CALL timestop(handle)
681 :
682 12816 : END SUBROUTINE tddfpt_compute_ritz_vects
683 :
684 : ! **************************************************************************************************
685 : !> \brief Expand Krylov space by computing residual vectors.
686 : !> \param residual_vects residual vectors (modified on exit)
687 : !> \param evals Ritz eigenvalues (modified on exit)
688 : !> \param ritz_vects Ritz eigenvectors
689 : !> \param Aop_ritz approximate action of TDDFPT operator on Ritz vectors
690 : !> \param gs_mos molecular orbitals optimised for the ground state
691 : !> \param matrix_s overlap matrix
692 : !> \param tddfpt_control ...
693 : !> \par History
694 : !> * 06.2016 created [Sergey Chulkov]
695 : !> * 03.2017 refactored to achieve significant performance gain [Sergey Chulkov]
696 : ! **************************************************************************************************
697 4968 : SUBROUTINE tddfpt_compute_residual_vects(residual_vects, evals, ritz_vects, Aop_ritz, gs_mos, &
698 : matrix_s, tddfpt_control)
699 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: residual_vects
700 : REAL(kind=dp), DIMENSION(:), INTENT(in) :: evals
701 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(in) :: ritz_vects, Aop_ritz
702 : TYPE(tddfpt_ground_state_mos), DIMENSION(:), &
703 : INTENT(in) :: gs_mos
704 : TYPE(dbcsr_type), POINTER :: matrix_s
705 : TYPE(tddfpt2_control_type), POINTER :: tddfpt_control
706 :
707 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_compute_residual_vects'
708 : REAL(kind=dp), PARAMETER :: eref_scale = 0.99_dp, threshold = 16.0_dp*EPSILON(1.0_dp)
709 :
710 : INTEGER :: handle, ica, icb, icol_local, &
711 : irow_local, irv, ispin, nao, &
712 : ncols_local, nrows_local, nrvs, &
713 : nspins, spin2, spinflip
714 4968 : INTEGER, DIMENSION(:), POINTER :: col_indices_local, row_indices_local
715 : INTEGER, DIMENSION(maxspins) :: nactive, nmo_virt
716 : REAL(kind=dp) :: e_occ_plus_lambda, eref, lambda
717 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :), &
718 4968 : POINTER :: weights_ldata
719 : TYPE(cp_fm_struct_type), POINTER :: ao_mo_struct, virt_mo_struct
720 4968 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: awork, vomat
721 :
722 4968 : CALL timeset(routineN, handle)
723 :
724 4968 : nspins = SIZE(residual_vects, 1)
725 4968 : nrvs = SIZE(residual_vects, 2)
726 4968 : spinflip = tddfpt_control%spinflip
727 :
728 4968 : IF (nrvs > 0) THEN
729 4968 : CALL dbcsr_get_info(matrix_s, nfullrows_total=nao)
730 35936 : ALLOCATE (awork(nspins), vomat(nspins))
731 10516 : DO ispin = 1, nspins
732 5548 : IF (spinflip == no_sf_tddfpt) THEN
733 : spin2 = ispin
734 : ELSE
735 74 : spin2 = 2
736 : END IF
737 5548 : nmo_virt(spin2) = SIZE(gs_mos(spin2)%evals_virt)
738 : !
739 : CALL cp_fm_get_info(matrix=ritz_vects(ispin, 1), matrix_struct=ao_mo_struct, &
740 5548 : ncol_global=nactive(ispin))
741 5548 : CALL cp_fm_create(awork(ispin), ao_mo_struct)
742 : !
743 : CALL cp_fm_struct_create(virt_mo_struct, nrow_global=nmo_virt(spin2), &
744 5548 : ncol_global=nactive(ispin), template_fmstruct=ao_mo_struct)
745 5548 : CALL cp_fm_create(vomat(ispin), virt_mo_struct)
746 10516 : CALL cp_fm_struct_release(virt_mo_struct)
747 : END DO
748 :
749 : ! *** actually compute residual vectors ***
750 16756 : DO irv = 1, nrvs
751 11788 : lambda = evals(irv)
752 :
753 30312 : DO ispin = 1, nspins
754 13556 : IF (spinflip == no_sf_tddfpt) THEN
755 : spin2 = ispin
756 : ELSE
757 338 : spin2 = 2
758 : END IF
759 : CALL cp_fm_get_info(vomat(ispin), nrow_local=nrows_local, &
760 : ncol_local=ncols_local, row_indices=row_indices_local, &
761 13556 : col_indices=col_indices_local, local_data=weights_ldata)
762 :
763 : ! awork := Ab(ispin, irv) - evals(irv) b(ispin, irv), where 'b' is a Ritz vector
764 : CALL cp_dbcsr_sm_fm_multiply(matrix_s, ritz_vects(ispin, irv), awork(ispin), &
765 13556 : ncol=nactive(ispin), alpha=-lambda, beta=0.0_dp)
766 13556 : CALL cp_fm_scale_and_add(1.0_dp, awork(ispin), 1.0_dp, Aop_ritz(ispin, irv))
767 : !
768 : CALL parallel_gemm('T', 'N', nmo_virt(spin2), nactive(ispin), nao, 1.0_dp, gs_mos(spin2)%mos_virt, &
769 13556 : awork(ispin), 0.0_dp, vomat(ispin))
770 :
771 : ! Apply Davidson preconditioner to the residue vectors vomat to obtain new directions
772 115124 : DO icol_local = 1, ncols_local
773 101568 : ica = col_indices_local(icol_local)
774 101568 : icb = gs_mos(ispin)%index_active(ica)
775 101568 : e_occ_plus_lambda = gs_mos(ispin)%evals_occ(icb) + lambda
776 :
777 3487703 : DO irow_local = 1, nrows_local
778 3372579 : eref = gs_mos(spin2)%evals_virt(row_indices_local(irow_local)) - e_occ_plus_lambda
779 :
780 : ! eref = e_virt - e_occ - lambda = e_virt - e_occ - (eref_scale*lambda + (1-eref_scale)*lambda);
781 : ! eref_new = e_virt - e_occ - eref_scale*lambda = eref + (1 - eref_scale)*lambda
782 3372579 : IF (ABS(eref) < threshold) &
783 74 : eref = eref + (1.0_dp - eref_scale)*lambda
784 :
785 3474147 : 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 38900 : vomat(ispin), 0.0_dp, residual_vects(ispin, irv))
791 : END DO
792 : END DO
793 : !
794 4968 : CALL cp_fm_release(awork)
795 9936 : CALL cp_fm_release(vomat)
796 : END IF
797 :
798 4968 : CALL timestop(handle)
799 :
800 9936 : 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 1448 : 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 1448 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_last
854 1448 : REAL(kind=dp), DIMENSION(:, :), POINTER :: Atilde
855 1448 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: Aop_krylov, Aop_ritz, krylov_vects, &
856 1448 : S_krylov
857 1448 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
858 :
859 1448 : CALL timeset(routineN, handle)
860 :
861 1448 : nspins = SIZE(evects, 1)
862 1448 : nstates = tddfpt_control%nstates
863 1448 : 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 1448 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
872 :
873 : ! adjust the number of Krylov vectors
874 1448 : max_krylov_vects = MIN(MAX(tddfpt_control%nkvs, nstates), INT(nstates_total))
875 :
876 13920 : ALLOCATE (Aop_ritz(nspins, nstates))
877 5310 : DO istate = 1, nstates
878 9576 : DO ispin = 1, nspins
879 8128 : 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 4344 : ALLOCATE (evals_last(max_krylov_vects))
884 : ALLOCATE (Aop_krylov(nspins, max_krylov_vects), krylov_vects(nspins, max_krylov_vects), &
885 976210 : S_krylov(nspins, max_krylov_vects))
886 :
887 5310 : DO istate = 1, nstates
888 9576 : DO ispin = 1, nspins
889 4266 : CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, krylov_vects(ispin, istate))
890 4266 : CALL cp_fm_to_fm(evects(ispin, istate), krylov_vects(ispin, istate))
891 :
892 4266 : CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, S_krylov(ispin, istate))
893 4266 : CALL cp_fm_to_fm(S_evects(ispin, istate), S_krylov(ispin, istate))
894 :
895 8128 : 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 1448 : nvects_exist = 0
900 1448 : nvects_new = nstates
901 :
902 1448 : t1 = m_walltime()
903 :
904 1448 : ALLOCATE (Atilde(1, 1))
905 :
906 6408 : DO
907 : ! davidson iteration
908 6408 : 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 6408 : 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 6408 : 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 6408 : logger=logger, tddfpt_print_section=tddfpt_print_section)
929 :
930 29190 : conv = MAXVAL(ABS(evals_last(1:nstates) - evals(1:nstates)))
931 :
932 6408 : nvects_exist = nvects_exist + nvects_new
933 6408 : IF (nvects_exist + nvects_new > max_krylov_vects) &
934 454 : nvects_new = max_krylov_vects - nvects_exist
935 6408 : IF (iter >= tddfpt_control%niters) nvects_new = 0
936 :
937 6408 : IF (conv > tddfpt_control%conv .AND. nvects_new > 0) THEN
938 : ! compute residual vectors for the next iteration
939 16756 : DO istate = 1, nvects_new
940 30312 : DO ispin = 1, nspins
941 : CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
942 13556 : krylov_vects(ispin, nvects_exist + istate))
943 : CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
944 13556 : S_krylov(ispin, nvects_exist + istate))
945 : CALL fm_pool_create_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, &
946 25344 : Aop_krylov(ispin, nvects_exist + istate))
947 : END DO
948 : END DO
949 :
950 : CALL tddfpt_compute_residual_vects(residual_vects=krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
951 : evals=evals_last(1:nvects_new), &
952 : ritz_vects=evects(:, 1:nvects_new), Aop_ritz=Aop_ritz(:, 1:nvects_new), &
953 4968 : gs_mos=gs_mos, matrix_s=matrix_s(1)%matrix, tddfpt_control=tddfpt_control)
954 :
955 : CALL tddfpt_orthogonalize_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
956 : work_matrices%S_C0_C0T, qs_env, &
957 4968 : gs_mos, evals(1:nstates), tddfpt_control, work_matrices%S_C0)
958 :
959 : CALL tddfpt_orthonormalize_psi1_psi1(krylov_vects(:, 1:nvects_exist + nvects_new), nvects_new, &
960 4968 : S_krylov(:, 1:nvects_exist + nvects_new), matrix_s(1)%matrix)
961 :
962 : is_nonortho = tddfpt_is_nonorthogonal_psi1_psi0(krylov_vects(:, nvects_exist + 1:nvects_exist + nvects_new), &
963 : work_matrices%S_C0, tddfpt_control%orthogonal_eps, &
964 4968 : tddfpt_control%spinflip)
965 : ELSE
966 : ! convergence or the maximum number of Krylov vectors have been achieved
967 1440 : nvects_new = 0
968 1440 : is_nonortho = .FALSE.
969 : END IF
970 :
971 6408 : t2 = m_walltime()
972 6408 : IF (energy_unit > 0) THEN
973 486 : WRITE (energy_unit, '(/,4X,A,T14,A,T36,A)') "State", "Exc. energy (eV)", "Convergence (eV)"
974 1243 : DO istate = 1, nstates
975 757 : WRITE (energy_unit, '(1X,I8,T12,F14.7,T38,ES11.4)') istate, &
976 2000 : evals_last(istate)*evolt, (evals_last(istate) - evals(istate))*evolt
977 : END DO
978 486 : WRITE (energy_unit, *)
979 486 : CALL m_flush(energy_unit)
980 : END IF
981 :
982 6408 : IF (iter_unit > 0) THEN
983 3204 : nstates_conv = 0
984 11391 : DO istate = 1, nstates
985 8187 : IF (ABS(evals_last(istate) - evals(istate)) <= tddfpt_control%conv) &
986 5903 : nstates_conv = nstates_conv + 1
987 : END DO
988 :
989 3204 : WRITE (iter_unit, '(T7,I8,T24,F7.1,T40,ES11.4,T66,I8)') iter, t2 - t1, conv, nstates_conv
990 3204 : CALL m_flush(iter_unit)
991 : END IF
992 :
993 22782 : t1 = t2
994 22782 : evals(1:nstates) = evals_last(1:nstates)
995 :
996 : ! nvects_new == 0 if iter >= tddfpt_control%niters
997 6408 : IF (nvects_new == 0 .OR. is_nonortho) THEN
998 : ! restart Davidson iterations
999 : CALL tddfpt_orthogonalize_psi1_psi0(evects, work_matrices%S_C0_C0T, qs_env, &
1000 : gs_mos, &
1001 1448 : evals(1:nstates), tddfpt_control, work_matrices%S_C0)
1002 1448 : CALL tddfpt_orthonormalize_psi1_psi1(evects, nstates, S_evects, matrix_s(1)%matrix)
1003 :
1004 : EXIT
1005 : END IF
1006 : END DO
1007 :
1008 1448 : DEALLOCATE (Atilde)
1009 :
1010 17098 : DO istate = nvects_exist + nvects_new, 1, -1
1011 34920 : DO ispin = nspins, 1, -1
1012 17822 : CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, Aop_krylov(ispin, istate))
1013 17822 : CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, S_krylov(ispin, istate))
1014 33472 : CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, krylov_vects(ispin, istate))
1015 : END DO
1016 : END DO
1017 1448 : DEALLOCATE (Aop_krylov, krylov_vects, S_krylov)
1018 1448 : DEALLOCATE (evals_last)
1019 :
1020 5310 : DO istate = nstates, 1, -1
1021 9576 : DO ispin = nspins, 1, -1
1022 8128 : CALL fm_pool_give_back_fm(work_matrices%fm_pool_ao_mo_active(ispin)%pool, Aop_ritz(ispin, istate))
1023 : END DO
1024 : END DO
1025 1448 : DEALLOCATE (Aop_ritz)
1026 :
1027 1448 : CALL timestop(handle)
1028 :
1029 1448 : END FUNCTION tddfpt_davidson_solver
1030 :
1031 : END MODULE qs_tddfpt2_eigensolver
|
