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