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             : ! **************************************************************************************************
       9             : !> \brief Methods related to (\cal S)^2 (i.e. spin)
      10             : !> \par History
      11             : !>      03.2006 copied compute_s_square from qs_scf_post (Joost VandeVondele)
      12             : !>      08.2021 revised (Matthias Krack, MK)
      13             : !> \author Joost VandeVondele
      14             : ! **************************************************************************************************
      15             : MODULE s_square_methods
      16             : 
      17             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      18             :    USE cp_control_types,                ONLY: s2_restraint_type
      19             :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply
      20             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      21             :                                               cp_fm_struct_release,&
      22             :                                               cp_fm_struct_type
      23             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      24             :                                               cp_fm_get_info,&
      25             :                                               cp_fm_release,&
      26             :                                               cp_fm_type
      27             :    USE dbcsr_api,                       ONLY: dbcsr_p_type
      28             :    USE input_constants,                 ONLY: do_s2_constraint,&
      29             :                                               do_s2_restraint
      30             :    USE kinds,                           ONLY: dp
      31             :    USE message_passing,                 ONLY: mp_para_env_type
      32             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      33             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      34             :                                               has_uniform_occupation,&
      35             :                                               mo_set_type
      36             : #include "./base/base_uses.f90"
      37             : 
      38             :    IMPLICIT NONE
      39             : 
      40             :    PRIVATE
      41             : 
      42             :    ! Global parameters
      43             : 
      44             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 's_square_methods'
      45             : 
      46             :    PUBLIC :: compute_s_square, s2_restraint
      47             : 
      48             : CONTAINS
      49             : 
      50             : ! **************************************************************************************************
      51             : !> \brief Compute the expectation value <(\cal S)^2> of the single determinant defined by the
      52             : !>        spin up (alpha) and spin down (beta) orbitals
      53             : !> \param mos [in] MO set with all MO information including the alpha and beta MO coefficients
      54             : !> \param matrix_s [in] AO overlap matrix S (do not mix with the spin operator (\cal S))
      55             : !> \param s_square [out] <(\cal S)^2> including potential spin contaminations
      56             : !> \param s_square_ideal [out] Ideal value for <(\cal S)^2> without any spin contaminations
      57             : !> \param mo_derivs [inout] If present, add the derivative of s_square wrt the MOs to mo_derivs
      58             : !> \param strength [in] Strength for constraining or restraining (\cal S)^2
      59             : !> \par History
      60             : !>      07.2004 created (Joost VandeVondele)
      61             : !>      08.2021 revised (Matthias Krack, MK)
      62             : !> \note
      63             : !>      See Eqs. 2.271 and 2.272 in Modern Quantum Chemistry by A. Szabo and N. S. Ostlund
      64             : ! **************************************************************************************************
      65        1810 :    SUBROUTINE compute_s_square(mos, matrix_s, s_square, s_square_ideal, mo_derivs, strength)
      66             : 
      67             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
      68             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
      69             :       REAL(KIND=dp), INTENT(OUT)                         :: s_square, s_square_ideal
      70             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT), &
      71             :          OPTIONAL                                        :: mo_derivs
      72             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: strength
      73             : 
      74             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'compute_s_square'
      75             : 
      76             :       INTEGER                                            :: handle, i, j, nalpha, nao, nao_beta, &
      77             :                                                             nbeta, ncol_local, nmo_alpha, &
      78             :                                                             nmo_beta, nrow_local
      79             :       LOGICAL                                            :: has_uniform_occupation_alpha, &
      80             :                                                             has_uniform_occupation_beta
      81             :       REAL(KIND=dp)                                      :: s2
      82             :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
      83        1810 :          POINTER                                         :: local_data
      84             :       TYPE(cp_blacs_env_type), POINTER                   :: context
      85             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
      86             :       TYPE(cp_fm_type)                                   :: catscb, sca, scb
      87             :       TYPE(cp_fm_type), POINTER                          :: c_alpha, c_beta
      88             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      89             : 
      90        1810 :       CALL timeset(routineN, handle)
      91             : 
      92        1810 :       NULLIFY (context)
      93        1810 :       NULLIFY (fm_struct_tmp)
      94        1810 :       NULLIFY (local_data)
      95        1810 :       NULLIFY (para_env)
      96             : 
      97        1810 :       SELECT CASE (SIZE(mos))
      98             :       CASE (1)
      99             :          ! Spin restricted case, i.e. nothing to do
     100           0 :          s_square = 0.0_dp
     101           0 :          s_square_ideal = 0.0_dp
     102             :          ! Restraining or constraining (\cal S) does not make sense
     103           0 :          CPASSERT(PRESENT(mo_derivs) .OR. PRESENT(strength))
     104             :       CASE (2)
     105        1810 :          CALL get_mo_set(mo_set=mos(1), mo_coeff=c_alpha, homo=nalpha, nmo=nmo_alpha, nao=nao)
     106        1810 :          CALL get_mo_set(mo_set=mos(2), mo_coeff=c_beta, homo=nbeta, nmo=nmo_beta, nao=nao_beta)
     107        1810 :          CPASSERT(nao == nao_beta)
     108        1810 :          has_uniform_occupation_alpha = has_uniform_occupation(mo_set=mos(1), last_mo=nalpha)
     109        1810 :          has_uniform_occupation_beta = has_uniform_occupation(mo_set=mos(2), last_mo=nbeta)
     110             :          ! This makes only sense if we have uniform occupations for the alpha and beta spin MOs while
     111             :          ! ignoring potentially added MOs with zero occupation
     112        1810 :          IF (has_uniform_occupation_alpha .AND. has_uniform_occupation_beta) THEN
     113             :             ! Eq. 2.272 in Modern Quantum Chemistry by A. Szabo and N. S. Ostlund
     114        1810 :             s_square_ideal = REAL((nalpha - nbeta)*(nalpha - nbeta + 2), KIND=dp)/4.0_dp
     115             :             ! Create overlap matrix
     116        1810 :             CALL cp_fm_get_info(c_alpha, para_env=para_env, context=context)
     117             :             CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
     118        1810 :                                      nrow_global=nalpha, ncol_global=nbeta)
     119             :             ! Prepare C(alpha)^T*S*C(beta)
     120        1810 :             CALL cp_fm_create(catscb, fm_struct_tmp, name="C(alpha)^T*S*C(beta)")
     121        1810 :             CALL cp_fm_struct_release(fm_struct_tmp)
     122             :             ! Create S*C(beta)
     123        1810 :             CALL cp_fm_get_info(c_beta, matrix_struct=fm_struct_tmp)
     124        1810 :             CALL cp_fm_create(scb, fm_struct_tmp, name="S*C(beta)")
     125             :             ! Compute S*C(beta)
     126        1810 :             CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, c_beta, scb, nbeta)
     127             :             ! Compute C(alpha)^T*S*C(beta)
     128        1810 :             CALL parallel_gemm('T', 'N', nalpha, nbeta, nao, 1.0_dp, c_alpha, scb, 0.0_dp, catscb)
     129             :             ! Eq. 2.271 in Modern Quantum Chemistry by A. Szabo and N. S. Ostlund
     130        1810 :             CALL cp_fm_get_info(catscb, local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
     131        1810 :             s2 = 0.0_dp
     132        8332 :             DO j = 1, ncol_local
     133       51069 :                DO i = 1, nrow_local
     134       49259 :                   s2 = s2 + local_data(i, j)**2
     135             :                END DO
     136             :             END DO
     137        1810 :             CALL para_env%sum(s2)
     138        1810 :             s_square = s_square_ideal + nbeta - s2
     139             :             ! Only needed for constraining or restraining (\cal S)
     140        1810 :             IF (PRESENT(mo_derivs)) THEN
     141           0 :                CPASSERT(SIZE(mo_derivs, 1) == 2)
     142             :                ! This gets really wrong for fractional occupations
     143           0 :                CALL get_mo_set(mo_set=mos(1), uniform_occupation=has_uniform_occupation_alpha)
     144           0 :                CPASSERT(has_uniform_occupation_alpha)
     145           0 :                CALL get_mo_set(mo_set=mos(2), uniform_occupation=has_uniform_occupation_beta)
     146           0 :                CPASSERT(has_uniform_occupation_beta)
     147             :                ! Add -strength*S*C(beta)*(C(alpha)^T*S*C(beta))^T to the alpha MO derivatives
     148           0 :                CALL parallel_gemm('N', 'T', nao, nalpha, nbeta, -strength, scb, catscb, 1.0_dp, mo_derivs(1))
     149             :                ! Create S*C(alpha)
     150           0 :                CALL cp_fm_get_info(c_alpha, matrix_struct=fm_struct_tmp)
     151           0 :                CALL cp_fm_create(sca, fm_struct_tmp, name="S*C(alpha)")
     152             :                ! Compute S*C(alpha)
     153           0 :                CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, c_alpha, sca, nalpha)
     154             :                ! Add -strength*S*C(alpha)*C(alpha)^T*S*C(beta) to the beta MO derivatives
     155           0 :                CALL parallel_gemm('N', 'N', nao, nbeta, nalpha, -strength, sca, catscb, 1.0_dp, mo_derivs(2))
     156           0 :                CALL cp_fm_release(sca)
     157             :             END IF
     158        1810 :             CALL cp_fm_release(scb)
     159        1810 :             CALL cp_fm_release(catscb)
     160             :          ELSE
     161           0 :             IF (.NOT. has_uniform_occupation_alpha) THEN
     162           0 :                CPHINT("The alpha orbitals (MOs) have a non-uniform occupation")
     163             :             END IF
     164           0 :             IF (.NOT. has_uniform_occupation_beta) THEN
     165           0 :                CPHINT("The beta orbitals (MOs) have a non-uniform occupation")
     166             :             END IF
     167           0 :             CPHINT("Skipping S**2 calculation")
     168             :          END IF
     169             :       CASE DEFAULT
     170             :          ! We should never get here
     171        1810 :          CPABORT("Alpha, beta, what else ?")
     172             :       END SELECT
     173             : 
     174        1810 :       CALL timestop(handle)
     175             : 
     176        1810 :    END SUBROUTINE compute_s_square
     177             : 
     178             : ! **************************************************************************************************
     179             : !> \brief restrains/constrains the value of s2 in a calculation
     180             : !> \param mos input
     181             : !> \param matrix_s input
     182             : !> \param mo_derivs inout if present, add the derivative of s_square wrt mos to mo_derivs
     183             : !> \param energy ...
     184             : !> \param s2_restraint_control ...
     185             : !> \param just_energy ...
     186             : !> \par History
     187             : !>      07.2004 created [ Joost VandeVondele ]
     188             : ! **************************************************************************************************
     189           0 :    SUBROUTINE s2_restraint(mos, matrix_s, mo_derivs, energy, &
     190             :                            s2_restraint_control, just_energy)
     191             : 
     192             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
     193             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     194             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT)      :: mo_derivs
     195             :       REAL(kind=dp)                                      :: energy
     196             :       TYPE(s2_restraint_type), POINTER                   :: s2_restraint_control
     197             :       LOGICAL                                            :: just_energy
     198             : 
     199             :       CHARACTER(len=*), PARAMETER                        :: routineN = 's2_restraint'
     200             : 
     201             :       INTEGER                                            :: handle
     202             :       REAL(kind=dp)                                      :: s_square, s_square_ideal
     203             : 
     204           0 :       CALL timeset(routineN, handle)
     205             : 
     206           0 :       SELECT CASE (s2_restraint_control%functional_form)
     207             :       CASE (do_s2_constraint)
     208           0 :          IF (just_energy) THEN
     209           0 :             CALL compute_s_square(mos, matrix_s, s_square, s_square_ideal)
     210             :          ELSE
     211             :             CALL compute_s_square(mos, matrix_s, s_square, s_square_ideal, &
     212           0 :                                   mo_derivs, s2_restraint_control%strength)
     213             :          END IF
     214           0 :          energy = s2_restraint_control%strength*(s_square - s2_restraint_control%target)
     215           0 :          s2_restraint_control%s2_order_p = s_square
     216             :       CASE (do_s2_restraint) ! not yet implemented
     217           0 :          CPABORT("")
     218             :       CASE DEFAULT
     219           0 :          CPABORT("")
     220             :       END SELECT
     221             : 
     222           0 :       CALL timestop(handle)
     223             : 
     224           0 :    END SUBROUTINE s2_restraint
     225             : 
     226             : END MODULE s_square_methods