pao_param_exp.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 Original matrix exponential parametrization
!> \author Ole Schuett
! **************************************************************************************************
MODULE pao_param_exp
   USE basis_set_types,                 ONLY: gto_basis_set_type
   USE cp_dbcsr_api,                    ONLY: &
        dbcsr_create, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
        dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
        dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type
   USE cp_dbcsr_contrib,                ONLY: dbcsr_reserve_diag_blocks
   USE dm_ls_scf_types,                 ONLY: ls_scf_env_type
   USE kinds,                           ONLY: dp
   USE mathlib,                         ONLY: diag_antisym,&
                                              diamat_all
   USE pao_param_methods,               ONLY: pao_calc_ab_from_u,&
                                              pao_calc_grad_lnv_wrt_u
   USE pao_potentials,                  ONLY: pao_guess_initial_potential
   USE pao_types,                       ONLY: pao_env_type
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
   USE qs_kind_types,                   ONLY: get_qs_kind,&
                                              qs_kind_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_param_exp'
 
   PUBLIC :: pao_param_init_exp, pao_param_finalize_exp, pao_calc_ab_exp
   PUBLIC :: pao_param_count_exp, pao_param_initguess_exp
 
CONTAINS
 
! **************************************************************************************************
!> \brief Initialize matrix exponential parametrization
!> \param pao ...
!> \param qs_env ...
! **************************************************************************************************
   SUBROUTINE pao_param_init_exp(pao, qs_env)
      TYPE(pao_env_type), POINTER                        :: pao
      TYPE(qs_environment_type), POINTER                 :: qs_env
 
      CHARACTER(len=*), PARAMETER :: routinen = 'pao_param_init_exp'
 
      INTEGER                                            :: acol, arow, handle, iatom, n
      LOGICAL                                            :: found
      REAL(dp), DIMENSION(:), POINTER                    :: h_evals
      REAL(dp), DIMENSION(:, :), POINTER                 :: block_h, block_h0, block_n, block_u0, &
                                                            block_v0, h_evecs
      TYPE(dbcsr_iterator_type)                          :: iter
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
 
      CALL timeset(routinen, handle)
 
      CALL get_qs_env(qs_env, matrix_s=matrix_s)
 
      ! allocate matrix_U0
      CALL dbcsr_create(pao%matrix_U0, &
                        name="PAO matrix_U0", &
                        matrix_type="N", &
                        dist=pao%diag_distribution, &
                        template=matrix_s(1)%matrix)
      CALL dbcsr_reserve_diag_blocks(pao%matrix_U0)
 
      ! diagonalize each block of H0 and store eigenvectors in U0
!$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env) &
!$OMP PRIVATE(iter,arow,acol,iatom,N,found,block_H0,block_V0,block_N,block_H,block_U0,H_evecs,H_evals)
      CALL dbcsr_iterator_start(iter, pao%matrix_U0)
      DO WHILE (dbcsr_iterator_blocks_left(iter))
         CALL dbcsr_iterator_next_block(iter, arow, acol, block_u0)
         iatom = arow; cpassert(arow == acol)
         CALL dbcsr_get_block_p(matrix=pao%matrix_H0, row=iatom, col=iatom, block=block_h0, found=found)
         CALL dbcsr_get_block_p(matrix=pao%matrix_N_diag, row=iatom, col=iatom, block=block_n, found=found)
         cpassert(ASSOCIATED(block_h0) .AND. ASSOCIATED(block_n))
         n = SIZE(block_u0, 1)
 
         ALLOCATE (block_v0(n, n))
         CALL pao_guess_initial_potential(qs_env, iatom, block_v0)
 
         ! construct H
         ALLOCATE (block_h(n, n))
         block_h = matmul(matmul(block_n, block_h0 + block_v0), block_n) ! transform into orthonormal basis
 
         ! diagonalize H
         ALLOCATE (h_evecs(n, n), h_evals(n))
         h_evecs = block_h
         CALL diamat_all(h_evecs, h_evals)
 
         ! use eigenvectors as initial guess
         block_u0 = h_evecs
 
         DEALLOCATE (block_h, h_evecs, h_evals, block_v0)
      END DO
      CALL dbcsr_iterator_stop(iter)
!$OMP END PARALLEL
 
      IF (pao%precondition) &
         cpabort("PAO preconditioning not supported for selected parametrization.")
 
      CALL timestop(handle)
   END SUBROUTINE pao_param_init_exp
 
! **************************************************************************************************
!> \brief Finalize exponential parametrization
!> \param pao ...
! **************************************************************************************************
   SUBROUTINE pao_param_finalize_exp(pao)
      TYPE(pao_env_type), POINTER                        :: pao
 
      CALL dbcsr_release(pao%matrix_U0)
 
   END SUBROUTINE pao_param_finalize_exp
 
