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 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_api, ONLY: dbcsr_p_type
20 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
21 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
22 : cp_fm_struct_release,&
23 : 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 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 1846 : 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 1846 : 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 1846 : CALL timeset(routineN, handle)
91 :
92 1846 : NULLIFY (context)
93 1846 : NULLIFY (fm_struct_tmp)
94 1846 : NULLIFY (local_data)
95 1846 : NULLIFY (para_env)
96 :
97 1846 : 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 1846 : CALL get_mo_set(mo_set=mos(1), mo_coeff=c_alpha, homo=nalpha, nmo=nmo_alpha, nao=nao)
106 1846 : CALL get_mo_set(mo_set=mos(2), mo_coeff=c_beta, homo=nbeta, nmo=nmo_beta, nao=nao_beta)
107 1846 : CPASSERT(nao == nao_beta)
108 1846 : has_uniform_occupation_alpha = has_uniform_occupation(mo_set=mos(1), last_mo=nalpha)
109 1846 : 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 1846 : 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 1846 : s_square_ideal = REAL((nalpha - nbeta)*(nalpha - nbeta + 2), KIND=dp)/4.0_dp
115 : ! Create overlap matrix
116 1846 : 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 1846 : nrow_global=nalpha, ncol_global=nbeta)
119 : ! Prepare C(alpha)^T*S*C(beta)
120 1846 : CALL cp_fm_create(catscb, fm_struct_tmp, name="C(alpha)^T*S*C(beta)")
121 1846 : CALL cp_fm_struct_release(fm_struct_tmp)
122 : ! Create S*C(beta)
123 1846 : CALL cp_fm_get_info(c_beta, matrix_struct=fm_struct_tmp)
124 1846 : CALL cp_fm_create(scb, fm_struct_tmp, name="S*C(beta)", set_zero=.TRUE.)
125 : ! Compute S*C(beta)
126 1846 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, c_beta, scb, nbeta)
127 : ! Compute C(alpha)^T*S*C(beta)
128 1846 : 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 1846 : CALL cp_fm_get_info(catscb, local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
131 1846 : s2 = 0.0_dp
132 8444 : DO j = 1, ncol_local
133 51297 : DO i = 1, nrow_local
134 49451 : s2 = s2 + local_data(i, j)**2
135 : END DO
136 : END DO
137 1846 : CALL para_env%sum(s2)
138 1846 : s_square = s_square_ideal + nbeta - s2
139 : ! Only needed for constraining or restraining (\cal S)
140 1846 : 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)", set_zero=.TRUE.)
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 1846 : CALL cp_fm_release(scb)
159 1846 : 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 1846 : CPABORT("Alpha, beta, what else?")
172 : END SELECT
173 :
174 1846 : CALL timestop(handle)
175 :
176 1846 : 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
|
