LCOV - code coverage report
Current view: top level - src - dm_ls_scf_curvy.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 340 349 97.4 %
Date: 2025-02-18 08:24:35 Functions: 13 13 100.0 %

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

Generated by: LCOV version 1.15