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

            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 density matrix optimization using exponential transformations
      10              : !> \par History
      11              : !>       2012.05 created [Florian Schiffmann]
      12              : !> \author Florian Schiffmann
      13              : ! **************************************************************************************************
      14              : 
      15              : MODULE dm_ls_scf_curvy
      16              :    USE bibliography,                    ONLY: Shao2003,&
      17              :                                               cite_reference
      18              :    USE cp_dbcsr_api,                    ONLY: &
      19              :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_filter, dbcsr_multiply, dbcsr_release, &
      20              :         dbcsr_scale, dbcsr_set, dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry
      21              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_dot,&
      22              :                                               dbcsr_frobenius_norm
      23              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      24              :                                               cp_logger_get_default_unit_nr,&
      25              :                                               cp_logger_type
      26              :    USE dm_ls_scf_types,                 ONLY: ls_scf_curvy_type,&
      27              :                                               ls_scf_env_type
      28              :    USE input_constants,                 ONLY: ls_scf_line_search_3point,&
      29              :                                               ls_scf_line_search_3point_2d
      30              :    USE iterate_matrix,                  ONLY: purify_mcweeny
      31              :    USE kinds,                           ONLY: dp
      32              :    USE machine,                         ONLY: m_flush
      33              :    USE mathconstants,                   ONLY: ifac
      34              :    USE mathlib,                         ONLY: invmat
      35              : #include "./base/base_uses.f90"
      36              : 
      37              :    IMPLICIT NONE
      38              : 
      39              :    PRIVATE
      40              : 
      41              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf_curvy'
      42              : 
      43              :    PUBLIC :: dm_ls_curvy_optimization, deallocate_curvy_data
      44              : 
      45              : CONTAINS
      46              : 
      47              : ! **************************************************************************************************
      48              : !> \brief driver routine for Head-Gordon curvy step approach
      49              : !> \param ls_scf_env ...
      50              : !> \param energy ...
      51              : !> \param check_conv ...
      52              : !> \par History
      53              : !>       2012.05 created [Florian Schiffmann]
      54              : !> \author Florian Schiffmann
      55              : ! **************************************************************************************************
      56              : 
      57           90 :    SUBROUTINE dm_ls_curvy_optimization(ls_scf_env, energy, check_conv)
      58              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
      59              :       REAL(KIND=dp)                                      :: energy
      60              :       LOGICAL                                            :: check_conv
      61              : 
      62              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'dm_ls_curvy_optimization'
      63              : 
      64              :       INTEGER                                            :: handle, i, lsstep
      65              : 
      66           90 :       CALL timeset(routineN, handle)
      67              : 
      68           90 :       CALL cite_reference(Shao2003)
      69              : 
      70              : ! Upon first call initialize all matrices needed curing optimization
      71              : ! In addition transform P into orthonormal basis. Will be scaled by 0.5 in closed shell case
      72              : ! Only to be done once as it will be stored and reused afterwards
      73              : ! TRS4 might yield a non-idempotent P therefore McWeeny purification is applied on initial P
      74              : 
      75           90 :       IF (.NOT. ALLOCATED(ls_scf_env%curvy_data%matrix_dp)) THEN
      76           18 :          CALL init_curvy(ls_scf_env%curvy_data, ls_scf_env%matrix_s, ls_scf_env%nspins)
      77           18 :          ls_scf_env%curvy_data%line_search_step = 1
      78              : 
      79           18 :          IF (ls_scf_env%curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
      80            6 :             DO i = 1, ls_scf_env%nspins
      81              :                CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_psave(i, 1), &
      82            6 :                                ls_scf_env%matrix_p(i))
      83              :             END DO
      84              :          END IF
      85           18 :          IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 0.5_dp)
      86              :          CALL transform_matrix_orth(ls_scf_env%matrix_p, ls_scf_env%matrix_s_sqrt, &
      87           18 :                                     ls_scf_env%eps_filter)
      88           18 :          CALL purify_mcweeny(ls_scf_env%matrix_p, ls_scf_env%eps_filter, 3)
      89           38 :          DO i = 1, ls_scf_env%nspins
      90           38 :             CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_p(i), ls_scf_env%matrix_p(i))
      91              :          END DO
      92              :       END IF
      93              : 
      94           90 :       lsstep = ls_scf_env%curvy_data%line_search_step
      95              : 
      96              : ! If new search direction has to be computed transform H into the orthnormal basis
      97              : 
      98           90 :       IF (ls_scf_env%curvy_data%line_search_step == 1) &
      99              :          CALL transform_matrix_orth(ls_scf_env%matrix_ks, ls_scf_env%matrix_s_sqrt_inv, &
     100           28 :                                     ls_scf_env%eps_filter)
     101              : 
     102              : ! Set the energies for the line search and make sure to give the correct energy back to scf_main
     103           90 :       ls_scf_env%curvy_data%energies(lsstep) = energy
     104           90 :       IF (lsstep .NE. 1) energy = ls_scf_env%curvy_data%energies(1)
     105              : 
     106              : ! start the optimization by calling the driver routine or simply combine saved P(2D line search)
     107           90 :       IF (lsstep .LE. 2) THEN
     108           56 :          CALL optimization_step(ls_scf_env%curvy_data, ls_scf_env)
     109           34 :       ELSE IF (lsstep == ls_scf_env%curvy_data%line_search_type) THEN
     110              : ! line_search type has the value appropriate to the number of energy calculations needed
     111           28 :          CALL optimization_step(ls_scf_env%curvy_data, ls_scf_env)
     112              :       ELSE
     113              :          CALL new_p_from_save(ls_scf_env%matrix_p, ls_scf_env%curvy_data%matrix_psave, lsstep, &
     114            6 :                               ls_scf_env%curvy_data%double_step_size)
     115            6 :          ls_scf_env%curvy_data%line_search_step = ls_scf_env%curvy_data%line_search_step + 1
     116            6 :          CALL timestop(handle)
     117            6 :          RETURN
     118              :       END IF
     119           84 :       lsstep = ls_scf_env%curvy_data%line_search_step
     120              : 
     121              : ! transform new density matrix back into nonorthonormal basis (again scaling might apply)
     122              : 
     123              :       CALL transform_matrix_orth(ls_scf_env%matrix_p, ls_scf_env%matrix_s_sqrt_inv, &
     124           84 :                                  ls_scf_env%eps_filter)
     125           84 :       IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(1), 2.0_dp)
     126              : 
     127              : ! P-matrices only need to be stored in case of 2D line search
     128           84 :       IF (lsstep .LE. 3 .AND. ls_scf_env%curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
     129           18 :          DO i = 1, ls_scf_env%nspins
     130              :             CALL dbcsr_copy(ls_scf_env%curvy_data%matrix_psave(i, lsstep), &
     131           18 :                             ls_scf_env%matrix_p(i))
     132              :          END DO
     133              :       END IF
     134           84 :       check_conv = lsstep == 1
     135              : 
     136           84 :       CALL timestop(handle)
     137              : 
     138              :    END SUBROUTINE dm_ls_curvy_optimization
     139              : 
     140              : ! **************************************************************************************************
     141              : !> \brief low level routine for Head-Gordons curvy step approach
     142              : !>        computes gradients, performs a cg and line search,
     143              : !>        and evaluates the BCH series to obtain the new P matrix
     144              : !> \param curvy_data ...
     145              : !> \param ls_scf_env ...
     146              : !> \par History
     147              : !>       2012.05 created [Florian Schiffmann]
     148              : !> \author Florian Schiffmann
     149              : ! **************************************************************************************************
     150              : 
     151           84 :    SUBROUTINE optimization_step(curvy_data, ls_scf_env)
     152              :       TYPE(ls_scf_curvy_type)                            :: curvy_data
     153              :       TYPE(ls_scf_env_type)                              :: ls_scf_env
     154              : 
     155              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'optimization_step'
     156              : 
     157              :       INTEGER                                            :: handle, ispin
     158              :       REAL(KIND=dp)                                      :: filter, step_size(2)
     159              : 
     160              : ! Upon first line search step compute new search direction and apply CG if required
     161              : 
     162           84 :       CALL timeset(routineN, handle)
     163              : 
     164           84 :       IF (curvy_data%line_search_step == 1) THEN
     165          196 :          curvy_data%step_size = MAXVAL(curvy_data%step_size)
     166           84 :          curvy_data%step_size = MIN(MAX(0.10_dp, 0.5_dp*ABS(curvy_data%step_size(1))), 0.5_dp)
     167              : ! Dynamic eps_filter for newton steps
     168              :          filter = MAX(ls_scf_env%eps_filter*curvy_data%min_filter, &
     169           28 :                       ls_scf_env%eps_filter*curvy_data%filter_factor)
     170              :          CALL compute_direction_newton(curvy_data%matrix_p, ls_scf_env%matrix_ks, &
     171              :                                        curvy_data%matrix_dp, filter, curvy_data%fix_shift, curvy_data%shift, &
     172           28 :                                        curvy_data%cg_numer, curvy_data%cg_denom, curvy_data%min_shift)
     173           28 :          curvy_data%filter_factor = curvy_data%scale_filter*curvy_data%filter_factor
     174           84 :          step_size = curvy_data%step_size
     175           84 :          curvy_data%BCH_saved = 0
     176           56 :       ELSE IF (curvy_data%line_search_step == 2) THEN
     177           84 :          step_size = curvy_data%step_size
     178           28 :          IF (curvy_data%energies(1) - curvy_data%energies(2) .GT. 0.0_dp) THEN
     179           84 :             curvy_data%step_size = curvy_data%step_size*2.0_dp
     180           28 :             curvy_data%double_step_size = .TRUE.
     181              :          ELSE
     182            0 :             curvy_data%step_size = curvy_data%step_size*0.5_dp
     183            0 :             curvy_data%double_step_size = .FALSE.
     184              :          END IF
     185           84 :          step_size = curvy_data%step_size
     186           28 :       ELSE IF (curvy_data%line_search_step == ls_scf_line_search_3point_2d) THEN
     187            2 :          CALL line_search_2d(curvy_data%energies, curvy_data%step_size)
     188            6 :          step_size = curvy_data%step_size
     189           26 :       ELSE IF (curvy_data%line_search_step == ls_scf_line_search_3point) THEN
     190           26 :          CALL line_search_3pnt(curvy_data%energies, curvy_data%step_size)
     191           78 :          step_size = curvy_data%step_size
     192              :       END IF
     193              : 
     194              :       CALL update_p_exp(curvy_data%matrix_p, ls_scf_env%matrix_p, curvy_data%matrix_dp, &
     195              :                         curvy_data%matrix_BCH, ls_scf_env%eps_filter, step_size, curvy_data%BCH_saved, &
     196           84 :                         curvy_data%n_bch_hist)
     197              : 
     198              : ! line_search type has the value appropriate to the numeber of energy calculations needed
     199           84 :       curvy_data%line_search_step = MOD(curvy_data%line_search_step, curvy_data%line_search_type) + 1
     200           84 :       IF (curvy_data%line_search_step == 1) THEN
     201           58 :          DO ispin = 1, SIZE(curvy_data%matrix_p)
     202           58 :             CALL dbcsr_copy(curvy_data%matrix_p(ispin), ls_scf_env%matrix_p(ispin))
     203              :          END DO
     204              :       END IF
     205           84 :       CALL timestop(handle)
     206              : 
     207           84 :    END SUBROUTINE optimization_step
     208              : 
     209              : ! **************************************************************************************************
     210              : !> \brief Perform a 6pnt-2D line search for spin polarized calculations.
     211              : !>        Fit a 2D parabolic function to 6 points
     212              : !> \param energies ...
     213              : !> \param step_size ...
     214              : !> \par History
     215              : !>       2012.05 created [Florian Schiffmann]
     216              : !> \author Florian Schiffmann
     217              : ! **************************************************************************************************
     218              : 
     219            2 :    SUBROUTINE line_search_2d(energies, step_size)
     220              :       REAL(KIND=dp)                                      :: energies(6), step_size(2)
     221              : 
     222              :       INTEGER                                            :: info, unit_nr
     223              :       REAL(KIND=dp)                                      :: e_pred, param(6), s1, s1sq, s2, s2sq, &
     224              :                                                             sys_lin_eq(6, 6), tmp_e, v1, v2
     225              :       TYPE(cp_logger_type), POINTER                      :: logger
     226              : 
     227            2 :       logger => cp_get_default_logger()
     228            2 :       IF (energies(1) - energies(2) .LT. 0._dp) THEN
     229            0 :          tmp_e = energies(2); energies(2) = energies(3); energies(3) = tmp_e
     230            0 :          step_size = step_size*2.0_dp
     231              :       END IF
     232            2 :       IF (logger%para_env%is_source()) THEN
     233            1 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     234              :       ELSE
     235              :          unit_nr = -1
     236              :       END IF
     237            2 :       s1 = 0.5_dp*step_size(1); s2 = step_size(1); s1sq = s1**2; s2sq = s2**2
     238           14 :       sys_lin_eq = 0.0_dp; sys_lin_eq(:, 6) = 1.0_dp
     239            2 :       sys_lin_eq(2, 1) = s1sq; sys_lin_eq(2, 2) = s1sq; sys_lin_eq(2, 3) = s1sq; sys_lin_eq(2, 4) = s1; sys_lin_eq(2, 5) = s1
     240            2 :       sys_lin_eq(3, 1) = s2sq; sys_lin_eq(3, 2) = s2sq; sys_lin_eq(3, 3) = s2sq; sys_lin_eq(3, 4) = s2; sys_lin_eq(3, 5) = s2
     241            2 :       sys_lin_eq(4, 3) = s1sq; sys_lin_eq(4, 5) = s1
     242            2 :       sys_lin_eq(5, 1) = s1sq; sys_lin_eq(5, 4) = s1
     243            2 :       sys_lin_eq(6, 3) = s2sq; sys_lin_eq(6, 5) = s2
     244              : 
     245            2 :       CALL invmat(sys_lin_eq, info)
     246           86 :       param = MATMUL(sys_lin_eq, energies)
     247            2 :       v1 = (param(2)*param(4))/(2.0_dp*param(1)) - param(5)
     248            2 :       v2 = -(param(2)**2)/(2.0_dp*param(1)) + 2.0_dp*param(3)
     249            2 :       step_size(2) = v1/v2
     250            2 :       step_size(1) = (-param(2)*step_size(2) - param(4))/(2.0_dp*param(1))
     251            2 :       IF (step_size(1) .LT. 0.0_dp) step_size(1) = 1.0_dp
     252            2 :       IF (step_size(2) .LT. 0.0_dp) step_size(2) = 1.0_dp
     253              : !    step_size(1)=MIN(step_size(1),2.0_dp)
     254              : !    step_size(2)=MIN(step_size(2),2.0_dp)
     255              :       e_pred = param(1)*step_size(1)**2 + param(2)*step_size(1)*step_size(2) + &
     256            2 :                param(3)*step_size(2)**2 + param(4)*step_size(1) + param(5)*step_size(2) + param(6)
     257            2 :       IF (unit_nr .GT. 0) WRITE (unit_nr, "(t3,a,F10.5,F10.5,A,F20.9)") &
     258            1 :          " Line Search: Step Size", step_size, " Predicted energy", e_pred
     259              :       e_pred = param(1)*s1**2 + param(2)*s2*s1*0.0_dp + &
     260              :                param(3)*s1**2*0.0_dp + param(4)*s1 + param(5)*s1*0.0_dp + param(6)
     261              : 
     262            2 :    END SUBROUTINE line_search_2d
     263              : 
     264              : ! **************************************************************************************************
     265              : !> \brief Perform a 3pnt line search
     266              : !> \param energies ...
     267              : !> \param step_size ...
     268              : !> \par History
     269              : !>       2012.05 created [Florian Schiffmann]
     270              : !> \author Florian Schiffmann
     271              : ! **************************************************************************************************
     272              : 
     273           26 :    SUBROUTINE line_search_3pnt(energies, step_size)
     274              :       REAL(KIND=dp)                                      :: energies(3), step_size(2)
     275              : 
     276              :       INTEGER                                            :: unit_nr
     277              :       REAL(KIND=dp)                                      :: a, b, c, e_pred, min_val, step1, tmp, &
     278              :                                                             tmp_e
     279              :       TYPE(cp_logger_type), POINTER                      :: logger
     280              : 
     281           26 :       logger => cp_get_default_logger()
     282           26 :       IF (energies(1) - energies(2) .LT. 0._dp) THEN
     283            0 :          tmp_e = energies(2); energies(2) = energies(3); energies(3) = tmp_e
     284            0 :          step_size = step_size*2.0_dp
     285              :       END IF
     286           26 :       IF (logger%para_env%is_source()) THEN
     287           13 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     288              :       ELSE
     289           13 :          unit_nr = -1
     290              :       END IF
     291           26 :       step1 = 0.5_dp*step_size(1)
     292           26 :       c = energies(1)
     293           26 :       a = (energies(3) + c - 2.0_dp*energies(2))/(2.0_dp*step1**2)
     294           26 :       b = (energies(2) - c - a*step1**2)/step1
     295           26 :       IF (a .LT. 1.0E-12_dp) a = -1.0E-12_dp
     296           26 :       min_val = -b/(2.0_dp*a)
     297           26 :       e_pred = a*min_val**2 + b*min_val + c
     298           26 :       tmp = step_size(1)
     299           26 :       IF (e_pred .LT. energies(1) .AND. e_pred .LT. energies(2)) THEN
     300              :          step_size = MAX(-1.0_dp, &
     301           54 :                          MIN(min_val, 10_dp*step_size))
     302              :       ELSE
     303           24 :          step_size = 1.0_dp
     304              :       END IF
     305           26 :       e_pred = a*(step_size(1))**2 + b*(step_size(1)) + c
     306           26 :       IF (unit_nr .GT. 0) THEN
     307           13 :          WRITE (unit_nr, "(t3,a,f16.8,a,F20.9)") "Line Search: Step Size", step_size(1), " Predicted energy", e_pred
     308           13 :          CALL m_flush(unit_nr)
     309              :       END IF
     310           26 :    END SUBROUTINE line_search_3pnt
     311              : 
     312              : ! **************************************************************************************************
     313              : !> \brief Get a new search direction. Iterate to obtain a Newton like step
     314              : !>        Refine with a CG update of the search direction
     315              : !> \param matrix_p ...
     316              : !> \param matrix_ks ...
     317              : !> \param matrix_dp ...
     318              : !> \param eps_filter ...
     319              : !> \param fix_shift ...
     320              : !> \param curvy_shift ...
     321              : !> \param cg_numer ...
     322              : !> \param cg_denom ...
     323              : !> \param min_shift ...
     324              : !> \par History
     325              : !>       2012.05 created [Florian Schiffmann]
     326              : !> \author Florian Schiffmann
     327              : ! **************************************************************************************************
     328              : 
     329           28 :    SUBROUTINE compute_direction_newton(matrix_p, matrix_ks, matrix_dp, eps_filter, fix_shift, &
     330              :                                        curvy_shift, cg_numer, cg_denom, min_shift)
     331              :       TYPE(dbcsr_type), DIMENSION(:)                     :: matrix_p, matrix_ks, matrix_dp
     332              :       REAL(KIND=dp)                                      :: eps_filter
     333              :       LOGICAL                                            :: fix_shift(2)
     334              :       REAL(KIND=dp)                                      :: curvy_shift(2), cg_numer(2), &
     335              :                                                             cg_denom(2), min_shift
     336              : 
     337              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_direction_newton'
     338              : 
     339              :       INTEGER                                            :: handle, i, ispin, ncyc, nspin, unit_nr
     340              :       LOGICAL                                            :: at_limit
     341              :       REAL(KIND=dp)                                      :: beta, conv_val, maxel, old_conv, shift
     342              :       TYPE(cp_logger_type), POINTER                      :: logger
     343              :       TYPE(dbcsr_type)                                   :: matrix_Ax, matrix_b, matrix_cg, &
     344              :                                                             matrix_dp_old, matrix_PKs, matrix_res, &
     345              :                                                             matrix_tmp, matrix_tmp1
     346              : 
     347           56 :       logger => cp_get_default_logger()
     348              : 
     349           28 :       IF (logger%para_env%is_source()) THEN
     350           14 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     351              :       ELSE
     352           14 :          unit_nr = -1
     353              :       END IF
     354           28 :       CALL timeset(routineN, handle)
     355           28 :       nspin = SIZE(matrix_p)
     356              : 
     357           28 :       CALL dbcsr_create(matrix_PKs, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
     358           28 :       CALL dbcsr_create(matrix_Ax, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
     359           28 :       CALL dbcsr_create(matrix_tmp, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
     360           28 :       CALL dbcsr_create(matrix_tmp1, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
     361           28 :       CALL dbcsr_create(matrix_res, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
     362           28 :       CALL dbcsr_create(matrix_cg, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
     363           28 :       CALL dbcsr_create(matrix_b, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
     364           28 :       CALL dbcsr_create(matrix_dp_old, template=matrix_dp(1), matrix_type=dbcsr_type_no_symmetry)
     365              : 
     366           58 :       DO ispin = 1, nspin
     367           30 :          CALL dbcsr_copy(matrix_dp_old, matrix_dp(ispin))
     368              : 
     369              : ! Precompute some matrices to save work during iterations
     370              :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p(ispin), matrix_ks(ispin), &
     371           30 :                              0.0_dp, matrix_PKs, filter_eps=eps_filter)
     372           30 :          CALL dbcsr_transposed(matrix_b, matrix_PKs)
     373           30 :          CALL dbcsr_copy(matrix_cg, matrix_b)
     374              : 
     375              : ! Starting CG with guess 0-matrix gives -2*gradient=[Ks*P-(Ks*P)T] for cg_matrix in second step
     376           30 :          CALL dbcsr_add(matrix_cg, matrix_PKs, 2.0_dp, -2.0_dp)
     377              : 
     378              : ! Residual matrix in first step=cg matrix. Keep Pks for later use in CG!
     379           30 :          CALL dbcsr_copy(matrix_res, matrix_cg)
     380              : 
     381              : ! Precompute -FP-[FP]T which will be used throughout the CG iterations
     382           30 :          CALL dbcsr_add(matrix_b, matrix_PKs, -1.0_dp, -1.0_dp)
     383              : 
     384              : ! Setup some values to check convergence and safety checks for eigenvalue shifting
     385           30 :          old_conv = dbcsr_frobenius_norm(matrix_res)
     386           30 :          shift = MIN(10.0_dp, MAX(min_shift, 0.05_dp*old_conv))
     387           30 :          conv_val = MAX(0.010_dp*old_conv, 100.0_dp*eps_filter)
     388           30 :          old_conv = 100.0_dp
     389           30 :          IF (fix_shift(ispin)) THEN
     390            0 :             shift = MAX(min_shift, MIN(10.0_dp, MAX(shift, curvy_shift(ispin) - 0.5_dp*curvy_shift(ispin))))
     391            0 :             curvy_shift(ispin) = shift
     392              :          END IF
     393              : 
     394              : ! Begin the real optimization loop
     395           30 :          CALL dbcsr_set(matrix_dp(ispin), 0.0_dp)
     396           30 :          ncyc = 10
     397          104 :          DO i = 1, ncyc
     398              : 
     399              : ! One step to compute: -FPD-DPF-DFP-PFD (not obvious but symmetry allows for some tricks)
     400          104 :             CALL commutator_symm(matrix_b, matrix_cg, matrix_Ax, eps_filter, 1.0_dp)
     401              : 
     402              : ! Compute the missing bits 2*(FDP+PDF) (again use symmetry to compute as a commutator)
     403              :             CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_cg, matrix_p(ispin), &
     404          104 :                                 0.0_dp, matrix_tmp, filter_eps=eps_filter)
     405          104 :             CALL commutator_symm(matrix_ks(ispin), matrix_tmp, matrix_tmp1, eps_filter, 2.0_dp)
     406          104 :             CALL dbcsr_add(matrix_Ax, matrix_tmp1, 1.0_dp, 1.0_dp)
     407              : 
     408              : ! Apply the shift and hope it's enough to stabilize the CG iterations
     409          104 :             CALL dbcsr_add(matrix_Ax, matrix_cg, 1.0_dp, shift)
     410              : 
     411              :             CALL compute_cg_matrices(matrix_Ax, matrix_res, matrix_cg, matrix_dp(ispin), &
     412          104 :                                      matrix_tmp, eps_filter, at_limit)
     413          104 :             CALL dbcsr_filter(matrix_cg, eps_filter)
     414              : 
     415              : ! check for convergence of the newton step
     416          104 :             maxel = dbcsr_frobenius_norm(matrix_res)
     417          104 :             IF (unit_nr .GT. 0) THEN
     418           52 :                WRITE (unit_nr, "(T3,A,F12.6)") "Convergence of Newton iteration ", maxel
     419           52 :                CALL m_flush(unit_nr)
     420              :             END IF
     421          104 :             at_limit = at_limit .OR. (old_conv/maxel .LT. 1.01_dp)
     422          104 :             old_conv = maxel
     423          104 :             IF (i == ncyc .AND. maxel/conv_val .GT. 5.0_dp) THEN
     424            0 :                fix_shift(ispin) = .TRUE.
     425            0 :                curvy_shift(ispin) = 4.0_dp*shift
     426              :             END IF
     427          104 :             IF (maxel .LT. conv_val .OR. at_limit) EXIT
     428              :          END DO
     429              : 
     430              : ! Refine the Newton like search direction with a preconditioned cg update
     431           30 :          CALL dbcsr_transposed(matrix_b, matrix_PKs)
     432              :          !compute b= -2*KsP+2*PKs=-(2*gradient)
     433           30 :          CALL dbcsr_copy(matrix_cg, matrix_b)
     434           30 :          CALL dbcsr_add(matrix_cg, matrix_PKs, 1.0_dp, -1.0_dp)
     435           30 :          cg_denom(ispin) = cg_numer(ispin)
     436           30 :          CALL dbcsr_dot(matrix_cg, matrix_dp(ispin), cg_numer(ispin))
     437           30 :          beta = cg_numer(ispin)/MAX(cg_denom(ispin), 1.0E-6_dp)
     438           30 :          IF (beta .LT. 1.0_dp) THEN
     439           28 :             beta = MAX(0.0_dp, beta)
     440           28 :             CALL dbcsr_add(matrix_dp(ispin), matrix_dp_old, 1.0_dp, beta)
     441              :          END IF
     442           58 :          IF (unit_nr .GT. 0) WRITE (unit_nr, "(A)") " "
     443              :       END DO
     444              : 
     445           28 :       CALL dbcsr_release(matrix_PKs)
     446           28 :       CALL dbcsr_release(matrix_dp_old)
     447           28 :       CALL dbcsr_release(matrix_b)
     448           28 :       CALL dbcsr_release(matrix_Ax)
     449           28 :       CALL dbcsr_release(matrix_tmp)
     450           28 :       CALL dbcsr_release(matrix_tmp1)
     451           28 :       CALL dbcsr_release(matrix_b)
     452           28 :       CALL dbcsr_release(matrix_res)
     453           28 :       CALL dbcsr_release(matrix_cg)
     454              : 
     455           28 :       IF (unit_nr .GT. 0) CALL m_flush(unit_nr)
     456           28 :       CALL timestop(handle)
     457           28 :    END SUBROUTINE compute_direction_newton
     458              : 
     459              : ! **************************************************************************************************
     460              : !> \brief compute the optimal step size of the current cycle and update the
     461              : !>        matrices needed to solve the system of linear equations
     462              : !> \param Ax ...
     463              : !> \param res ...
     464              : !> \param cg ...
     465              : !> \param deltp ...
     466              : !> \param tmp ...
     467              : !> \param eps_filter ...
     468              : !> \param at_limit ...
     469              : !> \par History
     470              : !>       2012.05 created [Florian Schiffmann]
     471              : !> \author Florian Schiffmann
     472              : ! **************************************************************************************************
     473              : 
     474          104 :    SUBROUTINE compute_cg_matrices(Ax, res, cg, deltp, tmp, eps_filter, at_limit)
     475              :       TYPE(dbcsr_type)                                   :: Ax, res, cg, deltp, tmp
     476              :       REAL(KIND=dp)                                      :: eps_filter
     477              :       LOGICAL                                            :: at_limit
     478              : 
     479              :       INTEGER                                            :: i, info
     480              :       REAL(KIND=dp)                                      :: alpha, beta, devi(3), fac, fac1, &
     481              :                                                             lin_eq(3, 3), new_norm, norm_cA, &
     482              :                                                             norm_rr, vec(3)
     483              : 
     484          104 :       at_limit = .FALSE.
     485          104 :       CALL dbcsr_dot(res, res, norm_rr)
     486          104 :       CALL dbcsr_dot(cg, Ax, norm_cA)
     487          104 :       lin_eq = 0.0_dp
     488          104 :       fac = norm_rr/norm_cA
     489          104 :       fac1 = fac
     490              : ! Use a 3point line search and a fit to a quadratic function to determine optimal step size
     491          416 :       DO i = 1, 3
     492          312 :          CALL dbcsr_copy(tmp, res)
     493          312 :          CALL dbcsr_add(tmp, Ax, 1.0_dp, -fac)
     494          312 :          devi(i) = dbcsr_frobenius_norm(tmp)
     495         1248 :          lin_eq(i, :) = (/fac**2, fac, 1.0_dp/)
     496          416 :          fac = fac1 + fac1*((-1)**i)*0.5_dp
     497              :       END DO
     498          104 :       CALL invmat(lin_eq, info)
     499         1352 :       vec = MATMUL(lin_eq, devi)
     500          104 :       alpha = -vec(2)/(2.0_dp*vec(1))
     501          104 :       fac = SQRT(norm_rr/(norm_cA*alpha))
     502              : !scale the previous matrices to match the step size
     503          104 :       CALL dbcsr_scale(Ax, fac)
     504          104 :       CALL dbcsr_scale(cg, fac)
     505          104 :       norm_cA = norm_cA*fac**2
     506              : 
     507              : ! USe CG to get the new matrices
     508          104 :       alpha = norm_rr/norm_cA
     509          104 :       CALL dbcsr_add(res, Ax, 1.0_dp, -alpha)
     510          104 :       CALL dbcsr_dot(res, res, new_norm)
     511          104 :       IF (norm_rr .LT. eps_filter*0.001_dp .OR. new_norm .LT. eps_filter*0.001_dp) THEN
     512              :          beta = 0.0_dp
     513           22 :          at_limit = .TRUE.
     514              :       ELSE
     515              :          beta = new_norm/norm_rr
     516           82 :          CALL dbcsr_add(deltp, cg, 1.0_dp, alpha)
     517              :       END IF
     518          104 :       beta = new_norm/norm_rr
     519          104 :       CALL dbcsr_add(cg, res, beta, 1.0_dp)
     520              : 
     521          208 :    END SUBROUTINE compute_cg_matrices
     522              : 
     523              : ! **************************************************************************************************
     524              : !> \brief Only for 2D line search. Use saved P-components to construct new
     525              : !>        test density matrix. Takes care as well, whether step_size
     526              : !>        increased or decreased during 2nd step and combines matrices accordingly
     527              : !> \param matrix_p ...
     528              : !> \param matrix_psave ...
     529              : !> \param lsstep ...
     530              : !> \param DOUBLE ...
     531              : !> \par History
     532              : !>       2012.05 created [Florian Schiffmann]
     533              : !> \author Florian Schiffmann
     534              : ! **************************************************************************************************
     535              : 
     536            6 :    SUBROUTINE new_p_from_save(matrix_p, matrix_psave, lsstep, DOUBLE)
     537              :       TYPE(dbcsr_type), DIMENSION(:)                     :: matrix_p
     538              :       TYPE(dbcsr_type), DIMENSION(:, :)                  :: matrix_psave
     539              :       INTEGER                                            :: lsstep
     540              :       LOGICAL                                            :: DOUBLE
     541              : 
     542            8 :       SELECT CASE (lsstep)
     543              :       CASE (3)
     544            2 :          CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 1))
     545            2 :          IF (DOUBLE) THEN
     546            2 :             CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 2))
     547              :          ELSE
     548            0 :             CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 3))
     549              :          END IF
     550              :       CASE (4)
     551            2 :          IF (DOUBLE) THEN
     552            2 :             CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 2))
     553              :          ELSE
     554            0 :             CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 3))
     555              :          END IF
     556            2 :          CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 1))
     557              :       CASE (5)
     558            2 :          CALL dbcsr_copy(matrix_p(1), matrix_psave(1, 1))
     559            8 :          IF (DOUBLE) THEN
     560            2 :             CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 3))
     561              :          ELSE
     562            0 :             CALL dbcsr_copy(matrix_p(2), matrix_psave(2, 2))
     563              :          END IF
     564              :       END SELECT
     565              : 
     566            6 :    END SUBROUTINE new_p_from_save
     567              : 
     568              : ! **************************************************************************************************
     569              : !> \brief computes a commutator exploiting symmetry RES=k*[A,B]=k*[AB-(AB)T]
     570              : !> \param a ...
     571              : !> \param b ...
     572              : !> \param res ...
     573              : !> \param eps_filter   filtering threshold for sparse matrices
     574              : !> \param prefac      prefactor k in above equation
     575              : !> \par History
     576              : !>       2012.05 created [Florian Schiffmann]
     577              : !> \author Florian Schiffmann
     578              : ! **************************************************************************************************
     579              : 
     580          208 :    SUBROUTINE commutator_symm(a, b, res, eps_filter, prefac)
     581              :       TYPE(dbcsr_type)                                   :: a, b, res
     582              :       REAL(KIND=dp)                                      :: eps_filter, prefac
     583              : 
     584              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'commutator_symm'
     585              : 
     586              :       INTEGER                                            :: handle
     587              :       TYPE(dbcsr_type)                                   :: work
     588              : 
     589          208 :       CALL timeset(routineN, handle)
     590              : 
     591          208 :       CALL dbcsr_create(work, template=a, matrix_type=dbcsr_type_no_symmetry)
     592              : 
     593          208 :       CALL dbcsr_multiply("N", "N", prefac, a, b, 0.0_dp, res, filter_eps=eps_filter)
     594          208 :       CALL dbcsr_transposed(work, res)
     595          208 :       CALL dbcsr_add(res, work, 1.0_dp, -1.0_dp)
     596              : 
     597          208 :       CALL dbcsr_release(work)
     598              : 
     599          208 :       CALL timestop(handle)
     600          208 :    END SUBROUTINE commutator_symm
     601              : 
     602              : ! **************************************************************************************************
     603              : !> \brief Use the BCH update to get the new idempotent P
     604              : !>        Numerics don't allow for perfect idempotency, therefore a mc weeny
     605              : !>        step is used to make sure we stay close to the idempotent surface
     606              : !> \param matrix_p_in ...
     607              : !> \param matrix_p_out ...
     608              : !> \param matrix_dp ...
     609              : !> \param matrix_BCH ...
     610              : !> \param threshold ...
     611              : !> \param step_size ...
     612              : !> \param BCH_saved ...
     613              : !> \param n_bch_hist ...
     614              : !> \par History
     615              : !>       2012.05 created [Florian Schiffmann]
     616              : !> \author Florian Schiffmann
     617              : ! **************************************************************************************************
     618              : 
     619           84 :    SUBROUTINE update_p_exp(matrix_p_in, matrix_p_out, matrix_dp, matrix_BCH, threshold, step_size, &
     620              :                            BCH_saved, n_bch_hist)
     621              :       TYPE(dbcsr_type), DIMENSION(:)                     :: matrix_p_in, matrix_p_out, matrix_dp
     622              :       TYPE(dbcsr_type), DIMENSION(:, :)                  :: matrix_BCH
     623              :       REAL(KIND=dp)                                      :: threshold, step_size(2)
     624              :       INTEGER                                            :: BCH_saved(2), n_bch_hist
     625              : 
     626              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'update_p_exp'
     627              : 
     628              :       INTEGER                                            :: handle, i, ispin, nsave, nspin, unit_nr
     629              :       LOGICAL                                            :: save_BCH
     630              :       REAL(KIND=dp)                                      :: frob_norm, step_fac
     631              :       TYPE(cp_logger_type), POINTER                      :: logger
     632              :       TYPE(dbcsr_type)                                   :: matrix, matrix_tmp
     633              : 
     634           84 :       CALL timeset(routineN, handle)
     635              : 
     636           84 :       logger => cp_get_default_logger()
     637           84 :       IF (logger%para_env%is_source()) THEN
     638           42 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     639              :       ELSE
     640           42 :          unit_nr = -1
     641              :       END IF
     642              : 
     643           84 :       CALL dbcsr_create(matrix, template=matrix_p_in(1), matrix_type=dbcsr_type_no_symmetry)
     644           84 :       CALL dbcsr_create(matrix_tmp, template=matrix_p_in(1), matrix_type=dbcsr_type_no_symmetry)
     645           84 :       nspin = SIZE(matrix_p_in)
     646              : 
     647          174 :       DO ispin = 1, nspin
     648           90 :          step_fac = 1.0_dp
     649           90 :          frob_norm = 1.0_dp
     650           90 :          nsave = 0
     651              : 
     652           90 :          CALL dbcsr_copy(matrix_tmp, matrix_p_in(ispin))
     653           90 :          CALL dbcsr_copy(matrix_p_out(ispin), matrix_p_in(ispin))
     654              : ! If a BCH history is used make good use of it and do a few steps as a copy and scale update of P
     655              : ! else BCH_saved will be 0 and loop is skipped
     656          130 :          DO i = 1, BCH_saved(ispin)
     657           86 :             step_fac = step_fac*step_size(ispin)
     658           86 :             CALL dbcsr_copy(matrix_tmp, matrix_p_out(ispin))
     659           86 :             CALL dbcsr_add(matrix_p_out(ispin), matrix_BCH(ispin, i), 1.0_dp, ifac(i)*step_fac)
     660           86 :             CALL dbcsr_add(matrix_tmp, matrix_p_out(ispin), 1.0_dp, -1.0_dp)
     661           86 :             frob_norm = dbcsr_frobenius_norm(matrix_tmp)
     662           86 :             IF (unit_nr .GT. 0) WRITE (unit_nr, "(t3,a,i3,a,f16.8)") "BCH: step", i, " Norm of P_old-Pnew:", frob_norm
     663          130 :             IF (frob_norm .LT. threshold) EXIT
     664              :          END DO
     665           90 :          IF (frob_norm .LT. threshold) CYCLE
     666              : 
     667              : ! If the copy and scale isn't enough compute a few more BCH steps. 20 seems high but except of the first step it will never be close
     668           44 :          save_BCH = BCH_saved(ispin) == 0 .AND. n_bch_hist .GT. 0
     669           86 :          DO i = BCH_saved(ispin) + 1, 20
     670           86 :             step_fac = step_fac*step_size(ispin)
     671              :             !allow for a bit of matrix magic here by exploiting matrix and matrix_tmp
     672              :             !matrix_tmp is alway the previous order of the BCH series
     673              :             CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_dp(ispin), &
     674           86 :                                 0.0_dp, matrix, filter_eps=threshold)
     675              : 
     676              :             !(anti)symmetry allows to sum the transposed instead of the full commutator, matrix becomes the latest result
     677              : 
     678           86 :             CALL dbcsr_transposed(matrix_tmp, matrix)
     679           86 :             CALL dbcsr_add(matrix, matrix_tmp, 1.0_dp, 1.0_dp)
     680              : 
     681              :             !Finally, add the new BCH order to P, but store the previous one for a convergence check
     682           86 :             CALL dbcsr_copy(matrix_tmp, matrix_p_out(ispin))
     683           86 :             CALL dbcsr_add(matrix_p_out(ispin), matrix, 1.0_dp, ifac(i)*step_fac)
     684           86 :             IF (save_BCH .AND. i .LE. n_bch_hist) THEN
     685           78 :                CALL dbcsr_copy(matrix_BCH(ispin, i), matrix)
     686           78 :                nsave = i
     687              :             END IF
     688              : 
     689           86 :             CALL dbcsr_add(matrix_tmp, matrix_p_out(ispin), 1.0_dp, -1.0_dp)
     690              : 
     691              :             !Stop the BCH-series if two successive P's differ by less the threshold
     692           86 :             frob_norm = dbcsr_frobenius_norm(matrix_tmp)
     693           86 :             IF (unit_nr .GT. 0) WRITE (unit_nr, "(t3,a,i3,a,f16.8)") "BCH: step", i, " Norm of P_old-Pnew:", frob_norm
     694           86 :             IF (frob_norm .LT. threshold) EXIT
     695              : 
     696              :             !Copy the latest BCH-matrix on matrix tmp, so we can cycle with all matrices in place
     697           42 :             CALL dbcsr_copy(matrix_tmp, matrix)
     698           86 :             CALL dbcsr_filter(matrix_tmp, threshold)
     699              :          END DO
     700           44 :          BCH_saved(ispin) = nsave
     701          128 :          IF (unit_nr .GT. 0) WRITE (unit_nr, "(A)") " "
     702              :       END DO
     703              : 
     704           84 :       CALL purify_mcweeny(matrix_p_out, threshold, 1)
     705           84 :       IF (unit_nr .GT. 0) CALL m_flush(unit_nr)
     706           84 :       CALL dbcsr_release(matrix_tmp)
     707           84 :       CALL dbcsr_release(matrix)
     708           84 :       CALL timestop(handle)
     709           84 :    END SUBROUTINE update_p_exp
     710              : 
     711              : ! **************************************************************************************************
     712              : !> \brief performs a transformation of a matrix back to/into orthonormal basis
     713              : !>        in case of P a scaling of 0.5 has to be applied for closed shell case
     714              : !> \param matrix       matrix to be transformed
     715              : !> \param matrix_trafo transformation matrix
     716              : !> \param eps_filter   filtering threshold for sparse matrices
     717              : !> \par History
     718              : !>       2012.05 created [Florian Schiffmann]
     719              : !> \author Florian Schiffmann
     720              : ! **************************************************************************************************
     721              : 
     722          130 :    SUBROUTINE transform_matrix_orth(matrix, matrix_trafo, eps_filter)
     723              :       TYPE(dbcsr_type), DIMENSION(:)                     :: matrix
     724              :       TYPE(dbcsr_type)                                   :: matrix_trafo
     725              :       REAL(KIND=dp)                                      :: eps_filter
     726              : 
     727              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'transform_matrix_orth'
     728              : 
     729              :       INTEGER                                            :: handle, ispin
     730              :       TYPE(dbcsr_type)                                   :: matrix_tmp, matrix_work
     731              : 
     732          130 :       CALL timeset(routineN, handle)
     733              : 
     734          130 :       CALL dbcsr_create(matrix_work, template=matrix(1), matrix_type=dbcsr_type_no_symmetry)
     735          130 :       CALL dbcsr_create(matrix_tmp, template=matrix(1), matrix_type=dbcsr_type_no_symmetry)
     736              : 
     737          270 :       DO ispin = 1, SIZE(matrix)
     738              :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix(ispin), matrix_trafo, &
     739          140 :                              0.0_dp, matrix_work, filter_eps=eps_filter)
     740              :          CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_trafo, matrix_work, &
     741          140 :                              0.0_dp, matrix_tmp, filter_eps=eps_filter)
     742              :          ! symmetrize results (this is again needed to make sure everything is stable)
     743          140 :          CALL dbcsr_transposed(matrix_work, matrix_tmp)
     744          140 :          CALL dbcsr_add(matrix_tmp, matrix_work, 0.5_dp, 0.5_dp)
     745          270 :          CALL dbcsr_copy(matrix(ispin), matrix_tmp)
     746              :       END DO
     747              : 
     748          130 :       CALL dbcsr_release(matrix_tmp)
     749          130 :       CALL dbcsr_release(matrix_work)
     750          130 :       CALL timestop(handle)
     751              : 
     752          130 :    END SUBROUTINE
     753              : 
     754              : ! **************************************************************************************************
     755              : !> \brief ...
     756              : !> \param curvy_data ...
     757              : ! **************************************************************************************************
     758          882 :    SUBROUTINE deallocate_curvy_data(curvy_data)
     759              :       TYPE(ls_scf_curvy_type)                            :: curvy_data
     760              : 
     761              :       INTEGER                                            :: i, j
     762              : 
     763          882 :       CALL release_dbcsr_array(curvy_data%matrix_dp)
     764          882 :       CALL release_dbcsr_array(curvy_data%matrix_p)
     765              : 
     766          882 :       IF (ALLOCATED(curvy_data%matrix_psave)) THEN
     767            6 :          DO i = 1, SIZE(curvy_data%matrix_psave, 1)
     768           18 :             DO j = 1, 3
     769           16 :                CALL dbcsr_release(curvy_data%matrix_psave(i, j))
     770              :             END DO
     771              :          END DO
     772            2 :          DEALLOCATE (curvy_data%matrix_psave)
     773              :       END IF
     774          882 :       IF (ALLOCATED(curvy_data%matrix_BCH)) THEN
     775           38 :          DO i = 1, SIZE(curvy_data%matrix_BCH, 1)
     776          178 :             DO j = 1, 7
     777          160 :                CALL dbcsr_release(curvy_data%matrix_BCH(i, j))
     778              :             END DO
     779              :          END DO
     780           18 :          DEALLOCATE (curvy_data%matrix_BCH)
     781              :       END IF
     782          882 :    END SUBROUTINE deallocate_curvy_data
     783              : 
     784              : ! **************************************************************************************************
     785              : !> \brief ...
     786              : !> \param matrix ...
     787              : ! **************************************************************************************************
     788         1764 :    SUBROUTINE release_dbcsr_array(matrix)
     789              :       TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:)        :: matrix
     790              : 
     791              :       INTEGER                                            :: i
     792              : 
     793         1764 :       IF (ALLOCATED(matrix)) THEN
     794           76 :          DO i = 1, SIZE(matrix)
     795           76 :             CALL dbcsr_release(matrix(i))
     796              :          END DO
     797           36 :          DEALLOCATE (matrix)
     798              :       END IF
     799         1764 :    END SUBROUTINE release_dbcsr_array
     800              : 
     801              : ! **************************************************************************************************
     802              : !> \brief ...
     803              : !> \param curvy_data ...
     804              : !> \param matrix_s ...
     805              : !> \param nspins ...
     806              : ! **************************************************************************************************
     807           18 :    SUBROUTINE init_curvy(curvy_data, matrix_s, nspins)
     808              :       TYPE(ls_scf_curvy_type)                            :: curvy_data
     809              :       TYPE(dbcsr_type)                                   :: matrix_s
     810              :       INTEGER                                            :: nspins
     811              : 
     812              :       INTEGER                                            :: ispin, j
     813              : 
     814           74 :       ALLOCATE (curvy_data%matrix_dp(nspins))
     815           74 :       ALLOCATE (curvy_data%matrix_p(nspins))
     816           38 :       DO ispin = 1, nspins
     817              :          CALL dbcsr_create(curvy_data%matrix_dp(ispin), template=matrix_s, &
     818           20 :                            matrix_type=dbcsr_type_no_symmetry)
     819           20 :          CALL dbcsr_set(curvy_data%matrix_dp(ispin), 0.0_dp)
     820              :          CALL dbcsr_create(curvy_data%matrix_p(ispin), template=matrix_s, &
     821           20 :                            matrix_type=dbcsr_type_no_symmetry)
     822           60 :          curvy_data%fix_shift = .FALSE.
     823           20 :          curvy_data%double_step_size = .TRUE.
     824           60 :          curvy_data%shift = 1.0_dp
     825           60 :          curvy_data%BCH_saved = 0
     826           60 :          curvy_data%step_size = 0.60_dp
     827           60 :          curvy_data%cg_numer = 0.00_dp
     828           78 :          curvy_data%cg_denom = 0.00_dp
     829              :       END DO
     830           18 :       IF (curvy_data%line_search_type == ls_scf_line_search_3point_2d) THEN
     831           24 :          ALLOCATE (curvy_data%matrix_psave(nspins, 3))
     832            6 :          DO ispin = 1, nspins
     833           18 :             DO j = 1, 3
     834              :                CALL dbcsr_create(curvy_data%matrix_psave(ispin, j), template=matrix_s, &
     835           16 :                                  matrix_type=dbcsr_type_no_symmetry)
     836              :             END DO
     837              :          END DO
     838              :       END IF
     839           18 :       IF (curvy_data%n_bch_hist .GT. 0) THEN
     840          338 :          ALLOCATE (curvy_data%matrix_BCH(nspins, curvy_data%n_bch_hist))
     841           38 :          DO ispin = 1, nspins
     842          178 :             DO j = 1, curvy_data%n_bch_hist
     843              :                CALL dbcsr_create(curvy_data%matrix_BCH(ispin, j), template=matrix_s, &
     844          160 :                                  matrix_type=dbcsr_type_no_symmetry)
     845              :             END DO
     846              :          END DO
     847              :       END IF
     848              : 
     849           18 :    END SUBROUTINE init_curvy
     850              : 
     851              : END MODULE dm_ls_scf_curvy
        

Generated by: LCOV version 2.0-1