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_api, ONLY: dbcsr_copy,&
17 : dbcsr_multiply,&
18 : dbcsr_release,&
19 : dbcsr_scale_by_vector,&
20 : dbcsr_set,&
21 : dbcsr_type
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_upper_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 184235 : 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 184235 : CALL timeset(routineN, handle)
95 :
96 184235 : my_use_dbcsr = .FALSE.
97 184235 : IF (PRESENT(use_dbcsr)) my_use_dbcsr = use_dbcsr
98 184235 : my_retain_sparsity = .TRUE.
99 184235 : IF (PRESENT(retain_sparsity)) my_retain_sparsity = retain_sparsity
100 184235 : IF (my_use_dbcsr) THEN
101 72710 : 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 184235 : CALL dbcsr_set(density_matrix, 0.0_dp)
107 :
108 184235 : IF (.NOT. mo_set%uniform_occupation) THEN ! not all orbitals 1..homo are equally occupied
109 14292 : 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 14292 : CALL cp_fm_create(fm_tmp, mo_set%mo_coeff%matrix_struct)
119 14292 : CALL cp_fm_to_fm(mo_set%mo_coeff, fm_tmp)
120 14292 : CALL cp_fm_column_scale(fm_tmp, mo_set%occupation_numbers(1:mo_set%homo))
121 14292 : 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 14292 : alpha=alpha)
127 14292 : CALL cp_fm_release(fm_tmp)
128 : END IF
129 : ELSE
130 169943 : 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 72710 : last_k=mo_set%homo)
134 : ELSE
135 97233 : 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 97233 : alpha=alpha)
140 : END IF
141 : END IF
142 :
143 184235 : CALL timestop(handle)
144 :
145 184235 : 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 3502 : 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 3502 : CALL timeset(routineN, handle)
171 :
172 3502 : CALL dbcsr_set(w_matrix, 0.0_dp)
173 3502 : CALL cp_fm_create(weighted_vectors, mo_set%mo_coeff%matrix_struct, "weighted_vectors")
174 3502 : CALL cp_fm_to_fm(mo_set%mo_coeff, weighted_vectors)
175 :
176 : ! scale every column with the occupation
177 10456 : ALLOCATE (eigocc(mo_set%homo))
178 :
179 42176 : DO imo = 1, mo_set%homo
180 42176 : eigocc(imo) = mo_set%eigenvalues(imo)*mo_set%occupation_numbers(imo)
181 : END DO
182 3502 : CALL cp_fm_column_scale(weighted_vectors, eigocc)
183 3502 : 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 3502 : ncol=mo_set%homo)
189 :
190 3502 : CALL cp_fm_release(weighted_vectors)
191 :
192 3502 : CALL timestop(handle)
193 :
194 7004 : 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 2471 : 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 2471 : 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 2471 : CALL timeset(routineN, handle)
226 2471 : NULLIFY (fm_struct_tmp)
227 :
228 : CALL cp_fm_get_info(matrix=mo_set%mo_coeff, &
229 : ncol_global=ncol_global, &
230 2471 : nrow_global=nrow_global)
231 :
232 2471 : 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 2471 : context=mo_set%mo_coeff%matrix_struct%context)
236 2471 : CALL cp_fm_create(h_block, fm_struct_tmp, name="h block")
237 : IF (do_symm) CALL cp_fm_create(h_block_t, fm_struct_tmp, name="h block t")
238 2471 : CALL cp_fm_struct_release(fm_struct_tmp)
239 :
240 2471 : CALL get_mo_set(mo_set=mo_set, occupation_numbers=occupation_numbers)
241 7377 : ALLOCATE (scaling_factor(SIZE(occupation_numbers)))
242 37674 : scaling_factor = 2.0_dp*occupation_numbers
243 2471 : CALL copy_dbcsr_to_fm(mo_deriv, weighted_vectors)
244 2471 : CALL cp_fm_column_scale(weighted_vectors, scaling_factor)
245 2471 : 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 2471 : 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 2471 : mo_set%mo_coeff, h_block, 0.0_dp, weighted_vectors)
262 :
263 2471 : 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 2471 : 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 2471 : CALL cp_fm_release(weighted_vectors)
291 2471 : CALL cp_fm_release(h_block)
292 :
293 2471 : CALL timestop(handle)
294 :
295 4942 : 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_upper_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 4112 : 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 1028 : CALL timeset(routineN, handle)
373 :
374 : ! CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow)
375 : ! CALL cp_fm_create(scrv, mo_set%mo_coeff%matrix_struct, "scr vectors")
376 : ! CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
377 : ! para_env=mo_set%mo_coeff%matrix_struct%para_env, &
378 : ! context=mo_set%mo_coeff%matrix_struct%context)
379 : ! CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
380 : ! CALL cp_fm_struct_release(fm_struct_tmp)
381 : ! CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, ncol)
382 : ! CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat)
383 : ! CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv)
384 : ! CALL dbcsr_set(w_matrix, 0.0_dp)
385 : ! CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, &
386 : ! ncol=mo_set%homo, symmetry_mode=1)
387 : ! CALL cp_fm_release(scrv)
388 : ! CALL cp_fm_release(ksmat)
389 1028 : CALL cp_fm_get_info(matrix=mo_set%mo_coeff, ncol_global=ncol, nrow_global=nrow)
390 1028 : nocc = mo_set%homo
391 1028 : CALL cp_fm_create(scrv, mo_set%mo_coeff%matrix_struct, "scr vectors")
392 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nocc, ncol_global=nocc, &
393 : para_env=mo_set%mo_coeff%matrix_struct%para_env, &
394 1028 : context=mo_set%mo_coeff%matrix_struct%context)
395 1028 : CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
396 1028 : CALL cp_fm_struct_release(fm_struct_tmp)
397 1028 : CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mo_set%mo_coeff, scrv, nocc)
398 1028 : CALL parallel_gemm("T", "N", nocc, nocc, nrow, 1.0_dp, mo_set%mo_coeff, scrv, 0.0_dp, ksmat)
399 1028 : CALL parallel_gemm("N", "N", nrow, nocc, nocc, 1.0_dp, mo_set%mo_coeff, ksmat, 0.0_dp, scrv)
400 1028 : CALL dbcsr_set(w_matrix, 0.0_dp)
401 1028 : CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=psi1, ncol=nocc, symmetry_mode=1)
402 1028 : CALL cp_fm_release(scrv)
403 1028 : CALL cp_fm_release(ksmat)
404 :
405 1028 : CALL timestop(handle)
406 :
407 1028 : END SUBROUTINE calculate_wz_matrix
408 :
409 : ! **************************************************************************************************
410 : !> \brief Calculate the Wz matrix from the MO eigenvectors, MO eigenvalues,
411 : !> and the MO occupation numbers. Only works if they are eigenstates
412 : !> \param c0vec ...
413 : !> \param hzm ...
414 : !> \param w_matrix sparse matrix
415 : !> \param focc ...
416 : !> \param nocc ...
417 : !> \par History
418 : !> adapted from calculate_w_matrix_1
419 : !> \author JGH
420 : ! **************************************************************************************************
421 6000 : SUBROUTINE calculate_whz_matrix(c0vec, hzm, w_matrix, focc, nocc)
422 :
423 : TYPE(cp_fm_type), INTENT(IN) :: c0vec
424 : TYPE(dbcsr_type), POINTER :: hzm, w_matrix
425 : REAL(KIND=dp), INTENT(IN) :: focc
426 : INTEGER, INTENT(IN) :: nocc
427 :
428 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_whz_matrix'
429 :
430 : INTEGER :: handle, nao, norb
431 : REAL(KIND=dp) :: falpha
432 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, fm_struct_mat
433 : TYPE(cp_fm_type) :: chcmat, hcvec
434 :
435 1500 : CALL timeset(routineN, handle)
436 :
437 1500 : falpha = focc
438 :
439 1500 : CALL cp_fm_create(hcvec, c0vec%matrix_struct, "hcvec")
440 1500 : CALL cp_fm_get_info(hcvec, matrix_struct=fm_struct, nrow_global=nao, ncol_global=norb)
441 1500 : CPASSERT(nocc <= norb .AND. nocc > 0)
442 1500 : norb = nocc
443 : CALL cp_fm_struct_create(fm_struct_mat, context=fm_struct%context, nrow_global=norb, &
444 1500 : ncol_global=norb, para_env=fm_struct%para_env)
445 1500 : CALL cp_fm_create(chcmat, fm_struct_mat)
446 1500 : CALL cp_fm_struct_release(fm_struct_mat)
447 :
448 1500 : CALL cp_dbcsr_sm_fm_multiply(hzm, c0vec, hcvec, norb)
449 1500 : CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, c0vec, hcvec, 0.0_dp, chcmat)
450 1500 : CALL parallel_gemm("N", "N", nao, norb, norb, 1.0_dp, c0vec, chcmat, 0.0_dp, hcvec)
451 :
452 1500 : CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=hcvec, matrix_g=c0vec, ncol=norb, alpha=falpha)
453 :
454 1500 : CALL cp_fm_release(hcvec)
455 1500 : CALL cp_fm_release(chcmat)
456 :
457 1500 : CALL timestop(handle)
458 :
459 1500 : END SUBROUTINE calculate_whz_matrix
460 :
461 : ! **************************************************************************************************
462 : !> \brief Calculate the excited state W matrix from the MO eigenvectors, KS matrix
463 : !> \param mos_occ ...
464 : !> \param xvec ...
465 : !> \param ks_matrix ...
466 : !> \param w_matrix ...
467 : !> \par History
468 : !> adapted from calculate_wz_matrix
469 : !> \author JGH
470 : ! **************************************************************************************************
471 2648 : SUBROUTINE calculate_wx_matrix(mos_occ, xvec, ks_matrix, w_matrix)
472 :
473 : TYPE(cp_fm_type), INTENT(IN) :: mos_occ, xvec
474 : TYPE(dbcsr_type), POINTER :: ks_matrix, w_matrix
475 :
476 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_wx_matrix'
477 :
478 : INTEGER :: handle, ncol, nrow
479 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
480 : TYPE(cp_fm_type) :: ksmat, scrv
481 :
482 662 : CALL timeset(routineN, handle)
483 :
484 662 : CALL cp_fm_get_info(matrix=mos_occ, ncol_global=ncol, nrow_global=nrow)
485 662 : CALL cp_fm_create(scrv, mos_occ%matrix_struct, "scr vectors")
486 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
487 : para_env=mos_occ%matrix_struct%para_env, &
488 662 : context=mos_occ%matrix_struct%context)
489 662 : CALL cp_fm_create(ksmat, fm_struct_tmp, name="KS")
490 662 : CALL cp_fm_struct_release(fm_struct_tmp)
491 662 : CALL cp_dbcsr_sm_fm_multiply(ks_matrix, mos_occ, scrv, ncol)
492 662 : CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, mos_occ, scrv, 0.0_dp, ksmat)
493 662 : CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, xvec, ksmat, 0.0_dp, scrv)
494 662 : CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=xvec, ncol=ncol, symmetry_mode=1)
495 662 : CALL cp_fm_release(scrv)
496 662 : CALL cp_fm_release(ksmat)
497 :
498 662 : CALL timestop(handle)
499 :
500 662 : END SUBROUTINE calculate_wx_matrix
501 :
502 : ! **************************************************************************************************
503 : !> \brief Calculate the excited state W matrix from the MO eigenvectors, KS matrix
504 : !> \param mos_occ ...
505 : !> \param xvec ...
506 : !> \param s_matrix ...
507 : !> \param ks_matrix ...
508 : !> \param w_matrix ...
509 : !> \param eval ...
510 : !> \par History
511 : !> adapted from calculate_wz_matrix
512 : !> \author JGH
513 : ! **************************************************************************************************
514 2648 : SUBROUTINE calculate_xwx_matrix(mos_occ, xvec, s_matrix, ks_matrix, w_matrix, eval)
515 :
516 : TYPE(cp_fm_type), INTENT(IN) :: mos_occ, xvec
517 : TYPE(dbcsr_type), POINTER :: s_matrix, ks_matrix, w_matrix
518 : REAL(KIND=dp), INTENT(IN) :: eval
519 :
520 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_xwx_matrix'
521 :
522 : INTEGER :: handle, ncol, nrow
523 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
524 : TYPE(cp_fm_type) :: scrv, xsxmat
525 :
526 662 : CALL timeset(routineN, handle)
527 :
528 662 : CALL cp_fm_get_info(matrix=mos_occ, ncol_global=ncol, nrow_global=nrow)
529 662 : CALL cp_fm_create(scrv, mos_occ%matrix_struct, "scr vectors")
530 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
531 : para_env=mos_occ%matrix_struct%para_env, &
532 662 : context=mos_occ%matrix_struct%context)
533 662 : CALL cp_fm_create(xsxmat, fm_struct_tmp, name="XSX")
534 662 : CALL cp_fm_struct_release(fm_struct_tmp)
535 :
536 662 : CALL cp_dbcsr_sm_fm_multiply(ks_matrix, xvec, scrv, ncol, 1.0_dp, 0.0_dp)
537 662 : CALL cp_dbcsr_sm_fm_multiply(s_matrix, xvec, scrv, ncol, eval, -1.0_dp)
538 662 : CALL parallel_gemm("T", "N", ncol, ncol, nrow, 1.0_dp, xvec, scrv, 0.0_dp, xsxmat)
539 :
540 662 : CALL parallel_gemm("N", "N", nrow, ncol, ncol, 1.0_dp, mos_occ, xsxmat, 0.0_dp, scrv)
541 662 : CALL cp_dbcsr_plus_fm_fm_t(w_matrix, matrix_v=scrv, matrix_g=mos_occ, ncol=ncol, symmetry_mode=1)
542 :
543 662 : CALL cp_fm_release(scrv)
544 662 : CALL cp_fm_release(xsxmat)
545 :
546 662 : CALL timestop(handle)
547 :
548 662 : END SUBROUTINE calculate_xwx_matrix
549 :
550 : END MODULE qs_density_matrices
|