! **************************************************************************************************
!> \brief Returns the number of parameters for given atomic kind
!> \param qs_env ...
!> \param ikind ...
!> \param nparams ...
! **************************************************************************************************
   SUBROUTINE pao_param_count_exp(qs_env, ikind, nparams)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      INTEGER, INTENT(IN)                                :: ikind
      INTEGER, INTENT(OUT)                               :: nparams
 
      INTEGER                                            :: cols, pao_basis_size, pri_basis_size, &
                                                            rows
      TYPE(gto_basis_set_type), POINTER                  :: basis_set
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
 
      CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
      CALL get_qs_kind(qs_kind_set(ikind), &
                       basis_set=basis_set, &
                       pao_basis_size=pao_basis_size)
      pri_basis_size = basis_set%nsgf
 
      ! we only consider rotations between occupied and virtuals
      rows = pao_basis_size
      cols = pri_basis_size - pao_basis_size
      nparams = rows*cols
 
   END SUBROUTINE pao_param_count_exp
 
! **************************************************************************************************
!> \brief Fills matrix_X with an initial guess
!> \param pao ...
! **************************************************************************************************
   SUBROUTINE pao_param_initguess_exp(pao)
      TYPE(pao_env_type), POINTER                        :: pao
 
      CALL dbcsr_set(pao%matrix_X, 0.0_dp) ! actual initial guess is matrix_U0
 
   END SUBROUTINE pao_param_initguess_exp
 
! **************************************************************************************************
!> \brief Takes current matrix_X and calculates the matrices A and B.
!> \param pao ...
!> \param qs_env ...
!> \param ls_scf_env ...
!> \param gradient ...
! **************************************************************************************************
   SUBROUTINE pao_calc_ab_exp(pao, qs_env, ls_scf_env, gradient)
      TYPE(pao_env_type), POINTER                        :: pao
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(ls_scf_env_type), TARGET                      :: ls_scf_env
      LOGICAL, INTENT(IN)                                :: gradient
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'pao_calc_AB_exp'
 
      INTEGER                                            :: handle
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
      TYPE(dbcsr_type)                                   :: matrix_m, matrix_u
 
      CALL timeset(routinen, handle)
      CALL get_qs_env(qs_env, matrix_s=matrix_s)
      CALL dbcsr_create(matrix_u, matrix_type="N", dist=pao%diag_distribution, template=matrix_s(1)%matrix)
      CALL dbcsr_reserve_diag_blocks(matrix_u)
 
      !TODO: move this condition into pao_calc_U, use matrix_N as template
      IF (gradient) THEN
         CALL pao_calc_grad_lnv_wrt_u(qs_env, ls_scf_env, matrix_m)
         CALL pao_calc_u_exp(pao, matrix_u, matrix_m, pao%matrix_G)
         CALL dbcsr_release(matrix_m)
      ELSE
         CALL pao_calc_u_exp(pao, matrix_u)
      END IF
 
      CALL pao_calc_ab_from_u(pao, qs_env, ls_scf_env, matrix_u)
      CALL dbcsr_release(matrix_u)
      CALL timestop(handle)
   END SUBROUTINE pao_calc_ab_exp
 
! **************************************************************************************************
!> \brief Calculate new matrix U and optionally its gradient G
!> \param pao ...
!> \param matrix_U ...
!> \param matrix_M ...
!> \param matrix_G ...
! **************************************************************************************************
   SUBROUTINE pao_calc_u_exp(pao, matrix_U, matrix_M, matrix_G)
      TYPE(pao_env_type), POINTER                        :: pao
      TYPE(dbcsr_type)                                   :: matrix_u
      TYPE(dbcsr_type), OPTIONAL                         :: matrix_m, matrix_g
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'pao_calc_U_exp'
 
      COMPLEX(dp)                                        :: denom
      COMPLEX(dp), DIMENSION(:), POINTER                 :: evals
      COMPLEX(dp), DIMENSION(:, :), POINTER              :: block_d, evecs
      INTEGER                                            :: acol, arow, handle, i, iatom, j, k, m, &
                                                            n, nparams
      INTEGER, DIMENSION(:), POINTER                     :: blk_sizes_pao, blk_sizes_pri
      LOGICAL                                            :: found
      REAL(dp), DIMENSION(:, :), POINTER                 :: block_g, block_g_full, block_m, &
                                                            block_tmp, block_u, block_u0, block_x, &
                                                            block_x_full
      TYPE(dbcsr_iterator_type)                          :: iter
 
      CALL timeset(routinen, handle)
 
      CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
 
