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 : !> \par History
10 : !> 09-JUL-2002, TCH, development started
11 : ! **************************************************************************************************
12 : MODULE qs_tddfpt_utils
13 : USE cp_control_types, ONLY: dft_control_type,&
14 : tddfpt_control_type
15 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
16 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
17 : cp_dbcsr_sm_fm_multiply
18 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale,&
19 : cp_fm_scale_and_add,&
20 : cp_fm_trace
21 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
22 : cp_fm_cholesky_invert
23 : USE cp_fm_types, ONLY: cp_fm_create,&
24 : cp_fm_get_submatrix,&
25 : cp_fm_init_random,&
26 : cp_fm_release,&
27 : cp_fm_set_all,&
28 : cp_fm_set_submatrix,&
29 : cp_fm_to_fm,&
30 : cp_fm_type
31 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
32 : USE kinds, ONLY: dp
33 : USE physcon, ONLY: evolt
34 : USE qs_environment_types, ONLY: get_qs_env,&
35 : qs_environment_type
36 : USE qs_mo_types, ONLY: get_mo_set
37 : USE qs_p_env_methods, ONLY: p_env_create,&
38 : p_env_psi0_changed
39 : USE qs_p_env_types, ONLY: p_env_release,&
40 : qs_p_env_type
41 : USE qs_tddfpt_types, ONLY: tddfpt_env_allocate,&
42 : tddfpt_env_deallocate,&
43 : tddfpt_env_type
44 : #include "./base/base_uses.f90"
45 :
46 : IMPLICIT NONE
47 :
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt_utils'
49 : LOGICAL, PARAMETER :: DEBUG_THIS_MODULE = .TRUE.
50 :
51 : ! **************************************************************************************************
52 : TYPE simple_solution_sorter
53 : INTEGER :: orbit = -1
54 : INTEGER :: lumo = -1
55 : REAL(KIND=DP) :: value = -1.0_dp
56 : TYPE(simple_solution_sorter), POINTER :: next => NULL()
57 : END TYPE simple_solution_sorter
58 :
59 : PRIVATE
60 :
61 : ! METHODS
62 : PUBLIC :: tddfpt_cleanup, &
63 : tddfpt_init, &
64 : co_initial_guess, &
65 : find_contributions, &
66 : normalize, &
67 : reorthogonalize
68 :
69 : CONTAINS
70 :
71 : ! **************************************************************************************************
72 : !> \brief Initialize some necessary structures for a tddfpt calculation.
73 : !> \param p_env perturbation environment to be initialized
74 : !> \param t_env tddfpt environment to be initialized
75 : !> \param qs_env Quickstep environment with the results of a
76 : !> ground state calcualtion
77 : ! **************************************************************************************************
78 12 : SUBROUTINE tddfpt_init(p_env, t_env, qs_env)
79 :
80 : TYPE(qs_p_env_type), INTENT(INOUT) :: p_env
81 : TYPE(tddfpt_env_type), INTENT(out) :: t_env
82 : TYPE(qs_environment_type), POINTER :: qs_env
83 :
84 : !------------------!
85 : ! create the p_env !
86 : !------------------!
87 :
88 12 : CALL p_env_create(p_env, qs_env, orthogonal_orbitals=.TRUE.)
89 12 : CALL p_env_psi0_changed(p_env, qs_env) ! update the m_epsilon matrix
90 :
91 : !------------------!
92 : ! create the t_env !
93 : !------------------!
94 12 : CALL tddfpt_env_allocate(t_env, p_env, qs_env)
95 12 : CALL tddfpt_env_init(t_env, qs_env)
96 :
97 12 : END SUBROUTINE tddfpt_init
98 :
99 : ! **************************************************************************************************
100 : !> \brief Initialize t_env with meaningfull values.
101 : !> \param t_env ...
102 : !> \param qs_env ...
103 : ! **************************************************************************************************
104 12 : SUBROUTINE tddfpt_env_init(t_env, qs_env)
105 :
106 : TYPE(tddfpt_env_type), INTENT(inout) :: t_env
107 : TYPE(qs_environment_type), POINTER :: qs_env
108 :
109 : INTEGER :: n_spins, spin
110 12 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
111 : TYPE(dft_control_type), POINTER :: dft_control
112 :
113 12 : NULLIFY (matrix_s, dft_control)
114 :
115 12 : CALL get_qs_env(qs_env, matrix_s=matrix_s, dft_control=dft_control)
116 12 : n_spins = dft_control%nspins
117 12 : IF (dft_control%tddfpt_control%invert_S) THEN
118 28 : DO spin = 1, n_spins
119 16 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, t_env%invS(spin))
120 16 : CALL cp_fm_cholesky_decompose(t_env%invS(spin))
121 28 : CALL cp_fm_cholesky_invert(t_env%invS(spin))
122 : END DO
123 : END IF
124 :
125 12 : END SUBROUTINE tddfpt_env_init
126 :
127 : ! **************************************************************************************************
128 : !> \brief ...
129 : !> \param t_env ...
130 : !> \param p_env ...
131 : ! **************************************************************************************************
132 12 : SUBROUTINE tddfpt_cleanup(t_env, p_env)
133 :
134 : TYPE(tddfpt_env_type), INTENT(inout) :: t_env
135 : TYPE(qs_p_env_type), INTENT(INOUT) :: p_env
136 :
137 12 : CALL tddfpt_env_deallocate(t_env)
138 12 : CALL p_env_release(p_env)
139 :
140 12 : END SUBROUTINE tddfpt_cleanup
141 :
142 : ! **************************************************************************************************
143 : !> \brief ...
144 : !> \param matrices ...
145 : !> \param energies ...
146 : !> \param n_v ...
147 : !> \param qs_env ...
148 : ! **************************************************************************************************
149 12 : SUBROUTINE co_initial_guess(matrices, energies, n_v, qs_env)
150 :
151 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: matrices
152 : REAL(kind=DP), DIMENSION(:), INTENT(OUT) :: energies
153 : INTEGER, INTENT(IN) :: n_v
154 : TYPE(qs_environment_type), POINTER :: qs_env
155 :
156 : INTEGER :: i, n_cols, n_lumos, n_orbits, n_rows, &
157 : n_spins, oo, spin, vo
158 : REAL(KIND=DP) :: evd
159 12 : REAL(KIND=DP), ALLOCATABLE, DIMENSION(:, :) :: guess, lumos
160 12 : REAL(KIND=DP), DIMENSION(:), POINTER :: orbital_eigenvalues
161 : TYPE(dft_control_type), POINTER :: dft_control
162 : TYPE(simple_solution_sorter), POINTER :: sorter_iterator, sorter_pointer, &
163 : sorter_start
164 : TYPE(tddfpt_control_type), POINTER :: tddfpt_control
165 :
166 : ! number of vectors to initialize
167 :
168 12 : NULLIFY (dft_control)
169 :
170 12 : CALL get_qs_env(qs_env, dft_control=dft_control)
171 12 : tddfpt_control => dft_control%tddfpt_control
172 12 : n_spins = dft_control%nspins
173 44 : energies = 0.0_dp
174 :
175 12 : IF (.NOT. ASSOCIATED(tddfpt_control%lumos)) THEN
176 0 : CPABORT("LUMOS missing")
177 : END IF
178 :
179 28 : DO spin = 1, n_spins
180 :
181 16 : n_cols = matrices(1, spin)%matrix_struct%ncol_global
182 16 : n_rows = matrices(1, spin)%matrix_struct%nrow_global
183 :
184 60 : DO i = 1, n_v
185 60 : CALL cp_fm_set_all(matrices(i, spin), 0.0_dp)
186 : END DO
187 :
188 16 : CALL get_mo_set(qs_env%mos(spin), eigenvalues=orbital_eigenvalues)
189 :
190 16 : n_lumos = tddfpt_control%lumos(spin)%matrix_struct%ncol_global
191 :
192 16 : n_orbits = SIZE(orbital_eigenvalues)
193 :
194 : !-----------------------------------------!
195 : ! create a SORTED list of initial guesses !
196 : !-----------------------------------------!
197 : ! first element
198 16 : evd = tddfpt_control%lumos_eigenvalues(1, spin) - orbital_eigenvalues(n_orbits)
199 16 : ALLOCATE (sorter_start)
200 16 : sorter_start%orbit = n_orbits
201 16 : sorter_start%lumo = 1
202 16 : sorter_start%value = evd
203 : NULLIFY (sorter_start%next)
204 : ! rest of the elements
205 80 : DO oo = n_orbits, 1, -1
206 376 : DO vo = 1, n_lumos
207 :
208 296 : IF (oo == n_orbits .AND. vo == 1) CYCLE ! already in list
209 :
210 280 : evd = tddfpt_control%lumos_eigenvalues(vo, spin) - orbital_eigenvalues(oo)
211 :
212 280 : sorter_iterator => sorter_start
213 280 : NULLIFY (sorter_pointer)
214 2314 : DO WHILE (ASSOCIATED(sorter_iterator%next))
215 2134 : IF (sorter_iterator%next%value > evd) THEN
216 : sorter_pointer => sorter_iterator%next
217 : EXIT
218 : END IF
219 180 : sorter_iterator => sorter_iterator%next
220 : END DO
221 :
222 280 : ALLOCATE (sorter_iterator%next)
223 280 : sorter_iterator%next%orbit = oo
224 280 : sorter_iterator%next%lumo = vo
225 280 : sorter_iterator%next%value = evd
226 360 : sorter_iterator%next%next => sorter_pointer
227 :
228 : END DO
229 : END DO
230 :
231 112 : ALLOCATE (lumos(n_rows, n_lumos), guess(n_rows, n_orbits))
232 : CALL cp_fm_get_submatrix(tddfpt_control%lumos(spin), lumos, &
233 16 : start_col=1, n_cols=n_lumos)
234 :
235 : !-------------------!
236 : ! fill the matrices !
237 : !-------------------!
238 16 : sorter_iterator => sorter_start
239 60 : DO i = 1, MIN(n_v, n_orbits*n_lumos)
240 3996 : guess(:, :) = 0.0_dp
241 : CALL dcopy(n_rows, lumos(:, sorter_iterator%lumo), 1, &
242 44 : guess(:, sorter_iterator%orbit), 1)
243 44 : CALL cp_fm_set_submatrix(matrices(i, spin), guess)
244 44 : energies(i) = energies(i) + sorter_iterator%value/REAL(n_spins, dp)
245 60 : sorter_iterator => sorter_iterator%next
246 : END DO
247 16 : IF (n_v > n_orbits*n_lumos) THEN
248 0 : DO i = n_orbits*n_lumos + 1, n_v
249 0 : CALL cp_fm_init_random(matrices(i, spin), n_orbits)
250 0 : energies(i) = 1.0E38_dp
251 : END DO
252 : END IF
253 :
254 : !--------------!
255 : ! some cleanup !
256 : !--------------!
257 16 : DEALLOCATE (lumos, guess)
258 16 : sorter_iterator => sorter_start
259 324 : DO WHILE (ASSOCIATED(sorter_iterator))
260 296 : sorter_pointer => sorter_iterator
261 296 : sorter_iterator => sorter_iterator%next
262 312 : DEALLOCATE (sorter_pointer)
263 : END DO
264 :
265 : END DO
266 :
267 24 : END SUBROUTINE co_initial_guess
268 :
269 : ! **************************************************************************************************
270 : !> \brief ...
271 : !> \param qs_env ...
272 : !> \param t_env ...
273 : ! **************************************************************************************************
274 12 : SUBROUTINE find_contributions(qs_env, t_env)
275 :
276 : TYPE(qs_environment_type), POINTER :: qs_env
277 : TYPE(tddfpt_env_type), INTENT(IN) :: t_env
278 :
279 : INTEGER :: i, j, n_ev, n_spins, occ, output_unit, &
280 : spin, virt
281 : INTEGER, DIMENSION(2) :: nhomos, nlumos, nrows
282 : REAL(KIND=dp) :: contribution, summed_contributions
283 12 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: homo_coeff_col, lumo_coeff_col
284 12 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: S_lumos
285 12 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
286 : TYPE(dft_control_type), POINTER :: dft_control
287 : TYPE(tddfpt_control_type) :: t_control
288 :
289 12 : NULLIFY (matrix_s, dft_control)
290 24 : output_unit = cp_logger_get_default_io_unit()
291 12 : CALL get_qs_env(qs_env, matrix_s=matrix_s, dft_control=dft_control)
292 :
293 12 : IF (output_unit > 0) WRITE (output_unit, *)
294 6 : IF (output_unit > 0) WRITE (output_unit, *)
295 :
296 12 : t_control = dft_control%tddfpt_control
297 12 : n_ev = t_control%n_ev
298 12 : n_spins = dft_control%nspins
299 :
300 52 : ALLOCATE (S_lumos(n_spins))
301 :
302 28 : DO spin = 1, n_spins
303 16 : nrows(spin) = t_control%lumos(spin)%matrix_struct%nrow_global
304 16 : nhomos(spin) = t_env%evecs(1, spin)%matrix_struct%ncol_global
305 16 : nlumos(spin) = t_control%lumos(spin)%matrix_struct%ncol_global
306 : CALL cp_fm_create(S_lumos(spin), t_control%lumos(spin)%matrix_struct, &
307 16 : "S times lumos")
308 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, t_control%lumos(spin), &
309 28 : S_lumos(spin), nlumos(spin), 1.0_dp, 0.0_dp)
310 : END DO
311 :
312 : ALLOCATE (homo_coeff_col(MAXVAL(nrows(1:n_spins)), 1), &
313 92 : lumo_coeff_col(MAXVAL(nrows(1:n_spins)), 1))
314 44 : DO i = 1, n_ev
315 32 : IF (output_unit > 0) THEN
316 16 : WRITE (output_unit, '(A,I3,5X,F15.6)') " excited state : ", i, t_env%evals(i)*evolt
317 16 : WRITE (output_unit, *)
318 : END IF
319 32 : summed_contributions = 0.0_dp
320 76 : DO spin = 1, n_spins
321 44 : IF (n_spins == 2) THEN
322 24 : IF (spin == 1) THEN
323 12 : IF (output_unit > 0) WRITE (output_unit, *) 'alpha:'
324 : ELSE
325 12 : IF (output_unit > 0) WRITE (output_unit, *) 'beta:'
326 : END IF
327 : END IF
328 252 : searchloop: DO occ = nhomos(spin), 1, -1
329 : CALL cp_fm_get_submatrix(t_env%evecs(i, spin), homo_coeff_col, &
330 176 : 1, occ, nrows(spin), 1)
331 1052 : DO virt = 1, nlumos(spin)
332 : CALL cp_fm_get_submatrix(S_lumos(spin), lumo_coeff_col, &
333 832 : 1, virt, nrows(spin), 1)
334 832 : contribution = 0.0_dp
335 19424 : DO j = 1, nrows(spin)
336 19424 : contribution = contribution + homo_coeff_col(j, 1)*lumo_coeff_col(j, 1)
337 : END DO
338 832 : summed_contributions = summed_contributions + (contribution)**2
339 832 : IF (ABS(contribution) > 5.0e-2_dp) THEN
340 244 : IF (output_unit > 0) WRITE (output_unit, '(14X,I5,A,I5,10X,F8.3,5X,F8.3)') &
341 122 : occ, " ->", nhomos(spin) + virt, ABS(contribution), summed_contributions
342 : END IF
343 1008 : IF (ABS(summed_contributions - 1.0_dp) < 1.0e-3_dp) CYCLE searchloop
344 : END DO
345 : END DO searchloop
346 : END DO
347 44 : IF (output_unit > 0) WRITE (output_unit, *)
348 : END DO
349 :
350 : !
351 : ! punch a checksum for the regs
352 12 : IF (output_unit > 0) THEN
353 22 : WRITE (output_unit, '(T2,A,E14.6)') ' TDDFPT : CheckSum =', SQRT(SUM(t_env%evals**2))
354 : END IF
355 :
356 12 : CALL cp_fm_release(S_lumos)
357 :
358 12 : DEALLOCATE (homo_coeff_col, lumo_coeff_col)
359 :
360 24 : END SUBROUTINE find_contributions
361 :
362 : ! **************************************************************************************************
363 : !> \brief ...
364 : !> \param X ...
365 : !> \param tmp_vec ...
366 : !> \param metric ...
367 : ! **************************************************************************************************
368 546 : SUBROUTINE normalize(X, tmp_vec, metric)
369 :
370 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: x, tmp_vec
371 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: metric
372 :
373 : INTEGER :: n_spins, spin
374 : REAL(KIND=dp) :: norm, tmp
375 :
376 546 : n_spins = SIZE(x)
377 546 : norm = 0.0_dp
378 :
379 1188 : DO spin = 1, n_spins
380 : tmp = 0.0_dp
381 : CALL cp_dbcsr_sm_fm_multiply(metric(1)%matrix, X(spin), &
382 : tmp_vec(spin), &
383 : X(spin)%matrix_struct%ncol_global, &
384 642 : 1.0_dp, 0.0_dp)
385 642 : CALL cp_fm_trace(X(spin), tmp_vec(spin), tmp)
386 1188 : norm = norm + tmp
387 : END DO
388 :
389 546 : norm = SQRT(norm)
390 1188 : DO spin = 1, n_spins
391 1188 : CALL cp_fm_scale((1.0_dp/norm), X(spin))
392 : END DO
393 :
394 546 : END SUBROUTINE normalize
395 :
396 : !---------------------------------------!
397 : ! x must not be changed in this routine !
398 : ! tmp_vec may be changed !
399 : !---------------------------------------!
400 : ! **************************************************************************************************
401 : !> \brief ...
402 : !> \param X ...
403 : !> \param V_set ...
404 : !> \param SV_set ...
405 : !> \param work ...
406 : !> \param n ...
407 : ! **************************************************************************************************
408 1092 : SUBROUTINE reorthogonalize(X, V_set, SV_set, work, n)
409 :
410 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: X
411 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: V_set, SV_set
412 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: work
413 : INTEGER, INTENT(IN) :: n
414 :
415 : CHARACTER(LEN=*), PARAMETER :: routineN = 'reorthogonalize'
416 :
417 : INTEGER :: handle, i, n_spins, spin
418 : REAL(DP) :: dot_product, tmp
419 :
420 1092 : CALL timeset(routineN, handle)
421 :
422 1092 : IF (n > 0) THEN
423 :
424 1020 : n_spins = SIZE(X)
425 2208 : DO spin = 1, n_spins
426 2208 : CALL cp_fm_to_fm(X(spin), work(spin))
427 : END DO
428 :
429 18844 : DO i = 1, n
430 : dot_product = 0.0_dp
431 36344 : DO spin = 1, n_spins
432 18520 : CALL cp_fm_trace(SV_set(i, spin), work(spin), tmp)
433 36344 : dot_product = dot_product + tmp
434 : END DO
435 37364 : DO spin = 1, n_spins
436 : CALL cp_fm_scale_and_add(1.0_dp, X(spin), &
437 36344 : -1.0_dp*dot_product, V_set(i, spin))
438 : END DO
439 : END DO
440 :
441 : END IF
442 :
443 1092 : CALL timestop(handle)
444 :
445 1092 : END SUBROUTINE reorthogonalize
446 :
447 : ! **************************************************************************************************
448 :
449 0 : END MODULE qs_tddfpt_utils
|
