LCOV - code coverage report
Current view: top level - src - pao_param_exp.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 98.6 % 70 69
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 7 7

            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 Original matrix exponential parametrization
      10              : !> \author Ole Schuett
      11              : ! **************************************************************************************************
      12              : MODULE pao_param_exp
      13              :    USE basis_set_types,                 ONLY: gto_basis_set_type
      14              :    USE cp_dbcsr_api,                    ONLY: &
      15              :         dbcsr_create, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, &
      16              :         dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
      17              :         dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type
      18              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_reserve_diag_blocks
      19              :    USE dm_ls_scf_types,                 ONLY: ls_scf_env_type
      20              :    USE kinds,                           ONLY: dp
      21              :    USE mathlib,                         ONLY: diag_antisym,&
      22              :                                               diamat_all
      23              :    USE pao_param_methods,               ONLY: pao_calc_AB_from_U,&
      24              :                                               pao_calc_grad_lnv_wrt_U
      25              :    USE pao_potentials,                  ONLY: pao_guess_initial_potential
      26              :    USE pao_types,                       ONLY: pao_env_type
      27              :    USE qs_environment_types,            ONLY: get_qs_env,&
      28              :                                               qs_environment_type
      29              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      30              :                                               qs_kind_type
      31              : #include "./base/base_uses.f90"
      32              : 
      33              :    IMPLICIT NONE
      34              : 
      35              :    PRIVATE
      36              : 
      37              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pao_param_exp'
      38              : 
      39              :    PUBLIC :: pao_param_init_exp, pao_param_finalize_exp, pao_calc_AB_exp
      40              :    PUBLIC :: pao_param_count_exp, pao_param_initguess_exp
      41              : 
      42              : CONTAINS
      43              : 
      44              : ! **************************************************************************************************
      45              : !> \brief Initialize matrix exponential parametrization
      46              : !> \param pao ...
      47              : !> \param qs_env ...
      48              : ! **************************************************************************************************
      49           24 :    SUBROUTINE pao_param_init_exp(pao, qs_env)
      50              :       TYPE(pao_env_type), POINTER                        :: pao
      51              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      52              : 
      53              :       CHARACTER(len=*), PARAMETER :: routineN = 'pao_param_init_exp'
      54              : 
      55              :       INTEGER                                            :: acol, arow, handle, iatom, N
      56              :       LOGICAL                                            :: found
      57           24 :       REAL(dp), DIMENSION(:), POINTER                    :: H_evals
      58           24 :       REAL(dp), DIMENSION(:, :), POINTER                 :: block_H, block_H0, block_N, block_U0, &
      59           24 :                                                             block_V0, H_evecs
      60              :       TYPE(dbcsr_iterator_type)                          :: iter
      61           24 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
      62              : 
      63           24 :       CALL timeset(routineN, handle)
      64              : 
      65           24 :       CALL get_qs_env(qs_env, matrix_s=matrix_s)
      66              : 
      67              :       ! allocate matrix_U0
      68              :       CALL dbcsr_create(pao%matrix_U0, &
      69              :                         name="PAO matrix_U0", &
      70              :                         matrix_type="N", &
      71              :                         dist=pao%diag_distribution, &
      72           24 :                         template=matrix_s(1)%matrix)
      73           24 :       CALL dbcsr_reserve_diag_blocks(pao%matrix_U0)
      74              : 
      75              :       ! diagonalize each block of H0 and store eigenvectors in U0
      76              : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,qs_env) &
      77           24 : !$OMP PRIVATE(iter,arow,acol,iatom,N,found,block_H0,block_V0,block_N,block_H,block_U0,H_evecs,H_evals)
      78              :       CALL dbcsr_iterator_start(iter, pao%matrix_U0)
      79              :       DO WHILE (dbcsr_iterator_blocks_left(iter))
      80              :          CALL dbcsr_iterator_next_block(iter, arow, acol, block_U0)
      81              :          iatom = arow; CPASSERT(arow == acol)
      82              :          CALL dbcsr_get_block_p(matrix=pao%matrix_H0, row=iatom, col=iatom, block=block_H0, found=found)
      83              :          CALL dbcsr_get_block_p(matrix=pao%matrix_N_diag, row=iatom, col=iatom, block=block_N, found=found)
      84              :          CPASSERT(ASSOCIATED(block_H0) .AND. ASSOCIATED(block_N))
      85              :          N = SIZE(block_U0, 1)
      86              : 
      87              :          ALLOCATE (block_V0(N, N))
      88              :          CALL pao_guess_initial_potential(qs_env, iatom, block_V0)
      89              : 
      90              :          ! construct H
      91              :          ALLOCATE (block_H(N, N))
      92              :          block_H = MATMUL(MATMUL(block_N, block_H0 + block_V0), block_N) ! transform into orthonormal basis
      93              : 
      94              :          ! diagonalize H
      95              :          ALLOCATE (H_evecs(N, N), H_evals(N))
      96              :          H_evecs = block_H
      97              :          CALL diamat_all(H_evecs, H_evals)
      98              : 
      99              :          ! use eigenvectors as initial guess
     100              :          block_U0 = H_evecs
     101              : 
     102              :          DEALLOCATE (block_H, H_evecs, H_evals, block_V0)
     103              :       END DO
     104              :       CALL dbcsr_iterator_stop(iter)
     105              : !$OMP END PARALLEL
     106              : 
     107           24 :       IF (pao%precondition) &
     108            0 :          CPABORT("PAO preconditioning not supported for selected parametrization.")
     109              : 
     110           24 :       CALL timestop(handle)
     111           24 :    END SUBROUTINE pao_param_init_exp
     112              : 
     113              : ! **************************************************************************************************
     114              : !> \brief Finalize exponential parametrization
     115              : !> \param pao ...
     116              : ! **************************************************************************************************
     117           24 :    SUBROUTINE pao_param_finalize_exp(pao)
     118              :       TYPE(pao_env_type), POINTER                        :: pao
     119              : 
     120           24 :       CALL dbcsr_release(pao%matrix_U0)
     121              : 
     122           24 :    END SUBROUTINE pao_param_finalize_exp
     123              : 
     124              : ! **************************************************************************************************
     125              : !> \brief Returns the number of parameters for given atomic kind
     126              : !> \param qs_env ...
     127              : !> \param ikind ...
     128              : !> \param nparams ...
     129              : ! **************************************************************************************************
     130          128 :    SUBROUTINE pao_param_count_exp(qs_env, ikind, nparams)
     131              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     132              :       INTEGER, INTENT(IN)                                :: ikind
     133              :       INTEGER, INTENT(OUT)                               :: nparams
     134              : 
     135              :       INTEGER                                            :: cols, pao_basis_size, pri_basis_size, &
     136              :                                                             rows
     137              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set
     138           64 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     139              : 
     140           64 :       CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set)
     141              :       CALL get_qs_kind(qs_kind_set(ikind), &
     142              :                        basis_set=basis_set, &
     143           64 :                        pao_basis_size=pao_basis_size)
     144           64 :       pri_basis_size = basis_set%nsgf
     145              : 
     146              :       ! we only consider rotations between occupied and virtuals
     147           64 :       rows = pao_basis_size
     148           64 :       cols = pri_basis_size - pao_basis_size
     149           64 :       nparams = rows*cols
     150              : 
     151           64 :    END SUBROUTINE pao_param_count_exp
     152              : 
     153              : ! **************************************************************************************************
     154              : !> \brief Fills matrix_X with an initial guess
     155              : !> \param pao ...
     156              : ! **************************************************************************************************
     157           14 :    SUBROUTINE pao_param_initguess_exp(pao)
     158              :       TYPE(pao_env_type), POINTER                        :: pao
     159              : 
     160           14 :       CALL dbcsr_set(pao%matrix_X, 0.0_dp) ! actual initial guess is matrix_U0
     161              : 
     162           14 :    END SUBROUTINE pao_param_initguess_exp
     163              : 
     164              : ! **************************************************************************************************
     165              : !> \brief Takes current matrix_X and calculates the matrices A and B.
     166              : !> \param pao ...
     167              : !> \param qs_env ...
     168              : !> \param ls_scf_env ...
     169              : !> \param gradient ...
     170              : ! **************************************************************************************************
     171         2710 :    SUBROUTINE pao_calc_AB_exp(pao, qs_env, ls_scf_env, gradient)
     172              :       TYPE(pao_env_type), POINTER                        :: pao
     173              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     174              :       TYPE(ls_scf_env_type), TARGET                      :: ls_scf_env
     175              :       LOGICAL, INTENT(IN)                                :: gradient
     176              : 
     177              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'pao_calc_AB_exp'
     178              : 
     179              :       INTEGER                                            :: handle
     180         2710 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     181              :       TYPE(dbcsr_type)                                   :: matrix_M, matrix_U
     182              : 
     183         2710 :       CALL timeset(routineN, handle)
     184         2710 :       CALL get_qs_env(qs_env, matrix_s=matrix_s)
     185         2710 :       CALL dbcsr_create(matrix_U, matrix_type="N", dist=pao%diag_distribution, template=matrix_s(1)%matrix)
     186         2710 :       CALL dbcsr_reserve_diag_blocks(matrix_U)
     187              : 
     188              :       !TODO: move this condition into pao_calc_U, use matrix_N as template
     189         2710 :       IF (gradient) THEN
     190          488 :          CALL pao_calc_grad_lnv_wrt_U(qs_env, ls_scf_env, matrix_M)
     191          488 :          CALL pao_calc_U_exp(pao, matrix_U, matrix_M, pao%matrix_G)
     192          488 :          CALL dbcsr_release(matrix_M)
     193              :       ELSE
     194         2222 :          CALL pao_calc_U_exp(pao, matrix_U)
     195              :       END IF
     196              : 
     197         2710 :       CALL pao_calc_AB_from_U(pao, qs_env, ls_scf_env, matrix_U)
     198         2710 :       CALL dbcsr_release(matrix_U)
     199         2710 :       CALL timestop(handle)
     200         2710 :    END SUBROUTINE pao_calc_AB_exp
     201              : 
     202              : ! **************************************************************************************************
     203              : !> \brief Calculate new matrix U and optionally its gradient G
     204              : !> \param pao ...
     205              : !> \param matrix_U ...
     206              : !> \param matrix_M ...
     207              : !> \param matrix_G ...
     208              : ! **************************************************************************************************
     209         2710 :    SUBROUTINE pao_calc_U_exp(pao, matrix_U, matrix_M, matrix_G)
     210              :       TYPE(pao_env_type), POINTER                        :: pao
     211              :       TYPE(dbcsr_type)                                   :: matrix_U
     212              :       TYPE(dbcsr_type), OPTIONAL                         :: matrix_M, matrix_G
     213              : 
     214              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'pao_calc_U_exp'
     215              : 
     216              :       COMPLEX(dp)                                        :: denom
     217         2710 :       COMPLEX(dp), DIMENSION(:), POINTER                 :: evals
     218         2710 :       COMPLEX(dp), DIMENSION(:, :), POINTER              :: block_D, evecs
     219              :       INTEGER                                            :: acol, arow, handle, i, iatom, j, k, M, &
     220              :                                                             N, nparams
     221         2710 :       INTEGER, DIMENSION(:), POINTER                     :: blk_sizes_pao, blk_sizes_pri
     222              :       LOGICAL                                            :: found
     223         2710 :       REAL(dp), DIMENSION(:, :), POINTER                 :: block_G, block_G_full, block_M, &
     224         2710 :                                                             block_tmp, block_U, block_U0, block_X, &
     225         2710 :                                                             block_X_full
     226              :       TYPE(dbcsr_iterator_type)                          :: iter
     227              : 
     228         2710 :       CALL timeset(routineN, handle)
     229              : 
     230         2710 :       CALL dbcsr_get_info(pao%matrix_Y, row_blk_size=blk_sizes_pri, col_blk_size=blk_sizes_pao)
     231              : 
     232              : !$OMP PARALLEL DEFAULT(NONE) SHARED(pao,matrix_U,matrix_M,matrix_G,blk_sizes_pri,blk_sizes_pao) &
     233              : !$OMP PRIVATE(iter,arow,acol,iatom,N,M,nparams,i,j,k,found) &
     234              : !$OMP PRIVATE(block_X,block_U,block_U0,block_X_full,evals,evecs) &
     235         2710 : !$OMP PRIVATE(block_M,block_G,block_D,block_tmp,block_G_full,denom)
     236              :       CALL dbcsr_iterator_start(iter, pao%matrix_X)
     237              :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     238              :          CALL dbcsr_iterator_next_block(iter, arow, acol, block_X)
     239              :          iatom = arow; CPASSERT(arow == acol)
     240              :          CALL dbcsr_get_block_p(matrix=matrix_U, row=iatom, col=iatom, block=block_U, found=found)
     241              :          CPASSERT(ASSOCIATED(block_U))
     242              :          CALL dbcsr_get_block_p(matrix=pao%matrix_U0, row=iatom, col=iatom, block=block_U0, found=found)
     243              :          CPASSERT(ASSOCIATED(block_U0))
     244              : 
     245              :          N = blk_sizes_pri(iatom) ! size of primary basis
     246              :          M = blk_sizes_pao(iatom) ! size of pao basis
     247              :          nparams = SIZE(block_X, 1)
     248              : 
     249              :          ! block_X stores only rotations between occupied and virtuals
     250              :          ! hence, we first have to build the full anti-symmetric exponent block
     251              :          ALLOCATE (block_X_full(N, N))
     252              :          block_X_full(:, :) = 0.0_dp
     253              :          DO i = 1, nparams
     254              :             block_X_full(MOD(i - 1, M) + 1, M + (i - 1)/M + 1) = +block_X(i, 1)
     255              :             block_X_full(M + (i - 1)/M + 1, MOD(i - 1, M) + 1) = -block_X(i, 1)
     256              :          END DO
     257              : 
     258              :          ! diagonalize block_X_full
     259              :          ALLOCATE (evals(N), evecs(N, N))
     260              :          CALL diag_antisym(block_X_full, evecs, evals)
     261              : 
     262              :          ! construct rotation matrix
     263              :          block_U(:, :) = 0.0_dp
     264              :          DO k = 1, N
     265              :             DO i = 1, N
     266              :                DO j = 1, N
     267              :                   block_U(i, j) = block_U(i, j) + REAL(EXP(evals(k))*evecs(i, k)*CONJG(evecs(j, k)), dp)
     268              :                END DO
     269              :             END DO
     270              :          END DO
     271              : 
     272              :          block_U = MATMUL(block_U0, block_U) ! prepend initial guess rotation
     273              : 
     274              :          ! TURNING POINT (if calc grad) ------------------------------------------
     275              :          IF (PRESENT(matrix_G)) THEN
     276              :             CPASSERT(PRESENT(matrix_M))
     277              : 
     278              :             CALL dbcsr_get_block_p(matrix=pao%matrix_G, row=iatom, col=iatom, block=block_G, found=found)
     279              :             CPASSERT(ASSOCIATED(block_G))
     280              :             CALL dbcsr_get_block_p(matrix=matrix_M, row=iatom, col=iatom, block=block_M, found=found)
     281              :             ! don't check ASSOCIATED(block_M), it might have been filtered out.
     282              : 
     283              :             ALLOCATE (block_D(N, N), block_tmp(N, N), block_G_full(N, N))
     284              :             DO i = 1, N
     285              :                DO j = 1, N
     286              :                   denom = evals(i) - evals(j)
     287              :                   IF (i == j) THEN
     288              :                      block_D(i, i) = EXP(evals(i)) ! diagonal elements
     289              :                   ELSE IF (ABS(denom) > 1e-10_dp) THEN
     290              :                      block_D(i, j) = (EXP(evals(i)) - EXP(evals(j)))/denom
     291              :                   ELSE
     292              :                      block_D(i, j) = 1.0_dp ! limit according to L'Hospital's rule
     293              :                   END IF
     294              :                END DO
     295              :             END DO
     296              : 
     297              :             IF (ASSOCIATED(block_M)) THEN
     298              :                block_tmp = MATMUL(TRANSPOSE(block_U0), block_M)
     299              :             ELSE
     300              :                block_tmp = 0.0_dp
     301              :             END IF
     302              :             block_G_full = fold_derivatives(block_tmp, block_D, evecs)
     303              : 
     304              :             ! return only gradient for rotations between occupied and virtuals
     305              :             DO i = 1, nparams
     306              :                block_G(i, 1) = 2.0_dp*block_G_full(MOD(i - 1, M) + 1, M + (i - 1)/M + 1)
     307              :             END DO
     308              : 
     309              :             DEALLOCATE (block_D, block_tmp, block_G_full)
     310              :          END IF
     311              : 
     312              :          DEALLOCATE (block_X_full, evals, evecs)
     313              : 
     314              :       END DO
     315              :       CALL dbcsr_iterator_stop(iter)
     316              : !$OMP END PARALLEL
     317              : 
     318         2710 :       CALL timestop(handle)
     319         2710 :    END SUBROUTINE pao_calc_U_exp
     320              : 
     321              : ! **************************************************************************************************
     322              : !> \brief Helper routine, for calculating derivatives
     323              : !> \param M ...
     324              : !> \param D ...
     325              : !> \param R ...
     326              : !> \return ...
     327              : ! **************************************************************************************************
     328          683 :    FUNCTION fold_derivatives(M, D, R) RESULT(G)
     329              :       REAL(dp), DIMENSION(:, :), INTENT(IN)              :: M
     330              :       COMPLEX(dp), DIMENSION(:, :), INTENT(IN)           :: D, R
     331              :       REAL(dp), DIMENSION(SIZE(M, 1), SIZE(M, 1))        :: G
     332              : 
     333          683 :       COMPLEX(dp), DIMENSION(:, :), POINTER              :: F, RF, RM, RMR
     334              :       INTEGER                                            :: n
     335          683 :       REAL(dp), DIMENSION(:, :), POINTER                 :: RFR
     336              : 
     337          683 :       n = SIZE(M, 1)
     338              : 
     339         8879 :       ALLOCATE (RM(n, n), RMR(n, n), F(n, n), RF(n, n), RFR(n, n))
     340              : 
     341       642537 :       RM = MATMUL(TRANSPOSE(CONJG(R)), TRANSPOSE(M))
     342      1133318 :       RMR = MATMUL(RM, R)
     343       101626 :       F = RMR*D !Hadamard product
     344       641854 :       RF = MATMUL(R, F)
     345      1133318 :       RFR = REAL(MATMUL(RF, TRANSPOSE(CONJG(R))), dp)
     346              : 
     347              :       ! gradient dE/dX has to be anti-symmetric
     348        50813 :       G = 0.5_dp*(TRANSPOSE(RFR) - RFR)
     349              : 
     350          683 :       DEALLOCATE (RM, RMR, F, RF, RFR)
     351          683 :    END FUNCTION fold_derivatives
     352              : 
     353          683 : END MODULE pao_param_exp
        

Generated by: LCOV version 2.0-1