LCOV - code coverage report
Current view: top level - src - dm_ls_scf.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 92.6 % 476 441
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 12 12

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Routines for a linear scaling quickstep SCF run based on the density
      10              : !>        matrix
      11              : !> \par History
      12              : !>       2010.10 created [Joost VandeVondele]
      13              : !> \author Joost VandeVondele
      14              : ! **************************************************************************************************
      15              : MODULE dm_ls_scf
      16              :    USE arnoldi_api,                     ONLY: arnoldi_extremal
      17              :    USE bibliography,                    ONLY: Kolafa2004,&
      18              :                                               Kuhne2007,&
      19              :                                               cite_reference
      20              :    USE cp_control_types,                ONLY: dft_control_type
      21              :    USE cp_dbcsr_api,                    ONLY: &
      22              :         dbcsr_add, dbcsr_binary_read, dbcsr_binary_write, dbcsr_copy, dbcsr_create, &
      23              :         dbcsr_distribution_type, dbcsr_filter, dbcsr_get_info, dbcsr_get_occupation, &
      24              :         dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, &
      25              :         dbcsr_type_no_symmetry
      26              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_checksum
      27              :    USE cp_external_control,             ONLY: external_control
      28              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      29              :                                               cp_logger_get_default_unit_nr,&
      30              :                                               cp_logger_type
      31              :    USE dm_ls_chebyshev,                 ONLY: compute_chebyshev
      32              :    USE dm_ls_scf_create,                ONLY: ls_scf_create
      33              :    USE dm_ls_scf_curvy,                 ONLY: deallocate_curvy_data,&
      34              :                                               dm_ls_curvy_optimization
      35              :    USE dm_ls_scf_methods,               ONLY: apply_matrix_preconditioner,&
      36              :                                               compute_homo_lumo,&
      37              :                                               density_matrix_sign,&
      38              :                                               density_matrix_sign_fixed_mu,&
      39              :                                               density_matrix_tc2,&
      40              :                                               density_matrix_trs4,&
      41              :                                               ls_scf_init_matrix_S
      42              :    USE dm_ls_scf_qs,                    ONLY: &
      43              :         ls_nonscf_energy, ls_nonscf_ks, ls_scf_dm_to_ks, ls_scf_init_qs, ls_scf_qs_atomic_guess, &
      44              :         matrix_ls_create, matrix_ls_to_qs, matrix_qs_to_ls, rho_mixing_ls_init
      45              :    USE dm_ls_scf_types,                 ONLY: ls_scf_env_type
      46              :    USE ec_env_types,                    ONLY: energy_correction_type
      47              :    USE input_constants,                 ONLY: ls_cluster_atomic,&
      48              :                                               ls_scf_pexsi,&
      49              :                                               ls_scf_sign,&
      50              :                                               ls_scf_tc2,&
      51              :                                               ls_scf_trs4,&
      52              :                                               transport_transmission
      53              :    USE input_section_types,             ONLY: section_vals_type
      54              :    USE iterate_matrix,                  ONLY: purify_mcweeny
      55              :    USE kinds,                           ONLY: default_path_length,&
      56              :                                               default_string_length,&
      57              :                                               dp
      58              :    USE machine,                         ONLY: m_flush,&
      59              :                                               m_walltime
      60              :    USE mathlib,                         ONLY: binomial
      61              :    USE molecule_types,                  ONLY: molecule_type
      62              :    USE pao_main,                        ONLY: pao_optimization_end,&
      63              :                                               pao_optimization_start,&
      64              :                                               pao_post_scf,&
      65              :                                               pao_update
      66              :    USE pexsi_methods,                   ONLY: density_matrix_pexsi,&
      67              :                                               pexsi_finalize_scf,&
      68              :                                               pexsi_init_scf,&
      69              :                                               pexsi_set_convergence_tolerance,&
      70              :                                               pexsi_to_qs
      71              :    USE qs_diis,                         ONLY: qs_diis_b_clear_sparse,&
      72              :                                               qs_diis_b_create_sparse,&
      73              :                                               qs_diis_b_step_4lscf
      74              :    USE qs_diis_types,                   ONLY: qs_diis_b_release_sparse,&
      75              :                                               qs_diis_buffer_type_sparse
      76              :    USE qs_environment_types,            ONLY: get_qs_env,&
      77              :                                               qs_environment_type
      78              :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      79              :    USE qs_nonscf_utils,                 ONLY: qs_nonscf_print_summary
      80              :    USE qs_scf_post_gpw,                 ONLY: qs_scf_post_moments,&
      81              :                                               write_mo_free_results
      82              :    USE qs_scf_post_tb,                  ONLY: scf_post_calculation_tb
      83              :    USE transport,                       ONLY: external_scf_method,&
      84              :                                               transport_initialize
      85              :    USE transport_env_types,             ONLY: transport_env_type
      86              : #include "./base/base_uses.f90"
      87              : 
      88              :    IMPLICIT NONE
      89              : 
      90              :    PRIVATE
      91              : 
      92              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf'
      93              : 
      94              :    PUBLIC :: calculate_w_matrix_ls, ls_scf, post_scf_sparsities
      95              : 
      96              : CONTAINS
      97              : 
      98              : ! **************************************************************************************************
      99              : !> \brief perform an linear scaling scf procedure: entry point
     100              : !>
     101              : !> \param qs_env ...
     102              : !> \param nonscf ...
     103              : !> \par History
     104              : !>       2010.10 created [Joost VandeVondele]
     105              : !> \author Joost VandeVondele
     106              : ! **************************************************************************************************
     107          664 :    SUBROUTINE ls_scf(qs_env, nonscf)
     108              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     109              :       LOGICAL, INTENT(IN), OPTIONAL                      :: nonscf
     110              : 
     111              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ls_scf'
     112              : 
     113              :       INTEGER                                            :: handle
     114              :       LOGICAL                                            :: do_scf, pao_is_done
     115              :       TYPE(ls_scf_env_type), POINTER                     :: ls_scf_env
     116              : 
     117          664 :       CALL timeset(routineN, handle)
     118          664 :       do_scf = .TRUE.
     119          664 :       IF (PRESENT(nonscf)) do_scf = .NOT. nonscf
     120              : 
     121              :       ! Moved here from qs_environment to remove dependencies
     122          664 :       CALL ls_scf_create(qs_env)
     123          664 :       CALL get_qs_env(qs_env, ls_scf_env=ls_scf_env)
     124              : 
     125          664 :       IF (do_scf) THEN
     126          598 :          CALL pao_optimization_start(qs_env, ls_scf_env)
     127          598 :          pao_is_done = .FALSE.
     128         1414 :          DO WHILE (.NOT. pao_is_done)
     129          816 :             CALL ls_scf_init_scf(qs_env, ls_scf_env, .FALSE.)
     130          816 :             CALL pao_update(qs_env, ls_scf_env, pao_is_done)
     131          816 :             CALL ls_scf_main(qs_env, ls_scf_env, .FALSE.)
     132          816 :             CALL pao_post_scf(qs_env, ls_scf_env, pao_is_done)
     133          816 :             CALL ls_scf_post(qs_env, ls_scf_env)
     134              :          END DO
     135          598 :          CALL pao_optimization_end(ls_scf_env)
     136              :       ELSE
     137           66 :          CALL ls_scf_init_scf(qs_env, ls_scf_env, .TRUE.)
     138           66 :          CALL ls_scf_main(qs_env, ls_scf_env, .TRUE.)
     139           66 :          CALL ls_scf_post(qs_env, ls_scf_env)
     140              :       END IF
     141              : 
     142          664 :       CALL timestop(handle)
     143              : 
     144          664 :    END SUBROUTINE ls_scf
     145              : 
     146              : ! **************************************************************************************************
     147              : !> \brief initialization needed for scf
     148              : !> \param qs_env ...
     149              : !> \param ls_scf_env ...
     150              : !> \param nonscf ...
     151              : !> \par History
     152              : !>       2010.10 created [Joost VandeVondele]
     153              : !> \author Joost VandeVondele
     154              : ! **************************************************************************************************
     155          882 :    SUBROUTINE ls_scf_init_scf(qs_env, ls_scf_env, nonscf)
     156              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     157              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
     158              :       LOGICAL, INTENT(IN)                                :: nonscf
     159              : 
     160              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ls_scf_init_scf'
     161              : 
     162              :       INTEGER                                            :: handle, ispin, nspin, unit_nr
     163              :       TYPE(cp_logger_type), POINTER                      :: logger
     164          882 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s, matrix_w
     165              :       TYPE(dft_control_type), POINTER                    :: dft_control
     166          882 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     167              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     168              :       TYPE(section_vals_type), POINTER                   :: input
     169              : 
     170          882 :       CALL timeset(routineN, handle)
     171              : 
     172              :       ! get a useful output_unit
     173          882 :       logger => cp_get_default_logger()
     174          882 :       IF (logger%para_env%is_source()) THEN
     175          441 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     176              :       ELSE
     177              :          unit_nr = -1
     178              :       END IF
     179              : 
     180              :       ! get basic quantities from the qs_env
     181              :       CALL get_qs_env(qs_env, nelectron_total=ls_scf_env%nelectron_total, &
     182              :                       matrix_s=matrix_s, &
     183              :                       matrix_w=matrix_w, &
     184              :                       ks_env=ks_env, &
     185              :                       dft_control=dft_control, &
     186              :                       molecule_set=molecule_set, &
     187              :                       input=input, &
     188              :                       has_unit_metric=ls_scf_env%has_unit_metric, &
     189              :                       para_env=ls_scf_env%para_env, &
     190          882 :                       nelectron_spin=ls_scf_env%nelectron_spin)
     191              : 
     192              :       ! needs forces ? There might be a better way to flag this
     193          882 :       ls_scf_env%calculate_forces = ASSOCIATED(matrix_w)
     194              : 
     195              :       ! some basic initialization of the QS side of things
     196          882 :       CALL ls_scf_init_qs(qs_env)
     197              : 
     198              :       ! create the matrix template for use in the ls procedures
     199              :       CALL matrix_ls_create(matrix_ls=ls_scf_env%matrix_s, matrix_qs=matrix_s(1)%matrix, &
     200          882 :                             ls_mstruct=ls_scf_env%ls_mstruct)
     201              : 
     202          882 :       nspin = ls_scf_env%nspins
     203          882 :       IF (ALLOCATED(ls_scf_env%matrix_p)) THEN
     204         1098 :          DO ispin = 1, SIZE(ls_scf_env%matrix_p)
     205         1098 :             CALL dbcsr_release(ls_scf_env%matrix_p(ispin))
     206              :          END DO
     207              :       ELSE
     208         1358 :          ALLOCATE (ls_scf_env%matrix_p(nspin))
     209              :       END IF
     210              : 
     211         1784 :       DO ispin = 1, nspin
     212              :          CALL dbcsr_create(ls_scf_env%matrix_p(ispin), template=ls_scf_env%matrix_s, &
     213         1784 :                            matrix_type=dbcsr_type_no_symmetry)
     214              :       END DO
     215              : 
     216         3548 :       ALLOCATE (ls_scf_env%matrix_ks(nspin))
     217         1784 :       DO ispin = 1, nspin
     218              :          CALL dbcsr_create(ls_scf_env%matrix_ks(ispin), template=ls_scf_env%matrix_s, &
     219         1784 :                            matrix_type=dbcsr_type_no_symmetry)
     220              :       END DO
     221              : 
     222              :       ! set up matrix S, and needed functions of S
     223          882 :       CALL ls_scf_init_matrix_s(matrix_s(1)%matrix, ls_scf_env)
     224              : 
     225              :       ! get the initial guess for the SCF
     226          882 :       CALL ls_scf_initial_guess(qs_env, ls_scf_env, nonscf)
     227              : 
     228          882 :       IF (ls_scf_env%do_rho_mixing) THEN
     229            0 :          CALL rho_mixing_ls_init(qs_env, ls_scf_env)
     230              :       END IF
     231              : 
     232          882 :       IF (ls_scf_env%do_pexsi) THEN
     233            0 :          CALL pexsi_init_scf(ks_env, ls_scf_env%pexsi, matrix_s(1)%matrix)
     234              :       END IF
     235              : 
     236          882 :       IF (qs_env%do_transport) THEN
     237            0 :          CALL transport_initialize(ks_env, qs_env%transport_env, matrix_s(1)%matrix)
     238              :       END IF
     239              : 
     240          882 :       CALL timestop(handle)
     241              : 
     242          882 :    END SUBROUTINE ls_scf_init_scf
     243              : 
     244              : ! **************************************************************************************************
     245              : !> \brief deal with the scf initial guess
     246              : !> \param qs_env ...
     247              : !> \param ls_scf_env ...
     248              : !> \param nonscf ...
     249              : !> \par History
     250              : !>       2012.11 created [Joost VandeVondele]
     251              : !> \author Joost VandeVondele
     252              : ! **************************************************************************************************
     253         1252 :    SUBROUTINE ls_scf_initial_guess(qs_env, ls_scf_env, nonscf)
     254              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     255              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
     256              :       LOGICAL, INTENT(IN)                                :: nonscf
     257              : 
     258              :       CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_initial_guess'
     259              :       INTEGER, PARAMETER                                 :: aspc_guess = 2, atomic_guess = 1, &
     260              :                                                             restart_guess = 3
     261              : 
     262              :       CHARACTER(LEN=default_path_length)                 :: file_name, project_name
     263              :       INTEGER                                            :: handle, iaspc, initial_guess_type, &
     264              :                                                             ispin, istore, naspc, unit_nr
     265              :       REAL(KIND=dp)                                      :: alpha, cs_pos
     266              :       TYPE(cp_logger_type), POINTER                      :: logger
     267              :       TYPE(dbcsr_distribution_type)                      :: dist
     268              :       TYPE(dbcsr_type)                                   :: matrix_tmp1
     269              : 
     270          512 :       IF (ls_scf_env%do_pao) RETURN ! pao has its own initial guess
     271              : 
     272          370 :       CALL timeset(routineN, handle)
     273              : 
     274              :       ! get a useful output_unit
     275          370 :       logger => cp_get_default_logger()
     276          370 :       IF (logger%para_env%is_source()) THEN
     277          185 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     278              :       ELSE
     279              :          unit_nr = -1
     280              :       END IF
     281              : 
     282          185 :       IF (unit_nr > 0) WRITE (unit_nr, '()')
     283              :       ! if there is no history go for the atomic guess, otherwise extrapolate the dm history
     284          370 :       IF (ls_scf_env%scf_history%istore == 0) THEN
     285          244 :          IF (ls_scf_env%restart_read) THEN
     286              :             initial_guess_type = restart_guess
     287              :          ELSE
     288              :             initial_guess_type = atomic_guess
     289              :          END IF
     290              :       ELSE
     291              :          initial_guess_type = aspc_guess
     292              :       END IF
     293              : 
     294              :       ! how to get the initial guess
     295              :       SELECT CASE (initial_guess_type)
     296              :       CASE (atomic_guess)
     297          240 :          CALL ls_scf_qs_atomic_guess(qs_env, ls_scf_env, ls_scf_env%energy_init, nonscf)
     298          240 :          IF (unit_nr > 0) WRITE (unit_nr, '()')
     299              :       CASE (restart_guess)
     300            4 :          project_name = logger%iter_info%project_name
     301            8 :          DO ispin = 1, SIZE(ls_scf_env%matrix_p)
     302            4 :             WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_", ispin, "_RESTART.dm"
     303            4 :             CALL dbcsr_get_info(ls_scf_env%matrix_p(1), distribution=dist)
     304            4 :             CALL dbcsr_binary_read(file_name, distribution=dist, matrix_new=ls_scf_env%matrix_p(ispin))
     305            4 :             cs_pos = dbcsr_checksum(ls_scf_env%matrix_p(ispin), pos=.TRUE.)
     306           12 :             IF (unit_nr > 0) THEN
     307            2 :                WRITE (unit_nr, '(T2,A,E20.8)') "Read restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
     308              :             END IF
     309              :          END DO
     310              : 
     311              :          ! directly go to computing the corresponding energy and ks matrix
     312            4 :          IF (nonscf) THEN
     313            0 :             CALL ls_nonscf_ks(qs_env, ls_scf_env, ls_scf_env%energy_init)
     314              :          ELSE
     315            4 :             CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, ls_scf_env%energy_init, iscf=0)
     316              :          END IF
     317              :       CASE (aspc_guess)
     318          126 :          CALL cite_reference(Kolafa2004)
     319          126 :          CALL cite_reference(Kuhne2007)
     320          126 :          naspc = MIN(ls_scf_env%scf_history%istore, ls_scf_env%scf_history%nstore)
     321          258 :          DO ispin = 1, SIZE(ls_scf_env%matrix_p)
     322              :             ! actual extrapolation
     323          132 :             CALL dbcsr_set(ls_scf_env%matrix_p(ispin), 0.0_dp)
     324          624 :             DO iaspc = 1, naspc
     325              :                alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* &
     326          366 :                        binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
     327          366 :                istore = MOD(ls_scf_env%scf_history%istore - iaspc, ls_scf_env%scf_history%nstore) + 1
     328          498 :                CALL dbcsr_add(ls_scf_env%matrix_p(ispin), ls_scf_env%scf_history%matrix(ispin, istore), 1.0_dp, alpha)
     329              :             END DO
     330              :          END DO
     331              :       END SELECT
     332              : 
     333              :       ! which cases need getting purified and non-orthogonal ?
     334          126 :       SELECT CASE (initial_guess_type)
     335              :       CASE (atomic_guess, restart_guess)
     336              :          ! do nothing
     337              :       CASE (aspc_guess)
     338              :          ! purification can't be done on the pexsi matrix, which is not necessarily idempotent,
     339              :          ! and not stored in an ortho basis form
     340          126 :          IF (.NOT. (ls_scf_env%do_pexsi)) THEN
     341          258 :             DO ispin = 1, SIZE(ls_scf_env%matrix_p)
     342              :                ! linear combination of P's is not idempotent. A bit of McWeeny is needed to ensure it is again
     343          132 :                IF (SIZE(ls_scf_env%matrix_p) == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 0.5_dp)
     344              :                ! to ensure that noisy blocks do not build up during MD (in particular with curvy) filter that guess a bit more
     345          132 :                CALL dbcsr_filter(ls_scf_env%matrix_p(ispin), ls_scf_env%eps_filter**(2.0_dp/3.0_dp))
     346          132 :                CALL purify_mcweeny(ls_scf_env%matrix_p(ispin:ispin), ls_scf_env%eps_filter, 3)
     347          132 :                IF (SIZE(ls_scf_env%matrix_p) == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 2.0_dp)
     348              : 
     349          258 :                IF (ls_scf_env%use_s_sqrt) THEN
     350              :                   ! need to get P in the non-orthogonal basis if it was stored differently
     351              :                   CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, &
     352          132 :                                     matrix_type=dbcsr_type_no_symmetry)
     353              :                   CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_p(ispin), &
     354          132 :                                       0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
     355              :                   CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt_inv, &
     356              :                                       0.0_dp, ls_scf_env%matrix_p(ispin), &
     357          132 :                                       filter_eps=ls_scf_env%eps_filter)
     358          132 :                   CALL dbcsr_release(matrix_tmp1)
     359              : 
     360          132 :                   IF (ls_scf_env%has_s_preconditioner) THEN
     361              :                      CALL apply_matrix_preconditioner(ls_scf_env%matrix_p(ispin), "forward", &
     362          126 :                                                       ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
     363              :                   END IF
     364              :                END IF
     365              :             END DO
     366              :          END IF
     367              : 
     368              :          ! compute corresponding energy and ks matrix
     369          496 :          IF (nonscf) THEN
     370           60 :             CALL ls_nonscf_ks(qs_env, ls_scf_env, ls_scf_env%energy_init)
     371              :          ELSE
     372           66 :             CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, ls_scf_env%energy_init, iscf=0)
     373              :          END IF
     374              :       END SELECT
     375              : 
     376          370 :       IF (unit_nr > 0) THEN
     377          185 :          WRITE (unit_nr, '(T2,A,F20.9)') "Energy with the initial guess:", ls_scf_env%energy_init
     378          185 :          WRITE (unit_nr, '()')
     379              :       END IF
     380              : 
     381          370 :       CALL timestop(handle)
     382              : 
     383          882 :    END SUBROUTINE ls_scf_initial_guess
     384              : 
     385              : ! **************************************************************************************************
     386              : !> \brief store a history of matrices for later use in ls_scf_initial_guess
     387              : !> \param ls_scf_env ...
     388              : !> \par History
     389              : !>       2012.11 created [Joost VandeVondele]
     390              : !> \author Joost VandeVondele
     391              : ! **************************************************************************************************
     392          370 :    SUBROUTINE ls_scf_store_result(ls_scf_env)
     393              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
     394              : 
     395              :       CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_store_result'
     396              : 
     397              :       CHARACTER(LEN=default_path_length)                 :: file_name, project_name
     398              :       INTEGER                                            :: handle, ispin, istore, unit_nr
     399              :       REAL(KIND=dp)                                      :: cs_pos
     400              :       TYPE(cp_logger_type), POINTER                      :: logger
     401              :       TYPE(dbcsr_type)                                   :: matrix_tmp1
     402              : 
     403          370 :       CALL timeset(routineN, handle)
     404              : 
     405              :       ! get a useful output_unit
     406          370 :       logger => cp_get_default_logger()
     407          370 :       IF (logger%para_env%is_source()) THEN
     408          185 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     409              :       ELSE
     410              :          unit_nr = -1
     411              :       END IF
     412              : 
     413          370 :       IF (ls_scf_env%restart_write) THEN
     414           12 :          DO ispin = 1, SIZE(ls_scf_env%matrix_p)
     415            6 :             project_name = logger%iter_info%project_name
     416            6 :             WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_", ispin, "_RESTART.dm"
     417            6 :             cs_pos = dbcsr_checksum(ls_scf_env%matrix_p(ispin), pos=.TRUE.)
     418            6 :             IF (unit_nr > 0) THEN
     419            3 :                WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
     420              :             END IF
     421            6 :             IF (ls_scf_env%do_transport .OR. ls_scf_env%do_pexsi) THEN
     422            0 :                IF (unit_nr > 0) THEN
     423            0 :                   WRITE (unit_nr, '(T6,A)') "The restart DM "//TRIM(file_name)//" has the sparsity of S, therefore,"
     424            0 :                   WRITE (unit_nr, '(T6,A)') "not compatible with methods that require a full DM! "
     425              :                END IF
     426              :             END IF
     427           12 :             CALL dbcsr_binary_write(ls_scf_env%matrix_p(ispin), file_name)
     428              :          END DO
     429              :       END IF
     430              : 
     431          370 :       IF (ls_scf_env%scf_history%nstore > 0) THEN
     432          362 :          ls_scf_env%scf_history%istore = ls_scf_env%scf_history%istore + 1
     433          744 :          DO ispin = 1, SIZE(ls_scf_env%matrix_p)
     434          382 :             istore = MOD(ls_scf_env%scf_history%istore - 1, ls_scf_env%scf_history%nstore) + 1
     435          382 :             CALL dbcsr_copy(ls_scf_env%scf_history%matrix(ispin, istore), ls_scf_env%matrix_p(ispin))
     436              : 
     437              :             ! if we have the sqrt around, we use it to go to the orthogonal basis
     438          744 :             IF (ls_scf_env%use_s_sqrt) THEN
     439              :                ! usually sqrt(S) * P * sqrt(S) should be available, or could be stored at least,
     440              :                ! so that the next multiplications could be saved.
     441              :                CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, &
     442          378 :                                  matrix_type=dbcsr_type_no_symmetry)
     443              : 
     444          378 :                IF (ls_scf_env%has_s_preconditioner) THEN
     445              :                   CALL apply_matrix_preconditioner(ls_scf_env%scf_history%matrix(ispin, istore), "backward", &
     446          332 :                                                    ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
     447              :                END IF
     448              :                CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt, ls_scf_env%scf_history%matrix(ispin, istore), &
     449          378 :                                    0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
     450              :                CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt, &
     451              :                                    0.0_dp, ls_scf_env%scf_history%matrix(ispin, istore), &
     452          378 :                                    filter_eps=ls_scf_env%eps_filter)
     453          378 :                CALL dbcsr_release(matrix_tmp1)
     454              :             END IF
     455              : 
     456              :          END DO
     457              :       END IF
     458              : 
     459          370 :       CALL timestop(handle)
     460              : 
     461          370 :    END SUBROUTINE ls_scf_store_result
     462              : 
     463              : ! **************************************************************************************************
     464              : !> \brief Main SCF routine. Can we keep it clean ?
     465              : !> \param qs_env ...
     466              : !> \param ls_scf_env ...
     467              : !> \param nonscf ...
     468              : !> \par History
     469              : !>       2010.10 created [Joost VandeVondele]
     470              : !> \author Joost VandeVondele
     471              : ! **************************************************************************************************
     472          882 :    SUBROUTINE ls_scf_main(qs_env, ls_scf_env, nonscf)
     473              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     474              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
     475              :       LOGICAL, INTENT(IN), OPTIONAL                      :: nonscf
     476              : 
     477              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ls_scf_main'
     478              : 
     479              :       INTEGER                                            :: handle, iscf, ispin, &
     480              :                                                             nelectron_spin_real, nmixing, nspin, &
     481              :                                                             unit_nr
     482              :       LOGICAL :: check_convergence, diis_step, do_transport, extra_scf, maxscf_reached, &
     483              :          scf_converged, should_stop, transm_maxscf_reached, transm_scf_converged
     484              :       REAL(KIND=dp)                                      :: energy_diff, energy_new, energy_old, &
     485              :                                                             eps_diis, t1, t2, tdiag
     486              :       TYPE(cp_logger_type), POINTER                      :: logger
     487          882 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_s
     488          882 :       TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:)        :: matrix_ks_deviation, matrix_mixing_old
     489              :       TYPE(energy_correction_type), POINTER              :: ec_env
     490              :       TYPE(qs_diis_buffer_type_sparse), POINTER          :: diis_buffer
     491              :       TYPE(transport_env_type), POINTER                  :: transport_env
     492              : 
     493          882 :       CALL timeset(routineN, handle)
     494              : 
     495              :       ! get a useful output_unit
     496          882 :       logger => cp_get_default_logger()
     497          882 :       IF (logger%para_env%is_source()) THEN
     498          441 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     499              :       ELSE
     500          441 :          unit_nr = -1
     501              :       END IF
     502              : 
     503          882 :       nspin = ls_scf_env%nspins
     504              : 
     505              :       ! old quantities, useful for mixing
     506         5332 :       ALLOCATE (matrix_mixing_old(nspin), matrix_ks_deviation(nspin))
     507         1784 :       DO ispin = 1, nspin
     508          902 :          CALL dbcsr_create(matrix_mixing_old(ispin), template=ls_scf_env%matrix_ks(ispin))
     509              : 
     510          902 :          CALL dbcsr_create(matrix_ks_deviation(ispin), template=ls_scf_env%matrix_ks(ispin))
     511         1784 :          CALL dbcsr_set(matrix_ks_deviation(ispin), 0.0_dp)
     512              :       END DO
     513         2646 :       ls_scf_env%homo_spin(:) = 0.0_dp
     514         2646 :       ls_scf_env%lumo_spin(:) = 0.0_dp
     515              : 
     516          882 :       transm_scf_converged = .FALSE.
     517          882 :       transm_maxscf_reached = .FALSE.
     518              : 
     519          882 :       energy_old = 0.0_dp
     520          882 :       IF (ls_scf_env%scf_history%istore > 0) energy_old = ls_scf_env%energy_init
     521          882 :       check_convergence = .TRUE.
     522          882 :       iscf = 0
     523          882 :       IF (ls_scf_env%ls_diis) THEN
     524            4 :          diis_step = .FALSE.
     525            4 :          eps_diis = ls_scf_env%eps_diis
     526            4 :          nmixing = ls_scf_env%nmixing
     527              :          NULLIFY (diis_buffer)
     528            4 :          ALLOCATE (diis_buffer)
     529              :          CALL qs_diis_b_create_sparse(diis_buffer, &
     530            4 :                                       nbuffer=ls_scf_env%max_diis)
     531            4 :          CALL qs_diis_b_clear_sparse(diis_buffer)
     532            4 :          CALL get_qs_env(qs_env, matrix_s=matrix_s)
     533              :       END IF
     534              : 
     535          882 :       CALL get_qs_env(qs_env, transport_env=transport_env, do_transport=do_transport)
     536              : 
     537              :       ! the real SCF loop
     538         2746 :       DO
     539              : 
     540              :          ! check on max SCF or timing/exit
     541         2746 :          CALL external_control(should_stop, "SCF", start_time=qs_env%start_time, target_time=qs_env%target_time)
     542         2746 :          IF (do_transport) THEN
     543            0 :             maxscf_reached = should_stop .OR. iscf >= ls_scf_env%max_scf
     544              :             ! one extra scf step for post-processing in transmission calculations
     545            0 :             IF (transport_env%params%method .EQ. transport_transmission) THEN
     546            0 :                IF (transm_maxscf_reached) THEN
     547            0 :                   IF (unit_nr > 0) WRITE (unit_nr, '(T2,A)') "SCF not converged! "
     548              :                   EXIT
     549              :                END IF
     550              :                transm_maxscf_reached = maxscf_reached
     551              :             ELSE
     552            0 :                IF (maxscf_reached) THEN
     553            0 :                   IF (unit_nr > 0) WRITE (unit_nr, '(T2,A)') "SCF not converged! "
     554              :                   EXIT
     555              :                END IF
     556              :             END IF
     557              :          ELSE
     558         2746 :             IF (should_stop .OR. iscf >= ls_scf_env%max_scf) THEN
     559           46 :                IF (unit_nr > 0) WRITE (unit_nr, '(T2,A)') "SCF not converged! "
     560              :                ! Skip Harris functional calculation if ground-state is NOT converged
     561           46 :                IF (qs_env%energy_correction) THEN
     562            0 :                   CALL get_qs_env(qs_env, ec_env=ec_env)
     563            0 :                   IF (ec_env%skip_ec) ec_env%do_skip = .TRUE.
     564              :                END IF
     565              :                EXIT
     566              :             END IF
     567              :          END IF
     568              : 
     569         2700 :          t1 = m_walltime()
     570         2700 :          iscf = iscf + 1
     571              : 
     572              :          ! first get a copy of the current KS matrix
     573         2700 :          CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
     574         5488 :          DO ispin = 1, nspin
     575              :             CALL matrix_qs_to_ls(ls_scf_env%matrix_ks(ispin), matrix_ks(ispin)%matrix, &
     576         2788 :                                  ls_scf_env%ls_mstruct, covariant=.TRUE.)
     577         2788 :             IF (ls_scf_env%has_s_preconditioner) THEN
     578              :                CALL apply_matrix_preconditioner(ls_scf_env%matrix_ks(ispin), "forward", &
     579         1302 :                                                 ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
     580              :             END IF
     581         5488 :             CALL dbcsr_filter(ls_scf_env%matrix_ks(ispin), ls_scf_env%eps_filter)
     582              :          END DO
     583              :          ! run curvy steps if required. Needs an idempotent DM (either perification or restart)
     584         2700 :          IF ((iscf > 1 .OR. ls_scf_env%scf_history%istore > 0) .AND. ls_scf_env%curvy_steps) THEN
     585           90 :             CALL dm_ls_curvy_optimization(ls_scf_env, energy_old, check_convergence)
     586              :          ELSE
     587              :             ! turn the KS matrix in a density matrix
     588         5296 :             DO ispin = 1, nspin
     589         2686 :                IF (nonscf) THEN
     590           66 :                   CALL dbcsr_copy(matrix_mixing_old(ispin), ls_scf_env%matrix_ks(ispin))
     591         2620 :                ELSEIF (ls_scf_env%do_rho_mixing) THEN
     592            0 :                   CALL dbcsr_copy(matrix_mixing_old(ispin), ls_scf_env%matrix_ks(ispin))
     593              :                ELSE
     594         2620 :                   IF (iscf == 1) THEN
     595              :                      ! initialize the mixing matrix with the current state if needed
     596          826 :                      CALL dbcsr_copy(matrix_mixing_old(ispin), ls_scf_env%matrix_ks(ispin))
     597              :                   ELSE
     598         1794 :                      IF (ls_scf_env%ls_diis) THEN ! ------- IF-DIIS+MIX--- START
     599           28 :                         IF (diis_step .AND. (iscf - 1) .GE. ls_scf_env%iter_ini_diis) THEN
     600           22 :                            IF (unit_nr > 0) THEN
     601              :                               WRITE (unit_nr, '(A61)') &
     602           11 :                                  '*************************************************************'
     603              :                               WRITE (unit_nr, '(A50,2(I3,A1),L1,A1)') &
     604           11 :                                  " Using DIIS mixed KS:  (iscf,INI_DIIS,DIIS_STEP)=(", &
     605           22 :                                  iscf, ",", ls_scf_env%iter_ini_diis, ",", diis_step, ")"
     606              :                               WRITE (unit_nr, '(A52)') &
     607           11 :                                  " KS_nw= DIIS-Linear-Combination-Previous KS matrices"
     608              :                               WRITE (unit_nr, '(61A)') &
     609           11 :                                  "*************************************************************"
     610              :                            END IF
     611              :                            CALL dbcsr_copy(matrix_mixing_old(ispin), & ! out
     612           22 :                                            ls_scf_env%matrix_ks(ispin)) ! in
     613              :                         ELSE
     614            6 :                            IF (unit_nr > 0) THEN
     615              :                               WRITE (unit_nr, '(A57)') &
     616            3 :                                  "*********************************************************"
     617              :                               WRITE (unit_nr, '(A23,F5.3,A25,I3)') &
     618            3 :                                  " Using MIXING_FRACTION=", ls_scf_env%mixing_fraction, &
     619            6 :                                  " to mix KS matrix:  iscf=", iscf
     620              :                               WRITE (unit_nr, '(A7,F5.3,A6,F5.3,A7)') &
     621            3 :                                  " KS_nw=", ls_scf_env%mixing_fraction, "*KS + ", &
     622            6 :                                  1.0_dp - ls_scf_env%mixing_fraction, "*KS_old"
     623              :                               WRITE (unit_nr, '(A57)') &
     624            3 :                                  "*********************************************************"
     625              :                            END IF
     626              :                            ! perform the mixing of ks matrices
     627              :                            CALL dbcsr_add(matrix_mixing_old(ispin), &
     628              :                                           ls_scf_env%matrix_ks(ispin), &
     629              :                                           1.0_dp - ls_scf_env%mixing_fraction, &
     630            6 :                                           ls_scf_env%mixing_fraction)
     631              :                         END IF
     632              :                      ELSE ! otherwise
     633         1766 :                         IF (unit_nr > 0) THEN
     634              :                            WRITE (unit_nr, '(A57)') &
     635          883 :                               "*********************************************************"
     636              :                            WRITE (unit_nr, '(A23,F5.3,A25,I3)') &
     637          883 :                               " Using MIXING_FRACTION=", ls_scf_env%mixing_fraction, &
     638         1766 :                               " to mix KS matrix:  iscf=", iscf
     639              :                            WRITE (unit_nr, '(A7,F5.3,A6,F5.3,A7)') &
     640          883 :                               " KS_nw=", ls_scf_env%mixing_fraction, "*KS + ", &
     641         1766 :                               1.0_dp - ls_scf_env%mixing_fraction, "*KS_old"
     642              :                            WRITE (unit_nr, '(A57)') &
     643          883 :                               "*********************************************************"
     644              :                         END IF
     645              :                         ! perform the mixing of ks matrices
     646              :                         CALL dbcsr_add(matrix_mixing_old(ispin), &
     647              :                                        ls_scf_env%matrix_ks(ispin), &
     648              :                                        1.0_dp - ls_scf_env%mixing_fraction, &
     649         1766 :                                        ls_scf_env%mixing_fraction)
     650              :                      END IF ! ------- IF-DIIS+MIX--- END
     651              :                   END IF
     652              :                END IF
     653              : 
     654              :                ! compute the density matrix that matches it
     655              :                ! we need the proper number of states
     656         2686 :                nelectron_spin_real = ls_scf_env%nelectron_spin(ispin)
     657         2686 :                IF (ls_scf_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
     658              : 
     659         2686 :                IF (do_transport) THEN
     660            0 :                   IF (ls_scf_env%has_s_preconditioner) &
     661            0 :                      CPABORT("NOT YET IMPLEMENTED with S preconditioner. ")
     662            0 :                   IF (ls_scf_env%ls_mstruct%cluster_type .NE. ls_cluster_atomic) &
     663            0 :                      CPABORT("NOT YET IMPLEMENTED with molecular clustering. ")
     664              : 
     665            0 :                   extra_scf = maxscf_reached .OR. scf_converged
     666              :                   ! get the current Kohn-Sham matrix (ks) and return matrix_p evaluated using an external C routine
     667              :                   CALL external_scf_method(transport_env, ls_scf_env%matrix_s, matrix_mixing_old(ispin), &
     668              :                                            ls_scf_env%matrix_p(ispin), nelectron_spin_real, ls_scf_env%natoms, &
     669            0 :                                            energy_diff, iscf, extra_scf)
     670              : 
     671              :                ELSE
     672         3594 :                   SELECT CASE (ls_scf_env%purification_method)
     673              :                   CASE (ls_scf_sign)
     674              :                      CALL density_matrix_sign(ls_scf_env%matrix_p(ispin), ls_scf_env%mu_spin(ispin), ls_scf_env%fixed_mu, &
     675              :                                               ls_scf_env%sign_method, ls_scf_env%sign_order, matrix_mixing_old(ispin), &
     676              :                                               ls_scf_env%matrix_s, ls_scf_env%matrix_s_inv, nelectron_spin_real, &
     677              :                                               ls_scf_env%eps_filter, ls_scf_env%sign_symmetric, ls_scf_env%submatrix_sign_method, &
     678          908 :                                               ls_scf_env%matrix_s_sqrt_inv)
     679              :                   CASE (ls_scf_tc2)
     680              :                      CALL density_matrix_tc2(ls_scf_env%matrix_p(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s_sqrt_inv, &
     681              :                                              nelectron_spin_real, ls_scf_env%eps_filter, ls_scf_env%homo_spin(ispin), &
     682              :                                              ls_scf_env%lumo_spin(ispin), non_monotonic=ls_scf_env%non_monotonic, &
     683              :                                              eps_lanczos=ls_scf_env%eps_lanczos, max_iter_lanczos=ls_scf_env%max_iter_lanczos, &
     684          118 :                                              iounit=-1)
     685              :                   CASE (ls_scf_trs4)
     686              :                      CALL density_matrix_trs4(ls_scf_env%matrix_p(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s_sqrt_inv, &
     687              :                                               nelectron_spin_real, ls_scf_env%eps_filter, ls_scf_env%homo_spin(ispin), &
     688              :                                               ls_scf_env%lumo_spin(ispin), ls_scf_env%mu_spin(ispin), &
     689              :                                               dynamic_threshold=ls_scf_env%dynamic_threshold, &
     690              :                                               matrix_ks_deviation=matrix_ks_deviation(ispin), &
     691              :                                               eps_lanczos=ls_scf_env%eps_lanczos, max_iter_lanczos=ls_scf_env%max_iter_lanczos, &
     692         1660 :                                               iounit=-1)
     693              :                   CASE (ls_scf_pexsi)
     694            0 :                      IF (ls_scf_env%has_s_preconditioner) &
     695            0 :                         CPABORT("S preconditioning not implemented in combination with the PEXSI library. ")
     696            0 :                      IF (ls_scf_env%ls_mstruct%cluster_type .NE. ls_cluster_atomic) &
     697              :                         CALL cp_abort(__LOCATION__, &
     698            0 :                                       "Molecular clustering not implemented in combination with the PEXSI library. ")
     699              :                      CALL density_matrix_pexsi(ls_scf_env%pexsi, ls_scf_env%matrix_p(ispin), ls_scf_env%pexsi%matrix_w(ispin), &
     700              :                                                ls_scf_env%pexsi%kTS(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s, &
     701         2686 :                                                nelectron_spin_real, ls_scf_env%mu_spin(ispin), iscf, ispin)
     702              :                   END SELECT
     703              :                END IF
     704              : 
     705         2686 :                IF (ls_scf_env%has_s_preconditioner) THEN
     706              :                   CALL apply_matrix_preconditioner(ls_scf_env%matrix_p(ispin), "forward", &
     707         1302 :                                                    ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
     708              :                END IF
     709         2686 :                CALL dbcsr_filter(ls_scf_env%matrix_p(ispin), ls_scf_env%eps_filter)
     710              : 
     711         5296 :                IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 2.0_dp)
     712              : 
     713              :             END DO
     714              :          END IF
     715              : 
     716              :          ! compute the corresponding new energy KS matrix and new energy
     717         2700 :          IF (nonscf) THEN
     718           66 :             CALL ls_nonscf_energy(qs_env, ls_scf_env)
     719              :          ELSE
     720         2634 :             CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, energy_new, iscf)
     721              :          END IF
     722              : 
     723         2700 :          IF (ls_scf_env%do_pexsi) THEN
     724            0 :             CALL pexsi_to_qs(ls_scf_env, qs_env, kTS=ls_scf_env%pexsi%kTS)
     725              :          END IF
     726              : 
     727         2700 :          t2 = m_walltime()
     728         2700 :          IF (nonscf) THEN
     729           66 :             tdiag = t2 - t1
     730           66 :             CALL qs_nonscf_print_summary(qs_env, tdiag, ls_scf_env%nelectron_total, unit_nr)
     731           66 :             EXIT
     732              :          ELSE
     733              :             ! report current SCF loop
     734         2634 :             energy_diff = energy_new - energy_old
     735         2634 :             energy_old = energy_new
     736         2634 :             IF (unit_nr > 0) THEN
     737         1317 :                WRITE (unit_nr, *)
     738         1317 :                WRITE (unit_nr, '(T2,A,I6,F20.9,F20.9,F12.6)') "SCF", iscf, energy_new, energy_diff, t2 - t1
     739         1317 :                WRITE (unit_nr, *)
     740         1317 :                CALL m_flush(unit_nr)
     741              :             END IF
     742              :          END IF
     743              : 
     744         2634 :          IF (do_transport) THEN
     745            0 :             scf_converged = check_convergence .AND. ABS(energy_diff) < ls_scf_env%eps_scf*ls_scf_env%nelectron_total
     746              :             ! one extra scf step for post-processing in transmission calculations
     747            0 :             IF (transport_env%params%method .EQ. transport_transmission) THEN
     748            0 :                IF (transm_scf_converged) EXIT
     749              :                transm_scf_converged = scf_converged
     750              :             ELSE
     751            0 :                IF (scf_converged) THEN
     752            0 :                   IF (unit_nr > 0) WRITE (unit_nr, '(/,T2,A,I5,A/)') "SCF run converged in ", iscf, " steps."
     753              :                   EXIT
     754              :                END IF
     755              :             END IF
     756              :          ELSE
     757              :             ! exit criterion on the energy only for the time being
     758         2634 :             IF (check_convergence .AND. ABS(energy_diff) < ls_scf_env%eps_scf*ls_scf_env%nelectron_total) THEN
     759          770 :                IF (unit_nr > 0) WRITE (unit_nr, '(/,T2,A,I5,A/)') "SCF run converged in ", iscf, " steps."
     760              :                ! Skip Harris functional calculation if ground-state is NOT converged
     761          770 :                IF (qs_env%energy_correction) THEN
     762           20 :                   CALL get_qs_env(qs_env, ec_env=ec_env)
     763           20 :                   IF (ec_env%skip_ec) ec_env%do_skip = .FALSE.
     764              :                END IF
     765              :                EXIT
     766              :             END IF
     767              :          END IF
     768              : 
     769         1864 :          IF (ls_scf_env%ls_diis) THEN
     770              : ! diis_buffer, buffer with 1) Kohn-Sham history matrix,
     771              : !                          2) KS error history matrix (f=KPS-SPK),
     772              : !                          3) B matrix (for finding DIIS weighting coefficients)
     773              :             CALL qs_diis_b_step_4lscf(diis_buffer, qs_env, ls_scf_env, unit_nr, &
     774              :                                       iscf, diis_step, eps_diis, nmixing, matrix_s(1)%matrix, &
     775           18 :                                       ls_scf_env%eps_filter)
     776              :          END IF
     777              : 
     778         1864 :          IF (ls_scf_env%do_pexsi) THEN
     779              :             CALL pexsi_set_convergence_tolerance(ls_scf_env%pexsi, energy_diff, &
     780              :                                                  ls_scf_env%eps_scf*ls_scf_env%nelectron_total, &
     781              :                                                  ! initialize in second scf step of first SCF cycle:
     782              :                                                  (iscf .EQ. 2) .AND. (ls_scf_env%scf_history%istore .EQ. 0), &
     783            0 :                                                  check_convergence)
     784              :          END IF
     785              : 
     786              :       END DO
     787              : 
     788              :       ! free storage
     789          882 :       IF (ls_scf_env%ls_diis) THEN
     790            4 :          CALL qs_diis_b_release_sparse(diis_buffer)
     791            4 :          DEALLOCATE (diis_buffer)
     792              :       END IF
     793         1784 :       DO ispin = 1, nspin
     794          902 :          CALL dbcsr_release(matrix_mixing_old(ispin))
     795         1784 :          CALL dbcsr_release(matrix_ks_deviation(ispin))
     796              :       END DO
     797          882 :       DEALLOCATE (matrix_mixing_old, matrix_ks_deviation)
     798              : 
     799          882 :       CALL timestop(handle)
     800              : 
     801          882 :    END SUBROUTINE ls_scf_main
     802              : 
     803              : ! **************************************************************************************************
     804              : !> \brief after SCF we have a density matrix, and the self consistent KS matrix
     805              : !>        analyze its properties.
     806              : !> \param qs_env ...
     807              : !> \param ls_scf_env ...
     808              : !> \par History
     809              : !>       2010.10 created [Joost VandeVondele]
     810              : !> \author Joost VandeVondele
     811              : ! **************************************************************************************************
     812          882 :    SUBROUTINE ls_scf_post(qs_env, ls_scf_env)
     813              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     814              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
     815              : 
     816              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ls_scf_post'
     817              : 
     818              :       INTEGER                                            :: handle, ispin, unit_nr
     819              :       REAL(KIND=dp)                                      :: occ
     820              :       TYPE(cp_logger_type), POINTER                      :: logger
     821          882 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_w
     822              :       TYPE(dft_control_type), POINTER                    :: dft_control
     823              : 
     824          882 :       CALL timeset(routineN, handle)
     825              : 
     826          882 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     827              : 
     828              :       ! get a useful output_unit
     829          882 :       logger => cp_get_default_logger()
     830          882 :       IF (logger%para_env%is_source()) THEN
     831          441 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     832              :       ELSE
     833          441 :          unit_nr = -1
     834              :       END IF
     835              : 
     836              :       ! store the matrix for a next scf run
     837          882 :       IF (.NOT. ls_scf_env%do_pao) &
     838          370 :          CALL ls_scf_store_result(ls_scf_env)
     839              : 
     840              :       ! write homo and lumo energy and occupation (if not already part of the output)
     841          882 :       IF (ls_scf_env%curvy_steps) THEN
     842           18 :          CALL post_scf_homo_lumo(ls_scf_env)
     843              : 
     844              :          ! always report P occ
     845           18 :          IF (unit_nr > 0) WRITE (unit_nr, *) ""
     846           38 :          DO ispin = 1, ls_scf_env%nspins
     847           20 :             occ = dbcsr_get_occupation(ls_scf_env%matrix_p(ispin))
     848           38 :             IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,F20.12)') "Density matrix (P) occupation ", occ
     849              :          END DO
     850              :       END IF
     851              : 
     852              :       ! compute the matrix_w if associated
     853          882 :       IF (ls_scf_env%calculate_forces) THEN
     854          170 :          CALL get_qs_env(qs_env, matrix_w=matrix_w)
     855          170 :          CPASSERT(ASSOCIATED(matrix_w))
     856          170 :          IF (ls_scf_env%do_pexsi) THEN
     857            0 :             CALL pexsi_to_qs(ls_scf_env, qs_env, matrix_w=ls_scf_env%pexsi%matrix_w)
     858              :          ELSE
     859          170 :             CALL calculate_w_matrix_ls(matrix_w, ls_scf_env)
     860              :          END IF
     861              :       END IF
     862              : 
     863              :       ! compute properties
     864              : 
     865          882 :       IF (ls_scf_env%perform_mu_scan) CALL post_scf_mu_scan(ls_scf_env)
     866              : 
     867          882 :       IF (ls_scf_env%report_all_sparsities) CALL post_scf_sparsities(ls_scf_env)
     868              : 
     869          882 :       IF (dft_control%qs_control%dftb) THEN
     870           44 :          CALL scf_post_calculation_tb(qs_env, "DFTB", .TRUE.)
     871          838 :       ELSE IF (dft_control%qs_control%xtb) THEN
     872           86 :          CALL scf_post_calculation_tb(qs_env, "xTB", .TRUE.)
     873              :       ELSE
     874          752 :          CALL write_mo_free_results(qs_env)
     875              :       END IF
     876              : 
     877          882 :       IF (ls_scf_env%chebyshev%compute_chebyshev) CALL compute_chebyshev(qs_env, ls_scf_env)
     878              : 
     879          882 :       IF (.TRUE.) CALL post_scf_experiment()
     880              : 
     881          882 :       IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
     882              :          !
     883              :       ELSE
     884          752 :          CALL qs_scf_post_moments(qs_env%input, logger, qs_env, unit_nr)
     885              :       END IF
     886              : 
     887              :       ! clean up used data
     888              : 
     889          882 :       CALL dbcsr_release(ls_scf_env%matrix_s)
     890          882 :       CALL deallocate_curvy_data(ls_scf_env%curvy_data)
     891              : 
     892          882 :       IF (ls_scf_env%has_s_preconditioner) THEN
     893          326 :          CALL dbcsr_release(ls_scf_env%matrix_bs_sqrt)
     894          326 :          CALL dbcsr_release(ls_scf_env%matrix_bs_sqrt_inv)
     895              :       END IF
     896              : 
     897          882 :       IF (ls_scf_env%needs_s_inv) THEN
     898          880 :          CALL dbcsr_release(ls_scf_env%matrix_s_inv)
     899              :       END IF
     900              : 
     901          882 :       IF (ls_scf_env%use_s_sqrt) THEN
     902          878 :          CALL dbcsr_release(ls_scf_env%matrix_s_sqrt)
     903          878 :          CALL dbcsr_release(ls_scf_env%matrix_s_sqrt_inv)
     904              :       END IF
     905              : 
     906         1784 :       DO ispin = 1, SIZE(ls_scf_env%matrix_ks)
     907         1784 :          CALL dbcsr_release(ls_scf_env%matrix_ks(ispin))
     908              :       END DO
     909          882 :       DEALLOCATE (ls_scf_env%matrix_ks)
     910              : 
     911          882 :       IF (ls_scf_env%do_pexsi) &
     912            0 :          CALL pexsi_finalize_scf(ls_scf_env%pexsi, ls_scf_env%mu_spin)
     913              : 
     914          882 :       CALL timestop(handle)
     915              : 
     916          882 :    END SUBROUTINE ls_scf_post
     917              : 
     918              : ! **************************************************************************************************
     919              : !> \brief Compute the HOMO LUMO energies post SCF
     920              : !> \param ls_scf_env ...
     921              : !> \par History
     922              : !>       2013.06 created [Joost VandeVondele]
     923              : !> \author Joost VandeVondele
     924              : ! **************************************************************************************************
     925           18 :    SUBROUTINE post_scf_homo_lumo(ls_scf_env)
     926              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
     927              : 
     928              :       CHARACTER(len=*), PARAMETER :: routineN = 'post_scf_homo_lumo'
     929              : 
     930              :       INTEGER                                            :: handle, ispin, nspin, unit_nr
     931              :       LOGICAL                                            :: converged
     932              :       REAL(KIND=dp)                                      :: eps_max, eps_min, homo, lumo
     933              :       TYPE(cp_logger_type), POINTER                      :: logger
     934              :       TYPE(dbcsr_type)                                   :: matrix_k, matrix_p, matrix_tmp
     935              : 
     936           18 :       CALL timeset(routineN, handle)
     937              : 
     938              :       ! get a useful output_unit
     939           18 :       logger => cp_get_default_logger()
     940           18 :       IF (logger%para_env%is_source()) THEN
     941            9 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     942              :       ELSE
     943            9 :          unit_nr = -1
     944              :       END IF
     945              : 
     946           18 :       IF (unit_nr > 0) WRITE (unit_nr, '(T2,A)') ""
     947              : 
     948              :       ! TODO: remove these limitations
     949           18 :       CPASSERT(.NOT. ls_scf_env%has_s_preconditioner)
     950           18 :       CPASSERT(ls_scf_env%use_s_sqrt)
     951              : 
     952           18 :       nspin = ls_scf_env%nspins
     953              : 
     954           18 :       CALL dbcsr_create(matrix_p, template=ls_scf_env%matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
     955              : 
     956           18 :       CALL dbcsr_create(matrix_k, template=ls_scf_env%matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
     957              : 
     958           18 :       CALL dbcsr_create(matrix_tmp, template=ls_scf_env%matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
     959              : 
     960           38 :       DO ispin = 1, nspin
     961              :          ! ortho basis ks
     962              :          CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_ks(ispin), &
     963           20 :                              0.0_dp, matrix_tmp, filter_eps=ls_scf_env%eps_filter)
     964              :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, ls_scf_env%matrix_s_sqrt_inv, &
     965           20 :                              0.0_dp, matrix_k, filter_eps=ls_scf_env%eps_filter)
     966              : 
     967              :          ! extremal eigenvalues ks
     968              :          CALL arnoldi_extremal(matrix_k, eps_max, eps_min, max_iter=ls_scf_env%max_iter_lanczos, &
     969           20 :                                threshold=ls_scf_env%eps_lanczos, converged=converged)
     970              : 
     971              :          ! ortho basis p
     972              :          CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt, ls_scf_env%matrix_p(ispin), &
     973           20 :                              0.0_dp, matrix_tmp, filter_eps=ls_scf_env%eps_filter)
     974              :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, ls_scf_env%matrix_s_sqrt, &
     975           20 :                              0.0_dp, matrix_p, filter_eps=ls_scf_env%eps_filter)
     976           20 :          IF (nspin == 1) CALL dbcsr_scale(matrix_p, 0.5_dp)
     977              : 
     978              :          ! go compute homo lumo
     979              :          CALL compute_homo_lumo(matrix_k, matrix_p, eps_min, eps_max, ls_scf_env%eps_filter, &
     980           58 :                                 ls_scf_env%max_iter_lanczos, ls_scf_env%eps_lanczos, homo, lumo, unit_nr)
     981              : 
     982              :       END DO
     983              : 
     984           18 :       CALL dbcsr_release(matrix_p)
     985           18 :       CALL dbcsr_release(matrix_k)
     986           18 :       CALL dbcsr_release(matrix_tmp)
     987              : 
     988           18 :       CALL timestop(handle)
     989              : 
     990           18 :    END SUBROUTINE post_scf_homo_lumo
     991              : 
     992              : ! **************************************************************************************************
     993              : !> \brief Compute the density matrix for various values of the chemical potential
     994              : !> \param ls_scf_env ...
     995              : !> \par History
     996              : !>       2010.10 created [Joost VandeVondele]
     997              : !> \author Joost VandeVondele
     998              : ! **************************************************************************************************
     999            2 :    SUBROUTINE post_scf_mu_scan(ls_scf_env)
    1000              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
    1001              : 
    1002              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'post_scf_mu_scan'
    1003              : 
    1004              :       INTEGER                                            :: handle, imu, ispin, nelectron_spin_real, &
    1005              :                                                             nmu, nspin, unit_nr
    1006              :       REAL(KIND=dp)                                      :: mu, t1, t2, trace
    1007              :       TYPE(cp_logger_type), POINTER                      :: logger
    1008              :       TYPE(dbcsr_type)                                   :: matrix_p
    1009              : 
    1010            2 :       CALL timeset(routineN, handle)
    1011              : 
    1012              :       ! get a useful output_unit
    1013            2 :       logger => cp_get_default_logger()
    1014            2 :       IF (logger%para_env%is_source()) THEN
    1015            1 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
    1016              :       ELSE
    1017              :          unit_nr = -1
    1018              :       END IF
    1019              : 
    1020            2 :       nspin = ls_scf_env%nspins
    1021              : 
    1022            2 :       CALL dbcsr_create(matrix_p, template=ls_scf_env%matrix_p(1))
    1023              : 
    1024            2 :       nmu = 10
    1025           24 :       DO imu = 0, nmu
    1026              : 
    1027           22 :          t1 = m_walltime()
    1028              : 
    1029           22 :          mu = -0.4_dp + imu*(0.1_dp + 0.4_dp)/nmu
    1030              : 
    1031           22 :          IF (unit_nr > 0) WRITE (unit_nr, *) "------- starting with mu ", mu
    1032              : 
    1033           44 :          DO ispin = 1, nspin
    1034              :             ! we need the proper number of states
    1035           22 :             nelectron_spin_real = ls_scf_env%nelectron_spin(ispin)
    1036           22 :             IF (ls_scf_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
    1037              : 
    1038              :             CALL density_matrix_sign_fixed_mu(matrix_p, trace, mu, ls_scf_env%sign_method, &
    1039              :                                               ls_scf_env%sign_order, ls_scf_env%matrix_ks(ispin), &
    1040              :                                               ls_scf_env%matrix_s, ls_scf_env%matrix_s_inv, &
    1041              :                                               ls_scf_env%eps_filter, ls_scf_env%sign_symmetric, &
    1042           44 :                                               ls_scf_env%submatrix_sign_method, ls_scf_env%matrix_s_sqrt_inv)
    1043              :          END DO
    1044              : 
    1045           22 :          t2 = m_walltime()
    1046              : 
    1047           24 :          IF (unit_nr > 0) WRITE (unit_nr, *) " obtained ", mu, trace, t2 - t1
    1048              : 
    1049              :       END DO
    1050              : 
    1051            2 :       CALL dbcsr_release(matrix_p)
    1052              : 
    1053            2 :       CALL timestop(handle)
    1054              : 
    1055            2 :    END SUBROUTINE post_scf_mu_scan
    1056              : 
    1057              : ! **************************************************************************************************
    1058              : !> \brief Report on the sparsities of various interesting matrices.
    1059              : !>
    1060              : !> \param ls_scf_env ...
    1061              : !> \par History
    1062              : !>       2010.10 created [Joost VandeVondele]
    1063              : !> \author Joost VandeVondele
    1064              : ! **************************************************************************************************
    1065          180 :    SUBROUTINE post_scf_sparsities(ls_scf_env)
    1066              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
    1067              : 
    1068              :       CHARACTER(len=*), PARAMETER :: routineN = 'post_scf_sparsities'
    1069              : 
    1070              :       CHARACTER(LEN=default_string_length)               :: title
    1071              :       INTEGER                                            :: handle, ispin, nspin, unit_nr
    1072              :       TYPE(cp_logger_type), POINTER                      :: logger
    1073              :       TYPE(dbcsr_type)                                   :: matrix_tmp1, matrix_tmp2
    1074              : 
    1075          180 :       CALL timeset(routineN, handle)
    1076              : 
    1077              :       ! get a useful output_unit
    1078          180 :       logger => cp_get_default_logger()
    1079          180 :       IF (logger%para_env%is_source()) THEN
    1080           90 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
    1081              :       ELSE
    1082           90 :          unit_nr = -1
    1083              :       END IF
    1084              : 
    1085          180 :       nspin = ls_scf_env%nspins
    1086              : 
    1087          180 :       IF (unit_nr > 0) THEN
    1088           90 :          WRITE (unit_nr, '()')
    1089           90 :          WRITE (unit_nr, '(T2,A,E17.3)') "Sparsity reports for eps_filter: ", ls_scf_env%eps_filter
    1090           90 :          WRITE (unit_nr, '()')
    1091              :       END IF
    1092              : 
    1093              :       CALL report_matrix_sparsity(ls_scf_env%matrix_s, unit_nr, "overlap matrix (S)", &
    1094          180 :                                   ls_scf_env%eps_filter)
    1095              : 
    1096          368 :       DO ispin = 1, nspin
    1097          188 :          WRITE (title, '(A,I3)') "Kohn-Sham matrix (H) for spin ", ispin
    1098              :          CALL report_matrix_sparsity(ls_scf_env%matrix_ks(ispin), unit_nr, title, &
    1099          368 :                                      ls_scf_env%eps_filter)
    1100              :       END DO
    1101              : 
    1102          180 :       CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
    1103          180 :       CALL dbcsr_create(matrix_tmp2, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
    1104              : 
    1105          368 :       DO ispin = 1, nspin
    1106          188 :          WRITE (title, '(A,I3)') "Density matrix (P) for spin ", ispin
    1107              :          CALL report_matrix_sparsity(ls_scf_env%matrix_p(ispin), unit_nr, title, &
    1108          188 :                                      ls_scf_env%eps_filter)
    1109              : 
    1110              :          CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s, ls_scf_env%matrix_p(ispin), &
    1111          188 :                              0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
    1112              : 
    1113          188 :          WRITE (title, '(A,I3,A)') "S * P(", ispin, ")"
    1114          188 :          CALL report_matrix_sparsity(matrix_tmp1, unit_nr, title, ls_scf_env%eps_filter)
    1115              : 
    1116              :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s, &
    1117          188 :                              0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
    1118          188 :          WRITE (title, '(A,I3,A)') "S * P(", ispin, ") * S"
    1119          368 :          CALL report_matrix_sparsity(matrix_tmp2, unit_nr, title, ls_scf_env%eps_filter)
    1120              :       END DO
    1121              : 
    1122          180 :       IF (ls_scf_env%needs_s_inv) THEN
    1123              :          CALL report_matrix_sparsity(ls_scf_env%matrix_s_inv, unit_nr, "inv(S)", &
    1124          180 :                                      ls_scf_env%eps_filter)
    1125          368 :          DO ispin = 1, nspin
    1126              :             CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_inv, ls_scf_env%matrix_ks(ispin), &
    1127          188 :                                 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
    1128              : 
    1129          188 :             WRITE (title, '(A,I3,A)') "inv(S) * H(", ispin, ")"
    1130          368 :             CALL report_matrix_sparsity(matrix_tmp1, unit_nr, title, ls_scf_env%eps_filter)
    1131              :          END DO
    1132              :       END IF
    1133              : 
    1134          180 :       IF (ls_scf_env%use_s_sqrt) THEN
    1135              : 
    1136              :          CALL report_matrix_sparsity(ls_scf_env%matrix_s_sqrt, unit_nr, "sqrt(S)", &
    1137          178 :                                      ls_scf_env%eps_filter)
    1138              :          CALL report_matrix_sparsity(ls_scf_env%matrix_s_sqrt_inv, unit_nr, "inv(sqrt(S))", &
    1139          178 :                                      ls_scf_env%eps_filter)
    1140              : 
    1141          364 :          DO ispin = 1, nspin
    1142              :             CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_ks(ispin), &
    1143          186 :                                 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
    1144              :             CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt_inv, &
    1145          186 :                                 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
    1146          186 :             WRITE (title, '(A,I3,A)') "inv(sqrt(S)) * H(", ispin, ") * inv(sqrt(S))"
    1147          364 :             CALL report_matrix_sparsity(matrix_tmp2, unit_nr, title, ls_scf_env%eps_filter)
    1148              :          END DO
    1149              : 
    1150          364 :          DO ispin = 1, nspin
    1151              :             CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt, ls_scf_env%matrix_p(ispin), &
    1152          186 :                                 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
    1153              :             CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt, &
    1154          186 :                                 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
    1155          186 :             WRITE (title, '(A,I3,A)') "sqrt(S) * P(", ispin, ") * sqrt(S)"
    1156          364 :             CALL report_matrix_sparsity(matrix_tmp2, unit_nr, title, ls_scf_env%eps_filter)
    1157              :          END DO
    1158              : 
    1159              :       END IF
    1160              : 
    1161          180 :       CALL dbcsr_release(matrix_tmp1)
    1162          180 :       CALL dbcsr_release(matrix_tmp2)
    1163              : 
    1164          180 :       CALL timestop(handle)
    1165              : 
    1166          180 :    END SUBROUTINE post_scf_sparsities
    1167              : 
    1168              : ! **************************************************************************************************
    1169              : !> \brief Helper routine to report on the sparsity of a single matrix,
    1170              : !>        for several filtering values
    1171              : !> \param matrix ...
    1172              : !> \param unit_nr ...
    1173              : !> \param title ...
    1174              : !> \param eps ...
    1175              : !> \par History
    1176              : !>       2010.10 created [Joost VandeVondele]
    1177              : !> \author Joost VandeVondele
    1178              : ! **************************************************************************************************
    1179         2028 :    SUBROUTINE report_matrix_sparsity(matrix, unit_nr, title, eps)
    1180              :       TYPE(dbcsr_type)                                   :: matrix
    1181              :       INTEGER                                            :: unit_nr
    1182              :       CHARACTER(LEN=*)                                   :: title
    1183              :       REAL(KIND=dp)                                      :: eps
    1184              : 
    1185              :       CHARACTER(len=*), PARAMETER :: routineN = 'report_matrix_sparsity'
    1186              : 
    1187              :       INTEGER                                            :: handle
    1188              :       REAL(KIND=dp)                                      :: eps_local, occ
    1189              :       TYPE(dbcsr_type)                                   :: matrix_tmp
    1190              : 
    1191         2028 :       CALL timeset(routineN, handle)
    1192         2028 :       CALL dbcsr_create(matrix_tmp, template=matrix, name=TRIM(title))
    1193         2028 :       CALL dbcsr_copy(matrix_tmp, matrix, name=TRIM(title))
    1194              : 
    1195         2028 :       IF (unit_nr > 0) THEN
    1196         1014 :          WRITE (unit_nr, '(T2,A)') "Sparsity for : "//TRIM(title)
    1197              :       END IF
    1198              : 
    1199         2028 :       eps_local = MAX(eps, 10E-14_dp)
    1200        14796 :       DO
    1201        16824 :          IF (eps_local > 1.1_dp) EXIT
    1202        14796 :          CALL dbcsr_filter(matrix_tmp, eps_local)
    1203        14796 :          occ = dbcsr_get_occupation(matrix_tmp)
    1204        14796 :          IF (unit_nr > 0) WRITE (unit_nr, '(T2,F16.12,A3,F16.12)') eps_local, " : ", occ
    1205        14796 :          eps_local = eps_local*10
    1206              :       END DO
    1207              : 
    1208         2028 :       CALL dbcsr_release(matrix_tmp)
    1209              : 
    1210         2028 :       CALL timestop(handle)
    1211              : 
    1212         2028 :    END SUBROUTINE report_matrix_sparsity
    1213              : 
    1214              : ! **************************************************************************************************
    1215              : !> \brief Compute matrix_w as needed for the forces
    1216              : !> \param matrix_w ...
    1217              : !> \param ls_scf_env ...
    1218              : !> \par History
    1219              : !>       2010.11 created [Joost VandeVondele]
    1220              : !> \author Joost VandeVondele
    1221              : ! **************************************************************************************************
    1222          200 :    SUBROUTINE calculate_w_matrix_ls(matrix_w, ls_scf_env)
    1223              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_w
    1224              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
    1225              : 
    1226              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_ls'
    1227              : 
    1228              :       INTEGER                                            :: handle, ispin
    1229              :       REAL(KIND=dp)                                      :: scaling
    1230              :       TYPE(dbcsr_type)                                   :: matrix_tmp1, matrix_tmp2, matrix_tmp3
    1231              : 
    1232          200 :       CALL timeset(routineN, handle)
    1233              : 
    1234          200 :       CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
    1235          200 :       CALL dbcsr_create(matrix_tmp2, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
    1236          200 :       CALL dbcsr_create(matrix_tmp3, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
    1237              : 
    1238          200 :       IF (ls_scf_env%nspins == 1) THEN
    1239          192 :          scaling = 0.5_dp
    1240              :       ELSE
    1241            8 :          scaling = 1.0_dp
    1242              :       END IF
    1243              : 
    1244          408 :       DO ispin = 1, ls_scf_env%nspins
    1245              : 
    1246          208 :          CALL dbcsr_copy(matrix_tmp3, ls_scf_env%matrix_ks(ispin))
    1247          208 :          IF (ls_scf_env%has_s_preconditioner) THEN
    1248              :             CALL apply_matrix_preconditioner(matrix_tmp3, "backward", &
    1249          136 :                                              ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
    1250              :          END IF
    1251          208 :          CALL dbcsr_filter(matrix_tmp3, ls_scf_env%eps_filter)
    1252              : 
    1253              :          CALL dbcsr_multiply("N", "N", scaling, ls_scf_env%matrix_p(ispin), matrix_tmp3, &
    1254          208 :                              0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
    1255              :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_p(ispin), &
    1256          208 :                              0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
    1257          408 :          CALL matrix_ls_to_qs(matrix_w(ispin)%matrix, matrix_tmp2, ls_scf_env%ls_mstruct, covariant=.FALSE.)
    1258              :       END DO
    1259              : 
    1260          200 :       CALL dbcsr_release(matrix_tmp1)
    1261          200 :       CALL dbcsr_release(matrix_tmp2)
    1262          200 :       CALL dbcsr_release(matrix_tmp3)
    1263              : 
    1264          200 :       CALL timestop(handle)
    1265              : 
    1266          200 :    END SUBROUTINE calculate_w_matrix_ls
    1267              : 
    1268              : ! **************************************************************************************************
    1269              : !> \brief a place for quick experiments
    1270              : !> \par History
    1271              : !>       2010.11 created [Joost VandeVondele]
    1272              : !> \author Joost VandeVondele
    1273              : ! **************************************************************************************************
    1274          882 :    SUBROUTINE post_scf_experiment()
    1275              : 
    1276              :       CHARACTER(len=*), PARAMETER :: routineN = 'post_scf_experiment'
    1277              : 
    1278              :       INTEGER                                            :: handle, unit_nr
    1279              :       TYPE(cp_logger_type), POINTER                      :: logger
    1280              : 
    1281          882 :       CALL timeset(routineN, handle)
    1282              : 
    1283              :       ! get a useful output_unit
    1284          882 :       logger => cp_get_default_logger()
    1285          882 :       IF (logger%para_env%is_source()) THEN
    1286          441 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
    1287              :       ELSE
    1288              :          unit_nr = -1
    1289              :       END IF
    1290              : 
    1291          882 :       CALL timestop(handle)
    1292          882 :    END SUBROUTINE post_scf_experiment
    1293              : 
    1294              : END MODULE dm_ls_scf
        

Generated by: LCOV version 2.0-1