Line data Source code
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 : ! **************************************************************************************************
14 : MODULE 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 :
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 199026 : 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 199026 : CALL timeset(routineN, handle)
95 :
96 199026 : my_use_dbcsr = .FALSE.
97 199026 : IF (PRESENT(use_dbcsr)) my_use_dbcsr = use_dbcsr
98 199026 : my_retain_sparsity = .TRUE.
99 199026 : IF (PRESENT(retain_sparsity)) my_retain_sparsity = retain_sparsity
100 199026 : IF (my_use_dbcsr) THEN
101 84209 : IF (.NOT. ASSOCIATED(mo_set%mo_coeff_b)) THEN
102 0 : CPABORT("mo_coeff_b NOT ASSOCIATED")
103 : END IF
104 : END IF
105 :
106 199026 : CALL dbcsr_set(density_matrix, 0.0_dp)
107 :
108 199026 : IF (.NOT. mo_set%uniform_occupation) THEN ! not all orbitals 1..homo are equally occupied
109 14100 : IF (my_use_dbcsr) THEN
110 0 : 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 0 : 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 0 : last_k=mo_set%homo)
116 0 : CALL dbcsr_release(dbcsr_tmp)
117 : ELSE
118 14100 : CALL cp_fm_create(fm_tmp, mo_set%mo_coeff%matrix_struct)
119 14100 : CALL cp_fm_to_fm(mo_set%mo_coeff, fm_tmp)
120 14100 : CALL cp_fm_column_scale(fm_tmp, mo_set%occupation_numbers(1:mo_set%homo))
121 14100 : 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 14100 : alpha=alpha)
127 14100 : CALL cp_fm_release(fm_tmp)
128 : END IF
129 : ELSE
130 184926 : 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 84209 : last_k=mo_set%homo)
134 : ELSE
135 100717 : 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 100717 : alpha=alpha)
140 : END IF
141 : END IF
142 :
143 199026 : CALL timestop(handle)
144 :
145 199026 : 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 3602 : 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 3602 : CALL timeset(routineN, handle)
171 :
172 3602 : CALL dbcsr_set(w_matrix, 0.0_dp)
173 3602 : CALL cp_fm_create(weighted_vectors, mo_set%mo_coeff%matrix_struct, "weighted_vectors")
174 3602 : CALL cp_fm_to_fm(mo_set%mo_coeff, weighted_vectors)
175 :
176 : ! scale every column with the occupation
177 10756 : ALLOCATE (eigocc(mo_set%homo))
178 :
179 43308 : DO imo = 1, mo_set%homo
180 43308 : eigocc(imo) = mo_set%eigenvalues(imo)*mo_set%occupation_numbers(imo)
181 : END DO
182 3602 : CALL cp_fm_column_scale(weighted_vectors, eigocc)
183 3602 : 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 3602 : ncol=mo_set%homo)
189 :
190 3602 : CALL cp_fm_release(weighted_vectors)
191 :
192 3602 : CALL timestop(handle)
193 :
194 7204 : 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 2705 : 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 2705 : 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 2705 : CALL timeset(routineN, handle)
226 2705 : NULLIFY (fm_struct_tmp)
227 :
228 : CALL cp_fm_get_info(matrix=mo_set%mo_coeff, &
229 : ncol_global=ncol_global, &
230 2705 : nrow_global=nrow_global)
231 :
232 2705 : 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 2705 : context=mo_set%mo_coeff%matrix_struct%context)
236 2705 : 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 2705 : CALL cp_fm_struct_release(fm_struct_tmp)
239 :
240 2705 : CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers)
241 8079 : ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
242 41758 : scaling_factor = 2.0_dp*occupation_numbers
243 2705 : CALL copy_dbcsr_to_fm(mo_deriv, weighted_vectors)
244 2705 : CALL cp_fm_column_scale(weighted_vectors, scaling_factor)
245 2705 : 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 2705 : 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 2705 : mo_set%mo_coeff, h_block, 0.0_dp, weighted_vectors)
262 :
263 2705 : 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 2705 : 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 2705 : CALL cp_fm_release(weighted_vectors)
291 2705 : CALL cp_fm_release(h_block)
292 :
293 2705 : CALL timestop(handle)
294 :
295 2705 : 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 260 : 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 52 : CALL timeset(routineN, handle)
320 :
321 52 : NULLIFY (context)
322 52 : NULLIFY (fm_struct)
323 52 : NULLIFY (para_env)
324 :
325 52 : CALL get_mo_set(mo_set=mo_set, mo_coeff=c)
326 52 : 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 52 : ncol_global=nao, para_env=para_env)
329 52 : CALL cp_fm_create(ks, fm_struct, name="Kohn-Sham matrix")
330 52 : CALL cp_fm_create(p, fm_struct, name="Density matrix")
331 52 : CALL cp_fm_create(work, fm_struct, name="Work matrix")
332 52 : CALL cp_fm_struct_release(fm_struct)
333 52 : CALL copy_dbcsr_to_fm(matrix_ks, ks)
334 52 : CALL copy_dbcsr_to_fm(matrix_p, p)
335 52 : CALL cp_fm_uplo_to_full(p, work)
336 52 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ks, p, 0.0_dp, work)
337 52 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, p, work, 0.0_dp, ks)
338 52 : CALL dbcsr_set(matrix_w, 0.0_dp)
339 52 : CALL copy_fm_to_dbcsr(ks, matrix_w, keep_sparsity=.TRUE.)
340 :
341 52 : CALL cp_fm_release(work)
342 52 : CALL cp_fm_release(p)
343 52 : CALL cp_fm_release(ks)
344 :
345 52 : CALL timestop(handle)
346 :
347 52 : 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 3258 : 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 1086 : CALL timeset(routineN, handle)
373 :
374 1086 : CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow)
375 1086 : nocc = mo_set%homo
376 1086 : 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 1086 : context=mo_set%mo_coeff%matrix_struct%context)
380 1086 : CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
381 1086 : CALL cp_fm_struct_release(fm_struct_tmp)
382 1086 : CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, nocc)
383 1086 : CALL parallel_gemm("T", "N", nocc, nocc, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat)
384 1086 : CALL parallel_gemm("N", "N", nrow, nocc, nocc, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv)
385 1086 : CALL dbcsr_set(w_matrix, 0.0_dp)
386 1086 : CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, ncol=nocc, symmetry_mode=1)
387 1086 : CALL cp_fm_release(scrv)
388 1086 : CALL cp_fm_release(ksmat)
389 :
390 1086 : CALL timestop(handle)
391 :
392 1086 : 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 4698 : 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 1566 : CALL timeset(routineN, handle)
421 :
422 1566 : falpha = focc
423 :
424 1566 : CALL cp_fm_create(hcvec, c0vec%matrix_struct, name="hcvec", set_zero=.TRUE.)
425 1566 : CALL cp_fm_get_info(hcvec, matrix_struct=fm_struct, nrow_global=nao, ncol_global=norb)
426 1566 : CPASSERT(nocc <= norb .AND. nocc > 0)
427 1566 : norb = nocc
428 : CALL cp_fm_struct_create(fm_struct_mat, context=fm_struct%context, nrow_global=norb, &
429 1566 : ncol_global=norb, para_env=fm_struct%para_env)
430 1566 : CALL cp_fm_create(chcmat, fm_struct_mat)
431 1566 : CALL cp_fm_struct_release(fm_struct_mat)
432 :
433 1566 : CALL cp_dbcsr_sm_fm_multiply(hzm, c0vec, hcvec, norb)
434 1566 : CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, c0vec, hcvec, 0.0_dp, chcmat)
435 1566 : CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, c0vec, chcmat, 0.0_dp, hcvec)
436 :
437 1566 : CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=hcvec, matrix_g=c0vec, ncol=norb, alpha=falpha)
438 :
439 1566 : CALL cp_fm_release(hcvec)
440 1566 : CALL cp_fm_release(chcmat)
441 :
442 1566 : CALL timestop(handle)
443 :
444 1566 : 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 2028 : 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 676 : CALL timeset(routineN, handle)
468 :
469 676 : CALL cp_fm_get_info(matrix=mos_active, ncol_global=ncol, nrow_global=nrow)
470 676 : 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 676 : context=mos_active%matrix_struct%context)
474 676 : CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
475 676 : CALL cp_fm_struct_release(fm_struct_tmp)
476 676 : CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mos_active, scrv, ncol)
477 676 : CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mos_active, scrv, 0.0_dp, ksmat)
478 676 : CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, xvec, ksmat, 0.0_dp, scrv)
479 676 : CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=xvec, ncol=ncol, symmetry_mode=1)
480 676 : CALL cp_fm_release(scrv)
481 676 : CALL cp_fm_release(ksmat)
482 :
483 676 : CALL timestop(handle)
484 :
485 676 : 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 2028 : 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 676 : CALL timeset(routineN, handle)
512 :
513 676 : CALL cp_fm_get_info(matrix=mos_active, ncol_global=ncol, nrow_global=nrow)
514 676 : 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 676 : context=mos_active%matrix_struct%context)
518 676 : CALL cp_fm_create(xsxmat, fm_struct_tmp, name="XSX")
519 676 : CALL cp_fm_struct_release(fm_struct_tmp)
520 :
521 676 : CALL cp_dbcsr_sm_fm_multiply(ks_matrix, xvec, scrv, ncol, 1.0_dp, 0.0_dp)
522 676 : CALL cp_dbcsr_sm_fm_multiply(s_matrix, xvec, scrv, ncol, eval, -1.0_dp)
523 676 : CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, xvec, scrv, 0.0_dp, xsxmat)
524 :
525 676 : CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mos_active, xsxmat, 0.0_dp, scrv)
526 676 : CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=mos_active, ncol=ncol, symmetry_mode=1)
527 :
528 676 : CALL cp_fm_release(scrv)
529 676 : CALL cp_fm_release(xsxmat)
530 :
531 676 : CALL timestop(handle)
532 :
533 676 : END SUBROUTINE calculate_xwx_matrix
534 :
535 : END MODULE qs_density_matrices
|
