LCOV - code coverage report
Current view: top level - src - qs_tddfpt2_assign.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.7 % 130 127
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 4 4

            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              : MODULE qs_tddfpt2_assign
       9              :    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type,&
      10              :                                               dbcsr_type
      11              :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply
      12              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_geadd,&
      13              :                                               cp_fm_scale
      14              :    USE cp_fm_diag,                      ONLY: cp_fm_power
      15              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      16              :                                               cp_fm_struct_release,&
      17              :                                               cp_fm_struct_type
      18              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      19              :                                               cp_fm_get_diag,&
      20              :                                               cp_fm_get_info,&
      21              :                                               cp_fm_release,&
      22              :                                               cp_fm_to_fm,&
      23              :                                               cp_fm_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
        

Generated by: LCOV version 2.0-1