LCOV - code coverage report
Current view: top level - src - qs_ot_types.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:ccc2433) Lines: 434 470 92.3 %
Date: 2024-04-25 07:09:54 Functions: 5 7 71.4 %

          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 orbital transformations
      10             : !> \par History
      11             : !>      Added Taylor expansion based computation of the matrix functions (01.2004)
      12             : !>      added additional rotation variables for non-equivalent occupied orbs (08.2004)
      13             : !> \author Joost VandeVondele (06.2002)
      14             : ! **************************************************************************************************
      15             : MODULE qs_ot_types
      16             :    USE bibliography,                    ONLY: VandeVondele2003,&
      17             :                                               Weber2008,&
      18             :                                               cite_reference
      19             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_release,&
      20             :                                               cp_blacs_env_type
      21             :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_m_by_n_from_row_template,&
      22             :                                               cp_dbcsr_m_by_n_from_template,&
      23             :                                               dbcsr_allocate_matrix_set,&
      24             :                                               dbcsr_deallocate_matrix_set
      25             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_get,&
      26             :                                               cp_fm_struct_type
      27             :    USE dbcsr_api,                       ONLY: dbcsr_init_p,&
      28             :                                               dbcsr_p_type,&
      29             :                                               dbcsr_release_p,&
      30             :                                               dbcsr_set,&
      31             :                                               dbcsr_type,&
      32             :                                               dbcsr_type_complex_default,&
      33             :                                               dbcsr_type_no_symmetry,&
      34             :                                               dbcsr_type_real_8
      35             :    USE input_constants,                 ONLY: &
      36             :         cholesky_reduce, ls_2pnt, ls_3pnt, ls_gold, ls_none, ot_algo_irac, ot_algo_taylor_or_diag, &
      37             :         ot_chol_irac, ot_lwdn_irac, ot_mini_broyden, ot_mini_cg, ot_mini_diis, ot_mini_sd, &
      38             :         ot_poly_irac, ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, &
      39             :         ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, &
      40             :         ot_precond_solver_default, ot_precond_solver_direct, ot_precond_solver_inv_chol, &
      41             :         ot_precond_solver_update
      42             :    USE input_section_types,             ONLY: section_vals_type,&
      43             :                                               section_vals_val_get
      44             :    USE kinds,                           ONLY: dp
      45             :    USE message_passing,                 ONLY: mp_para_env_release,&
      46             :                                               mp_para_env_type
      47             :    USE preconditioner_types,            ONLY: preconditioner_type
      48             : #include "./base/base_uses.f90"
      49             : 
      50             :    IMPLICIT NONE
      51             : 
      52             :    PRIVATE
      53             : 
      54             :    PUBLIC  :: qs_ot_type
      55             :    PUBLIC  :: qs_ot_settings_type
      56             :    PUBLIC  :: qs_ot_destroy
      57             :    PUBLIC  :: qs_ot_allocate
      58             :    PUBLIC  :: qs_ot_init
      59             :    PUBLIC  :: qs_ot_settings_init
      60             :    PUBLIC  :: ot_readwrite_input
      61             : 
      62             :    ! notice, this variable needs to be copyable !
      63             :    ! needed for spins as e.g. in qs_ot_scf      !
      64             : ! **************************************************************************************************
      65             :    TYPE qs_ot_settings_type
      66             :       LOGICAL           :: do_rotation, do_ener
      67             :       LOGICAL           :: ks
      68             :       CHARACTER(LEN=4)  :: ot_method
      69             :       CHARACTER(LEN=3)  :: ot_algorithm
      70             :       CHARACTER(LEN=4)  :: line_search_method
      71             :       CHARACTER(LEN=20) :: preconditioner_name
      72             :       INTEGER           :: preconditioner_type
      73             :       INTEGER           :: cholesky_type
      74             :       CHARACTER(LEN=20) :: precond_solver_name
      75             :       INTEGER           :: precond_solver_type
      76             :       LOGICAL           :: safer_diis
      77             :       REAL(KIND=dp)   :: ds_min
      78             :       REAL(KIND=dp)   :: energy_gap
      79             :       INTEGER           :: diis_m
      80             :       REAL(KIND=dp)   :: gold_target
      81             :       REAL(KIND=dp)   :: eps_taylor ! minimum accuracy of the taylor expansion
      82             :       INTEGER           :: max_taylor ! maximum order of the taylor expansion before switching to diagonalization
      83             :       INTEGER           :: irac_degree ! this is used to control the refinement polynomial degree
      84             :       INTEGER           :: max_irac ! maximum number of iteration for the refinement
      85             :       REAL(KIND=dp)   :: eps_irac ! target accuracy for the refinement
      86             :       REAL(KIND=dp)   :: eps_irac_quick_exit
      87             :       REAL(dp)          :: eps_irac_filter_matrix
      88             :       REAL(KIND=dp)   :: eps_irac_switch
      89             :       LOGICAL           :: on_the_fly_loc
      90             :       CHARACTER(LEN=4)  :: ortho_irac
      91             :       LOGICAL           :: occupation_preconditioner, add_nondiag_energy
      92             :       REAL(KIND=dp)   :: nondiag_energy_strength
      93             :       REAL(KIND=dp)   :: broyden_beta, broyden_gamma, broyden_sigma
      94             :       REAL(KIND=dp)   :: broyden_eta, broyden_omega, broyden_sigma_decrease
      95             :       REAL(KIND=dp)   :: broyden_sigma_min
      96             :       LOGICAL           :: broyden_forget_history, broyden_adaptive_sigma
      97             :       LOGICAL           :: broyden_enable_flip
      98             :    END TYPE qs_ot_settings_type
      99             : 
     100             : ! **************************************************************************************************
     101             :    TYPE qs_ot_type
     102             :       ! this sets the method to be used
     103             :       TYPE(qs_ot_settings_type) :: settings
     104             :       LOGICAL                   :: restricted
     105             : 
     106             :       ! first part of the variables, for occupied subspace invariant optimisation
     107             : 
     108             :       ! add a preconditioner matrix. should be symmetric and positive definite
     109             :       ! the type of this matrix might change in the future
     110             :       TYPE(preconditioner_type), POINTER :: preconditioner
     111             : 
     112             :       ! these will/might change during iterations
     113             : 
     114             :       ! OT / TOD
     115             :       TYPE(dbcsr_type), POINTER :: matrix_p
     116             :       TYPE(dbcsr_type), POINTER :: matrix_r
     117             :       TYPE(dbcsr_type), POINTER :: matrix_sinp
     118             :       TYPE(dbcsr_type), POINTER :: matrix_cosp
     119             :       TYPE(dbcsr_type), POINTER :: matrix_sinp_b
     120             :       TYPE(dbcsr_type), POINTER :: matrix_cosp_b
     121             :       TYPE(dbcsr_type), POINTER :: matrix_buf1
     122             :       TYPE(dbcsr_type), POINTER :: matrix_buf2
     123             :       TYPE(dbcsr_type), POINTER :: matrix_buf3
     124             :       TYPE(dbcsr_type), POINTER :: matrix_buf4
     125             :       TYPE(dbcsr_type), POINTER :: matrix_os
     126             :       TYPE(dbcsr_type), POINTER :: matrix_buf1_ortho
     127             :       TYPE(dbcsr_type), POINTER :: matrix_buf2_ortho
     128             : 
     129             :       REAL(KIND=dp), DIMENSION(:), POINTER :: evals
     130             :       REAL(KIND=dp), DIMENSION(:), POINTER :: dum
     131             : 
     132             :       ! matrix os valid
     133             :       LOGICAL os_valid
     134             : 
     135             :       ! for efficient/parallel writing to the blacs_matrix
     136             :       TYPE(mp_para_env_type), POINTER :: para_env
     137             :       TYPE(cp_blacs_env_type), POINTER :: blacs_env
     138             : 
     139             :       ! mo-like vectors
     140             :       TYPE(dbcsr_type), POINTER :: matrix_c0, matrix_sc0, matrix_psc0
     141             : 
     142             :       ! OT / IR
     143             :       TYPE(dbcsr_type), POINTER :: buf1_k_k_nosym, buf2_k_k_nosym, &
     144             :                                    buf3_k_k_nosym, buf4_k_k_nosym, &
     145             :                                    buf1_k_k_sym, buf2_k_k_sym, buf3_k_k_sym, buf4_k_k_sym, &
     146             :                                    p_k_k_sym, buf1_n_k, buf1_n_k_dp
     147             : 
     148             :       ! only here for the ease of programming. These will have to be supplied
     149             :       ! explicitly at all times
     150             :       TYPE(dbcsr_type), POINTER :: matrix_x, matrix_sx, matrix_gx
     151             :       TYPE(dbcsr_type), POINTER :: matrix_dx, matrix_gx_old
     152             : 
     153             :       LOGICAL :: use_gx_old, use_dx
     154             : 
     155             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h_e, matrix_h_x
     156             : 
     157             :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: ls_diis
     158             :       REAL(KIND=dp), DIMENSION(:, :), POINTER  :: lss_diis
     159             :       REAL(KIND=dp), DIMENSION(:), POINTER  :: c_diis
     160             :       REAL(KIND=dp), DIMENSION(:), POINTER  :: c_broy
     161             :       REAL(KIND=dp), DIMENSION(:), POINTER  :: energy_h
     162             :       INTEGER, DIMENSION(:), POINTER         :: ipivot
     163             : 
     164             :       REAL(KIND=dp)  :: ot_pos(53), ot_energy(53), ot_grad(53) ! HARD LIMIT FOR THE LS
     165             :       INTEGER          :: line_search_left, line_search_right, line_search_mid
     166             :       INTEGER          :: line_search_count
     167             :       LOGICAL          :: line_search_might_be_done
     168             :       REAL(KIND=dp)  :: delta, gnorm, gnorm_old, etotal, gradient
     169             :       LOGICAL          :: energy_only
     170             :       INTEGER          :: diis_iter
     171             :       CHARACTER(LEN=8) :: OT_METHOD_FULL
     172             :       INTEGER          :: OT_count
     173             :       REAL(KIND=dp)  :: ds_min
     174             :       REAL(KIND=dp)  :: broyden_adaptive_sigma
     175             : 
     176             :       LOGICAL          :: do_taylor
     177             :       INTEGER          :: taylor_order
     178             :       REAL(KIND=dp)  :: largest_eval_upper_bound
     179             : 
     180             :       ! second part of the variables, if an explicit rotation is required as well
     181             :       TYPE(dbcsr_type), POINTER :: rot_mat_u ! rotation matrix
     182             :       TYPE(dbcsr_type), POINTER :: rot_mat_x ! antisymmetric matrix that parametrises rot_matrix_u
     183             :       TYPE(dbcsr_type), POINTER :: rot_mat_dedu ! derivative of the total energy wrt to u
     184             :       TYPE(dbcsr_type), POINTER :: rot_mat_chc ! for convencience, the matrix c^T H c
     185             : 
     186             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rot_mat_h_e
     187             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rot_mat_h_x
     188             :       TYPE(dbcsr_type), POINTER :: rot_mat_gx
     189             :       TYPE(dbcsr_type), POINTER :: rot_mat_gx_old
     190             :       TYPE(dbcsr_type), POINTER :: rot_mat_dx
     191             : 
     192             :       REAL(KIND=dp), DIMENSION(:), POINTER :: rot_mat_evals
     193             :       TYPE(dbcsr_type), POINTER :: rot_mat_evec
     194             : 
     195             :       ! third part of the variables, if we need to optimize orbital energies
     196             :       REAL(KIND=dp), POINTER, DIMENSION(:) :: ener_x
     197             :       REAL(KIND=dp), POINTER, DIMENSION(:) :: ener_dx
     198             :       REAL(KIND=dp), POINTER, DIMENSION(:) :: ener_gx
     199             :       REAL(KIND=dp), POINTER, DIMENSION(:) :: ener_gx_old
     200             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: ener_h_e
     201             :       REAL(KIND=dp), POINTER, DIMENSION(:, :) :: ener_h_x
     202             :    END TYPE qs_ot_type
     203             : 
     204             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ot_types'
     205             : 
     206             : CONTAINS
     207             : 
     208             : ! **************************************************************************************************
     209             : !> \brief sets default values for the settings type
     210             : !> \param settings ...
     211             : !> \par History
     212             : !>      10.2004 created [Joost VandeVondele]
     213             : ! **************************************************************************************************
     214       13830 :    SUBROUTINE qs_ot_settings_init(settings)
     215             :       TYPE(qs_ot_settings_type)                          :: settings
     216             : 
     217       13830 :       settings%ot_method = "CG"
     218       13830 :       settings%ot_algorithm = "TOD"
     219       13830 :       settings%diis_m = 7
     220       13830 :       settings%preconditioner_name = "FULL_KINETIC"
     221       13830 :       settings%preconditioner_type = ot_precond_full_kinetic
     222       13830 :       settings%cholesky_type = cholesky_reduce
     223       13830 :       settings%precond_solver_name = "CHOLESKY_INVERSE"
     224       13830 :       settings%precond_solver_type = ot_precond_solver_inv_chol
     225       13830 :       settings%line_search_method = "2PNT"
     226       13830 :       settings%ds_min = 0.15_dp
     227       13830 :       settings%safer_diis = .TRUE.
     228       13830 :       settings%energy_gap = 0.2_dp
     229       13830 :       settings%eps_taylor = 1.0E-16_dp
     230       13830 :       settings%max_taylor = 4
     231       13830 :       settings%gold_target = 0.01_dp
     232       13830 :       settings%do_rotation = .FALSE.
     233       13830 :       settings%do_ener = .FALSE.
     234       13830 :       settings%irac_degree = 4
     235       13830 :       settings%max_irac = 50
     236       13830 :       settings%eps_irac = 1.0E-10_dp
     237       13830 :       settings%eps_irac_quick_exit = 1.0E-5_dp
     238       13830 :       settings%eps_irac_switch = 1.0E-2
     239       13830 :       settings%eps_irac_filter_matrix = 0.0_dp
     240       13830 :       settings%on_the_fly_loc = .FALSE.
     241       13830 :       settings%ortho_irac = "CHOL"
     242       13830 :       settings%ks = .TRUE.
     243       13830 :       settings%occupation_preconditioner = .FALSE.
     244       13830 :       settings%add_nondiag_energy = .FALSE.
     245       13830 :       settings%nondiag_energy_strength = 0.0_dp
     246       13830 :    END SUBROUTINE qs_ot_settings_init
     247             : 
     248             :    ! init matrices, needs c0 and sc0 so that c0*sc0=1
     249             : ! **************************************************************************************************
     250             : !> \brief ...
     251             : !> \param qs_ot_env ...
     252             : ! **************************************************************************************************
     253        8561 :    SUBROUTINE qs_ot_init(qs_ot_env)
     254             :       TYPE(qs_ot_type)                                   :: qs_ot_env
     255             : 
     256      462294 :       qs_ot_env%OT_energy(:) = 0.0_dp
     257      462294 :       qs_ot_env%OT_pos(:) = 0.0_dp
     258      462294 :       qs_ot_env%OT_grad(:) = 0.0_dp
     259        8561 :       qs_ot_env%line_search_count = 0
     260             : 
     261        8561 :       qs_ot_env%energy_only = .FALSE.
     262        8561 :       qs_ot_env%gnorm_old = 1.0_dp
     263        8561 :       qs_ot_env%diis_iter = 0
     264        8561 :       qs_ot_env%ds_min = qs_ot_env%settings%ds_min
     265        8561 :       qs_ot_env%os_valid = .FALSE.
     266             : 
     267        8561 :       CALL dbcsr_set(qs_ot_env%matrix_gx, 0.0_dp)
     268             : 
     269        8561 :       IF (qs_ot_env%use_dx) &
     270        3930 :          CALL dbcsr_set(qs_ot_env%matrix_dx, 0.0_dp)
     271             : 
     272        8561 :       IF (qs_ot_env%use_gx_old) &
     273        3930 :          CALL dbcsr_set(qs_ot_env%matrix_gx_old, 0.0_dp)
     274             : 
     275        8561 :       IF (qs_ot_env%settings%do_rotation) THEN
     276         232 :          CALL dbcsr_set(qs_ot_env%rot_mat_gx, 0.0_dp)
     277         232 :          IF (qs_ot_env%use_dx) &
     278         202 :             CALL dbcsr_set(qs_ot_env%rot_mat_dx, 0.0_dp)
     279         232 :          IF (qs_ot_env%use_gx_old) &
     280         202 :             CALL dbcsr_set(qs_ot_env%rot_mat_gx_old, 0.0_dp)
     281             :       END IF
     282        8561 :       IF (qs_ot_env%settings%do_ener) THEN
     283           0 :          qs_ot_env%ener_gx(:) = 0.0_dp
     284           0 :          IF (qs_ot_env%use_dx) &
     285           0 :             qs_ot_env%ener_dx(:) = 0.0_dp
     286           0 :          IF (qs_ot_env%use_gx_old) &
     287           0 :             qs_ot_env%ener_gx_old(:) = 0.0_dp
     288             :       END IF
     289             : 
     290        8561 :    END SUBROUTINE qs_ot_init
     291             : 
     292             :    ! allocates the data in qs_ot_env, for a calculation with fm_struct_ref
     293             :    ! ortho_k allows for specifying an additional orthogonal subspace (i.e. c will
     294             :    ! be kept orthogonal provided c0 was, used in qs_ot_eigensolver)
     295             : ! **************************************************************************************************
     296             : !> \brief ...
     297             : !> \param qs_ot_env ...
     298             : !> \param matrix_s ...
     299             : !> \param fm_struct_ref ...
     300             : !> \param ortho_k ...
     301             : ! **************************************************************************************************
     302        8561 :    SUBROUTINE qs_ot_allocate(qs_ot_env, matrix_s, fm_struct_ref, ortho_k)
     303             :       TYPE(qs_ot_type)                                   :: qs_ot_env
     304             :       TYPE(dbcsr_type), POINTER                          :: matrix_s
     305             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_ref
     306             :       INTEGER, OPTIONAL                                  :: ortho_k
     307             : 
     308             :       INTEGER                                            :: i, k, m_diis, my_ortho_k, n, ncoef
     309             :       TYPE(cp_blacs_env_type), POINTER                   :: context
     310             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     311             : 
     312        8561 :       CALL cite_reference(VandeVondele2003)
     313             : 
     314        8561 :       NULLIFY (qs_ot_env%preconditioner)
     315        8561 :       NULLIFY (qs_ot_env%matrix_psc0)
     316        8561 :       NULLIFY (qs_ot_env%para_env)
     317        8561 :       NULLIFY (qs_ot_env%blacs_env)
     318             : 
     319             :       CALL cp_fm_struct_get(fm_struct_ref, nrow_global=n, ncol_global=k, &
     320        8561 :                             para_env=para_env, context=context)
     321             : 
     322        8561 :       qs_ot_env%para_env => para_env
     323        8561 :       qs_ot_env%blacs_env => context
     324        8561 :       CALL para_env%retain()
     325        8561 :       CALL context%retain()
     326             : 
     327        8561 :       IF (PRESENT(ortho_k)) THEN
     328         654 :          my_ortho_k = ortho_k
     329             :       ELSE
     330        7907 :          my_ortho_k = k
     331             :       END IF
     332             : 
     333        8561 :       m_diis = qs_ot_env%settings%diis_m
     334             : 
     335        8561 :       qs_ot_env%use_gx_old = .FALSE.
     336        8561 :       qs_ot_env%use_dx = .FALSE.
     337             : 
     338        3930 :       SELECT CASE (qs_ot_env%settings%ot_method)
     339             :       CASE ("SD")
     340             :          ! nothing
     341             :       CASE ("CG")
     342        3930 :          qs_ot_env%use_gx_old = .TRUE.
     343        3930 :          qs_ot_env%use_dx = .TRUE.
     344             :       CASE ("DIIS", "BROY")
     345        4621 :          IF (m_diis .LT. 1) CPABORT("m_diis less than one")
     346             :       CASE DEFAULT
     347        8561 :          CPABORT("Unknown option")
     348             :       END SELECT
     349             : 
     350        8561 :       IF (qs_ot_env%settings%ot_method .EQ. "DIIS" .OR. &
     351             :           qs_ot_env%settings%ot_method .EQ. "BROY") THEN
     352       18484 :          ALLOCATE (qs_ot_env%ls_diis(m_diis + 1, m_diis + 1))
     353      348543 :          qs_ot_env%ls_diis = 0.0_dp
     354       18484 :          ALLOCATE (qs_ot_env%lss_diis(m_diis + 1, m_diis + 1))
     355       13863 :          ALLOCATE (qs_ot_env%c_diis(m_diis + 1))
     356       13863 :          ALLOCATE (qs_ot_env%c_broy(m_diis))
     357       13863 :          ALLOCATE (qs_ot_env%energy_h(m_diis))
     358       13863 :          ALLOCATE (qs_ot_env%ipivot(m_diis + 1))
     359             :       END IF
     360             : 
     361       25427 :       ALLOCATE (qs_ot_env%evals(k))
     362       25427 :       ALLOCATE (qs_ot_env%dum(k))
     363             : 
     364        8561 :       NULLIFY (qs_ot_env%matrix_os)
     365        8561 :       NULLIFY (qs_ot_env%matrix_buf1_ortho)
     366        8561 :       NULLIFY (qs_ot_env%matrix_buf2_ortho)
     367        8561 :       NULLIFY (qs_ot_env%matrix_p)
     368        8561 :       NULLIFY (qs_ot_env%matrix_r)
     369        8561 :       NULLIFY (qs_ot_env%matrix_sinp)
     370        8561 :       NULLIFY (qs_ot_env%matrix_cosp)
     371        8561 :       NULLIFY (qs_ot_env%matrix_sinp_b)
     372        8561 :       NULLIFY (qs_ot_env%matrix_cosp_b)
     373        8561 :       NULLIFY (qs_ot_env%matrix_buf1)
     374        8561 :       NULLIFY (qs_ot_env%matrix_buf2)
     375        8561 :       NULLIFY (qs_ot_env%matrix_buf3)
     376        8561 :       NULLIFY (qs_ot_env%matrix_buf4)
     377        8561 :       NULLIFY (qs_ot_env%matrix_c0)
     378        8561 :       NULLIFY (qs_ot_env%matrix_sc0)
     379        8561 :       NULLIFY (qs_ot_env%matrix_x)
     380        8561 :       NULLIFY (qs_ot_env%matrix_sx)
     381        8561 :       NULLIFY (qs_ot_env%matrix_gx)
     382        8561 :       NULLIFY (qs_ot_env%matrix_gx_old)
     383        8561 :       NULLIFY (qs_ot_env%matrix_dx)
     384        8561 :       NULLIFY (qs_ot_env%buf1_k_k_nosym)
     385        8561 :       NULLIFY (qs_ot_env%buf2_k_k_nosym)
     386        8561 :       NULLIFY (qs_ot_env%buf3_k_k_nosym)
     387        8561 :       NULLIFY (qs_ot_env%buf4_k_k_nosym)
     388        8561 :       NULLIFY (qs_ot_env%buf1_k_k_sym)
     389        8561 :       NULLIFY (qs_ot_env%buf2_k_k_sym)
     390        8561 :       NULLIFY (qs_ot_env%buf3_k_k_sym)
     391        8561 :       NULLIFY (qs_ot_env%buf4_k_k_sym)
     392        8561 :       NULLIFY (qs_ot_env%buf1_n_k)
     393        8561 :       NULLIFY (qs_ot_env%buf1_n_k_dp)
     394        8561 :       NULLIFY (qs_ot_env%p_k_k_sym)
     395             : 
     396             :       ! COMMON MATRICES
     397        8561 :       CALL dbcsr_init_p(qs_ot_env%matrix_c0)
     398             :       CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_c0, template=matrix_s, n=k, &
     399        8561 :                                              sym=dbcsr_type_no_symmetry)
     400             : 
     401        8561 :       CALL dbcsr_init_p(qs_ot_env%matrix_sc0)
     402             :       CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_sc0, template=matrix_s, n=my_ortho_k, &
     403        8561 :                                              sym=dbcsr_type_no_symmetry)
     404             : 
     405        8561 :       CALL dbcsr_init_p(qs_ot_env%matrix_x)
     406             :       CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_x, template=matrix_s, n=k, &
     407        8561 :                                              sym=dbcsr_type_no_symmetry)
     408             : 
     409        8561 :       CALL dbcsr_init_p(qs_ot_env%matrix_sx)
     410             :       CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_sx, template=matrix_s, n=k, &
     411        8561 :                                              sym=dbcsr_type_no_symmetry)
     412             : 
     413        8561 :       CALL dbcsr_init_p(qs_ot_env%matrix_gx)
     414             :       CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_gx, template=matrix_s, n=k, &
     415        8561 :                                              sym=dbcsr_type_no_symmetry, data_type=dbcsr_type_real_8)
     416             : 
     417        8561 :       IF (qs_ot_env%use_dx) THEN
     418        3930 :          CALL dbcsr_init_p(qs_ot_env%matrix_dx)
     419             :          CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_dx, template=matrix_s, n=k, &
     420        3930 :                                                 sym=dbcsr_type_no_symmetry, data_type=dbcsr_type_real_8)
     421             :       END IF
     422             : 
     423        8561 :       IF (qs_ot_env%use_gx_old) THEN
     424        3930 :          CALL dbcsr_init_p(qs_ot_env%matrix_gx_old)
     425             :          CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_gx_old, template=matrix_s, n=k, &
     426        3930 :                                                 sym=dbcsr_type_no_symmetry, data_type=dbcsr_type_real_8)
     427             :       END IF
     428             : 
     429        7485 :       SELECT CASE (qs_ot_env%settings%ot_algorithm)
     430             :       CASE ("TOD")
     431        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_p)
     432             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_p, template=matrix_s, m=k, n=k, &
     433        7485 :                                             sym=dbcsr_type_no_symmetry)
     434             : 
     435        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_r)
     436             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_r, template=matrix_s, m=k, n=k, &
     437        7485 :                                             sym=dbcsr_type_no_symmetry)
     438             : 
     439        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_sinp)
     440             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_sinp, template=matrix_s, m=k, n=k, &
     441        7485 :                                             sym=dbcsr_type_no_symmetry)
     442             : 
     443        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_cosp)
     444             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_cosp, template=matrix_s, m=k, n=k, &
     445        7485 :                                             sym=dbcsr_type_no_symmetry)
     446             : 
     447        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_sinp_b)
     448             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_sinp_b, template=matrix_s, m=k, n=k, &
     449        7485 :                                             sym=dbcsr_type_no_symmetry)
     450             : 
     451        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_cosp_b)
     452             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_cosp_b, template=matrix_s, m=k, n=k, &
     453        7485 :                                             sym=dbcsr_type_no_symmetry)
     454             : 
     455        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_buf1)
     456             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf1, template=matrix_s, m=k, n=k, &
     457        7485 :                                             sym=dbcsr_type_no_symmetry)
     458             : 
     459        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_buf2)
     460             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf2, template=matrix_s, m=k, n=k, &
     461        7485 :                                             sym=dbcsr_type_no_symmetry)
     462             : 
     463        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_buf3)
     464             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf3, template=matrix_s, m=k, n=k, &
     465        7485 :                                             sym=dbcsr_type_no_symmetry)
     466             : 
     467        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_buf4)
     468             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf4, template=matrix_s, m=k, n=k, &
     469        7485 :                                             sym=dbcsr_type_no_symmetry)
     470             : 
     471        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_os)
     472             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_os, template=matrix_s, m=my_ortho_k, n=my_ortho_k, &
     473        7485 :                                             sym=dbcsr_type_no_symmetry)
     474             : 
     475        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_buf1_ortho)
     476             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf1_ortho, template=matrix_s, m=my_ortho_k, n=k, &
     477        7485 :                                             sym=dbcsr_type_no_symmetry)
     478             : 
     479        7485 :          CALL dbcsr_init_p(qs_ot_env%matrix_buf2_ortho)
     480             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf2_ortho, template=matrix_s, m=my_ortho_k, n=k, &
     481        7485 :                                             sym=dbcsr_type_no_symmetry)
     482             : 
     483             :       CASE ("REF")
     484        1076 :          CALL dbcsr_init_p(qs_ot_env%buf1_k_k_nosym)
     485             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf1_k_k_nosym, template=matrix_s, m=k, n=k, &
     486        1076 :                                             sym=dbcsr_type_no_symmetry)
     487             : 
     488        1076 :          CALL dbcsr_init_p(qs_ot_env%buf2_k_k_nosym)
     489             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf2_k_k_nosym, template=matrix_s, m=k, n=k, &
     490        1076 :                                             sym=dbcsr_type_no_symmetry)
     491             : 
     492        1076 :          CALL dbcsr_init_p(qs_ot_env%buf3_k_k_nosym)
     493             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf3_k_k_nosym, template=matrix_s, m=k, n=k, &
     494        1076 :                                             sym=dbcsr_type_no_symmetry)
     495             : 
     496        1076 :          CALL dbcsr_init_p(qs_ot_env%buf4_k_k_nosym)
     497             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf4_k_k_nosym, template=matrix_s, m=k, n=k, &
     498        1076 :                                             sym=dbcsr_type_no_symmetry)
     499             : 
     500             : ! It claims to be symmetric but to avoid dbcsr conusion nonsymmetric is kept
     501        1076 :          CALL dbcsr_init_p(qs_ot_env%buf1_k_k_sym)
     502             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf1_k_k_sym, template=matrix_s, m=k, n=k, &
     503        1076 :                                             sym=dbcsr_type_no_symmetry)
     504             : 
     505        1076 :          CALL dbcsr_init_p(qs_ot_env%buf2_k_k_sym)
     506             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf2_k_k_sym, template=matrix_s, m=k, n=k, &
     507        1076 :                                             sym=dbcsr_type_no_symmetry)
     508             : 
     509        1076 :          CALL dbcsr_init_p(qs_ot_env%buf3_k_k_sym)
     510             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf3_k_k_sym, template=matrix_s, m=k, n=k, &
     511        1076 :                                             sym=dbcsr_type_no_symmetry)
     512             :          !
     513        1076 :          CALL dbcsr_init_p(qs_ot_env%buf4_k_k_sym)
     514             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%buf4_k_k_sym, template=matrix_s, m=k, n=k, &
     515        1076 :                                             sym=dbcsr_type_no_symmetry)
     516             :          !
     517        1076 :          CALL dbcsr_init_p(qs_ot_env%p_k_k_sym)
     518             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%p_k_k_sym, template=matrix_s, m=k, n=k, &
     519        1076 :                                             sym=dbcsr_type_no_symmetry)
     520             :          !
     521        1076 :          CALL dbcsr_init_p(qs_ot_env%buf1_n_k)
     522             :          CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%buf1_n_k, template=matrix_s, n=k, &
     523        1076 :                                                 sym=dbcsr_type_no_symmetry)
     524             :          !
     525        1076 :          CALL dbcsr_init_p(qs_ot_env%matrix_buf1)
     526             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%matrix_buf1, template=matrix_s, m=k, n=k, &
     527        9637 :                                             sym=dbcsr_type_no_symmetry)
     528             : 
     529             :       END SELECT
     530             : 
     531        8561 :       IF (qs_ot_env%settings%ot_method .EQ. "DIIS" .OR. &
     532             :           qs_ot_env%settings%ot_method .EQ. "BROY") THEN
     533        4621 :          NULLIFY (qs_ot_env%matrix_h_e)
     534        4621 :          NULLIFY (qs_ot_env%matrix_h_x)
     535        4621 :          CALL dbcsr_allocate_matrix_set(qs_ot_env%matrix_h_e, m_diis)
     536        4621 :          CALL dbcsr_allocate_matrix_set(qs_ot_env%matrix_h_x, m_diis)
     537       37558 :          DO i = 1, m_diis
     538       32937 :             CALL dbcsr_init_p(qs_ot_env%matrix_h_x(i)%matrix)
     539             :             CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_h_x(i)%matrix, template=matrix_s, n=k, &
     540       32937 :                                                    sym=dbcsr_type_no_symmetry, data_type=dbcsr_type_real_8)
     541             : 
     542       32937 :             CALL dbcsr_init_p(qs_ot_env%matrix_h_e(i)%matrix)
     543             :             CALL cp_dbcsr_m_by_n_from_row_template(qs_ot_env%matrix_h_e(i)%matrix, template=matrix_s, n=k, &
     544       41498 :                                                    sym=dbcsr_type_no_symmetry, data_type=dbcsr_type_real_8)
     545             :          END DO
     546             :       END IF
     547             : 
     548        8561 :       NULLIFY (qs_ot_env%rot_mat_u, qs_ot_env%rot_mat_x, qs_ot_env%rot_mat_h_e, qs_ot_env%rot_mat_h_x, &
     549        8561 :                qs_ot_env%rot_mat_gx, qs_ot_env%rot_mat_gx_old, qs_ot_env%rot_mat_dx, &
     550        8561 :                qs_ot_env%rot_mat_evals, qs_ot_env%rot_mat_evec, qs_ot_env%rot_mat_dedu, qs_ot_env%rot_mat_chc)
     551             : 
     552        8561 :       IF (qs_ot_env%settings%do_rotation) THEN
     553         232 :          CALL dbcsr_init_p(qs_ot_env%rot_mat_u)
     554             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_u, template=matrix_s, m=k, n=k, &
     555         232 :                                             sym=dbcsr_type_no_symmetry)
     556             : 
     557         232 :          CALL dbcsr_init_p(qs_ot_env%rot_mat_x)
     558             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_x, template=matrix_s, m=k, n=k, &
     559         232 :                                             sym=dbcsr_type_no_symmetry)
     560             : 
     561         232 :          CALL dbcsr_init_p(qs_ot_env%rot_mat_dedu)
     562             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_dedu, template=matrix_s, m=k, n=k, &
     563         232 :                                             sym=dbcsr_type_no_symmetry)
     564             : 
     565         232 :          CALL dbcsr_init_p(qs_ot_env%rot_mat_chc)
     566             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_chc, template=matrix_s, m=k, n=k, &
     567         232 :                                             sym=dbcsr_type_no_symmetry)
     568             : 
     569         232 :          IF (qs_ot_env%settings%ot_method .EQ. "DIIS") THEN
     570          30 :             CALL dbcsr_allocate_matrix_set(qs_ot_env%rot_mat_h_e, m_diis)
     571          30 :             CALL dbcsr_allocate_matrix_set(qs_ot_env%rot_mat_h_x, m_diis)
     572         240 :             DO i = 1, m_diis
     573         210 :                CALL dbcsr_init_p(qs_ot_env%rot_mat_h_e(i)%matrix)
     574             :                CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_h_e(i)%matrix, template=matrix_s, m=k, n=k, &
     575         210 :                                                   sym=dbcsr_type_no_symmetry)
     576             : 
     577         210 :                CALL dbcsr_init_p(qs_ot_env%rot_mat_h_x(i)%matrix)
     578             :                CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_h_x(i)%matrix, template=matrix_s, m=k, n=k, &
     579         240 :                                                   sym=dbcsr_type_no_symmetry)
     580             :             END DO
     581             :          END IF
     582             : 
     583         692 :          ALLOCATE (qs_ot_env%rot_mat_evals(k))
     584             : 
     585         232 :          CALL dbcsr_init_p(qs_ot_env%rot_mat_evec)
     586             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_evec, template=matrix_s, m=k, n=k, &
     587         232 :                                             sym=dbcsr_type_no_symmetry, data_type=dbcsr_type_complex_default)
     588             : 
     589         232 :          CALL dbcsr_init_p(qs_ot_env%rot_mat_gx)
     590             :          CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_gx, template=matrix_s, m=k, n=k, &
     591         232 :                                             sym=dbcsr_type_no_symmetry)
     592             : 
     593         232 :          IF (qs_ot_env%use_gx_old) THEN
     594         202 :             CALL dbcsr_init_p(qs_ot_env%rot_mat_gx_old)
     595             :             CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_gx_old, template=matrix_s, m=k, n=k, &
     596         202 :                                                sym=dbcsr_type_no_symmetry)
     597             :          END IF
     598             : 
     599         232 :          IF (qs_ot_env%use_dx) THEN
     600         202 :             CALL dbcsr_init_p(qs_ot_env%rot_mat_dx)
     601             :             CALL cp_dbcsr_m_by_n_from_template(qs_ot_env%rot_mat_dx, template=matrix_s, m=k, n=k, &
     602         202 :                                                sym=dbcsr_type_no_symmetry)
     603             :          END IF
     604             : 
     605             :       END IF
     606             : 
     607        8561 :       IF (qs_ot_env%settings%do_ener) THEN
     608           0 :          ncoef = k
     609           0 :          ALLOCATE (qs_ot_env%ener_x(ncoef))
     610             : 
     611           0 :          IF (qs_ot_env%settings%ot_method .EQ. "DIIS") THEN
     612           0 :             ALLOCATE (qs_ot_env%ener_h_e(m_diis, ncoef))
     613           0 :             ALLOCATE (qs_ot_env%ener_h_x(m_diis, ncoef))
     614             :          END IF
     615             : 
     616           0 :          ALLOCATE (qs_ot_env%ener_gx(ncoef))
     617             : 
     618           0 :          IF (qs_ot_env%use_gx_old) THEN
     619           0 :             ALLOCATE (qs_ot_env%ener_gx_old(ncoef))
     620             :          END IF
     621             : 
     622           0 :          IF (qs_ot_env%use_dx) THEN
     623           0 :             ALLOCATE (qs_ot_env%ener_dx(ncoef))
     624           0 :             qs_ot_env%ener_dx = 0.0_dp
     625             :          END IF
     626             :       END IF
     627             : 
     628        8561 :    END SUBROUTINE qs_ot_allocate
     629             : 
     630             :    ! deallocates data
     631             : ! **************************************************************************************************
     632             : !> \brief ...
     633             : !> \param qs_ot_env ...
     634             : ! **************************************************************************************************
     635        8561 :    SUBROUTINE qs_ot_destroy(qs_ot_env)
     636             :       TYPE(qs_ot_type)                                   :: qs_ot_env
     637             : 
     638        8561 :       CALL mp_para_env_release(qs_ot_env%para_env)
     639        8561 :       CALL cp_blacs_env_release(qs_ot_env%blacs_env)
     640             : 
     641        8561 :       DEALLOCATE (qs_ot_env%evals)
     642        8561 :       DEALLOCATE (qs_ot_env%dum)
     643             : 
     644        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_os)) CALL dbcsr_release_p(qs_ot_env%matrix_os)
     645        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_p)) CALL dbcsr_release_p(qs_ot_env%matrix_p)
     646        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_cosp)) CALL dbcsr_release_p(qs_ot_env%matrix_cosp)
     647        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_sinp)) CALL dbcsr_release_p(qs_ot_env%matrix_sinp)
     648        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_r)) CALL dbcsr_release_p(qs_ot_env%matrix_r)
     649        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_cosp_b)) CALL dbcsr_release_p(qs_ot_env%matrix_cosp_b)
     650        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_sinp_b)) CALL dbcsr_release_p(qs_ot_env%matrix_sinp_b)
     651        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_buf1)) CALL dbcsr_release_p(qs_ot_env%matrix_buf1)
     652        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_buf2)) CALL dbcsr_release_p(qs_ot_env%matrix_buf2)
     653        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_buf3)) CALL dbcsr_release_p(qs_ot_env%matrix_buf3)
     654        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_buf4)) CALL dbcsr_release_p(qs_ot_env%matrix_buf4)
     655        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_buf1_ortho)) CALL dbcsr_release_p(qs_ot_env%matrix_buf1_ortho)
     656        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_buf2_ortho)) CALL dbcsr_release_p(qs_ot_env%matrix_buf2_ortho)
     657        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_c0)) CALL dbcsr_release_p(qs_ot_env%matrix_c0)
     658        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_sc0)) CALL dbcsr_release_p(qs_ot_env%matrix_sc0)
     659        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_psc0)) CALL dbcsr_release_p(qs_ot_env%matrix_psc0)
     660        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_x)) CALL dbcsr_release_p(qs_ot_env%matrix_x)
     661        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_sx)) CALL dbcsr_release_p(qs_ot_env%matrix_sx)
     662        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_gx)) CALL dbcsr_release_p(qs_ot_env%matrix_gx)
     663        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_dx)) CALL dbcsr_release_p(qs_ot_env%matrix_dx)
     664        8561 :       IF (ASSOCIATED(qs_ot_env%matrix_gx_old)) CALL dbcsr_release_p(qs_ot_env%matrix_gx_old)
     665        8561 :       IF (ASSOCIATED(qs_ot_env%buf1_k_k_nosym)) CALL dbcsr_release_p(qs_ot_env%buf1_k_k_nosym)
     666        8561 :       IF (ASSOCIATED(qs_ot_env%buf2_k_k_nosym)) CALL dbcsr_release_p(qs_ot_env%buf2_k_k_nosym)
     667        8561 :       IF (ASSOCIATED(qs_ot_env%buf3_k_k_nosym)) CALL dbcsr_release_p(qs_ot_env%buf3_k_k_nosym)
     668        8561 :       IF (ASSOCIATED(qs_ot_env%buf4_k_k_nosym)) CALL dbcsr_release_p(qs_ot_env%buf4_k_k_nosym)
     669        8561 :       IF (ASSOCIATED(qs_ot_env%p_k_k_sym)) CALL dbcsr_release_p(qs_ot_env%p_k_k_sym)
     670        8561 :       IF (ASSOCIATED(qs_ot_env%buf1_k_k_sym)) CALL dbcsr_release_p(qs_ot_env%buf1_k_k_sym)
     671        8561 :       IF (ASSOCIATED(qs_ot_env%buf2_k_k_sym)) CALL dbcsr_release_p(qs_ot_env%buf2_k_k_sym)
     672        8561 :       IF (ASSOCIATED(qs_ot_env%buf3_k_k_sym)) CALL dbcsr_release_p(qs_ot_env%buf3_k_k_sym)
     673        8561 :       IF (ASSOCIATED(qs_ot_env%buf4_k_k_sym)) CALL dbcsr_release_p(qs_ot_env%buf4_k_k_sym)
     674        8561 :       IF (ASSOCIATED(qs_ot_env%buf1_n_k)) CALL dbcsr_release_p(qs_ot_env%buf1_n_k)
     675        8561 :       IF (ASSOCIATED(qs_ot_env%buf1_n_k_dp)) CALL dbcsr_release_p(qs_ot_env%buf1_n_k_dp)
     676             : 
     677        8561 :       IF (qs_ot_env%settings%ot_method .EQ. "DIIS" .OR. &
     678             :           qs_ot_env%settings%ot_method .EQ. "BROY") THEN
     679        4621 :          CALL dbcsr_deallocate_matrix_set(qs_ot_env%matrix_h_x)
     680        4621 :          CALL dbcsr_deallocate_matrix_set(qs_ot_env%matrix_h_e)
     681        4621 :          DEALLOCATE (qs_ot_env%ls_diis)
     682        4621 :          DEALLOCATE (qs_ot_env%lss_diis)
     683        4621 :          DEALLOCATE (qs_ot_env%c_diis)
     684        4621 :          DEALLOCATE (qs_ot_env%c_broy)
     685        4621 :          DEALLOCATE (qs_ot_env%energy_h)
     686        4621 :          DEALLOCATE (qs_ot_env%ipivot)
     687             :       END IF
     688             : 
     689        8561 :       IF (qs_ot_env%settings%do_rotation) THEN
     690             : 
     691         232 :          IF (ASSOCIATED(qs_ot_env%rot_mat_u)) CALL dbcsr_release_p(qs_ot_env%rot_mat_u)
     692         232 :          IF (ASSOCIATED(qs_ot_env%rot_mat_x)) CALL dbcsr_release_p(qs_ot_env%rot_mat_x)
     693         232 :          IF (ASSOCIATED(qs_ot_env%rot_mat_dedu)) CALL dbcsr_release_p(qs_ot_env%rot_mat_dedu)
     694         232 :          IF (ASSOCIATED(qs_ot_env%rot_mat_chc)) CALL dbcsr_release_p(qs_ot_env%rot_mat_chc)
     695             : 
     696         232 :          IF (qs_ot_env%settings%ot_method .EQ. "DIIS") THEN
     697          30 :             CALL dbcsr_deallocate_matrix_set(qs_ot_env%rot_mat_h_x)
     698          30 :             CALL dbcsr_deallocate_matrix_set(qs_ot_env%rot_mat_h_e)
     699             :          END IF
     700             : 
     701         232 :          DEALLOCATE (qs_ot_env%rot_mat_evals)
     702             : 
     703         232 :          IF (ASSOCIATED(qs_ot_env%rot_mat_evec)) CALL dbcsr_release_p(qs_ot_env%rot_mat_evec)
     704         232 :          IF (ASSOCIATED(qs_ot_env%rot_mat_gx)) CALL dbcsr_release_p(qs_ot_env%rot_mat_gx)
     705         232 :          IF (ASSOCIATED(qs_ot_env%rot_mat_gx_old)) CALL dbcsr_release_p(qs_ot_env%rot_mat_gx_old)
     706         232 :          IF (ASSOCIATED(qs_ot_env%rot_mat_dx)) CALL dbcsr_release_p(qs_ot_env%rot_mat_dx)
     707             :       END IF
     708             : 
     709        8561 :       IF (qs_ot_env%settings%do_ener) THEN
     710           0 :          DEALLOCATE (qs_ot_env%ener_x)
     711           0 :          DEALLOCATE (qs_ot_env%ener_gx)
     712           0 :          IF (qs_ot_env%settings%ot_method .EQ. "DIIS") THEN
     713           0 :             DEALLOCATE (qs_ot_env%ener_h_x)
     714           0 :             DEALLOCATE (qs_ot_env%ener_h_e)
     715             :          END IF
     716           0 :          IF (qs_ot_env%use_dx) THEN
     717           0 :             DEALLOCATE (qs_ot_env%ener_dx)
     718             :          END IF
     719           0 :          IF (qs_ot_env%use_gx_old) THEN
     720           0 :             DEALLOCATE (qs_ot_env%ener_gx_old)
     721             :          END IF
     722             :       END IF
     723             : 
     724        8561 :    END SUBROUTINE qs_ot_destroy
     725             : 
     726             : ! **************************************************************************************************
     727             : !> \brief ...
     728             : !> \param settings ...
     729             : !> \param ot_section ...
     730             : !> \param output_unit ...
     731             : ! **************************************************************************************************
     732       33270 :    SUBROUTINE ot_readwrite_input(settings, ot_section, output_unit)
     733             :       TYPE(qs_ot_settings_type)                          :: settings
     734             :       TYPE(section_vals_type), POINTER                   :: ot_section
     735             :       INTEGER, INTENT(IN)                                :: output_unit
     736             : 
     737             :       CHARACTER(len=*), PARAMETER :: routineN = 'ot_readwrite_input'
     738             : 
     739             :       INTEGER                                            :: handle, ls_method, ot_algorithm, &
     740             :                                                             ot_method, ot_ortho_irac
     741             : 
     742        6654 :       CALL timeset(routineN, handle)
     743             : 
     744             :       ! choose algorithm
     745        6654 :       CALL section_vals_val_get(ot_section, "ALGORITHM", i_val=ot_algorithm)
     746        6084 :       SELECT CASE (ot_algorithm)
     747             :       CASE (ot_algo_taylor_or_diag)
     748        6084 :          settings%ot_algorithm = "TOD"
     749             :       CASE (ot_algo_irac)
     750         570 :          CALL cite_reference(Weber2008)
     751         570 :          settings%ot_algorithm = "REF"
     752             :       CASE DEFAULT
     753        6654 :          CPABORT("Value unknown")
     754             :       END SELECT
     755             : 
     756             :       ! irac input
     757        6654 :       CALL section_vals_val_get(ot_section, "IRAC_DEGREE", i_val=settings%irac_degree)
     758        6654 :       IF (settings%irac_degree < 2 .OR. settings%irac_degree > 4) THEN
     759           0 :          CPABORT("READ OT IRAC_DEGREE: Value unknown")
     760             :       END IF
     761        6654 :       CALL section_vals_val_get(ot_section, "MAX_IRAC", i_val=settings%max_irac)
     762        6654 :       IF (settings%max_irac < 1) THEN
     763           0 :          CPABORT("READ OT MAX_IRAC: VALUE MUST BE GREATER THAN ZERO")
     764             :       END IF
     765        6654 :       CALL section_vals_val_get(ot_section, "EPS_IRAC_FILTER_MATRIX", r_val=settings%eps_irac_filter_matrix)
     766        6654 :       CALL section_vals_val_get(ot_section, "EPS_IRAC", r_val=settings%eps_irac)
     767        6654 :       IF (settings%eps_irac < 0.0_dp) THEN
     768           0 :          CPABORT("READ OT EPS_IRAC: VALUE MUST BE GREATER THAN ZERO")
     769             :       END IF
     770        6654 :       CALL section_vals_val_get(ot_section, "EPS_IRAC_QUICK_EXIT", r_val=settings%eps_irac_quick_exit)
     771        6654 :       IF (settings%eps_irac_quick_exit < 0.0_dp) THEN
     772           0 :          CPABORT("READ OT EPS_IRAC_QUICK_EXIT: VALUE MUST BE GREATER THAN ZERO")
     773             :       END IF
     774             : 
     775        6654 :       CALL section_vals_val_get(ot_section, "EPS_IRAC_SWITCH", r_val=settings%eps_irac_switch)
     776        6654 :       IF (settings%eps_irac_switch < 0.0_dp) THEN
     777           0 :          CPABORT("READ OT EPS_IRAC_SWITCH: VALUE MUST BE GREATER THAN ZERO")
     778             :       END IF
     779             : 
     780        6654 :       CALL section_vals_val_get(ot_section, "ORTHO_IRAC", i_val=ot_ortho_irac)
     781        6570 :       SELECT CASE (ot_ortho_irac)
     782             :       CASE (ot_chol_irac)
     783        6570 :          settings%ortho_irac = "CHOL"
     784             :       CASE (ot_poly_irac)
     785          34 :          settings%ortho_irac = "POLY"
     786             :       CASE (ot_lwdn_irac)
     787          50 :          settings%ortho_irac = "LWDN"
     788             :       CASE DEFAULT
     789        6654 :          CPABORT("READ OT ORTHO_IRAC: Value unknown")
     790             :       END SELECT
     791             : 
     792        6654 :       CALL section_vals_val_get(ot_section, "ON_THE_FLY_LOC", l_val=settings%on_the_fly_loc)
     793             : 
     794        6654 :       CALL section_vals_val_get(ot_section, "MINIMIZER", i_val=ot_method)
     795             :       ! compatibility
     796          10 :       SELECT CASE (ot_method)
     797             :       CASE (ot_mini_sd)
     798          10 :          settings%ot_method = "SD"
     799             :       CASE (ot_mini_cg)
     800        2301 :          settings%ot_method = "CG"
     801             :       CASE (ot_mini_diis)
     802        4329 :          settings%ot_method = "DIIS"
     803        4329 :          CALL section_vals_val_get(ot_section, "N_HISTORY_VEC", i_val=settings%diis_m)
     804             :       CASE (ot_mini_broyden)
     805          14 :          CALL section_vals_val_get(ot_section, "N_HISTORY_VEC", i_val=settings%diis_m)
     806          14 :          CALL section_vals_val_get(ot_section, "BROYDEN_BETA", r_val=settings%broyden_beta)
     807          14 :          CALL section_vals_val_get(ot_section, "BROYDEN_GAMMA", r_val=settings%broyden_gamma)
     808          14 :          CALL section_vals_val_get(ot_section, "BROYDEN_SIGMA", r_val=settings%broyden_sigma)
     809          14 :          CALL section_vals_val_get(ot_section, "BROYDEN_ETA", r_val=settings%broyden_eta)
     810          14 :          CALL section_vals_val_get(ot_section, "BROYDEN_OMEGA", r_val=settings%broyden_omega)
     811          14 :          CALL section_vals_val_get(ot_section, "BROYDEN_SIGMA_DECREASE", r_val=settings%broyden_sigma_decrease)
     812          14 :          CALL section_vals_val_get(ot_section, "BROYDEN_SIGMA_MIN", r_val=settings%broyden_sigma_min)
     813          14 :          CALL section_vals_val_get(ot_section, "BROYDEN_FORGET_HISTORY", l_val=settings%broyden_forget_history)
     814          14 :          CALL section_vals_val_get(ot_section, "BROYDEN_ADAPTIVE_SIGMA", l_val=settings%broyden_adaptive_sigma)
     815          14 :          CALL section_vals_val_get(ot_section, "BROYDEN_ENABLE_FLIP", l_val=settings%broyden_enable_flip)
     816          14 :          settings%ot_method = "BROY"
     817             :       CASE DEFAULT
     818        6654 :          CPABORT("READ OTSCF MINIMIZER: Value unknown")
     819             :       END SELECT
     820        6654 :       CALL section_vals_val_get(ot_section, "SAFER_DIIS", l_val=settings%safer_diis)
     821        6654 :       CALL section_vals_val_get(ot_section, "LINESEARCH", i_val=ls_method)
     822           2 :       SELECT CASE (ls_method)
     823             :       CASE (ls_none)
     824           2 :          settings%line_search_method = "NONE"
     825             :       CASE (ls_2pnt)
     826        6576 :          settings%line_search_method = "2PNT"
     827             :       CASE (ls_3pnt)
     828          64 :          settings%line_search_method = "3PNT"
     829             :       CASE (ls_gold)
     830          12 :          settings%line_search_method = "GOLD"
     831          12 :          CALL section_vals_val_get(ot_section, "GOLD_TARGET", r_val=settings%gold_target)
     832             :       CASE DEFAULT
     833        6654 :          CPABORT("READ OTSCF LS: Value unknown")
     834             :       END SELECT
     835             : 
     836        6654 :       CALL section_vals_val_get(ot_section, "PRECOND_SOLVER", i_val=settings%precond_solver_type)
     837       13286 :       SELECT CASE (settings%precond_solver_type)
     838             :       CASE (ot_precond_solver_default)
     839        6632 :          settings%precond_solver_name = "DEFAULT"
     840             :       CASE (ot_precond_solver_inv_chol)
     841          16 :          settings%precond_solver_name = "INVERSE_CHOLESKY"
     842             :       CASE (ot_precond_solver_direct)
     843           0 :          settings%precond_solver_name = "DIRECT"
     844             :       CASE (ot_precond_solver_update)
     845           6 :          settings%precond_solver_name = "INVERSE_UPDATE"
     846             :       CASE DEFAULT
     847        6654 :          CPABORT("READ OTSCF SOLVER: Value unknown")
     848             :       END SELECT
     849             : 
     850             :       !If these values are negative we will set them "optimal" for a given precondtioner below
     851        6654 :       CALL section_vals_val_get(ot_section, "STEPSIZE", r_val=settings%ds_min)
     852        6654 :       CALL section_vals_val_get(ot_section, "ENERGY_GAP", r_val=settings%energy_gap)
     853             : 
     854        6654 :       CALL section_vals_val_get(ot_section, "PRECONDITIONER", i_val=settings%preconditioner_type)
     855        7442 :       SELECT CASE (settings%preconditioner_type)
     856             :       CASE (ot_precond_none)
     857         788 :          settings%preconditioner_name = "NONE"
     858         788 :          IF (settings%ds_min < 0.0_dp) settings%ds_min = 0.15_dp
     859         788 :          IF (settings%energy_gap < 0.0_dp) settings%energy_gap = 0.2_dp
     860             :       CASE (ot_precond_full_single)
     861          40 :          settings%preconditioner_name = "FULL_SINGLE"
     862          40 :          IF (settings%ds_min < 0.0_dp) settings%ds_min = 0.15_dp
     863          40 :          IF (settings%energy_gap < 0.0_dp) settings%energy_gap = 0.2_dp
     864             :       CASE (ot_precond_full_single_inverse)
     865        2737 :          settings%preconditioner_name = "FULL_SINGLE_INVERSE"
     866        2737 :          IF (settings%ds_min < 0.0_dp) settings%ds_min = 0.08_dp
     867        2737 :          IF (settings%energy_gap < 0.0_dp) settings%energy_gap = 0.08_dp
     868             :       CASE (ot_precond_full_all)
     869        1906 :          settings%preconditioner_name = "FULL_ALL"
     870        1906 :          IF (settings%ds_min < 0.0_dp) settings%ds_min = 0.15_dp
     871        1906 :          IF (settings%energy_gap < 0.0_dp) settings%energy_gap = 0.08_dp
     872             :       CASE (ot_precond_full_kinetic)
     873        1147 :          settings%preconditioner_name = "FULL_KINETIC"
     874        1147 :          IF (settings%ds_min < 0.0_dp) settings%ds_min = 0.15_dp
     875        1147 :          IF (settings%energy_gap < 0.0_dp) settings%energy_gap = 0.2_dp
     876             :       CASE (ot_precond_s_inverse)
     877          36 :          settings%preconditioner_name = "FULL_S_INVERSE"
     878          36 :          IF (settings%ds_min < 0.0_dp) settings%ds_min = 0.15_dp
     879          36 :          IF (settings%energy_gap < 0.0_dp) settings%energy_gap = 0.2_dp
     880             :       CASE DEFAULT
     881        6654 :          CPABORT("READ OTSCF PRECONDITIONER: Value unknown")
     882             :       END SELECT
     883        6654 :       CALL section_vals_val_get(ot_section, "CHOLESKY", i_val=settings%cholesky_type)
     884        6654 :       CALL section_vals_val_get(ot_section, "EPS_TAYLOR", r_val=settings%eps_taylor)
     885        6654 :       CALL section_vals_val_get(ot_section, "MAX_TAYLOR", i_val=settings%max_taylor)
     886        6654 :       CALL section_vals_val_get(ot_section, "ROTATION", l_val=settings%do_rotation)
     887        6654 :       CALL section_vals_val_get(ot_section, "ENERGIES", l_val=settings%do_ener)
     888             :       CALL section_vals_val_get(ot_section, "OCCUPATION_PRECONDITIONER", &
     889        6654 :                                 l_val=settings%occupation_preconditioner)
     890        6654 :       CALL section_vals_val_get(ot_section, "NONDIAG_ENERGY", l_val=settings%add_nondiag_energy)
     891             :       CALL section_vals_val_get(ot_section, "NONDIAG_ENERGY_STRENGTH", &
     892        6654 :                                 r_val=settings%nondiag_energy_strength)
     893             :       ! not yet fully implemented
     894        6654 :       CPASSERT(.NOT. settings%do_ener)
     895             : 
     896             :       ! write OT output
     897             : 
     898        6654 :       IF (output_unit > 0) THEN
     899        3444 :          WRITE (output_unit, '(/,A)') "  ----------------------------------- OT ---------------------------------------"
     900        3444 :          IF (settings%do_rotation) THEN
     901         110 :             WRITE (output_unit, '(A)') "  Allowing for rotations "
     902             :          END IF
     903        3444 :          IF (settings%do_ener) THEN
     904           0 :             WRITE (output_unit, '(A,L2)') "  Optimizing orbital energies "
     905             :          END IF
     906           5 :          SELECT CASE (settings%OT_METHOD)
     907             :          CASE ("SD")
     908           5 :             WRITE (output_unit, '(A)') "  Minimizer      : SD                  : steepest descent"
     909             :          CASE ("CG")
     910        1180 :             WRITE (output_unit, '(A)') "  Minimizer      : CG                  : conjugate gradient"
     911             :          CASE ("DIIS")
     912        2252 :             WRITE (output_unit, '(A)') "  Minimizer      : DIIS                : direct inversion"
     913        2252 :             WRITE (output_unit, '(A)') "                                         in the iterative subspace"
     914        2252 :             WRITE (output_unit, '(A,I3,A)') "                                         using ", settings%diis_m, " DIIS vectors"
     915        2252 :             IF (settings%safer_diis) THEN
     916        2252 :                WRITE (output_unit, '(A,I3,A)') "                                         safer DIIS on"
     917             :             ELSE
     918           0 :                WRITE (output_unit, '(A,I3,A)') "                                         safer DIIS off"
     919             :             END IF
     920             :          CASE ("BROY")
     921           7 :             WRITE (output_unit, '(A)') "  Minimizer      : BROYDEN             : Broyden "
     922           7 :             WRITE (output_unit, '(A,F16.8)') "                   BETA                : ", settings%broyden_beta
     923           7 :             WRITE (output_unit, '(A,F16.8)') "                   GAMMA               : ", settings%broyden_gamma
     924           7 :             WRITE (output_unit, '(A,F16.8)') "                   SIGMA               : ", settings%broyden_sigma
     925           7 :             WRITE (output_unit, '(A,I3,A)') "                            using      : - ", &
     926          14 :                settings%diis_m, " BROYDEN vectors"
     927             :          CASE DEFAULT
     928        3444 :             WRITE (output_unit, '(3A)') "  Minimizer      :      ", settings%OT_METHOD, " : UNKNOWN"
     929             :          END SELECT
     930          20 :          SELECT CASE (settings%preconditioner_name)
     931             :          CASE ("FULL_SINGLE")
     932          20 :             WRITE (output_unit, '(A)') "  Preconditioner : FULL_SINGLE         : diagonalization based"
     933             :          CASE ("FULL_SINGLE_INVERSE")
     934        1436 :             WRITE (output_unit, '(A,/,A)') "  Preconditioner : FULL_SINGLE_INVERSE : inversion of ", &
     935        2872 :                "                                         H + eS - 2*(Sc)(c^T*H*c+const)(Sc)^T"
     936             :          CASE ("FULL_ALL")
     937         983 :             WRITE (output_unit, '(A)') "  Preconditioner : FULL_ALL            : diagonalization, state selective"
     938             :          CASE ("FULL_KINETIC")
     939         593 :             WRITE (output_unit, '(A)') "  Preconditioner : FULL_KINETIC        : inversion of T + eS"
     940             :          CASE ("FULL_S_INVERSE")
     941          18 :             WRITE (output_unit, '(A)') "  Preconditioner : FULL_S_INVERSE      : cholesky inversion of S"
     942             :          CASE ("SPARSE_DIAG")
     943             :             WRITE (output_unit, '(A)') &
     944           0 :                "  Preconditioner : SPARSE_DIAG    : diagonal atomic block diagonalization"
     945             :          CASE ("SPARSE_KINETIC")
     946           0 :             WRITE (output_unit, '(A)') "  Preconditioner : SPARSE_KINETIC      : sparse linear solver for T + eS"
     947             :          CASE ("NONE")
     948         394 :             WRITE (output_unit, '(A)') "  Preconditioner : NONE"
     949             :          CASE DEFAULT
     950        3444 :             WRITE (output_unit, '(3A)') "  Preconditioner : ", settings%preconditioner_name, " : UNKNOWN"
     951             :          END SELECT
     952             : 
     953        3444 :          WRITE (output_unit, '(A)') "  Precond_solver : "//TRIM(settings%precond_solver_name)
     954             : 
     955        3444 :          IF (settings%OT_METHOD .EQ. "SD" .OR. settings%OT_METHOD .EQ. "CG") THEN
     956        1148 :             SELECT CASE (settings%line_search_method)
     957             :             CASE ("2PNT")
     958        1148 :                WRITE (output_unit, '(A)') "  Line search    : 2PNT                : 2 energies, one gradient"
     959             :             CASE ("3PNT")
     960          30 :                WRITE (output_unit, '(A)') "  Line search    : 3PNT                : 3 energies"
     961             :             CASE ("GOLD")
     962           6 :                WRITE (output_unit, '(A)') "  Line search    : GOLD                : bracketing and golden section search"
     963           6 :                WRITE (output_unit, '(A,F14.8)') "                   target rel accuracy : ", settings%gold_target
     964             :             CASE ("NONE")
     965           1 :                WRITE (output_unit, '(A)') "  Line search    : NONE"
     966             :             CASE DEFAULT
     967        3444 :                WRITE (output_unit, '(3A)') "  Line search : ", settings%line_search_method, " : UNKNOWN"
     968             :             END SELECT
     969             :          END IF
     970        3444 :          WRITE (output_unit, '(A,F14.8,T49,A,F14.8)') "  stepsize       :", settings%ds_min, &
     971        6888 :             "  energy_gap     :", settings%energy_gap
     972        3444 :          IF (settings%ot_algorithm .EQ. 'TOD') THEN
     973        3141 :             WRITE (output_unit, '(A,E14.5,T49,A,I14)') "  eps_taylor     :", settings%eps_taylor, &
     974        6282 :                "  max_taylor     :", settings%max_taylor
     975             :          END IF
     976        3444 :          IF (settings%ot_algorithm .EQ. 'REF') THEN
     977         303 :             WRITE (output_unit, '(A,1X,A,T49,A,I14)') "  ortho_irac     :", settings%ortho_irac, &
     978         606 :                "  irac_degree    :", settings%irac_degree
     979         303 :             WRITE (output_unit, '(A,I14,T49,A,E14.5)') "  max_irac       :", settings%max_irac, &
     980         606 :                "  eps_irac       :", settings%eps_irac
     981         303 :             WRITE (output_unit, '(A,E14.5,T49,A,E10.3)') "  eps_irac_switch:", settings%eps_irac_switch, &
     982         606 :                "  eps_irac_quick_exit:", settings%eps_irac_quick_exit
     983         303 :             WRITE (output_unit, '(A,L2)') "  on_the_fly_loc :", settings%on_the_fly_loc
     984             :          END IF
     985        3444 :          WRITE (output_unit, '(A)') "  ----------------------------------- OT ---------------------------------------"
     986             :          WRITE (UNIT=output_unit, &
     987             :                 FMT="(/,T3,A,T12,A,T31,A,T39,A,T59,A,T75,A,/,T3,A)") &
     988        3444 :             "Step", "Update method", "Time", "Convergence", "Total energy", "Change", &
     989        6888 :             REPEAT("-", 78)
     990             :       END IF
     991             : 
     992        6654 :       CALL timestop(handle)
     993             : 
     994        6654 :    END SUBROUTINE ot_readwrite_input
     995             : ! **************************************************************************************************
     996             : 
     997           0 : END MODULE qs_ot_types

Generated by: LCOV version 1.15