(git:339af6c)
Loading...
Searching...
No Matches
qs_density_matrices.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief collects routines that calculate density matrices
10!> \note
11!> first version : most routines imported
12!> \author JGH (2020-01)
13! **************************************************************************************************
14MODULE qs_density_matrices
15 USE cp_blacs_env, ONLY: cp_blacs_env_type
16 USE cp_dbcsr_api, ONLY: dbcsr_copy,&
17 dbcsr_multiply,&
18 dbcsr_release,&
19 dbcsr_set,&
20 dbcsr_type
21 USE cp_dbcsr_contrib, ONLY: dbcsr_scale_by_vector
22 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
23 copy_fm_to_dbcsr,&
24 cp_dbcsr_plus_fm_fm_t,&
25 cp_dbcsr_sm_fm_multiply
26 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
27 cp_fm_scale_and_add,&
28 cp_fm_symm,&
29 cp_fm_transpose,&
30 cp_fm_uplo_to_full
31 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
32 cp_fm_struct_release,&
33 cp_fm_struct_type
34 USE cp_fm_types, ONLY: cp_fm_create,&
35 cp_fm_get_info,&
36 cp_fm_release,&
37 cp_fm_to_fm,&
38 cp_fm_type
39 USE cp_log_handling, ONLY: cp_get_default_logger,&
40 cp_logger_get_default_unit_nr,&
41 cp_logger_type
42 USE kinds, ONLY: dp
43 USE message_passing, ONLY: mp_para_env_type
44 USE parallel_gemm_api, ONLY: parallel_gemm
45 USE qs_mo_types, ONLY: get_mo_set,&
46 mo_set_type
47#include "./base/base_uses.f90"
48
49 IMPLICIT NONE
50 PRIVATE
51 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_density_matrices'
52
53 PUBLIC :: calculate_density_matrix
54 PUBLIC :: calculate_w_matrix, calculate_w_matrix_ot
55 PUBLIC :: calculate_wz_matrix, calculate_whz_matrix
56 PUBLIC :: calculate_wx_matrix, calculate_xwx_matrix
57
58 INTERFACE calculate_density_matrix
59 MODULE PROCEDURE calculate_dm_sparse
60 END INTERFACE
61
62 INTERFACE calculate_w_matrix
63 MODULE PROCEDURE calculate_w_matrix_1, calculate_w_matrix_roks
64 END INTERFACE
65
66CONTAINS
67
68! **************************************************************************************************
69!> \brief Calculate the density matrix
70!> \param mo_set ...
71!> \param density_matrix ...
72!> \param use_dbcsr ...
73!> \param retain_sparsity ...
74!> \date 06.2002
75!> \par History
76!> - Fractional occupied orbitals (MK)
77!> \author Joost VandeVondele
78!> \version 1.0
79! **************************************************************************************************
80 SUBROUTINE calculate_dm_sparse(mo_set, density_matrix, use_dbcsr, retain_sparsity)
81
82 TYPE(mo_set_type), INTENT(IN) :: mo_set
83 TYPE(dbcsr_type), POINTER :: density_matrix
84 LOGICAL, INTENT(IN), OPTIONAL :: use_dbcsr, retain_sparsity
85
86 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_dm_sparse'
87
88 INTEGER :: handle
89 LOGICAL :: my_retain_sparsity, my_use_dbcsr
90 REAL(KIND=dp) :: alpha
91 TYPE(cp_fm_type) :: fm_tmp
92 TYPE(dbcsr_type) :: dbcsr_tmp
93
94 CALL timeset(routinen, handle)
95
96 my_use_dbcsr = .false.
97 IF (PRESENT(use_dbcsr)) my_use_dbcsr = use_dbcsr
98 my_retain_sparsity = .true.
99 IF (PRESENT(retain_sparsity)) my_retain_sparsity = retain_sparsity
100 IF (my_use_dbcsr) THEN
101 IF (.NOT. ASSOCIATED(mo_set%mo_coeff_b)) THEN
102 cpabort("mo_coeff_b NOT ASSOCIATED")
103 END IF
104 END IF
105
106 CALL dbcsr_set(density_matrix, 0.0_dp)
107
108 IF (.NOT. mo_set%uniform_occupation) THEN ! not all orbitals 1..homo are equally occupied
109 IF (my_use_dbcsr) THEN
110 CALL dbcsr_copy(dbcsr_tmp, mo_set%mo_coeff_b)
111 CALL dbcsr_scale_by_vector(dbcsr_tmp, mo_set%occupation_numbers(1:mo_set%homo), &
112 side='right')
113 CALL dbcsr_multiply("N", "T", 1.0_dp, mo_set%mo_coeff_b, dbcsr_tmp, &
114 1.0_dp, density_matrix, retain_sparsity=my_retain_sparsity, &
115 last_k=mo_set%homo)
116 CALL dbcsr_release(dbcsr_tmp)
117 ELSE
118 CALL cp_fm_create(fm_tmp, mo_set%mo_coeff%matrix_struct)
119 CALL cp_fm_to_fm(mo_set%mo_coeff, fm_tmp)
120 CALL cp_fm_column_scale(fm_tmp, mo_set%occupation_numbers(1:mo_set%homo))
121 alpha = 1.0_dp
122 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix, &
123 matrix_v=mo_set%mo_coeff, &
124 matrix_g=fm_tmp, &
125 ncol=mo_set%homo, &
126 alpha=alpha)
127 CALL cp_fm_release(fm_tmp)
128 END IF
129 ELSE
130 IF (my_use_dbcsr) THEN
131 CALL dbcsr_multiply("N", "T", mo_set%maxocc, mo_set%mo_coeff_b, mo_set%mo_coeff_b, &
132 1.0_dp, density_matrix, retain_sparsity=my_retain_sparsity, &
133 last_k=mo_set%homo)
134 ELSE
135 alpha = mo_set%maxocc
136 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix, &
137 matrix_v=mo_set%mo_coeff, &
138 ncol=mo_set%homo, &
139 alpha=alpha)
140 END IF
141 END IF
142
143 CALL timestop(handle)
144
145 END SUBROUTINE calculate_dm_sparse
146
147! **************************************************************************************************
148!> \brief Calculate the W matrix from the MO eigenvectors, MO eigenvalues,
149!> and the MO occupation numbers. Only works if they are eigenstates
150!> \param mo_set type containing the full matrix of the MO and the eigenvalues
151!> \param w_matrix sparse matrix
152!> error
153!> \par History
154!> Creation (03.03.03,MK)
155!> Modification that computes it as a full block, several times (e.g. 20)
156!> faster at the cost of some additional memory
157!> \author MK
158! **************************************************************************************************
159 SUBROUTINE calculate_w_matrix_1(mo_set, w_matrix)
160
161 TYPE(mo_set_type), INTENT(IN) :: mo_set
162 TYPE(dbcsr_type), POINTER :: w_matrix
163
164 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_1'
165
166 INTEGER :: handle, imo
167 REAL(KIND=dp), DIMENSION(:), POINTER :: eigocc
168 TYPE(cp_fm_type) :: weighted_vectors
169
170 CALL timeset(routinen, handle)
171
172 CALL dbcsr_set(w_matrix, 0.0_dp)
173 CALL cp_fm_create(weighted_vectors, mo_set%mo_coeff%matrix_struct, "weighted_vectors")
174 CALL cp_fm_to_fm(mo_set%mo_coeff, weighted_vectors)
175
176 ! scale every column with the occupation
177 ALLOCATE (eigocc(mo_set%homo))
178
179 DO imo = 1, mo_set%homo
180 eigocc(imo) = mo_set%eigenvalues(imo)*mo_set%occupation_numbers(imo)
181 END DO
182 CALL cp_fm_column_scale(weighted_vectors, eigocc)
183 DEALLOCATE (eigocc)
184
185 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=w_matrix, &
186 matrix_v=mo_set%mo_coeff, &
187 matrix_g=weighted_vectors, &
188 ncol=mo_set%homo)
189
190 CALL cp_fm_release(weighted_vectors)
191
192 CALL timestop(handle)
193
194 END SUBROUTINE calculate_w_matrix_1
195
196! **************************************************************************************************
197!> \brief Calculate the W matrix from the MO coefs, MO derivs
198!> could overwrite the mo_derivs for increased memory efficiency
199!> \param mo_set type containing the full matrix of the MO coefs
200!> mo_deriv:
201!> \param mo_deriv ...
202!> \param w_matrix sparse matrix
203!> \param s_matrix sparse matrix for the overlap
204!> error
205!> \par History
206!> Creation (JV)
207!> \author MK
208! **************************************************************************************************
209 SUBROUTINE calculate_w_matrix_ot(mo_set, mo_deriv, w_matrix, s_matrix)
210
211 TYPE(mo_set_type), INTENT(IN) :: mo_set
212 TYPE(dbcsr_type), POINTER :: mo_deriv, w_matrix, s_matrix
213
214 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_w_matrix_ot'
215 LOGICAL, PARAMETER :: check_gradient = .false., &
216 do_symm = .false.
217
218 INTEGER :: handle, iounit, ncol_global, nrow_global
219 REAL(kind=dp), DIMENSION(:), POINTER :: occupation_numbers, scaling_factor
220 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
221 TYPE(cp_fm_type) :: gradient, h_block, h_block_t, &
222 weighted_vectors
223 TYPE(cp_logger_type), POINTER :: logger
224
225 CALL timeset(routinen, handle)
226 NULLIFY (fm_struct_tmp)
227
228 CALL cp_fm_get_info(matrix=mo_set%mo_coeff, &
229 ncol_global=ncol_global, &
230 nrow_global=nrow_global)
231
232 CALL cp_fm_create(weighted_vectors, mo_set%mo_coeff%matrix_struct, "weighted_vectors")
233 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol_global, ncol_global=ncol_global, &
234 para_env=mo_set%mo_coeff%matrix_struct%para_env, &
235 context=mo_set%mo_coeff%matrix_struct%context)
236 CALL cp_fm_create(h_block, fm_struct_tmp, name="h block", set_zero=.true.)
237 IF (do_symm) CALL cp_fm_create(h_block_t, fm_struct_tmp, name="h block t")
238 CALL cp_fm_struct_release(fm_struct_tmp)
239
240 CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers)
241 ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
242 scaling_factor = 2.0_dp*occupation_numbers
243 CALL copy_dbcsr_to_fm(mo_deriv, weighted_vectors)
244 CALL cp_fm_column_scale(weighted_vectors, scaling_factor)
245 DEALLOCATE (scaling_factor)
246
247 ! the convention seems to require the half here, the factor of two is presumably taken care of
248 ! internally in qs_core_hamiltonian
249 CALL parallel_gemm('T', 'N', ncol_global, ncol_global, nrow_global, 0.5_dp, &
250 mo_set%mo_coeff, weighted_vectors, 0.0_dp, h_block)
251
252 IF (do_symm) THEN
253 ! at the minimum things are anyway symmetric, but numerically it might not be the case
254 ! needs some investigation to find out if using this is better
255 CALL cp_fm_transpose(h_block, h_block_t)
256 CALL cp_fm_scale_and_add(0.5_dp, h_block, 0.5_dp, h_block_t)
257 END IF
258
259 ! this could overwrite the mo_derivs to save the weighted_vectors
260 CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
261 mo_set%mo_coeff, h_block, 0.0_dp, weighted_vectors)
262
263 CALL dbcsr_set(w_matrix, 0.0_dp)
264 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=w_matrix, &
265 matrix_v=mo_set%mo_coeff, &
266 matrix_g=weighted_vectors, &
267 ncol=mo_set%homo)
268
269 IF (check_gradient) THEN
270 CALL cp_fm_create(gradient, mo_set%mo_coeff%matrix_struct, "gradient")
271 CALL cp_dbcsr_sm_fm_multiply(s_matrix, weighted_vectors, &
272 gradient, ncol_global)
273
274 ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
275 scaling_factor = 2.0_dp*occupation_numbers
276 CALL copy_dbcsr_to_fm(mo_deriv, weighted_vectors)
277 CALL cp_fm_column_scale(weighted_vectors, scaling_factor)
278 DEALLOCATE (scaling_factor)
279
280 logger => cp_get_default_logger()
281 IF (logger%para_env%is_source()) THEN
282 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
283 WRITE (iounit, *) " maxabs difference ", &
284 maxval(abs(weighted_vectors%local_data - 2.0_dp*gradient%local_data))
285 END IF
286 CALL cp_fm_release(gradient)
287 END IF
288
289 IF (do_symm) CALL cp_fm_release(h_block_t)
290 CALL cp_fm_release(weighted_vectors)
291 CALL cp_fm_release(h_block)
292
293 CALL timestop(handle)
294
295 END SUBROUTINE calculate_w_matrix_ot
296
297! **************************************************************************************************
298!> \brief Calculate the energy-weighted density matrix W if ROKS is active.
299!> The W matrix is returned in matrix_w.
300!> \param mo_set ...
301!> \param matrix_ks ...
302!> \param matrix_p ...
303!> \param matrix_w ...
304!> \author 04.05.06,MK
305! **************************************************************************************************
306 SUBROUTINE calculate_w_matrix_roks(mo_set, matrix_ks, matrix_p, matrix_w)
307 TYPE(mo_set_type), INTENT(IN) :: mo_set
308 TYPE(dbcsr_type), POINTER :: matrix_ks, matrix_p, matrix_w
309
310 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_roks'
311
312 INTEGER :: handle, nao
313 TYPE(cp_blacs_env_type), POINTER :: context
314 TYPE(cp_fm_struct_type), POINTER :: fm_struct
315 TYPE(cp_fm_type) :: ks, p, work
316 TYPE(cp_fm_type), POINTER :: c
317 TYPE(mp_para_env_type), POINTER :: para_env
318
319 CALL timeset(routinen, handle)
320
321 NULLIFY (context)
322 NULLIFY (fm_struct)
323 NULLIFY (para_env)
324
325 CALL get_mo_set(mo_set=mo_set, mo_coeff=c)
326 CALL cp_fm_get_info(c, context=context, nrow_global=nao, para_env=para_env)
327 CALL cp_fm_struct_create(fm_struct, context=context, nrow_global=nao, &
328 ncol_global=nao, para_env=para_env)
329 CALL cp_fm_create(ks, fm_struct, name="Kohn-Sham matrix")
330 CALL cp_fm_create(p, fm_struct, name="Density matrix")
331 CALL cp_fm_create(work, fm_struct, name="Work matrix")
332 CALL cp_fm_struct_release(fm_struct)
333 CALL copy_dbcsr_to_fm(matrix_ks, ks)
334 CALL copy_dbcsr_to_fm(matrix_p, p)
335 CALL cp_fm_uplo_to_full(p, work)
336 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ks, p, 0.0_dp, work)
337 CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, p, work, 0.0_dp, ks)
338 CALL dbcsr_set(matrix_w, 0.0_dp)
339 CALL copy_fm_to_dbcsr(ks, matrix_w, keep_sparsity=.true.)
340
341 CALL cp_fm_release(work)
342 CALL cp_fm_release(p)
343 CALL cp_fm_release(ks)
344
345 CALL timestop(handle)
346
347 END SUBROUTINE calculate_w_matrix_roks
348
349! **************************************************************************************************
350!> \brief Calculate the response W matrix from the MO eigenvectors, MO eigenvalues,
351!> and the MO occupation numbers. Only works if they are eigenstates
352!> \param mo_set type containing the full matrix of the MO and the eigenvalues
353!> \param psi1 response orbitals
354!> \param ks_matrix Kohn-Sham sparse matrix
355!> \param w_matrix sparse matrix
356!> \par History
357!> adapted from calculate_w_matrix_1
358!> \author JGH
359! **************************************************************************************************
360 SUBROUTINE calculate_wz_matrix(mo_set, psi1, ks_matrix, w_matrix)
361
362 TYPE(mo_set_type), INTENT(IN) :: mo_set
363 TYPE(cp_fm_type), INTENT(IN) :: psi1
364 TYPE(dbcsr_type), POINTER :: ks_matrix, w_matrix
365
366 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_wz_matrix'
367
368 INTEGER :: handle, ncol, nocc, nrow
369 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
370 TYPE(cp_fm_type) :: ksmat, scrv
371
372 CALL timeset(routinen, handle)
373
374 CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow)
375 nocc = mo_set%homo
376 CALL cp_fm_create(scrv, mo_set%mo_coeff%matrix_struct, "scr vectors", set_zero=.true.)
377 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nocc, ncol_global=nocc, &
378 para_env=mo_set%mo_coeff%matrix_struct%para_env, &
379 context=mo_set%mo_coeff%matrix_struct%context)
380 CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
381 CALL cp_fm_struct_release(fm_struct_tmp)
382 CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, nocc)
383 CALL parallel_gemm("T", "N", nocc, nocc, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat)
384 CALL parallel_gemm("N", "N", nrow, nocc, nocc, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv)
385 CALL dbcsr_set(w_matrix, 0.0_dp)
386 CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, ncol=nocc, symmetry_mode=1)
387 CALL cp_fm_release(scrv)
388 CALL cp_fm_release(ksmat)
389
390 CALL timestop(handle)
391
392 END SUBROUTINE calculate_wz_matrix
393
394! **************************************************************************************************
395!> \brief Calculate the Wz matrix from the MO eigenvectors, MO eigenvalues,
396!> and the MO occupation numbers. Only works if they are eigenstates
397!> \param c0vec ...
398!> \param hzm ...
399!> \param w_matrix sparse matrix
400!> \param focc ...
401!> \param nocc ...
402!> \par History
403!> adapted from calculate_w_matrix_1
404!> \author JGH
405! **************************************************************************************************
406 SUBROUTINE calculate_whz_matrix(c0vec, hzm, w_matrix, focc, nocc)
407
408 TYPE(cp_fm_type), INTENT(IN) :: c0vec
409 TYPE(dbcsr_type), POINTER :: hzm, w_matrix
410 REAL(kind=dp), INTENT(IN) :: focc
411 INTEGER, INTENT(IN) :: nocc
412
413 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_whz_matrix'
414
415 INTEGER :: handle, nao, norb
416 REAL(kind=dp) :: falpha
417 TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_struct_mat
418 TYPE(cp_fm_type) :: chcmat, hcvec
419
420 CALL timeset(routinen, handle)
421
422 falpha = focc
423
424 CALL cp_fm_create(hcvec, c0vec%matrix_struct, name="hcvec", set_zero=.true.)
425 CALL cp_fm_get_info(hcvec, matrix_struct=fm_struct, nrow_global=nao, ncol_global=norb)
426 cpassert(nocc <= norb .AND. nocc > 0)
427 norb = nocc
428 CALL cp_fm_struct_create(fm_struct_mat, context=fm_struct%context, nrow_global=norb, &
429 ncol_global=norb, para_env=fm_struct%para_env)
430 CALL cp_fm_create(chcmat, fm_struct_mat)
431 CALL cp_fm_struct_release(fm_struct_mat)
432
433 CALL cp_dbcsr_sm_fm_multiply(hzm, c0vec, hcvec, norb)
434 CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, c0vec, hcvec, 0.0_dp, chcmat)
435 CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, c0vec, chcmat, 0.0_dp, hcvec)
436
437 CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=hcvec, matrix_g=c0vec, ncol=norb, alpha=falpha)
438
439 CALL cp_fm_release(hcvec)
440 CALL cp_fm_release(chcmat)
441
442 CALL timestop(handle)
443
444 END SUBROUTINE calculate_whz_matrix
445
446! **************************************************************************************************
447!> \brief Calculate the excited state W matrix from the MO eigenvectors, KS matrix
448!> \param mos_active ...
449!> \param xvec ...
450!> \param ks_matrix ...
451!> \param w_matrix ...
452!> \par History
453!> adapted from calculate_wz_matrix
454!> \author JGH
455! **************************************************************************************************
456 SUBROUTINE calculate_wx_matrix(mos_active, xvec, ks_matrix, w_matrix)
457
458 TYPE(cp_fm_type), INTENT(IN) :: mos_active, xvec
459 TYPE(dbcsr_type), POINTER :: ks_matrix, w_matrix
460
461 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_wx_matrix'
462
463 INTEGER :: handle, ncol, nrow
464 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
465 TYPE(cp_fm_type) :: ksmat, scrv
466
467 CALL timeset(routinen, handle)
468
469 CALL cp_fm_get_info(matrix=mos_active, ncol_global=ncol, nrow_global=nrow)
470 CALL cp_fm_create(scrv, mos_active%matrix_struct, name="scr vectors", set_zero=.true.)
471 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
472 para_env=mos_active%matrix_struct%para_env, &
473 context=mos_active%matrix_struct%context)
474 CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
475 CALL cp_fm_struct_release(fm_struct_tmp)
476 CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mos_active, scrv, ncol)
477 CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mos_active, scrv, 0.0_dp, ksmat)
478 CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, xvec, ksmat, 0.0_dp, scrv)
479 CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=xvec, ncol=ncol, symmetry_mode=1)
480 CALL cp_fm_release(scrv)
481 CALL cp_fm_release(ksmat)
482
483 CALL timestop(handle)
484
485 END SUBROUTINE calculate_wx_matrix
486
487! **************************************************************************************************
488!> \brief Calculate the excited state W matrix from the MO eigenvectors, KS matrix
489!> \param mos_active ...
490!> \param xvec ...
491!> \param s_matrix ...
492!> \param ks_matrix ...
493!> \param w_matrix ...
494!> \param eval ...
495!> \par History
496!> adapted from calculate_wz_matrix
497!> \author JGH
498! **************************************************************************************************
499 SUBROUTINE calculate_xwx_matrix(mos_active, xvec, s_matrix, ks_matrix, w_matrix, eval)
500
501 TYPE(cp_fm_type), INTENT(IN) :: mos_active, xvec
502 TYPE(dbcsr_type), POINTER :: s_matrix, ks_matrix, w_matrix
503 REAL(kind=dp), INTENT(IN) :: eval
504
505 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_xwx_matrix'
506
507 INTEGER :: handle, ncol, nrow
508 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
509 TYPE(cp_fm_type) :: scrv, xsxmat
510
511 CALL timeset(routinen, handle)
512
513 CALL cp_fm_get_info(matrix=mos_active, ncol_global=ncol, nrow_global=nrow)
514 CALL cp_fm_create(scrv, mos_active%matrix_struct, name="scr vectors", set_zero=.true.)
515 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
516 para_env=mos_active%matrix_struct%para_env, &
517 context=mos_active%matrix_struct%context)
518 CALL cp_fm_create(xsxmat, fm_struct_tmp, name="XSX")
519 CALL cp_fm_struct_release(fm_struct_tmp)
520
521 CALL cp_dbcsr_sm_fm_multiply(ks_matrix, xvec, scrv, ncol, 1.0_dp, 0.0_dp)
522 CALL cp_dbcsr_sm_fm_multiply(s_matrix, xvec, scrv, ncol, eval, -1.0_dp)
523 CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, xvec, scrv, 0.0_dp, xsxmat)
524
525 CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mos_active, xsxmat, 0.0_dp, scrv)
526 CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=mos_active, ncol=ncol, symmetry_mode=1)
527
528 CALL cp_fm_release(scrv)
529 CALL cp_fm_release(xsxmat)
530
531 CALL timestop(handle)
532
533 END SUBROUTINE calculate_xwx_matrix
534
535END MODULE qs_density_matrices
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:84
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
qs_density_matrices::calculate_density_matrix
Definition qs_density_matrices.F:58
qs_density_matrices::calculate_w_matrix
Definition qs_density_matrices.F:62
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1086
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1194
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1132
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_scale_by_vector
subroutine, public dbcsr_scale_by_vector(matrix, alpha, side)
Scales the rows/columns of given matrix.
Definition cp_dbcsr_contrib.F:486
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:552
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:214
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
Definition cp_dbcsr_operations.F:695
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:111
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1937
cp_fm_basic_linalg::cp_fm_transpose
subroutine, public cp_fm_transpose(matrix, matrixt)
transposes a matrix matrixt = matrix ^ T
Definition cp_fm_basic_linalg.F:1749
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:167
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1801
cp_fm_basic_linalg::cp_fm_symm
subroutine, public cp_fm_symm(side, uplo, m, n, alpha, matrix_a, matrix_b, beta, matrix_c)
computes matrix_c = beta * matrix_c + alpha * matrix_a * matrix_b computes matrix_c = beta * matrix_c...
Definition cp_fm_basic_linalg.F:620
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1087
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:169
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
qs_density_matrices
collects routines that calculate density matrices
Definition qs_density_matrices.F:14
qs_density_matrices::calculate_whz_matrix
subroutine, public calculate_whz_matrix(c0vec, hzm, w_matrix, focc, nocc)
Calculate the Wz matrix from the MO eigenvectors, MO eigenvalues, and the MO occupation numbers....
Definition qs_density_matrices.F:407
qs_density_matrices::calculate_xwx_matrix
subroutine, public calculate_xwx_matrix(mos_active, xvec, s_matrix, ks_matrix, w_matrix, eval)
Calculate the excited state W matrix from the MO eigenvectors, KS matrix.
Definition qs_density_matrices.F:500
qs_density_matrices::calculate_w_matrix_ot
subroutine, public calculate_w_matrix_ot(mo_set, mo_deriv, w_matrix, s_matrix)
Calculate the W matrix from the MO coefs, MO derivs could overwrite the mo_derivs for increased memor...
Definition qs_density_matrices.F:210
qs_density_matrices::calculate_wz_matrix
subroutine, public calculate_wz_matrix(mo_set, psi1, ks_matrix, w_matrix)
Calculate the response W matrix from the MO eigenvectors, MO eigenvalues, and the MO occupation numbe...
Definition qs_density_matrices.F:361
qs_density_matrices::calculate_wx_matrix
subroutine, public calculate_wx_matrix(mos_active, xvec, ks_matrix, w_matrix)
Calculate the excited state W matrix from the MO eigenvectors, KS matrix.
Definition qs_density_matrices.F:457
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:399
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:115
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:729
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47