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 Calculate MAO's and analyze wavefunctions
      10             : !> \par History
      11             : !>      03.2016 created [JGH]
      12             : !>      12.2016 split into four modules [JGH]
      13             : !> \author JGH
      14             : ! **************************************************************************************************
      15             : MODULE mao_optimizer
      16             :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
      17             :                                               dbcsr_deallocate_matrix_set
      18             :    USE dbcsr_api,                       ONLY: &
      19             :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_get_info, dbcsr_norm, &
      20             :         dbcsr_norm_maxabsnorm, dbcsr_p_type, dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_type, &
      21             :         dbcsr_type_symmetric
      22             :    USE kinds,                           ONLY: dp
      23             :    USE mao_methods,                     ONLY: mao_function,&
      24             :                                               mao_function_gradient,&
      25             :                                               mao_initialization,&
      26             :                                               mao_orthogonalization,&
      27             :                                               mao_scalar_product
      28             : #include "./base/base_uses.f90"
      29             : 
      30             :    IMPLICIT NONE
      31             :    PRIVATE
      32             : 
      33             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_optimizer'
      34             : 
      35             :    PUBLIC ::  mao_optimize
      36             : 
      37             : ! **************************************************************************************************
      38             : 
      39             : CONTAINS
      40             : 
      41             : ! **************************************************************************************************
      42             : !> \brief ...
      43             : !> \param mao_coef ...
      44             : !> \param matrix_q ...
      45             : !> \param matrix_smm ...
      46             : !> \param electra ...
      47             : !> \param max_iter ...
      48             : !> \param eps_grad ...
      49             : !> \param eps1 ...
      50             : !> \param iolevel ...
      51             : !> \param iw ...
      52             : ! **************************************************************************************************
      53          14 :    SUBROUTINE mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, eps1, &
      54             :                            iolevel, iw)
      55             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mao_coef, matrix_q, matrix_smm
      56             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: electra
      57             :       INTEGER, INTENT(IN)                                :: max_iter
      58             :       REAL(KIND=dp), INTENT(IN)                          :: eps_grad, eps1
      59             :       INTEGER, INTENT(IN)                                :: iolevel, iw
      60             : 
      61             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'mao_optimize'
      62             : 
      63             :       INTEGER                                            :: handle, i, ispin, iter, nspin
      64          14 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_sizes, mao_sizes_a, mao_sizes_b
      65             :       LOGICAL                                            :: spin_warning
      66             :       REAL(KIND=dp)                                      :: a1, a2, alpha, an, beta, eps_fun, fa1, &
      67             :                                                             fa2, fnnew, fnold, fval, grad_norm
      68             :       TYPE(dbcsr_distribution_type)                      :: dbcsr_dist
      69          14 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mao_grad
      70             :       TYPE(dbcsr_type)                                   :: amat, binv, cgmat
      71             : 
      72          14 :       CALL timeset(routineN, handle)
      73             : 
      74          14 :       eps_fun = 10._dp*eps_grad
      75             : 
      76          14 :       nspin = SIZE(mao_coef, 1)
      77             : 
      78          14 :       IF (iw > 0) THEN
      79           7 :          WRITE (iw, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
      80           7 :          WRITE (UNIT=iw, FMT="(T32,A)") "MAO BASIS GENERATION"
      81           7 :          WRITE (iw, '(T2,A)') '!-----------------------------------------------------------------------------!'
      82             :       END IF
      83             : 
      84             :       ! initialize MAO coeficients from diagonal blocks of the Q matrix
      85          30 :       DO ispin = 1, nspin
      86             :          CALL mao_initialization(mao_coef(ispin)%matrix, &
      87          30 :                                  matrix_q(ispin)%matrix, matrix_smm(1)%matrix, eps1, iolevel, iw)
      88             :       END DO
      89             :       ! check for mao sizes
      90          14 :       spin_warning = .FALSE.
      91          14 :       CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=mao_sizes_a, distribution=dbcsr_dist)
      92          14 :       IF (nspin == 2) THEN
      93           2 :          CALL dbcsr_get_info(mao_coef(2)%matrix, col_blk_size=mao_sizes_b, distribution=dbcsr_dist)
      94           8 :          DO i = 1, SIZE(mao_sizes_a)
      95           8 :             IF (mao_sizes_a(i) /= mao_sizes_b(i)) spin_warning = .TRUE.
      96             :          END DO
      97             :       END IF
      98          14 :       IF (iw > 0 .AND. spin_warning) THEN
      99           0 :          WRITE (UNIT=iw, FMT="(/,A)") ">>>> WARNING: Different number of MAOs for different spins"
     100             :       END IF
     101          14 :       IF (iw > 0) THEN
     102          15 :          DO ispin = 1, nspin
     103          15 :             IF (ispin == 1) THEN
     104          43 :                WRITE (UNIT=iw, FMT="(/,A,I2,T70,I11)") " MAO's for SPIN ", ispin, SUM(mao_sizes_a)
     105           7 :                WRITE (UNIT=iw, FMT="(20I4)") mao_sizes_a(:)
     106           1 :             ELSE IF (ispin == 2) THEN
     107           4 :                WRITE (UNIT=iw, FMT="(/,A,I2,T70,I11)") " MAO's for SPIN ", ispin, SUM(mao_sizes_b)
     108           1 :                WRITE (UNIT=iw, FMT="(20I4)") mao_sizes_b(:)
     109             :             END IF
     110             :          END DO
     111           7 :          WRITE (UNIT=iw, FMT="(A)")
     112             :       END IF
     113             : 
     114          14 :       IF (max_iter < 1) THEN
     115             :          ! projection only
     116          12 :          DO ispin = 1, nspin
     117           6 :             CALL dbcsr_get_info(mao_coef(ispin)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
     118             :             CALL dbcsr_create(binv, name="Binv", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
     119           6 :                               row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes, nze=0)
     120             :             CALL mao_function(mao_coef(ispin)%matrix, fval, matrix_q(ispin)%matrix, &
     121           6 :                               matrix_smm(1)%matrix, binv, .FALSE.)
     122           6 :             IF (iw > 0) THEN
     123             :                WRITE (UNIT=iw, FMT="(T2,A,T31,A,I2,T56,A,F12.8)") &
     124           3 :                   "MAO Projection", "Spin:", ispin, "Completeness:", fval/electra(ispin)
     125             :             END IF
     126          12 :             CALL dbcsr_release(binv)
     127             :          END DO
     128           6 :          IF (iw > 0) WRITE (UNIT=iw, FMT="(A)")
     129             :       ELSE
     130             :          ! optimize MAOs
     131           8 :          NULLIFY (mao_grad)
     132           8 :          CALL dbcsr_allocate_matrix_set(mao_grad, nspin)
     133          18 :          DO ispin = 1, nspin
     134          10 :             ALLOCATE (mao_grad(ispin)%matrix)
     135          10 :             CALL dbcsr_create(matrix=mao_grad(ispin)%matrix, name="MAO_GRAD", template=mao_coef(1)%matrix)
     136          18 :             CALL dbcsr_reserve_diag_blocks(matrix=mao_grad(ispin)%matrix)
     137             :          END DO
     138           8 :          CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist)
     139             :          CALL dbcsr_create(binv, name="Binv", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
     140           8 :                            row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes, nze=0)
     141           8 :          CALL dbcsr_create(cgmat, template=mao_grad(1)%matrix)
     142           8 :          CALL dbcsr_create(amat, template=mao_coef(1)%matrix)
     143          18 :          DO ispin = 1, nspin
     144          10 :             alpha = 0.25_dp
     145          10 :             beta = 0.0_dp
     146             :             CALL mao_function_gradient(mao_coef(ispin)%matrix, fval, mao_grad(ispin)%matrix, &
     147          10 :                                        matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .FALSE.)
     148          10 :             CALL dbcsr_copy(cgmat, mao_grad(ispin)%matrix)
     149          10 :             CALL dbcsr_norm(mao_grad(ispin)%matrix, dbcsr_norm_maxabsnorm, norm_scalar=grad_norm)
     150          10 :             fnold = mao_scalar_product(mao_grad(ispin)%matrix, mao_grad(ispin)%matrix)
     151          10 :             IF (iw > 0) THEN
     152           5 :                WRITE (UNIT=iw, FMT="(/,T2,A,T73,A,I2)") "MAO OPTIMIZATION", "Spin =", ispin
     153             :                WRITE (UNIT=iw, FMT="(T2,A,T24,A,F11.8,T48,A,F11.8,T69,A,F6.3)") &
     154           5 :                   "Initialization", "fval =", fval/electra(ispin), "grad =", grad_norm, "step =", alpha
     155             :             END IF
     156         174 :             DO iter = 1, max_iter
     157         160 :                IF (grad_norm < eps_grad) EXIT
     158         160 :                IF ((1.0_dp - fval/electra(ispin)) < eps_fun) EXIT
     159         156 :                CALL dbcsr_add(mao_coef(ispin)%matrix, cgmat, 1.0_dp, alpha)
     160         156 :                CALL mao_orthogonalization(mao_coef(ispin)%matrix, matrix_smm(1)%matrix)
     161             :                CALL mao_function_gradient(mao_coef(ispin)%matrix, fval, mao_grad(ispin)%matrix, &
     162         156 :                                           matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .TRUE.)
     163         156 :                CALL dbcsr_norm(mao_grad(ispin)%matrix, dbcsr_norm_maxabsnorm, norm_scalar=grad_norm)
     164         156 :                IF (iw > 0) THEN
     165             :                   WRITE (UNIT=iw, FMT="(T2,A,i8,T24,A,F11.8,T48,A,F11.8,T69,A,F6.3)") &
     166          78 :                      "iter=", iter, "fval =", fval/electra(ispin), "grad =", grad_norm, "step =", alpha
     167             :                END IF
     168         156 :                fnnew = mao_scalar_product(mao_grad(ispin)%matrix, mao_grad(ispin)%matrix)
     169         156 :                beta = fnnew/fnold
     170         156 :                CALL dbcsr_add(cgmat, mao_grad(ispin)%matrix, beta, 1.0_dp)
     171         156 :                fnold = fnnew
     172             :                ! line search, update alpha
     173         156 :                CALL dbcsr_copy(amat, mao_coef(ispin)%matrix)
     174         156 :                CALL dbcsr_add(amat, cgmat, 1.0_dp, 0.5_dp*alpha)
     175         156 :                CALL mao_orthogonalization(amat, matrix_smm(1)%matrix)
     176         156 :                CALL mao_function(amat, fa1, matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .TRUE.)
     177         156 :                CALL dbcsr_copy(amat, mao_coef(ispin)%matrix)
     178         156 :                CALL dbcsr_add(amat, cgmat, 1.0_dp, alpha)
     179         156 :                CALL mao_orthogonalization(amat, matrix_smm(1)%matrix)
     180         156 :                CALL mao_function(amat, fa2, matrix_q(ispin)%matrix, matrix_smm(1)%matrix, binv, .TRUE.)
     181         156 :                a2 = (4._dp*fa1 - fa2 - 3._dp*fval)/alpha
     182         156 :                a1 = (fa2 - fval - a2*alpha)/(alpha*alpha)
     183         156 :                IF (ABS(a1) > 1.e-14_dp) THEN
     184         156 :                   an = -a2/(2._dp*a1)
     185         156 :                   an = MIN(an, 2.0_dp*alpha)
     186             :                ELSE
     187           0 :                   an = 2.0_dp*alpha
     188             :                END IF
     189         166 :                IF (an < 0.05_dp .OR. a1 > 0.0_dp) THEN
     190           0 :                   CALL dbcsr_copy(cgmat, mao_grad(ispin)%matrix)
     191           0 :                   alpha = 0.25_dp
     192             :                ELSE
     193         156 :                   alpha = an
     194             :                END IF
     195             :             END DO
     196             :          END DO
     197           8 :          CALL dbcsr_release(binv)
     198           8 :          CALL dbcsr_release(cgmat)
     199           8 :          CALL dbcsr_release(amat)
     200           8 :          IF (ASSOCIATED(mao_grad)) CALL dbcsr_deallocate_matrix_set(mao_grad)
     201             :       END IF
     202             : 
     203          14 :       CALL timestop(handle)
     204             : 
     205          14 :    END SUBROUTINE mao_optimize
     206             : 
     207             : END MODULE mao_optimizer