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 Common framework for using eigenvectors of a Fock matrix as PAO basis.
      10             : !> \author Ole Schuett
      11             : ! **************************************************************************************************
      12             : MODULE pao_param_fock
      13             :    USE dbcsr_api,                       ONLY: dbcsr_get_block_p,&
      14             :                                               dbcsr_get_info
      15             :    USE kinds,                           ONLY: dp
      16             :    USE mathlib,                         ONLY: diamat_all
      17             :    USE pao_types,                       ONLY: pao_env_type
      18             : #include "./base/base_uses.f90"
      19             : 
      20             :    IMPLICIT NONE
      21             : 
      22             :    PRIVATE
      23             : 
      24             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_param_fock'
      25             : 
      26             :    PUBLIC :: pao_calc_U_block_fock
      27             : 
      28             : CONTAINS
      29             : 
      30             : ! **************************************************************************************************
      31             : !> \brief Calculate new matrix U and optinally its gradient G
      32             : !> \param pao ...
      33             : !> \param iatom ...
      34             : !> \param V ...
      35             : !> \param U ...
      36             : !> \param penalty ...
      37             : !> \param gap ...
      38             : !> \param evals ...
      39             : !> \param M1 ...
      40             : !> \param G ...
      41             : ! **************************************************************************************************
      42       15520 :    SUBROUTINE pao_calc_U_block_fock(pao, iatom, V, U, penalty, gap, evals, M1, G)
      43             :       TYPE(pao_env_type), POINTER                        :: pao
      44             :       INTEGER, INTENT(IN)                                :: iatom
      45             :       REAL(dp), DIMENSION(:, :), POINTER                 :: V, U
      46             :       REAL(dp), INTENT(INOUT), OPTIONAL                  :: penalty
      47             :       REAL(dp), INTENT(OUT)                              :: gap
      48             :       REAL(dp), DIMENSION(:), INTENT(OUT), OPTIONAL      :: evals
      49             :       REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER       :: M1, G
      50             : 
      51             :       CHARACTER(len=*), PARAMETER :: routineN = 'pao_calc_U_block_fock'
      52             : 
      53             :       INTEGER                                            :: handle, i, j, m, n
      54       15520 :       INTEGER, DIMENSION(:), POINTER                     :: blk_sizes_pao, blk_sizes_pri
      55             :       LOGICAL                                            :: found
      56             :       REAL(dp)                                           :: alpha, beta, denom, diff
      57       15520 :       REAL(dp), DIMENSION(:), POINTER                    :: H_evals
      58       15520 :       REAL(dp), DIMENSION(:, :), POINTER                 :: block_N, D1, D2, H, H0, H_evecs, M2, M3, &
      59       15520 :                                                             M4, M5
      60             : 
      61       15520 :       CALL timeset(routineN, handle)
      62             : 
      63       15520 :       CALL dbcsr_get_block_p(matrix=pao%matrix_H0, row=iatom, col=iatom, block=H0, found=found)
      64       15520 :       CPASSERT(ASSOCIATED(H0))
      65       15520 :       CALL dbcsr_get_block_p(matrix=pao%matrix_N_diag, row=iatom, col=iatom, block=block_N, found=found)
      66       15520 :       CPASSERT(ASSOCIATED(block_N))
      67     1019058 :       IF (MAXVAL(ABS(V - TRANSPOSE(V))) > 1e-14_dp) CPABORT("Expect symmetric matrix")
      68             : 
      69             :       ! figure out basis sizes
      70       15520 :       CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
      71       15520 :       n = blk_sizes_pri(iatom) ! size of primary basis
      72       15520 :       m = blk_sizes_pao(iatom) ! size of pao basis
      73             : 
      74             :       ! calculate H in the orthonormal basis
      75       62080 :       ALLOCATE (H(n, n))
      76    61749582 :       H = MATMUL(MATMUL(block_N, H0 + V), block_N)
      77             : 
      78             :       ! diagonalize H
      79       93120 :       ALLOCATE (H_evals(n), H_evecs(n, n))
      80     2022596 :       H_evecs = H
      81       15520 :       CALL diamat_all(H_evecs, H_evals)
      82             : 
      83             :       ! the eigenvectors of H become the rotation matrix U
      84     2022596 :       U = H_evecs
      85             : 
      86             :       ! copy eigenvectors around the gap from H_evals into evals array
      87       15520 :       IF (PRESENT(evals)) THEN
      88       12786 :          CPASSERT(MOD(SIZE(evals), 2) == 0) ! gap will be exactely in the middle
      89       12786 :          i = SIZE(evals)/2
      90       12786 :          j = MIN(m, i)
      91       53320 :          evals(1 + i - j:i) = H_evals(1 + m - j:m) ! eigenvalues below gap
      92       12786 :          j = MIN(n - m, i)
      93       58783 :          evals(i:i + j) = H_evals(m:m + j) ! eigenvalues above gap
      94             :       END IF
      95             : 
      96             :       ! calculate homo-lumo gap (it's useful for detecting numerical issues)
      97       15520 :       gap = HUGE(dp)
      98       15520 :       IF (m < n) & ! catch special case n==m
      99       14978 :          gap = H_evals(m + 1) - H_evals(m)
     100             : 
     101       15520 :       IF (PRESENT(penalty)) THEN
     102             :          ! penalty terms: occupied and virtual eigenvalues repel each other
     103       15200 :          alpha = pao%penalty_strength
     104       15200 :          beta = pao%penalty_dist
     105       67564 :          DO i = 1, m
     106      222228 :          DO j = m + 1, n
     107      154664 :             diff = H_evals(i) - H_evals(j)
     108      207028 :             penalty = penalty + alpha*EXP(-(diff/beta)**2)
     109             :          END DO
     110             :          END DO
     111             : 
     112             :          ! regularization energy
     113     1007324 :          penalty = penalty + pao%regularization*SUM(V**2)
     114             :       END IF
     115             : 
     116       15520 :       IF (PRESENT(G)) THEN ! TURNING POINT (if calc grad) -------------------------
     117             : 
     118        2471 :          CPASSERT(PRESENT(M1))
     119             : 
     120             :          ! calculate derivatives between eigenvectors of H
     121       32123 :          ALLOCATE (D1(n, n), M2(n, n), M3(n, n), M4(n, n))
     122       19060 :          DO i = 1, n
     123      157033 :          DO j = 1, n
     124             :             ! ignore changes among occupied or virtual eigenvectors
     125             :             ! They will get filtered out by M2*D1 anyways, however this early
     126             :             ! intervention might stabilize numerics in the case of level-crossings.
     127      154562 :             IF (i <= m .EQV. j <= m) THEN
     128       84705 :                D1(i, j) = 0.0_dp
     129             :             ELSE
     130       53268 :                denom = H_evals(i) - H_evals(j)
     131       53268 :                IF (ABS(denom) > 1e-9_dp) THEN ! avoid division by zero
     132       53268 :                   D1(i, j) = 1.0_dp/denom
     133             :                ELSE
     134           0 :                   D1(i, j) = SIGN(1e+9_dp, denom)
     135             :                END IF
     136             :             END IF
     137             :          END DO
     138             :          END DO
     139        2471 :          IF (ASSOCIATED(M1)) THEN
     140     3433480 :             M2 = MATMUL(TRANSPOSE(M1), H_evecs)
     141             :          ELSE
     142           0 :             M2 = 0.0_dp
     143             :          END IF
     144      311595 :          M3 = M2*D1 ! Hadamard product
     145     6403274 :          M4 = MATMUL(MATMUL(H_evecs, M3), TRANSPOSE(H_evecs))
     146             : 
     147             :          ! gradient contribution from penalty terms
     148        2471 :          IF (PRESENT(penalty)) THEN
     149        9748 :             ALLOCATE (D2(n, n))
     150      155955 :             D2 = 0.0_dp
     151       18854 :             DO i = 1, n
     152      155955 :             DO j = 1, n
     153      137101 :                IF (i <= m .EQV. j <= m) CYCLE
     154       52992 :                diff = H_evals(i) - H_evals(j)
     155      153518 :                D2(i, i) = D2(i, i) - 2.0_dp*alpha*diff/beta**2*EXP(-(diff/beta)**2)
     156             :             END DO
     157             :             END DO
     158     9188927 :             M4 = M4 + MATMUL(MATMUL(H_evecs, D2), TRANSPOSE(H_evecs))
     159        2437 :             DEALLOCATE (D2)
     160             :          END IF
     161             : 
     162             :          ! dH / dV, return to non-orthonormal basis
     163        9884 :          ALLOCATE (M5(n, n))
     164    12028796 :          M5 = MATMUL(MATMUL(block_N, M4), block_N)
     165             : 
     166             :          ! add regularization gradient
     167        2471 :          IF (PRESENT(penalty)) &
     168      309473 :             M5 = M5 + 2.0_dp*pao%regularization*V
     169             : 
     170             :          ! symmetrize
     171      311595 :          G = 0.5_dp*(M5 + TRANSPOSE(M5)) ! the final gradient
     172             : 
     173        2471 :          DEALLOCATE (D1, M2, M3, M4, M5)
     174             :       END IF
     175             : 
     176       15520 :       DEALLOCATE (H, H_evals, H_evecs)
     177             : 
     178       15520 :       CALL timestop(handle)
     179       31040 :    END SUBROUTINE pao_calc_U_block_fock
     180             : 
     181       25370 : END MODULE pao_param_fock