Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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_cfm_basic_linalg, ONLY: cp_cfm_schur_product_cc
19 : USE cp_cfm_types, ONLY: cp_cfm_create,&
20 : cp_cfm_get_info,&
21 : cp_cfm_release,&
22 : cp_cfm_to_fm,&
23 : cp_cfm_type
24 : USE cp_control_types, ONLY: s2_restraint_type
25 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
26 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
27 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
28 : cp_fm_row_scale,&
29 : cp_fm_schur_product
30 : USE cp_fm_pool_types, ONLY: cp_fm_pool_type,&
31 : fm_pool_create_fm,&
32 : fm_pool_give_back_fm
33 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
34 : cp_fm_struct_release,&
35 : cp_fm_struct_type
36 : USE cp_fm_types, ONLY: cp_fm_create,&
37 : cp_fm_get_info,&
38 : cp_fm_release,&
39 : cp_fm_sumall,&
40 : cp_fm_type
41 : USE input_constants, ONLY: do_s2_constraint,&
42 : do_s2_restraint
43 : USE kinds, ONLY: dp
44 : USE kpoint_types, ONLY: get_kpoint_env,&
45 : get_kpoint_info,&
46 : kpoint_env_type,&
47 : kpoint_type
48 : USE mathconstants, ONLY: one,&
49 : z_one,&
50 : z_zero,&
51 : zero
52 : USE message_passing, ONLY: mp_para_env_type
53 : USE parallel_gemm_api, ONLY: parallel_gemm
54 : USE qs_kpoint_operators, ONLY: kpoint_operator_context_create,&
55 : kpoint_operator_context_release,&
56 : kpoint_operator_context_type,&
57 : kpoint_operator_finish,&
58 : kpoint_operator_prepare,&
59 : kpoint_operator_start,&
60 : kpoint_slot_s,&
61 : kpoint_spin_free
62 : USE qs_matrix_pools, ONLY: mpools_get,&
63 : qs_matrix_pools_type
64 : USE qs_mo_types, ONLY: get_mo_set,&
65 : has_uniform_occupation,&
66 : mo_set_type
67 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
68 : #include "./base/base_uses.f90"
69 :
70 : IMPLICIT NONE
71 :
72 : PRIVATE
73 :
74 : ! Global parameters
75 :
76 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 's_square_methods'
77 :
78 : PUBLIC :: compute_s_square, compute_s_square_kp
79 : PUBLIC :: s2_restraint
80 :
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief Compute the expectation value <(\cal S)^2> of the single determinant defined by the
85 : !> spin up (alpha) and spin down (beta) orbitals
86 : !> \param mos [in] MO set with all MO information including the alpha and beta MO coefficients
87 : !> \param matrix_s [in] AO overlap matrix S (do not mix with the spin operator (\cal S))
88 : !> \param s_square [out] <(\cal S)^2> including potential spin contaminations
89 : !> \param s_square_ideal [out] Ideal value for <(\cal S)^2> without any spin contaminations
90 : !> \par History
91 : !> 07.2004 created (Joost VandeVondele)
92 : !> 08.2021 revised (Matthias Krack, MK)
93 : !> 10.2026 adapted and generalized to partial occupation and k-points
94 : !> \note
95 : !> See Eqs. 2.271 and 2.272 in Modern Quantum Chemistry by A. Szabo and N. S. Ostlund
96 : ! **************************************************************************************************
97 2472 : SUBROUTINE compute_s_square(mos, matrix_s, s_square, s_square_ideal)
98 :
99 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
100 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
101 : REAL(KIND=dp), INTENT(OUT) :: s_square, s_square_ideal
102 :
103 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_s_square'
104 :
105 : INTEGER :: handle, nalpha, nao, nao_beta, nbeta
106 : REAL(KIND=dp) :: nela, nelb, s2
107 2472 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:), POINTER :: occa, occb
108 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
109 : TYPE(cp_fm_type) :: asb, asb2, scb
110 : TYPE(cp_fm_type), POINTER :: c_alpha, c_beta
111 :
112 2472 : CALL timeset(routineN, handle)
113 :
114 2472 : CALL get_mo_set(mo_set=mos(1), mo_coeff=c_alpha, nmo=nalpha, nao=nao)
115 2472 : CALL get_mo_set(mo_set=mos(2), mo_coeff=c_beta, nmo=nbeta, nao=nao_beta)
116 2472 : CPASSERT(nao == nao_beta)
117 2472 : CALL get_mo_set(mo_set=mos(1), n_el_f=nela, occupation_numbers=occa)
118 2472 : CALL get_mo_set(mo_set=mos(2), n_el_f=nelb, occupation_numbers=occb)
119 : ! Eq. 2.272 in Modern Quantum Chemistry by A. Szabo and N. S. Ostlund
120 2472 : s_square_ideal = REAL((nela - nelb)*(nela - nelb + 2), KIND=dp)/4.0_dp
121 : ! Create MO overlap matrix
122 2472 : CALL cp_fm_get_info(c_beta, matrix_struct=fm_struct)
123 : ! Prepare C(alpha)^T*S*C(beta)
124 2472 : CALL cp_fm_create(asb, fm_struct, nrow=nalpha, ncol=nbeta)
125 2472 : CALL cp_fm_create(asb2, fm_struct, nrow=nalpha, ncol=nbeta)
126 2472 : CALL cp_fm_create(scb, fm_struct)
127 : ! Compute S*C(beta)
128 2472 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, c_beta, scb, nbeta)
129 : ! Compute C(alpha)^T*S*C(beta)
130 2472 : CALL parallel_gemm('T', 'N', nalpha, nbeta, nao, one, c_alpha, scb, zero, asb)
131 : ! square the matrix elements
132 2472 : CALL cp_fm_schur_product(asb, asb, asb2)
133 : ! occupation number weighting
134 2472 : CALL cp_fm_row_scale(asb2, occa)
135 2472 : CALL cp_fm_column_scale(asb2, occb)
136 2472 : s2 = cp_fm_sumall(asb2)
137 2472 : s_square = s_square_ideal + nelb - s2
138 2472 : CALL cp_fm_release(scb)
139 2472 : CALL cp_fm_release(asb)
140 2472 : CALL cp_fm_release(asb2)
141 :
142 2472 : CALL timestop(handle)
143 :
144 2472 : END SUBROUTINE compute_s_square
145 :
146 : ! **************************************************************************************************
147 : !> \brief Compute the expectation value <(\cal S)^2> of the single determinant defined by the
148 : !> spin up (alpha) and spin down (beta) orbitals
149 : !> \param mosin ...
150 : !> \param kpoints ...
151 : !> \param matrix_s ...
152 : !> \param s_square [out] <(\cal S)^2> including potential spin contaminations
153 : !> \param s_square_ideal [out] Ideal value for <(\cal S)^2> without any spin contaminations
154 : !> \par History
155 : !> 10.2026 adapted and generalized to partial occupation and k-points
156 : ! **************************************************************************************************
157 226 : SUBROUTINE compute_s_square_kp(mosin, kpoints, matrix_s, s_square, s_square_ideal)
158 :
159 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mosin
160 : TYPE(kpoint_type), INTENT(IN), POINTER :: kpoints
161 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
162 : REAL(KIND=dp), INTENT(OUT) :: s_square, s_square_ideal
163 :
164 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_s_square_kp'
165 :
166 : INTEGER :: handle, ikp, kplocal, kpmax, nalpha, &
167 : nao, nbeta
168 : INTEGER, DIMENSION(2) :: kp_range
169 226 : INTEGER, DIMENSION(:, :), POINTER :: kp_dist
170 : LOGICAL :: wfn_real
171 : REAL(KIND=dp) :: nela, nelb, s2, wkp
172 226 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:), POINTER :: occa, occb
173 : TYPE(cp_cfm_type) :: casb, casb2, cscb, csmat
174 : TYPE(cp_cfm_type), POINTER :: c_alpha, c_beta
175 : TYPE(cp_fm_pool_type), POINTER :: ao_ao_fm_pool
176 : TYPE(cp_fm_struct_type), POINTER :: fm_struct, matrix_struct
177 : TYPE(cp_fm_type) :: asb, asb2, fmglobal, fmlocal, rsmat, scb
178 : TYPE(cp_fm_type), POINTER :: moc, r_alpha, r_beta
179 : TYPE(kpoint_env_type), POINTER :: kp
180 226 : TYPE(kpoint_operator_context_type) :: op_ctx
181 226 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
182 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
183 226 : POINTER :: sab_nl
184 : TYPE(qs_matrix_pools_type), POINTER :: mpools
185 :
186 226 : CALL timeset(routineN, handle)
187 :
188 : CALL get_kpoint_info(kpoints, use_real_wfn=wfn_real, kp_range=kp_range, &
189 226 : kp_dist=kp_dist, sab_nl=sab_nl)
190 226 : CPASSERT(ASSOCIATED(sab_nl))
191 :
192 226 : kplocal = kp_range(2) - kp_range(1) + 1
193 634 : kpmax = MAXVAL(kp_dist(2, :) - kp_dist(1, :) + 1)
194 :
195 : ! fm pools to be used within a kpoint group
196 226 : CALL get_kpoint_info(kpoints, mpools=mpools)
197 226 : CALL mpools_get(mpools, ao_ao_fm_pool=ao_ao_fm_pool)
198 226 : CALL fm_pool_create_fm(ao_ao_fm_pool, fmlocal)
199 226 : CALL cp_fm_get_info(fmlocal, matrix_struct=matrix_struct)
200 226 : IF (wfn_real) THEN
201 2 : CALL cp_fm_create(rsmat, matrix_struct)
202 : ELSE
203 224 : CALL cp_cfm_create(csmat, matrix_struct)
204 : END IF
205 226 : CALL fm_pool_give_back_fm(ao_ao_fm_pool, fmlocal)
206 :
207 : ! fm matrix_structure for the global environment
208 226 : CALL get_mo_set(mo_set=mosin(1), mo_coeff=moc, nao=nao)
209 226 : CALL cp_fm_create(fmglobal, moc%matrix_struct, nrow=nao, ncol=nao)
210 : ! The overlap is spin free: the routine drives every instance under the
211 : ! single spin-free key.
212 226 : CALL kpoint_operator_context_create(op_ctx, kpoints, fmglobal, 1)
213 226 : CALL kpoint_operator_prepare(op_ctx, matrix_s, kpoint_slot_s, matrix_row=kpoint_spin_free)
214 992 : DO ikp = 1, kpmax
215 : CALL kpoint_operator_start(op_ctx, ikp, kpoint_spin_free, matrix_s, kpoint_slot_s, &
216 992 : use_grid=.FALSE., matrix_row=kpoint_spin_free)
217 : END DO
218 226 : CALL cp_fm_release(fmglobal)
219 :
220 226 : s2 = 0.0_dp
221 992 : DO ikp = 1, kpmax
222 766 : IF (wfn_real) THEN
223 2 : CALL kpoint_operator_finish(op_ctx, ikp, kpoint_spin_free, kpoint_slot_s, fm_re=rsmat)
224 : ELSE
225 764 : CALL kpoint_operator_finish(op_ctx, ikp, kpoint_spin_free, kpoint_slot_s, cmat=csmat)
226 : END IF
227 :
228 766 : IF (ikp > kplocal) CYCLE
229 :
230 762 : kp => kpoints%kp_env(ikp)%kpoint_env
231 762 : CALL get_kpoint_env(kp, wkp=wkp, mos=mos)
232 762 : CALL get_mo_set(mo_set=mos(1), occupation_numbers=occa)
233 762 : CALL get_mo_set(mo_set=mos(2), occupation_numbers=occb)
234 :
235 762 : IF (wfn_real) THEN
236 2 : CALL get_mo_set(mo_set=mos(1), mo_coeff=r_alpha, nmo=nalpha)
237 2 : CALL get_mo_set(mo_set=mos(2), mo_coeff=r_beta, nmo=nbeta)
238 2 : CALL cp_fm_get_info(r_beta, matrix_struct=fm_struct)
239 2 : CALL cp_fm_create(scb, fm_struct)
240 2 : CALL cp_fm_create(asb, fm_struct, nrow=nalpha, ncol=nbeta)
241 2 : CALL cp_fm_create(asb2, fm_struct, nrow=nalpha, ncol=nbeta)
242 2 : CALL parallel_gemm('N', 'N', nao, nbeta, nao, one, rsmat, r_beta, zero, scb)
243 2 : CALL parallel_gemm('T', 'N', nalpha, nbeta, nao, one, r_alpha, scb, zero, asb)
244 2 : CALL cp_fm_schur_product(asb, asb, asb2)
245 2 : CALL cp_fm_release(scb)
246 2 : CALL cp_fm_release(asb)
247 : ELSE
248 760 : CALL get_mo_set(mo_set=mos(1), cmo_coeff=c_alpha, nmo=nalpha)
249 760 : CALL get_mo_set(mo_set=mos(2), cmo_coeff=c_beta, nmo=nbeta)
250 760 : CALL cp_cfm_get_info(c_beta, matrix_struct=fm_struct)
251 760 : CALL cp_cfm_create(cscb, fm_struct)
252 760 : CALL cp_cfm_create(casb, fm_struct, nrow=nalpha, ncol=nbeta)
253 760 : CALL cp_cfm_create(casb2, fm_struct, nrow=nalpha, ncol=nbeta)
254 760 : CALL cp_fm_create(asb2, fm_struct, nrow=nalpha, ncol=nbeta)
255 760 : CALL parallel_gemm('N', 'N', nao, nbeta, nao, z_one, csmat, c_beta, z_zero, cscb)
256 760 : CALL parallel_gemm('C', 'N', nalpha, nbeta, nao, z_one, c_alpha, cscb, z_zero, casb)
257 760 : CALL cp_cfm_schur_product_cc(casb, casb, casb2)
258 760 : CALL cp_cfm_to_fm(casb2, asb2)
259 760 : CALL cp_cfm_release(cscb)
260 760 : CALL cp_cfm_release(casb)
261 760 : CALL cp_cfm_release(casb2)
262 : END IF
263 762 : CALL cp_fm_row_scale(asb2, occa)
264 762 : CALL cp_fm_column_scale(asb2, occb)
265 762 : s2 = s2 + wkp*cp_fm_sumall(asb2)
266 992 : CALL cp_fm_release(asb2)
267 :
268 : END DO
269 :
270 : !sum s2 over kp group
271 226 : CALL kpoints%para_env_inter_kp%sum(s2)
272 :
273 226 : ikp = 1
274 226 : IF (ikp <= kplocal) THEN
275 226 : kp => kpoints%kp_env(ikp)%kpoint_env
276 226 : CALL get_kpoint_env(kp, mos=mos)
277 226 : CALL get_mo_set(mo_set=mos(1), occupation_numbers=occa)
278 226 : CALL get_mo_set(mo_set=mos(2), occupation_numbers=occb)
279 1482 : nela = SUM(occa)
280 1482 : nelb = SUM(occb)
281 : ELSE
282 0 : nela = 0
283 0 : nelb = 0
284 : END IF
285 : !
286 226 : CALL kpoints%para_env_kp%max(nela)
287 226 : CALL kpoints%para_env_kp%max(nelb)
288 : !
289 226 : s_square_ideal = REAL((nela - nelb)*(nela - nelb + 2), KIND=dp)/4.0_dp
290 226 : s_square = s_square_ideal + nelb - s2
291 :
292 226 : CALL kpoint_operator_context_release(op_ctx)
293 :
294 226 : IF (wfn_real) THEN
295 2 : CALL cp_fm_release(rsmat)
296 : ELSE
297 224 : CALL cp_cfm_release(csmat)
298 : END IF
299 :
300 226 : CALL timestop(handle)
301 :
302 452 : END SUBROUTINE compute_s_square_kp
303 :
304 : ! **************************************************************************************************
305 : !> \brief restrains/constrains the value of s2 in a calculation
306 : !> \param mos input
307 : !> \param matrix_s input
308 : !> \param mo_derivs inout if present, add the derivative of s_square wrt mos to mo_derivs
309 : !> \param energy ...
310 : !> \param s2_restraint_control ...
311 : !> \param just_energy ...
312 : !> \par History
313 : !> 07.2004 created [ Joost VandeVondele ]
314 : ! **************************************************************************************************
315 0 : SUBROUTINE s2_restraint(mos, matrix_s, mo_derivs, energy, &
316 : s2_restraint_control, just_energy)
317 :
318 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
319 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
320 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: mo_derivs
321 : REAL(kind=dp) :: energy
322 : TYPE(s2_restraint_type), POINTER :: s2_restraint_control
323 : LOGICAL :: just_energy
324 :
325 : CHARACTER(len=*), PARAMETER :: routineN = 's2_restraint'
326 :
327 : INTEGER :: handle
328 : REAL(kind=dp) :: s_square, s_square_ideal
329 :
330 0 : CALL timeset(routineN, handle)
331 :
332 0 : SELECT CASE (s2_restraint_control%functional_form)
333 : CASE (do_s2_constraint)
334 0 : IF (just_energy) THEN
335 0 : CALL compute_s2_energyforce(mos, matrix_s, s_square, s_square_ideal)
336 : ELSE
337 : CALL compute_s2_energyforce(mos, matrix_s, s_square, s_square_ideal, &
338 0 : mo_derivs, s2_restraint_control%strength)
339 : END IF
340 0 : energy = s2_restraint_control%strength*(s_square - s2_restraint_control%target)
341 0 : s2_restraint_control%s2_order_p = s_square
342 : CASE (do_s2_restraint) ! not yet implemented
343 0 : CPABORT("do_s2_restraint for <S**2> not yet implemented")
344 : CASE DEFAULT
345 0 : CPABORT("Unknown <S**2> constraint/restraint form")
346 : END SELECT
347 :
348 0 : CALL timestop(handle)
349 :
350 0 : END SUBROUTINE s2_restraint
351 :
352 : ! **************************************************************************************************
353 : !> \brief Compute the expectation value <(\cal S)^2> of the single determinant defined by the
354 : !> spin up (alpha) and spin down (beta) orbitals
355 : !> \param mos [in] MO set with all MO information including the alpha and beta MO coefficients
356 : !> \param matrix_s [in] AO overlap matrix S (do not mix with the spin operator (\cal S))
357 : !> \param s_square [out] <(\cal S)^2> including potential spin contaminations
358 : !> \param s_square_ideal [out] Ideal value for <(\cal S)^2> without any spin contaminations
359 : !> \param mo_derivs [inout] If present, add the derivative of s_square wrt the MOs to mo_derivs
360 : !> \param strength [in] Strength for constraining or restraining (\cal S)^2
361 : !> \par History
362 : !> 07.2004 created (Joost VandeVondele)
363 : !> 08.2021 revised (Matthias Krack, MK)
364 : !> \note
365 : !> See Eqs. 2.271 and 2.272 in Modern Quantum Chemistry by A. Szabo and N. S. Ostlund
366 : ! **************************************************************************************************
367 0 : SUBROUTINE compute_s2_energyforce(mos, matrix_s, s_square, s_square_ideal, mo_derivs, strength)
368 :
369 : TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
370 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
371 : REAL(KIND=dp), INTENT(OUT) :: s_square, s_square_ideal
372 : TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT), &
373 : OPTIONAL :: mo_derivs
374 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: strength
375 :
376 : CHARACTER(len=*), PARAMETER :: routineN = 'compute_s2_energyforce'
377 :
378 : INTEGER :: handle, i, j, nalpha, nao, nao_beta, &
379 : nbeta, ncol_local, nmo_alpha, &
380 : nmo_beta, nrow_local
381 : LOGICAL :: has_uniform_occupation_alpha, &
382 : has_uniform_occupation_beta
383 : REAL(KIND=dp) :: s2
384 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
385 0 : POINTER :: local_data
386 : TYPE(cp_blacs_env_type), POINTER :: context
387 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
388 : TYPE(cp_fm_type) :: catscb, sca, scb
389 : TYPE(cp_fm_type), POINTER :: c_alpha, c_beta
390 : TYPE(mp_para_env_type), POINTER :: para_env
391 :
392 0 : CALL timeset(routineN, handle)
393 :
394 0 : NULLIFY (context)
395 0 : NULLIFY (fm_struct_tmp)
396 0 : NULLIFY (local_data)
397 0 : NULLIFY (para_env)
398 :
399 0 : SELECT CASE (SIZE(mos))
400 : CASE (1)
401 : ! Spin restricted case, i.e. nothing to do
402 0 : s_square = 0.0_dp
403 0 : s_square_ideal = 0.0_dp
404 : ! Restraining or constraining (\cal S) does not make sense
405 0 : CPASSERT(PRESENT(mo_derivs) .OR. PRESENT(strength))
406 : CASE (2)
407 0 : CALL get_mo_set(mo_set=mos(1), mo_coeff=c_alpha, homo=nalpha, nmo=nmo_alpha, nao=nao)
408 0 : CALL get_mo_set(mo_set=mos(2), mo_coeff=c_beta, homo=nbeta, nmo=nmo_beta, nao=nao_beta)
409 0 : CPASSERT(nao == nao_beta)
410 0 : has_uniform_occupation_alpha = has_uniform_occupation(mo_set=mos(1), last_mo=nalpha)
411 0 : has_uniform_occupation_beta = has_uniform_occupation(mo_set=mos(2), last_mo=nbeta)
412 : ! This makes only sense if we have uniform occupations for the alpha and beta spin MOs while
413 : ! ignoring potentially added MOs with zero occupation
414 0 : IF (has_uniform_occupation_alpha .AND. has_uniform_occupation_beta) THEN
415 : ! Eq. 2.272 in Modern Quantum Chemistry by A. Szabo and N. S. Ostlund
416 0 : s_square_ideal = REAL((nalpha - nbeta)*(nalpha - nbeta + 2), KIND=dp)/4.0_dp
417 : ! Create overlap matrix
418 0 : CALL cp_fm_get_info(c_alpha, para_env=para_env, context=context)
419 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
420 0 : nrow_global=nalpha, ncol_global=nbeta)
421 : ! Prepare C(alpha)^T*S*C(beta)
422 0 : CALL cp_fm_create(catscb, fm_struct_tmp, name="C(alpha)^T*S*C(beta)")
423 0 : CALL cp_fm_struct_release(fm_struct_tmp)
424 : ! Create S*C(beta)
425 0 : CALL cp_fm_get_info(c_beta, matrix_struct=fm_struct_tmp)
426 0 : CALL cp_fm_create(scb, fm_struct_tmp, name="S*C(beta)", set_zero=.TRUE.)
427 : ! Compute S*C(beta)
428 0 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, c_beta, scb, nbeta)
429 : ! Compute C(alpha)^T*S*C(beta)
430 0 : CALL parallel_gemm('T', 'N', nalpha, nbeta, nao, 1.0_dp, c_alpha, scb, 0.0_dp, catscb)
431 : ! Eq. 2.271 in Modern Quantum Chemistry by A. Szabo and N. S. Ostlund
432 0 : CALL cp_fm_get_info(catscb, local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
433 0 : s2 = 0.0_dp
434 0 : DO j = 1, ncol_local
435 0 : DO i = 1, nrow_local
436 0 : s2 = s2 + local_data(i, j)**2
437 : END DO
438 : END DO
439 0 : CALL para_env%sum(s2)
440 0 : s_square = s_square_ideal + nbeta - s2
441 : ! Only needed for constraining or restraining (\cal S)
442 0 : IF (PRESENT(mo_derivs)) THEN
443 0 : CPASSERT(SIZE(mo_derivs, 1) == 2)
444 : ! This gets really wrong for fractional occupations
445 0 : CALL get_mo_set(mo_set=mos(1), uniform_occupation=has_uniform_occupation_alpha)
446 0 : CPASSERT(has_uniform_occupation_alpha)
447 0 : CALL get_mo_set(mo_set=mos(2), uniform_occupation=has_uniform_occupation_beta)
448 0 : CPASSERT(has_uniform_occupation_beta)
449 : ! Add -strength*S*C(beta)*(C(alpha)^T*S*C(beta))^T to the alpha MO derivatives
450 0 : CALL parallel_gemm('N', 'T', nao, nalpha, nbeta, -strength, scb, catscb, 1.0_dp, mo_derivs(1))
451 : ! Create S*C(alpha)
452 0 : CALL cp_fm_get_info(c_alpha, matrix_struct=fm_struct_tmp)
453 0 : CALL cp_fm_create(sca, fm_struct_tmp, name="S*C(alpha)", set_zero=.TRUE.)
454 : ! Compute S*C(alpha)
455 0 : CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, c_alpha, sca, nalpha)
456 : ! Add -strength*S*C(alpha)*C(alpha)^T*S*C(beta) to the beta MO derivatives
457 0 : CALL parallel_gemm('N', 'N', nao, nbeta, nalpha, -strength, sca, catscb, 1.0_dp, mo_derivs(2))
458 0 : CALL cp_fm_release(sca)
459 : END IF
460 0 : CALL cp_fm_release(scb)
461 0 : CALL cp_fm_release(catscb)
462 : ELSE
463 0 : IF (.NOT. has_uniform_occupation_alpha) THEN
464 0 : CPHINT("The alpha orbitals (MOs) have a non-uniform occupation")
465 : END IF
466 0 : IF (.NOT. has_uniform_occupation_beta) THEN
467 0 : CPHINT("The beta orbitals (MOs) have a non-uniform occupation")
468 : END IF
469 0 : CPHINT("Skipping S**2 calculation")
470 : END IF
471 : CASE DEFAULT
472 : ! We should never get here
473 0 : CPABORT("Alpha, beta, what else?")
474 : END SELECT
475 :
476 0 : CALL timestop(handle)
477 :
478 0 : END SUBROUTINE compute_s2_energyforce
479 :
480 : END MODULE s_square_methods
|
