(git:374b731)
Loading...
Searching...
No Matches
qs_scf_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief 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! **************************************************************************************************
19MODULE qs_scf_methods
20
21 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
22 cp_dbcsr_sm_fm_multiply
23 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
24 cp_fm_symm,&
25 cp_fm_triangular_multiply,&
26 cp_fm_upper_to_full
27 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_reduce,&
28 cp_fm_cholesky_restore
29 USE cp_fm_diag, ONLY: choose_eigv_solver,&
30 cp_fm_block_jacobi
31 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
32 cp_fm_struct_equivalent,&
33 cp_fm_struct_release,&
34 cp_fm_struct_type
35 USE cp_fm_types, ONLY: cp_fm_create,&
36 cp_fm_get_info,&
37 cp_fm_release,&
38 cp_fm_to_fm,&
39 cp_fm_type
40 USE dbcsr_api, ONLY: &
41 dbcsr_add, dbcsr_desymmetrize, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
42 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
43 dbcsr_multiply, dbcsr_p_type, dbcsr_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_dbcsr, &
67 eigensolver_symm, &
68 eigensolver_simple, &
69 scf_env_density_mixing
70
71 PRIVATE :: correct_mo_eigenvalues, shift_unocc_mos
72
73 INTERFACE combine_ks_matrices
74 MODULE PROCEDURE combine_ks_matrices_1, &
75 combine_ks_matrices_2
76 END INTERFACE combine_ks_matrices
77
78CONTAINS
79
80! **************************************************************************************************
81!> \brief perform (if requested) a density mixing
82!> \param p_mix_new New density matrices
83!> \param mixing_store ...
84!> \param rho_ao Density environment
85!> \param para_env ...
86!> \param iter_delta ...
87!> \param iter_count ...
88!> \param diis ...
89!> \param invert Invert mixing
90!> \par History
91!> 02.2003 created [fawzi]
92!> 08.2014 adapted for kpoints [JGH]
93!> \author fawzi
94! **************************************************************************************************
95 SUBROUTINE scf_env_density_mixing(p_mix_new, mixing_store, rho_ao, para_env, &
96 iter_delta, iter_count, diis, invert)
97 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: p_mix_new
98 TYPE(mixing_storage_type), POINTER :: mixing_store
99 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
100 TYPE(mp_para_env_type), POINTER :: para_env
101 REAL(kind=dp), INTENT(INOUT) :: iter_delta
102 INTEGER, INTENT(IN) :: iter_count
103 LOGICAL, INTENT(in), OPTIONAL :: diis, invert
104
105 CHARACTER(len=*), PARAMETER :: routinen = 'scf_env_density_mixing'
106
107 INTEGER :: handle, ic, ispin
108 LOGICAL :: my_diis, my_invert
109 REAL(kind=dp) :: my_p_mix, tmp
110
111 CALL timeset(routinen, handle)
112
113 my_diis = .false.
114 IF (PRESENT(diis)) my_diis = diis
115 my_invert = .false.
116 IF (PRESENT(invert)) my_invert = invert
117 my_p_mix = mixing_store%alpha
118 IF (my_diis .OR. iter_count < mixing_store%nskip_mixing) THEN
119 my_p_mix = 1.0_dp
120 END IF
121
122 iter_delta = 0.0_dp
123 cpassert(ASSOCIATED(p_mix_new))
124 DO ic = 1, SIZE(p_mix_new, 2)
125 DO ispin = 1, SIZE(p_mix_new, 1)
126 IF (my_invert) THEN
127 cpassert(my_p_mix /= 0.0_dp)
128 IF (my_p_mix /= 1.0_dp) THEN
129 CALL dbcsr_add(matrix_a=p_mix_new(ispin, ic)%matrix, &
130 alpha_scalar=1.0_dp/my_p_mix, &
131 matrix_b=rho_ao(ispin, ic)%matrix, &
132 beta_scalar=(my_p_mix - 1.0_dp)/my_p_mix)
133 END IF
134 ELSE
135 CALL cp_sm_mix(m1=p_mix_new(ispin, ic)%matrix, &
136 m2=rho_ao(ispin, ic)%matrix, &
137 p_mix=my_p_mix, &
138 delta=tmp, &
139 para_env=para_env)
140 iter_delta = max(iter_delta, tmp)
141 END IF
142 END DO
143 END DO
144
145 CALL timestop(handle)
146
147 END SUBROUTINE scf_env_density_mixing
148
149! **************************************************************************************************
150!> \brief Diagonalise the Kohn-Sham matrix to get a new set of MO eigen-
151!> vectors and MO eigenvalues. ks will be modified
152!> \param matrix_ks_fm ...
153!> \param mo_set ...
154!> \param ortho ...
155!> \param work ...
156!> \param cholesky_method ...
157!> \param do_level_shift activate the level shifting technique
158!> \param level_shift amount of shift applied (in a.u.)
159!> \param matrix_u_fm matrix U : S (overlap matrix) = U^T * U
160!> \param use_jacobi ...
161!> \date 01.05.2001
162!> \author Matthias Krack
163!> \version 1.0
164! **************************************************************************************************
165 SUBROUTINE eigensolver(matrix_ks_fm, mo_set, ortho, work, &
166 cholesky_method, do_level_shift, &
167 level_shift, matrix_u_fm, use_jacobi)
168 TYPE(cp_fm_type), INTENT(IN) :: matrix_ks_fm
169 TYPE(mo_set_type), INTENT(IN) :: mo_set
170 TYPE(cp_fm_type), INTENT(IN) :: ortho, work
171 INTEGER, INTENT(inout) :: cholesky_method
172 LOGICAL, INTENT(in) :: do_level_shift
173 REAL(kind=dp), INTENT(in) :: level_shift
174 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_u_fm
175 LOGICAL, INTENT(in) :: use_jacobi
176
177 CHARACTER(len=*), PARAMETER :: routinen = 'eigensolver'
178
179 INTEGER :: handle, homo, nao, nmo
180 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
181 TYPE(cp_fm_type), POINTER :: mo_coeff
182
183 CALL timeset(routinen, handle)
184
185 NULLIFY (mo_coeff)
186 NULLIFY (mo_eigenvalues)
187
188 ! Diagonalise the Kohn-Sham matrix
189
190 CALL get_mo_set(mo_set=mo_set, &
191 nao=nao, &
192 nmo=nmo, &
193 homo=homo, &
194 eigenvalues=mo_eigenvalues, &
195 mo_coeff=mo_coeff)
196
197 SELECT CASE (cholesky_method)
198 CASE (cholesky_reduce)
199 CALL cp_fm_cholesky_reduce(matrix_ks_fm, ortho)
200
201 IF (do_level_shift) &
202 CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, nmo=nmo, nao=nao, &
203 level_shift=level_shift, is_triangular=.true., matrix_u_fm=matrix_u_fm)
204
205 CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
206 CALL cp_fm_cholesky_restore(work, nmo, ortho, mo_coeff, "SOLVE")
207 IF (do_level_shift) &
208 CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
209
210 CASE (cholesky_restore)
211 CALL cp_fm_upper_to_full(matrix_ks_fm, work)
212 CALL cp_fm_cholesky_restore(matrix_ks_fm, nao, ortho, work, &
213 "SOLVE", pos="RIGHT")
214 CALL cp_fm_cholesky_restore(work, nao, ortho, matrix_ks_fm, &
215 "SOLVE", pos="LEFT", transa="T")
216
217 IF (do_level_shift) &
218 CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, nmo=nmo, nao=nao, &
219 level_shift=level_shift, is_triangular=.true., matrix_u_fm=matrix_u_fm)
220
221 CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
222 CALL cp_fm_cholesky_restore(work, nmo, ortho, mo_coeff, "SOLVE")
223
224 IF (do_level_shift) &
225 CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
226
227 CASE (cholesky_inverse)
228 CALL cp_fm_upper_to_full(matrix_ks_fm, work)
229
230 CALL cp_fm_triangular_multiply(ortho, matrix_ks_fm, side="R", transpose_tr=.false., &
231 invert_tr=.false., uplo_tr="U", n_rows=nao, n_cols=nao, alpha=1.0_dp)
232 CALL cp_fm_triangular_multiply(ortho, matrix_ks_fm, side="L", transpose_tr=.true., &
233 invert_tr=.false., uplo_tr="U", n_rows=nao, n_cols=nao, alpha=1.0_dp)
234
235 IF (do_level_shift) &
236 CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, nmo=nmo, nao=nao, &
237 level_shift=level_shift, is_triangular=.true., matrix_u_fm=matrix_u_fm)
238
239 CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
240 CALL cp_fm_triangular_multiply(ortho, work, side="L", transpose_tr=.false., &
241 invert_tr=.false., uplo_tr="U", n_rows=nao, n_cols=nmo, alpha=1.0_dp)
242 CALL cp_fm_to_fm(work, mo_coeff, nmo, 1, 1)
243
244 IF (do_level_shift) &
245 CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
246
247 END SELECT
248
249 IF (use_jacobi) THEN
250 CALL cp_fm_to_fm(mo_coeff, ortho)
251 cholesky_method = cholesky_off
252 END IF
253
254 CALL timestop(handle)
255
256 END SUBROUTINE eigensolver
257
258! **************************************************************************************************
259!> \brief ...
260!> \param matrix_ks ...
261!> \param matrix_ks_fm ...
262!> \param mo_set ...
263!> \param ortho_dbcsr ...
264!> \param ksbuf1 ...
265!> \param ksbuf2 ...
266! **************************************************************************************************
267 SUBROUTINE eigensolver_dbcsr(matrix_ks, matrix_ks_fm, mo_set, ortho_dbcsr, ksbuf1, ksbuf2)
268 TYPE(dbcsr_type), POINTER :: matrix_ks
269 TYPE(cp_fm_type), INTENT(IN) :: matrix_ks_fm
270 TYPE(mo_set_type), INTENT(IN) :: mo_set
271 TYPE(dbcsr_type), POINTER :: ortho_dbcsr, ksbuf1, ksbuf2
272
273 CHARACTER(len=*), PARAMETER :: routinen = 'eigensolver_dbcsr'
274
275 INTEGER :: handle, nao, nmo
276 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
277 TYPE(cp_fm_type) :: all_evecs, nmo_evecs
278 TYPE(cp_fm_type), POINTER :: mo_coeff
279
280 CALL timeset(routinen, handle)
281
282 NULLIFY (mo_coeff)
283 NULLIFY (mo_eigenvalues)
284
285 CALL get_mo_set(mo_set=mo_set, &
286 nao=nao, &
287 nmo=nmo, &
288 eigenvalues=mo_eigenvalues, &
289 mo_coeff=mo_coeff)
290
291! Reduce KS matrix
292 CALL dbcsr_desymmetrize(matrix_ks, ksbuf2)
293 CALL dbcsr_multiply('N', 'N', 1.0_dp, ksbuf2, ortho_dbcsr, 0.0_dp, ksbuf1)
294 CALL dbcsr_multiply('T', 'N', 1.0_dp, ortho_dbcsr, ksbuf1, 0.0_dp, ksbuf2)
295
296! Solve the eigenvalue problem
297 CALL copy_dbcsr_to_fm(ksbuf2, matrix_ks_fm)
298 CALL cp_fm_create(all_evecs, matrix_ks_fm%matrix_struct)
299 CALL choose_eigv_solver(matrix_ks_fm, all_evecs, mo_eigenvalues)
300
301 ! select first nmo eigenvectors
302 CALL cp_fm_create(nmo_evecs, mo_coeff%matrix_struct)
303 CALL cp_fm_to_fm(msource=all_evecs, mtarget=nmo_evecs, ncol=nmo)
304 CALL cp_fm_release(all_evecs)
305
306! Restore the eigenvector of the general eig. problem
307 CALL cp_dbcsr_sm_fm_multiply(ortho_dbcsr, nmo_evecs, mo_coeff, nmo)
308
309 CALL cp_fm_release(nmo_evecs)
310 CALL timestop(handle)
311
312 END SUBROUTINE eigensolver_dbcsr
313
314! **************************************************************************************************
315!> \brief ...
316!> \param matrix_ks_fm ...
317!> \param mo_set ...
318!> \param ortho ...
319!> \param work ...
320!> \param do_level_shift activate the level shifting technique
321!> \param level_shift amount of shift applied (in a.u.)
322!> \param matrix_u_fm matrix U : S (overlap matrix) = U^T * U
323!> \param use_jacobi ...
324!> \param jacobi_threshold ...
325! **************************************************************************************************
326 SUBROUTINE eigensolver_symm(matrix_ks_fm, mo_set, ortho, work, do_level_shift, &
327 level_shift, matrix_u_fm, use_jacobi, jacobi_threshold)
328 TYPE(cp_fm_type), INTENT(IN) :: matrix_ks_fm
329 TYPE(mo_set_type), INTENT(IN) :: mo_set
330 TYPE(cp_fm_type), INTENT(IN) :: ortho, work
331 LOGICAL, INTENT(IN) :: do_level_shift
332 REAL(kind=dp), INTENT(IN) :: level_shift
333 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_u_fm
334 LOGICAL, INTENT(IN) :: use_jacobi
335 REAL(kind=dp), INTENT(IN) :: jacobi_threshold
336
337 CHARACTER(len=*), PARAMETER :: routinen = 'eigensolver_symm'
338
339 INTEGER :: handle, homo, nao, nelectron, nmo
340 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
341 TYPE(cp_fm_type), POINTER :: mo_coeff
342
343 CALL timeset(routinen, handle)
344
345 NULLIFY (mo_coeff)
346 NULLIFY (mo_eigenvalues)
347
348 ! Diagonalise the Kohn-Sham matrix
349
350 CALL get_mo_set(mo_set=mo_set, &
351 nao=nao, &
352 nmo=nmo, &
353 homo=homo, &
354 nelectron=nelectron, &
355 eigenvalues=mo_eigenvalues, &
356 mo_coeff=mo_coeff)
357
358 IF (use_jacobi) THEN
359
360 CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, matrix_ks_fm, mo_coeff, 0.0_dp, work)
361 CALL parallel_gemm("T", "N", homo, nao - homo, nao, 1.0_dp, work, mo_coeff, &
362 0.0_dp, matrix_ks_fm, b_first_col=homo + 1, c_first_col=homo + 1)
363
364 ! Block Jacobi (pseudo-diagonalization, only one sweep)
365 CALL cp_fm_block_jacobi(matrix_ks_fm, mo_coeff, mo_eigenvalues, &
366 jacobi_threshold, homo + 1)
367
368 ELSE ! full S^(-1/2) has been computed
369
370 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, matrix_ks_fm, ortho, 0.0_dp, work)
371 CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, ortho, work, 0.0_dp, matrix_ks_fm)
372
373 IF (do_level_shift) &
374 CALL shift_unocc_mos(matrix_ks_fm=matrix_ks_fm, mo_coeff=mo_coeff, homo=homo, nmo=nmo, nao=nao, &
375 level_shift=level_shift, is_triangular=.false., matrix_u_fm=matrix_u_fm)
376
377 CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
378
379 CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, &
380 mo_coeff)
381
382 IF (do_level_shift) &
383 CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
384
385 END IF
386
387 CALL timestop(handle)
388
389 END SUBROUTINE eigensolver_symm
390
391! **************************************************************************************************
392
393! **************************************************************************************************
394!> \brief ...
395!> \param matrix_ks ...
396!> \param mo_set ...
397!> \param work ...
398!> \param do_level_shift activate the level shifting technique
399!> \param level_shift amount of shift applied (in a.u.)
400!> \param use_jacobi ...
401!> \param jacobi_threshold ...
402! **************************************************************************************************
403 SUBROUTINE eigensolver_simple(matrix_ks, mo_set, work, do_level_shift, &
404 level_shift, use_jacobi, jacobi_threshold)
405
406 TYPE(cp_fm_type), INTENT(IN) :: matrix_ks
407 TYPE(mo_set_type), INTENT(IN) :: mo_set
408 TYPE(cp_fm_type), INTENT(IN) :: work
409 LOGICAL, INTENT(IN) :: do_level_shift
410 REAL(kind=dp), INTENT(IN) :: level_shift
411 LOGICAL, INTENT(IN) :: use_jacobi
412 REAL(kind=dp), INTENT(IN) :: jacobi_threshold
413
414 CHARACTER(len=*), PARAMETER :: routinen = 'eigensolver_simple'
415
416 INTEGER :: handle, homo, nao, nelectron, nmo
417 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
418 TYPE(cp_fm_type), POINTER :: mo_coeff
419
420 CALL timeset(routinen, handle)
421
422 NULLIFY (mo_coeff)
423 NULLIFY (mo_eigenvalues)
424
425 CALL get_mo_set(mo_set=mo_set, &
426 nao=nao, &
427 nmo=nmo, &
428 homo=homo, &
429 nelectron=nelectron, &
430 eigenvalues=mo_eigenvalues, &
431 mo_coeff=mo_coeff)
432
433 IF (do_level_shift) THEN
434 ! matrix_u_fm is simply an identity matrix, so we omit it here
435 CALL shift_unocc_mos(matrix_ks_fm=matrix_ks, mo_coeff=mo_coeff, homo=homo, nmo=nmo, nao=nao, &
436 level_shift=level_shift, is_triangular=.false.)
437 END IF
438
439 IF (use_jacobi) THEN
440 CALL cp_fm_symm("L", "U", nao, homo, 1.0_dp, matrix_ks, mo_coeff, 0.0_dp, work)
441 CALL parallel_gemm("T", "N", homo, nao - homo, nao, 1.0_dp, work, mo_coeff, &
442 0.0_dp, matrix_ks, b_first_col=homo + 1, c_first_col=homo + 1)
443 ! Block Jacobi (pseudo-diagonalization, only one sweep)
444 CALL cp_fm_block_jacobi(matrix_ks, mo_coeff, mo_eigenvalues, jacobi_threshold, homo + 1)
445 ELSE
446
447 CALL choose_eigv_solver(matrix_ks, work, mo_eigenvalues)
448
449 CALL cp_fm_to_fm(work, mo_coeff, nmo, 1, 1)
450
451 END IF
452
453 IF (do_level_shift) &
454 CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
455
456 CALL timestop(handle)
457
458 END SUBROUTINE eigensolver_simple
459
460! **************************************************************************************************
461!> \brief Perform a mixing of the given matrixes into the first matrix
462!> m1 = m2 + p_mix (m1-m2)
463!> \param m1 first (new) matrix, is modified
464!> \param m2 the second (old) matrix
465!> \param p_mix how much m1 is conserved (0: none, 1: all)
466!> \param delta maximum norm of m1-m2
467!> \param para_env ...
468!> \param m3 ...
469!> \par History
470!> 02.2003 rewamped [fawzi]
471!> \author fawzi
472!> \note
473!> if you what to store the result in m2 swap m1 and m2 an use
474!> (1-pmix) as pmix
475!> para_env should be removed (embedded in matrix)
476! **************************************************************************************************
477 SUBROUTINE cp_sm_mix(m1, m2, p_mix, delta, para_env, m3)
478
479 TYPE(dbcsr_type), POINTER :: m1, m2
480 REAL(kind=dp), INTENT(IN) :: p_mix
481 REAL(kind=dp), INTENT(OUT) :: delta
482 TYPE(mp_para_env_type), POINTER :: para_env
483 TYPE(dbcsr_type), OPTIONAL, POINTER :: m3
484
485 CHARACTER(len=*), PARAMETER :: routinen = 'cp_sm_mix'
486
487 INTEGER :: blk, handle, i, iblock_col, iblock_row, j
488 LOGICAL :: found
489 REAL(kind=dp), DIMENSION(:, :), POINTER :: p_delta_block, p_new_block, p_old_block
490 TYPE(dbcsr_iterator_type) :: iter
491
492 CALL timeset(routinen, handle)
493 delta = 0.0_dp
494
495 CALL dbcsr_iterator_start(iter, m1)
496 DO WHILE (dbcsr_iterator_blocks_left(iter))
497 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, p_new_block, blk)
498 CALL dbcsr_get_block_p(matrix=m2, row=iblock_row, col=iblock_col, &
499 block=p_old_block, found=found)
500 cpassert(ASSOCIATED(p_old_block))
501 IF (PRESENT(m3)) THEN
502 CALL dbcsr_get_block_p(matrix=m3, row=iblock_row, col=iblock_col, &
503 block=p_delta_block, found=found)
504 cpassert(ASSOCIATED(p_delta_block))
505
506 DO j = 1, SIZE(p_new_block, 2)
507 DO i = 1, SIZE(p_new_block, 1)
508 p_delta_block(i, j) = p_new_block(i, j) - p_old_block(i, j)
509 delta = max(delta, abs(p_delta_block(i, j)))
510 END DO
511 END DO
512 ELSE
513 DO j = 1, SIZE(p_new_block, 2)
514 DO i = 1, SIZE(p_new_block, 1)
515 p_new_block(i, j) = p_new_block(i, j) - p_old_block(i, j)
516 delta = max(delta, abs(p_new_block(i, j)))
517 p_new_block(i, j) = p_old_block(i, j) + p_mix*p_new_block(i, j)
518 END DO
519 END DO
520 END IF
521 END DO
522 CALL dbcsr_iterator_stop(iter)
523
524 CALL para_env%max(delta)
525
526 CALL timestop(handle)
527
528 END SUBROUTINE cp_sm_mix
529
530! **************************************************************************************************
531!> \brief ...
532!> \param ksa ...
533!> \param ksb ...
534!> \param occa ...
535!> \param occb ...
536!> \param roks_parameter ...
537! **************************************************************************************************
538 SUBROUTINE combine_ks_matrices_1(ksa, ksb, occa, occb, roks_parameter)
539
540 ! Combine the alpha and beta Kohn-Sham matrices during a restricted open
541 ! Kohn-Sham (ROKS) calculation
542 ! On input ksa and ksb contain the alpha and beta Kohn-Sham matrices,
543 ! respectively. occa and occb contain the corresponding MO occupation
544 ! numbers. On output the combined ROKS operator matrix is returned in ksa.
545
546 ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
547 ! - M. F. Guest and V. R. Saunders, Mol. Phys. 28(3), 819 (1974)
548
549 TYPE(cp_fm_type), INTENT(IN) :: ksa, ksb
550 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occa, occb
551 REAL(KIND=dp), DIMENSION(0:2, 0:2, 1:2), &
552 INTENT(IN) :: roks_parameter
553
554 CHARACTER(LEN=*), PARAMETER :: routineN = 'combine_ks_matrices_1'
555
556 INTEGER :: handle, i, icol_global, icol_local, &
557 irow_global, irow_local, j, &
558 ncol_local, nrow_local
559 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
560 LOGICAL :: compatible_matrices
561 REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
562 POINTER :: fa, fb
563 TYPE(cp_fm_struct_type), POINTER :: ksa_struct, ksb_struct
564
565! -------------------------------------------------------------------------
566
567 CALL timeset(routinen, handle)
568
569 CALL cp_fm_get_info(matrix=ksa, &
570 matrix_struct=ksa_struct, &
571 nrow_local=nrow_local, &
572 ncol_local=ncol_local, &
573 row_indices=row_indices, &
574 col_indices=col_indices, &
575 local_data=fa)
576
577 CALL cp_fm_get_info(matrix=ksb, &
578 matrix_struct=ksb_struct, &
579 local_data=fb)
580
581 compatible_matrices = cp_fm_struct_equivalent(ksa_struct, ksb_struct)
582 cpassert(compatible_matrices)
583
584 IF (sum(occb) == 0.0_dp) fb = 0.0_dp
585
586 DO icol_local = 1, ncol_local
587 icol_global = col_indices(icol_local)
588 j = int(occa(icol_global)) + int(occb(icol_global))
589 DO irow_local = 1, nrow_local
590 irow_global = row_indices(irow_local)
591 i = int(occa(irow_global)) + int(occb(irow_global))
592 fa(irow_local, icol_local) = &
593 roks_parameter(i, j, 1)*fa(irow_local, icol_local) + &
594 roks_parameter(i, j, 2)*fb(irow_local, icol_local)
595 END DO
596 END DO
597
598 CALL timestop(handle)
599
600 END SUBROUTINE combine_ks_matrices_1
601
602! **************************************************************************************************
603!> \brief ...
604!> \param ksa ...
605!> \param ksb ...
606!> \param occa ...
607!> \param occb ...
608!> \param f ...
609!> \param nalpha ...
610!> \param nbeta ...
611! **************************************************************************************************
612 SUBROUTINE combine_ks_matrices_2(ksa, ksb, occa, occb, f, nalpha, nbeta)
613
614 ! Combine the alpha and beta Kohn-Sham matrices during a restricted open
615 ! Kohn-Sham (ROKS) calculation
616 ! On input ksa and ksb contain the alpha and beta Kohn-Sham matrices,
617 ! respectively. occa and occb contain the corresponding MO occupation
618 ! numbers. On output the combined ROKS operator matrix is returned in ksa.
619
620 ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
621 ! - M. Filatov and S. Shaik, Chem. Phys. Lett. 288, 689 (1998)
622
623 TYPE(cp_fm_type), INTENT(IN) :: ksa, ksb
624 REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occa, occb
625 REAL(KIND=dp), INTENT(IN) :: f
626 INTEGER, INTENT(IN) :: nalpha, nbeta
627
628 CHARACTER(LEN=*), PARAMETER :: routineN = 'combine_ks_matrices_2'
629
630 INTEGER :: handle, icol_global, icol_local, &
631 irow_global, irow_local, ncol_local, &
632 nrow_local
633 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
634 LOGICAL :: compatible_matrices
635 REAL(KIND=dp) :: beta, t1, t2, ta, tb
636 REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
637 POINTER :: fa, fb
638 TYPE(cp_fm_struct_type), POINTER :: ksa_struct, ksb_struct
639
640! -------------------------------------------------------------------------
641
642 CALL timeset(routinen, handle)
643
644 CALL cp_fm_get_info(matrix=ksa, &
645 matrix_struct=ksa_struct, &
646 nrow_local=nrow_local, &
647 ncol_local=ncol_local, &
648 row_indices=row_indices, &
649 col_indices=col_indices, &
650 local_data=fa)
651
652 CALL cp_fm_get_info(matrix=ksb, &
653 matrix_struct=ksb_struct, &
654 local_data=fb)
655
656 compatible_matrices = cp_fm_struct_equivalent(ksa_struct, ksb_struct)
657 cpassert(compatible_matrices)
658
659 beta = 1.0_dp/(1.0_dp - f)
660
661 DO icol_local = 1, ncol_local
662
663 icol_global = col_indices(icol_local)
664
665 DO irow_local = 1, nrow_local
666
667 irow_global = row_indices(irow_local)
668
669 t1 = 0.5_dp*(fa(irow_local, icol_local) + fb(irow_local, icol_local))
670
671 IF ((0 < irow_global) .AND. (irow_global <= nbeta)) THEN
672 IF ((0 < icol_global) .AND. (icol_global <= nbeta)) THEN
673 ! closed-closed
674 fa(irow_local, icol_local) = t1
675 ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
676 ! closed-open
677 ta = 0.5_dp*(f - real(occa(icol_global), kind=dp))/f
678 tb = 0.5_dp*(f - real(occb(icol_global), kind=dp))/f
679 t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
680 fa(irow_local, icol_local) = t1 + (beta - 1.0_dp)*t2
681 ELSE
682 ! closed-virtual
683 fa(irow_local, icol_local) = t1
684 END IF
685 ELSE IF ((nbeta < irow_global) .AND. (irow_global <= nalpha)) THEN
686 IF ((0 < irow_global) .AND. (irow_global <= nbeta)) THEN
687 ! open-closed
688 ta = 0.5_dp*(f - real(occa(irow_global), kind=dp))/f
689 tb = 0.5_dp*(f - real(occb(irow_global), kind=dp))/f
690 t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
691 fa(irow_local, icol_local) = t1 + (beta - 1.0_dp)*t2
692 ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
693 ! open-open
694 ta = 0.5_dp*(f - real(occa(icol_global), kind=dp))/f
695 tb = 0.5_dp*(f - real(occb(icol_global), kind=dp))/f
696 t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
697 IF (irow_global == icol_global) THEN
698 fa(irow_local, icol_local) = t1 - t2
699 ELSE
700 fa(irow_local, icol_local) = t1 - 0.5_dp*t2
701 END IF
702 ELSE
703 ! open-virtual
704 ta = 0.5_dp*(f - real(occa(irow_global), kind=dp))/f
705 tb = 0.5_dp*(f - real(occb(irow_global), kind=dp))/f
706 t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
707 fa(irow_local, icol_local) = t1 - t2
708 END IF
709 ELSE
710 IF ((0 < irow_global) .AND. (irow_global < nbeta)) THEN
711 ! virtual-closed
712 fa(irow_local, icol_local) = t1
713 ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
714 ! virtual-open
715 ta = 0.5_dp*(f - real(occa(icol_global), kind=dp))/f
716 tb = 0.5_dp*(f - real(occb(icol_global), kind=dp))/f
717 t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
718 fa(irow_local, icol_local) = t1 - t2
719 ELSE
720 ! virtual-virtual
721 fa(irow_local, icol_local) = t1
722 END IF
723 END IF
724
725 END DO
726 END DO
727
728 CALL timestop(handle)
729
730 END SUBROUTINE combine_ks_matrices_2
731
732! **************************************************************************************************
733!> \brief Correct MO eigenvalues after MO level shifting.
734!> \param mo_eigenvalues vector of eigenvalues
735!> \param homo index of the highest occupied molecular orbital
736!> \param nmo number of molecular orbitals
737!> \param level_shift amount of applied level shifting (in a.u.)
738!> \date 19.04.2002
739!> \par History
740!> - correct_mo_eigenvalues added (18.04.02,MK)
741!> - moved from module qs_mo_types, revised interface (03.2016, Sergey Chulkov)
742!> \author MK
743!> \version 1.0
744! **************************************************************************************************
745 PURE SUBROUTINE correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
746
747 REAL(kind=dp), DIMENSION(:), INTENT(inout) :: mo_eigenvalues
748 INTEGER, INTENT(in) :: homo, nmo
749 REAL(kind=dp), INTENT(in) :: level_shift
750
751 INTEGER :: imo
752
753 DO imo = homo + 1, nmo
754 mo_eigenvalues(imo) = mo_eigenvalues(imo) - level_shift
755 END DO
756
757 END SUBROUTINE correct_mo_eigenvalues
758
759! **************************************************************************************************
760!> \brief Adjust the Kohn-Sham matrix by shifting the orbital energies of all
761!> unoccupied molecular orbitals
762!> \param matrix_ks_fm transformed Kohn-Sham matrix = U^{-1,T} * KS * U^{-1}
763!> \param mo_coeff matrix of molecular orbitals (C)
764!> \param homo number of occupied molecular orbitals
765!> \param nmo number of molecular orbitals
766!> \param nao number of atomic orbitals
767!> \param level_shift amount of shift applying (in a.u.)
768!> \param is_triangular indicates that matrix_u_fm contains an upper triangular matrix
769!> \param matrix_u_fm matrix U: S (overlap matrix) = U^T * U;
770!> assume an identity matrix if omitted
771!> \par History
772!> 03.2016 created [Sergey Chulkov]
773! **************************************************************************************************
774 SUBROUTINE shift_unocc_mos(matrix_ks_fm, mo_coeff, homo, nmo, nao, &
775 level_shift, is_triangular, matrix_u_fm)
776
777 TYPE(cp_fm_type), INTENT(IN) :: matrix_ks_fm, mo_coeff
778 INTEGER, INTENT(in) :: homo, nmo, nao
779 REAL(kind=dp), INTENT(in) :: level_shift
780 LOGICAL, INTENT(in) :: is_triangular
781 TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: matrix_u_fm
782
783 CHARACTER(len=*), PARAMETER :: routineN = 'shift_unocc_mos'
784
785 INTEGER :: handle
786 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: weights
787 TYPE(cp_fm_struct_type), POINTER :: ao_mo_fmstruct
788 TYPE(cp_fm_type) :: u_mo, u_mo_scaled
789
790 CALL timeset(routinen, handle)
791
792 NULLIFY (ao_mo_fmstruct)
793 CALL cp_fm_struct_create(ao_mo_fmstruct, nrow_global=nao, ncol_global=nmo, &
794 para_env=mo_coeff%matrix_struct%para_env, context=mo_coeff%matrix_struct%context)
795
796 CALL cp_fm_create(u_mo, ao_mo_fmstruct)
797 CALL cp_fm_create(u_mo_scaled, ao_mo_fmstruct)
798
799 CALL cp_fm_struct_release(ao_mo_fmstruct)
800
801 ! U * C
802 IF (PRESENT(matrix_u_fm)) THEN
803 IF (is_triangular) THEN
804 CALL cp_fm_to_fm(mo_coeff, u_mo)
805 CALL cp_fm_triangular_multiply(matrix_u_fm, u_mo, side="L", transpose_tr=.false., &
806 invert_tr=.false., uplo_tr="U", n_rows=nao, n_cols=nmo, alpha=1.0_dp)
807 ELSE
808 CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, matrix_u_fm, mo_coeff, 0.0_dp, u_mo)
809 END IF
810 ELSE
811 ! assume U is an identity matrix
812 CALL cp_fm_to_fm(mo_coeff, u_mo)
813 END IF
814
815 CALL cp_fm_to_fm(u_mo, u_mo_scaled)
816
817 ! up-shift all unoccupied molecular orbitals by the amount of 'level_shift'
818 ! weight = diag(DELTA) = (0, ... 0, level_shift, ..., level_shift)
819 ! MO index : 1 .. homo homo+1 ... nmo
820 ALLOCATE (weights(nmo))
821 weights(1:homo) = 0.0_dp
822 weights(homo + 1:nmo) = level_shift
823 ! DELTA * U * C
824 ! DELTA is a diagonal matrix, so simply scale all the columns of (U * C) by weights(:)
825 CALL cp_fm_column_scale(u_mo_scaled, weights)
826 DEALLOCATE (weights)
827
828 ! NewKS = U^{-1,T} * KS * U^{-1} + (U * C) * DELTA * (U * C)^T
829 CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, u_mo, u_mo_scaled, 1.0_dp, matrix_ks_fm)
830
831 CALL cp_fm_release(u_mo_scaled)
832 CALL cp_fm_release(u_mo)
833
834 CALL timestop(handle)
835
836 END SUBROUTINE shift_unocc_mos
837
838END MODULE qs_scf_methods
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:85
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
qs_scf_methods::combine_ks_matrices
Definition qs_scf_methods.F:73
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:744
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:208
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_column_scale
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
Definition cp_fm_basic_linalg.F:1882
cp_fm_basic_linalg::cp_fm_upper_to_full
subroutine, public cp_fm_upper_to_full(matrix, work)
given an upper triangular matrix computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1759
cp_fm_basic_linalg::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
...
Definition cp_fm_basic_linalg.F:3162
cp_fm_basic_linalg::cp_fm_symm
subroutine, public cp_fm_symm(side, uplo, m, n, alpha, matrix_a, matrix_b, beta, matrix_c)
computes matrix_c = beta * matrix_c + alpha * matrix_a * matrix_b computes matrix_c = beta * matrix_c...
Definition cp_fm_basic_linalg.F:576
cp_fm_basic_linalg::cp_fm_triangular_multiply
subroutine, public cp_fm_triangular_multiply(triangular_matrix, matrix_b, side, transpose_tr, invert_tr, uplo_tr, unit_diag_tr, n_rows, n_cols, alpha)
multiplies in place by a triangular matrix: matrix_b = alpha op(triangular_matrix) matrix_b or (if si...
Definition cp_fm_basic_linalg.F:1592
cp_fm_cholesky
various cholesky decomposition related routines
Definition cp_fm_cholesky.F:14
cp_fm_cholesky::cp_fm_cholesky_reduce
subroutine, public cp_fm_cholesky_reduce(matrix, matrixb, itype)
reduce a matrix pencil A,B to normal form B has to be cholesky decomposed with cp_fm_cholesky_decompo...
Definition cp_fm_cholesky.F:192
cp_fm_diag
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
cp_fm_diag::cp_fm_block_jacobi
subroutine, public cp_fm_block_jacobi(matrix, eigenvectors, eigval, thresh, start_sec_block)
...
Definition cp_fm_diag.F:1075
cp_fm_diag::choose_eigv_solver
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:208
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_equivalent
logical function, public cp_fm_struct_equivalent(fmstruct1, fmstruct2)
returns true if the two matrix structures are equivalent, false otherwise.
Definition cp_fm_struct.F:357
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1016
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::cholesky_restore
integer, parameter, public cholesky_restore
Definition input_constants.F:195
input_constants::cholesky_off
integer, parameter, public cholesky_off
Definition input_constants.F:195
input_constants::cholesky_reduce
integer, parameter, public cholesky_reduce
Definition input_constants.F:195
input_constants::cholesky_inverse
integer, parameter, public cholesky_inverse
Definition input_constants.F:195
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
qs_density_mixing_types
module that contains the definitions of the scf types
Definition qs_density_mixing_types.F:14
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:397
qs_scf_methods
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
Definition qs_scf_methods.F:19
qs_scf_methods::eigensolver_simple
subroutine, public eigensolver_simple(matrix_ks, mo_set, work, do_level_shift, level_shift, use_jacobi, jacobi_threshold)
...
Definition qs_scf_methods.F:405
qs_scf_methods::eigensolver_dbcsr
subroutine, public eigensolver_dbcsr(matrix_ks, matrix_ks_fm, mo_set, ortho_dbcsr, ksbuf1, ksbuf2)
...
Definition qs_scf_methods.F:268
qs_scf_methods::scf_env_density_mixing
subroutine, public scf_env_density_mixing(p_mix_new, mixing_store, rho_ao, para_env, iter_delta, iter_count, diis, invert)
perform (if requested) a density mixing
Definition qs_scf_methods.F:97
qs_scf_methods::eigensolver_symm
subroutine, public eigensolver_symm(matrix_ks_fm, mo_set, ortho, work, do_level_shift, level_shift, matrix_u_fm, use_jacobi, jacobi_threshold)
...
Definition qs_scf_methods.F:328
qs_scf_methods::eigensolver
subroutine, public eigensolver(matrix_ks_fm, mo_set, ortho, work, cholesky_method, do_level_shift, level_shift, matrix_u_fm, use_jacobi)
Diagonalise the Kohn-Sham matrix to get a new set of MO eigen- vectors and MO eigenvalues....
Definition qs_scf_methods.F:168
qs_scf_methods::cp_sm_mix
subroutine, public cp_sm_mix(m1, m2, p_mix, delta, para_env, m3)
Perform a mixing of the given matrixes into the first matrix m1 = m2 + p_mix (m1-m2)
Definition qs_scf_methods.F:478
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
qs_density_mixing_types::mixing_storage_type
Definition qs_density_mixing_types.F:55
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46