(git:ccc2433)
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 ! **************************************************************************************************
19 MODULE 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 
78 CONTAINS
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 
838 END MODULE qs_scf_methods
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_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_cholesky::cp_fm_cholesky_restore
subroutine, public cp_fm_cholesky_restore(matrix, neig, matrixb, matrixout, op, pos, transa)
...
Definition: cp_fm_cholesky.F:254
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