Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief groups fairly general SCF methods, so that modules other than qs_scf can use them too
10 : !> split off from qs_scf to reduce dependencies
11 : !> \par History
12 : !> - Joost VandeVondele (03.2006)
13 : !> - combine_ks_matrices added (05.04.06,MK)
14 : !> - second ROKS scheme added (15.04.06,MK)
15 : !> - MO occupation management moved (29.08.2008,MK)
16 : !> - correct_mo_eigenvalues was moved from qs_mo_types;
17 : !> new subroutine shift_unocc_mos (03.2016, Sergey Chulkov)
18 : ! **************************************************************************************************
19 : MODULE qs_scf_methods
20 : USE cp_dbcsr_api, ONLY: &
21 : dbcsr_add, dbcsr_desymmetrize, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
22 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
23 : dbcsr_multiply, dbcsr_p_type, dbcsr_type
24 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
25 : cp_dbcsr_sm_fm_multiply
26 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
27 : cp_fm_symm,&
28 : cp_fm_triangular_multiply,&
29 : cp_fm_uplo_to_full
30 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_reduce,&
31 : cp_fm_cholesky_restore
32 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
33 : cp_fm_block_jacobi,&
34 : cp_fm_geeig
35 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
36 : cp_fm_struct_equivalent,&
37 : cp_fm_struct_release,&
38 : cp_fm_struct_type
39 : USE cp_fm_types, ONLY: cp_fm_create,&
40 : cp_fm_get_info,&
41 : cp_fm_release,&
42 : cp_fm_to_fm,&
43 : cp_fm_type
44 : USE input_constants, ONLY: cholesky_inverse,&
45 : cholesky_off,&
46 : cholesky_reduce,&
47 : cholesky_restore
48 : USE kinds, ONLY: dp
49 : USE message_passing, ONLY: mp_para_env_type
50 : USE parallel_gemm_api, ONLY: parallel_gemm
51 : USE qs_density_mixing_types, ONLY: mixing_storage_type
52 : USE qs_mo_types, ONLY: get_mo_set,&
53 : mo_set_type
54 : #include "./base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 :
58 : PRIVATE
59 :
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_methods'
61 : REAL(KIND=dp), PARAMETER :: ratio = 0.25_dp
62 :
63 : PUBLIC :: combine_ks_matrices, &
64 : cp_sm_mix, &
65 : eigensolver, &
66 : eigensolver_generalized, &
67 : eigensolver_dbcsr, &
68 : eigensolver_symm, &
69 : eigensolver_simple, &
70 : scf_env_density_mixing
71 :
72 : INTERFACE combine_ks_matrices
73 : MODULE PROCEDURE combine_ks_matrices_1, &
74 : combine_ks_matrices_2
75 : END INTERFACE combine_ks_matrices
76 :
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief perform (if requested) a density mixing
81 : !> \param p_mix_new New density matrices
82 : !> \param mixing_store ...
83 : !> \param rho_ao Density environment
84 : !> \param para_env ...
85 : !> \param iter_delta ...
86 : !> \param iter_count ...
87 : !> \param diis ...
88 : !> \param invert Invert mixing
89 : !> \par History
90 : !> 02.2003 created [fawzi]
91 : !> 08.2014 adapted for kpoints [JGH]
92 : !> \author fawzi
93 : ! **************************************************************************************************
94 161841 : SUBROUTINE scf_env_density_mixing(p_mix_new, mixing_store, rho_ao, para_env, &
95 : iter_delta, iter_count, diis, invert)
96 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_mix_new
97 : TYPE(mixing_storage_type), POINTER :: mixing_store
98 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
99 : TYPE(mp_para_env_type), POINTER :: para_env
100 : REAL(KIND=dp), INTENT(INOUT) :: iter_delta
101 : INTEGER, INTENT(IN) :: iter_count
102 : LOGICAL, INTENT(in), OPTIONAL :: diis, invert
103 :
104 : CHARACTER(len=*), PARAMETER :: routineN = 'scf_env_density_mixing'
105 :
106 : INTEGER :: handle, ic, ispin
107 : LOGICAL :: my_diis, my_invert
108 : REAL(KIND=dp) :: my_p_mix, tmp
109 :
110 161841 : CALL timeset(routineN, handle)
111 :
112 161841 : my_diis = .FALSE.
113 161841 : IF (PRESENT(diis)) my_diis = diis
114 161841 : my_invert = .FALSE.
115 161841 : IF (PRESENT(invert)) my_invert = invert
116 161841 : my_p_mix = mixing_store%alpha
117 161841 : IF (my_diis .OR. iter_count < mixing_store%nskip_mixing) THEN
118 112235 : my_p_mix = 1.0_dp
119 : END IF
120 :
121 161841 : iter_delta = 0.0_dp
122 161841 : CPASSERT(ASSOCIATED(p_mix_new))
123 1487270 : DO ic = 1, SIZE(p_mix_new, 2)
124 2967299 : DO ispin = 1, SIZE(p_mix_new, 1)
125 2805458 : IF (my_invert) THEN
126 162664 : CPASSERT(my_p_mix /= 0.0_dp)
127 162664 : IF (my_p_mix /= 1.0_dp) THEN
128 : CALL dbcsr_add(matrix_a=p_mix_new(ispin, ic)%matrix, &
129 : alpha_scalar=1.0_dp/my_p_mix, &
130 : matrix_b=rho_ao(ispin, ic)%matrix, &
131 24914 : beta_scalar=(my_p_mix - 1.0_dp)/my_p_mix)
132 : END IF
133 : ELSE
134 : CALL cp_sm_mix(m1=p_mix_new(ispin, ic)%matrix, &
135 : m2=rho_ao(ispin, ic)%matrix, &
136 : p_mix=my_p_mix, &
137 : delta=tmp, &
138 1317365 : para_env=para_env)
139 1317365 : iter_delta = MAX(iter_delta, tmp)
140 : END IF
141 : END DO
142 : END DO
143 :
144 161841 : CALL timestop(handle)
145 :
146 161841 : END SUBROUTINE scf_env_density_mixing
147 :
148 : ! **************************************************************************************************
149 : !> \brief Diagonalise the Kohn-Sham matrix to get a new set of MO eigen-
150 : !> vectors and MO eigenvalues. ks will be modified
151 : !> \param matrix_ks_fm ...
152 : !> \param mo_set ...
153 : !> \param ortho ...
154 : !> \param work ...
155 : !> \param cholesky_method ...
156 : !> \param do_level_shift activate the level shifting technique
157 : !> \param level_shift amount of shift applied (in a.u.)
158 : !> \param matrix_u_fm matrix U : S (overlap matrix) = U^T * U
159 : !> \param use_jacobi ...
160 : !> \date 01.05.2001
161 : !> \author Matthias Krack
162 : !> \version 1.0
163 : ! **************************************************************************************************
164 243562 : SUBROUTINE eigensolver(matrix_ks_fm, mo_set, ortho, work, &
165 : cholesky_method, do_level_shift, &
166 : level_shift, matrix_u_fm, use_jacobi)
167 : TYPE(cp_fm_type), INTENT(IN) :: matrix_ks_fm
168 : TYPE(mo_set_type), INTENT(IN) :: mo_set
169 : TYPE(cp_fm_type), INTENT(IN) :: ortho, work
170 : INTEGER, INTENT(inout) :: cholesky_method
171 : LOGICAL, INTENT(in) :: do_level_shift
172 : REAL(KIND=dp), INTENT(in) :: level_shift
173 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_u_fm
174 : LOGICAL, INTENT(in) :: use_jacobi
175 :
176 : CHARACTER(len=*), PARAMETER :: routineN = 'eigensolver'
177 :
178 : INTEGER :: handle, homo, nao, nmo
179 121781 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
180 : TYPE(cp_fm_type), POINTER :: mo_coeff
181 :
182 121781 : CALL timeset(routineN, handle)
183 :
184 121781 : NULLIFY (mo_coeff)
185 121781 : NULLIFY (mo_eigenvalues)
186 :
187 : ! Diagonalise the Kohn-Sham matrix
188 :
189 : CALL get_mo_set(mo_set=mo_set, &
190 : nao=nao, &
191 : nmo=nmo, &
192 : homo=homo, &
193 : eigenvalues=mo_eigenvalues, &
194 121781 : mo_coeff=mo_coeff)
195 :
196 121837 : SELECT CASE (cholesky_method)
197 : CASE (cholesky_reduce)
198 56 : CALL cp_fm_cholesky_reduce(matrix_ks_fm, ortho)
199 :
200 56 : IF (do_level_shift) THEN
201 : CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
202 34 : level_shift=level_shift, is_triangular=.TRUE., matrix_u_fm=matrix_u_fm)
203 : END IF
204 :
205 56 : CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
206 56 : CALL cp_fm_cholesky_restore(work, nmo, ortho, mo_coeff, "SOLVE")
207 56 : IF (do_level_shift) THEN
208 34 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
209 : END IF
210 :
211 : CASE (cholesky_restore)
212 120001 : CALL cp_fm_uplo_to_full(matrix_ks_fm, work)
213 : CALL cp_fm_cholesky_restore(matrix_ks_fm, nao, ortho, work, &
214 120001 : "SOLVE", pos="RIGHT")
215 : CALL cp_fm_cholesky_restore(work, nao, ortho, matrix_ks_fm, &
216 120001 : "SOLVE", pos="LEFT", transa="T")
217 :
218 120001 : IF (do_level_shift) THEN
219 : CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
220 82 : level_shift=level_shift, is_triangular=.TRUE., matrix_u_fm=matrix_u_fm)
221 : END IF
222 :
223 120001 : CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
224 120001 : CALL cp_fm_cholesky_restore(work, nmo, ortho, mo_coeff, "SOLVE")
225 :
226 120001 : IF (do_level_shift) THEN
227 82 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
228 : END IF
229 :
230 : CASE (cholesky_inverse)
231 1724 : CALL cp_fm_uplo_to_full(matrix_ks_fm, work)
232 :
233 : CALL cp_fm_triangular_multiply(ortho, matrix_ks_fm, side="R", transpose_tr=.FALSE., &
234 1724 : invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nao, alpha=1.0_dp)
235 : CALL cp_fm_triangular_multiply(ortho, matrix_ks_fm, side="L", transpose_tr=.TRUE., &
236 1724 : invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nao, alpha=1.0_dp)
237 :
238 1724 : IF (do_level_shift) THEN
239 : CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
240 34 : level_shift=level_shift, is_triangular=.TRUE., matrix_u_fm=matrix_u_fm)
241 : END IF
242 :
243 1724 : CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
244 : CALL cp_fm_triangular_multiply(ortho, work, side="L", transpose_tr=.FALSE., &
245 1724 : invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nmo, alpha=1.0_dp)
246 1724 : CALL cp_fm_to_fm(work, mo_coeff, nmo, 1, 1)
247 :
248 123505 : IF (do_level_shift) THEN
249 34 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
250 : END IF
251 :
252 : END SELECT
253 :
254 121781 : IF (use_jacobi) THEN
255 0 : CALL cp_fm_to_fm(mo_coeff, ortho)
256 0 : cholesky_method = cholesky_off
257 : END IF
258 :
259 121781 : CALL timestop(handle)
260 :
261 121781 : END SUBROUTINE eigensolver
262 :
263 : ! **************************************************************************************************
264 : !> \brief Solve the generalized eigenvalue problem
265 : !> \param matrix_ks_fm Kohn-Sham matrix in FM format
266 : !> \param matrix_s Overlap matrix (DBCSR)
267 : !> \param mo_set Molecular orbital set
268 : !> \param work Work matrix (used as eigenvector buffer)
269 : ! **************************************************************************************************
270 0 : SUBROUTINE eigensolver_generalized(matrix_ks_fm, matrix_s, mo_set, work)
271 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_ks_fm
272 : TYPE(dbcsr_type), INTENT(IN) :: matrix_s
273 : TYPE(mo_set_type), INTENT(IN) :: mo_set
274 : TYPE(cp_fm_type), INTENT(INOUT) :: work
275 :
276 : CHARACTER(len=*), PARAMETER :: routineN = 'eigensolver_generalized'
277 :
278 : INTEGER :: handle, nmo
279 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
280 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
281 : TYPE(cp_fm_type) :: s_fm
282 : TYPE(cp_fm_type), POINTER :: mo_coeff
283 :
284 0 : CALL timeset(routineN, handle)
285 :
286 0 : NULLIFY (mo_coeff)
287 0 : NULLIFY (mo_eigenvalues)
288 :
289 0 : CALL get_mo_set(mo_set=mo_set, nmo=nmo, eigenvalues=mo_eigenvalues, mo_coeff=mo_coeff)
290 0 : CALL cp_fm_get_info(matrix_ks_fm, matrix_struct=fm_struct)
291 :
292 : ! Convert S matrix from DBCSR to FM (required for cuSOLVERMp)
293 0 : CALL cp_fm_create(s_fm, fm_struct)
294 0 : CALL copy_dbcsr_to_fm(matrix_s, s_fm)
295 :
296 : ! Solve generalized eigenvalue problem - eigenvectors output to work buffer
297 0 : CALL cp_fm_geeig(matrix_ks_fm, s_fm, mo_coeff, mo_eigenvalues, work)
298 :
299 0 : CALL cp_fm_release(s_fm)
300 :
301 0 : CALL timestop(handle)
302 :
303 0 : END SUBROUTINE eigensolver_generalized
304 :
305 : ! **************************************************************************************************
306 : !> \brief ...
307 : !> \param matrix_ks ...
308 : !> \param matrix_ks_fm ...
309 : !> \param mo_set ...
310 : !> \param ortho_dbcsr ...
311 : !> \param ksbuf1 ...
312 : !> \param ksbuf2 ...
313 : ! **************************************************************************************************
314 8504 : SUBROUTINE eigensolver_dbcsr(matrix_ks, matrix_ks_fm, mo_set, ortho_dbcsr, ksbuf1, ksbuf2)
315 : TYPE(dbcsr_type), INTENT(IN) :: matrix_ks
316 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_ks_fm
317 : TYPE(mo_set_type), INTENT(IN) :: mo_set
318 : TYPE(dbcsr_type), INTENT(IN) :: ortho_dbcsr
319 : TYPE(dbcsr_type), INTENT(INOUT) :: ksbuf1, ksbuf2
320 :
321 : CHARACTER(len=*), PARAMETER :: routineN = 'eigensolver_dbcsr'
322 :
323 : INTEGER :: handle, nao, nmo
324 2126 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
325 : TYPE(cp_fm_type) :: all_evecs, nmo_evecs
326 : TYPE(cp_fm_type), POINTER :: mo_coeff
327 :
328 2126 : CALL timeset(routineN, handle)
329 :
330 2126 : NULLIFY (mo_coeff)
331 2126 : NULLIFY (mo_eigenvalues)
332 :
333 : CALL get_mo_set(mo_set=mo_set, &
334 : nao=nao, &
335 : nmo=nmo, &
336 : eigenvalues=mo_eigenvalues, &
337 2126 : mo_coeff=mo_coeff)
338 :
339 : ! Reduce KS matrix
340 2126 : CALL dbcsr_desymmetrize(matrix_ks, ksbuf2)
341 2126 : CALL dbcsr_multiply('N', 'N', 1.0_dp, ksbuf2, ortho_dbcsr, 0.0_dp, ksbuf1)
342 2126 : CALL dbcsr_multiply('T', 'N', 1.0_dp, ortho_dbcsr, ksbuf1, 0.0_dp, ksbuf2)
343 :
344 : ! Solve the eigenvalue problem
345 2126 : CALL copy_dbcsr_to_fm(ksbuf2, matrix_ks_fm)
346 2126 : CALL cp_fm_create(all_evecs, matrix_ks_fm%matrix_struct)
347 2126 : CALL choose_eigv_solver(matrix_ks_fm, all_evecs, mo_eigenvalues)
348 :
349 : ! select first nmo eigenvectors
350 2126 : CALL cp_fm_create(nmo_evecs, mo_coeff%matrix_struct)
351 2126 : CALL cp_fm_to_fm(msource=all_evecs, mtarget=nmo_evecs, ncol=nmo)
352 2126 : CALL cp_fm_release(all_evecs)
353 :
354 : ! Restore the eigenvector of the general eig. problem
355 2126 : CALL cp_dbcsr_sm_fm_multiply(ortho_dbcsr, nmo_evecs, mo_coeff, nmo)
356 :
357 2126 : CALL cp_fm_release(nmo_evecs)
358 2126 : CALL timestop(handle)
359 :
360 2126 : END SUBROUTINE eigensolver_dbcsr
361 :
362 : ! **************************************************************************************************
363 : !> \brief ...
364 : !> \param matrix_ks_fm ...
365 : !> \param mo_set ...
366 : !> \param ortho ...
367 : !> \param work ...
368 : !> \param do_level_shift activate the level shifting technique
369 : !> \param level_shift amount of shift applied (in a.u.)
370 : !> \param matrix_u_fm matrix U : S (overlap matrix) = U^T * U
371 : !> \param use_jacobi ...
372 : !> \param jacobi_threshold ...
373 : !> \param ortho_red ...
374 : !> \param work_red ...
375 : !> \param matrix_ks_fm_red ...
376 : !> \param matrix_u_fm_red ...
377 : ! **************************************************************************************************
378 532 : SUBROUTINE eigensolver_symm(matrix_ks_fm, mo_set, ortho, work, do_level_shift, &
379 : level_shift, matrix_u_fm, use_jacobi, jacobi_threshold, &
380 : ortho_red, work_red, matrix_ks_fm_red, matrix_u_fm_red)
381 : TYPE(cp_fm_type), INTENT(IN) :: matrix_ks_fm
382 : TYPE(mo_set_type), INTENT(IN) :: mo_set
383 : TYPE(cp_fm_type), INTENT(IN) :: ortho, work
384 : LOGICAL, INTENT(IN) :: do_level_shift
385 : REAL(KIND=dp), INTENT(IN) :: level_shift
386 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_u_fm
387 : LOGICAL, INTENT(IN) :: use_jacobi
388 : REAL(KIND=dp), INTENT(IN) :: jacobi_threshold
389 : TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: ortho_red, work_red, matrix_ks_fm_red, &
390 : matrix_u_fm_red
391 :
392 : CHARACTER(len=*), PARAMETER :: routineN = 'eigensolver_symm'
393 :
394 : INTEGER :: handle, homo, nao, nao_red, nelectron, &
395 : nmo
396 532 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
397 532 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
398 : TYPE(cp_fm_type) :: work_red2
399 : TYPE(cp_fm_type), POINTER :: mo_coeff
400 :
401 532 : CALL timeset(routineN, handle)
402 :
403 532 : NULLIFY (mo_coeff)
404 532 : NULLIFY (mo_eigenvalues)
405 :
406 : ! Diagonalise the Kohn-Sham matrix
407 :
408 : CALL get_mo_set(mo_set=mo_set, &
409 : nao=nao, &
410 : nmo=nmo, &
411 : homo=homo, &
412 : nelectron=nelectron, &
413 : eigenvalues=mo_eigenvalues, &
414 532 : mo_coeff=mo_coeff)
415 :
416 532 : IF (use_jacobi) THEN
417 :
418 0 : CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, matrix_ks_fm, mo_coeff, 0.0_dp, work)
419 : CALL parallel_gemm("T", "N", homo, nao - homo, nao, 1.0_dp, work, mo_coeff, &
420 0 : 0.0_dp, matrix_ks_fm, b_first_col=homo + 1, c_first_col=homo + 1)
421 :
422 : ! Block Jacobi (pseudo-diagonalization, only one sweep)
423 : CALL cp_fm_block_jacobi(matrix_ks_fm, mo_coeff, mo_eigenvalues, &
424 0 : jacobi_threshold, homo + 1)
425 :
426 : ELSE ! full S^(-1/2) has been computed
427 532 : IF (PRESENT(work_red) .AND. PRESENT(ortho_red) .AND. PRESENT(matrix_ks_fm_red)) THEN
428 532 : CALL cp_fm_get_info(ortho_red, ncol_global=nao_red)
429 532 : CALL cp_fm_symm("L", "U", nao, nao_red, 1.0_dp, matrix_ks_fm, ortho_red, 0.0_dp, work_red)
430 532 : CALL parallel_gemm("T", "N", nao_red, nao_red, nao, 1.0_dp, ortho_red, work_red, 0.0_dp, matrix_ks_fm_red)
431 :
432 532 : IF (do_level_shift) THEN
433 : CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm_red, mo_coeff=mo_coeff, homo=homo, &
434 86 : level_shift=level_shift, is_triangular=.FALSE., matrix_u_fm=matrix_u_fm_red)
435 : END IF
436 :
437 532 : CALL cp_fm_create(work_red2, matrix_ks_fm_red%matrix_struct)
438 1596 : ALLOCATE (eigenvalues(nao_red))
439 532 : CALL choose_eigv_solver(matrix_ks_fm_red, work_red2, eigenvalues)
440 6244 : mo_eigenvalues(1:MIN(nao_red, nmo)) = eigenvalues(1:MIN(nao_red, nmo))
441 : CALL parallel_gemm("N", "N", nao, nmo, nao_red, 1.0_dp, ortho_red, work_red2, 0.0_dp, &
442 532 : mo_coeff)
443 1596 : CALL cp_fm_release(work_red2)
444 : ELSE
445 0 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, matrix_ks_fm, ortho, 0.0_dp, work)
446 0 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, ortho, work, 0.0_dp, matrix_ks_fm)
447 0 : IF (do_level_shift) THEN
448 : CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
449 0 : level_shift=level_shift, is_triangular=.FALSE., matrix_u_fm=matrix_u_fm)
450 : END IF
451 0 : CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
452 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, &
453 0 : mo_coeff)
454 : END IF
455 :
456 532 : IF (do_level_shift) THEN
457 86 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
458 : END IF
459 :
460 : END IF
461 :
462 532 : CALL timestop(handle)
463 :
464 1064 : END SUBROUTINE eigensolver_symm
465 :
466 : ! **************************************************************************************************
467 :
468 : ! **************************************************************************************************
469 : !> \brief ...
470 : !> \param matrix_ks ...
471 : !> \param mo_set ...
472 : !> \param work ...
473 : !> \param do_level_shift activate the level shifting technique
474 : !> \param level_shift amount of shift applied (in a.u.)
475 : !> \param use_jacobi ...
476 : !> \param jacobi_threshold ...
477 : ! **************************************************************************************************
478 37296 : SUBROUTINE eigensolver_simple(matrix_ks, mo_set, work, do_level_shift, &
479 : level_shift, use_jacobi, jacobi_threshold)
480 :
481 : TYPE(cp_fm_type), INTENT(IN) :: matrix_ks
482 : TYPE(mo_set_type), INTENT(IN) :: mo_set
483 : TYPE(cp_fm_type), INTENT(IN) :: work
484 : LOGICAL, INTENT(IN) :: do_level_shift
485 : REAL(KIND=dp), INTENT(IN) :: level_shift
486 : LOGICAL, INTENT(IN) :: use_jacobi
487 : REAL(KIND=dp), INTENT(IN) :: jacobi_threshold
488 :
489 : CHARACTER(len=*), PARAMETER :: routineN = 'eigensolver_simple'
490 :
491 : INTEGER :: handle, homo, nao, nelectron, nmo
492 18648 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
493 : TYPE(cp_fm_type), POINTER :: mo_coeff
494 :
495 18648 : CALL timeset(routineN, handle)
496 :
497 18648 : NULLIFY (mo_coeff)
498 18648 : NULLIFY (mo_eigenvalues)
499 :
500 : CALL get_mo_set(mo_set=mo_set, &
501 : nao=nao, &
502 : nmo=nmo, &
503 : homo=homo, &
504 : nelectron=nelectron, &
505 : eigenvalues=mo_eigenvalues, &
506 18648 : mo_coeff=mo_coeff)
507 :
508 18648 : IF (do_level_shift) THEN
509 : ! matrix_u_fm is simply an identity matrix, so we omit it here
510 : CALL shift_unocc_mos(matrix_ks_fm=matrix_ks, mo_coeff=mo_coeff, homo=homo, &
511 0 : level_shift=level_shift, is_triangular=.FALSE.)
512 : END IF
513 :
514 18648 : IF (use_jacobi) THEN
515 18 : CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, matrix_ks, mo_coeff, 0.0_dp, work)
516 : CALL parallel_gemm("T", "N", homo, nao - homo, nao, 1.0_dp, work, mo_coeff, &
517 18 : 0.0_dp, matrix_ks, b_first_col=homo + 1, c_first_col=homo + 1)
518 : ! Block Jacobi (pseudo-diagonalization, only one sweep)
519 18 : CALL cp_fm_block_jacobi(matrix_ks, mo_coeff, mo_eigenvalues, jacobi_threshold, homo + 1)
520 : ELSE
521 :
522 18630 : CALL choose_eigv_solver(matrix_ks, work, mo_eigenvalues)
523 :
524 18630 : CALL cp_fm_to_fm(work, mo_coeff, nmo, 1, 1)
525 :
526 : END IF
527 :
528 18648 : IF (do_level_shift) THEN
529 0 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
530 : END IF
531 :
532 18648 : CALL timestop(handle)
533 :
534 18648 : END SUBROUTINE eigensolver_simple
535 :
536 : ! **************************************************************************************************
537 : !> \brief Perform a mixing of the given matrixes into the first matrix
538 : !> m1 = m2 + p_mix (m1-m2)
539 : !> \param m1 first (new) matrix, is modified
540 : !> \param m2 the second (old) matrix
541 : !> \param p_mix how much m1 is conserved (0: none, 1: all)
542 : !> \param delta maximum norm of m1-m2
543 : !> \param para_env ...
544 : !> \param m3 ...
545 : !> \par History
546 : !> 02.2003 rewamped [fawzi]
547 : !> \author fawzi
548 : !> \note
549 : !> if you what to store the result in m2 swap m1 and m2 an use
550 : !> (1-pmix) as pmix
551 : !> para_env should be removed (embedded in matrix)
552 : ! **************************************************************************************************
553 1979509 : SUBROUTINE cp_sm_mix(m1, m2, p_mix, delta, para_env, m3)
554 :
555 : TYPE(dbcsr_type), POINTER :: m1, m2
556 : REAL(KIND=dp), INTENT(IN) :: p_mix
557 : REAL(KIND=dp), INTENT(OUT) :: delta
558 : TYPE(mp_para_env_type), POINTER :: para_env
559 : TYPE(dbcsr_type), OPTIONAL, POINTER :: m3
560 :
561 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_sm_mix'
562 :
563 : INTEGER :: handle, i, iblock_col, iblock_row, j
564 : LOGICAL :: found
565 1979509 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_delta_block, p_new_block, p_old_block
566 : TYPE(dbcsr_iterator_type) :: iter
567 :
568 1979509 : CALL timeset(routineN, handle)
569 1979509 : delta = 0.0_dp
570 :
571 : !$OMP PARALLEL DEFAULT(NONE) SHARED(m1, m2, m3, p_mix) &
572 : !$OMP PRIVATE(iter, iblock_row, iblock_col, p_new_block, p_old_block, p_delta_block, found, i, j) &
573 1979509 : !$OMP REDUCTION(MAX:delta)
574 : CALL dbcsr_iterator_start(iter, m1)
575 : DO WHILE (dbcsr_iterator_blocks_left(iter))
576 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, p_new_block)
577 : CALL dbcsr_get_block_p(matrix=m2, row=iblock_row, col=iblock_col, &
578 : BLOCK=p_old_block, found=found)
579 : CPASSERT(ASSOCIATED(p_old_block))
580 : IF (PRESENT(m3)) THEN
581 : CALL dbcsr_get_block_p(matrix=m3, row=iblock_row, col=iblock_col, &
582 : BLOCK=p_delta_block, found=found)
583 : CPASSERT(ASSOCIATED(p_delta_block))
584 :
585 : DO j = 1, SIZE(p_new_block, 2)
586 : DO i = 1, SIZE(p_new_block, 1)
587 : p_delta_block(i, j) = p_new_block(i, j) - p_old_block(i, j)
588 : delta = MAX(delta, ABS(p_delta_block(i, j)))
589 : END DO
590 : END DO
591 : ELSE
592 : DO j = 1, SIZE(p_new_block, 2)
593 : DO i = 1, SIZE(p_new_block, 1)
594 : p_new_block(i, j) = p_new_block(i, j) - p_old_block(i, j)
595 : delta = MAX(delta, ABS(p_new_block(i, j)))
596 : p_new_block(i, j) = p_old_block(i, j) + p_mix*p_new_block(i, j)
597 : END DO
598 : END DO
599 : END IF
600 : END DO
601 : CALL dbcsr_iterator_stop(iter)
602 : !$OMP END PARALLEL
603 :
604 1979509 : CALL para_env%max(delta)
605 :
606 1979509 : CALL timestop(handle)
607 :
608 1979509 : END SUBROUTINE cp_sm_mix
609 :
610 : ! **************************************************************************************************
611 : !> \brief ...
612 : !> \param ksa ...
613 : !> \param ksb ...
614 : !> \param occa ...
615 : !> \param occb ...
616 : !> \param roks_parameter ...
617 : ! **************************************************************************************************
618 966 : SUBROUTINE combine_ks_matrices_1(ksa, ksb, occa, occb, roks_parameter)
619 :
620 : ! Combine the alpha and beta Kohn-Sham matrices during a restricted open
621 : ! Kohn-Sham (ROKS) calculation
622 : ! On input ksa and ksb contain the alpha and beta Kohn-Sham matrices,
623 : ! respectively. occa and occb contain the corresponding MO occupation
624 : ! numbers. On output the combined ROKS operator matrix is returned in ksa.
625 :
626 : ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
627 : ! - M. F. Guest and V. R. Saunders, Mol. Phys. 28(3), 819 (1974)
628 :
629 : TYPE(cp_fm_type), INTENT(IN) :: ksa, ksb
630 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occa, occb
631 : REAL(KIND=dp), DIMENSION(0:2, 0:2, 1:2), &
632 : INTENT(IN) :: roks_parameter
633 :
634 : CHARACTER(LEN=*), PARAMETER :: routineN = 'combine_ks_matrices_1'
635 :
636 : INTEGER :: handle, i, icol_global, icol_local, &
637 : irow_global, irow_local, j, &
638 : ncol_local, nrow_local
639 966 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
640 : LOGICAL :: compatible_matrices
641 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
642 966 : POINTER :: fa, fb
643 : TYPE(cp_fm_struct_type), POINTER :: ksa_struct, ksb_struct
644 :
645 : ! -------------------------------------------------------------------------
646 :
647 966 : CALL timeset(routineN, handle)
648 :
649 : CALL cp_fm_get_info(matrix=ksa, &
650 : matrix_struct=ksa_struct, &
651 : nrow_local=nrow_local, &
652 : ncol_local=ncol_local, &
653 : row_indices=row_indices, &
654 : col_indices=col_indices, &
655 966 : local_data=fa)
656 :
657 : CALL cp_fm_get_info(matrix=ksb, &
658 : matrix_struct=ksb_struct, &
659 966 : local_data=fb)
660 :
661 966 : compatible_matrices = cp_fm_struct_equivalent(ksa_struct, ksb_struct)
662 966 : CPASSERT(compatible_matrices)
663 :
664 15528 : IF (SUM(occb) == 0.0_dp) fb = 0.0_dp
665 :
666 15528 : DO icol_local = 1, ncol_local
667 14562 : icol_global = col_indices(icol_local)
668 14562 : j = INT(occa(icol_global)) + INT(occb(icol_global))
669 174995 : DO irow_local = 1, nrow_local
670 159467 : irow_global = row_indices(irow_local)
671 159467 : i = INT(occa(irow_global)) + INT(occb(irow_global))
672 : fa(irow_local, icol_local) = &
673 : roks_parameter(i, j, 1)*fa(irow_local, icol_local) + &
674 174029 : roks_parameter(i, j, 2)*fb(irow_local, icol_local)
675 : END DO
676 : END DO
677 :
678 966 : CALL timestop(handle)
679 :
680 966 : END SUBROUTINE combine_ks_matrices_1
681 :
682 : ! **************************************************************************************************
683 : !> \brief ...
684 : !> \param ksa ...
685 : !> \param ksb ...
686 : !> \param occa ...
687 : !> \param occb ...
688 : !> \param f ...
689 : !> \param nalpha ...
690 : !> \param nbeta ...
691 : ! **************************************************************************************************
692 0 : SUBROUTINE combine_ks_matrices_2(ksa, ksb, occa, occb, f, nalpha, nbeta)
693 :
694 : ! Combine the alpha and beta Kohn-Sham matrices during a restricted open
695 : ! Kohn-Sham (ROKS) calculation
696 : ! On input ksa and ksb contain the alpha and beta Kohn-Sham matrices,
697 : ! respectively. occa and occb contain the corresponding MO occupation
698 : ! numbers. On output the combined ROKS operator matrix is returned in ksa.
699 :
700 : ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
701 : ! - M. Filatov and S. Shaik, Chem. Phys. Lett. 288, 689 (1998)
702 :
703 : TYPE(cp_fm_type), INTENT(IN) :: ksa, ksb
704 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occa, occb
705 : REAL(KIND=dp), INTENT(IN) :: f
706 : INTEGER, INTENT(IN) :: nalpha, nbeta
707 :
708 : CHARACTER(LEN=*), PARAMETER :: routineN = 'combine_ks_matrices_2'
709 :
710 : INTEGER :: handle, icol_global, icol_local, &
711 : irow_global, irow_local, ncol_local, &
712 : nrow_local
713 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
714 : LOGICAL :: compatible_matrices
715 : REAL(KIND=dp) :: beta, t1, t2, ta, tb
716 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
717 0 : POINTER :: fa, fb
718 : TYPE(cp_fm_struct_type), POINTER :: ksa_struct, ksb_struct
719 :
720 : ! -------------------------------------------------------------------------
721 :
722 0 : CALL timeset(routineN, handle)
723 :
724 : CALL cp_fm_get_info(matrix=ksa, &
725 : matrix_struct=ksa_struct, &
726 : nrow_local=nrow_local, &
727 : ncol_local=ncol_local, &
728 : row_indices=row_indices, &
729 : col_indices=col_indices, &
730 0 : local_data=fa)
731 :
732 : CALL cp_fm_get_info(matrix=ksb, &
733 : matrix_struct=ksb_struct, &
734 0 : local_data=fb)
735 :
736 0 : compatible_matrices = cp_fm_struct_equivalent(ksa_struct, ksb_struct)
737 0 : CPASSERT(compatible_matrices)
738 :
739 0 : beta = 1.0_dp/(1.0_dp - f)
740 :
741 0 : DO icol_local = 1, ncol_local
742 :
743 0 : icol_global = col_indices(icol_local)
744 :
745 0 : DO irow_local = 1, nrow_local
746 :
747 0 : irow_global = row_indices(irow_local)
748 :
749 0 : t1 = 0.5_dp*(fa(irow_local, icol_local) + fb(irow_local, icol_local))
750 :
751 0 : IF ((0 < irow_global) .AND. (irow_global <= nbeta)) THEN
752 0 : IF ((0 < icol_global) .AND. (icol_global <= nbeta)) THEN
753 : ! closed-closed
754 0 : fa(irow_local, icol_local) = t1
755 0 : ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
756 : ! closed-open
757 0 : ta = 0.5_dp*(f - REAL(occa(icol_global), KIND=dp))/f
758 0 : tb = 0.5_dp*(f - REAL(occb(icol_global), KIND=dp))/f
759 0 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
760 0 : fa(irow_local, icol_local) = t1 + (beta - 1.0_dp)*t2
761 : ELSE
762 : ! closed-virtual
763 0 : fa(irow_local, icol_local) = t1
764 : END IF
765 0 : ELSE IF ((nbeta < irow_global) .AND. (irow_global <= nalpha)) THEN
766 : IF ((0 < irow_global) .AND. (irow_global <= nbeta)) THEN
767 : ! open-closed
768 : ta = 0.5_dp*(f - REAL(occa(irow_global), KIND=dp))/f
769 : tb = 0.5_dp*(f - REAL(occb(irow_global), KIND=dp))/f
770 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
771 : fa(irow_local, icol_local) = t1 + (beta - 1.0_dp)*t2
772 0 : ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
773 : ! open-open
774 0 : ta = 0.5_dp*(f - REAL(occa(icol_global), KIND=dp))/f
775 0 : tb = 0.5_dp*(f - REAL(occb(icol_global), KIND=dp))/f
776 0 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
777 0 : IF (irow_global == icol_global) THEN
778 0 : fa(irow_local, icol_local) = t1 - t2
779 : ELSE
780 0 : fa(irow_local, icol_local) = t1 - 0.5_dp*t2
781 : END IF
782 : ELSE
783 : ! open-virtual
784 0 : ta = 0.5_dp*(f - REAL(occa(irow_global), KIND=dp))/f
785 0 : tb = 0.5_dp*(f - REAL(occb(irow_global), KIND=dp))/f
786 0 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
787 0 : fa(irow_local, icol_local) = t1 - t2
788 : END IF
789 : ELSE
790 0 : IF ((0 < irow_global) .AND. (irow_global < nbeta)) THEN
791 : ! virtual-closed
792 0 : fa(irow_local, icol_local) = t1
793 0 : ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
794 : ! virtual-open
795 0 : ta = 0.5_dp*(f - REAL(occa(icol_global), KIND=dp))/f
796 0 : tb = 0.5_dp*(f - REAL(occb(icol_global), KIND=dp))/f
797 0 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
798 0 : fa(irow_local, icol_local) = t1 - t2
799 : ELSE
800 : ! virtual-virtual
801 0 : fa(irow_local, icol_local) = t1
802 : END IF
803 : END IF
804 :
805 : END DO
806 : END DO
807 :
808 0 : CALL timestop(handle)
809 :
810 0 : END SUBROUTINE combine_ks_matrices_2
811 :
812 : ! **************************************************************************************************
813 : !> \brief Correct MO eigenvalues after MO level shifting.
814 : !> \param mo_eigenvalues vector of eigenvalues
815 : !> \param homo index of the highest occupied molecular orbital
816 : !> \param nmo number of molecular orbitals
817 : !> \param level_shift amount of applied level shifting (in a.u.)
818 : !> \date 19.04.2002
819 : !> \par History
820 : !> - correct_mo_eigenvalues added (18.04.02,MK)
821 : !> - moved from module qs_mo_types, revised interface (03.2016, Sergey Chulkov)
822 : !> \author MK
823 : !> \version 1.0
824 : ! **************************************************************************************************
825 236 : PURE SUBROUTINE correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
826 :
827 : REAL(kind=dp), DIMENSION(:), INTENT(inout) :: mo_eigenvalues
828 : INTEGER, INTENT(in) :: homo, nmo
829 : REAL(kind=dp), INTENT(in) :: level_shift
830 :
831 : INTEGER :: imo
832 :
833 6182 : DO imo = homo + 1, nmo
834 6182 : mo_eigenvalues(imo) = mo_eigenvalues(imo) - level_shift
835 : END DO
836 :
837 236 : END SUBROUTINE correct_mo_eigenvalues
838 :
839 : ! **************************************************************************************************
840 : !> \brief Adjust the Kohn-Sham matrix by shifting the orbital energies of all
841 : !> unoccupied molecular orbitals
842 : !> \param matrix_ks_fm transformed Kohn-Sham matrix = U^{-1,T} * KS * U^{-1}
843 : !> \param mo_coeff matrix of molecular orbitals (C)
844 : !> \param homo number of occupied molecular orbitals
845 : !> \param level_shift amount of shift applying (in a.u.)
846 : !> \param is_triangular indicates that matrix_u_fm contains an upper triangular matrix
847 : !> \param matrix_u_fm matrix U: S (overlap matrix) = U^T * U;
848 : !> assume an identity matrix if omitted
849 : !> \par History
850 : !> 03.2016 created [Sergey Chulkov]
851 : ! **************************************************************************************************
852 236 : SUBROUTINE shift_unocc_mos(matrix_ks_fm, mo_coeff, homo, &
853 : level_shift, is_triangular, matrix_u_fm)
854 :
855 : TYPE(cp_fm_type), INTENT(IN) :: matrix_ks_fm, mo_coeff
856 : INTEGER, INTENT(in) :: homo
857 : REAL(kind=dp), INTENT(in) :: level_shift
858 : LOGICAL, INTENT(in) :: is_triangular
859 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_u_fm
860 :
861 : CHARACTER(len=*), PARAMETER :: routineN = 'shift_unocc_mos'
862 :
863 : INTEGER :: handle, nao, nao_red, nmo
864 236 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights
865 : TYPE(cp_fm_struct_type), POINTER :: ao_mo_fmstruct
866 : TYPE(cp_fm_type) :: u_mo, u_mo_scaled
867 :
868 236 : CALL timeset(routineN, handle)
869 :
870 236 : IF (PRESENT(matrix_u_fm)) THEN
871 236 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
872 236 : CALL cp_fm_get_info(matrix_u_fm, nrow_global=nao_red, ncol_global=nao)
873 : ELSE
874 0 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo, nrow_global=nao)
875 0 : nao_red = nao
876 : END IF
877 :
878 236 : NULLIFY (ao_mo_fmstruct)
879 : CALL cp_fm_struct_create(ao_mo_fmstruct, nrow_global=nao_red, ncol_global=nmo, &
880 236 : para_env=mo_coeff%matrix_struct%para_env, context=mo_coeff%matrix_struct%context)
881 :
882 236 : CALL cp_fm_create(u_mo, ao_mo_fmstruct)
883 236 : CALL cp_fm_create(u_mo_scaled, ao_mo_fmstruct)
884 :
885 236 : CALL cp_fm_struct_release(ao_mo_fmstruct)
886 :
887 : ! U * C
888 236 : IF (PRESENT(matrix_u_fm)) THEN
889 236 : IF (is_triangular) THEN
890 150 : CALL cp_fm_to_fm(mo_coeff, u_mo)
891 : CALL cp_fm_triangular_multiply(matrix_u_fm, u_mo, side="L", transpose_tr=.FALSE., &
892 150 : invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nmo, alpha=1.0_dp)
893 : ELSE
894 86 : CALL parallel_gemm("N", "N", nao_red, nmo, nao, 1.0_dp, matrix_u_fm, mo_coeff, 0.0_dp, u_mo)
895 : END IF
896 : ELSE
897 : ! assume U is an identity matrix
898 0 : CALL cp_fm_to_fm(mo_coeff, u_mo)
899 : END IF
900 :
901 236 : CALL cp_fm_to_fm(u_mo, u_mo_scaled)
902 :
903 : ! up-shift all unoccupied molecular orbitals by the amount of 'level_shift'
904 : ! weight = diag(DELTA) = (0, ... 0, level_shift, ..., level_shift)
905 : ! MO index : 1 .. homo homo+1 ... nmo
906 708 : ALLOCATE (weights(nmo))
907 1248 : weights(1:homo) = 0.0_dp
908 6472 : weights(homo + 1:nmo) = level_shift
909 : ! DELTA * U * C
910 : ! DELTA is a diagonal matrix, so simply scale all the columns of (U * C) by weights(:)
911 236 : CALL cp_fm_column_scale(u_mo_scaled, weights)
912 236 : DEALLOCATE (weights)
913 :
914 : ! NewKS = U^{-1,T} * KS * U^{-1} + (U * C) * DELTA * (U * C)^T
915 236 : CALL parallel_gemm("N", "T", nao_red, nao_red, nmo, 1.0_dp, u_mo, u_mo_scaled, 1.0_dp, matrix_ks_fm)
916 :
917 236 : CALL cp_fm_release(u_mo_scaled)
918 236 : CALL cp_fm_release(u_mo)
919 :
920 236 : CALL timestop(handle)
921 :
922 472 : END SUBROUTINE shift_unocc_mos
923 :
924 : END MODULE qs_scf_methods
|
