(git:1f285aa)
s_square_methods.F
Go to the documentation of this file.
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  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  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  CALL timeset(routinen, handle)
91 
92  NULLIFY (context)
93  NULLIFY (fm_struct_tmp)
94  NULLIFY (local_data)
95  NULLIFY (para_env)
96 
97  SELECT CASE (SIZE(mos))
98  CASE (1)
99  ! Spin restricted case, i.e. nothing to do
100  s_square = 0.0_dp
101  s_square_ideal = 0.0_dp
102  ! Restraining or constraining (\cal S) does not make sense
103  cpassert(PRESENT(mo_derivs) .OR. PRESENT(strength))
104  CASE (2)
105  CALL get_mo_set(mo_set=mos(1), mo_coeff=c_alpha, homo=nalpha, nmo=nmo_alpha, nao=nao)
106  CALL get_mo_set(mo_set=mos(2), mo_coeff=c_beta, homo=nbeta, nmo=nmo_beta, nao=nao_beta)
107  cpassert(nao == nao_beta)
108  has_uniform_occupation_alpha = has_uniform_occupation(mo_set=mos(1), last_mo=nalpha)
109  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  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  s_square_ideal = real((nalpha - nbeta)*(nalpha - nbeta + 2), kind=dp)/4.0_dp
115  ! Create overlap matrix
116  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  nrow_global=nalpha, ncol_global=nbeta)
119  ! Prepare C(alpha)^T*S*C(beta)
120  CALL cp_fm_create(catscb, fm_struct_tmp, name="C(alpha)^T*S*C(beta)")
121  CALL cp_fm_struct_release(fm_struct_tmp)
122  ! Create S*C(beta)
123  CALL cp_fm_get_info(c_beta, matrix_struct=fm_struct_tmp)
124  CALL cp_fm_create(scb, fm_struct_tmp, name="S*C(beta)")
125  ! Compute S*C(beta)
126  CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, c_beta, scb, nbeta)
127  ! Compute C(alpha)^T*S*C(beta)
128  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  CALL cp_fm_get_info(catscb, local_data=local_data, nrow_local=nrow_local, ncol_local=ncol_local)
131  s2 = 0.0_dp
132  DO j = 1, ncol_local
133  DO i = 1, nrow_local
134  s2 = s2 + local_data(i, j)**2
135  END DO
136  END DO
137  CALL para_env%sum(s2)
138  s_square = s_square_ideal + nbeta - s2
139  ! Only needed for constraining or restraining (\cal S)
140  IF (PRESENT(mo_derivs)) THEN
141  cpassert(SIZE(mo_derivs, 1) == 2)
142  ! This gets really wrong for fractional occupations
143  CALL get_mo_set(mo_set=mos(1), uniform_occupation=has_uniform_occupation_alpha)
144  cpassert(has_uniform_occupation_alpha)
145  CALL get_mo_set(mo_set=mos(2), uniform_occupation=has_uniform_occupation_beta)
146  cpassert(has_uniform_occupation_beta)
147  ! Add -strength*S*C(beta)*(C(alpha)^T*S*C(beta))^T to the alpha MO derivatives
148  CALL parallel_gemm('N', 'T', nao, nalpha, nbeta, -strength, scb, catscb, 1.0_dp, mo_derivs(1))
149  ! Create S*C(alpha)
150  CALL cp_fm_get_info(c_alpha, matrix_struct=fm_struct_tmp)
151  CALL cp_fm_create(sca, fm_struct_tmp, name="S*C(alpha)")
152  ! Compute S*C(alpha)
153  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  CALL parallel_gemm('N', 'N', nao, nbeta, nalpha, -strength, sca, catscb, 1.0_dp, mo_derivs(2))
156  CALL cp_fm_release(sca)
157  END IF
158  CALL cp_fm_release(scb)
159  CALL cp_fm_release(catscb)
160  ELSE
161  IF (.NOT. has_uniform_occupation_alpha) THEN
162  cphint("The alpha orbitals (MOs) have a non-uniform occupation")
163  END IF
164  IF (.NOT. has_uniform_occupation_beta) THEN
165  cphint("The beta orbitals (MOs) have a non-uniform occupation")
166  END IF
167  cphint("Skipping S**2 calculation")
168  END IF
169  CASE DEFAULT
170  ! We should never get here
171  cpabort("Alpha, beta, what else ?")
172  END SELECT
173 
174  CALL timestop(handle)
175 
176  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  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  CALL timeset(routinen, handle)
205 
206  SELECT CASE (s2_restraint_control%functional_form)
207  CASE (do_s2_constraint)
208  IF (just_energy) THEN
209  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  mo_derivs, s2_restraint_control%strength)
213  END IF
214  energy = s2_restraint_control%strength*(s_square - s2_restraint_control%target)
215  s2_restraint_control%s2_order_p = s_square
216  CASE (do_s2_restraint) ! not yet implemented
217  cpabort("")
218  CASE DEFAULT
219  cpabort("")
220  END SELECT
221 
222  CALL timestop(handle)
223 
224  END SUBROUTINE s2_restraint
225 
226 END MODULE s_square_methods
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_s2_restraint
integer, parameter, public do_s2_restraint
Definition: input_constants.F:754
input_constants::do_s2_constraint
integer, parameter, public do_s2_constraint
Definition: input_constants.F:754
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_mo_types::has_uniform_occupation
logical function, public has_uniform_occupation(mo_set, first_mo, last_mo, occupation, tolerance)
Check if the set of MOs in mo_set specifed by the MO index range [first_mo,last_mo] an integer occupa...
Definition: qs_mo_types.F:503
s_square_methods
Methods related to (\cal S)^2 (i.e. spin)
Definition: s_square_methods.F:15
s_square_methods::s2_restraint
subroutine, public s2_restraint(mos, matrix_s, mo_derivs, energy, s2_restraint_control, just_energy)
restrains/constrains the value of s2 in a calculation
Definition: s_square_methods.F:191
s_square_methods::compute_s_square
subroutine, public compute_s_square(mos, matrix_s, s_square, s_square_ideal, mo_derivs, strength)
Compute the expectation value <(\cal S)^2> of the single determinant defined by the spin up (alpha) a...
Definition: s_square_methods.F:66