Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \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 183687 : 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 183687 : CALL timeset(routineN, handle)
95 :
96 183687 : my_use_dbcsr = .FALSE.
97 183687 : IF (PRESENT(use_dbcsr)) my_use_dbcsr = use_dbcsr
98 183687 : my_retain_sparsity = .TRUE.
99 183687 : IF (PRESENT(retain_sparsity)) my_retain_sparsity = retain_sparsity
100 183687 : IF (my_use_dbcsr) THEN
101 74466 : 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 183687 : CALL dbcsr_set(density_matrix, 0.0_dp)
107 :
108 183687 : IF (.NOT. mo_set%uniform_occupation) THEN ! not all orbitals 1..homo are equally occupied
109 13020 : 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 13020 : CALL cp_fm_create(fm_tmp, mo_set%mo_coeff%matrix_struct)
119 13020 : CALL cp_fm_to_fm(mo_set%mo_coeff, fm_tmp)
120 13020 : CALL cp_fm_column_scale(fm_tmp, mo_set%occupation_numbers(1:mo_set%homo))
121 13020 : 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 13020 : alpha=alpha)
127 13020 : CALL cp_fm_release(fm_tmp)
128 : END IF
129 : ELSE
130 170667 : 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 74466 : last_k=mo_set%homo)
134 : ELSE
135 96201 : 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 96201 : alpha=alpha)
140 : END IF
141 : END IF
142 :
143 183687 : CALL timestop(handle)
144 :
145 183687 : 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 3466 : 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 3466 : CALL timeset(routineN, handle)
171 :
172 3466 : CALL dbcsr_set(w_matrix, 0.0_dp)
173 3466 : CALL cp_fm_create(weighted_vectors, mo_set%mo_coeff%matrix_struct, "weighted_vectors")
174 3466 : CALL cp_fm_to_fm(mo_set%mo_coeff, weighted_vectors)
175 :
176 : ! scale every column with the occupation
177 10348 : ALLOCATE (eigocc(mo_set%homo))
178 :
179 41842 : DO imo = 1, mo_set%homo
180 41842 : eigocc(imo) = mo_set%eigenvalues(imo)*mo_set%occupation_numbers(imo)
181 : END DO
182 3466 : CALL cp_fm_column_scale(weighted_vectors, eigocc)
183 3466 : 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 3466 : ncol=mo_set%homo)
189 :
190 3466 : CALL cp_fm_release(weighted_vectors)
191 :
192 3466 : CALL timestop(handle)
193 :
194 3466 : 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 2457 : 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 2457 : 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 2457 : CALL timeset(routineN, handle)
227 2457 : 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 2457 : ncol_block=ncol_block)
234 :
235 2457 : 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 2457 : context=mo_set%mo_coeff%matrix_struct%context)
239 2457 : 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 2457 : CALL cp_fm_struct_release(fm_struct_tmp)
242 :
243 2457 : CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers)
244 7335 : ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
245 37530 : scaling_factor = 2.0_dp*occupation_numbers
246 2457 : CALL copy_dbcsr_to_fm(mo_deriv, weighted_vectors)
247 2457 : CALL cp_fm_column_scale(weighted_vectors, scaling_factor)
248 2457 : 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 2457 : 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 2457 : mo_set%mo_coeff, h_block, 0.0_dp, weighted_vectors)
265 :
266 2457 : 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 2457 : 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 2457 : CALL cp_fm_release(weighted_vectors)
294 2457 : CALL cp_fm_release(h_block)
295 :
296 2457 : CALL timestop(handle)
297 :
298 2457 : 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 104 : 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 52 : CALL timeset(routineN, handle)
323 :
324 52 : NULLIFY (context)
325 52 : NULLIFY (fm_struct)
326 52 : NULLIFY (para_env)
327 :
328 52 : CALL get_mo_set(mo_set=mo_set, mo_coeff=c)
329 52 : 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 52 : ncol_global=nao, para_env=para_env)
332 52 : CALL cp_fm_create(ks, fm_struct, name="Kohn-Sham matrix")
333 52 : CALL cp_fm_create(p, fm_struct, name="Density matrix")
334 52 : CALL cp_fm_create(work, fm_struct, name="Work matrix")
335 52 : CALL cp_fm_struct_release(fm_struct)
336 52 : CALL copy_dbcsr_to_fm(matrix_ks, ks)
337 52 : CALL copy_dbcsr_to_fm(matrix_p, p)
338 52 : CALL cp_fm_upper_to_full(p, work)
339 52 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ks, p, 0.0_dp, work)
340 52 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, p, work, 0.0_dp, ks)
341 52 : CALL dbcsr_set(matrix_w, 0.0_dp)
342 52 : CALL copy_fm_to_dbcsr(ks, matrix_w, keep_sparsity=.TRUE.)
343 :
344 52 : CALL cp_fm_release(work)
345 52 : CALL cp_fm_release(p)
346 52 : CALL cp_fm_release(ks)
347 :
348 52 : CALL timestop(handle)
349 :
350 52 : 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 2048 : 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 1024 : 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 1024 : CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow)
393 1024 : nocc = mo_set%homo
394 1024 : 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 1024 : context=mo_set%mo_coeff%matrix_struct%context)
398 1024 : CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
399 1024 : CALL cp_fm_struct_release(fm_struct_tmp)
400 1024 : CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, nocc)
401 1024 : CALL parallel_gemm("T", "N", nocc, nocc, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat)
402 1024 : CALL parallel_gemm("N", "N", nrow, nocc, nocc, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv)
403 1024 : CALL dbcsr_set(w_matrix, 0.0_dp)
404 1024 : CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, ncol=nocc, symmetry_mode=1)
405 1024 : CALL cp_fm_release(scrv)
406 1024 : CALL cp_fm_release(ksmat)
407 :
408 1024 : CALL timestop(handle)
409 :
410 1024 : 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 2980 : 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 1490 : CALL timeset(routineN, handle)
439 :
440 1490 : falpha = focc
441 :
442 1490 : CALL cp_fm_create(hcvec, c0vec%matrix_struct, "hcvec")
443 1490 : CALL cp_fm_get_info(hcvec, matrix_struct=fm_struct, nrow_global=nao, ncol_global=norb)
444 1490 : CPASSERT(nocc <= norb .AND. nocc > 0)
445 1490 : norb = nocc
446 : CALL cp_fm_struct_create(fm_struct_mat, context=fm_struct%context, nrow_global=norb, &
447 1490 : ncol_global=norb, para_env=fm_struct%para_env)
448 1490 : CALL cp_fm_create(chcmat, fm_struct_mat)
449 1490 : CALL cp_fm_struct_release(fm_struct_mat)
450 :
451 1490 : CALL cp_dbcsr_sm_fm_multiply(hzm, c0vec, hcvec, norb)
452 1490 : CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, c0vec, hcvec, 0.0_dp, chcmat)
453 1490 : CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, c0vec, chcmat, 0.0_dp, hcvec)
454 :
455 1490 : CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=hcvec, matrix_g=c0vec, ncol=norb, alpha=falpha)
456 :
457 1490 : CALL cp_fm_release(hcvec)
458 1490 : CALL cp_fm_release(chcmat)
459 :
460 1490 : CALL timestop(handle)
461 :
462 1490 : 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 1324 : 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 662 : CALL timeset(routineN, handle)
486 :
487 662 : CALL cp_fm_get_info(matrix=mos_occ, ncol_global=ncol, nrow_global=nrow)
488 662 : 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 662 : context=mos_occ%matrix_struct%context)
492 662 : CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
493 662 : CALL cp_fm_struct_release(fm_struct_tmp)
494 662 : CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mos_occ, scrv, ncol)
495 662 : CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mos_occ, scrv, 0.0_dp, ksmat)
496 662 : CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, xvec, ksmat, 0.0_dp, scrv)
497 662 : CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=xvec, ncol=ncol, symmetry_mode=1)
498 662 : CALL cp_fm_release(scrv)
499 662 : CALL cp_fm_release(ksmat)
500 :
501 662 : CALL timestop(handle)
502 :
503 662 : 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 1324 : 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 662 : CALL timeset(routineN, handle)
530 :
531 662 : CALL cp_fm_get_info(matrix=mos_occ, ncol_global=ncol, nrow_global=nrow)
532 662 : 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 662 : context=mos_occ%matrix_struct%context)
536 662 : CALL cp_fm_create(xsxmat, fm_struct_tmp, name="XSX")
537 662 : CALL cp_fm_struct_release(fm_struct_tmp)
538 :
539 662 : CALL cp_dbcsr_sm_fm_multiply(ks_matrix, xvec, scrv, ncol, 1.0_dp, 0.0_dp)
540 662 : CALL cp_dbcsr_sm_fm_multiply(s_matrix, xvec, scrv, ncol, eval, -1.0_dp)
541 662 : CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, xvec, scrv, 0.0_dp, xsxmat)
542 :
543 662 : CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mos_occ, xsxmat, 0.0_dp, scrv)
544 662 : CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=mos_occ, ncol=ncol, symmetry_mode=1)
545 :
546 662 : CALL cp_fm_release(scrv)
547 662 : CALL cp_fm_release(xsxmat)
548 :
549 662 : CALL timestop(handle)
550 :
551 662 : END SUBROUTINE calculate_xwx_matrix
552 :
553 : END MODULE qs_density_matrices
|
