(git:374b731)
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-2024 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_operations, ONLY: copy_dbcsr_to_fm,&
17 copy_fm_to_dbcsr,&
18 cp_dbcsr_plus_fm_fm_t,&
19 cp_dbcsr_sm_fm_multiply
20 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
21 cp_fm_scale_and_add,&
22 cp_fm_symm,&
23 cp_fm_transpose,&
24 cp_fm_upper_to_full
25 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
26 cp_fm_struct_release,&
27 cp_fm_struct_type
28 USE cp_fm_types, ONLY: cp_fm_create,&
29 cp_fm_get_info,&
30 cp_fm_release,&
31 cp_fm_to_fm,&
32 cp_fm_type
33 USE cp_log_handling, ONLY: cp_get_default_logger,&
34 cp_logger_get_default_unit_nr,&
35 cp_logger_type
36 USE dbcsr_api, ONLY: dbcsr_copy,&
37 dbcsr_multiply,&
38 dbcsr_release,&
39 dbcsr_scale_by_vector,&
40 dbcsr_set,&
41 dbcsr_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_block, ncol_global, &
219 nrow_block, nrow_global
220 REAL(kind=dp), DIMENSION(:), POINTER :: occupation_numbers, scaling_factor
221 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
222 TYPE(cp_fm_type) :: gradient, h_block, h_block_t, &
223 weighted_vectors
224 TYPE(cp_logger_type), POINTER :: logger
225
226 CALL timeset(routinen, handle)
227 NULLIFY (fm_struct_tmp)
228
229 CALL cp_fm_get_info(matrix=mo_set%mo_coeff, &
230 ncol_global=ncol_global, &
231 nrow_global=nrow_global, &
232 nrow_block=nrow_block, &
233 ncol_block=ncol_block)
234
235 CALL cp_fm_create(weighted_vectors, mo_set%mo_coeff%matrix_struct, "weighted_vectors")
236 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol_global, ncol_global=ncol_global, &
237 para_env=mo_set%mo_coeff%matrix_struct%para_env, &
238 context=mo_set%mo_coeff%matrix_struct%context)
239 CALL cp_fm_create(h_block, fm_struct_tmp, name="h block")
240 IF (do_symm) CALL cp_fm_create(h_block_t, fm_struct_tmp, name="h block t")
241 CALL cp_fm_struct_release(fm_struct_tmp)
242
243 CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers)
244 ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
245 scaling_factor = 2.0_dp*occupation_numbers
246 CALL copy_dbcsr_to_fm(mo_deriv, weighted_vectors)
247 CALL cp_fm_column_scale(weighted_vectors, scaling_factor)
248 DEALLOCATE (scaling_factor)
249
250 ! the convention seems to require the half here, the factor of two is presumably taken care of
251 ! internally in qs_core_hamiltonian
252 CALL parallel_gemm('T', 'N', ncol_global, ncol_global, nrow_global, 0.5_dp, &
253 mo_set%mo_coeff, weighted_vectors, 0.0_dp, h_block)
254
255 IF (do_symm) THEN
256 ! at the minimum things are anyway symmetric, but numerically it might not be the case
257 ! needs some investigation to find out if using this is better
258 CALL cp_fm_transpose(h_block, h_block_t)
259 CALL cp_fm_scale_and_add(0.5_dp, h_block, 0.5_dp, h_block_t)
260 END IF
261
262 ! this could overwrite the mo_derivs to save the weighted_vectors
263 CALL parallel_gemm('N', 'N', nrow_global, ncol_global, ncol_global, 1.0_dp, &
264 mo_set%mo_coeff, h_block, 0.0_dp, weighted_vectors)
265
266 CALL dbcsr_set(w_matrix, 0.0_dp)
267 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=w_matrix, &
268 matrix_v=mo_set%mo_coeff, &
269 matrix_g=weighted_vectors, &
270 ncol=mo_set%homo)
271
272 IF (check_gradient) THEN
273 CALL cp_fm_create(gradient, mo_set%mo_coeff%matrix_struct, "gradient")
274 CALL cp_dbcsr_sm_fm_multiply(s_matrix, weighted_vectors, &
275 gradient, ncol_global)
276
277 ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
278 scaling_factor = 2.0_dp*occupation_numbers
279 CALL copy_dbcsr_to_fm(mo_deriv, weighted_vectors)
280 CALL cp_fm_column_scale(weighted_vectors, scaling_factor)
281 DEALLOCATE (scaling_factor)
282
283 logger => cp_get_default_logger()
284 IF (logger%para_env%is_source()) THEN
285 iounit = cp_logger_get_default_unit_nr(logger, local=.true.)
286 WRITE (iounit, *) " maxabs difference ", &
287 maxval(abs(weighted_vectors%local_data - 2.0_dp*gradient%local_data))
288 END IF
289 CALL cp_fm_release(gradient)
290 END IF
291
292 IF (do_symm) CALL cp_fm_release(h_block_t)
293 CALL cp_fm_release(weighted_vectors)
294 CALL cp_fm_release(h_block)
295
296 CALL timestop(handle)
297
298 END SUBROUTINE calculate_w_matrix_ot
299
300! **************************************************************************************************
301!> \brief Calculate the energy-weighted density matrix W if ROKS is active.
302!> The W matrix is returned in matrix_w.
303!> \param mo_set ...
304!> \param matrix_ks ...
305!> \param matrix_p ...
306!> \param matrix_w ...
307!> \author 04.05.06,MK
308! **************************************************************************************************
309 SUBROUTINE calculate_w_matrix_roks(mo_set, matrix_ks, matrix_p, matrix_w)
310 TYPE(mo_set_type), INTENT(IN) :: mo_set
311 TYPE(dbcsr_type), POINTER :: matrix_ks, matrix_p, matrix_w
312
313 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_roks'
314
315 INTEGER :: handle, nao
316 TYPE(cp_blacs_env_type), POINTER :: context
317 TYPE(cp_fm_struct_type), POINTER :: fm_struct
318 TYPE(cp_fm_type) :: ks, p, work
319 TYPE(cp_fm_type), POINTER :: c
320 TYPE(mp_para_env_type), POINTER :: para_env
321
322 CALL timeset(routinen, handle)
323
324 NULLIFY (context)
325 NULLIFY (fm_struct)
326 NULLIFY (para_env)
327
328 CALL get_mo_set(mo_set=mo_set, mo_coeff=c)
329 CALL cp_fm_get_info(c, context=context, nrow_global=nao, para_env=para_env)
330 CALL cp_fm_struct_create(fm_struct, context=context, nrow_global=nao, &
331 ncol_global=nao, para_env=para_env)
332 CALL cp_fm_create(ks, fm_struct, name="Kohn-Sham matrix")
333 CALL cp_fm_create(p, fm_struct, name="Density matrix")
334 CALL cp_fm_create(work, fm_struct, name="Work matrix")
335 CALL cp_fm_struct_release(fm_struct)
336 CALL copy_dbcsr_to_fm(matrix_ks, ks)
337 CALL copy_dbcsr_to_fm(matrix_p, p)
338 CALL cp_fm_upper_to_full(p, work)
339 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ks, p, 0.0_dp, work)
340 CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, p, work, 0.0_dp, ks)
341 CALL dbcsr_set(matrix_w, 0.0_dp)
342 CALL copy_fm_to_dbcsr(ks, matrix_w, keep_sparsity=.true.)
343
344 CALL cp_fm_release(work)
345 CALL cp_fm_release(p)
346 CALL cp_fm_release(ks)
347
348 CALL timestop(handle)
349
350 END SUBROUTINE calculate_w_matrix_roks
351
352! **************************************************************************************************
353!> \brief Calculate the response W matrix from the MO eigenvectors, MO eigenvalues,
354!> and the MO occupation numbers. Only works if they are eigenstates
355!> \param mo_set type containing the full matrix of the MO and the eigenvalues
356!> \param psi1 response orbitals
357!> \param ks_matrix Kohn-Sham sparse matrix
358!> \param w_matrix sparse matrix
359!> \par History
360!> adapted from calculate_w_matrix_1
361!> \author JGH
362! **************************************************************************************************
363 SUBROUTINE calculate_wz_matrix(mo_set, psi1, ks_matrix, w_matrix)
364
365 TYPE(mo_set_type), INTENT(IN) :: mo_set
366 TYPE(cp_fm_type), INTENT(IN) :: psi1
367 TYPE(dbcsr_type), POINTER :: ks_matrix, w_matrix
368
369 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_wz_matrix'
370
371 INTEGER :: handle, ncol, nocc, nrow
372 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
373 TYPE(cp_fm_type) :: ksmat, scrv
374
375 CALL timeset(routinen, handle)
376
377! CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow)
378! CALL cp_fm_create(scrv, mo_set%mo_coeff%matrix_struct, "scr vectors")
379! CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
380! para_env=mo_set%mo_coeff%matrix_struct%para_env, &
381! context=mo_set%mo_coeff%matrix_struct%context)
382! CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
383! CALL cp_fm_struct_release(fm_struct_tmp)
384! CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, ncol)
385! CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat)
386! CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv)
387! CALL dbcsr_set(w_matrix, 0.0_dp)
388! CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, &
389! ncol=mo_set%homo, symmetry_mode=1)
390! CALL cp_fm_release(scrv)
391! CALL cp_fm_release(ksmat)
392 CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow)
393 nocc = mo_set%homo
394 CALL cp_fm_create(scrv, mo_set%mo_coeff%matrix_struct, "scr vectors")
395 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nocc, ncol_global=nocc, &
396 para_env=mo_set%mo_coeff%matrix_struct%para_env, &
397 context=mo_set%mo_coeff%matrix_struct%context)
398 CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
399 CALL cp_fm_struct_release(fm_struct_tmp)
400 CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, nocc)
401 CALL parallel_gemm("T", "N", nocc, nocc, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat)
402 CALL parallel_gemm("N", "N", nrow, nocc, nocc, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv)
403 CALL dbcsr_set(w_matrix, 0.0_dp)
404 CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, ncol=nocc, symmetry_mode=1)
405 CALL cp_fm_release(scrv)
406 CALL cp_fm_release(ksmat)
407
408 CALL timestop(handle)
409
410 END SUBROUTINE calculate_wz_matrix
411
412! **************************************************************************************************
413!> \brief Calculate the Wz matrix from the MO eigenvectors, MO eigenvalues,
414!> and the MO occupation numbers. Only works if they are eigenstates
415!> \param c0vec ...
416!> \param hzm ...
417!> \param w_matrix sparse matrix
418!> \param focc ...
419!> \param nocc ...
420!> \par History
421!> adapted from calculate_w_matrix_1
422!> \author JGH
423! **************************************************************************************************
424 SUBROUTINE calculate_whz_matrix(c0vec, hzm, w_matrix, focc, nocc)
425
426 TYPE(cp_fm_type), INTENT(IN) :: c0vec
427 TYPE(dbcsr_type), POINTER :: hzm, w_matrix
428 REAL(kind=dp), INTENT(IN) :: focc
429 INTEGER, INTENT(IN) :: nocc
430
431 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_whz_matrix'
432
433 INTEGER :: handle, nao, norb
434 REAL(kind=dp) :: falpha
435 TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_struct_mat
436 TYPE(cp_fm_type) :: chcmat, hcvec
437
438 CALL timeset(routinen, handle)
439
440 falpha = focc
441
442 CALL cp_fm_create(hcvec, c0vec%matrix_struct, "hcvec")
443 CALL cp_fm_get_info(hcvec, matrix_struct=fm_struct, nrow_global=nao, ncol_global=norb)
444 cpassert(nocc <= norb .AND. nocc > 0)
445 norb = nocc
446 CALL cp_fm_struct_create(fm_struct_mat, context=fm_struct%context, nrow_global=norb, &
447 ncol_global=norb, para_env=fm_struct%para_env)
448 CALL cp_fm_create(chcmat, fm_struct_mat)
449 CALL cp_fm_struct_release(fm_struct_mat)
450
451 CALL cp_dbcsr_sm_fm_multiply(hzm, c0vec, hcvec, norb)
452 CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, c0vec, hcvec, 0.0_dp, chcmat)
453 CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, c0vec, chcmat, 0.0_dp, hcvec)
454
455 CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=hcvec, matrix_g=c0vec, ncol=norb, alpha=falpha)
456
457 CALL cp_fm_release(hcvec)
458 CALL cp_fm_release(chcmat)
459
460 CALL timestop(handle)
461
462 END SUBROUTINE calculate_whz_matrix
463
464! **************************************************************************************************
465!> \brief Calculate the excited state W matrix from the MO eigenvectors, KS matrix
466!> \param mos_occ ...
467!> \param xvec ...
468!> \param ks_matrix ...
469!> \param w_matrix ...
470!> \par History
471!> adapted from calculate_wz_matrix
472!> \author JGH
473! **************************************************************************************************
474 SUBROUTINE calculate_wx_matrix(mos_occ, xvec, ks_matrix, w_matrix)
475
476 TYPE(cp_fm_type), INTENT(IN) :: mos_occ, xvec
477 TYPE(dbcsr_type), POINTER :: ks_matrix, w_matrix
478
479 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_wx_matrix'
480
481 INTEGER :: handle, ncol, nrow
482 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
483 TYPE(cp_fm_type) :: ksmat, scrv
484
485 CALL timeset(routinen, handle)
486
487 CALL cp_fm_get_info(matrix=mos_occ, ncol_global=ncol, nrow_global=nrow)
488 CALL cp_fm_create(scrv, mos_occ%matrix_struct, "scr vectors")
489 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
490 para_env=mos_occ%matrix_struct%para_env, &
491 context=mos_occ%matrix_struct%context)
492 CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
493 CALL cp_fm_struct_release(fm_struct_tmp)
494 CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mos_occ, scrv, ncol)
495 CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mos_occ, scrv, 0.0_dp, ksmat)
496 CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, xvec, ksmat, 0.0_dp, scrv)
497 CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=xvec, ncol=ncol, symmetry_mode=1)
498 CALL cp_fm_release(scrv)
499 CALL cp_fm_release(ksmat)
500
501 CALL timestop(handle)
502
503 END SUBROUTINE calculate_wx_matrix
504
505! **************************************************************************************************
506!> \brief Calculate the excited state W matrix from the MO eigenvectors, KS matrix
507!> \param mos_occ ...
508!> \param xvec ...
509!> \param s_matrix ...
510!> \param ks_matrix ...
511!> \param w_matrix ...
512!> \param eval ...
513!> \par History
514!> adapted from calculate_wz_matrix
515!> \author JGH
516! **************************************************************************************************
517 SUBROUTINE calculate_xwx_matrix(mos_occ, xvec, s_matrix, ks_matrix, w_matrix, eval)
518
519 TYPE(cp_fm_type), INTENT(IN) :: mos_occ, xvec
520 TYPE(dbcsr_type), POINTER :: s_matrix, ks_matrix, w_matrix
521 REAL(kind=dp), INTENT(IN) :: eval
522
523 CHARACTER(len=*), PARAMETER :: routinen = 'calculate_xwx_matrix'
524
525 INTEGER :: handle, ncol, nrow
526 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
527 TYPE(cp_fm_type) :: scrv, xsxmat
528
529 CALL timeset(routinen, handle)
530
531 CALL cp_fm_get_info(matrix=mos_occ, ncol_global=ncol, nrow_global=nrow)
532 CALL cp_fm_create(scrv, mos_occ%matrix_struct, "scr vectors")
533 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
534 para_env=mos_occ%matrix_struct%para_env, &
535 context=mos_occ%matrix_struct%context)
536 CALL cp_fm_create(xsxmat, fm_struct_tmp, name="XSX")
537 CALL cp_fm_struct_release(fm_struct_tmp)
538
539 CALL cp_dbcsr_sm_fm_multiply(ks_matrix, xvec, scrv, ncol, 1.0_dp, 0.0_dp)
540 CALL cp_dbcsr_sm_fm_multiply(s_matrix, xvec, scrv, ncol, eval, -1.0_dp)
541 CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, xvec, scrv, 0.0_dp, xsxmat)
542
543 CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mos_occ, xsxmat, 0.0_dp, scrv)
544 CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=mos_occ, ncol=ncol, symmetry_mode=1)
545
546 CALL cp_fm_release(scrv)
547 CALL cp_fm_release(xsxmat)
548
549 CALL timestop(handle)
550
551 END SUBROUTINE calculate_xwx_matrix
552
553END MODULE qs_density_matrices
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
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_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:744
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:208
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:887
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:118
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:1882
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:1710
cp_fm_basic_linalg::cp_fm_upper_to_full
subroutine, public cp_fm_upper_to_full(matrix, work)
given an upper triangular matrix computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1759
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_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:576
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:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:320
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:1016
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
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:425
qs_density_matrices::calculate_wx_matrix
subroutine, public calculate_wx_matrix(mos_occ, xvec, ks_matrix, w_matrix)
Calculate the excited state W matrix from the MO eigenvectors, KS matrix.
Definition qs_density_matrices.F:475
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:364
qs_density_matrices::calculate_xwx_matrix
subroutine, public calculate_xwx_matrix(mos_occ, 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:518
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:397
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_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:116
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:763
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46