LCOV - code coverage report
Current view: top level - src - qs_ot_scf.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:4dc10b3) Lines: 113 151 74.8 %
Date: 2024-11-21 06:45:46 Functions: 4 4 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief basic functionality for using ot in the scf routines.
      10             : !> \par History
      11             : !>      01.2003 : Joost VandeVondele : adapted for LSD
      12             : !> \author Joost VandeVondele (25.08.2002)
      13             : ! **************************************************************************************************
      14             : MODULE qs_ot_scf
      15             :    USE cp_array_utils,                  ONLY: cp_1d_r_p_type
      16             :    USE cp_dbcsr_api,                    ONLY: &
      17             :         dbcsr_copy, dbcsr_dot, dbcsr_get_diag, dbcsr_get_info, dbcsr_init_p, dbcsr_multiply, &
      18             :         dbcsr_p_type, dbcsr_release, dbcsr_scale_by_vector, dbcsr_set, dbcsr_set_diag, dbcsr_type, &
      19             :         dbcsr_type_no_symmetry
      20             :    USE cp_dbcsr_operations,             ONLY: copy_fm_to_dbcsr,&
      21             :                                               cp_dbcsr_m_by_n_from_row_template
      22             :    USE cp_fm_types,                     ONLY: cp_fm_type
      23             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      24             :                                               cp_logger_type
      25             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      26             :                                               cp_print_key_unit_nr
      27             :    USE input_section_types,             ONLY: section_vals_get,&
      28             :                                               section_vals_get_subs_vals,&
      29             :                                               section_vals_type
      30             :    USE kinds,                           ONLY: dp
      31             :    USE qs_mo_occupation,                ONLY: set_mo_occupation
      32             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      33             :                                               mo_set_restrict,&
      34             :                                               mo_set_type
      35             :    USE qs_ot,                           ONLY: qs_ot_get_orbitals,&
      36             :                                               qs_ot_get_orbitals_ref,&
      37             :                                               qs_ot_get_p
      38             :    USE qs_ot_minimizer,                 ONLY: ot_mini
      39             :    USE qs_ot_types,                     ONLY: ot_readwrite_input,&
      40             :                                               qs_ot_allocate,&
      41             :                                               qs_ot_destroy,&
      42             :                                               qs_ot_init,&
      43             :                                               qs_ot_settings_init,&
      44             :                                               qs_ot_type
      45             :    USE scf_control_types,               ONLY: smear_type
      46             : #include "./base/base_uses.f90"
      47             : 
      48             :    IMPLICIT NONE
      49             : 
      50             :    PRIVATE
      51             : 
      52             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot_scf'
      53             :    ! *** Public subroutines ***
      54             : 
      55             :    PUBLIC :: ot_scf_init
      56             :    PUBLIC :: ot_scf_mini
      57             :    PUBLIC :: ot_scf_destroy
      58             :    PUBLIC :: ot_scf_read_input
      59             : 
      60             : CONTAINS
      61             : 
      62             : ! **************************************************************************************************
      63             : !> \brief ...
      64             : !> \param qs_ot_env ...
      65             : !> \param scf_section ...
      66             : ! **************************************************************************************************
      67       13148 :    SUBROUTINE ot_scf_read_input(qs_ot_env, scf_section)
      68             :       TYPE(qs_ot_type), DIMENSION(:), POINTER            :: qs_ot_env
      69             :       TYPE(section_vals_type), POINTER                   :: scf_section
      70             : 
      71             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ot_scf_read_input'
      72             : 
      73             :       INTEGER                                            :: handle, ispin, nspin, output_unit
      74             :       LOGICAL                                            :: explicit
      75             :       TYPE(cp_logger_type), POINTER                      :: logger
      76             :       TYPE(section_vals_type), POINTER                   :: ot_section
      77             : 
      78        6574 :       CALL timeset(routineN, handle)
      79             : 
      80        6574 :       logger => cp_get_default_logger()
      81             :       output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%PROGRAM_RUN_INFO", &
      82        6574 :                                          extension=".log")
      83             : 
      84             :       ! decide default settings
      85        6574 :       CALL qs_ot_settings_init(qs_ot_env(1)%settings)
      86             : 
      87             :       ! use ot input new style
      88        6574 :       ot_section => section_vals_get_subs_vals(scf_section, "OT")
      89        6574 :       CALL section_vals_get(ot_section, explicit=explicit)
      90             : 
      91        6574 :       CALL ot_readwrite_input(qs_ot_env(1)%settings, ot_section, output_unit)
      92             : 
      93             :       CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
      94        6574 :                                         "PRINT%PROGRAM_RUN_INFO")
      95             : 
      96             :       ! copy the ot settings type so it is identical
      97        6574 :       nspin = SIZE(qs_ot_env)
      98        7837 :       DO ispin = 2, nspin
      99        7837 :          qs_ot_env(ispin)%settings = qs_ot_env(1)%settings
     100             :       END DO
     101             : 
     102        6574 :       CALL timestop(handle)
     103             : 
     104        6574 :    END SUBROUTINE ot_scf_read_input
     105             : ! **************************************************************************************************
     106             :    !
     107             :    ! performs the actual minimisation, needs only limited info
     108             :    ! updated for restricted calculations
     109             :    ! matrix_dedc is the derivative of the energy with respect to the orbitals (except for a factor 2*fi)
     110             :    ! a null pointer for matrix_s implies that matrix_s is the unit matrix
     111             :    !
     112             :    !
     113             : ! **************************************************************************************************
     114             : !> \brief ...
     115             : !> \param mo_array ...
     116             : !> \param matrix_dedc ...
     117             : !> \param smear ...
     118             : !> \param matrix_s ...
     119             : !> \param energy ...
     120             : !> \param energy_only ...
     121             : !> \param delta ...
     122             : !> \param qs_ot_env ...
     123             : ! **************************************************************************************************
     124       62196 :    SUBROUTINE ot_scf_mini(mo_array, matrix_dedc, smear, matrix_s, energy, &
     125             :                           energy_only, delta, qs_ot_env)
     126             : 
     127             :       TYPE(mo_set_type), DIMENSION(:), INTENT(INOUT)     :: mo_array
     128             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_dedc
     129             :       TYPE(smear_type), POINTER                          :: smear
     130             :       TYPE(dbcsr_type), POINTER                          :: matrix_s
     131             :       REAL(KIND=dp)                                      :: energy
     132             :       LOGICAL, INTENT(INOUT)                             :: energy_only
     133             :       REAL(KIND=dp)                                      :: delta
     134             :       TYPE(qs_ot_type), DIMENSION(:), POINTER            :: qs_ot_env
     135             : 
     136             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ot_scf_mini'
     137             : 
     138             :       INTEGER                                            :: handle, ispin, k, n, nspin
     139             :       REAL(KIND=dp)                                      :: ener_nondiag, trace
     140       62196 :       TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:)    :: expectation_values, occupation_numbers, &
     141       62196 :                                                             scaling_factor
     142             :       TYPE(cp_logger_type), POINTER                      :: logger
     143       62196 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_dedc_scaled
     144             :       TYPE(dbcsr_type), POINTER                          :: mo_coeff
     145             : 
     146       62196 :       CALL timeset(routineN, handle)
     147             : 
     148       62196 :       NULLIFY (logger)
     149       62196 :       logger => cp_get_default_logger()
     150             : 
     151       62196 :       nspin = SIZE(mo_array)
     152             : 
     153      259268 :       ALLOCATE (occupation_numbers(nspin))
     154      197072 :       ALLOCATE (scaling_factor(nspin))
     155             : 
     156       62196 :       IF (qs_ot_env(1)%settings%do_ener) THEN
     157           0 :          ALLOCATE (expectation_values(nspin))
     158             :       END IF
     159             : 
     160      134876 :       DO ispin = 1, nspin
     161       72680 :          CALL get_mo_set(mo_set=mo_array(ispin), occupation_numbers=occupation_numbers(ispin)%array)
     162      217468 :          ALLOCATE (scaling_factor(ispin)%array(SIZE(occupation_numbers(ispin)%array)))
     163      847840 :          scaling_factor(ispin)%array = 2.0_dp*occupation_numbers(ispin)%array
     164      134876 :          IF (qs_ot_env(1)%settings%do_ener) THEN
     165           0 :             ALLOCATE (expectation_values(ispin)%array(SIZE(occupation_numbers(ispin)%array)))
     166             :          END IF
     167             :       END DO
     168             : 
     169             :       ! optimizing orbital energies somehow implies non-equivalent orbitals
     170       62196 :       IF (qs_ot_env(1)%settings%do_ener) THEN
     171           0 :          CPASSERT(qs_ot_env(1)%settings%do_rotation)
     172             :       END IF
     173             :       ! add_nondiag_energy requires do_ener
     174       62196 :       IF (qs_ot_env(1)%settings%add_nondiag_energy) THEN
     175           0 :          CPASSERT(qs_ot_env(1)%settings%do_ener)
     176             :       END IF
     177             : 
     178             :       ! get a rotational force
     179       62196 :       IF (.NOT. energy_only) THEN
     180       51803 :          IF (qs_ot_env(1)%settings%do_rotation) THEN
     181        2724 :             DO ispin = 1, SIZE(qs_ot_env)
     182        1388 :                CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
     183        1388 :                CALL dbcsr_get_info(mo_coeff, nfullrows_total=n, nfullcols_total=k)
     184             :                CALL dbcsr_multiply('T', 'N', 1.0_dp, mo_coeff, matrix_dedc(ispin)%matrix, &
     185        1388 :                                    0.0_dp, qs_ot_env(ispin)%rot_mat_chc)
     186        1388 :                CALL dbcsr_copy(qs_ot_env(ispin)%matrix_buf1, qs_ot_env(ispin)%rot_mat_chc)
     187             : 
     188        1388 :                CALL dbcsr_scale_by_vector(qs_ot_env(ispin)%matrix_buf1, alpha=scaling_factor(ispin)%array, side='right')
     189             :                ! create the derivative of the energy wrt to rot_mat_u
     190             :                CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env(ispin)%rot_mat_u, qs_ot_env(ispin)%matrix_buf1, &
     191        4112 :                                    0.0_dp, qs_ot_env(ispin)%rot_mat_dedu)
     192             :             END DO
     193             : 
     194             :             ! here we construct the derivative of the free energy with respect to the evals
     195             :             ! (note that this requires the diagonal elements of chc)
     196             :             ! the mo occupations should in principle remain unaltered
     197        1336 :             IF (qs_ot_env(1)%settings%do_ener) THEN
     198           0 :                DO ispin = 1, SIZE(mo_array)
     199           0 :                   CALL dbcsr_get_diag(qs_ot_env(ispin)%rot_mat_chc, expectation_values(ispin)%array)
     200           0 :                   qs_ot_env(ispin)%ener_gx = expectation_values(ispin)%array
     201             :                   CALL set_mo_occupation(mo_set=mo_array(ispin), &
     202           0 :                                          smear=smear, eval_deriv=qs_ot_env(ispin)%ener_gx)
     203             :                END DO
     204             :             END IF
     205             : 
     206             :             ! chc only needs to be stored in u independent form if we require add_nondiag_energy,
     207             :             ! which will use it in non-selfconsistent form for e.g. the linesearch
     208             :             ! transform C^T H C -> U C^T H C U ^ T
     209        1336 :             IF (qs_ot_env(1)%settings%add_nondiag_energy) THEN
     210           0 :                DO ispin = 1, SIZE(qs_ot_env)
     211           0 :                   CALL dbcsr_get_info(qs_ot_env(ispin)%rot_mat_u, nfullcols_total=k)
     212             :                   CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env(ispin)%rot_mat_u, qs_ot_env(ispin)%rot_mat_chc, &
     213           0 :                                       0.0_dp, qs_ot_env(ispin)%matrix_buf1)
     214             :                   CALL dbcsr_multiply('N', 'T', 1.0_dp, qs_ot_env(ispin)%matrix_buf1, qs_ot_env(ispin)%rot_mat_u, &
     215           0 :                                       0.0_dp, qs_ot_env(ispin)%rot_mat_chc)
     216             :                END DO
     217             :             END IF
     218             :          END IF
     219             :       END IF
     220             : 
     221             :       ! evaluate non-diagonal energy contribution
     222       62196 :       ener_nondiag = 0.0_dp
     223       62196 :       IF (qs_ot_env(1)%settings%add_nondiag_energy) THEN
     224           0 :          DO ispin = 1, SIZE(qs_ot_env)
     225             :             ! transform \tilde H to the current basis of C (assuming non-selfconsistent H)
     226           0 :             CALL dbcsr_get_info(qs_ot_env(ispin)%rot_mat_u, nfullcols_total=k)
     227             :             CALL dbcsr_multiply('T', 'N', 1.0_dp, qs_ot_env(ispin)%rot_mat_u, qs_ot_env(ispin)%rot_mat_chc, &
     228           0 :                                 0.0_dp, qs_ot_env(ispin)%matrix_buf1)
     229             :             CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env(ispin)%matrix_buf1, qs_ot_env(ispin)%rot_mat_u, &
     230           0 :                                 0.0_dp, qs_ot_env(ispin)%matrix_buf2)
     231             : 
     232             :             ! subtract the current ener_x from the diagonal
     233           0 :             CALL dbcsr_get_diag(qs_ot_env(ispin)%matrix_buf2, expectation_values(ispin)%array)
     234           0 :             expectation_values(ispin)%array = expectation_values(ispin)%array - qs_ot_env(ispin)%ener_x
     235           0 :             CALL dbcsr_set_diag(qs_ot_env(ispin)%matrix_buf2, expectation_values(ispin)%array)
     236             : 
     237             :             ! get nondiag energy trace (D^T D)
     238           0 :             CALL dbcsr_dot(qs_ot_env(ispin)%matrix_buf2, qs_ot_env(ispin)%matrix_buf2, trace)
     239           0 :             ener_nondiag = ener_nondiag + 0.5_dp*qs_ot_env(1)%settings%nondiag_energy_strength*trace
     240             : 
     241             :             ! get gradient (again ignoring dependencies of H)
     242           0 :             IF (.NOT. energy_only) THEN
     243             :                ! first for the ener_x (-2*(diag(C^T H C)-ener_x))
     244             :                qs_ot_env(ispin)%ener_gx = qs_ot_env(ispin)%ener_gx - &
     245           0 :                                           qs_ot_env(1)%settings%nondiag_energy_strength*expectation_values(ispin)%array
     246             : 
     247             :                ! next for the rot_mat_u derivative (2 * k * \tilde H U D)
     248             :                CALL dbcsr_multiply('N', 'N', 1.0_dp, qs_ot_env(ispin)%rot_mat_chc, qs_ot_env(ispin)%rot_mat_u, &
     249           0 :                                    0.0_dp, qs_ot_env(ispin)%matrix_buf1)
     250             :                CALL dbcsr_multiply('N', 'N', 2.0_dp*qs_ot_env(1)%settings%nondiag_energy_strength, &
     251             :                                    qs_ot_env(ispin)%matrix_buf1, qs_ot_env(ispin)%matrix_buf2, &
     252           0 :                                    1.0_dp, qs_ot_env(ispin)%rot_mat_dedu)
     253             :             END IF
     254             :          END DO
     255             :       END IF
     256             : 
     257             :       ! this is kind of a hack so far (costly memory wise), we locally recreate the scaled matrix_hc, and
     258             :       ! use it in the following, eventually, as occupations numbers get more integrated, it should become possible
     259             :       ! to remove this.
     260      258630 :       ALLOCATE (matrix_dedc_scaled(SIZE(matrix_dedc)))
     261      134238 :       DO ispin = 1, SIZE(matrix_dedc)
     262       72042 :          ALLOCATE (matrix_dedc_scaled(ispin)%matrix)
     263       72042 :          CALL dbcsr_copy(matrix_dedc_scaled(ispin)%matrix, matrix_dedc(ispin)%matrix)
     264             : 
     265             :          ! as a preconditioner, one might want to scale only with a constant, not with f(i)
     266             :          ! for the convergence criterion, maybe take it back out
     267       72042 :          IF (qs_ot_env(1)%settings%occupation_preconditioner) THEN
     268           0 :             scaling_factor(ispin)%array = 2.0_dp
     269             :          END IF
     270      134238 :          CALL dbcsr_scale_by_vector(matrix_dedc_scaled(ispin)%matrix, alpha=scaling_factor(ispin)%array, side='right')
     271             :       END DO
     272             : 
     273             :       ! notice we use qs_ot_env(1) for driving all output and the minimization in case of LSD
     274       62196 :       qs_ot_env(1)%etotal = energy + ener_nondiag
     275             : 
     276       62196 :       CALL ot_mini(qs_ot_env, matrix_dedc_scaled)
     277             : 
     278       62196 :       delta = qs_ot_env(1)%delta
     279       62196 :       energy_only = qs_ot_env(1)%energy_only
     280             : 
     281             :       ! generate the orbitals using the new matrix_x
     282      134238 :       DO ispin = 1, SIZE(qs_ot_env)
     283       72042 :          CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff)
     284       72042 :          CALL dbcsr_get_info(mo_coeff, nfullrows_total=n, nfullcols_total=k)
     285      134238 :          SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
     286             :          CASE ("TOD")
     287       67670 :             IF (ASSOCIATED(matrix_s)) THEN
     288             :                CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(ispin)%matrix_x, &
     289       50360 :                                    0.0_dp, qs_ot_env(ispin)%matrix_sx)
     290             :             ELSE
     291       17310 :                CALL dbcsr_copy(qs_ot_env(ispin)%matrix_sx, qs_ot_env(ispin)%matrix_x)
     292             :             END IF
     293       67670 :             CALL qs_ot_get_p(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_sx, qs_ot_env(ispin))
     294       67670 :             CALL qs_ot_get_orbitals(mo_coeff, qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin))
     295             :          CASE ("REF")
     296             :             CALL qs_ot_get_orbitals_ref(mo_coeff, matrix_s, qs_ot_env(ispin)%matrix_x, &
     297             :                                         qs_ot_env(ispin)%matrix_sx, qs_ot_env(ispin)%matrix_gx_old, &
     298        4372 :                                         qs_ot_env(ispin)%matrix_dx, qs_ot_env(ispin), qs_ot_env(1))
     299             :          CASE DEFAULT
     300       72042 :             CPABORT("Algorithm not yet implemented")
     301             :          END SELECT
     302             :       END DO
     303             : 
     304       62196 :       IF (qs_ot_env(1)%restricted) THEN
     305         638 :          CALL mo_set_restrict(mo_array, convert_dbcsr=.TRUE.)
     306             :       END IF
     307             :       !
     308             :       ! obtain the new set of OT eigenvalues and set the occupations accordingly
     309             :       !
     310       62196 :       IF (qs_ot_env(1)%settings%do_ener) THEN
     311           0 :          DO ispin = 1, SIZE(mo_array)
     312           0 :             mo_array(ispin)%eigenvalues = qs_ot_env(ispin)%ener_x
     313             :             CALL set_mo_occupation(mo_set=mo_array(ispin), &
     314           0 :                                    smear=smear)
     315             :          END DO
     316             :       END IF
     317             : 
     318             :       ! cleanup
     319      134876 :       DO ispin = 1, SIZE(scaling_factor)
     320      134876 :          DEALLOCATE (scaling_factor(ispin)%array)
     321             :       END DO
     322       62196 :       DEALLOCATE (scaling_factor)
     323       62196 :       IF (qs_ot_env(1)%settings%do_ener) THEN
     324           0 :          DO ispin = 1, SIZE(expectation_values)
     325           0 :             DEALLOCATE (expectation_values(ispin)%array)
     326             :          END DO
     327           0 :          DEALLOCATE (expectation_values)
     328             :       END IF
     329       62196 :       DEALLOCATE (occupation_numbers)
     330      134238 :       DO ispin = 1, SIZE(matrix_dedc_scaled)
     331       72042 :          CALL dbcsr_release(matrix_dedc_scaled(ispin)%matrix)
     332      134238 :          DEALLOCATE (matrix_dedc_scaled(ispin)%matrix)
     333             :       END DO
     334       62196 :       DEALLOCATE (matrix_dedc_scaled)
     335             : 
     336       62196 :       CALL timestop(handle)
     337             : 
     338      124392 :    END SUBROUTINE ot_scf_mini
     339             :    !
     340             :    ! initialises qs_ot_env so that mo_coeff is the current point
     341             :    ! and that the mimizization can be started.
     342             :    !
     343             : ! **************************************************************************************************
     344             : !> \brief ...
     345             : !> \param mo_array ...
     346             : !> \param matrix_s ...
     347             : !> \param qs_ot_env ...
     348             : !> \param matrix_ks ...
     349             : !> \param broyden_adaptive_sigma ...
     350             : ! **************************************************************************************************
     351        6574 :    SUBROUTINE ot_scf_init(mo_array, matrix_s, qs_ot_env, matrix_ks, broyden_adaptive_sigma)
     352             : 
     353             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mo_array
     354             :       TYPE(dbcsr_type), POINTER                          :: matrix_s
     355             :       TYPE(qs_ot_type), DIMENSION(:), POINTER            :: qs_ot_env
     356             :       TYPE(dbcsr_type), POINTER                          :: matrix_ks
     357             :       REAL(KIND=dp)                                      :: broyden_adaptive_sigma
     358             : 
     359             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'ot_scf_init'
     360             : 
     361             :       INTEGER                                            :: handle, ispin, k, n, nspin
     362             :       LOGICAL                                            :: is_equal
     363             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff_fm
     364             :       TYPE(dbcsr_type), POINTER                          :: mo_coeff
     365             : 
     366        6574 :       CALL timeset(routineN, handle)
     367             : 
     368       14503 :       DO ispin = 1, SIZE(mo_array)
     369       14503 :          IF (.NOT. ASSOCIATED(mo_array(ispin)%mo_coeff_b)) THEN
     370           0 :             CPABORT("Shouldn't get there")
     371             :             ! we do ot then copy fm to dbcsr
     372             :             ! allocate that somewhere else ! fm -> dbcsr
     373           0 :             CALL dbcsr_init_p(mo_array(ispin)%mo_coeff_b)
     374             :             CALL cp_dbcsr_m_by_n_from_row_template(mo_array(ispin)%mo_coeff_b, template=matrix_ks, &
     375             :                                                    n=mo_array(ispin)%nmo, &
     376           0 :                                                    sym=dbcsr_type_no_symmetry)
     377             :          END IF
     378             :       END DO
     379             : 
     380             :       ! *** set a history for broyden
     381       14411 :       DO ispin = 1, SIZE(qs_ot_env)
     382       14411 :          qs_ot_env(ispin)%broyden_adaptive_sigma = broyden_adaptive_sigma
     383             :       END DO
     384             : 
     385             :       ! **** SCP
     386             :       ! **** SCP
     387             :       ! adapted for work with the restricted keyword
     388        6574 :       nspin = SIZE(qs_ot_env)
     389             : 
     390       14411 :       DO ispin = 1, nspin
     391             : 
     392        7837 :          NULLIFY (mo_coeff)
     393        7837 :          CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff_b=mo_coeff, mo_coeff=mo_coeff_fm)
     394        7837 :          CALL copy_fm_to_dbcsr(mo_coeff_fm, mo_coeff) !fm -> dbcsr
     395             : 
     396        7837 :          CALL dbcsr_get_info(mo_coeff, nfullrows_total=n, nfullcols_total=k)
     397             : 
     398             :          ! allocate
     399        7837 :          CALL qs_ot_allocate(qs_ot_env(ispin), matrix_ks, mo_coeff_fm%matrix_struct)
     400             : 
     401             :          ! set c0,sc0
     402        7837 :          CALL dbcsr_copy(qs_ot_env(ispin)%matrix_c0, mo_coeff)
     403        7837 :          IF (ASSOCIATED(matrix_s)) THEN
     404             :             CALL dbcsr_multiply('N', 'N', 1.0_dp, matrix_s, qs_ot_env(ispin)%matrix_c0, &
     405        6627 :                                 0.0_dp, qs_ot_env(ispin)%matrix_sc0)
     406             :          ELSE
     407        1210 :             CALL dbcsr_copy(qs_ot_env(ispin)%matrix_sc0, qs_ot_env(ispin)%matrix_c0)
     408             :          END IF
     409             : 
     410             :          ! init
     411        7837 :          CALL qs_ot_init(qs_ot_env(ispin))
     412             : 
     413             :          ! set x
     414        7837 :          CALL dbcsr_set(qs_ot_env(ispin)%matrix_x, 0.0_dp)
     415        7837 :          CALL dbcsr_set(qs_ot_env(ispin)%matrix_sx, 0.0_dp)
     416             : 
     417        7837 :          IF (qs_ot_env(ispin)%settings%do_rotation) THEN
     418         224 :             CALL dbcsr_set(qs_ot_env(ispin)%rot_mat_x, 0.0_dp)
     419             :          END IF
     420             : 
     421        7837 :          IF (qs_ot_env(ispin)%settings%do_ener) THEN
     422           0 :             is_equal = SIZE(qs_ot_env(ispin)%ener_x) == SIZE(mo_array(ispin)%eigenvalues)
     423           0 :             CPASSERT(is_equal)
     424           0 :             qs_ot_env(ispin)%ener_x = mo_array(ispin)%eigenvalues
     425             :          END IF
     426             : 
     427       14411 :          SELECT CASE (qs_ot_env(1)%settings%ot_algorithm)
     428             :          CASE ("TOD")
     429             :             ! get c
     430        6761 :             CALL qs_ot_get_p(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_sx, qs_ot_env(ispin))
     431             :          CASE ("REF")
     432        1076 :             CALL dbcsr_copy(qs_ot_env(ispin)%matrix_x, qs_ot_env(ispin)%matrix_c0)
     433        1076 :             CALL dbcsr_copy(qs_ot_env(ispin)%matrix_sx, qs_ot_env(ispin)%matrix_sc0)
     434             :          CASE DEFAULT
     435        7837 :             CPABORT("Algorithm not yet implemented")
     436             :          END SELECT
     437             : 
     438             :       END DO
     439        6574 :       CALL timestop(handle)
     440        6574 :    END SUBROUTINE ot_scf_init
     441             : 
     442             : ! **************************************************************************************************
     443             : !> \brief ...
     444             : !> \param qs_ot_env ...
     445             : ! **************************************************************************************************
     446        7837 :    SUBROUTINE ot_scf_destroy(qs_ot_env)
     447             : 
     448             :       TYPE(qs_ot_type)                                   :: qs_ot_env
     449             : 
     450        7837 :       CALL qs_ot_destroy(qs_ot_env)
     451             : 
     452        7837 :    END SUBROUTINE ot_scf_destroy
     453             : 
     454             : END MODULE qs_ot_scf
     455             : 

Generated by: LCOV version 1.15