LCOV - code coverage report
Current view: top level - src - dm_ls_chebyshev.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:71c3ab0) Lines: 99.4 % 174 173
Test Date: 2026-07-25 06:35:44 Functions: 100.0 % 3 3

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Routines using linear scaling chebyshev methods
      10              : !> \par History
      11              : !>       2012.10 created [Jinwoong Cha]
      12              : !> \author Jinwoong Cha
      13              : ! **************************************************************************************************
      14              : MODULE dm_ls_chebyshev
      15              :    USE arnoldi_api,                     ONLY: arnoldi_extremal
      16              :    USE cp_dbcsr_api,                    ONLY: &
      17              :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_get_info, dbcsr_get_occupation, dbcsr_multiply, &
      18              :         dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
      19              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_add_on_diag,&
      20              :                                               dbcsr_frobenius_norm,&
      21              :                                               dbcsr_trace
      22              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      23              :                                               cp_logger_get_default_unit_nr,&
      24              :                                               cp_logger_type
      25              :    USE cp_output_handling,              ONLY: cp_p_file,&
      26              :                                               cp_print_key_finished_output,&
      27              :                                               cp_print_key_should_output,&
      28              :                                               cp_print_key_unit_nr
      29              :    USE dm_ls_scf_qs,                    ONLY: write_matrix_to_cube
      30              :    USE dm_ls_scf_types,                 ONLY: ls_scf_env_type
      31              :    USE input_section_types,             ONLY: section_get_ivals,&
      32              :                                               section_vals_val_get
      33              :    USE kinds,                           ONLY: default_string_length,&
      34              :                                               dp
      35              :    USE machine,                         ONLY: m_flush,&
      36              :                                               m_walltime
      37              :    USE mathconstants,                   ONLY: pi
      38              :    USE qs_environment_types,            ONLY: qs_environment_type
      39              : #include "./base/base_uses.f90"
      40              : 
      41              :    IMPLICIT NONE
      42              : 
      43              :    PRIVATE
      44              : 
      45              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_chebyshev'
      46              : 
      47              :    PUBLIC :: compute_chebyshev
      48              : 
      49              : CONTAINS
      50              : 
      51              : ! **************************************************************************************************
      52              : !> \brief compute chebyshev polynomials up to order n for a given value of x
      53              : !> \param value ...
      54              : !> \param x ...
      55              : !> \param n ...
      56              : !> \par History
      57              : !>       2012.11 created [Jinwoong Cha]
      58              : !> \author Jinwoong Cha
      59              : ! **************************************************************************************************
      60      3729300 :    SUBROUTINE chebyshev_poly(value, x, n)
      61              :       REAL(KIND=dp), INTENT(OUT)                         :: value
      62              :       REAL(KIND=dp), INTENT(IN)                          :: x
      63              :       INTEGER, INTENT(IN)                                :: n
      64              : 
      65              : !polynomial values
      66              : !number of chev polynomials
      67              : 
      68      3729300 :       value = COS((n - 1)*ACOS(x))
      69              : 
      70      3729300 :    END SUBROUTINE chebyshev_poly
      71              : 
      72              : ! **************************************************************************************************
      73              : !> \brief kernel for chebyshev polynomials expansion (Jackson kernel)
      74              : !> \param value ...
      75              : !> \param n ...
      76              : !> \param nc ...
      77              : !> \par History
      78              : !>       2012.11 created [Jinwoong Cha]
      79              : !> \author Jinwoong Cha
      80              : ! **************************************************************************************************
      81         3000 :    SUBROUTINE kernel(value, n, nc)
      82              :       REAL(KIND=dp), INTENT(OUT)                         :: value
      83              :       INTEGER, INTENT(IN)                                :: n, nc
      84              : 
      85              : !kernel at n
      86              : !n-1 order of chebyshev polynomials
      87              : !number of total chebyshev polynomials
      88              : !Kernel define
      89              : 
      90              :       value = 1.0_dp/(nc + 1.0_dp)*((nc - (n - 1) + 1.0_dp)* &
      91         3000 :                                     COS(pi*(n - 1)/(nc + 1.0_dp)) + SIN(pi*(n - 1)/(nc + 1.0_dp))*1.0_dp/TAN(pi/(nc + 1.0_dp)))
      92              : 
      93         3000 :    END SUBROUTINE kernel
      94              : 
      95              : ! **************************************************************************************************
      96              : !> \brief compute properties based on chebyshev expansion
      97              : !> \param qs_env ...
      98              : !> \param ls_scf_env ...
      99              : !> \par History
     100              : !>       2012.10 created [Jinwoong Cha]
     101              : !> \author Jinwoong Cha
     102              : ! **************************************************************************************************
     103            6 :    SUBROUTINE compute_chebyshev(qs_env, ls_scf_env)
     104              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     105              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
     106              : 
     107              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'compute_chebyshev'
     108              :       REAL(KIND=dp), PARAMETER                           :: scale_evals = 1.01_dp
     109              : 
     110              :       CHARACTER(LEN=30)                                  :: middle_name
     111              :       CHARACTER(LEN=default_string_length)               :: title
     112              :       INTEGER                                            :: handle, icheb, igrid, iinte, ispin, &
     113              :                                                             iwindow, n_gridpoint_dos, ncheb, &
     114              :                                                             ninte, Nrows, nwindow, unit_cube, &
     115              :                                                             unit_dos, unit_nr
     116              :       LOGICAL                                            :: converged, write_cubes
     117              :       REAL(KIND=dp) :: chev_T, chev_T_dos, dummy1, final, frob_matrix, initial, interval_a, &
     118              :          interval_b, max_ev, min_ev, occ, orbital_occ, summa, t1, t2
     119            6 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: chev_E, chev_Es_dos, dos, dummy2, ev1, &
     120            6 :                                                             ev2, kernel_g, mu, sev1, sev2, trace_dm
     121            6 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: aitchev_T, E_inte, gdensity, sqrt_vec
     122            6 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tmp_r
     123              :       TYPE(cp_logger_type), POINTER                      :: logger
     124              :       TYPE(dbcsr_type)                                   :: matrix_dummy1, matrix_F, matrix_tmp1, &
     125              :                                                             matrix_tmp2, matrix_tmp3
     126            6 :       TYPE(dbcsr_type), DIMENSION(:), POINTER            :: matrix_dummy2
     127              : 
     128            6 :       IF (.NOT. ls_scf_env%chebyshev%compute_chebyshev) RETURN
     129              : 
     130            6 :       CALL timeset(routineN, handle)
     131              : 
     132              :       ! get a useful output_unit
     133            6 :       logger => cp_get_default_logger()
     134            6 :       IF (logger%para_env%is_source()) THEN
     135            3 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     136              :       ELSE
     137            3 :          unit_nr = -1
     138              :       END IF
     139              : 
     140            6 :       ncheb = ls_scf_env%chebyshev%n_chebyshev
     141            6 :       ninte = 2*ncheb
     142            6 :       n_gridpoint_dos = ls_scf_env%chebyshev%n_gridpoint_dos
     143              : 
     144            6 :       write_cubes = BTEST(cp_print_key_should_output(logger%iter_info, ls_scf_env%chebyshev%print_key_cube), cp_p_file)
     145            6 :       IF (write_cubes) THEN
     146            2 :          IF (ASSOCIATED(ls_scf_env%chebyshev%min_energy)) DEALLOCATE (ls_scf_env%chebyshev%min_energy)
     147            2 :          CALL section_vals_val_get(ls_scf_env%chebyshev%print_key_cube, "MIN_ENERGY", r_vals=tmp_r)
     148            6 :          ALLOCATE (ls_scf_env%chebyshev%min_energy(SIZE(tmp_r)))
     149            8 :          ls_scf_env%chebyshev%min_energy = tmp_r
     150              : 
     151            2 :          IF (ASSOCIATED(ls_scf_env%chebyshev%max_energy)) DEALLOCATE (ls_scf_env%chebyshev%max_energy)
     152            2 :          CALL section_vals_val_get(ls_scf_env%chebyshev%print_key_cube, "MAX_ENERGY", r_vals=tmp_r)
     153            6 :          ALLOCATE (ls_scf_env%chebyshev%max_energy(SIZE(tmp_r)))
     154            8 :          ls_scf_env%chebyshev%max_energy = tmp_r
     155              : 
     156            2 :          nwindow = SIZE(ls_scf_env%chebyshev%min_energy)
     157              :       ELSE
     158              :          nwindow = 0
     159              :       END IF
     160              : 
     161           14 :       ALLOCATE (ev1(1:nwindow))
     162            8 :       ALLOCATE (ev2(1:nwindow))
     163            8 :       ALLOCATE (sev1(1:nwindow))
     164            8 :       ALLOCATE (sev2(1:nwindow))
     165            8 :       ALLOCATE (trace_dm(1:nwindow))
     166           20 :       ALLOCATE (matrix_dummy2(1:nwindow))
     167              : 
     168           12 :       DO iwindow = 1, nwindow
     169            6 :          ev1(iwindow) = ls_scf_env%chebyshev%min_energy(iwindow)
     170           12 :          ev2(iwindow) = ls_scf_env%chebyshev%max_energy(iwindow)
     171              :       END DO
     172              : 
     173            6 :       IF (unit_nr > 0) THEN
     174            3 :          WRITE (unit_nr, '()')
     175            3 :          WRITE (unit_nr, '(T2,A)') "STARTING CHEBYSHEV CALCULATION"
     176              :       END IF
     177              : 
     178              :       ! create 3 temporary matrices
     179            6 :       CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
     180            6 :       CALL dbcsr_create(matrix_tmp2, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
     181            6 :       CALL dbcsr_create(matrix_tmp3, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
     182            6 :       CALL dbcsr_create(matrix_F, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
     183            6 :       CALL dbcsr_create(matrix_dummy1, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
     184              : 
     185           12 :       DO iwindow = 1, nwindow
     186              :          CALL dbcsr_create(matrix_dummy2(iwindow), template=ls_scf_env%matrix_s, &
     187           12 :                            matrix_type=dbcsr_type_no_symmetry)
     188              :       END DO
     189              : 
     190           12 :       DO ispin = 1, SIZE(ls_scf_env%matrix_ks)
     191              :          ! create matrix_F=inv(sqrt(S))*H*inv(sqrt(S))
     192              :          CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_ks(ispin), &
     193            6 :                              0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
     194              :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt_inv, &
     195            6 :                              0.0_dp, matrix_F, filter_eps=ls_scf_env%eps_filter)
     196              : 
     197              :          ! find largest and smallest eigenvalues
     198              :          CALL arnoldi_extremal(matrix_F, max_ev, min_ev, converged=converged, max_iter=ls_scf_env%max_iter_lanczos, &
     199            6 :                                threshold=ls_scf_env%eps_lanczos) !Lanczos algorithm to calculate eigenvalue
     200            6 :          IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,2F16.8,A,L2)') &
     201            3 :             "smallest largest eigenvalue", min_ev, max_ev, " converged ", converged
     202            6 :          IF (nwindow > 0) THEN
     203            2 :             IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,1000F16.8)') "requested interval-min_energy", ev1(:)
     204            2 :             IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,1000F16.8)') "requested interval-max_energy", ev2(:)
     205              :          END IF
     206            6 :          interval_a = (max_ev - min_ev)*scale_evals/2
     207            6 :          interval_b = (max_ev + min_ev)/2
     208              : 
     209           12 :          sev1(:) = (ev1(:) - interval_b)/interval_a !scaled ev1 vector
     210           12 :          sev2(:) = (ev2(:) - interval_b)/interval_a !scaled ev2 vector
     211              : 
     212              :          !chebyshev domain,pi*sqrt(1-x^2) vector construction and chebyshev polynomials for integration (for g(E))
     213           20 :          ALLOCATE (E_inte(1:ninte + 1, 1:nwindow))
     214           14 :          ALLOCATE (sqrt_vec(1:ninte + 1, 1:nwindow))
     215              : 
     216           12 :          DO iwindow = 1, nwindow
     217         4818 :             DO iinte = 1, ninte + 1
     218         4806 :                E_inte(iinte, iwindow) = sev1(iwindow) + ((sev2(iwindow) - sev1(iwindow))/ninte)*(iinte - 1)
     219         4812 :                sqrt_vec(iinte, iwindow) = pi*SQRT(1.0_dp - E_inte(iinte, iwindow)*E_inte(iinte, iwindow))
     220              :             END DO
     221              :          END DO
     222              : 
     223              :          !integral.. (identical to the coefficient for g(E))
     224              : 
     225           20 :          ALLOCATE (aitchev_T(1:ncheb, 1:nwindow)) !after intergral. =>ainte
     226              : 
     227           12 :          DO iwindow = 1, nwindow
     228         2412 :             DO icheb = 1, ncheb
     229         2400 :                CALL chebyshev_poly(initial, E_inte(1, iwindow), icheb)
     230         2400 :                CALL chebyshev_poly(final, E_inte(1, iwindow), icheb)
     231         2400 :           summa = (sev2(iwindow) - sev1(iwindow))/(2.0_dp*ninte)*(initial/sqrt_vec(1, iwindow) + final/sqrt_vec(ninte + 1, iwindow))
     232      1920000 :                DO iinte = 2, ninte
     233      1917600 :                   CALL chebyshev_poly(chev_T, E_inte(iinte, iwindow), icheb)
     234      1920000 :                   summa = summa + ((sev2(iwindow) - sev1(iwindow))/ninte)*(chev_T/sqrt_vec(iinte, iwindow))
     235              :                END DO
     236         2400 :                aitchev_T(icheb, iwindow) = summa
     237         4806 :                summa = 0
     238              :             END DO
     239              :          END DO
     240              : 
     241              :          ! scale the matrix to get evals in the interval -1,1
     242            6 :          CALL dbcsr_add_on_diag(matrix_F, -interval_b)
     243            6 :          CALL dbcsr_scale(matrix_F, 1/interval_a)
     244              : 
     245              :          ! compute chebyshev matrix recursion
     246            6 :          CALL dbcsr_get_info(matrix=matrix_F, nfullrows_total=Nrows) !get information about a matrix
     247            6 :          CALL dbcsr_set(matrix_dummy1, 0.0_dp) !empty matrix creation(for density matrix)
     248              : 
     249           12 :          DO iwindow = 1, nwindow
     250           12 :             CALL dbcsr_set(matrix_dummy2(iwindow), 0.0_dp) !empty matrix creation(for density matrix)
     251              :          END DO
     252              : 
     253           18 :          ALLOCATE (mu(1:ncheb))
     254           12 :          ALLOCATE (kernel_g(1:ncheb))
     255            6 :          CALL kernel(kernel_g(1), 1, ncheb)
     256            6 :          CALL kernel(kernel_g(2), 2, ncheb)
     257              : 
     258            6 :          CALL dbcsr_set(matrix_tmp1, 0.0_dp) !matrix creation
     259            6 :          CALL dbcsr_add_on_diag(matrix_tmp1, 1.0_dp) !add a only number to diagonal elements
     260            6 :          CALL dbcsr_trace(matrix_tmp1, trace=mu(1))
     261            6 :          CALL dbcsr_copy(matrix_tmp2, matrix_F) !make matrix_tmp2 = matrix_F
     262            6 :          CALL dbcsr_trace(matrix_tmp2, trace=mu(2))
     263              : 
     264           12 :          DO iwindow = 1, nwindow
     265            6 :             CALL dbcsr_copy(matrix_dummy1, matrix_tmp1)
     266              :             !matrix_dummy2=
     267            6 :             CALL dbcsr_copy(matrix_dummy2(iwindow), matrix_tmp2)
     268              :             !first term of chebyshev poly(matrix)
     269            6 :             CALL dbcsr_scale(matrix_dummy1, kernel_g(1)*aitchev_T(1, iwindow))
     270              :             !second term of chebyshev poly(matrix)
     271            6 :             CALL dbcsr_scale(matrix_dummy2(iwindow), 2.0_dp*kernel_g(2)*aitchev_T(2, iwindow))
     272              : 
     273           12 :             CALL dbcsr_add(matrix_dummy2(iwindow), matrix_dummy1, 1.0_dp, 1.0_dp)
     274              :          END DO
     275              : 
     276         2994 :          DO icheb = 2, ncheb - 1
     277         2988 :             t1 = m_walltime()
     278              :             !matrix multiplication(Recursion)
     279              :             CALL dbcsr_multiply("N", "N", 2.0_dp, matrix_F, matrix_tmp2, &
     280         2988 :                                 -1.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
     281         2988 :             CALL dbcsr_copy(matrix_tmp3, matrix_tmp1)
     282         2988 :             CALL dbcsr_copy(matrix_tmp1, matrix_tmp2)
     283         2988 :             CALL dbcsr_copy(matrix_tmp2, matrix_tmp3)
     284              :             !icheb+1 th coefficient
     285         2988 :             CALL dbcsr_trace(matrix_tmp2, trace=mu(icheb + 1))
     286         2988 :             CALL kernel(kernel_g(icheb + 1), icheb + 1, ncheb)
     287              : 
     288         5376 :             DO iwindow = 1, nwindow
     289              : 
     290         2388 :                CALL dbcsr_copy(matrix_dummy1, matrix_tmp2)
     291              :                !second term of chebyshev poly(matrix)
     292         2388 :                CALL dbcsr_scale(matrix_dummy1, 2.0_dp*kernel_g(icheb + 1)*aitchev_T(icheb + 1, iwindow))
     293         2388 :                CALL dbcsr_add(matrix_dummy2(iwindow), matrix_dummy1, 1.0_dp, 1.0_dp)
     294         5376 :                CALL dbcsr_trace(matrix_dummy2(iwindow), trace=trace_dm(iwindow)) !icheb+1 th coefficient
     295              : 
     296              :             END DO
     297              : 
     298         2988 :             occ = dbcsr_get_occupation(matrix_tmp1)
     299         2988 :             t2 = m_walltime()
     300         2994 :             IF (unit_nr > 0 .AND. MOD(icheb, 20) == 0) THEN
     301           72 :                CALL m_flush(unit_nr)
     302           72 :                IF (nwindow > 0) THEN
     303              :                   WRITE (unit_nr, '(T2,A,I5,1X,A,1X,F8.3,1X,A,1X,F8.6,1X,A,1X,1000F16.8)') &
     304           19 :                      "Iter.", icheb, "time=", t2 - t1, "occ=", occ, "traces=", trace_dm(:)
     305              :                ELSE
     306              :                   WRITE (unit_nr, '(T2,A,I5,1X,A,1X,F8.3,1X,A,1X,F8.6)') &
     307           53 :                      "Iter.", icheb, "time=", t2 - t1, "occ=", occ
     308              :                END IF
     309              :             END IF
     310              :          END DO
     311              : 
     312           12 :          DO iwindow = 1, nwindow
     313            6 :             IF (SIZE(ls_scf_env%matrix_ks) == 1) THEN
     314            6 :                orbital_occ = 2.0_dp
     315              :             ELSE
     316            0 :                orbital_occ = 1.0_dp
     317              :             END IF
     318              :             CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, matrix_dummy2(iwindow), &
     319            6 :                                 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
     320              :             CALL dbcsr_multiply("N", "N", orbital_occ, matrix_tmp1, ls_scf_env%matrix_s_sqrt_inv, &
     321            6 :                                 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
     322            6 :             CALL dbcsr_copy(matrix_dummy2(iwindow), matrix_tmp2)
     323              : 
     324              :             ! look at the difference with the density matrix from the ls routines
     325            6 :             IF (.FALSE.) THEN
     326              :                CALL dbcsr_copy(matrix_tmp1, matrix_tmp2)
     327              :                CALL dbcsr_add(matrix_tmp1, ls_scf_env%matrix_p(ispin), 1.0_dp, -1.0_dp) !comparison
     328              :                frob_matrix = dbcsr_frobenius_norm(matrix_tmp1)
     329              :                IF (unit_nr > 0) WRITE (unit_nr, *) "Difference between Chebyshev DM and LS DM", frob_matrix
     330              :             END IF
     331              :          END DO
     332              : 
     333              :          write_cubes = BTEST(cp_print_key_should_output(logger%iter_info, &
     334            6 :                                                         ls_scf_env%chebyshev%print_key_cube), cp_p_file)
     335           24 :          IF (write_cubes) THEN
     336            8 :             DO iwindow = 1, nwindow
     337            6 :                WRITE (middle_name, "(A,I0)") "E_DENSITY_WINDOW_", iwindow
     338            6 :                WRITE (title, "(A,1X,F16.8,1X,A,1X,F16.8)") "Energy range : ", ev1(iwindow), "to", ev2(iwindow)
     339              :                unit_cube = cp_print_key_unit_nr(logger, ls_scf_env%chebyshev%print_key_cube, &
     340              :                                                 "", extension=".cube", & !added 01/22/2012
     341            6 :                                                 middle_name=TRIM(middle_name), log_filename=.FALSE.)
     342              :                CALL write_matrix_to_cube(qs_env, ls_scf_env, matrix_dummy2(iwindow), unit_cube, title, &
     343            6 :                                          section_get_ivals(ls_scf_env%chebyshev%print_key_cube, "STRIDE"))
     344            8 :                CALL cp_print_key_finished_output(unit_cube, logger, ls_scf_env%chebyshev%print_key_cube, "")
     345              :             END DO
     346              :          END IF
     347              : 
     348              :       END DO
     349              : 
     350              :       ! Chebyshev expansion with calculated coefficient
     351              :       ! grid construction and rescaling (by J)
     352              :       unit_dos = cp_print_key_unit_nr(logger, ls_scf_env%chebyshev%print_key_dos, "", extension=".xy", &
     353            6 :                                       middle_name="DOS", log_filename=.FALSE.)
     354              : 
     355            6 :       IF (unit_dos > 0) THEN
     356            9 :          ALLOCATE (dos(1:n_gridpoint_dos))
     357           10 :          ALLOCATE (gdensity(1:n_gridpoint_dos, 1:nwindow))
     358            6 :          ALLOCATE (chev_E(1:n_gridpoint_dos))
     359            6 :          ALLOCATE (chev_Es_dos(1:n_gridpoint_dos))
     360            4 :          ALLOCATE (dummy2(1:nwindow))
     361         3103 :          DO igrid = 1, n_gridpoint_dos
     362         3100 :             chev_E(igrid) = min_ev + (igrid - 1)*(max_ev - min_ev)/(n_gridpoint_dos - 1)
     363         3103 :             chev_Es_dos(igrid) = (chev_E(igrid) - interval_b)/interval_a
     364              :          END DO
     365         3103 :          DO igrid = 1, n_gridpoint_dos
     366         3100 :             dummy1 = 0.0_dp !summation of polynomials
     367         3100 :             dummy2(:) = 0.0_dp !summation of polynomials
     368      1810000 :             DO icheb = 2, ncheb
     369      1806900 :                CALL chebyshev_poly(chev_T_dos, chev_Es_dos(igrid), icheb)
     370      1806900 :                dummy1 = dummy1 + kernel_g(icheb)*mu(icheb)*chev_T_dos
     371      4204000 :                DO iwindow = 1, nwindow
     372      4200900 :                   dummy2(iwindow) = dummy2(iwindow) + kernel_g(icheb)*aitchev_T(icheb, iwindow)*chev_T_dos
     373              :                END DO
     374              :             END DO
     375              :             dos(igrid) = 1.0_dp/(interval_a*Nrows* &
     376         3100 :                                  (pi*SQRT(1.0_dp - chev_Es_dos(igrid)*chev_Es_dos(igrid))))*(kernel_g(1)*mu(1) + 2.0_dp*dummy1)
     377         9100 :             DO iwindow = 1, nwindow
     378         9100 :                gdensity(igrid, iwindow) = kernel_g(1)*aitchev_T(1, iwindow) + 2.0_dp*dummy2(iwindow)
     379              :             END DO
     380         3103 :             WRITE (unit_dos, '(1000F16.8)') chev_E(igrid), dos(igrid), gdensity(igrid, :)
     381              :          END DO
     382            3 :          DEALLOCATE (chev_Es_dos, chev_E, dos, gdensity)
     383              :       END IF
     384            6 :       CALL cp_print_key_finished_output(unit_dos, logger, ls_scf_env%chebyshev%print_key_dos, "")
     385              : 
     386              :       ! free the matrices
     387            6 :       CALL dbcsr_release(matrix_tmp1)
     388            6 :       CALL dbcsr_release(matrix_tmp2)
     389            6 :       CALL dbcsr_release(matrix_tmp3)
     390            6 :       CALL dbcsr_release(matrix_F)
     391            6 :       CALL dbcsr_release(matrix_dummy1)
     392              : 
     393           12 :       DO iwindow = 1, nwindow
     394           12 :          CALL dbcsr_release(matrix_dummy2(iwindow))
     395              :       END DO
     396              : 
     397            6 :       DEALLOCATE (ev1, ev2, sev1, sev2, matrix_dummy2)
     398              : 
     399              :       !Need deallocation
     400            6 :       DEALLOCATE (mu, kernel_g, aitchev_T, E_inte, sqrt_vec)
     401              : 
     402            6 :       IF (unit_nr > 0) WRITE (unit_nr, '(T2,A)') "ENDING CHEBYSHEV CALCULATION"
     403              : 
     404            6 :       CALL timestop(handle)
     405              : 
     406           12 :    END SUBROUTINE compute_chebyshev
     407              : 
     408              : END MODULE dm_ls_chebyshev
        

Generated by: LCOV version 2.0-1