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_cusolver_api, ONLY: cp_fm_general_cusolver
33 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
34 : cp_fm_block_jacobi
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 157557 : 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 157557 : CALL timeset(routineN, handle)
111 :
112 157557 : my_diis = .FALSE.
113 157557 : IF (PRESENT(diis)) my_diis = diis
114 157557 : my_invert = .FALSE.
115 157557 : IF (PRESENT(invert)) my_invert = invert
116 157557 : my_p_mix = mixing_store%alpha
117 157557 : IF (my_diis .OR. iter_count < mixing_store%nskip_mixing) THEN
118 110115 : my_p_mix = 1.0_dp
119 : END IF
120 :
121 157557 : iter_delta = 0.0_dp
122 157557 : CPASSERT(ASSOCIATED(p_mix_new))
123 1450150 : DO ic = 1, SIZE(p_mix_new, 2)
124 2893271 : DO ispin = 1, SIZE(p_mix_new, 1)
125 2735714 : IF (my_invert) THEN
126 159974 : CPASSERT(my_p_mix /= 0.0_dp)
127 159974 : 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 22958 : 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 1283147 : para_env=para_env)
139 1283147 : iter_delta = MAX(iter_delta, tmp)
140 : END IF
141 : END DO
142 : END DO
143 :
144 157557 : CALL timestop(handle)
145 :
146 157557 : 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 226926 : 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 113463 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
180 : TYPE(cp_fm_type), POINTER :: mo_coeff
181 :
182 113463 : CALL timeset(routineN, handle)
183 :
184 113463 : NULLIFY (mo_coeff)
185 113463 : 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 113463 : mo_coeff=mo_coeff)
195 :
196 113511 : SELECT CASE (cholesky_method)
197 : CASE (cholesky_reduce)
198 48 : CALL cp_fm_cholesky_reduce(matrix_ks_fm, ortho)
199 :
200 48 : IF (do_level_shift) THEN
201 : CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
202 28 : level_shift=level_shift, is_triangular=.TRUE., matrix_u_fm=matrix_u_fm)
203 : END IF
204 :
205 48 : CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
206 48 : CALL cp_fm_cholesky_restore(work, nmo, ortho, mo_coeff, "SOLVE")
207 48 : IF (do_level_shift) THEN
208 28 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
209 : END IF
210 :
211 : CASE (cholesky_restore)
212 112199 : CALL cp_fm_uplo_to_full(matrix_ks_fm, work)
213 : CALL cp_fm_cholesky_restore(matrix_ks_fm, nao, ortho, work, &
214 112199 : "SOLVE", pos="RIGHT")
215 : CALL cp_fm_cholesky_restore(work, nao, ortho, matrix_ks_fm, &
216 112199 : "SOLVE", pos="LEFT", transa="T")
217 :
218 112199 : IF (do_level_shift) THEN
219 : CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
220 88 : level_shift=level_shift, is_triangular=.TRUE., matrix_u_fm=matrix_u_fm)
221 : END IF
222 :
223 112199 : CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
224 112199 : CALL cp_fm_cholesky_restore(work, nmo, ortho, mo_coeff, "SOLVE")
225 :
226 112199 : IF (do_level_shift) THEN
227 88 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
228 : END IF
229 :
230 : CASE (cholesky_inverse)
231 1216 : 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 1216 : 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 1216 : invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nao, alpha=1.0_dp)
237 :
238 1216 : IF (do_level_shift) THEN
239 : CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
240 28 : level_shift=level_shift, is_triangular=.TRUE., matrix_u_fm=matrix_u_fm)
241 : END IF
242 :
243 1216 : CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
244 : CALL cp_fm_triangular_multiply(ortho, work, side="L", transpose_tr=.FALSE., &
245 1216 : invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nmo, alpha=1.0_dp)
246 1216 : CALL cp_fm_to_fm(work, mo_coeff, nmo, 1, 1)
247 :
248 114679 : IF (do_level_shift) THEN
249 28 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
250 : END IF
251 :
252 : END SELECT
253 :
254 113463 : 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 113463 : CALL timestop(handle)
260 :
261 113463 : END SUBROUTINE eigensolver
262 :
263 : ! **************************************************************************************************
264 : !> \brief Solve the generalized eigenvalue problem using cusolverMpSygvd
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_general_cusolver(matrix_ks_fm, s_fm, work, mo_eigenvalues)
298 :
299 : ! Copy only the occupied MOs to mo_coeff
300 0 : CALL cp_fm_to_fm(work, mo_coeff, nmo)
301 :
302 0 : CALL cp_fm_release(s_fm)
303 :
304 0 : CALL timestop(handle)
305 :
306 0 : END SUBROUTINE eigensolver_generalized
307 :
308 : ! **************************************************************************************************
309 : !> \brief ...
310 : !> \param matrix_ks ...
311 : !> \param matrix_ks_fm ...
312 : !> \param mo_set ...
313 : !> \param ortho_dbcsr ...
314 : !> \param ksbuf1 ...
315 : !> \param ksbuf2 ...
316 : ! **************************************************************************************************
317 8504 : SUBROUTINE eigensolver_dbcsr(matrix_ks, matrix_ks_fm, mo_set, ortho_dbcsr, ksbuf1, ksbuf2)
318 : TYPE(dbcsr_type), INTENT(IN) :: matrix_ks
319 : TYPE(cp_fm_type), INTENT(INOUT) :: matrix_ks_fm
320 : TYPE(mo_set_type), INTENT(IN) :: mo_set
321 : TYPE(dbcsr_type), INTENT(IN) :: ortho_dbcsr
322 : TYPE(dbcsr_type), INTENT(INOUT) :: ksbuf1, ksbuf2
323 :
324 : CHARACTER(len=*), PARAMETER :: routineN = 'eigensolver_dbcsr'
325 :
326 : INTEGER :: handle, nao, nmo
327 2126 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
328 : TYPE(cp_fm_type) :: all_evecs, nmo_evecs
329 : TYPE(cp_fm_type), POINTER :: mo_coeff
330 :
331 2126 : CALL timeset(routineN, handle)
332 :
333 2126 : NULLIFY (mo_coeff)
334 2126 : NULLIFY (mo_eigenvalues)
335 :
336 : CALL get_mo_set(mo_set=mo_set, &
337 : nao=nao, &
338 : nmo=nmo, &
339 : eigenvalues=mo_eigenvalues, &
340 2126 : mo_coeff=mo_coeff)
341 :
342 : ! Reduce KS matrix
343 2126 : CALL dbcsr_desymmetrize(matrix_ks, ksbuf2)
344 2126 : CALL dbcsr_multiply('N', 'N', 1.0_dp, ksbuf2, ortho_dbcsr, 0.0_dp, ksbuf1)
345 2126 : CALL dbcsr_multiply('T', 'N', 1.0_dp, ortho_dbcsr, ksbuf1, 0.0_dp, ksbuf2)
346 :
347 : ! Solve the eigenvalue problem
348 2126 : CALL copy_dbcsr_to_fm(ksbuf2, matrix_ks_fm)
349 2126 : CALL cp_fm_create(all_evecs, matrix_ks_fm%matrix_struct)
350 2126 : CALL choose_eigv_solver(matrix_ks_fm, all_evecs, mo_eigenvalues)
351 :
352 : ! select first nmo eigenvectors
353 2126 : CALL cp_fm_create(nmo_evecs, mo_coeff%matrix_struct)
354 2126 : CALL cp_fm_to_fm(msource=all_evecs, mtarget=nmo_evecs, ncol=nmo)
355 2126 : CALL cp_fm_release(all_evecs)
356 :
357 : ! Restore the eigenvector of the general eig. problem
358 2126 : CALL cp_dbcsr_sm_fm_multiply(ortho_dbcsr, nmo_evecs, mo_coeff, nmo)
359 :
360 2126 : CALL cp_fm_release(nmo_evecs)
361 2126 : CALL timestop(handle)
362 :
363 2126 : END SUBROUTINE eigensolver_dbcsr
364 :
365 : ! **************************************************************************************************
366 : !> \brief ...
367 : !> \param matrix_ks_fm ...
368 : !> \param mo_set ...
369 : !> \param ortho ...
370 : !> \param work ...
371 : !> \param do_level_shift activate the level shifting technique
372 : !> \param level_shift amount of shift applied (in a.u.)
373 : !> \param matrix_u_fm matrix U : S (overlap matrix) = U^T * U
374 : !> \param use_jacobi ...
375 : !> \param jacobi_threshold ...
376 : !> \param ortho_red ...
377 : !> \param work_red ...
378 : !> \param matrix_ks_fm_red ...
379 : !> \param matrix_u_fm_red ...
380 : ! **************************************************************************************************
381 532 : SUBROUTINE eigensolver_symm(matrix_ks_fm, mo_set, ortho, work, do_level_shift, &
382 : level_shift, matrix_u_fm, use_jacobi, jacobi_threshold, &
383 : ortho_red, work_red, matrix_ks_fm_red, matrix_u_fm_red)
384 : TYPE(cp_fm_type), INTENT(IN) :: matrix_ks_fm
385 : TYPE(mo_set_type), INTENT(IN) :: mo_set
386 : TYPE(cp_fm_type), INTENT(IN) :: ortho, work
387 : LOGICAL, INTENT(IN) :: do_level_shift
388 : REAL(KIND=dp), INTENT(IN) :: level_shift
389 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_u_fm
390 : LOGICAL, INTENT(IN) :: use_jacobi
391 : REAL(KIND=dp), INTENT(IN) :: jacobi_threshold
392 : TYPE(cp_fm_type), INTENT(INOUT), OPTIONAL :: ortho_red, work_red, matrix_ks_fm_red, &
393 : matrix_u_fm_red
394 :
395 : CHARACTER(len=*), PARAMETER :: routineN = 'eigensolver_symm'
396 :
397 : INTEGER :: handle, homo, nao, nao_red, nelectron, &
398 : nmo
399 532 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
400 532 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
401 : TYPE(cp_fm_type) :: work_red2
402 : TYPE(cp_fm_type), POINTER :: mo_coeff
403 :
404 532 : CALL timeset(routineN, handle)
405 :
406 532 : NULLIFY (mo_coeff)
407 532 : NULLIFY (mo_eigenvalues)
408 :
409 : ! Diagonalise the Kohn-Sham matrix
410 :
411 : CALL get_mo_set(mo_set=mo_set, &
412 : nao=nao, &
413 : nmo=nmo, &
414 : homo=homo, &
415 : nelectron=nelectron, &
416 : eigenvalues=mo_eigenvalues, &
417 532 : mo_coeff=mo_coeff)
418 :
419 532 : IF (use_jacobi) THEN
420 :
421 0 : CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, matrix_ks_fm, mo_coeff, 0.0_dp, work)
422 : CALL parallel_gemm("T", "N", homo, nao - homo, nao, 1.0_dp, work, mo_coeff, &
423 0 : 0.0_dp, matrix_ks_fm, b_first_col=homo + 1, c_first_col=homo + 1)
424 :
425 : ! Block Jacobi (pseudo-diagonalization, only one sweep)
426 : CALL cp_fm_block_jacobi(matrix_ks_fm, mo_coeff, mo_eigenvalues, &
427 0 : jacobi_threshold, homo + 1)
428 :
429 : ELSE ! full S^(-1/2) has been computed
430 532 : IF (PRESENT(work_red) .AND. PRESENT(ortho_red) .AND. PRESENT(matrix_ks_fm_red)) THEN
431 532 : CALL cp_fm_get_info(ortho_red, ncol_global=nao_red)
432 532 : CALL cp_fm_symm("L", "U", nao, nao_red, 1.0_dp, matrix_ks_fm, ortho_red, 0.0_dp, work_red)
433 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)
434 :
435 532 : IF (do_level_shift) THEN
436 : CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm_red, mo_coeff=mo_coeff, homo=homo, &
437 86 : level_shift=level_shift, is_triangular=.FALSE., matrix_u_fm=matrix_u_fm_red)
438 : END IF
439 :
440 532 : CALL cp_fm_create(work_red2, matrix_ks_fm_red%matrix_struct)
441 1596 : ALLOCATE (eigenvalues(nao_red))
442 532 : CALL choose_eigv_solver(matrix_ks_fm_red, work_red2, eigenvalues)
443 7108 : mo_eigenvalues(1:MIN(nao_red, nmo)) = eigenvalues(1:MIN(nao_red, nmo))
444 : CALL parallel_gemm("N", "N", nao, nmo, nao_red, 1.0_dp, ortho_red, work_red2, 0.0_dp, &
445 532 : mo_coeff)
446 1596 : CALL cp_fm_release(work_red2)
447 : ELSE
448 0 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, matrix_ks_fm, ortho, 0.0_dp, work)
449 0 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, ortho, work, 0.0_dp, matrix_ks_fm)
450 0 : IF (do_level_shift) THEN
451 : CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, &
452 0 : level_shift=level_shift, is_triangular=.FALSE., matrix_u_fm=matrix_u_fm)
453 : END IF
454 0 : CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
455 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, &
456 0 : mo_coeff)
457 : END IF
458 :
459 532 : IF (do_level_shift) THEN
460 86 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
461 : END IF
462 :
463 : END IF
464 :
465 532 : CALL timestop(handle)
466 :
467 1064 : END SUBROUTINE eigensolver_symm
468 :
469 : ! **************************************************************************************************
470 :
471 : ! **************************************************************************************************
472 : !> \brief ...
473 : !> \param matrix_ks ...
474 : !> \param mo_set ...
475 : !> \param work ...
476 : !> \param do_level_shift activate the level shifting technique
477 : !> \param level_shift amount of shift applied (in a.u.)
478 : !> \param use_jacobi ...
479 : !> \param jacobi_threshold ...
480 : ! **************************************************************************************************
481 37296 : SUBROUTINE eigensolver_simple(matrix_ks, mo_set, work, do_level_shift, &
482 : level_shift, use_jacobi, jacobi_threshold)
483 :
484 : TYPE(cp_fm_type), INTENT(IN) :: matrix_ks
485 : TYPE(mo_set_type), INTENT(IN) :: mo_set
486 : TYPE(cp_fm_type), INTENT(IN) :: work
487 : LOGICAL, INTENT(IN) :: do_level_shift
488 : REAL(KIND=dp), INTENT(IN) :: level_shift
489 : LOGICAL, INTENT(IN) :: use_jacobi
490 : REAL(KIND=dp), INTENT(IN) :: jacobi_threshold
491 :
492 : CHARACTER(len=*), PARAMETER :: routineN = 'eigensolver_simple'
493 :
494 : INTEGER :: handle, homo, nao, nelectron, nmo
495 18648 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
496 : TYPE(cp_fm_type), POINTER :: mo_coeff
497 :
498 18648 : CALL timeset(routineN, handle)
499 :
500 18648 : NULLIFY (mo_coeff)
501 18648 : NULLIFY (mo_eigenvalues)
502 :
503 : CALL get_mo_set(mo_set=mo_set, &
504 : nao=nao, &
505 : nmo=nmo, &
506 : homo=homo, &
507 : nelectron=nelectron, &
508 : eigenvalues=mo_eigenvalues, &
509 18648 : mo_coeff=mo_coeff)
510 :
511 18648 : IF (do_level_shift) THEN
512 : ! matrix_u_fm is simply an identity matrix, so we omit it here
513 : CALL shift_unocc_mos(matrix_ks_fm=matrix_ks, mo_coeff=mo_coeff, homo=homo, &
514 0 : level_shift=level_shift, is_triangular=.FALSE.)
515 : END IF
516 :
517 18648 : IF (use_jacobi) THEN
518 18 : CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, matrix_ks, mo_coeff, 0.0_dp, work)
519 : CALL parallel_gemm("T", "N", homo, nao - homo, nao, 1.0_dp, work, mo_coeff, &
520 18 : 0.0_dp, matrix_ks, b_first_col=homo + 1, c_first_col=homo + 1)
521 : ! Block Jacobi (pseudo-diagonalization, only one sweep)
522 18 : CALL cp_fm_block_jacobi(matrix_ks, mo_coeff, mo_eigenvalues, jacobi_threshold, homo + 1)
523 : ELSE
524 :
525 18630 : CALL choose_eigv_solver(matrix_ks, work, mo_eigenvalues)
526 :
527 18630 : CALL cp_fm_to_fm(work, mo_coeff, nmo, 1, 1)
528 :
529 : END IF
530 :
531 18648 : IF (do_level_shift) THEN
532 0 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
533 : END IF
534 :
535 18648 : CALL timestop(handle)
536 :
537 18648 : END SUBROUTINE eigensolver_simple
538 :
539 : ! **************************************************************************************************
540 : !> \brief Perform a mixing of the given matrixes into the first matrix
541 : !> m1 = m2 + p_mix (m1-m2)
542 : !> \param m1 first (new) matrix, is modified
543 : !> \param m2 the second (old) matrix
544 : !> \param p_mix how much m1 is conserved (0: none, 1: all)
545 : !> \param delta maximum norm of m1-m2
546 : !> \param para_env ...
547 : !> \param m3 ...
548 : !> \par History
549 : !> 02.2003 rewamped [fawzi]
550 : !> \author fawzi
551 : !> \note
552 : !> if you what to store the result in m2 swap m1 and m2 an use
553 : !> (1-pmix) as pmix
554 : !> para_env should be removed (embedded in matrix)
555 : ! **************************************************************************************************
556 3338326 : SUBROUTINE cp_sm_mix(m1, m2, p_mix, delta, para_env, m3)
557 :
558 : TYPE(dbcsr_type), POINTER :: m1, m2
559 : REAL(KIND=dp), INTENT(IN) :: p_mix
560 : REAL(KIND=dp), INTENT(OUT) :: delta
561 : TYPE(mp_para_env_type), POINTER :: para_env
562 : TYPE(dbcsr_type), OPTIONAL, POINTER :: m3
563 :
564 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_sm_mix'
565 :
566 : INTEGER :: handle, i, iblock_col, iblock_row, j
567 : LOGICAL :: found
568 1669163 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_delta_block, p_new_block, p_old_block
569 : TYPE(dbcsr_iterator_type) :: iter
570 :
571 1669163 : CALL timeset(routineN, handle)
572 1669163 : delta = 0.0_dp
573 :
574 1669163 : CALL dbcsr_iterator_start(iter, m1)
575 24965237 : DO WHILE (dbcsr_iterator_blocks_left(iter))
576 23296074 : 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 23296074 : BLOCK=p_old_block, found=found)
579 23296074 : CPASSERT(ASSOCIATED(p_old_block))
580 24965237 : IF (PRESENT(m3)) THEN
581 : CALL dbcsr_get_block_p(matrix=m3, row=iblock_row, col=iblock_col, &
582 6180132 : BLOCK=p_delta_block, found=found)
583 6180132 : CPASSERT(ASSOCIATED(p_delta_block))
584 :
585 39167683 : DO j = 1, SIZE(p_new_block, 2)
586 291813408 : DO i = 1, SIZE(p_new_block, 1)
587 252645725 : p_delta_block(i, j) = p_new_block(i, j) - p_old_block(i, j)
588 285633276 : delta = MAX(delta, ABS(p_delta_block(i, j)))
589 : END DO
590 : END DO
591 : ELSE
592 101745873 : DO j = 1, SIZE(p_new_block, 2)
593 765373580 : DO i = 1, SIZE(p_new_block, 1)
594 663627707 : p_new_block(i, j) = p_new_block(i, j) - p_old_block(i, j)
595 663627707 : delta = MAX(delta, ABS(p_new_block(i, j)))
596 748257638 : 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 1669163 : CALL dbcsr_iterator_stop(iter)
602 :
603 1669163 : CALL para_env%max(delta)
604 :
605 1669163 : CALL timestop(handle)
606 :
607 1669163 : END SUBROUTINE cp_sm_mix
608 :
609 : ! **************************************************************************************************
610 : !> \brief ...
611 : !> \param ksa ...
612 : !> \param ksb ...
613 : !> \param occa ...
614 : !> \param occb ...
615 : !> \param roks_parameter ...
616 : ! **************************************************************************************************
617 1024 : SUBROUTINE combine_ks_matrices_1(ksa, ksb, occa, occb, roks_parameter)
618 :
619 : ! Combine the alpha and beta Kohn-Sham matrices during a restricted open
620 : ! Kohn-Sham (ROKS) calculation
621 : ! On input ksa and ksb contain the alpha and beta Kohn-Sham matrices,
622 : ! respectively. occa and occb contain the corresponding MO occupation
623 : ! numbers. On output the combined ROKS operator matrix is returned in ksa.
624 :
625 : ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
626 : ! - M. F. Guest and V. R. Saunders, Mol. Phys. 28(3), 819 (1974)
627 :
628 : TYPE(cp_fm_type), INTENT(IN) :: ksa, ksb
629 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occa, occb
630 : REAL(KIND=dp), DIMENSION(0:2, 0:2, 1:2), &
631 : INTENT(IN) :: roks_parameter
632 :
633 : CHARACTER(LEN=*), PARAMETER :: routineN = 'combine_ks_matrices_1'
634 :
635 : INTEGER :: handle, i, icol_global, icol_local, &
636 : irow_global, irow_local, j, &
637 : ncol_local, nrow_local
638 1024 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
639 : LOGICAL :: compatible_matrices
640 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
641 1024 : POINTER :: fa, fb
642 : TYPE(cp_fm_struct_type), POINTER :: ksa_struct, ksb_struct
643 :
644 : ! -------------------------------------------------------------------------
645 :
646 1024 : CALL timeset(routineN, handle)
647 :
648 : CALL cp_fm_get_info(matrix=ksa, &
649 : matrix_struct=ksa_struct, &
650 : nrow_local=nrow_local, &
651 : ncol_local=ncol_local, &
652 : row_indices=row_indices, &
653 : col_indices=col_indices, &
654 1024 : local_data=fa)
655 :
656 : CALL cp_fm_get_info(matrix=ksb, &
657 : matrix_struct=ksb_struct, &
658 1024 : local_data=fb)
659 :
660 1024 : compatible_matrices = cp_fm_struct_equivalent(ksa_struct, ksb_struct)
661 1024 : CPASSERT(compatible_matrices)
662 :
663 16398 : IF (SUM(occb) == 0.0_dp) fb = 0.0_dp
664 :
665 16398 : DO icol_local = 1, ncol_local
666 15374 : icol_global = col_indices(icol_local)
667 15374 : j = INT(occa(icol_global)) + INT(occb(icol_global))
668 181549 : DO irow_local = 1, nrow_local
669 165151 : irow_global = row_indices(irow_local)
670 165151 : i = INT(occa(irow_global)) + INT(occb(irow_global))
671 : fa(irow_local, icol_local) = &
672 : roks_parameter(i, j, 1)*fa(irow_local, icol_local) + &
673 180525 : roks_parameter(i, j, 2)*fb(irow_local, icol_local)
674 : END DO
675 : END DO
676 :
677 1024 : CALL timestop(handle)
678 :
679 1024 : END SUBROUTINE combine_ks_matrices_1
680 :
681 : ! **************************************************************************************************
682 : !> \brief ...
683 : !> \param ksa ...
684 : !> \param ksb ...
685 : !> \param occa ...
686 : !> \param occb ...
687 : !> \param f ...
688 : !> \param nalpha ...
689 : !> \param nbeta ...
690 : ! **************************************************************************************************
691 0 : SUBROUTINE combine_ks_matrices_2(ksa, ksb, occa, occb, f, nalpha, nbeta)
692 :
693 : ! Combine the alpha and beta Kohn-Sham matrices during a restricted open
694 : ! Kohn-Sham (ROKS) calculation
695 : ! On input ksa and ksb contain the alpha and beta Kohn-Sham matrices,
696 : ! respectively. occa and occb contain the corresponding MO occupation
697 : ! numbers. On output the combined ROKS operator matrix is returned in ksa.
698 :
699 : ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
700 : ! - M. Filatov and S. Shaik, Chem. Phys. Lett. 288, 689 (1998)
701 :
702 : TYPE(cp_fm_type), INTENT(IN) :: ksa, ksb
703 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occa, occb
704 : REAL(KIND=dp), INTENT(IN) :: f
705 : INTEGER, INTENT(IN) :: nalpha, nbeta
706 :
707 : CHARACTER(LEN=*), PARAMETER :: routineN = 'combine_ks_matrices_2'
708 :
709 : INTEGER :: handle, icol_global, icol_local, &
710 : irow_global, irow_local, ncol_local, &
711 : nrow_local
712 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
713 : LOGICAL :: compatible_matrices
714 : REAL(KIND=dp) :: beta, t1, t2, ta, tb
715 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
716 0 : POINTER :: fa, fb
717 : TYPE(cp_fm_struct_type), POINTER :: ksa_struct, ksb_struct
718 :
719 : ! -------------------------------------------------------------------------
720 :
721 0 : CALL timeset(routineN, handle)
722 :
723 : CALL cp_fm_get_info(matrix=ksa, &
724 : matrix_struct=ksa_struct, &
725 : nrow_local=nrow_local, &
726 : ncol_local=ncol_local, &
727 : row_indices=row_indices, &
728 : col_indices=col_indices, &
729 0 : local_data=fa)
730 :
731 : CALL cp_fm_get_info(matrix=ksb, &
732 : matrix_struct=ksb_struct, &
733 0 : local_data=fb)
734 :
735 0 : compatible_matrices = cp_fm_struct_equivalent(ksa_struct, ksb_struct)
736 0 : CPASSERT(compatible_matrices)
737 :
738 0 : beta = 1.0_dp/(1.0_dp - f)
739 :
740 0 : DO icol_local = 1, ncol_local
741 :
742 0 : icol_global = col_indices(icol_local)
743 :
744 0 : DO irow_local = 1, nrow_local
745 :
746 0 : irow_global = row_indices(irow_local)
747 :
748 0 : t1 = 0.5_dp*(fa(irow_local, icol_local) + fb(irow_local, icol_local))
749 :
750 0 : IF ((0 < irow_global) .AND. (irow_global <= nbeta)) THEN
751 0 : IF ((0 < icol_global) .AND. (icol_global <= nbeta)) THEN
752 : ! closed-closed
753 0 : fa(irow_local, icol_local) = t1
754 0 : ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
755 : ! closed-open
756 0 : ta = 0.5_dp*(f - REAL(occa(icol_global), KIND=dp))/f
757 0 : tb = 0.5_dp*(f - REAL(occb(icol_global), KIND=dp))/f
758 0 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
759 0 : fa(irow_local, icol_local) = t1 + (beta - 1.0_dp)*t2
760 : ELSE
761 : ! closed-virtual
762 0 : fa(irow_local, icol_local) = t1
763 : END IF
764 0 : ELSE IF ((nbeta < irow_global) .AND. (irow_global <= nalpha)) THEN
765 : IF ((0 < irow_global) .AND. (irow_global <= nbeta)) THEN
766 : ! open-closed
767 : ta = 0.5_dp*(f - REAL(occa(irow_global), KIND=dp))/f
768 : tb = 0.5_dp*(f - REAL(occb(irow_global), KIND=dp))/f
769 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
770 : fa(irow_local, icol_local) = t1 + (beta - 1.0_dp)*t2
771 0 : ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
772 : ! open-open
773 0 : ta = 0.5_dp*(f - REAL(occa(icol_global), KIND=dp))/f
774 0 : tb = 0.5_dp*(f - REAL(occb(icol_global), KIND=dp))/f
775 0 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
776 0 : IF (irow_global == icol_global) THEN
777 0 : fa(irow_local, icol_local) = t1 - t2
778 : ELSE
779 0 : fa(irow_local, icol_local) = t1 - 0.5_dp*t2
780 : END IF
781 : ELSE
782 : ! open-virtual
783 0 : ta = 0.5_dp*(f - REAL(occa(irow_global), KIND=dp))/f
784 0 : tb = 0.5_dp*(f - REAL(occb(irow_global), KIND=dp))/f
785 0 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
786 0 : fa(irow_local, icol_local) = t1 - t2
787 : END IF
788 : ELSE
789 0 : IF ((0 < irow_global) .AND. (irow_global < nbeta)) THEN
790 : ! virtual-closed
791 0 : fa(irow_local, icol_local) = t1
792 0 : ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
793 : ! virtual-open
794 0 : ta = 0.5_dp*(f - REAL(occa(icol_global), KIND=dp))/f
795 0 : tb = 0.5_dp*(f - REAL(occb(icol_global), KIND=dp))/f
796 0 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
797 0 : fa(irow_local, icol_local) = t1 - t2
798 : ELSE
799 : ! virtual-virtual
800 0 : fa(irow_local, icol_local) = t1
801 : END IF
802 : END IF
803 :
804 : END DO
805 : END DO
806 :
807 0 : CALL timestop(handle)
808 :
809 0 : END SUBROUTINE combine_ks_matrices_2
810 :
811 : ! **************************************************************************************************
812 : !> \brief Correct MO eigenvalues after MO level shifting.
813 : !> \param mo_eigenvalues vector of eigenvalues
814 : !> \param homo index of the highest occupied molecular orbital
815 : !> \param nmo number of molecular orbitals
816 : !> \param level_shift amount of applied level shifting (in a.u.)
817 : !> \date 19.04.2002
818 : !> \par History
819 : !> - correct_mo_eigenvalues added (18.04.02,MK)
820 : !> - moved from module qs_mo_types, revised interface (03.2016, Sergey Chulkov)
821 : !> \author MK
822 : !> \version 1.0
823 : ! **************************************************************************************************
824 230 : PURE SUBROUTINE correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
825 :
826 : REAL(kind=dp), DIMENSION(:), INTENT(inout) :: mo_eigenvalues
827 : INTEGER, INTENT(in) :: homo, nmo
828 : REAL(kind=dp), INTENT(in) :: level_shift
829 :
830 : INTEGER :: imo
831 :
832 5720 : DO imo = homo + 1, nmo
833 5720 : mo_eigenvalues(imo) = mo_eigenvalues(imo) - level_shift
834 : END DO
835 :
836 230 : END SUBROUTINE correct_mo_eigenvalues
837 :
838 : ! **************************************************************************************************
839 : !> \brief Adjust the Kohn-Sham matrix by shifting the orbital energies of all
840 : !> unoccupied molecular orbitals
841 : !> \param matrix_ks_fm transformed Kohn-Sham matrix = U^{-1,T} * KS * U^{-1}
842 : !> \param mo_coeff matrix of molecular orbitals (C)
843 : !> \param homo number of occupied molecular orbitals
844 : !> \param level_shift amount of shift applying (in a.u.)
845 : !> \param is_triangular indicates that matrix_u_fm contains an upper triangular matrix
846 : !> \param matrix_u_fm matrix U: S (overlap matrix) = U^T * U;
847 : !> assume an identity matrix if omitted
848 : !> \par History
849 : !> 03.2016 created [Sergey Chulkov]
850 : ! **************************************************************************************************
851 230 : SUBROUTINE shift_unocc_mos(matrix_ks_fm, mo_coeff, homo, &
852 : level_shift, is_triangular, matrix_u_fm)
853 :
854 : TYPE(cp_fm_type), INTENT(IN) :: matrix_ks_fm, mo_coeff
855 : INTEGER, INTENT(in) :: homo
856 : REAL(kind=dp), INTENT(in) :: level_shift
857 : LOGICAL, INTENT(in) :: is_triangular
858 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_u_fm
859 :
860 : CHARACTER(len=*), PARAMETER :: routineN = 'shift_unocc_mos'
861 :
862 : INTEGER :: handle, nao, nao_red, nmo
863 230 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights
864 : TYPE(cp_fm_struct_type), POINTER :: ao_mo_fmstruct
865 : TYPE(cp_fm_type) :: u_mo, u_mo_scaled
866 :
867 230 : CALL timeset(routineN, handle)
868 :
869 230 : IF (PRESENT(matrix_u_fm)) THEN
870 230 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
871 230 : CALL cp_fm_get_info(matrix_u_fm, nrow_global=nao_red, ncol_global=nao)
872 : ELSE
873 0 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo, nrow_global=nao)
874 0 : nao_red = nao
875 : END IF
876 :
877 230 : NULLIFY (ao_mo_fmstruct)
878 : CALL cp_fm_struct_create(ao_mo_fmstruct, nrow_global=nao_red, ncol_global=nmo, &
879 230 : para_env=mo_coeff%matrix_struct%para_env, context=mo_coeff%matrix_struct%context)
880 :
881 230 : CALL cp_fm_create(u_mo, ao_mo_fmstruct)
882 230 : CALL cp_fm_create(u_mo_scaled, ao_mo_fmstruct)
883 :
884 230 : CALL cp_fm_struct_release(ao_mo_fmstruct)
885 :
886 : ! U * C
887 230 : IF (PRESENT(matrix_u_fm)) THEN
888 230 : IF (is_triangular) THEN
889 144 : CALL cp_fm_to_fm(mo_coeff, u_mo)
890 : CALL cp_fm_triangular_multiply(matrix_u_fm, u_mo, side="L", transpose_tr=.FALSE., &
891 144 : invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nmo, alpha=1.0_dp)
892 : ELSE
893 86 : CALL parallel_gemm("N", "N", nao_red, nmo, nao, 1.0_dp, matrix_u_fm, mo_coeff, 0.0_dp, u_mo)
894 : END IF
895 : ELSE
896 : ! assume U is an identity matrix
897 0 : CALL cp_fm_to_fm(mo_coeff, u_mo)
898 : END IF
899 :
900 230 : CALL cp_fm_to_fm(u_mo, u_mo_scaled)
901 :
902 : ! up-shift all unoccupied molecular orbitals by the amount of 'level_shift'
903 : ! weight = diag(DELTA) = (0, ... 0, level_shift, ..., level_shift)
904 : ! MO index : 1 .. homo homo+1 ... nmo
905 690 : ALLOCATE (weights(nmo))
906 1170 : weights(1:homo) = 0.0_dp
907 6010 : weights(homo + 1:nmo) = level_shift
908 : ! DELTA * U * C
909 : ! DELTA is a diagonal matrix, so simply scale all the columns of (U * C) by weights(:)
910 230 : CALL cp_fm_column_scale(u_mo_scaled, weights)
911 230 : DEALLOCATE (weights)
912 :
913 : ! NewKS = U^{-1,T} * KS * U^{-1} + (U * C) * DELTA * (U * C)^T
914 230 : CALL parallel_gemm("N", "T", nao_red, nao_red, nmo, 1.0_dp, u_mo, u_mo_scaled, 1.0_dp, matrix_ks_fm)
915 :
916 230 : CALL cp_fm_release(u_mo_scaled)
917 230 : CALL cp_fm_release(u_mo)
918 :
919 230 : CALL timestop(handle)
920 :
921 460 : END SUBROUTINE shift_unocc_mos
922 :
923 : END MODULE qs_scf_methods
|
