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             : MODULE qs_tddfpt2_assign
       9             :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply
      10             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_geadd,&
      11             :                                               cp_fm_scale
      12             :    USE cp_fm_diag,                      ONLY: cp_fm_power
      13             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      14             :                                               cp_fm_struct_release,&
      15             :                                               cp_fm_struct_type
      16             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      17             :                                               cp_fm_get_diag,&
      18             :                                               cp_fm_get_info,&
      19             :                                               cp_fm_release,&
      20             :                                               cp_fm_to_fm,&
      21             :                                               cp_fm_type
      22             :    USE dbcsr_api,                       ONLY: dbcsr_p_type,&
      23             :                                               dbcsr_type
      24             :    USE exstates_types,                  ONLY: wfn_history_type
      25             :    USE kinds,                           ONLY: dp
      26             :    USE message_passing,                 ONLY: mp_para_env_type
      27             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      28             :    USE qs_environment_types,            ONLY: get_qs_env,&
      29             :                                               qs_environment_type
      30             : #include "./base/base_uses.f90"
      31             : 
      32             :    IMPLICIT NONE
      33             : 
      34             :    PRIVATE
      35             : 
      36             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_assign'
      37             : 
      38             :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
      39             : 
      40             :    PUBLIC :: assign_state
      41             : 
      42             : ! **************************************************************************************************
      43             : 
      44             : CONTAINS
      45             : 
      46             : ! **************************************************************************************************
      47             : !> \brief ...
      48             : !> \param qs_env ...
      49             : !> \param matrix_s ...
      50             : !> \param evects ...
      51             : !> \param psi0 ...
      52             : !> \param wfn_history ...
      53             : !> \param my_state ...
      54             : ! **************************************************************************************************
      55           8 :    SUBROUTINE assign_state(qs_env, matrix_s, evects, psi0, wfn_history, my_state)
      56             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      57             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
      58             :       TYPE(cp_fm_type), DIMENSION(:, :)                  :: evects
      59             :       TYPE(cp_fm_type), DIMENSION(:)                     :: psi0
      60             :       TYPE(wfn_history_type)                             :: wfn_history
      61             :       INTEGER, INTENT(INOUT)                             :: my_state
      62             : 
      63             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'assign_state'
      64             : 
      65             :       INTEGER                                            :: handle, is, ispin, natom, ncol, nspins, &
      66             :                                                             nstate
      67             :       REAL(KIND=dp)                                      :: xsum
      68           8 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: dv, rdiag
      69             :       TYPE(dbcsr_type), POINTER                          :: smat
      70             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      71             : 
      72           8 :       CALL timeset(routineN, handle)
      73             : 
      74           8 :       CALL get_qs_env(qs_env, natom=natom, para_env=para_env)
      75           8 :       nspins = SIZE(psi0)
      76           8 :       nstate = SIZE(evects, 2)
      77             :       !
      78           8 :       smat => matrix_s(1)%matrix
      79             :       !
      80           8 :       IF (ASSOCIATED(wfn_history%evect)) THEN
      81          18 :          ALLOCATE (dv(nstate))
      82             :          !
      83           6 :          wfn_history%gsval = 0.0_dp
      84           6 :          wfn_history%gsmin = 1.0_dp
      85          12 :          DO ispin = 1, nspins
      86           6 :             CALL cp_fm_get_info(wfn_history%evect(ispin), ncol_global=ncol)
      87             :             CALL lowdin_orthogonalization(wfn_history%cpmos(ispin), wfn_history%evect(ispin), &
      88           6 :                                           ncol, smat)
      89          18 :             ALLOCATE (rdiag(ncol))
      90             :             CALL wfn_align(psi0(ispin), wfn_history%cpmos(ispin), wfn_history%evect(ispin), &
      91           6 :                            rdiag, smat)
      92          30 :             wfn_history%gsval = wfn_history%gsval + SUM(rdiag)/REAL(ncol*nspins, KIND=dp)
      93          36 :             wfn_history%gsmin = MIN(wfn_history%gsmin, MINVAL(rdiag))
      94          18 :             DEALLOCATE (rdiag)
      95             :          END DO
      96          36 :          DO is = 1, nstate
      97             :             xsum = 0.0_dp
      98          60 :             DO ispin = 1, nspins
      99          30 :                CALL cp_fm_get_info(wfn_history%evect(ispin), ncol_global=ncol)
     100          90 :                ALLOCATE (rdiag(ncol))
     101          30 :                CALL xvec_ovlp(evects(ispin, is), wfn_history%evect(ispin), rdiag, smat)
     102         150 :                xsum = xsum + SUM(rdiag)
     103          90 :                DEALLOCATE (rdiag)
     104             :             END DO
     105          36 :             dv(is) = ABS(xsum)/SQRT(REAL(nspins, dp))
     106             :          END DO
     107          42 :          my_state = MAXVAL(MAXLOC(dv))
     108           6 :          wfn_history%xsval = dv(my_state)
     109           6 :          IF (wfn_history%xsval < 0.75_dp) THEN
     110           0 :             dv(my_state) = 0.0_dp
     111           0 :             IF (wfn_history%xsval/MAXVAL(dv) < 0.5_dp) THEN
     112             :                CALL cp_warn(__LOCATION__, "Uncertain assignment for State following."// &
     113             :                             " Reduce trust radius in Geometry Optimization or timestep"// &
     114           0 :                             " in MD runs.")
     115             :             END IF
     116             :          END IF
     117          12 :          DO ispin = 1, nspins
     118           6 :             CALL cp_fm_get_info(wfn_history%evect(ispin), ncol_global=ncol)
     119           6 :             CALL cp_fm_to_fm(evects(ispin, my_state), wfn_history%evect(ispin))
     120          18 :             CALL cp_fm_to_fm(psi0(ispin), wfn_history%cpmos(ispin), ncol, 1, 1)
     121             :          END DO
     122             :          !
     123           6 :          DEALLOCATE (dv)
     124             :       ELSE
     125             :          !
     126           8 :          ALLOCATE (wfn_history%evect(nspins))
     127           6 :          ALLOCATE (wfn_history%cpmos(nspins))
     128           4 :          DO ispin = 1, nspins
     129           2 :             CALL cp_fm_create(wfn_history%evect(ispin), evects(ispin, 1)%matrix_struct, "Xvec")
     130           4 :             CALL cp_fm_create(wfn_history%cpmos(ispin), evects(ispin, 1)%matrix_struct, "Cvec")
     131             :          END DO
     132           4 :          DO ispin = 1, nspins
     133           2 :             CALL cp_fm_get_info(wfn_history%evect(ispin), ncol_global=ncol)
     134           2 :             CALL cp_fm_to_fm(evects(ispin, my_state), wfn_history%evect(ispin))
     135           6 :             CALL cp_fm_to_fm(psi0(ispin), wfn_history%cpmos(ispin), ncol, 1, 1)
     136             :          END DO
     137           2 :          wfn_history%xsval = 1.0_dp
     138           2 :          wfn_history%gsval = 1.0_dp
     139           2 :          wfn_history%gsmin = 1.0_dp
     140             :       END IF
     141             : 
     142           8 :       CALL timestop(handle)
     143             : 
     144           8 :    END SUBROUTINE assign_state
     145             : 
     146             : ! **************************************************************************************************
     147             : !> \brief return a set of S orthonormal vectors (C^T S C == 1) where
     148             : !>      a Lodwin transformation is applied to keep the rotated vectors as close
     149             : !>      as possible to the original ones
     150             : !> \param vmatrix ...
     151             : !> \param xmatrix ...
     152             : !> \param ncol ...
     153             : !> \param matrix_s ...
     154             : !> \param
     155             : !> \par History
     156             : !>      05.2009 created [MI]
     157             : !>      06.2023 adapted to include a second set of vectors [JGH]
     158             : !> \note
     159             : ! **************************************************************************************************
     160          12 :    SUBROUTINE lowdin_orthogonalization(vmatrix, xmatrix, ncol, matrix_s)
     161             : 
     162             :       TYPE(cp_fm_type), INTENT(IN)                       :: vmatrix, xmatrix
     163             :       INTEGER, INTENT(IN)                                :: ncol
     164             :       TYPE(dbcsr_type)                                   :: matrix_s
     165             : 
     166             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'lowdin_orthogonalization'
     167             :       REAL(KIND=dp), PARAMETER                           :: rone = 1.0_dp, rzero = 0.0_dp
     168             : 
     169             :       INTEGER                                            :: handle, n, ncol_global, ndep
     170             :       REAL(dp)                                           :: threshold, xsum
     171          12 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: rdiag
     172             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     173             :       TYPE(cp_fm_type)                                   :: csc, sc, work
     174             : 
     175          12 :       IF (ncol .EQ. 0) RETURN
     176             : 
     177          12 :       CALL timeset(routineN, handle)
     178             : 
     179          12 :       threshold = 1.0E-7_dp
     180          12 :       CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
     181          12 :       IF (ncol .GT. ncol_global) CPABORT("Wrong ncol value")
     182             : 
     183          12 :       CALL cp_fm_create(sc, xmatrix%matrix_struct, "SC")
     184          12 :       CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, sc, ncol)
     185             : 
     186          12 :       NULLIFY (fm_struct_tmp)
     187             :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
     188             :                                para_env=vmatrix%matrix_struct%para_env, &
     189          12 :                                context=vmatrix%matrix_struct%context)
     190          12 :       CALL cp_fm_create(csc, fm_struct_tmp, "csc")
     191          12 :       CALL cp_fm_create(work, fm_struct_tmp, "work")
     192          12 :       CALL cp_fm_struct_release(fm_struct_tmp)
     193             : 
     194          12 :       CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, sc, rzero, csc)
     195          12 :       CALL cp_fm_power(csc, work, -0.5_dp, threshold, ndep)
     196          12 :       CALL parallel_gemm('N', 'N', n, ncol, ncol, rone, vmatrix, csc, rzero, sc)
     197          12 :       CALL cp_fm_to_fm(sc, vmatrix, ncol, 1, 1)
     198             :       !
     199          12 :       CALL parallel_gemm('N', 'N', n, ncol, ncol, rone, xmatrix, csc, rzero, sc)
     200          12 :       CALL cp_fm_to_fm(sc, xmatrix, ncol, 1, 1)
     201             :       ! projecton for xSv = 0
     202          12 :       CALL cp_dbcsr_sm_fm_multiply(matrix_s, xmatrix, sc, ncol)
     203          12 :       CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, vmatrix, sc, rzero, csc)
     204          12 :       CALL parallel_gemm('N', 'N', n, ncol, ncol, rone, vmatrix, csc, rzero, sc)
     205          12 :       CALL cp_fm_geadd(-1.0_dp, 'N', sc, 1.0_dp, xmatrix)
     206             :       ! normalisation
     207          12 :       CALL cp_dbcsr_sm_fm_multiply(matrix_s, xmatrix, sc, ncol)
     208          12 :       CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, xmatrix, sc, rzero, csc)
     209          36 :       ALLOCATE (rdiag(ncol))
     210          12 :       CALL cp_fm_get_diag(csc, rdiag)
     211          60 :       xsum = SUM(rdiag)
     212          12 :       DEALLOCATE (rdiag)
     213          12 :       xsum = 1._dp/SQRT(xsum)
     214          12 :       CALL cp_fm_scale(xsum, xmatrix)
     215             : 
     216          12 :       CALL cp_fm_release(csc)
     217          12 :       CALL cp_fm_release(sc)
     218          12 :       CALL cp_fm_release(work)
     219             : 
     220          12 :       CALL timestop(handle)
     221             : 
     222          36 :    END SUBROUTINE lowdin_orthogonalization
     223             : 
     224             : ! **************************************************************************************************
     225             : !> \brief ...
     226             : !> \param gmatrix ...
     227             : !> \param vmatrix ...
     228             : !> \param xmatrix ...
     229             : !> \param rdiag ...
     230             : !> \param matrix_s ...
     231             : ! **************************************************************************************************
     232           6 :    SUBROUTINE wfn_align(gmatrix, vmatrix, xmatrix, rdiag, matrix_s)
     233             : 
     234             :       TYPE(cp_fm_type), INTENT(IN)                       :: gmatrix, vmatrix, xmatrix
     235             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: rdiag
     236             :       TYPE(dbcsr_type)                                   :: matrix_s
     237             : 
     238             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'wfn_align'
     239             :       REAL(KIND=dp), PARAMETER                           :: rone = 1.0_dp, rzero = 0.0_dp
     240             : 
     241             :       INTEGER                                            :: handle, n, ncol, ncol_global
     242             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     243             :       TYPE(cp_fm_type)                                   :: csc, sc
     244             : 
     245           6 :       CALL timeset(routineN, handle)
     246             : 
     247           6 :       ncol = SIZE(rdiag)
     248           6 :       CALL cp_fm_get_info(matrix=vmatrix, nrow_global=n, ncol_global=ncol_global)
     249           6 :       IF (ncol .GT. ncol_global) CPABORT("Wrong ncol value")
     250             : 
     251           6 :       CALL cp_fm_create(sc, vmatrix%matrix_struct, "SC")
     252           6 :       CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, sc, ncol)
     253             : 
     254           6 :       NULLIFY (fm_struct_tmp)
     255             :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
     256             :                                para_env=vmatrix%matrix_struct%para_env, &
     257           6 :                                context=vmatrix%matrix_struct%context)
     258           6 :       CALL cp_fm_create(csc, fm_struct_tmp, "csc")
     259           6 :       CALL cp_fm_struct_release(fm_struct_tmp)
     260             : 
     261           6 :       CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, gmatrix, sc, rzero, csc)
     262           6 :       CALL parallel_gemm('N', 'T', n, ncol, ncol, rone, vmatrix, csc, rzero, sc)
     263           6 :       CALL cp_fm_to_fm(sc, vmatrix, ncol, 1, 1)
     264           6 :       CALL parallel_gemm('N', 'T', n, ncol, ncol, rone, xmatrix, csc, rzero, sc)
     265           6 :       CALL cp_fm_to_fm(sc, xmatrix, ncol, 1, 1)
     266             :       !
     267           6 :       CALL lowdin_orthogonalization(vmatrix, xmatrix, ncol, matrix_s)
     268             :       !
     269           6 :       CALL cp_dbcsr_sm_fm_multiply(matrix_s, vmatrix, sc, ncol)
     270           6 :       CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, gmatrix, sc, rzero, csc)
     271           6 :       CALL cp_fm_get_diag(csc, rdiag)
     272             : 
     273           6 :       CALL cp_fm_release(csc)
     274           6 :       CALL cp_fm_release(sc)
     275             : 
     276           6 :       CALL timestop(handle)
     277             : 
     278           6 :    END SUBROUTINE wfn_align
     279             : 
     280             : ! **************************************************************************************************
     281             : !> \brief ...
     282             : !> \param ematrix ...
     283             : !> \param xmatrix ...
     284             : !> \param rdiag ...
     285             : !> \param matrix_s ...
     286             : ! **************************************************************************************************
     287          30 :    SUBROUTINE xvec_ovlp(ematrix, xmatrix, rdiag, matrix_s)
     288             : 
     289             :       TYPE(cp_fm_type), INTENT(IN)                       :: ematrix, xmatrix
     290             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: rdiag
     291             :       TYPE(dbcsr_type)                                   :: matrix_s
     292             : 
     293             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'xvec_ovlp'
     294             :       REAL(KIND=dp), PARAMETER                           :: rone = 1.0_dp, rzero = 0.0_dp
     295             : 
     296             :       INTEGER                                            :: handle, n, ncol, ncol_global
     297             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     298             :       TYPE(cp_fm_type)                                   :: csc, sc
     299             : 
     300          30 :       CALL timeset(routineN, handle)
     301             : 
     302          30 :       ncol = SIZE(rdiag)
     303          30 :       CALL cp_fm_get_info(matrix=xmatrix, nrow_global=n, ncol_global=ncol_global)
     304          30 :       IF (ncol .GT. ncol_global) CPABORT("Wrong ncol value")
     305             : 
     306          30 :       CALL cp_fm_create(sc, xmatrix%matrix_struct, "SC")
     307          30 :       CALL cp_dbcsr_sm_fm_multiply(matrix_s, xmatrix, sc, ncol)
     308             : 
     309          30 :       NULLIFY (fm_struct_tmp)
     310             :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol, ncol_global=ncol, &
     311             :                                para_env=xmatrix%matrix_struct%para_env, &
     312          30 :                                context=xmatrix%matrix_struct%context)
     313          30 :       CALL cp_fm_create(csc, fm_struct_tmp, "csc")
     314          30 :       CALL cp_fm_struct_release(fm_struct_tmp)
     315             : 
     316          30 :       CALL parallel_gemm('T', 'N', ncol, ncol, n, rone, ematrix, sc, rzero, csc)
     317          30 :       CALL cp_fm_get_diag(csc, rdiag)
     318             : 
     319          30 :       CALL cp_fm_release(csc)
     320          30 :       CALL cp_fm_release(sc)
     321             : 
     322          30 :       CALL timestop(handle)
     323             : 
     324          30 :    END SUBROUTINE xvec_ovlp
     325             : 
     326             : END MODULE qs_tddfpt2_assign