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