Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : MODULE qs_tddfpt_eigensolver
10 : USE cp_blacs_env, ONLY: cp_blacs_env_type
11 : USE cp_control_types, ONLY: dft_control_type,&
12 : tddfpt_control_type
13 : USE cp_dbcsr_api, ONLY: dbcsr_p_type,&
14 : dbcsr_set
15 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,&
16 : cp_dbcsr_sm_fm_multiply
17 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
18 : cp_fm_symm,&
19 : cp_fm_trace
20 : USE cp_fm_diag, ONLY: cp_fm_syevd
21 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
22 : fm_pools_create_fm_vect,&
23 : fm_pools_give_back_fm_vect
24 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
25 : cp_fm_struct_p_type,&
26 : cp_fm_struct_release,&
27 : cp_fm_struct_type
28 : USE cp_fm_types, ONLY: cp_fm_create,&
29 : cp_fm_get_element,&
30 : cp_fm_release,&
31 : cp_fm_set_all,&
32 : cp_fm_set_element,&
33 : cp_fm_to_fm,&
34 : cp_fm_type
35 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,&
36 : cp_to_string
37 : USE input_constants, ONLY: tddfpt_davidson,&
38 : tddfpt_lanczos
39 : USE kinds, ONLY: default_string_length,&
40 : dp
41 : USE message_passing, ONLY: mp_para_env_type
42 : USE physcon, ONLY: evolt
43 : USE qs_environment_types, ONLY: get_qs_env,&
44 : qs_environment_type
45 : USE qs_matrix_pools, ONLY: mpools_get
46 : USE qs_p_env_methods, ONLY: p_op_l1,&
47 : p_op_l2,&
48 : p_postortho,&
49 : p_preortho
50 : USE qs_p_env_types, ONLY: qs_p_env_type
51 : USE qs_tddfpt_types, ONLY: tddfpt_env_type
52 : USE qs_tddfpt_utils, ONLY: co_initial_guess,&
53 : normalize,&
54 : reorthogonalize
55 : #include "./base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 :
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt_eigensolver'
60 :
61 : PRIVATE
62 :
63 : PUBLIC :: eigensolver
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief ...
69 : !> \param p_env ...
70 : !> \param qs_env ...
71 : !> \param t_env ...
72 : ! **************************************************************************************************
73 12 : SUBROUTINE eigensolver(p_env, qs_env, t_env)
74 :
75 : TYPE(qs_p_env_type) :: p_env
76 : TYPE(qs_environment_type), POINTER :: qs_env
77 : TYPE(tddfpt_env_type), INTENT(INOUT) :: t_env
78 :
79 : CHARACTER(len=*), PARAMETER :: routineN = 'eigensolver'
80 :
81 : INTEGER :: handle, n_ev, nspins, output_unit, &
82 : restarts
83 : LOGICAL :: do_kernel_save
84 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ievals
85 : TYPE(dft_control_type), POINTER :: dft_control
86 :
87 12 : CALL timeset(routineN, handle)
88 :
89 12 : NULLIFY (dft_control)
90 :
91 12 : output_unit = cp_logger_get_default_io_unit()
92 :
93 12 : CALL get_qs_env(qs_env, dft_control=dft_control)
94 12 : n_ev = dft_control%tddfpt_control%n_ev
95 12 : nspins = dft_control%nspins
96 :
97 36 : ALLOCATE (ievals(n_ev))
98 :
99 : !---------------!
100 : ! initial guess !
101 : !---------------!
102 12 : do_kernel_save = dft_control%tddfpt_control%do_kernel
103 12 : dft_control%tddfpt_control%do_kernel = .FALSE.
104 12 : IF (output_unit > 0) THEN
105 6 : WRITE (output_unit, *) " Generating initial guess"
106 6 : WRITE (output_unit, *)
107 : END IF
108 12 : IF (ASSOCIATED(dft_control%tddfpt_control%lumos)) THEN
109 12 : CALL co_initial_guess(t_env%evecs, ievals, n_ev, qs_env)
110 : ELSE
111 0 : IF (output_unit > 0) WRITE (output_unit, *) "LUMOS are needed in TDDFPT!"
112 0 : CPABORT("")
113 : END IF
114 12 : DO restarts = 1, dft_control%tddfpt_control%n_restarts
115 12 : IF (iterative_solver(ievals, t_env, p_env, qs_env, ievals)) EXIT
116 12 : IF (output_unit > 0) THEN
117 0 : WRITE (output_unit, *) " Restarting"
118 0 : WRITE (output_unit, *)
119 : END IF
120 : END DO
121 12 : dft_control%tddfpt_control%do_kernel = do_kernel_save
122 :
123 : !-----------------!
124 : ! call the solver !
125 : !-----------------!
126 12 : IF (output_unit > 0) THEN
127 6 : WRITE (output_unit, *)
128 6 : WRITE (output_unit, *) " Doing TDDFPT calculation"
129 6 : WRITE (output_unit, *)
130 : END IF
131 36 : DO restarts = 1, dft_control%tddfpt_control%n_restarts
132 24 : IF (iterative_solver(ievals, t_env, p_env, qs_env, t_env%evals)) EXIT
133 36 : IF (output_unit > 0) THEN
134 12 : WRITE (output_unit, *) " Restarting"
135 12 : WRITE (output_unit, *)
136 : END IF
137 : END DO
138 :
139 : !---------!
140 : ! cleanup !
141 : !---------!
142 12 : DEALLOCATE (ievals)
143 :
144 12 : CALL timestop(handle)
145 :
146 12 : END SUBROUTINE eigensolver
147 :
148 : ! in_evals : approximations to the eigenvalues for the preconditioner
149 : ! t_env : TD-DFT environment values
150 : ! p_env : perturbation environment values
151 : ! qs_env : general Quickstep environment values
152 : ! out_evals : the resulting eigenvalues
153 : ! error : used for error handling
154 : !
155 : ! res : the function will return wheter the eigenvalues are converged or not
156 :
157 : ! **************************************************************************************************
158 : !> \brief ...
159 : !> \param in_evals ...
160 : !> \param t_env ...
161 : !> \param p_env ...
162 : !> \param qs_env ...
163 : !> \param out_evals ...
164 : !> \return ...
165 : ! **************************************************************************************************
166 36 : FUNCTION iterative_solver(in_evals, &
167 : t_env, p_env, qs_env, &
168 36 : out_evals) RESULT(res)
169 :
170 : REAL(KIND=dp), DIMENSION(:) :: in_evals
171 : TYPE(tddfpt_env_type), INTENT(INOUT) :: t_env
172 : TYPE(qs_p_env_type) :: p_env
173 : TYPE(qs_environment_type), POINTER :: qs_env
174 : REAL(kind=dp), DIMENSION(:), OPTIONAL :: out_evals
175 : LOGICAL :: res
176 :
177 : CHARACTER(len=*), PARAMETER :: routineN = 'iterative_solver', &
178 : routineP = moduleN//':'//routineN
179 :
180 : CHARACTER :: mode
181 : INTEGER :: col, handle, i, iev, iter, j, k, &
182 : max_krylovspace_dim, max_kv, n_ev, &
183 : n_kv, nspins, output_unit, row, spin
184 36 : INTEGER, ALLOCATABLE, DIMENSION(:) :: must_improve
185 : REAL(dp) :: Atilde_ij, convergence, tmp, tmp2
186 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals_difference, evals_tmp
187 36 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: evals
188 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
189 36 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools
190 : TYPE(cp_fm_struct_p_type), ALLOCATABLE, &
191 36 : DIMENSION(:) :: kv_fm_struct
192 : TYPE(cp_fm_struct_type), POINTER :: tilde_fm_struct
193 : TYPE(cp_fm_type) :: Atilde, Us
194 36 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: R, X
195 36 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: Ab, b, Sb
196 36 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
197 : TYPE(dft_control_type), POINTER :: dft_control
198 : TYPE(mp_para_env_type), POINTER :: para_env
199 : TYPE(tddfpt_control_type), POINTER :: tddfpt_control
200 :
201 36 : res = .FALSE.
202 :
203 36 : CALL timeset(routineN, handle)
204 :
205 36 : NULLIFY (ao_mo_fm_pools, tddfpt_control, &
206 36 : tilde_fm_struct, matrix_s, dft_control, &
207 36 : para_env, blacs_env)
208 :
209 : CALL get_qs_env(qs_env, &
210 : matrix_s=matrix_s, &
211 : dft_control=dft_control, &
212 : para_env=para_env, &
213 36 : blacs_env=blacs_env)
214 :
215 36 : tddfpt_control => dft_control%tddfpt_control
216 36 : output_unit = cp_logger_get_default_io_unit()
217 36 : n_ev = tddfpt_control%n_ev
218 36 : nspins = dft_control%nspins
219 :
220 36 : IF (dft_control%tddfpt_control%diag_method == tddfpt_lanczos) THEN
221 : mode = 'L'
222 36 : ELSE IF (dft_control%tddfpt_control%diag_method == tddfpt_davidson) THEN
223 36 : mode = 'D'
224 : END IF
225 :
226 : !-----------------------------------------!
227 : ! determine the size of the problem !
228 : ! and how many krylov space vetors to use !
229 : !-----------------------------------------!
230 84 : max_krylovspace_dim = SUM(p_env%n_ao(1:nspins)*p_env%n_mo(1:nspins))
231 36 : max_kv = tddfpt_control%max_kv
232 36 : IF (max_krylovspace_dim <= max_kv) THEN
233 6 : max_kv = max_krylovspace_dim
234 6 : IF (output_unit > 0) THEN
235 3 : WRITE (output_unit, *) " Setting the maximum number of krylov vectors to ", max_kv, "!"
236 : END IF
237 : END IF
238 :
239 : !----------------------!
240 : ! allocate the vectors !
241 : !----------------------!
242 36 : CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools)
243 36 : CALL fm_pools_create_fm_vect(ao_mo_fm_pools, X, name=routineP//":X")
244 36 : CALL fm_pools_create_fm_vect(ao_mo_fm_pools, R, name=routineP//":R")
245 :
246 108 : ALLOCATE (evals_difference(n_ev))
247 :
248 108 : ALLOCATE (must_improve(n_ev))
249 :
250 144 : ALLOCATE (evals(max_kv, 0:max_kv))
251 108 : ALLOCATE (evals_tmp(max_kv))
252 :
253 0 : ALLOCATE (b(max_kv, nspins), Ab(max_kv, nspins), &
254 3024 : Sb(max_kv, nspins))
255 :
256 156 : ALLOCATE (kv_fm_struct(nspins))
257 :
258 84 : DO spin = 1, nspins
259 : CALL cp_fm_struct_create(kv_fm_struct(spin)%struct, para_env, blacs_env, &
260 84 : p_env%n_ao(spin), p_env%n_mo(spin))
261 : END DO
262 :
263 36 : IF (output_unit > 0) THEN
264 : WRITE (output_unit, '(2X,A,T69,A)') &
265 18 : "nvec", "Convergence"
266 : WRITE (output_unit, '(2X,A)') &
267 18 : "-----------------------------------------------------------------------------"
268 : END IF
269 :
270 36 : iter = 1
271 36 : k = 0
272 36 : n_kv = n_ev
273 496 : iteration: DO
274 :
275 248 : CALL allocate_krylov_vectors(b, "b-", k + 1, n_kv, nspins, kv_fm_struct)
276 248 : CALL allocate_krylov_vectors(Ab, "Ab-", k + 1, n_kv, nspins, kv_fm_struct)
277 248 : CALL allocate_krylov_vectors(Sb, "Sb-", k + 1, n_kv, nspins, kv_fm_struct)
278 :
279 794 : DO i = 1, n_kv
280 546 : k = k + 1
281 :
282 546 : IF (k <= SIZE(t_env%evecs, 1)) THEN ! the first iteration
283 :
284 : ! take the initial guess
285 228 : DO spin = 1, nspins
286 228 : CALL cp_fm_to_fm(t_env%evecs(k, spin), b(k, spin))
287 : END DO
288 :
289 : ELSE
290 :
291 : ! create a new vector
292 450 : IF (mode == 'L') THEN
293 :
294 0 : DO spin = 1, nspins
295 0 : IF (tddfpt_control%invert_S) THEN
296 : CALL cp_fm_symm('L', 'U', p_env%n_ao(spin), p_env%n_mo(spin), &
297 : 1.0_dp, t_env%invS(spin), Ab(k - 1, spin), &
298 0 : 0.0_dp, b(k, spin))
299 : ELSE
300 0 : CALL cp_fm_to_fm(Ab(k - 1, spin), b(k, spin))
301 : END IF
302 : END DO
303 :
304 450 : ELSE IF (mode == 'D') THEN
305 :
306 450 : iev = must_improve(i)
307 : ! create the new davidson vector
308 960 : DO spin = 1, nspins
309 :
310 510 : CALL cp_fm_set_all(R(spin), 0.0_dp)
311 9312 : DO j = 1, k - i
312 8802 : CALL cp_fm_to_fm(Ab(j, spin), X(spin))
313 : CALL cp_fm_scale_and_add(1.0_dp, X(spin), &
314 8802 : -evals(iev, iter - 1), Sb(j, spin))
315 8802 : CALL cp_fm_get_element(Us, j, iev, tmp)
316 : CALL cp_fm_scale_and_add(1.0_dp, R(spin), &
317 9312 : tmp, X(spin))
318 : END DO
319 :
320 510 : IF (tddfpt_control%invert_S) THEN
321 : CALL cp_fm_symm('L', 'U', p_env%n_ao(spin), p_env%n_mo(spin), &
322 : 1.0_dp, t_env%invS(spin), R(spin), &
323 510 : 0.0_dp, X(spin))
324 : ELSE
325 0 : CALL cp_fm_to_fm(R(spin), X(spin))
326 : END IF
327 :
328 : !----------------!
329 : ! preconditioner !
330 : !----------------!
331 960 : IF (dft_control%tddfpt_control%precond) THEN
332 900 : DO col = 1, p_env%n_mo(spin)
333 720 : IF (col <= n_ev) THEN
334 540 : tmp2 = ABS(evals(iev, iter - 1) - in_evals(col))
335 : ELSE
336 180 : tmp2 = ABS(evals(iev, iter - 1) - (in_evals(n_ev) + 10.0_dp))
337 : END IF
338 : ! protect against division by 0 by a introducing a cutoff.
339 720 : tmp2 = MAX(tmp2, 100*EPSILON(1.0_dp))
340 17460 : DO row = 1, p_env%n_ao(spin)
341 16560 : CALL cp_fm_get_element(X(spin), row, col, tmp)
342 17280 : CALL cp_fm_set_element(b(k, spin), row, col, tmp/tmp2)
343 : END DO
344 : END DO
345 : ELSE
346 330 : CALL cp_fm_to_fm(X(spin), b(k, spin))
347 : END IF
348 :
349 : END DO
350 :
351 : ELSE
352 0 : IF (output_unit > 0) WRITE (output_unit, *) "unknown mode"
353 0 : CPABORT("")
354 : END IF
355 :
356 : END IF
357 :
358 546 : CALL p_preortho(p_env, qs_env, b(k, :))
359 1638 : DO j = 1, tddfpt_control%n_reortho
360 1638 : CALL reorthogonalize(b(k, :), b, Sb, R, k - 1) ! R is temp
361 : END DO
362 546 : CALL normalize(b(k, :), R, matrix_s) ! R is temp
363 1188 : DO spin = 1, nspins
364 1188 : CALL cp_fm_to_fm(b(k, spin), X(spin))
365 : END DO
366 : CALL apply_op(X, Ab(k, :), p_env, qs_env, &
367 546 : dft_control%tddfpt_control%do_kernel)
368 546 : CALL p_postortho(p_env, qs_env, Ab(k, :))
369 1436 : DO spin = 1, nspins
370 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, &
371 : b(k, spin), &
372 : Sb(k, spin), &
373 1188 : p_env%n_mo(spin))
374 : END DO
375 : END DO
376 :
377 : !--------------------------------------------!
378 : ! deallocate memory for the reduced matrices !
379 : !--------------------------------------------!
380 248 : CALL cp_fm_release(Atilde)
381 248 : CALL cp_fm_release(Us)
382 248 : IF (ASSOCIATED(tilde_fm_struct)) CALL cp_fm_struct_release(tilde_fm_struct)
383 :
384 : !------------------------------------------!
385 : ! allocate memory for the reduced matrices !
386 : !------------------------------------------!
387 248 : CALL cp_fm_struct_create(tilde_fm_struct, para_env, blacs_env, k, k)
388 : CALL cp_fm_create(Atilde, &
389 : tilde_fm_struct, &
390 248 : routineP//"Atilde")
391 : CALL cp_fm_create(Us, &
392 : tilde_fm_struct, &
393 248 : routineP//"Us")
394 :
395 : !---------------------------------------!
396 : ! calc the matrix Atilde = transp(b)*Ab !
397 : !---------------------------------------!
398 5088 : DO i = 1, k
399 171024 : DO j = 1, k
400 165936 : Atilde_ij = 0.0_dp
401 333060 : DO spin = 1, nspins
402 167124 : CALL cp_fm_trace(b(i, spin), Ab(j, spin), tmp)
403 333060 : Atilde_ij = Atilde_ij + tmp
404 : END DO
405 170776 : CALL cp_fm_set_element(Atilde, i, j, Atilde_ij)
406 : END DO
407 : END DO
408 :
409 : !--------------------!
410 : ! diagonalize Atilde !
411 : !--------------------!
412 9984 : evals_tmp(:) = evals(:, iter)
413 248 : CALL cp_fm_syevd(Atilde, Us, evals_tmp(:))
414 9984 : evals(:, iter) = evals_tmp(:)
415 :
416 : !-------------------!
417 : ! check convergence !
418 : !-------------------!
419 808 : evals_difference = 1.0_dp
420 248 : IF (iter /= 1) THEN
421 :
422 676 : evals_difference(:) = ABS((evals(1:n_ev, iter - 1) - evals(1:n_ev, iter)))
423 : ! For debugging
424 212 : IF (output_unit > 0) THEN
425 106 : WRITE (output_unit, *)
426 338 : DO i = 1, n_ev
427 338 : WRITE (output_unit, '(2X,F10.7,T69,ES11.4)') evals(i, iter)*evolt, evals_difference(i)
428 : END DO
429 106 : WRITE (output_unit, *)
430 : END IF
431 :
432 888 : convergence = MAXVAL(evals_difference)
433 212 : IF (output_unit > 0) WRITE (output_unit, '(2X,I4,T69,ES11.4)') k, convergence
434 :
435 212 : IF (convergence < tddfpt_control%tolerance) THEN
436 : res = .TRUE.
437 : EXIT iteration
438 : END IF
439 : END IF
440 :
441 236 : IF (mode == 'L') THEN
442 0 : n_kv = 1
443 : ELSE
444 764 : must_improve = 0
445 764 : DO i = 1, n_ev
446 764 : IF (evals_difference(i) > tddfpt_control%tolerance) must_improve(i) = 1
447 : END DO
448 : !! Set must_improve to 1 if all the vectors should
449 : !! be updated in one iteration.
450 : !! must_improve = 1
451 764 : n_kv = SUM(must_improve)
452 236 : j = 1
453 764 : DO i = 1, n_ev
454 764 : IF (must_improve(i) == 1) THEN
455 514 : must_improve(j) = i
456 514 : j = j + 1
457 : END IF
458 : END DO
459 : END IF
460 :
461 236 : IF (k + n_kv > max_kv) EXIT iteration
462 :
463 212 : iter = iter + 1
464 :
465 : END DO iteration
466 :
467 36 : IF (PRESENT(out_evals)) THEN
468 132 : out_evals(1:n_ev) = evals(1:n_ev, iter)
469 : END IF
470 :
471 84 : DO spin = 1, nspins
472 216 : DO j = 1, n_ev
473 132 : CALL cp_fm_set_all(t_env%evecs(j, spin), 0.0_dp)
474 1752 : DO i = 1, k
475 1572 : CALL cp_fm_get_element(Us, i, j, tmp)
476 : CALL cp_fm_scale_and_add(1.0_dp, t_env%evecs(j, spin), &
477 1704 : tmp, b(i, spin))
478 : END DO
479 : END DO
480 : END DO
481 :
482 : !----------!
483 : ! clean up !
484 : !----------!
485 36 : CALL cp_fm_release(Atilde)
486 36 : CALL cp_fm_release(Us)
487 36 : IF (ASSOCIATED(tilde_fm_struct)) CALL cp_fm_struct_release(tilde_fm_struct)
488 36 : CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools, X)
489 36 : CALL fm_pools_give_back_fm_vect(ao_mo_fm_pools, R)
490 84 : DO spin = 1, nspins
491 84 : CALL cp_fm_struct_release(kv_fm_struct(spin)%struct)
492 : END DO
493 36 : CALL cp_fm_release(b)
494 36 : CALL cp_fm_release(Ab)
495 36 : CALL cp_fm_release(Sb)
496 36 : DEALLOCATE (evals, evals_tmp, evals_difference, must_improve, kv_fm_struct)
497 :
498 36 : CALL timestop(handle)
499 :
500 72 : END FUNCTION iterative_solver
501 :
502 : ! X : the vector on which to apply the op
503 : ! R : the result
504 : ! t_env : td-dft environment (mainly control information)
505 : ! p_env : perturbation environment (variables)
506 : ! both of these carry info for the tddfpt calculation
507 : ! qs_env : info about a quickstep ground state calculation
508 :
509 : ! **************************************************************************************************
510 : !> \brief ...
511 : !> \param X ...
512 : !> \param R ...
513 : !> \param p_env ...
514 : !> \param qs_env ...
515 : !> \param do_kernel ...
516 : ! **************************************************************************************************
517 546 : SUBROUTINE apply_op(X, R, p_env, qs_env, do_kernel)
518 :
519 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: X, R
520 : TYPE(qs_p_env_type) :: p_env
521 : TYPE(qs_environment_type), POINTER :: qs_env
522 : LOGICAL, INTENT(IN) :: do_kernel
523 :
524 : CHARACTER(LEN=*), PARAMETER :: routineN = 'apply_op'
525 :
526 : INTEGER :: handle, nspins, spin
527 : INTEGER, SAVE :: counter = 0
528 : TYPE(dft_control_type), POINTER :: dft_control
529 :
530 546 : NULLIFY (dft_control)
531 :
532 546 : CALL timeset(routineN, handle)
533 :
534 546 : counter = counter + 1
535 546 : CALL get_qs_env(qs_env, dft_control=dft_control)
536 546 : nspins = dft_control%nspins
537 :
538 : !------------!
539 : ! R = HX-SXL !
540 : !------------!
541 546 : CALL p_op_l1(p_env, qs_env, X, R) ! acts on both spins, result in R
542 :
543 : !-----------------!
544 : ! calc P1 and !
545 : ! R = R + K(P1)*C !
546 : !-----------------!
547 546 : IF (do_kernel) THEN
548 1032 : DO spin = 1, nspins
549 552 : CALL dbcsr_set(p_env%p1(spin)%matrix, 0.0_dp) ! optimize?
550 : CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(spin)%matrix, &
551 : matrix_v=p_env%psi0d(spin), &
552 : matrix_g=X(spin), &
553 : ncol=p_env%n_mo(spin), &
554 1032 : symmetry_mode=1)
555 : END DO
556 1032 : DO spin = 1, nspins
557 1032 : CALL cp_fm_set_all(X(spin), 0.0_dp)
558 : END DO
559 : CALL p_op_l2(p_env, qs_env, p_env%p1, X, &
560 480 : alpha=1.0_dp, beta=0.0_dp) ! X = beta*X + alpha*K(P1)*C
561 1032 : DO spin = 1, nspins
562 : CALL cp_fm_scale_and_add(1.0_dp, R(spin), &
563 1032 : 1.0_dp, X(spin)) ! add X to R
564 : END DO
565 : END IF
566 :
567 546 : CALL timestop(handle)
568 :
569 546 : END SUBROUTINE apply_op
570 :
571 : ! **************************************************************************************************
572 : !> \brief ...
573 : !> \param vectors ...
574 : !> \param vectors_name ...
575 : !> \param startv ...
576 : !> \param n_v ...
577 : !> \param nspins ...
578 : !> \param fm_struct ...
579 : ! **************************************************************************************************
580 744 : SUBROUTINE allocate_krylov_vectors(vectors, vectors_name, &
581 744 : startv, n_v, nspins, fm_struct)
582 :
583 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: vectors
584 : CHARACTER(LEN=*), INTENT(IN) :: vectors_name
585 : INTEGER, INTENT(IN) :: startv, n_v, nspins
586 : TYPE(cp_fm_struct_p_type), DIMENSION(:), &
587 : INTENT(IN) :: fm_struct
588 :
589 : CHARACTER(LEN=*), PARAMETER :: routineN = 'allocate_krylov_vectors', &
590 : routineP = moduleN//':'//routineN
591 :
592 : CHARACTER(LEN=default_string_length) :: mat_name
593 : INTEGER :: index, spin
594 :
595 1584 : DO spin = 1, nspins
596 3510 : DO index = startv, startv + n_v - 1
597 : mat_name = routineP//vectors_name//TRIM(cp_to_string(index)) &
598 1926 : //","//TRIM(cp_to_string(spin))
599 : CALL cp_fm_create(vectors(index, spin), &
600 1926 : fm_struct(spin)%struct, mat_name)
601 1926 : IF (.NOT. ASSOCIATED(vectors(index, spin)%matrix_struct)) &
602 840 : CPABORT("Could not allocate "//TRIM(mat_name)//".")
603 : END DO
604 : END DO
605 :
606 744 : END SUBROUTINE allocate_krylov_vectors
607 :
608 : END MODULE qs_tddfpt_eigensolver
|
