LCOV - code coverage report
Current view: top level - src - et_coupling.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 0.0 % 102 0
Test Date: 2025-07-25 12:55:17 Functions: 0.0 % 1 0

            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 calculates the electron transfer coupling elements
      10              : !>      Wu, Van Voorhis, JCP 125, 164105 (2006)
      11              : !> \author fschiff (01.2007)
      12              : ! **************************************************************************************************
      13              : MODULE et_coupling
      14              :    USE cp_control_types,                ONLY: dft_control_type
      15              :    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type
      16              :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply,&
      17              :                                               dbcsr_deallocate_matrix_set
      18              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_det,&
      19              :                                               cp_fm_invert,&
      20              :                                               cp_fm_transpose
      21              :    USE cp_fm_pool_types,                ONLY: cp_fm_pool_p_type,&
      22              :                                               fm_pool_get_el_struct
      23              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_type
      24              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      25              :                                               cp_fm_get_info,&
      26              :                                               cp_fm_release,&
      27              :                                               cp_fm_type
      28              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      29              :                                               cp_logger_type,&
      30              :                                               cp_to_string
      31              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      32              :                                               cp_print_key_unit_nr
      33              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      34              :                                               section_vals_type
      35              :    USE kinds,                           ONLY: dp
      36              :    USE mathlib,                         ONLY: diamat_all
      37              :    USE message_passing,                 ONLY: mp_para_env_type
      38              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      39              :    USE qs_energy_types,                 ONLY: qs_energy_type
      40              :    USE qs_environment_types,            ONLY: get_qs_env,&
      41              :                                               qs_environment_type
      42              :    USE qs_matrix_pools,                 ONLY: mpools_get
      43              :    USE qs_mo_types,                     ONLY: get_mo_set
      44              : #include "./base/base_uses.f90"
      45              : 
      46              :    IMPLICIT NONE
      47              : 
      48              :    PRIVATE
      49              : 
      50              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'et_coupling'
      51              :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
      52              : 
      53              : ! *** Public subroutines ***
      54              : 
      55              :    PUBLIC :: calc_et_coupling
      56              : 
      57              : CONTAINS
      58              : ! **************************************************************************************************
      59              : !> \brief ...
      60              : !> \param qs_env ...
      61              : ! **************************************************************************************************
      62            0 :    SUBROUTINE calc_et_coupling(qs_env)
      63              : 
      64              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      65              : 
      66              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'calc_et_coupling'
      67              : 
      68              :       INTEGER                                            :: handle, i, iw, j, k, nao, ncol_local, &
      69              :                                                             nmo, nrow_local
      70            0 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
      71              :       LOGICAL                                            :: is_spin_constraint
      72              :       REAL(KIND=dp)                                      :: Sda, strength, Waa, Wbb, Wda
      73            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: a, b, S_det
      74              :       REAL(KIND=dp), DIMENSION(2)                        :: eigenv
      75              :       REAL(KIND=dp), DIMENSION(2, 2)                     :: S_mat, tmp_mat, U, W_mat
      76            0 :       TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER     :: mo_mo_fm_pools
      77              :       TYPE(cp_fm_struct_type), POINTER                   :: mo_mo_fmstruct
      78              :       TYPE(cp_fm_type)                                   :: inverse_mat, SMO, Tinverse, tmp2
      79            0 :       TYPE(cp_fm_type), DIMENSION(2)                     :: rest_MO
      80              :       TYPE(cp_logger_type), POINTER                      :: logger
      81            0 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
      82              :       TYPE(dft_control_type), POINTER                    :: dft_control
      83              :       TYPE(mp_para_env_type), POINTER                    :: para_env
      84              :       TYPE(qs_energy_type), POINTER                      :: energy
      85              :       TYPE(section_vals_type), POINTER                   :: et_coupling_section
      86              : 
      87            0 :       NULLIFY (mo_mo_fmstruct, energy, matrix_s, dft_control, para_env)
      88              : 
      89            0 :       CALL timeset(routineN, handle)
      90              : 
      91            0 :       logger => cp_get_default_logger()
      92              :       et_coupling_section => section_vals_get_subs_vals(qs_env%input, &
      93            0 :                                                         "PROPERTIES%ET_COUPLING")
      94              : 
      95            0 :       CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env)
      96              : 
      97              :       iw = cp_print_key_unit_nr(logger, et_coupling_section, "PROGRAM_RUN_INFO", &
      98            0 :                                 extension=".log")
      99              : 
     100            0 :       is_spin_constraint = .FALSE.
     101            0 :       ALLOCATE (a(dft_control%nspins))
     102            0 :       ALLOCATE (b(dft_control%nspins))
     103            0 :       ALLOCATE (S_det(dft_control%nspins))
     104              : 
     105            0 :       CALL mpools_get(qs_env%mpools, mo_mo_fm_pools=mo_mo_fm_pools)
     106            0 :       mo_mo_fmstruct => fm_pool_get_el_struct(mo_mo_fm_pools(1)%pool)
     107            0 :       DO i = 1, dft_control%nspins
     108            0 :          mo_mo_fmstruct => fm_pool_get_el_struct(mo_mo_fm_pools(i)%pool)
     109              : 
     110              :          CALL get_mo_set(mo_set=qs_env%mos(i), &
     111              :                          nao=nao, &
     112            0 :                          nmo=nmo)
     113              : 
     114              :          CALL cp_fm_create(matrix=tmp2, &
     115              :                            matrix_struct=qs_env%mos(i)%mo_coeff%matrix_struct, &
     116            0 :                            name="ET_TMP"//TRIM(ADJUSTL(cp_to_string(2)))//"MATRIX")
     117              :          CALL cp_fm_create(matrix=inverse_mat, &
     118              :                            matrix_struct=mo_mo_fmstruct, &
     119            0 :                            name="INVERSE"//TRIM(ADJUSTL(cp_to_string(2)))//"MATRIX")
     120              :          CALL cp_fm_create(matrix=Tinverse, &
     121              :                            matrix_struct=mo_mo_fmstruct, &
     122            0 :                            name="T_INVERSE"//TRIM(ADJUSTL(cp_to_string(2)))//"MATRIX")
     123              :          CALL cp_fm_create(matrix=SMO, &
     124              :                            matrix_struct=mo_mo_fmstruct, &
     125            0 :                            name="ET_SMO"//TRIM(ADJUSTL(cp_to_string(1)))//"MATRIX")
     126            0 :          DO j = 1, 2
     127              :             CALL cp_fm_create(matrix=rest_MO(j), &
     128              :                               matrix_struct=mo_mo_fmstruct, &
     129            0 :                               name="ET_rest_MO"//TRIM(ADJUSTL(cp_to_string(j)))//"MATRIX")
     130              :          END DO
     131              : 
     132              : !   calculate MO-overlap
     133              : 
     134            0 :          CALL get_qs_env(qs_env, matrix_s=matrix_s)
     135              :          CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, qs_env%et_coupling%et_mo_coeff(i), &
     136            0 :                                       tmp2, nmo, 1.0_dp, 0.0_dp)
     137              :          CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
     138              :                             qs_env%mos(i)%mo_coeff, &
     139            0 :                             tmp2, 0.0_dp, SMO)
     140              : 
     141              : !    calculate the MO-representation of the restraint matrix A
     142              : 
     143              :          CALL cp_dbcsr_sm_fm_multiply(qs_env%et_coupling%rest_mat(1)%matrix, &
     144              :                                       qs_env%et_coupling%et_mo_coeff(i), &
     145            0 :                                       tmp2, nmo, 1.0_dp, 0.0_dp)
     146              : 
     147              :          CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
     148              :                             qs_env%mos(i)%mo_coeff, &
     149            0 :                             tmp2, 0.0_dp, rest_MO(1))
     150              : 
     151              : !    calculate the MO-representation of the restraint matrix D
     152              : 
     153              :          CALL cp_dbcsr_sm_fm_multiply(qs_env%et_coupling%rest_mat(2)%matrix, &
     154              :                                       qs_env%mos(i)%mo_coeff, &
     155            0 :                                       tmp2, nmo, 1.0_dp, 0.0_dp)
     156              : 
     157              :          CALL parallel_gemm('T', 'N', nmo, nmo, nao, 1.0_dp, &
     158              :                             qs_env%et_coupling%et_mo_coeff(i), &
     159            0 :                             tmp2, 0.0_dp, rest_MO(2))
     160              : ! TODO: could fix cp_fm_invert/LU determinant pivoting instead of calling cp_fm_det to save a pdgetrf call
     161            0 :          CALL cp_fm_det(SMO, S_det(i))
     162            0 :          CALL cp_fm_invert(SMO, inverse_mat)
     163              : 
     164              :          CALL cp_fm_get_info(inverse_mat, nrow_local=nrow_local, ncol_local=ncol_local, &
     165            0 :                              row_indices=row_indices, col_indices=col_indices)
     166            0 :          b(i) = 0.0_dp
     167              : 
     168            0 :          DO j = 1, ncol_local
     169            0 :             DO k = 1, nrow_local
     170            0 :                b(i) = b(i) + rest_MO(2)%local_data(k, j)*inverse_mat%local_data(k, j)
     171              :             END DO
     172              :          END DO
     173              : 
     174            0 :          CALL cp_fm_transpose(inverse_mat, Tinverse)
     175            0 :          a(i) = 0.0_dp
     176            0 :          DO j = 1, ncol_local
     177            0 :             DO k = 1, nrow_local
     178            0 :                a(i) = a(i) + rest_MO(1)%local_data(k, j)*Tinverse%local_data(k, j)
     179              :             END DO
     180              :          END DO
     181              :          IF (is_spin_constraint) THEN
     182              :             a(i) = -a(i)
     183              :             b(i) = -b(i)
     184              :          END IF
     185            0 :          CALL para_env%sum(a(i))
     186              : 
     187            0 :          CALL para_env%sum(b(i))
     188              : 
     189            0 :          CALL cp_fm_release(tmp2)
     190            0 :          DO j = 1, 2
     191            0 :             CALL cp_fm_release(rest_MO(j))
     192              :          END DO
     193            0 :          CALL cp_fm_release(SMO)
     194            0 :          CALL cp_fm_release(Tinverse)
     195            0 :          CALL cp_fm_release(inverse_mat)
     196              :       END DO
     197              : 
     198              : !    solve eigenstates for the projector matrix
     199              : 
     200            0 :       IF (dft_control%nspins == 2) THEN
     201            0 :          Sda = S_det(1)*S_det(2)
     202            0 :          Wda = ((a(1) + a(2)) + (b(1) + b(2)))*0.5_dp*Sda
     203              :       ELSE
     204            0 :          Sda = S_det(1)**2
     205            0 :          Wda = (a(1) + b(1))*Sda
     206              :       END IF
     207              : 
     208            0 :       IF (dft_control%qs_control%ddapc_restraint) THEN
     209            0 :          Waa = qs_env%et_coupling%order_p
     210            0 :          Wbb = dft_control%qs_control%ddapc_restraint_control(1)%ddapc_order_p
     211            0 :          strength = dft_control%qs_control%ddapc_restraint_control(1)%strength
     212              :       END IF
     213              : 
     214              : !!   construct S and W   !!!
     215            0 :       S_mat(1, 1) = 1.0_dp
     216            0 :       S_mat(2, 2) = 1.0_dp
     217            0 :       S_mat(2, 1) = Sda
     218            0 :       S_mat(1, 2) = Sda
     219              : 
     220            0 :       W_mat(1, 1) = Wbb
     221            0 :       W_mat(2, 2) = Waa
     222            0 :       W_mat(2, 1) = Wda
     223            0 :       W_mat(1, 2) = Wda
     224              : 
     225              : !!  solve WC=SCN
     226            0 :       CALL diamat_all(S_mat, eigenv, .TRUE.)
     227              :       ! U = S**(-1/2)
     228            0 :       U = 0.0_dp
     229            0 :       U(1, 1) = 1.0_dp/SQRT(eigenv(1))
     230            0 :       U(2, 2) = 1.0_dp/SQRT(eigenv(2))
     231            0 :       tmp_mat = MATMUL(U, TRANSPOSE(S_mat))
     232            0 :       U = MATMUL(S_mat, tmp_mat)
     233            0 :       tmp_mat = MATMUL(W_mat, U)
     234            0 :       W_mat = MATMUL(U, tmp_mat)
     235            0 :       CALL diamat_all(W_mat, eigenv, .TRUE.)
     236            0 :       tmp_mat = MATMUL(U, W_mat)
     237              : 
     238            0 :       CALL get_qs_env(qs_env, energy=energy)
     239            0 :       W_mat(1, 1) = energy%total
     240            0 :       W_mat(2, 2) = qs_env%et_coupling%energy
     241            0 :       a(1) = (energy%total + strength*Wbb)*Sda - strength*Wda
     242            0 :       a(2) = (qs_env%et_coupling%energy + qs_env%et_coupling%e1*Waa)*Sda - qs_env%et_coupling%e1*Wda
     243            0 :       W_mat(1, 2) = (a(1) + a(2))*0.5_dp
     244            0 :       W_mat(2, 1) = W_mat(1, 2)
     245              : 
     246            0 :       S_mat = MATMUL(W_mat, (tmp_mat))
     247            0 :       W_mat = MATMUL(TRANSPOSE(tmp_mat), S_mat)
     248              : 
     249            0 :       IF (iw > 0) THEN
     250            0 :          WRITE (iw, *)
     251            0 :          WRITE (iw, '(T3,A,T60,(3X,F12.6))') 'Strength of constraint A          :', qs_env%et_coupling%e1
     252            0 :          WRITE (iw, '(T3,A,T60,(3X,F12.6))') 'Strength of constraint B          :', strength
     253            0 :          WRITE (iw, '(T3,A,T60,(3X,F12.6))') 'Final target value of constraint A:', Waa
     254            0 :          WRITE (iw, '(T3,A,T60,(3X,F12.6))') 'Final target value of constraint B:', Wbb
     255            0 :          WRITE (iw, *)
     256              :          WRITE (iw, '(T3,A,T60,(3X,F12.6))') &
     257            0 :             'Diabatic electronic coupling matrix element(mHartree):', ABS(W_mat(1, 2)*1000.0_dp)
     258              : 
     259              :       END IF
     260              : 
     261            0 :       CALL dbcsr_deallocate_matrix_set(qs_env%et_coupling%rest_mat)
     262              : 
     263              :       CALL cp_print_key_finished_output(iw, logger, et_coupling_section, &
     264            0 :                                         "PROGRAM_RUN_INFO")
     265            0 :       CALL timestop(handle)
     266            0 :    END SUBROUTINE calc_et_coupling
     267              : 
     268              : END MODULE et_coupling
     269              : 
        

Generated by: LCOV version 2.0-1