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