!$OMP PARALLEL DEFAULT(NONE) SHARED(pao,matrix_U,matrix_M,matrix_G,blk_sizes_pri,blk_sizes_pao) &
!$OMP PRIVATE(iter,arow,acol,iatom,N,M,nparams,i,j,k,found) &
!$OMP PRIVATE(block_X,block_U,block_U0,block_X_full,evals,evecs) &
!$OMP PRIVATE(block_M,block_G,block_D,block_tmp,block_G_full,denom)
      CALL dbcsr_iterator_start(iter, pao%matrix_X)
      DO WHILE (dbcsr_iterator_blocks_left(iter))
         CALL dbcsr_iterator_next_block(iter, arow, acol, block_x)
         iatom = arow; cpassert(arow == acol)
         CALL dbcsr_get_block_p(matrix=matrix_u, row=iatom, col=iatom, block=block_u, found=found)
         cpassert(ASSOCIATED(block_u))
         CALL dbcsr_get_block_p(matrix=pao%matrix_U0, row=iatom, col=iatom, block=block_u0, found=found)
         cpassert(ASSOCIATED(block_u0))
 
         n = blk_sizes_pri(iatom) ! size of primary basis
         m = blk_sizes_pao(iatom) ! size of pao basis
         nparams = SIZE(block_x, 1)
 
         ! block_X stores only rotations between occupied and virtuals
         ! hence, we first have to build the full anti-symmetric exponent block
         ALLOCATE (block_x_full(n, n))
         block_x_full(:, :) = 0.0_dp
         DO i = 1, nparams
            block_x_full(mod(i - 1, m) + 1, m + (i - 1)/m + 1) = +block_x(i, 1)
            block_x_full(m + (i - 1)/m + 1, mod(i - 1, m) + 1) = -block_x(i, 1)
         END DO
 
         ! diagonalize block_X_full
         ALLOCATE (evals(n), evecs(n, n))
         CALL diag_antisym(block_x_full, evecs, evals)
 
         ! construct rotation matrix
         block_u(:, :) = 0.0_dp
         DO k = 1, n
            DO i = 1, n
               DO j = 1, n
                  block_u(i, j) = block_u(i, j) + real(exp(evals(k))*evecs(i, k)*conjg(evecs(j, k)), dp)
               END DO
            END DO
         END DO
 
         block_u = matmul(block_u0, block_u) ! prepend initial guess rotation
 
         ! TURNING POINT (if calc grad) ------------------------------------------
         IF (PRESENT(matrix_g)) THEN
            cpassert(PRESENT(matrix_m))
 
            CALL dbcsr_get_block_p(matrix=pao%matrix_G, row=iatom, col=iatom, block=block_g, found=found)
            cpassert(ASSOCIATED(block_g))
            CALL dbcsr_get_block_p(matrix=matrix_m, row=iatom, col=iatom, block=block_m, found=found)
            ! don't check ASSOCIATED(block_M), it might have been filtered out.
 
            ALLOCATE (block_d(n, n), block_tmp(n, n), block_g_full(n, n))
            DO i = 1, n
               DO j = 1, n
                  denom = evals(i) - evals(j)
                  IF (i == j) THEN
                     block_d(i, i) = exp(evals(i)) ! diagonal elements
                  ELSE IF (abs(denom) > 1e-10_dp) THEN
                     block_d(i, j) = (exp(evals(i)) - exp(evals(j)))/denom
                  ELSE
                     block_d(i, j) = 1.0_dp ! limit according to L'Hospital's rule
                  END IF
               END DO
            END DO
 
            IF (ASSOCIATED(block_m)) THEN
               block_tmp = matmul(transpose(block_u0), block_m)
            ELSE
               block_tmp = 0.0_dp
            END IF
            block_g_full = fold_derivatives(block_tmp, block_d, evecs)
 
            ! return only gradient for rotations between occupied and virtuals
            DO i = 1, nparams
               block_g(i, 1) = 2.0_dp*block_g_full(mod(i - 1, m) + 1, m + (i - 1)/m + 1)
            END DO
 
            DEALLOCATE (block_d, block_tmp, block_g_full)
         END IF
 
         DEALLOCATE (block_x_full, evals, evecs)
 
      END DO
      CALL dbcsr_iterator_stop(iter)
!$OMP END PARALLEL
 
      CALL timestop(handle)
   END SUBROUTINE pao_calc_u_exp
 
! **************************************************************************************************
!> \brief Helper routine, for calculating derivatives
!> \param M ...
!> \param D ...
!> \param R ...
!> \return ...
! **************************************************************************************************
   FUNCTION fold_derivatives(M, D, R) RESULT(G)
      REAL(dp), DIMENSION(:, :), INTENT(IN)              :: m
      COMPLEX(dp), DIMENSION(:, :), INTENT(IN)           :: d, r
      REAL(dp), DIMENSION(SIZE(M, 1), SIZE(M, 1))        :: g
 
      COMPLEX(dp), DIMENSION(:, :), POINTER              :: f, rf, rm, rmr
      INTEGER                                            :: n
      REAL(dp), DIMENSION(:, :), POINTER                 :: rfr
 
      n = SIZE(m, 1)
 
      ALLOCATE (rm(n, n), rmr(n, n), f(n, n), rf(n, n), rfr(n, n))
 
      rm = matmul(transpose(conjg(r)), transpose(m))
      rmr = matmul(rm, r)
      f = rmr*d !Hadamard product
      rf = matmul(r, f)
      rfr = real(matmul(rf, transpose(conjg(r))), dp)
 
      ! gradient dE/dX has to be anti-symmetric
      g = 0.5_dp*(transpose(rfr) - rfr)
 
      DEALLOCATE (rm, rmr, f, rf, rfr)
   END FUNCTION fold_derivatives
 
END MODULE pao_param_exp