Line data    Source code

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