(git:ccc2433)
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 ! **************************************************************************************************
14 MODULE 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 
66 CONTAINS
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 
553 END MODULE qs_density_matrices
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