(git:e5b1968)
Loading...
Searching...
No Matches
qs_scf_diagonalization.F
Go to the documentation of this file.
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 Different diagonalization schemes that can be used
10!> for the iterative solution of the eigenvalue problem
11!> \par History
12!> started from routines previously located in the qs_scf module
13!> 05.2009
14! **************************************************************************************************
20 USE cp_cfm_diag, ONLY: cp_cfm_geeig,&
22 USE cp_cfm_types, ONLY: cp_cfm_create,&
29 USE cp_dbcsr_api, ONLY: &
31 dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, &
32 dbcsr_type_symmetric
51 USE cp_fm_types, ONLY: &
59 USE input_constants, ONLY: &
64 USE kinds, ONLY: dp
69 USE kpoint_types, ONLY: get_kpoint_info,&
72 USE machine, ONLY: m_flush,&
90 USE qs_ks_types, ONLY: qs_ks_did_change,&
92 USE qs_matrix_pools, ONLY: mpools_get,&
100 USE qs_mo_types, ONLY: get_mo_set,&
106 USE qs_rho_types, ONLY: qs_rho_get,&
118 USE qs_scf_types, ONLY: qs_scf_env_type,&
121#include "./base/base_uses.f90"
122
123 IMPLICIT NONE
124
125 PRIVATE
126
127 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_diagonalization'
128
132
133CONTAINS
134
135! **************************************************************************************************
136!> \brief the inner loop of scf, specific to diagonalization with S matrix
137!> basically, in goes the ks matrix out goes a new p matrix
138!> \param scf_env ...
139!> \param mos ...
140!> \param matrix_ks ...
141!> \param matrix_s ...
142!> \param scf_control ...
143!> \param scf_section ...
144!> \param diis_step ...
145!> \par History
146!> 03.2006 created [Joost VandeVondele]
147! **************************************************************************************************
148
149 SUBROUTINE general_eigenproblem(scf_env, mos, matrix_ks, &
150 matrix_s, scf_control, scf_section, &
151 diis_step)
152
153 TYPE(qs_scf_env_type), POINTER :: scf_env
154 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
155 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
156 TYPE(scf_control_type), POINTER :: scf_control
157 TYPE(section_vals_type), POINTER :: scf_section
158 LOGICAL, INTENT(INOUT) :: diis_step
159
160 INTEGER :: ispin, nspin
161 LOGICAL :: do_level_shift, owns_ortho, use_jacobi
162 REAL(kind=dp) :: diis_error, eps_diis
163 TYPE(cp_fm_type), POINTER :: ortho
164 TYPE(dbcsr_type), POINTER :: ortho_dbcsr
165
166 nspin = SIZE(matrix_ks)
167 NULLIFY (ortho, ortho_dbcsr)
168
169 DO ispin = 1, nspin
170 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
171 scf_env%scf_work1(ispin))
172 END DO
173
174 eps_diis = scf_control%eps_diis
175
176 IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis) THEN
177 CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
178 scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
179 eps_diis, scf_control%nmixing, &
180 s_matrix=matrix_s, &
181 scf_section=scf_section)
182 ELSE
183 diis_step = .false.
184 END IF
185
186 do_level_shift = ((scf_control%level_shift /= 0.0_dp) .AND. &
187 ((scf_control%density_guess == core_guess) .OR. &
188 (scf_env%iter_count > 1)))
189
190 IF ((scf_env%iter_count > 1) .AND. &
191 (scf_env%iter_delta < scf_control%diagonalization%eps_jacobi)) THEN
192 use_jacobi = .true.
193 ELSE
194 use_jacobi = .false.
195 END IF
196
197 IF (diis_step) THEN
198 scf_env%iter_param = diis_error
199 IF (use_jacobi) THEN
200 scf_env%iter_method = "DIIS/Jacobi"
201 ELSE
202 scf_env%iter_method = "DIIS/Diag."
203 END IF
204 ELSE
205 IF (scf_env%mixing_method == 0) THEN
206 scf_env%iter_method = "NoMix/Diag."
207 ELSE IF (scf_env%mixing_method == 1) THEN
208 scf_env%iter_param = scf_env%p_mix_alpha
209 IF (use_jacobi) THEN
210 scf_env%iter_method = "P_Mix/Jacobi"
211 ELSE
212 scf_env%iter_method = "P_Mix/Diag."
213 END IF
214 ELSEIF (scf_env%mixing_method > 1) THEN
215 scf_env%iter_param = scf_env%mixing_store%alpha
216 IF (use_jacobi) THEN
217 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Jacobi"
218 ELSE
219 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Diag."
220 END IF
221 END IF
222 END IF
223
224 IF (scf_env%cholesky_method == cholesky_dbcsr) THEN
225 ortho_dbcsr => scf_env%ortho_dbcsr
226 DO ispin = 1, nspin
227 CALL eigensolver_dbcsr(matrix_ks=matrix_ks(ispin)%matrix, matrix_ks_fm=scf_env%scf_work1(ispin), &
228 mo_set=mos(ispin), &
229 ortho_dbcsr=ortho_dbcsr, &
230 ksbuf1=scf_env%buf1_dbcsr, ksbuf2=scf_env%buf2_dbcsr)
231 END DO
232
233 ELSE IF (scf_env%cholesky_method > cholesky_off) THEN
234 IF (scf_env%cholesky_method == cholesky_inverse) THEN
235 ortho => scf_env%ortho_m1
236 ELSE
237 ortho => scf_env%ortho
238 END IF
239
240 owns_ortho = .false.
241 IF (.NOT. ASSOCIATED(ortho)) THEN
242 ALLOCATE (ortho)
243 owns_ortho = .true.
244 END IF
245
246 DO ispin = 1, nspin
247 IF (do_level_shift) THEN
248 CALL eigensolver(matrix_ks_fm=scf_env%scf_work1(ispin), &
249 mo_set=mos(ispin), &
250 ortho=ortho, &
251 work=scf_env%scf_work2, &
252 cholesky_method=scf_env%cholesky_method, &
253 do_level_shift=do_level_shift, &
254 level_shift=scf_control%level_shift, &
255 matrix_u_fm=scf_env%ortho, &
256 use_jacobi=use_jacobi)
257 ELSE
258 CALL eigensolver(matrix_ks_fm=scf_env%scf_work1(ispin), &
259 mo_set=mos(ispin), &
260 ortho=ortho, &
261 work=scf_env%scf_work2, &
262 cholesky_method=scf_env%cholesky_method, &
263 do_level_shift=do_level_shift, &
264 level_shift=scf_control%level_shift, &
265 use_jacobi=use_jacobi)
266 END IF
267 END DO
268
269 IF (owns_ortho) DEALLOCATE (ortho)
270 ELSE
271 ortho => scf_env%ortho
272
273 owns_ortho = .false.
274 IF (.NOT. ASSOCIATED(ortho)) THEN
275 ALLOCATE (ortho)
276 owns_ortho = .true.
277 END IF
278
279 IF (do_level_shift) THEN
280 DO ispin = 1, nspin
281 IF (ASSOCIATED(scf_env%scf_work1_red) .AND. ASSOCIATED(scf_env%scf_work2_red) &
282 .AND. ASSOCIATED(scf_env%ortho_red) .AND. ASSOCIATED(scf_env%ortho_m1_red)) THEN
283 CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
284 mo_set=mos(ispin), &
285 ortho=ortho, &
286 work=scf_env%scf_work2, &
287 do_level_shift=do_level_shift, &
288 level_shift=scf_control%level_shift, &
289 matrix_u_fm=scf_env%ortho_m1, &
290 use_jacobi=use_jacobi, &
291 jacobi_threshold=scf_control%diagonalization%jacobi_threshold, &
292 matrix_ks_fm_red=scf_env%scf_work1_red(ispin), &
293 ortho_red=scf_env%ortho_red, &
294 work_red=scf_env%scf_work2_red, &
295 matrix_u_fm_red=scf_env%ortho_m1_red)
296 ELSE
297 CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
298 mo_set=mos(ispin), &
299 ortho=ortho, &
300 work=scf_env%scf_work2, &
301 do_level_shift=do_level_shift, &
302 level_shift=scf_control%level_shift, &
303 matrix_u_fm=scf_env%ortho_m1, &
304 use_jacobi=use_jacobi, &
305 jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
306 END IF
307 END DO
308 ELSE
309 DO ispin = 1, nspin
310 IF (ASSOCIATED(scf_env%scf_work1_red) .AND. ASSOCIATED(scf_env%scf_work2_red) &
311 .AND. ASSOCIATED(scf_env%ortho_red)) THEN
312 CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
313 mo_set=mos(ispin), &
314 ortho=ortho, &
315 work=scf_env%scf_work2, &
316 do_level_shift=do_level_shift, &
317 level_shift=scf_control%level_shift, &
318 use_jacobi=use_jacobi, &
319 jacobi_threshold=scf_control%diagonalization%jacobi_threshold, &
320 matrix_ks_fm_red=scf_env%scf_work1_red(ispin), &
321 ortho_red=scf_env%ortho_red, &
322 work_red=scf_env%scf_work2_red)
323 ELSE
324 CALL eigensolver_symm(matrix_ks_fm=scf_env%scf_work1(ispin), &
325 mo_set=mos(ispin), &
326 ortho=ortho, &
327 work=scf_env%scf_work2, &
328 do_level_shift=do_level_shift, &
329 level_shift=scf_control%level_shift, &
330 use_jacobi=use_jacobi, &
331 jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
332 END IF
333 END DO
334 END IF
335
336 IF (owns_ortho) DEALLOCATE (ortho)
337 END IF
338
339 END SUBROUTINE general_eigenproblem
340
341! **************************************************************************************************
342!> \brief ...
343!> \param scf_env ...
344!> \param mos ...
345!> \param matrix_ks ...
346!> \param matrix_s ...
347!> \param scf_control ...
348!> \param scf_section ...
349!> \param diis_step ...
350!> \param probe ...
351! **************************************************************************************************
352 SUBROUTINE do_general_diag(scf_env, mos, matrix_ks, &
353 matrix_s, scf_control, scf_section, &
354 diis_step, probe)
355
356 TYPE(qs_scf_env_type), POINTER :: scf_env
357 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
358 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
359 TYPE(scf_control_type), POINTER :: scf_control
360 TYPE(section_vals_type), POINTER :: scf_section
361 LOGICAL, INTENT(INOUT) :: diis_step
362 TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
363 POINTER :: probe
364
365 INTEGER :: ispin, nspin
366 REAL(kind=dp) :: total_zeff_corr
367
368 nspin = SIZE(matrix_ks)
369
370 CALL general_eigenproblem(scf_env, mos, matrix_ks, &
371 matrix_s, scf_control, scf_section, diis_step)
372
373 total_zeff_corr = 0.0_dp
374 total_zeff_corr = scf_env%sum_zeff_corr
375
376 IF (abs(total_zeff_corr) > 0.0_dp) THEN
377 CALL set_mo_occupation(mo_array=mos, &
378 smear=scf_control%smear, tot_zeff_corr=total_zeff_corr)
379 ELSE
380 IF (PRESENT(probe) .EQV. .true.) THEN
381 scf_control%smear%do_smear = .false.
382 CALL set_mo_occupation(mo_array=mos, &
383 smear=scf_control%smear, &
384 probe=probe)
385 ELSE
386 CALL set_mo_occupation(mo_array=mos, &
387 smear=scf_control%smear)
388 END IF
389 END IF
390
391 DO ispin = 1, nspin
392 CALL calculate_density_matrix(mos(ispin), &
393 scf_env%p_mix_new(ispin, 1)%matrix)
394 END DO
395
396 END SUBROUTINE do_general_diag
397
398! **************************************************************************************************
399!> \brief Kpoint diagonalization routine
400!> Transforms matrices to kpoint, distributes kpoint groups, performs
401!> general diagonalization (no storgae of overlap decomposition), stores
402!> MOs, calculates occupation numbers, calculates density matrices
403!> in kpoint representation, transforms density matrices to real space
404!> \param matrix_ks Kohn-sham matrices (RS indices, global)
405!> \param matrix_s Overlap matrices (RS indices, global)
406!> \param kpoints Kpoint environment
407!> \param scf_env SCF environment
408!> \param scf_control SCF control variables
409!> \param update_p ...
410!> \param diis_step ...
411!> \param diis_error ...
412!> \param qs_env ...
413!> \param probe ...
414!> \par History
415!> 08.2014 created [JGH]
416! **************************************************************************************************
417 SUBROUTINE do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, update_p, &
418 diis_step, diis_error, qs_env, probe)
419
420 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s
421 TYPE(kpoint_type), POINTER :: kpoints
422 TYPE(qs_scf_env_type), POINTER :: scf_env
423 TYPE(scf_control_type), POINTER :: scf_control
424 LOGICAL, INTENT(IN) :: update_p
425 LOGICAL, INTENT(INOUT) :: diis_step
426 REAL(dp), INTENT(INOUT), OPTIONAL :: diis_error
427 TYPE(qs_environment_type), OPTIONAL, POINTER :: qs_env
428 TYPE(hairy_probes_type), DIMENSION(:), OPTIONAL, &
429 POINTER :: probe
430
431 CHARACTER(len=*), PARAMETER :: routinen = 'do_general_diag_kp'
432 COMPLEX(KIND=dp), PARAMETER :: cone = cmplx(1.0_dp, 0.0_dp, kind=dp), &
433 czero = cmplx(0.0_dp, 0.0_dp, kind=dp), ione = cmplx(0.0_dp, 1.0_dp, kind=dp)
434
435 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coeffs
436 INTEGER :: handle, ib, igroup, ik, ikp, indx, &
437 ispin, jb, kplocal, nb, nkp, &
438 nkp_groups, nspin
439 INTEGER, DIMENSION(2) :: kp_range
440 INTEGER, DIMENSION(:, :), POINTER :: kp_dist
441 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
442 LOGICAL :: do_diis, my_kpgrp, use_real_wfn
443 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
444 REAL(kind=dp), DIMENSION(:, :), POINTER :: xkp
445 TYPE(copy_info_type), ALLOCATABLE, DIMENSION(:, :) :: info
446 TYPE(cp_cfm_type) :: cksmat, cmos, csmat, cwork
447 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
448 TYPE(cp_fm_struct_type), POINTER :: matrix_struct, mo_struct
449 TYPE(cp_fm_type) :: fmdummy, fmlocal, rksmat, rsmat
450 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fmwork
451 TYPE(cp_fm_type), POINTER :: imos, mo_coeff, rmos
452 TYPE(dbcsr_type), POINTER :: cmatrix, rmatrix, tmpmat
453 TYPE(kpoint_env_type), POINTER :: kp
454 TYPE(mp_para_env_type), POINTER :: para_env, para_env_global
455 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
456 POINTER :: sab_nl
457 TYPE(qs_matrix_pools_type), POINTER :: mpools
458 TYPE(section_vals_type), POINTER :: scf_section
459
460 CALL timeset(routinen, handle)
461
462 NULLIFY (sab_nl)
463 CALL get_kpoint_info(kpoints, nkp=nkp, xkp=xkp, use_real_wfn=use_real_wfn, kp_range=kp_range, &
464 nkp_groups=nkp_groups, kp_dist=kp_dist, sab_nl=sab_nl, &
465 cell_to_index=cell_to_index)
466 cpassert(ASSOCIATED(sab_nl))
467 kplocal = kp_range(2) - kp_range(1) + 1
468
469 !Whether we use DIIS for k-points
470 do_diis = .false.
471 IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis .AND. .NOT. use_real_wfn &
472 .AND. PRESENT(diis_error) .AND. PRESENT(qs_env)) do_diis = .true.
473
474 ! allocate some work matrices
475 ALLOCATE (rmatrix, cmatrix, tmpmat)
476 CALL dbcsr_create(rmatrix, template=matrix_ks(1, 1)%matrix, &
477 matrix_type=dbcsr_type_symmetric)
478 CALL dbcsr_create(cmatrix, template=matrix_ks(1, 1)%matrix, &
479 matrix_type=dbcsr_type_antisymmetric)
480 CALL dbcsr_create(tmpmat, template=matrix_ks(1, 1)%matrix, &
481 matrix_type=dbcsr_type_no_symmetry)
482 CALL cp_dbcsr_alloc_block_from_nbl(rmatrix, sab_nl)
483 CALL cp_dbcsr_alloc_block_from_nbl(cmatrix, sab_nl)
484
485 fmwork => scf_env%scf_work1
486
487 ! fm pools to be used within a kpoint group
488 CALL get_kpoint_info(kpoints, mpools=mpools)
489 CALL mpools_get(mpools, ao_ao_fm_pools=ao_ao_fm_pools)
490
491 CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, fmlocal)
492 CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
493
494 IF (use_real_wfn) THEN
495 CALL cp_fm_create(rksmat, matrix_struct)
496 CALL cp_fm_create(rsmat, matrix_struct)
497 ELSE
498 CALL cp_cfm_create(cksmat, matrix_struct)
499 CALL cp_cfm_create(csmat, matrix_struct)
500 CALL cp_cfm_create(cwork, matrix_struct)
501 kp => kpoints%kp_env(1)%kpoint_env
502 CALL get_mo_set(kp%mos(1, 1), mo_coeff=mo_coeff)
503 CALL cp_fm_get_info(mo_coeff, matrix_struct=mo_struct)
504 CALL cp_cfm_create(cmos, mo_struct)
505 END IF
506
507 para_env => kpoints%blacs_env_all%para_env
508 nspin = SIZE(matrix_ks, 1)
509 ALLOCATE (info(kplocal*nspin*nkp_groups, 4))
510
511 ! Setup and start all the communication
512 indx = 0
513 DO ikp = 1, kplocal
514 DO ispin = 1, nspin
515 DO igroup = 1, nkp_groups
516 ! number of current kpoint
517 ik = kp_dist(1, igroup) + ikp - 1
518 my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
519 indx = indx + 1
520 IF (use_real_wfn) THEN
521 ! FT of matrices KS and S, then transfer to FM type
522 CALL dbcsr_set(rmatrix, 0.0_dp)
523 CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_ks, ispin=ispin, &
524 xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
525 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
526 CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
527 ! s matrix is not spin dependent
528 CALL dbcsr_set(rmatrix, 0.0_dp)
529 CALL rskp_transform(rmatrix=rmatrix, rsmat=matrix_s, ispin=1, &
530 xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
531 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
532 CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
533 ELSE
534 ! FT of matrices KS and S, then transfer to FM type
535 CALL dbcsr_set(rmatrix, 0.0_dp)
536 CALL dbcsr_set(cmatrix, 0.0_dp)
537 CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_ks, ispin=ispin, &
538 xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
539 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
540 CALL copy_dbcsr_to_fm(tmpmat, fmwork(1))
541 CALL dbcsr_desymmetrize(cmatrix, tmpmat)
542 CALL copy_dbcsr_to_fm(tmpmat, fmwork(2))
543 ! s matrix is not spin dependent, double the work
544 CALL dbcsr_set(rmatrix, 0.0_dp)
545 CALL dbcsr_set(cmatrix, 0.0_dp)
546 CALL rskp_transform(rmatrix=rmatrix, cmatrix=cmatrix, rsmat=matrix_s, ispin=1, &
547 xkp=xkp(1:3, ik), cell_to_index=cell_to_index, sab_nl=sab_nl)
548 CALL dbcsr_desymmetrize(rmatrix, tmpmat)
549 CALL copy_dbcsr_to_fm(tmpmat, fmwork(3))
550 CALL dbcsr_desymmetrize(cmatrix, tmpmat)
551 CALL copy_dbcsr_to_fm(tmpmat, fmwork(4))
552 END IF
553 ! transfer to kpoint group
554 ! redistribution of matrices, new blacs environment
555 ! fmwork -> fmlocal -> rksmat/cksmat
556 ! fmwork -> fmlocal -> rsmat/csmat
557 IF (use_real_wfn) THEN
558 IF (my_kpgrp) THEN
559 CALL cp_fm_start_copy_general(fmwork(1), rksmat, para_env, info(indx, 1))
560 CALL cp_fm_start_copy_general(fmwork(3), rsmat, para_env, info(indx, 2))
561 ELSE
562 CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
563 CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 2))
564 END IF
565 ELSE
566 IF (my_kpgrp) THEN
567 CALL cp_fm_start_copy_general(fmwork(1), fmlocal, para_env, info(indx, 1))
568 CALL cp_fm_start_copy_general(fmwork(2), fmlocal, para_env, info(indx, 2))
569 CALL cp_fm_start_copy_general(fmwork(3), fmlocal, para_env, info(indx, 3))
570 CALL cp_fm_start_copy_general(fmwork(4), fmlocal, para_env, info(indx, 4))
571 ELSE
572 CALL cp_fm_start_copy_general(fmwork(1), fmdummy, para_env, info(indx, 1))
573 CALL cp_fm_start_copy_general(fmwork(2), fmdummy, para_env, info(indx, 2))
574 CALL cp_fm_start_copy_general(fmwork(3), fmdummy, para_env, info(indx, 3))
575 CALL cp_fm_start_copy_general(fmwork(4), fmdummy, para_env, info(indx, 4))
576 END IF
577 END IF
578 END DO
579 END DO
580 END DO
581
582 ! Finish communication then diagonalise in each group
583 IF (do_diis) THEN
584 CALL get_qs_env(qs_env, para_env=para_env_global)
585 scf_section => section_vals_get_subs_vals(qs_env%input, "DFT%SCF")
586 CALL qs_diis_b_info_kp(kpoints%scf_diis_buffer, ib, nb)
587 indx = 0
588 DO ikp = 1, kplocal
589 DO ispin = 1, nspin
590 DO igroup = 1, nkp_groups
591 ! number of current kpoint
592 ik = kp_dist(1, igroup) + ikp - 1
593 my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
594 indx = indx + 1
595 IF (my_kpgrp) THEN
596 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
597 CALL cp_cfm_scale_and_add_fm(czero, cksmat, cone, fmlocal)
598 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
599 CALL cp_cfm_scale_and_add_fm(cone, cksmat, ione, fmlocal)
600 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 3))
601 CALL cp_cfm_scale_and_add_fm(czero, csmat, cone, fmlocal)
602 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 4))
603 CALL cp_cfm_scale_and_add_fm(cone, csmat, ione, fmlocal)
604 END IF
605 END DO !igroup
606
607 kp => kpoints%kp_env(ikp)%kpoint_env
608 CALL qs_diis_b_calc_err_kp(kpoints%scf_diis_buffer, ib, kp%mos, cksmat, csmat, &
609 ispin, ikp, kplocal, scf_section)
610
611 END DO !ispin
612 END DO !ikp
613
614 ALLOCATE (coeffs(nb))
615 CALL qs_diis_b_step_kp(kpoints%scf_diis_buffer, coeffs, ib, nb, scf_env%iter_delta, diis_error, &
616 diis_step, scf_control%eps_diis, nspin, nkp, kplocal, scf_control%nmixing, &
617 scf_section, para_env_global)
618
619 !build the ks matrices and idagonalize
620 DO ikp = 1, kplocal
621 DO ispin = 1, nspin
622 kp => kpoints%kp_env(ikp)%kpoint_env
623 CALL cp_cfm_to_cfm(kpoints%scf_diis_buffer%smat(ikp), csmat)
624
625 CALL cp_cfm_scale(czero, cksmat)
626 DO jb = 1, nb
627 CALL cp_cfm_scale_and_add(cone, cksmat, coeffs(jb), kpoints%scf_diis_buffer%param(jb, ispin, ikp))
628 END DO
629
630 CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
631 CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
632 IF (scf_env%cholesky_method == cholesky_off) THEN
633 CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, &
634 scf_control%eps_eigval)
635 ELSE
636 CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
637 END IF
638 ! copy eigenvalues to imag set (keep them in sync)
639 kp%mos(2, ispin)%eigenvalues = eigenvalues
640 ! split real and imaginary part of mos
641 CALL cp_cfm_to_fm(cmos, rmos, imos)
642 END DO
643 END DO
644
645 ELSE !no DIIS
646 diis_step = .false.
647 indx = 0
648 DO ikp = 1, kplocal
649 DO ispin = 1, nspin
650 DO igroup = 1, nkp_groups
651 ! number of current kpoint
652 ik = kp_dist(1, igroup) + ikp - 1
653 my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
654 indx = indx + 1
655 IF (my_kpgrp) THEN
656 IF (use_real_wfn) THEN
657 CALL cp_fm_finish_copy_general(rksmat, info(indx, 1))
658 CALL cp_fm_finish_copy_general(rsmat, info(indx, 2))
659 ELSE
660 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 1))
661 CALL cp_cfm_scale_and_add_fm(czero, cksmat, cone, fmlocal)
662 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 2))
663 CALL cp_cfm_scale_and_add_fm(cone, cksmat, ione, fmlocal)
664 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 3))
665 CALL cp_cfm_scale_and_add_fm(czero, csmat, cone, fmlocal)
666 CALL cp_fm_finish_copy_general(fmlocal, info(indx, 4))
667 CALL cp_cfm_scale_and_add_fm(cone, csmat, ione, fmlocal)
668 END IF
669 END IF
670 END DO
671
672 ! Each kpoint group has now information on a kpoint to be diagonalized
673 ! General eigensolver Hermite or Symmetric
674 kp => kpoints%kp_env(ikp)%kpoint_env
675 IF (use_real_wfn) THEN
676 CALL get_mo_set(kp%mos(1, ispin), mo_coeff=mo_coeff, eigenvalues=eigenvalues)
677 IF (scf_env%cholesky_method == cholesky_off) THEN
678 CALL cp_fm_geeig_canon(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal, &
679 scf_control%eps_eigval)
680 ELSE
681 CALL cp_fm_geeig(rksmat, rsmat, mo_coeff, eigenvalues, fmlocal)
682 END IF
683 ELSE
684 CALL get_mo_set(kp%mos(1, ispin), mo_coeff=rmos, eigenvalues=eigenvalues)
685 CALL get_mo_set(kp%mos(2, ispin), mo_coeff=imos)
686 IF (scf_env%cholesky_method == cholesky_off) THEN
687 CALL cp_cfm_geeig_canon(cksmat, csmat, cmos, eigenvalues, cwork, &
688 scf_control%eps_eigval)
689 ELSE
690 CALL cp_cfm_geeig(cksmat, csmat, cmos, eigenvalues, cwork)
691 END IF
692 ! copy eigenvalues to imag set (keep them in sync)
693 kp%mos(2, ispin)%eigenvalues = eigenvalues
694 ! split real and imaginary part of mos
695 CALL cp_cfm_to_fm(cmos, rmos, imos)
696 END IF
697 END DO
698 END DO
699 END IF
700
701 ! Clean up communication
702 indx = 0
703 DO ikp = 1, kplocal
704 DO ispin = 1, nspin
705 DO igroup = 1, nkp_groups
706 ! number of current kpoint
707 ik = kp_dist(1, igroup) + ikp - 1
708 my_kpgrp = (ik >= kpoints%kp_range(1) .AND. ik <= kpoints%kp_range(2))
709 indx = indx + 1
710 IF (use_real_wfn) THEN
711 CALL cp_fm_cleanup_copy_general(info(indx, 1))
712 CALL cp_fm_cleanup_copy_general(info(indx, 2))
713 ELSE
714 CALL cp_fm_cleanup_copy_general(info(indx, 1))
715 CALL cp_fm_cleanup_copy_general(info(indx, 2))
716 CALL cp_fm_cleanup_copy_general(info(indx, 3))
717 CALL cp_fm_cleanup_copy_general(info(indx, 4))
718 END IF
719 END DO
720 END DO
721 END DO
722
723 ! All done
724 DEALLOCATE (info)
725
726 IF (update_p) THEN
727 ! MO occupations
728 IF (PRESENT(probe) .EQV. .true.) THEN
729 scf_control%smear%do_smear = .false.
730 CALL kpoint_set_mo_occupation(kpoints, scf_control%smear, &
731 probe=probe)
732 ELSE
733 CALL kpoint_set_mo_occupation(kpoints, scf_control%smear)
734 END IF
735 ! density matrices
736 CALL kpoint_density_matrices(kpoints)
737 ! density matrices in real space
738 CALL kpoint_density_transform(kpoints, scf_env%p_mix_new, .false., &
739 matrix_s(1, 1)%matrix, sab_nl, fmwork)
740 END IF
741
742 CALL dbcsr_deallocate_matrix(rmatrix)
743 CALL dbcsr_deallocate_matrix(cmatrix)
744 CALL dbcsr_deallocate_matrix(tmpmat)
745
746 IF (use_real_wfn) THEN
747 CALL cp_fm_release(rksmat)
748 CALL cp_fm_release(rsmat)
749 ELSE
750 CALL cp_cfm_release(cksmat)
751 CALL cp_cfm_release(csmat)
752 CALL cp_cfm_release(cwork)
753 CALL cp_cfm_release(cmos)
754 END IF
755 CALL fm_pool_give_back_fm(ao_ao_fm_pools(1)%pool, fmlocal)
756
757 CALL timestop(handle)
758
759 END SUBROUTINE do_general_diag_kp
760
761! **************************************************************************************************
762!> \brief inner loop within MOS subspace, to refine occupation and density,
763!> before next diagonalization of the Hamiltonian
764!> \param qs_env ...
765!> \param scf_env ...
766!> \param subspace_env ...
767!> \param mos ...
768!> \param rho ...
769!> \param ks_env ...
770!> \param scf_section ...
771!> \param scf_control ...
772!> \par History
773!> 09.2009 created [MI]
774!> \note it is assumed that when diagonalization is used, also some mixing procedure is active
775! **************************************************************************************************
776 SUBROUTINE do_scf_diag_subspace(qs_env, scf_env, subspace_env, mos, rho, &
777 ks_env, scf_section, scf_control)
778
779 TYPE(qs_environment_type), POINTER :: qs_env
780 TYPE(qs_scf_env_type), POINTER :: scf_env
781 TYPE(subspace_env_type), POINTER :: subspace_env
782 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
783 TYPE(qs_rho_type), POINTER :: rho
784 TYPE(qs_ks_env_type), POINTER :: ks_env
785 TYPE(section_vals_type), POINTER :: scf_section
786 TYPE(scf_control_type), POINTER :: scf_control
787
788 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_scf_diag_subspace'
789 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
790
791 INTEGER :: handle, i, iloop, ispin, nao, nmo, &
792 nspin, output_unit
793 LOGICAL :: converged
794 REAL(dp) :: ene_diff, ene_old, iter_delta, max_val, &
795 sum_band, sum_val, t1, t2
796 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues, mo_occupations
797 TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: eval_first, occ_first
798 TYPE(cp_fm_type) :: work
799 TYPE(cp_fm_type), POINTER :: c0, chc, evec, mo_coeff
800 TYPE(cp_logger_type), POINTER :: logger
801 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
802 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp
803 TYPE(dft_control_type), POINTER :: dft_control
804 TYPE(mp_para_env_type), POINTER :: para_env
805 TYPE(qs_energy_type), POINTER :: energy
806 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom
807
808 CALL timeset(routinen, handle)
809 NULLIFY (c0, chc, energy, evec, matrix_ks, mo_coeff, mo_eigenvalues, &
810 mo_occupations, dft_control, rho_ao, rho_ao_kp)
811
812 logger => cp_get_default_logger()
813 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DIAG_SUB_SCF", &
814 extension=".scfLog")
815
816 !Extra loop keeping mos unchanged and refining the subspace occupation
817 nspin = SIZE(mos)
818 CALL qs_rho_get(rho, rho_ao=rho_ao, rho_ao_kp=rho_ao_kp)
819
820 ALLOCATE (eval_first(nspin))
821 ALLOCATE (occ_first(nspin))
822 DO ispin = 1, nspin
823 CALL get_mo_set(mo_set=mos(ispin), &
824 nmo=nmo, &
825 eigenvalues=mo_eigenvalues, &
826 occupation_numbers=mo_occupations)
827 ALLOCATE (eval_first(ispin)%array(nmo))
828 ALLOCATE (occ_first(ispin)%array(nmo))
829 eval_first(ispin)%array(1:nmo) = mo_eigenvalues(1:nmo)
830 occ_first(ispin)%array(1:nmo) = mo_occupations(1:nmo)
831 END DO
832
833 DO ispin = 1, nspin
834 ! does not yet handle k-points
835 CALL dbcsr_copy(subspace_env%p_matrix_store(ispin)%matrix, rho_ao(ispin)%matrix)
836 CALL dbcsr_copy(rho_ao(ispin)%matrix, scf_env%p_mix_new(ispin, 1)%matrix)
837 END DO
838
839 subspace_env%p_matrix_mix => scf_env%p_mix_new
840
841 NULLIFY (matrix_ks, energy, para_env, matrix_s)
842 CALL get_qs_env(qs_env, &
843 matrix_ks=matrix_ks, &
844 energy=energy, &
845 matrix_s=matrix_s, &
846 para_env=para_env, &
847 dft_control=dft_control)
848
849 ! mixing storage allocation
850 IF (subspace_env%mixing_method >= gspace_mixing_nr) THEN
851 CALL mixing_allocate(qs_env, subspace_env%mixing_method, scf_env%p_mix_new, &
852 scf_env%p_delta, nspin, subspace_env%mixing_store)
853 IF (dft_control%qs_control%gapw) THEN
854 CALL get_qs_env(qs_env=qs_env, rho_atom_set=rho_atom)
855 CALL mixing_init(subspace_env%mixing_method, rho, subspace_env%mixing_store, &
856 para_env, rho_atom=rho_atom)
857 ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
858 CALL charge_mixing_init(subspace_env%mixing_store)
859 ELSEIF (dft_control%qs_control%semi_empirical) THEN
860 cpabort('SE Code not possible')
861 ELSE
862 CALL mixing_init(subspace_env%mixing_method, rho, subspace_env%mixing_store, para_env)
863 END IF
864 END IF
865
866 ene_old = 0.0_dp
867 ene_diff = 0.0_dp
868 IF (output_unit > 0) THEN
869 WRITE (output_unit, "(/T19,A)") '<<<<<<<<< SUBSPACE ROTATION <<<<<<<<<<'
870 WRITE (output_unit, "(T4,A,T13,A,T21,A,T38,A,T51,A,T65,A/,T4,A)") &
871 "In-step", "Time", "Convergence", "Band ene.", "Total ene.", "Energy diff.", repeat("-", 74)
872 END IF
873
874 ! recalculate density matrix here
875
876 ! update of density
877 CALL qs_rho_update_rho(rho, qs_env=qs_env)
878
879 DO iloop = 1, subspace_env%max_iter
880 t1 = m_walltime()
881 converged = .false.
882 ene_old = energy%total
883
884 CALL qs_ks_did_change(ks_env, rho_changed=.true.)
885 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.false., &
886 just_energy=.false., print_active=.false.)
887
888 max_val = 0.0_dp
889 sum_val = 0.0_dp
890 sum_band = 0.0_dp
891 DO ispin = 1, SIZE(matrix_ks)
892 CALL get_mo_set(mo_set=mos(ispin), &
893 nao=nao, &
894 nmo=nmo, &
895 eigenvalues=mo_eigenvalues, &
896 occupation_numbers=mo_occupations, &
897 mo_coeff=mo_coeff)
898
899 !compute C'HC
900 chc => subspace_env%chc_mat(ispin)
901 evec => subspace_env%c_vec(ispin)
902 c0 => subspace_env%c0(ispin)
903 CALL cp_fm_to_fm(mo_coeff, c0)
904 CALL cp_fm_create(work, c0%matrix_struct)
905 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, c0, work, nmo)
906 CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, work, rzero, chc)
907 CALL cp_fm_release(work)
908 !diagonalize C'HC
909 CALL choose_eigv_solver(chc, evec, mo_eigenvalues)
910
911 !rotate the mos by the eigenvectors of C'HC
912 CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, mo_coeff)
913
914 CALL set_mo_occupation(mo_set=mos(ispin), &
915 smear=scf_control%smear)
916
917 ! does not yet handle k-points
918 CALL calculate_density_matrix(mos(ispin), &
919 subspace_env%p_matrix_mix(ispin, 1)%matrix)
920
921 DO i = 1, nmo
922 sum_band = sum_band + mo_eigenvalues(i)*mo_occupations(i)
923 END DO
924
925 !check for self consistency
926 END DO
927
928 IF (subspace_env%mixing_method == direct_mixing_nr) THEN
929 CALL scf_env_density_mixing(subspace_env%p_matrix_mix, &
930 scf_env%mixing_store, rho_ao_kp, para_env, iter_delta, iloop)
931 ELSE
932 CALL self_consistency_check(rho_ao_kp, scf_env%p_delta, para_env, &
933 subspace_env%p_matrix_mix, delta=iter_delta)
934 END IF
935
936 DO ispin = 1, nspin
937 ! does not yet handle k-points
938 CALL dbcsr_copy(rho_ao(ispin)%matrix, subspace_env%p_matrix_mix(ispin, 1)%matrix)
939 END DO
940 ! update of density
941 CALL qs_rho_update_rho(rho, qs_env=qs_env)
942 ! Mixing in reciprocal space
943 IF (subspace_env%mixing_method >= gspace_mixing_nr) THEN
944 CALL gspace_mixing(qs_env, scf_env%mixing_method, subspace_env%mixing_store, &
945 rho, para_env, scf_env%iter_count)
946 END IF
947
948 ene_diff = energy%total - ene_old
949 converged = (abs(ene_diff) < subspace_env%eps_ene .AND. &
950 iter_delta < subspace_env%eps_adapt*scf_env%iter_delta)
951 t2 = m_walltime()
952 IF (output_unit > 0) THEN
953 WRITE (output_unit, "(T4,I5,T11,F8.3,T18,E14.4,T34,F12.5,T46,F16.8,T62,E14.4)") &
954 iloop, t2 - t1, iter_delta, sum_band, energy%total, ene_diff
955 CALL m_flush(output_unit)
956 END IF
957 IF (converged) THEN
958 IF (output_unit > 0) WRITE (output_unit, "(T10,A,I6,A,/)") &
959 " Reached convergence in ", iloop, " iterations "
960 EXIT
961 END IF
962
963 END DO ! iloop
964
965 NULLIFY (subspace_env%p_matrix_mix)
966 DO ispin = 1, nspin
967 ! does not yet handle k-points
968 CALL dbcsr_copy(scf_env%p_mix_new(ispin, 1)%matrix, rho_ao(ispin)%matrix)
969 CALL dbcsr_copy(rho_ao(ispin)%matrix, subspace_env%p_matrix_store(ispin)%matrix)
970
971 DEALLOCATE (eval_first(ispin)%array, occ_first(ispin)%array)
972 END DO
973 DEALLOCATE (eval_first, occ_first)
974
975 CALL timestop(handle)
976
977 END SUBROUTINE do_scf_diag_subspace
978
979! **************************************************************************************************
980!> \brief ...
981!> \param subspace_env ...
982!> \param qs_env ...
983!> \param mos ...
984! **************************************************************************************************
985 SUBROUTINE diag_subspace_allocate(subspace_env, qs_env, mos)
986
987 TYPE(subspace_env_type), POINTER :: subspace_env
988 TYPE(qs_environment_type), POINTER :: qs_env
989 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
990
991 CHARACTER(LEN=*), PARAMETER :: routinen = 'diag_subspace_allocate'
992
993 INTEGER :: handle, i, ispin, nmo, nspin
994 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
995 TYPE(cp_fm_type), POINTER :: mo_coeff
996 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
997 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
998 POINTER :: sab_orb
999
1000 CALL timeset(routinen, handle)
1001
1002 NULLIFY (sab_orb, matrix_s)
1003 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb, &
1004 matrix_s=matrix_s)
1005
1006 nspin = SIZE(mos)
1007! *** allocate p_atrix_store ***
1008 IF (.NOT. ASSOCIATED(subspace_env%p_matrix_store)) THEN
1009 CALL dbcsr_allocate_matrix_set(subspace_env%p_matrix_store, nspin)
1010
1011 DO i = 1, nspin
1012 ALLOCATE (subspace_env%p_matrix_store(i)%matrix)
1013 CALL dbcsr_create(matrix=subspace_env%p_matrix_store(i)%matrix, template=matrix_s(1)%matrix, &
1014 name="DENSITY_STORE", matrix_type=dbcsr_type_symmetric)
1015 CALL cp_dbcsr_alloc_block_from_nbl(subspace_env%p_matrix_store(i)%matrix, &
1016 sab_orb)
1017 CALL dbcsr_set(subspace_env%p_matrix_store(i)%matrix, 0.0_dp)
1018 END DO
1019
1020 END IF
1021
1022 ALLOCATE (subspace_env%chc_mat(nspin))
1023 ALLOCATE (subspace_env%c_vec(nspin))
1024 ALLOCATE (subspace_env%c0(nspin))
1025
1026 DO ispin = 1, nspin
1027 CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1028 CALL cp_fm_create(subspace_env%c0(ispin), mo_coeff%matrix_struct)
1029 NULLIFY (fm_struct_tmp)
1030 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nmo, ncol_global=nmo, &
1031 para_env=mo_coeff%matrix_struct%para_env, &
1032 context=mo_coeff%matrix_struct%context)
1033 CALL cp_fm_create(subspace_env%chc_mat(ispin), fm_struct_tmp, "chc")
1034 CALL cp_fm_create(subspace_env%c_vec(ispin), fm_struct_tmp, "vec")
1035 CALL cp_fm_struct_release(fm_struct_tmp)
1036 END DO
1037
1038 CALL timestop(handle)
1039
1040 END SUBROUTINE diag_subspace_allocate
1041
1042! **************************************************************************************************
1043!> \brief the inner loop of scf, specific to diagonalization without S matrix
1044!> basically, in goes the ks matrix out goes a new p matrix
1045!> \param scf_env ...
1046!> \param mos ...
1047!> \param matrix_ks ...
1048!> \param scf_control ...
1049!> \param scf_section ...
1050!> \param diis_step ...
1051!> \par History
1052!> 03.2006 created [Joost VandeVondele]
1053! **************************************************************************************************
1054 SUBROUTINE do_special_diag(scf_env, mos, matrix_ks, scf_control, &
1055 scf_section, diis_step)
1056
1057 TYPE(qs_scf_env_type), POINTER :: scf_env
1058 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1059 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1060 TYPE(scf_control_type), POINTER :: scf_control
1061 TYPE(section_vals_type), POINTER :: scf_section
1062 LOGICAL, INTENT(INOUT) :: diis_step
1063
1064 INTEGER :: ispin, nspin
1065 LOGICAL :: do_level_shift, use_jacobi
1066 REAL(kind=dp) :: diis_error
1067
1068 nspin = SIZE(matrix_ks)
1069
1070 DO ispin = 1, nspin
1071 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, scf_env%scf_work1(ispin))
1072 END DO
1073 IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis) THEN
1074 CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
1075 scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
1076 scf_control%eps_diis, scf_control%nmixing, &
1077 scf_section=scf_section)
1078 ELSE
1079 diis_step = .false.
1080 END IF
1081
1082 IF ((scf_env%iter_count > 1) .AND. (scf_env%iter_delta < scf_control%diagonalization%eps_jacobi)) THEN
1083 use_jacobi = .true.
1084 ELSE
1085 use_jacobi = .false.
1086 END IF
1087
1088 do_level_shift = ((scf_control%level_shift /= 0.0_dp) .AND. &
1089 ((scf_control%density_guess == core_guess) .OR. (scf_env%iter_count > 1)))
1090 IF (diis_step) THEN
1091 scf_env%iter_param = diis_error
1092 IF (use_jacobi) THEN
1093 scf_env%iter_method = "DIIS/Jacobi"
1094 ELSE
1095 scf_env%iter_method = "DIIS/Diag."
1096 END IF
1097 ELSE
1098 IF (scf_env%mixing_method == 1) THEN
1099 scf_env%iter_param = scf_env%p_mix_alpha
1100 IF (use_jacobi) THEN
1101 scf_env%iter_method = "P_Mix/Jacobi"
1102 ELSE
1103 scf_env%iter_method = "P_Mix/Diag."
1104 END IF
1105 ELSEIF (scf_env%mixing_method > 1) THEN
1106 scf_env%iter_param = scf_env%mixing_store%alpha
1107 IF (use_jacobi) THEN
1108 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Jacobi"
1109 ELSE
1110 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Diag."
1111 END IF
1112 END IF
1113 END IF
1114 scf_env%iter_delta = 0.0_dp
1115
1116 DO ispin = 1, nspin
1117 CALL eigensolver_simple(matrix_ks=scf_env%scf_work1(ispin), &
1118 mo_set=mos(ispin), &
1119 work=scf_env%scf_work2, &
1120 do_level_shift=do_level_shift, &
1121 level_shift=scf_control%level_shift, &
1122 use_jacobi=use_jacobi, &
1123 jacobi_threshold=scf_control%diagonalization%jacobi_threshold)
1124 END DO
1125
1126 CALL set_mo_occupation(mo_array=mos, &
1127 smear=scf_control%smear)
1128
1129 DO ispin = 1, nspin
1130 ! does not yet handle k-points
1131 CALL calculate_density_matrix(mos(ispin), &
1132 scf_env%p_mix_new(ispin, 1)%matrix)
1133 END DO
1134
1135 END SUBROUTINE do_special_diag
1136
1137! **************************************************************************************************
1138!> \brief the inner loop of scf, specific to iterative diagonalization using OT
1139!> with S matrix; basically, in goes the ks matrix out goes a new p matrix
1140!> \param scf_env ...
1141!> \param mos ...
1142!> \param matrix_ks ...
1143!> \param matrix_s ...
1144!> \param scf_control ...
1145!> \param scf_section ...
1146!> \param diis_step ...
1147!> \par History
1148!> 10.2008 created [JGH]
1149! **************************************************************************************************
1150 SUBROUTINE do_ot_diag(scf_env, mos, matrix_ks, matrix_s, &
1151 scf_control, scf_section, diis_step)
1152
1153 TYPE(qs_scf_env_type), POINTER :: scf_env
1154 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1155 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
1156 TYPE(scf_control_type), POINTER :: scf_control
1157 TYPE(section_vals_type), POINTER :: scf_section
1158 LOGICAL, INTENT(INOUT) :: diis_step
1159
1160 INTEGER :: homo, ispin, nmo, nspin
1161 REAL(kind=dp) :: diis_error, eps_iter
1162 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
1163 TYPE(cp_fm_type), POINTER :: mo_coeff
1164
1165 NULLIFY (eigenvalues)
1166
1167 nspin = SIZE(matrix_ks)
1168
1169 DO ispin = 1, nspin
1170 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
1171 scf_env%scf_work1(ispin))
1172 END DO
1173
1174 IF ((scf_env%iter_count > 1) .AND. (.NOT. scf_env%skip_diis)) THEN
1175 CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
1176 scf_env%scf_work2, scf_env%iter_delta, diis_error, diis_step, &
1177 scf_control%eps_diis, scf_control%nmixing, &
1178 s_matrix=matrix_s, &
1179 scf_section=scf_section)
1180 ELSE
1181 diis_step = .false.
1182 END IF
1183
1184 eps_iter = scf_control%diagonalization%eps_iter
1185 IF (diis_step) THEN
1186 scf_env%iter_param = diis_error
1187 scf_env%iter_method = "DIIS/OTdiag"
1188 DO ispin = 1, nspin
1189 CALL copy_fm_to_dbcsr(scf_env%scf_work1(ispin), &
1190 matrix_ks(ispin)%matrix, keep_sparsity=.true.)
1191 END DO
1192 eps_iter = max(eps_iter, scf_control%diagonalization%eps_adapt*diis_error)
1193 ELSE
1194 IF (scf_env%mixing_method == 1) THEN
1195 scf_env%iter_param = scf_env%p_mix_alpha
1196 scf_env%iter_method = "P_Mix/OTdiag."
1197 ELSEIF (scf_env%mixing_method > 1) THEN
1198 scf_env%iter_param = scf_env%mixing_store%alpha
1199 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/OTdiag."
1200 END IF
1201 END IF
1202
1203 scf_env%iter_delta = 0.0_dp
1204
1205 DO ispin = 1, nspin
1206 CALL get_mo_set(mos(ispin), &
1207 mo_coeff=mo_coeff, &
1208 eigenvalues=eigenvalues, &
1209 nmo=nmo, &
1210 homo=homo)
1211 CALL ot_eigensolver(matrix_h=matrix_ks(ispin)%matrix, &
1212 matrix_s=matrix_s(1)%matrix, &
1213 matrix_c_fm=mo_coeff, &
1214 preconditioner=scf_env%ot_preconditioner(1)%preconditioner, &
1215 eps_gradient=eps_iter, &
1216 iter_max=scf_control%diagonalization%max_iter, &
1217 silent=.true., &
1218 ot_settings=scf_control%diagonalization%ot_settings)
1219 CALL calculate_subspace_eigenvalues(mo_coeff, matrix_ks(ispin)%matrix, &
1220 evals_arg=eigenvalues, &
1221 do_rotation=.true.)
1222 CALL copy_fm_to_dbcsr(mos(ispin)%mo_coeff, &
1223 mos(ispin)%mo_coeff_b)
1224 !fm->dbcsr
1225 END DO
1226
1227 CALL set_mo_occupation(mo_array=mos, &
1228 smear=scf_control%smear)
1229
1230 DO ispin = 1, nspin
1231 ! does not yet handle k-points
1232 CALL calculate_density_matrix(mos(ispin), &
1233 scf_env%p_mix_new(ispin, 1)%matrix)
1234 END DO
1235
1236 END SUBROUTINE do_ot_diag
1237
1238! **************************************************************************************************
1239!> \brief Solve a set restricted open Kohn-Sham (ROKS) equations based on the
1240!> alpha and beta Kohn-Sham matrices from unrestricted Kohn-Sham.
1241!> \param scf_env ...
1242!> \param mos ...
1243!> \param matrix_ks ...
1244!> \param matrix_s ...
1245!> \param scf_control ...
1246!> \param scf_section ...
1247!> \param diis_step ...
1248!> \param orthogonal_basis ...
1249!> \par History
1250!> 04.2006 created [MK]
1251!> Revised (01.05.06,MK)
1252!> \note
1253!> this is only a high-spin ROKS.
1254! **************************************************************************************************
1255 SUBROUTINE do_roks_diag(scf_env, mos, matrix_ks, matrix_s, &
1256 scf_control, scf_section, diis_step, &
1257 orthogonal_basis)
1258
1259 ! Literature: - C. C. J. Roothaan, Rev. Mod. Phys. 32, 179 (1960)
1260 ! - M. F. Guest and V. R. Saunders, Mol. Phys. 28(3), 819 (1974)
1261 ! - M. Filatov and S. Shaik, Chem. Phys. Lett. 288, 689 (1998)
1262
1263 TYPE(qs_scf_env_type), POINTER :: scf_env
1264 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
1265 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
1266 TYPE(scf_control_type), POINTER :: scf_control
1267 TYPE(section_vals_type), POINTER :: scf_section
1268 LOGICAL, INTENT(INOUT) :: diis_step
1269 LOGICAL, INTENT(IN) :: orthogonal_basis
1270
1271 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_roks_diag'
1272
1273 INTEGER :: handle, homoa, homob, imo, nalpha, nao, &
1274 nbeta, nmo
1275 REAL(kind=dp) :: diis_error, level_shift_loc
1276 REAL(kind=dp), DIMENSION(:), POINTER :: eiga, eigb, occa, occb
1277 TYPE(cp_fm_type), POINTER :: ksa, ksb, mo2ao, moa, mob, ortho, work
1278
1279! -------------------------------------------------------------------------
1280
1281 CALL timeset(routinen, handle)
1282
1283 IF (scf_env%cholesky_method == cholesky_inverse) THEN
1284 ortho => scf_env%ortho_m1
1285 ELSE
1286 ortho => scf_env%ortho
1287 END IF
1288 work => scf_env%scf_work2
1289
1290 ksa => scf_env%scf_work1(1)
1291 ksb => scf_env%scf_work1(2)
1292
1293 CALL copy_dbcsr_to_fm(matrix_ks(1)%matrix, ksa)
1294 CALL copy_dbcsr_to_fm(matrix_ks(2)%matrix, ksb)
1295
1296 ! Get MO information
1297
1298 CALL get_mo_set(mo_set=mos(1), &
1299 nao=nao, &
1300 nmo=nmo, &
1301 nelectron=nalpha, &
1302 homo=homoa, &
1303 eigenvalues=eiga, &
1304 occupation_numbers=occa, &
1305 mo_coeff=moa)
1306
1307 CALL get_mo_set(mo_set=mos(2), &
1308 nelectron=nbeta, &
1309 homo=homob, &
1310 eigenvalues=eigb, &
1311 occupation_numbers=occb, &
1312 mo_coeff=mob)
1313
1314 ! Define the amount of level-shifting
1315
1316 IF ((scf_control%level_shift /= 0.0_dp) .AND. &
1317 ((scf_control%density_guess == core_guess) .OR. &
1318 (scf_control%density_guess == restart_guess) .OR. &
1319 (scf_env%iter_count > 1))) THEN
1320 level_shift_loc = scf_control%level_shift
1321 ELSE
1322 level_shift_loc = 0.0_dp
1323 END IF
1324
1325 IF ((scf_env%iter_count > 1) .OR. &
1326 (scf_control%density_guess == core_guess) .OR. &
1327 (scf_control%density_guess == restart_guess)) THEN
1328
1329 ! Transform the spin unrestricted alpha and beta Kohn-Sham matrices
1330 ! from AO basis to MO basis: K(MO) = C(T)*K(AO)*C
1331
1332 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, moa, 0.0_dp, work)
1333 CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksa)
1334
1335 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksb, moa, 0.0_dp, work)
1336 CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksb)
1337
1338 ! Combine the spin unrestricted alpha and beta Kohn-Sham matrices
1339 ! in the MO basis
1340
1341 IF (scf_control%roks_scheme == general_roks) THEN
1342 CALL combine_ks_matrices(ksa, ksb, occa, occb, scf_control%roks_f, &
1343 nalpha, nbeta)
1344 ELSE IF (scf_control%roks_scheme == high_spin_roks) THEN
1345 CALL combine_ks_matrices(ksa, ksb, occa, occb, scf_control%roks_parameter)
1346 ELSE
1347 cpabort("Unknown ROKS scheme requested")
1348 END IF
1349
1350 ! Back-transform the restricted open Kohn-Sham matrix from MO basis
1351 ! to AO basis
1352
1353 IF (orthogonal_basis) THEN
1354 ! Q = C
1355 mo2ao => moa
1356 ELSE
1357 ! Q = S*C
1358 mo2ao => mob
1359!MK CALL copy_sm_to_fm(matrix_s(1)%matrix,work)
1360!MK CALL cp_fm_symm("L", "U",nao, nao, 1.0_dp, work, moa, 0.0_dp, mo2ao)
1361 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, moa, mo2ao, nao)
1362 END IF
1363
1364 ! K(AO) = Q*K(MO)*Q(T)
1365
1366 CALL parallel_gemm("N", "T", nao, nao, nao, 1.0_dp, ksa, mo2ao, 0.0_dp, work)
1367 CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, mo2ao, work, 0.0_dp, ksa)
1368
1369 ELSE
1370
1371 ! No transformation matrix available, yet. The closed shell part,
1372 ! i.e. the beta Kohn-Sham matrix in AO basis, is taken.
1373 ! There might be better choices, anyhow.
1374
1375 CALL cp_fm_to_fm(ksb, ksa)
1376
1377 END IF
1378
1379 ! Update DIIS buffer and possibly perform DIIS extrapolation step
1380
1381 IF (scf_env%iter_count > 1) THEN
1382 IF (orthogonal_basis) THEN
1383 CALL qs_diis_b_step(diis_buffer=scf_env%scf_diis_buffer, &
1384 mo_array=mos, &
1385 kc=scf_env%scf_work1, &
1386 sc=work, &
1387 delta=scf_env%iter_delta, &
1388 error_max=diis_error, &
1389 diis_step=diis_step, &
1390 eps_diis=scf_control%eps_diis, &
1391 scf_section=scf_section, &
1392 roks=.true.)
1393 cpassert(scf_env%iter_delta == scf_env%iter_delta)
1394 ELSE
1395 CALL qs_diis_b_step(diis_buffer=scf_env%scf_diis_buffer, &
1396 mo_array=mos, &
1397 kc=scf_env%scf_work1, &
1398 sc=work, &
1399 delta=scf_env%iter_delta, &
1400 error_max=diis_error, &
1401 diis_step=diis_step, &
1402 eps_diis=scf_control%eps_diis, &
1403 scf_section=scf_section, &
1404 s_matrix=matrix_s, &
1405 roks=.true.)
1406 END IF
1407 END IF
1408
1409 IF (diis_step) THEN
1410 scf_env%iter_param = diis_error
1411 scf_env%iter_method = "DIIS/Diag."
1412 ELSE
1413 IF (scf_env%mixing_method == 1) THEN
1414 scf_env%iter_param = scf_env%p_mix_alpha
1415 scf_env%iter_method = "P_Mix/Diag."
1416 ELSEIF (scf_env%mixing_method > 1) THEN
1417 scf_env%iter_param = scf_env%mixing_store%alpha
1418 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Diag."
1419 END IF
1420 END IF
1421
1422 scf_env%iter_delta = 0.0_dp
1423
1424 IF (level_shift_loc /= 0.0_dp) THEN
1425
1426 ! Transform the current Kohn-Sham matrix from AO to MO basis
1427 ! for level-shifting using the current MO set
1428
1429 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, moa, 0.0_dp, work)
1430 CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, moa, work, 0.0_dp, ksa)
1431
1432 ! Apply level-shifting using 50:50 split of the shift (could be relaxed)
1433
1434 DO imo = homob + 1, homoa
1435 CALL cp_fm_add_to_element(ksa, imo, imo, 0.5_dp*level_shift_loc)
1436 END DO
1437 DO imo = homoa + 1, nmo
1438 CALL cp_fm_add_to_element(ksa, imo, imo, level_shift_loc)
1439 END DO
1440
1441 ELSE IF (.NOT. orthogonal_basis) THEN
1442
1443 ! Transform the current Kohn-Sham matrix to an orthogonal basis
1444 SELECT CASE (scf_env%cholesky_method)
1445 CASE (cholesky_reduce)
1446 CALL cp_fm_cholesky_reduce(ksa, ortho)
1447 CASE (cholesky_restore)
1448 CALL cp_fm_uplo_to_full(ksa, work)
1449 CALL cp_fm_cholesky_restore(ksa, nao, ortho, work, &
1450 "SOLVE", pos="RIGHT")
1451 CALL cp_fm_cholesky_restore(work, nao, ortho, ksa, &
1452 "SOLVE", pos="LEFT", transa="T")
1453 CASE (cholesky_inverse)
1454 CALL cp_fm_uplo_to_full(ksa, work)
1455 CALL cp_fm_cholesky_restore(ksa, nao, ortho, work, &
1456 "MULTIPLY", pos="RIGHT")
1457 CALL cp_fm_cholesky_restore(work, nao, ortho, ksa, &
1458 "MULTIPLY", pos="LEFT", transa="T")
1459 CASE (cholesky_off)
1460 CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, ksa, ortho, 0.0_dp, work)
1461 CALL parallel_gemm("N", "N", nao, nao, nao, 1.0_dp, ortho, work, 0.0_dp, ksa)
1462 END SELECT
1463
1464 END IF
1465
1466 ! Diagonalization of the ROKS operator matrix
1467
1468 CALL choose_eigv_solver(ksa, work, eiga)
1469
1470 ! Back-transformation of the orthonormal eigenvectors if needed
1471
1472 IF (level_shift_loc /= 0.0_dp) THEN
1473 ! Use old MO set for back-transformation if level-shifting was applied
1474 CALL cp_fm_to_fm(moa, ortho)
1475 CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, moa)
1476 ELSE
1477 IF (orthogonal_basis) THEN
1478 CALL cp_fm_to_fm(work, moa)
1479 ELSE
1480 SELECT CASE (scf_env%cholesky_method)
1482 CALL cp_fm_cholesky_restore(work, nmo, ortho, moa, "SOLVE")
1483 CASE (cholesky_inverse)
1484 CALL cp_fm_cholesky_restore(work, nmo, ortho, moa, "MULTIPLY")
1485 CASE (cholesky_off)
1486 CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, ortho, work, 0.0_dp, moa)
1487 END SELECT
1488 END IF
1489 END IF
1490
1491 ! Correct MO eigenvalues, if level-shifting was applied
1492
1493 IF (level_shift_loc /= 0.0_dp) THEN
1494 DO imo = homob + 1, homoa
1495 eiga(imo) = eiga(imo) - 0.5_dp*level_shift_loc
1496 END DO
1497 DO imo = homoa + 1, nmo
1498 eiga(imo) = eiga(imo) - level_shift_loc
1499 END DO
1500 END IF
1501
1502 ! Update also the beta MO set
1503
1504 eigb(:) = eiga(:)
1505 CALL cp_fm_to_fm(moa, mob)
1506
1507 ! Calculate the new alpha and beta density matrix
1508
1509 ! does not yet handle k-points
1510 CALL calculate_density_matrix(mos(1), scf_env%p_mix_new(1, 1)%matrix)
1511 CALL calculate_density_matrix(mos(2), scf_env%p_mix_new(2, 1)%matrix)
1512
1513 CALL timestop(handle)
1514
1515 END SUBROUTINE do_roks_diag
1516
1517! **************************************************************************************************
1518!> \brief iterative diagonalization using the block Krylov-space approach
1519!> \param scf_env ...
1520!> \param mos ...
1521!> \param matrix_ks ...
1522!> \param scf_control ...
1523!> \param scf_section ...
1524!> \param check_moconv_only ...
1525!> \param
1526!> \par History
1527!> 05.2009 created [MI]
1528! **************************************************************************************************
1529
1530 SUBROUTINE do_block_krylov_diag(scf_env, mos, matrix_ks, &
1531 scf_control, scf_section, check_moconv_only)
1532
1533 TYPE(qs_scf_env_type), POINTER :: scf_env
1534 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1535 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1536 TYPE(scf_control_type), POINTER :: scf_control
1537 TYPE(section_vals_type), POINTER :: scf_section
1538 LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
1539
1540 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_block_krylov_diag'
1541 REAL(kind=dp), PARAMETER :: rone = 1.0_dp, rzero = 0.0_dp
1542
1543 INTEGER :: handle, homo, ispin, iter, nao, nmo, &
1544 output_unit
1545 LOGICAL :: converged, my_check_moconv_only
1546 REAL(dp) :: eps_iter, t1, t2
1547 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
1548 TYPE(cp_fm_type), POINTER :: c0, c1, chc, evec, ks, mo_coeff, ortho, &
1549 work
1550 TYPE(cp_logger_type), POINTER :: logger
1551
1552 logger => cp_get_default_logger()
1553 CALL timeset(routinen, handle)
1554
1555 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%LANCZOS", &
1556 extension=".scfLog")
1557
1558 my_check_moconv_only = .false.
1559 IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
1560
1561 NULLIFY (mo_coeff, ortho, work, ks)
1562 NULLIFY (mo_eigenvalues)
1563 NULLIFY (c0, c1)
1564
1565 IF (scf_env%cholesky_method == cholesky_inverse) THEN
1566 ortho => scf_env%ortho_m1
1567 ELSE
1568 ortho => scf_env%ortho
1569 END IF
1570 work => scf_env%scf_work2
1571
1572 DO ispin = 1, SIZE(matrix_ks)
1573 CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
1574 scf_env%scf_work1(ispin))
1575 END DO
1576
1577 IF (scf_env%mixing_method == 1) THEN
1578 scf_env%iter_param = scf_env%p_mix_alpha
1579 scf_env%iter_method = "P_Mix/Lanczos"
1580 ELSE
1581! scf_env%iter_param = scf_env%mixing_store%alpha
1582 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Lanc."
1583 END IF
1584
1585 DO ispin = 1, SIZE(matrix_ks)
1586
1587 ks => scf_env%scf_work1(ispin)
1588 CALL cp_fm_uplo_to_full(ks, work)
1589
1590 CALL get_mo_set(mo_set=mos(ispin), &
1591 nao=nao, &
1592 nmo=nmo, &
1593 homo=homo, &
1594 eigenvalues=mo_eigenvalues, &
1595 mo_coeff=mo_coeff)
1596
1597 NULLIFY (c0, c1)
1598 c0 => scf_env%krylov_space%mo_conv(ispin)
1599 c1 => scf_env%krylov_space%mo_refine(ispin)
1600 SELECT CASE (scf_env%cholesky_method)
1601 CASE (cholesky_reduce)
1602 CALL cp_fm_cholesky_reduce(ks, ortho)
1603 CALL cp_fm_uplo_to_full(ks, work)
1604 CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "MULTIPLY")
1605 CASE (cholesky_restore)
1606 CALL cp_fm_cholesky_restore(ks, nao, ortho, work, &
1607 "SOLVE", pos="RIGHT")
1608 CALL cp_fm_cholesky_restore(work, nao, ortho, ks, &
1609 "SOLVE", pos="LEFT", transa="T")
1610 CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "MULTIPLY")
1611 CASE (cholesky_inverse)
1612 CALL cp_fm_cholesky_restore(ks, nao, ortho, work, &
1613 "MULTIPLY", pos="RIGHT")
1614 CALL cp_fm_cholesky_restore(work, nao, ortho, ks, &
1615 "MULTIPLY", pos="LEFT", transa="T")
1616 CALL cp_fm_cholesky_restore(mo_coeff, nmo, ortho, c0, "SOLVE")
1617 END SELECT
1618
1619 scf_env%krylov_space%nmo_nc = nmo
1620 scf_env%krylov_space%nmo_conv = 0
1621
1622 t1 = m_walltime()
1623 IF (output_unit > 0) THEN
1624 WRITE (output_unit, "(/T15,A)") '<<<<<<<<< LANCZOS REFINEMENT <<<<<<<<<<'
1625 WRITE (output_unit, "(T8,A,T15,A,T23,A,T36,A,T49,A,T60,A,/,T8,A)") &
1626 " Spin ", " Cycle ", &
1627 " conv. MOS ", " B2MAX ", " B2MIN ", " Time", repeat("-", 60)
1628 END IF
1629 eps_iter = max(scf_env%krylov_space%eps_conv, scf_env%krylov_space%eps_adapt*scf_env%iter_delta)
1630 iter = 0
1631 converged = .false.
1632 !Check convergence of MOS
1633 IF (my_check_moconv_only) THEN
1634
1635 CALL lanczos_refinement(scf_env%krylov_space, ks, c0, c1, mo_eigenvalues, &
1636 nao, eps_iter, ispin, check_moconv_only=my_check_moconv_only)
1637 t2 = m_walltime()
1638 IF (output_unit > 0) &
1639 WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
1640 ispin, iter, scf_env%krylov_space%nmo_conv, &
1641 scf_env%krylov_space%max_res_norm, scf_env%krylov_space%min_res_norm, t2 - t1
1642
1643 cycle
1644 ELSE
1645 !Block Lanczos refinement
1646 DO iter = 1, scf_env%krylov_space%max_iter
1647 CALL lanczos_refinement_2v(scf_env%krylov_space, ks, c0, c1, mo_eigenvalues, &
1648 nao, eps_iter, ispin)
1649 t2 = m_walltime()
1650 IF (output_unit > 0) THEN
1651 WRITE (output_unit, '(T8,I3,T16,I5,T24,I6,T33,E12.4,2x,E12.4,T60,F8.3)') &
1652 ispin, iter, scf_env%krylov_space%nmo_conv, &
1653 scf_env%krylov_space%max_res_norm, scf_env%krylov_space%min_res_norm, t2 - t1
1654 END IF
1655 t1 = m_walltime()
1656 IF (scf_env%krylov_space%max_res_norm < eps_iter) THEN
1657 converged = .true.
1658 IF (output_unit > 0) WRITE (output_unit, *) &
1659 " Reached convergence in ", iter, " iterations "
1660 EXIT
1661 END IF
1662 END DO
1663
1664 IF (.NOT. converged .AND. output_unit > 0) THEN
1665 WRITE (output_unit, "(T4, A)") " WARNING Lanczos refinement could "// &
1666 "not converge all the mos:"
1667 WRITE (output_unit, "(T40,A,T70,I10)") " number of not converged mos ", &
1668 scf_env%krylov_space%nmo_nc
1669 WRITE (output_unit, "(T40,A,T70,E10.2)") " max norm of the residual ", &
1670 scf_env%krylov_space%max_res_norm
1671
1672 END IF
1673
1674 ! For the moment skip the re-orthogonalization
1675 IF (.false.) THEN
1676 !Re-orthogonalization
1677 NULLIFY (chc, evec)
1678 chc => scf_env%krylov_space%chc_mat(ispin)
1679 evec => scf_env%krylov_space%c_vec(ispin)
1680 CALL parallel_gemm('N', 'N', nao, nmo, nao, rone, ks, c0, rzero, work)
1681 CALL parallel_gemm('T', 'N', nmo, nmo, nao, rone, c0, work, rzero, chc)
1682 !Diagonalize (C^t)HC
1683 CALL choose_eigv_solver(chc, evec, mo_eigenvalues)
1684 !Rotate the C vectors
1685 CALL parallel_gemm('N', 'N', nao, nmo, nmo, rone, c0, evec, rzero, c1)
1686 c0 => scf_env%krylov_space%mo_refine(ispin)
1687 END IF
1688
1689 IF (scf_env%cholesky_method == cholesky_inverse) THEN
1690 CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "MULTIPLY")
1691 ELSE
1692 CALL cp_fm_cholesky_restore(c0, nmo, ortho, mo_coeff, "SOLVE")
1693 END IF
1694
1695 CALL set_mo_occupation(mo_set=mos(ispin), &
1696 smear=scf_control%smear)
1697
1698 ! does not yet handle k-points
1699 CALL calculate_density_matrix(mos(ispin), &
1700 scf_env%p_mix_new(ispin, 1)%matrix)
1701 END IF
1702 END DO ! ispin
1703
1704 IF (output_unit > 0) THEN
1705 WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END LANCZOS REFINEMENT <<<<<<<<<<'
1706 END IF
1707
1708 CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1709 "PRINT%LANCZOS")
1710
1711 CALL timestop(handle)
1712
1713 END SUBROUTINE do_block_krylov_diag
1714
1715! **************************************************************************************************
1716!> \brief iterative diagonalization using the block davidson space approach
1717!> \param qs_env ...
1718!> \param scf_env ...
1719!> \param mos ...
1720!> \param matrix_ks ...
1721!> \param matrix_s ...
1722!> \param scf_control ...
1723!> \param scf_section ...
1724!> \param check_moconv_only ...
1725!> \param
1726!> \par History
1727!> 05.2011 created [MI]
1728! **************************************************************************************************
1729
1730 SUBROUTINE do_block_davidson_diag(qs_env, scf_env, mos, matrix_ks, matrix_s, &
1731 scf_control, scf_section, check_moconv_only)
1732
1733 TYPE(qs_environment_type), POINTER :: qs_env
1734 TYPE(qs_scf_env_type), POINTER :: scf_env
1735 TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT) :: mos
1736 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
1737 TYPE(scf_control_type), POINTER :: scf_control
1738 TYPE(section_vals_type), POINTER :: scf_section
1739 LOGICAL, INTENT(IN), OPTIONAL :: check_moconv_only
1740
1741 CHARACTER(LEN=*), PARAMETER :: routinen = 'do_block_davidson_diag'
1742
1743 INTEGER :: handle, ispin, nspins, output_unit
1744 LOGICAL :: do_prec, my_check_moconv_only
1745 TYPE(cp_logger_type), POINTER :: logger
1746
1747 logger => cp_get_default_logger()
1748 CALL timeset(routinen, handle)
1749
1750 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%DAVIDSON", &
1751 extension=".scfLog")
1752
1753 IF (output_unit > 0) &
1754 WRITE (output_unit, "(/T15,A)") '<<<<<<<<< DAVIDSON ITERATIONS <<<<<<<<<<'
1755
1756 IF (scf_env%mixing_method == 1) THEN
1757 scf_env%iter_param = scf_env%p_mix_alpha
1758 scf_env%iter_method = "P_Mix/Dav."
1759 ELSE
1760 scf_env%iter_param = scf_env%mixing_store%alpha
1761 scf_env%iter_method = trim(scf_env%mixing_store%iter_method)//"/Dav."
1762 END IF
1763
1764 my_check_moconv_only = .false.
1765 IF (PRESENT(check_moconv_only)) my_check_moconv_only = check_moconv_only
1766 do_prec = .false.
1767 IF (scf_env%block_davidson_env(1)%prec_type /= 0 .AND. &
1768 scf_env%iter_count >= scf_env%block_davidson_env(1)%first_prec) THEN
1769 do_prec = .true.
1770 END IF
1771
1772 nspins = SIZE(matrix_ks)
1773
1774 IF (do_prec .AND. (scf_env%iter_count == scf_env%block_davidson_env(1)%first_prec .OR. &
1775 modulo(scf_env%iter_count, scf_env%block_davidson_env(1)%niter_new_prec) == 0)) THEN
1776 CALL restart_preconditioner(qs_env, scf_env%ot_preconditioner, &
1777 prec_type=scf_env%block_davidson_env(1)%prec_type, nspins=nspins)
1778 CALL prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, scf_env%ot_preconditioner, &
1779 scf_env%block_davidson_env(1)%prec_type, &
1780 scf_env%block_davidson_env(1)%solver_type, &
1781 scf_env%block_davidson_env(1)%energy_gap, nspins, &
1782 convert_to_dbcsr=scf_env%block_davidson_env(1)%use_sparse_mos, &
1783 full_mo_set=.true.)
1784 END IF
1785
1786 DO ispin = 1, nspins
1787 IF (scf_env%block_davidson_env(ispin)%use_sparse_mos) THEN
1788 IF (.NOT. do_prec) THEN
1789 CALL generate_extended_space_sparse(scf_env%block_davidson_env(ispin), mos(ispin), &
1790 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit)
1791 ELSE
1792 CALL generate_extended_space_sparse(scf_env%block_davidson_env(ispin), mos(ispin), &
1793 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit, &
1794 scf_env%ot_preconditioner(ispin)%preconditioner)
1795 END IF
1796
1797 ELSE
1798 IF (.NOT. do_prec) THEN
1799 CALL generate_extended_space(scf_env%block_davidson_env(ispin), mos(ispin), &
1800 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit)
1801 ELSE
1802 CALL generate_extended_space(scf_env%block_davidson_env(ispin), mos(ispin), &
1803 matrix_ks(ispin)%matrix, matrix_s(1)%matrix, output_unit, &
1804 scf_env%ot_preconditioner(ispin)%preconditioner)
1805 END IF
1806 END IF
1807 END DO !ispin
1808
1809 CALL set_mo_occupation(mo_array=mos, &
1810 smear=scf_control%smear)
1811
1812 DO ispin = 1, nspins
1813 ! does not yet handle k-points
1814 CALL calculate_density_matrix(mos(ispin), &
1815 scf_env%p_mix_new(ispin, 1)%matrix)
1816 END DO
1817
1818 IF (output_unit > 0) THEN
1819 WRITE (output_unit, "(T15,A/)") '<<<<<<<<< END DAVIDSON ITERATION <<<<<<<<<<'
1820 END IF
1821
1822 CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1823 "PRINT%DAVIDSON")
1824
1825 CALL timestop(handle)
1826
1827 END SUBROUTINE do_block_davidson_diag
1828
1829END MODULE qs_scf_diagonalization
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Basic linear algebra operations for complex full matrices.
subroutine, public cp_cfm_scale_and_add(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b).
subroutine, public cp_cfm_scale_and_add_fm(alpha, matrix_a, beta, matrix_b)
Scale and add two BLACS matrices (a = alpha*a + beta*b). where b is a real matrix (adapted from cp_cf...
used for collecting diagonalization schemes available for cp_cfm_type
Definition cp_cfm_diag.F:14
subroutine, public cp_cfm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
subroutine, public cp_cfm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
General Eigenvalue Problem AX = BXE Use canonical orthogonalization.
Represents a complex full matrix distributed on many processors.
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
subroutine, public cp_cfm_to_fm(msource, mtargetr, mtargeti)
Copy real and imaginary parts of a complex full matrix into separate real-value full matrices.
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_deallocate_matrix(matrix)
...
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_set(matrix, alpha)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Basic linear algebra operations for full matrices.
subroutine, public cp_fm_cholesky_restore(fm_matrix, neig, fm_matrixb, fm_matrixout, op, pos, transa)
...
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
subroutine, public cp_fm_symm(side, uplo, m, n, alpha, matrix_a, matrix_b, beta, matrix_c)
computes matrix_c = beta * matrix_c + alpha * matrix_a * matrix_b computes matrix_c = beta * matrix_c...
various cholesky decomposition related routines
subroutine, public cp_fm_cholesky_reduce(matrix, matrixb, itype)
reduce a matrix pencil A,B to normal form B has to be cholesky decomposed with cp_fm_cholesky_decompo...
used for collecting some of the diagonalization schemes available for cp_fm_type. cp_fm_power also mo...
Definition cp_fm_diag.F:17
subroutine, public cp_fm_geeig(amatrix, bmatrix, eigenvectors, eigenvalues, work)
General Eigenvalue Problem AX = BXE Single option version: Cholesky decomposition of B.
subroutine, public choose_eigv_solver(matrix, eigenvectors, eigenvalues, info)
Choose the Eigensolver depending on which library is available ELPA seems to be unstable for small sy...
Definition cp_fm_diag.F:238
subroutine, public cp_fm_geeig_canon(amatrix, bmatrix, eigenvectors, eigenvalues, work, epseig)
General Eigenvalue Problem AX = BXE Use canonical diagonalization : U*s**(-1/2)
pool for for elements that are retained and released
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
subroutine, public fm_pool_give_back_fm(pool, element)
returns the element to the pool
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_start_copy_general(source, destination, para_env, info)
Initiates the copy operation: get distribution data, post MPI isend and irecvs.
subroutine, public cp_fm_cleanup_copy_general(info)
Completes the copy operation: wait for comms clean up MPI state.
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_add_to_element(matrix, irow_global, icol_global, alpha)
...
subroutine, public cp_fm_finish_copy_general(destination, info)
Completes the copy operation: wait for comms, unpack, clean up MPI state.
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public core_guess
integer, parameter, public cholesky_restore
integer, parameter, public cholesky_dbcsr
integer, parameter, public cholesky_off
integer, parameter, public cholesky_reduce
integer, parameter, public high_spin_roks
integer, parameter, public cholesky_inverse
integer, parameter, public general_roks
integer, parameter, public restart_guess
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Routines needed for kpoint calculation.
subroutine, public rskp_transform(rmatrix, cmatrix, rsmat, ispin, xkp, cell_to_index, sab_nl, is_complex, rs_sign)
Transformation of real space matrices to a kpoint.
subroutine, public kpoint_density_transform(kpoint, denmat, wtype, tempmat, sab_nl, fmwork, for_aux_fit, pmat_ext)
generate real space density matrices in DBCSR format
subroutine, public kpoint_density_matrices(kpoint, energy_weighted, for_aux_fit)
Calculate kpoint density matrices (rho(k), owned by kpoint groups)
subroutine, public kpoint_set_mo_occupation(kpoint, smear, probe)
Given the eigenvalues of all kpoints, calculates the occupation numbers.
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:130
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:147
Interface to the message passing library MPI.
basic linear algebra operations for full matrixes
computes preconditioners, and implements methods to apply them currently used in qs_ot
subroutine, public restart_preconditioner(qs_env, preconditioner, prec_type, nspins)
Allows for a restart of the preconditioner depending on the method it purges all arrays or keeps them...
subroutine, public prepare_preconditioner(qs_env, mos, matrix_ks, matrix_s, ot_preconditioner, prec_type, solver_type, energy_gap, nspins, has_unit_metric, convert_to_dbcsr, chol_type, full_mo_set)
...
collects routines that calculate density matrices
module that contains the definitions of the scf types
integer, parameter, public direct_mixing_nr
integer, parameter, public gspace_mixing_nr
Apply the direct inversion in the iterative subspace (DIIS) of Pulay in the framework of an SCF itera...
Definition qs_diis.F:21
subroutine, public qs_diis_b_info_kp(diis_buffer, ib, nb)
Update info about the current buffer step ib and the current number of buffers nb.
Definition qs_diis.F:989
subroutine, public qs_diis_b_calc_err_kp(diis_buffer, ib, mos, kc, sc, ispin, ikp, nkp_local, scf_section)
Calculate and store the error for a given k-point.
Definition qs_diis.F:1013
subroutine, public qs_diis_b_step_kp(diis_buffer, coeffs, ib, nb, delta, error_max, diis_step, eps_diis, nspin, nkp, nkp_local, nmixing, scf_section, para_env)
Update the SCF DIIS buffer, and if appropriate does a diis step, for k-points.
Definition qs_diis.F:1100
subroutine, public qs_diis_b_step(diis_buffer, mo_array, kc, sc, delta, error_max, diis_step, eps_diis, nmixing, s_matrix, scf_section, roks)
Update the SCF DIIS buffer, and if appropriate does a diis step.
Definition qs_diis.F:230
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, tb_tblite)
Get the QUICKSTEP environment.
subroutine, public gspace_mixing(qs_env, mixing_method, mixing_store, rho, para_env, iter_count)
Driver for the g-space mixing, calls the proper routine given the requested method.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
wrapper for the pools of matrixes
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
elemental subroutine, public charge_mixing_init(mixing_store)
initialiation needed when charge mixing is used
subroutine, public mixing_init(mixing_method, rho, mixing_store, para_env, rho_atom)
initialiation needed when gspace mixing is used
subroutine, public mixing_allocate(qs_env, mixing_method, p_mix_new, p_delta, nspins, mixing_store)
allocation needed when density mixing is used
subroutine, public self_consistency_check(rho_ao, p_delta, para_env, p_out, delta)
...
collects routines that perform operations directly related to MOs
Set occupation of molecular orbitals.
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Define the neighbor list data types and the corresponding functionality.
an eigen-space solver for the generalised symmetric eigenvalue problem for sparse matrices,...
subroutine, public ot_eigensolver(matrix_h, matrix_s, matrix_orthogonal_space_fm, matrix_c_fm, preconditioner, eps_gradient, iter_max, size_ortho_space, silent, ot_settings)
...
methods of the rho structure (defined in qs_rho_types)
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
module that contains the algorithms to perform an itrative diagonalization by the block-Davidson appr...
subroutine, public generate_extended_space_sparse(bdav_env, mo_set, matrix_h, matrix_s, output_unit, preconditioner)
...
subroutine, public generate_extended_space(bdav_env, mo_set, matrix_h, matrix_s, output_unit, preconditioner)
...
Different diagonalization schemes that can be used for the iterative solution of the eigenvalue probl...
subroutine, public general_eigenproblem(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step)
the inner loop of scf, specific to diagonalization with S matrix basically, in goes the ks matrix out...
subroutine, public diag_subspace_allocate(subspace_env, qs_env, mos)
...
subroutine, public do_ot_diag(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step)
the inner loop of scf, specific to iterative diagonalization using OT with S matrix; basically,...
subroutine, public do_block_davidson_diag(qs_env, scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, check_moconv_only)
iterative diagonalization using the block davidson space approach
subroutine, public do_roks_diag(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step, orthogonal_basis)
Solve a set restricted open Kohn-Sham (ROKS) equations based on the alpha and beta Kohn-Sham matrices...
subroutine, public do_scf_diag_subspace(qs_env, scf_env, subspace_env, mos, rho, ks_env, scf_section, scf_control)
inner loop within MOS subspace, to refine occupation and density, before next diagonalization of the ...
subroutine, public do_block_krylov_diag(scf_env, mos, matrix_ks, scf_control, scf_section, check_moconv_only)
iterative diagonalization using the block Krylov-space approach
subroutine, public do_special_diag(scf_env, mos, matrix_ks, scf_control, scf_section, diis_step)
the inner loop of scf, specific to diagonalization without S matrix basically, in goes the ks matrix ...
subroutine, public do_general_diag(scf_env, mos, matrix_ks, matrix_s, scf_control, scf_section, diis_step, probe)
...
subroutine, public do_general_diag_kp(matrix_ks, matrix_s, kpoints, scf_env, scf_control, update_p, diis_step, diis_error, qs_env, probe)
Kpoint diagonalization routine Transforms matrices to kpoint, distributes kpoint groups,...
module that contains the algorithms to perform an itrative diagonalization by the block-Lanczos appro...
subroutine, public lanczos_refinement(krylov_space, ks, c0, c1, eval, nao, eps_iter, ispin, check_moconv_only)
lanczos refinement by blocks of not-converged MOs
subroutine, public lanczos_refinement_2v(krylov_space, ks, c0, c1, eval, nao, eps_iter, ispin, check_moconv_only)
...
groups fairly general SCF methods, so that modules other than qs_scf can use them too split off from ...
subroutine, public eigensolver_simple(matrix_ks, mo_set, work, do_level_shift, level_shift, use_jacobi, jacobi_threshold)
...
subroutine, public eigensolver_dbcsr(matrix_ks, matrix_ks_fm, mo_set, ortho_dbcsr, ksbuf1, ksbuf2)
...
subroutine, public scf_env_density_mixing(p_mix_new, mixing_store, rho_ao, para_env, iter_delta, iter_count, diis, invert)
perform (if requested) a density mixing
subroutine, public eigensolver(matrix_ks_fm, mo_set, ortho, work, cholesky_method, do_level_shift, level_shift, matrix_u_fm, use_jacobi)
Diagonalise the Kohn-Sham matrix to get a new set of MO eigen- vectors and MO eigenvalues....
subroutine, public eigensolver_symm(matrix_ks_fm, mo_set, ortho, work, do_level_shift, level_shift, matrix_u_fm, use_jacobi, jacobi_threshold, ortho_red, work_red, matrix_ks_fm_red, matrix_u_fm_red)
...
module that contains the definitions of the scf types
parameters that control an scf iteration
represent a pointer to a 1d array
Represent a complex full matrix.
keeps the information about the structure of a full matrix
Stores the state of a copy between cp_fm_start_copy_general and cp_fm_finish_copy_general.
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Keeps information about a specific k-point.
Contains information about kpoints.
stores all the informations relevant to an mpi environment
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
container for the pools of matrixes used by qs
keeps the density in various representations, keeping track of which ones are valid.