Line data Source code
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_api, ONLY: &
22 : dbcsr_add, dbcsr_desymmetrize, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
23 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
24 : dbcsr_multiply, dbcsr_p_type, dbcsr_type
25 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
26 : cp_dbcsr_sm_fm_multiply
27 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
28 : cp_fm_symm,&
29 : cp_fm_triangular_multiply,&
30 : cp_fm_upper_to_full
31 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_reduce,&
32 : cp_fm_cholesky_restore
33 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
34 : cp_fm_block_jacobi
35 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
36 : cp_fm_struct_equivalent,&
37 : cp_fm_struct_release,&
38 : cp_fm_struct_type
39 : USE cp_fm_types, ONLY: cp_fm_create,&
40 : cp_fm_get_info,&
41 : cp_fm_release,&
42 : cp_fm_to_fm,&
43 : cp_fm_type
44 : USE input_constants, ONLY: cholesky_inverse,&
45 : cholesky_off,&
46 : cholesky_reduce,&
47 : cholesky_restore
48 : USE kinds, ONLY: dp
49 : USE message_passing, ONLY: mp_para_env_type
50 : USE parallel_gemm_api, ONLY: parallel_gemm
51 : USE qs_density_mixing_types, ONLY: mixing_storage_type
52 : USE qs_mo_types, ONLY: get_mo_set,&
53 : mo_set_type
54 : #include "./base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 :
58 : PRIVATE
59 :
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_methods'
61 : REAL(KIND=dp), PARAMETER :: ratio = 0.25_dp
62 :
63 : PUBLIC :: combine_ks_matrices, &
64 : cp_sm_mix, &
65 : eigensolver, &
66 : eigensolver_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 96213 : 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 96213 : CALL timeset(routineN, handle)
112 :
113 96213 : my_diis = .FALSE.
114 96213 : IF (PRESENT(diis)) my_diis = diis
115 96213 : my_invert = .FALSE.
116 96213 : IF (PRESENT(invert)) my_invert = invert
117 96213 : my_p_mix = mixing_store%alpha
118 96213 : IF (my_diis .OR. iter_count < mixing_store%nskip_mixing) THEN
119 61197 : my_p_mix = 1.0_dp
120 : END IF
121 :
122 96213 : iter_delta = 0.0_dp
123 96213 : CPASSERT(ASSOCIATED(p_mix_new))
124 496364 : DO ic = 1, SIZE(p_mix_new, 2)
125 968961 : DO ispin = 1, SIZE(p_mix_new, 1)
126 872748 : IF (my_invert) THEN
127 65984 : CPASSERT(my_p_mix /= 0.0_dp)
128 65984 : 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 15490 : 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 406613 : para_env=para_env)
140 406613 : iter_delta = MAX(iter_delta, tmp)
141 : END IF
142 : END DO
143 : END DO
144 :
145 96213 : CALL timestop(handle)
146 :
147 96213 : 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 147458 : 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 73729 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
181 : TYPE(cp_fm_type), POINTER :: mo_coeff
182 :
183 73729 : CALL timeset(routineN, handle)
184 :
185 73729 : NULLIFY (mo_coeff)
186 73729 : 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 73729 : mo_coeff=mo_coeff)
196 :
197 73799 : SELECT CASE (cholesky_method)
198 : CASE (cholesky_reduce)
199 70 : CALL cp_fm_cholesky_reduce(matrix_ks_fm, ortho)
200 :
201 70 : 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 56 : level_shift=level_shift, is_triangular=.TRUE., matrix_u_fm=matrix_u_fm)
204 :
205 70 : CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
206 70 : CALL cp_fm_cholesky_restore(work, nmo, ortho, mo_coeff, "SOLVE")
207 70 : IF (do_level_shift) &
208 56 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
209 :
210 : CASE (cholesky_restore)
211 72647 : CALL cp_fm_upper_to_full(matrix_ks_fm, work)
212 : CALL cp_fm_cholesky_restore(matrix_ks_fm, nao, ortho, work, &
213 72647 : "SOLVE", pos="RIGHT")
214 : CALL cp_fm_cholesky_restore(work, nao, ortho, matrix_ks_fm, &
215 72647 : "SOLVE", pos="LEFT", transa="T")
216 :
217 72647 : 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 100 : level_shift=level_shift, is_triangular=.TRUE., matrix_u_fm=matrix_u_fm)
220 :
221 72647 : CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
222 72647 : CALL cp_fm_cholesky_restore(work, nmo, ortho, mo_coeff, "SOLVE")
223 :
224 72647 : IF (do_level_shift) &
225 100 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
226 :
227 : CASE (cholesky_inverse)
228 1012 : 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 1012 : 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 1012 : invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nao, alpha=1.0_dp)
234 :
235 1012 : 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 56 : level_shift=level_shift, is_triangular=.TRUE., matrix_u_fm=matrix_u_fm)
238 :
239 1012 : CALL choose_eigv_solver(matrix_ks_fm, work, mo_eigenvalues)
240 : CALL cp_fm_triangular_multiply(ortho, work, side="L", transpose_tr=.FALSE., &
241 1012 : invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nmo, alpha=1.0_dp)
242 1012 : CALL cp_fm_to_fm(work, mo_coeff, nmo, 1, 1)
243 :
244 1012 : IF (do_level_shift) &
245 73785 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
246 :
247 : END SELECT
248 :
249 73729 : IF (use_jacobi) THEN
250 0 : CALL cp_fm_to_fm(mo_coeff, ortho)
251 0 : cholesky_method = cholesky_off
252 : END IF
253 :
254 73729 : CALL timestop(handle)
255 :
256 73729 : 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 8472 : 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 2118 : 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 2118 : CALL timeset(routineN, handle)
281 :
282 2118 : NULLIFY (mo_coeff)
283 2118 : NULLIFY (mo_eigenvalues)
284 :
285 : CALL get_mo_set(mo_set=mo_set, &
286 : nao=nao, &
287 : nmo=nmo, &
288 : eigenvalues=mo_eigenvalues, &
289 2118 : mo_coeff=mo_coeff)
290 :
291 : ! Reduce KS matrix
292 2118 : CALL dbcsr_desymmetrize(matrix_ks, ksbuf2)
293 2118 : CALL dbcsr_multiply('N', 'N', 1.0_dp, ksbuf2, ortho_dbcsr, 0.0_dp, ksbuf1)
294 2118 : CALL dbcsr_multiply('T', 'N', 1.0_dp, ortho_dbcsr, ksbuf1, 0.0_dp, ksbuf2)
295 :
296 : ! Solve the eigenvalue problem
297 2118 : CALL copy_dbcsr_to_fm(ksbuf2, matrix_ks_fm)
298 2118 : CALL cp_fm_create(all_evecs, matrix_ks_fm%matrix_struct)
299 2118 : CALL choose_eigv_solver(matrix_ks_fm, all_evecs, mo_eigenvalues)
300 :
301 : ! select first nmo eigenvectors
302 2118 : CALL cp_fm_create(nmo_evecs, mo_coeff%matrix_struct)
303 2118 : CALL cp_fm_to_fm(msource=all_evecs, mtarget=nmo_evecs, ncol=nmo)
304 2118 : CALL cp_fm_release(all_evecs)
305 :
306 : ! Restore the eigenvector of the general eig. problem
307 2118 : CALL cp_dbcsr_sm_fm_multiply(ortho_dbcsr, nmo_evecs, mo_coeff, nmo)
308 :
309 2118 : CALL cp_fm_release(nmo_evecs)
310 2118 : CALL timestop(handle)
311 :
312 2118 : 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 412 : 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 206 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
341 : TYPE(cp_fm_type), POINTER :: mo_coeff
342 :
343 206 : CALL timeset(routineN, handle)
344 :
345 206 : NULLIFY (mo_coeff)
346 206 : 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 206 : mo_coeff=mo_coeff)
357 :
358 206 : IF (use_jacobi) THEN
359 :
360 0 : 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.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 0 : jacobi_threshold, homo + 1)
367 :
368 : ELSE ! full S^(-1/2) has been computed
369 :
370 206 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, matrix_ks_fm, ortho, 0.0_dp, work)
371 206 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, ortho, work, 0.0_dp, matrix_ks_fm)
372 :
373 206 : 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 56 : level_shift=level_shift, is_triangular=.FALSE., matrix_u_fm=matrix_u_fm)
376 :
377 206 : 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 206 : mo_coeff)
381 :
382 206 : IF (do_level_shift) &
383 56 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
384 :
385 : END IF
386 :
387 206 : CALL timestop(handle)
388 :
389 206 : 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 37356 : 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 18678 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
418 : TYPE(cp_fm_type), POINTER :: mo_coeff
419 :
420 18678 : CALL timeset(routineN, handle)
421 :
422 18678 : NULLIFY (mo_coeff)
423 18678 : 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 18678 : mo_coeff=mo_coeff)
432 :
433 18678 : 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 0 : level_shift=level_shift, is_triangular=.FALSE.)
437 : END IF
438 :
439 18678 : IF (use_jacobi) THEN
440 18 : 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 18 : 0.0_dp, matrix_ks, b_first_col=homo + 1, c_first_col=homo + 1)
443 : ! Block Jacobi (pseudo-diagonalization, only one sweep)
444 18 : CALL cp_fm_block_jacobi(matrix_ks, mo_coeff, mo_eigenvalues, jacobi_threshold, homo + 1)
445 : ELSE
446 :
447 18660 : CALL choose_eigv_solver(matrix_ks, work, mo_eigenvalues)
448 :
449 18660 : CALL cp_fm_to_fm(work, mo_coeff, nmo, 1, 1)
450 :
451 : END IF
452 :
453 18678 : IF (do_level_shift) &
454 0 : CALL correct_mo_eigenvalues(mo_eigenvalues, homo, nmo, level_shift)
455 :
456 18678 : CALL timestop(handle)
457 :
458 18678 : 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 1097298 : 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 548649 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_delta_block, p_new_block, p_old_block
490 : TYPE(dbcsr_iterator_type) :: iter
491 :
492 548649 : CALL timeset(routineN, handle)
493 548649 : delta = 0.0_dp
494 :
495 548649 : CALL dbcsr_iterator_start(iter, m1)
496 11924131 : DO WHILE (dbcsr_iterator_blocks_left(iter))
497 11375482 : 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 11375482 : BLOCK=p_old_block, found=found)
500 11375482 : CPASSERT(ASSOCIATED(p_old_block))
501 11924131 : IF (PRESENT(m3)) THEN
502 : CALL dbcsr_get_block_p(matrix=m3, row=iblock_row, col=iblock_col, &
503 1950282 : BLOCK=p_delta_block, found=found)
504 1950282 : CPASSERT(ASSOCIATED(p_delta_block))
505 :
506 18953598 : DO j = 1, SIZE(p_new_block, 2)
507 174749482 : DO i = 1, SIZE(p_new_block, 1)
508 155795884 : p_delta_block(i, j) = p_new_block(i, j) - p_old_block(i, j)
509 172799200 : delta = MAX(delta, ABS(p_delta_block(i, j)))
510 : END DO
511 : END DO
512 : ELSE
513 41884705 : DO j = 1, SIZE(p_new_block, 2)
514 263112058 : DO i = 1, SIZE(p_new_block, 1)
515 221227353 : p_new_block(i, j) = p_new_block(i, j) - p_old_block(i, j)
516 221227353 : delta = MAX(delta, ABS(p_new_block(i, j)))
517 253686858 : 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 548649 : CALL dbcsr_iterator_stop(iter)
523 :
524 548649 : CALL para_env%max(delta)
525 :
526 548649 : CALL timestop(handle)
527 :
528 548649 : 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 940 : 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 940 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
560 : LOGICAL :: compatible_matrices
561 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
562 940 : POINTER :: fa, fb
563 : TYPE(cp_fm_struct_type), POINTER :: ksa_struct, ksb_struct
564 :
565 : ! -------------------------------------------------------------------------
566 :
567 940 : 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 940 : local_data=fa)
576 :
577 : CALL cp_fm_get_info(matrix=ksb, &
578 : matrix_struct=ksb_struct, &
579 940 : local_data=fb)
580 :
581 940 : compatible_matrices = cp_fm_struct_equivalent(ksa_struct, ksb_struct)
582 940 : CPASSERT(compatible_matrices)
583 :
584 40865 : IF (SUM(occb) == 0.0_dp) fb = 0.0_dp
585 :
586 15168 : DO icol_local = 1, ncol_local
587 14228 : icol_global = col_indices(icol_local)
588 14228 : j = INT(occa(icol_global)) + INT(occb(icol_global))
589 169604 : DO irow_local = 1, nrow_local
590 154436 : irow_global = row_indices(irow_local)
591 154436 : 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 168664 : roks_parameter(i, j, 2)*fb(irow_local, icol_local)
595 : END DO
596 : END DO
597 :
598 940 : CALL timestop(handle)
599 :
600 940 : 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 0 : 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 0 : 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 0 : POINTER :: fa, fb
638 : TYPE(cp_fm_struct_type), POINTER :: ksa_struct, ksb_struct
639 :
640 : ! -------------------------------------------------------------------------
641 :
642 0 : 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 0 : local_data=fa)
651 :
652 : CALL cp_fm_get_info(matrix=ksb, &
653 : matrix_struct=ksb_struct, &
654 0 : local_data=fb)
655 :
656 0 : compatible_matrices = cp_fm_struct_equivalent(ksa_struct, ksb_struct)
657 0 : CPASSERT(compatible_matrices)
658 :
659 0 : beta = 1.0_dp/(1.0_dp - f)
660 :
661 0 : DO icol_local = 1, ncol_local
662 :
663 0 : icol_global = col_indices(icol_local)
664 :
665 0 : DO irow_local = 1, nrow_local
666 :
667 0 : irow_global = row_indices(irow_local)
668 :
669 0 : t1 = 0.5_dp*(fa(irow_local, icol_local) + fb(irow_local, icol_local))
670 :
671 0 : IF ((0 < irow_global) .AND. (irow_global <= nbeta)) THEN
672 0 : IF ((0 < icol_global) .AND. (icol_global <= nbeta)) THEN
673 : ! closed-closed
674 0 : fa(irow_local, icol_local) = t1
675 0 : ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
676 : ! closed-open
677 0 : ta = 0.5_dp*(f - REAL(occa(icol_global), KIND=dp))/f
678 0 : tb = 0.5_dp*(f - REAL(occb(icol_global), KIND=dp))/f
679 0 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
680 0 : fa(irow_local, icol_local) = t1 + (beta - 1.0_dp)*t2
681 : ELSE
682 : ! closed-virtual
683 0 : fa(irow_local, icol_local) = t1
684 : END IF
685 0 : 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 0 : ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
693 : ! open-open
694 0 : ta = 0.5_dp*(f - REAL(occa(icol_global), KIND=dp))/f
695 0 : tb = 0.5_dp*(f - REAL(occb(icol_global), KIND=dp))/f
696 0 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
697 0 : IF (irow_global == icol_global) THEN
698 0 : fa(irow_local, icol_local) = t1 - t2
699 : ELSE
700 0 : fa(irow_local, icol_local) = t1 - 0.5_dp*t2
701 : END IF
702 : ELSE
703 : ! open-virtual
704 0 : ta = 0.5_dp*(f - REAL(occa(irow_global), KIND=dp))/f
705 0 : tb = 0.5_dp*(f - REAL(occb(irow_global), KIND=dp))/f
706 0 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
707 0 : fa(irow_local, icol_local) = t1 - t2
708 : END IF
709 : ELSE
710 0 : IF ((0 < irow_global) .AND. (irow_global < nbeta)) THEN
711 : ! virtual-closed
712 0 : fa(irow_local, icol_local) = t1
713 0 : ELSE IF ((nbeta < icol_global) .AND. (icol_global <= nalpha)) THEN
714 : ! virtual-open
715 0 : ta = 0.5_dp*(f - REAL(occa(icol_global), KIND=dp))/f
716 0 : tb = 0.5_dp*(f - REAL(occb(icol_global), KIND=dp))/f
717 0 : t2 = ta*fa(irow_local, icol_local) + tb*fb(irow_local, icol_local)
718 0 : fa(irow_local, icol_local) = t1 - t2
719 : ELSE
720 : ! virtual-virtual
721 0 : fa(irow_local, icol_local) = t1
722 : END IF
723 : END IF
724 :
725 : END DO
726 : END DO
727 :
728 0 : CALL timestop(handle)
729 :
730 0 : 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 268 : 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 7586 : DO imo = homo + 1, nmo
754 7586 : mo_eigenvalues(imo) = mo_eigenvalues(imo) - level_shift
755 : END DO
756 :
757 268 : 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 268 : 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 268 : 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 268 : CALL timeset(routineN, handle)
791 :
792 268 : NULLIFY (ao_mo_fmstruct)
793 : CALL cp_fm_struct_create(ao_mo_fmstruct, nrow_global=nao, ncol_global=nmo, &
794 268 : para_env=mo_coeff%matrix_struct%para_env, context=mo_coeff%matrix_struct%context)
795 :
796 268 : CALL cp_fm_create(u_mo, ao_mo_fmstruct)
797 268 : CALL cp_fm_create(u_mo_scaled, ao_mo_fmstruct)
798 :
799 268 : CALL cp_fm_struct_release(ao_mo_fmstruct)
800 :
801 : ! U * C
802 268 : IF (PRESENT(matrix_u_fm)) THEN
803 268 : IF (is_triangular) THEN
804 212 : 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 212 : invert_tr=.FALSE., uplo_tr="U", n_rows=nao, n_cols=nmo, alpha=1.0_dp)
807 : ELSE
808 56 : 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 0 : CALL cp_fm_to_fm(mo_coeff, u_mo)
813 : END IF
814 :
815 268 : 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 804 : ALLOCATE (weights(nmo))
821 1454 : weights(1:homo) = 0.0_dp
822 7586 : 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 268 : CALL cp_fm_column_scale(u_mo_scaled, weights)
826 268 : DEALLOCATE (weights)
827 :
828 : ! NewKS = U^{-1,T} * KS * U^{-1} + (U * C) * DELTA * (U * C)^T
829 268 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, u_mo, u_mo_scaled, 1.0_dp, matrix_ks_fm)
830 :
831 268 : CALL cp_fm_release(u_mo_scaled)
832 268 : CALL cp_fm_release(u_mo)
833 :
834 268 : CALL timestop(handle)
835 :
836 536 : END SUBROUTINE shift_unocc_mos
837 :
838 : END MODULE qs_scf_methods
|
