pao_param_fock.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Common framework for using eigenvectors of a Fock matrix as PAO basis.
!> \author Ole Schuett
! **************************************************************************************************
MODULE pao_param_fock
   USE cp_dbcsr_api,                    ONLY: dbcsr_get_block_p,&
                                              dbcsr_get_info
   USE kinds,                           ONLY: dp
   USE mathlib,                         ONLY: diamat_all
   USE pao_types,                       ONLY: pao_env_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_param_fock'
 
   PUBLIC :: pao_calc_u_block_fock
 
CONTAINS
 
! **************************************************************************************************
!> \brief Calculate new matrix U and optinally its gradient G
!> \param pao ...
!> \param iatom ...
!> \param V ...
!> \param U ...
!> \param penalty ...
!> \param gap ...
!> \param evals ...
!> \param M1 ...
!> \param G ...
! **************************************************************************************************
   SUBROUTINE pao_calc_u_block_fock(pao, iatom, V, U, penalty, gap, evals, M1, G)
      TYPE(pao_env_type), POINTER                        :: pao
      INTEGER, INTENT(IN)                                :: iatom
      REAL(dp), DIMENSION(:, :), POINTER                 :: v, u
      REAL(dp), INTENT(INOUT), OPTIONAL                  :: penalty
      REAL(dp), INTENT(OUT)                              :: gap
      REAL(dp), DIMENSION(:), INTENT(OUT), OPTIONAL      :: evals
      REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER       :: m1, g
 
      CHARACTER(len=*), PARAMETER :: routinen = 'pao_calc_U_block_fock'
 
      INTEGER                                            :: handle, i, j, m, n
      INTEGER, DIMENSION(:), POINTER                     :: blk_sizes_pao, blk_sizes_pri
      LOGICAL                                            :: found
      REAL(dp)                                           :: alpha, beta, denom, diff
      REAL(dp), DIMENSION(:), POINTER                    :: h_evals
      REAL(dp), DIMENSION(:, :), POINTER                 :: block_n, d1, d2, h, h0, h_evecs, m2, m3, &
                                                            m4, m5
 
      CALL timeset(routinen, handle)
 
      CALL dbcsr_get_block_p(matrix=pao%matrix_H0, row=iatom, col=iatom, block=h0, found=found)
      cpassert(ASSOCIATED(h0))
      CALL dbcsr_get_block_p(matrix=pao%matrix_N_diag, row=iatom, col=iatom, block=block_n, found=found)
      cpassert(ASSOCIATED(block_n))
      IF (maxval(abs(v - transpose(v))) > 1e-14_dp) cpabort("Expect symmetric matrix")
 
      ! figure out basis sizes
      CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
      n = blk_sizes_pri(iatom) ! size of primary basis
      m = blk_sizes_pao(iatom) ! size of pao basis
 
      ! calculate H in the orthonormal basis
      ALLOCATE (h(n, n))
      h = matmul(matmul(block_n, h0 + v), block_n)
 
      ! diagonalize H
      ALLOCATE (h_evals(n), h_evecs(n, n))
      h_evecs = h
      CALL diamat_all(h_evecs, h_evals)
 
      ! the eigenvectors of H become the rotation matrix U
      u = h_evecs
 
      ! copy eigenvectors around the gap from H_evals into evals array
      IF (PRESENT(evals)) THEN
         cpassert(mod(SIZE(evals), 2) == 0) ! gap will be exactely in the middle
         i = SIZE(evals)/2
         j = min(m, i)
         evals(1 + i - j:i) = h_evals(1 + m - j:m) ! eigenvalues below gap
         j = min(n - m, i)
         evals(i:i + j) = h_evals(m:m + j) ! eigenvalues above gap
      END IF
 
      ! calculate homo-lumo gap (it's useful for detecting numerical issues)
      gap = huge(dp)
      IF (m < n) & ! catch special case n==m
         gap = h_evals(m + 1) - h_evals(m)
 
      IF (PRESENT(penalty)) THEN
         ! penalty terms: occupied and virtual eigenvalues repel each other
         alpha = pao%penalty_strength
         beta = pao%penalty_dist
         DO i = 1, m
         DO j = m + 1, n
            diff = h_evals(i) - h_evals(j)
            penalty = penalty + alpha*exp(-(diff/beta)**2)
         END DO
         END DO
 
         ! regularization energy
         penalty = penalty + pao%regularization*sum(v**2)
      END IF
 
      IF (PRESENT(g)) THEN ! TURNING POINT (if calc grad) -------------------------
 
         cpassert(PRESENT(m1))
 
         ! calculate derivatives between eigenvectors of H
         ALLOCATE (d1(n, n), m2(n, n), m3(n, n), m4(n, n))
         DO i = 1, n
         DO j = 1, n
            ! ignore changes among occupied or virtual eigenvectors
            ! They will get filtered out by M2*D1 anyways, however this early
            ! intervention might stabilize numerics in the case of level-crossings.
            IF (i <= m .EQV. j <= m) THEN
               d1(i, j) = 0.0_dp
            ELSE
               denom = h_evals(i) - h_evals(j)
               IF (abs(denom) > 1e-9_dp) THEN ! avoid division by zero
                  d1(i, j) = 1.0_dp/denom
               ELSE
                  d1(i, j) = sign(1e+9_dp, denom)
               END IF
            END IF
         END DO
         END DO
         IF (ASSOCIATED(m1)) THEN
            m2 = matmul(transpose(m1), h_evecs)
         ELSE
            m2 = 0.0_dp
         END IF
         m3 = m2*d1 ! Hadamard product
         m4 = matmul(matmul(h_evecs, m3), transpose(h_evecs))
 
         ! gradient contribution from penalty terms
         IF (PRESENT(penalty)) THEN
            ALLOCATE (d2(n, n))
            d2 = 0.0_dp
            DO i = 1, n
            DO j = 1, n
               IF (i <= m .EQV. j <= m) cycle
               diff = h_evals(i) - h_evals(j)
               d2(i, i) = d2(i, i) - 2.0_dp*alpha*diff/beta**2*exp(-(diff/beta)**2)
            END DO
            END DO
            m4 = m4 + matmul(matmul(h_evecs, d2), transpose(h_evecs))
            DEALLOCATE (d2)
         END IF
 
         ! dH / dV, return to non-orthonormal basis
         ALLOCATE (m5(n, n))
         m5 = matmul(matmul(block_n, m4), block_n)
 
         ! add regularization gradient
         IF (PRESENT(penalty)) &
            m5 = m5 + 2.0_dp*pao%regularization*v
 
         ! symmetrize
         g = 0.5_dp*(m5 + transpose(m5)) ! the final gradient
 
         DEALLOCATE (d1, m2, m3, m4, m5)
      END IF
 
      DEALLOCATE (h, h_evals, h_evecs)
 
      CALL timestop(handle)
   END SUBROUTINE pao_calc_u_block_fock
 
END MODULE pao_param_fock