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 Utilities for transforming k-point MO coefficients under symmetry operations
10 : ! **************************************************************************************************
11 : MODULE kpoint_mo_symmetry_methods
12 : USE cp_cfm_types, ONLY: cp_cfm_get_info,&
13 : cp_cfm_get_submatrix,&
14 : cp_cfm_set_submatrix,&
15 : cp_cfm_to_cfm,&
16 : cp_cfm_type
17 : USE kinds, ONLY: dp
18 : USE kpoint_types, ONLY: find_kpoint_rotation_slot,&
19 : kind_rotmat_type,&
20 : kpoint_sym_type,&
21 : kpoint_type
22 : USE mathconstants, ONLY: twopi
23 : #include "./base/base_uses.f90"
24 :
25 : IMPLICIT NONE
26 : PRIVATE
27 :
28 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kpoint_mo_symmetry_methods'
29 :
30 : PUBLIC :: kpoint_same_periodic, &
31 : kpoint_transform_scf_mo, kpoint_cell_gauge_action
32 :
33 : CONTAINS
34 :
35 : ! **************************************************************************************************
36 : !> \brief AO action in the cell Fourier gauge H(k)=sum_R exp(+ikR) H(R).
37 : !> \param kpoints Symmetry geometry and Gaussian basis rotation blocks
38 : !> \param representative Source irreducible point
39 : !> \param operation Orbit operation; zero identity, minus one pure time reversal
40 : !> \param TARGET Target fractional k point, modulo reciprocal vectors
41 : !> \param action Scalar AO action, excluding conjugation and physical spin rotation
42 : !> \param rotation Cartesian spatial operation
43 : !> \param antiunitary Whether the source coefficients must be conjugated
44 : !> \param status Zero on success, negative for unavailable or inconsistent mappings
45 : !> \note An atom shifted to image L contributes exp(-i*k_target*L). No atomic-gauge
46 : !> reciprocal-vector phase or legacy SCF phase calibration belongs in this cell gauge.
47 : ! **************************************************************************************************
48 16 : SUBROUTINE kpoint_cell_gauge_action(kpoints, representative, operation, TARGET, action, &
49 : rotation, antiunitary, status)
50 : TYPE(kpoint_type), POINTER :: kpoints
51 : INTEGER, INTENT(IN) :: representative, operation
52 : REAL(dp), INTENT(IN) :: target(3)
53 : COMPLEX(dp), INTENT(OUT) :: action(:, :)
54 : REAL(dp), INTENT(OUT) :: rotation(3, 3)
55 : LOGICAL, INTENT(OUT) :: antiunitary
56 : INTEGER, INTENT(OUT) :: status
57 :
58 : INTEGER :: a, b, first, i, kind, n, na, nb, slot
59 16 : INTEGER, ALLOCATABLE :: offset(:), seen(:), width(:)
60 : REAL(dp) :: angle
61 : TYPE(kpoint_sym_type), POINTER :: sym
62 :
63 8128 : action = 0.0_dp
64 16 : rotation = 0.0_dp
65 16 : antiunitary = .FALSE.
66 16 : status = -1
67 16 : n = SIZE(action, 1)
68 16 : IF (n < 1 .OR. SIZE(action, 2) /= n .OR. .NOT. ASSOCIATED(kpoints)) RETURN
69 16 : IF (representative < 1 .OR. representative > kpoints%nkp) RETURN
70 16 : IF (operation == 0 .OR. operation == -1) THEN
71 2 : antiunitary = operation == -1
72 2 : IF (operation == 0) THEN
73 0 : IF (.NOT. kpoint_same_periodic(TARGET, kpoints%xkp(:, representative))) RETURN
74 : ELSE
75 8 : IF (.NOT. kpoint_same_periodic(TARGET, -kpoints%xkp(:, representative))) RETURN
76 : END IF
77 28 : DO i = 1, n
78 28 : action(i, i) = 1.0_dp
79 : END DO
80 8 : DO i = 1, 3
81 8 : rotation(i, i) = 1.0_dp
82 : END DO
83 2 : status = 0
84 2 : RETURN
85 : END IF
86 14 : IF (operation < 1 .OR. .NOT. ASSOCIATED(kpoints%kp_sym)) RETURN
87 14 : sym => kpoints%kp_sym(representative)%kpoint_sym
88 14 : IF (.NOT. ASSOCIATED(sym)) RETURN
89 14 : IF (.NOT. sym%apply_symmetry .OR. operation > sym%nwred) RETURN
90 14 : IF (.NOT. kpoint_same_periodic(TARGET, sym%xkp(:, operation))) RETURN
91 14 : IF (.NOT. ASSOCIATED(kpoints%atype) .OR. .NOT. ASSOCIATED(kpoints%kind_rotmat)) RETURN
92 14 : IF (.NOT. ASSOCIATED(sym%fcell_gauge)) RETURN
93 14 : slot = find_kpoint_rotation_slot(kpoints, sym%rotp(operation))
94 14 : IF (slot == 0) RETURN
95 14 : na = SIZE(kpoints%atype)
96 70 : ALLOCATE (offset(na), width(na), seen(na))
97 14 : first = 1
98 14 : seen = 0
99 38 : DO a = 1, na
100 24 : kind = kpoints%atype(a)
101 24 : IF (kind < 1 .OR. kind > SIZE(kpoints%kind_rotmat, 2)) RETURN
102 24 : IF (.NOT. ASSOCIATED(kpoints%kind_rotmat(slot, kind)%rmat)) RETURN
103 24 : offset(a) = first
104 24 : width(a) = SIZE(kpoints%kind_rotmat(slot, kind)%rmat, 2)
105 38 : first = first + width(a)
106 : END DO
107 14 : IF (first /= n + 1) RETURN
108 38 : DO a = 1, na
109 24 : b = sym%f0(a, operation)
110 24 : IF (b < 1 .OR. b > na) RETURN
111 24 : seen(b) = seen(b) + 1
112 24 : IF (seen(b) /= 1 .OR. kpoints%atype(b) /= kpoints%atype(a)) RETURN
113 24 : kind = kpoints%atype(a)
114 24 : nb = width(b)
115 24 : IF (SIZE(kpoints%kind_rotmat(slot, kind)%rmat, 1) /= nb) RETURN
116 96 : angle = -twopi*DOT_PRODUCT(TARGET, REAL(sym%fcell_gauge(:, a, operation), dp))
117 : action(offset(b):offset(b) + nb - 1, offset(a):offset(a) + width(a) - 1) = &
118 4406 : CMPLX(COS(angle), SIN(angle), dp)*kpoints%kind_rotmat(slot, kind)%rmat
119 : END DO
120 182 : rotation = sym%rot(:, :, operation)
121 14 : antiunitary = sym%rotp(operation) < 0
122 14 : status = 0
123 16 : END SUBROUTINE kpoint_cell_gauge_action
124 :
125 : ! **************************************************************************************************
126 : !> \brief Compare two fractional k-points modulo reciprocal lattice vectors.
127 : !> \param xkp_a first k-point
128 : !> \param xkp_b second k-point
129 : !> \return true if the k-points are equivalent
130 : ! **************************************************************************************************
131 155900 : LOGICAL FUNCTION kpoint_same_periodic(xkp_a, xkp_b) RESULT(same)
132 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp_a, xkp_b
133 :
134 : REAL(KIND=dp), DIMENSION(3) :: diff
135 :
136 623600 : diff(1:3) = xkp_a(1:3) - xkp_b(1:3)
137 623600 : diff(1:3) = diff(1:3) - REAL(NINT(diff(1:3)), KIND=dp)
138 623600 : same = SUM(ABS(diff(1:3))) < 1.0e-8_dp
139 :
140 155900 : END FUNCTION kpoint_same_periodic
141 :
142 : ! **************************************************************************************************
143 : !> \brief Transform one SCF MO coefficient matrix to an equivalent full-mesh k-point.
144 : !> \param src source MO coefficients
145 : !> \param dst transformed MO coefficients
146 : !> \param qs_kpoint SCF k-point object containing symmetry operations
147 : !> \param ikred representative k-point index
148 : !> \param isym symmetry operation index; 0 direct copy, -1 time reversal
149 : !> \param success true if the transformation was performed
150 : !> \param reason diagnostic message
151 : ! **************************************************************************************************
152 2728 : SUBROUTINE kpoint_transform_scf_mo(src, dst, qs_kpoint, ikred, &
153 : isym, success, reason)
154 : TYPE(cp_cfm_type), INTENT(IN) :: src, dst
155 : TYPE(kpoint_type), POINTER :: qs_kpoint
156 : INTEGER, INTENT(IN) :: ikred, isym
157 : LOGICAL, INTENT(OUT) :: success
158 : CHARACTER(LEN=*), INTENT(OUT) :: reason
159 :
160 : COMPLEX(KIND=dp) :: coeff, phase
161 2728 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dst_coeff, src_coeff
162 : INTEGER :: iao, iatom, ikind, imo, irow, irow_source, irow_target, irow_work, nao, natom, &
163 : nmo, rot_slot, source_atom, source_dim, target_atom, target_dim
164 2728 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_size, ao_start
165 : LOGICAL :: reverse_phase, time_reversal
166 : REAL(KIND=dp) :: arg, coskl, sinkl
167 : REAL(KIND=dp), DIMENSION(3) :: xkp_phase
168 2728 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rotmat
169 2728 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kind_rot
170 : TYPE(kpoint_sym_type), POINTER :: kpsym
171 :
172 2728 : success = .FALSE.
173 2728 : reason = ""
174 2728 : IF (isym == 0) THEN
175 126 : CALL cp_cfm_to_cfm(src, dst)
176 126 : success = .TRUE.
177 126 : RETURN
178 : END IF
179 :
180 2602 : CALL cp_cfm_get_info(src, nrow_global=nao, ncol_global=nmo)
181 15612 : ALLOCATE (src_coeff(nao, nmo), dst_coeff(nao, nmo))
182 2602 : CALL cp_cfm_get_submatrix(src, src_coeff)
183 2602 : dst_coeff(:, :) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
184 :
185 2602 : IF (isym == -1) THEN
186 10874 : dst_coeff(:, :) = CONJG(src_coeff)
187 62 : CALL cp_cfm_set_submatrix(dst, dst_coeff)
188 62 : success = .TRUE.
189 62 : RETURN
190 : END IF
191 :
192 2540 : IF (.NOT. ASSOCIATED(qs_kpoint%kp_sym)) THEN
193 0 : reason = "SCF symmetry operation data are not available"
194 0 : RETURN
195 : END IF
196 2540 : kpsym => qs_kpoint%kp_sym(ikred)%kpoint_sym
197 2540 : IF (.NOT. ASSOCIATED(kpsym)) THEN
198 0 : reason = "SCF k-point symmetry operation is not available"
199 0 : RETURN
200 : END IF
201 2540 : IF (isym > kpsym%nwred) THEN
202 0 : reason = "SCF k-point symmetry operation index is out of range"
203 0 : RETURN
204 : END IF
205 2540 : IF (.NOT. ASSOCIATED(qs_kpoint%atype) .OR. .NOT. ASSOCIATED(qs_kpoint%kind_rotmat)) THEN
206 0 : reason = "SCF atom mappings or basis rotation matrices are not available"
207 0 : RETURN
208 : END IF
209 :
210 2540 : rot_slot = find_kpoint_rotation_slot(qs_kpoint, kpsym%rotp(isym))
211 2540 : IF (rot_slot == 0) THEN
212 0 : reason = "could not match the SCF symmetry operation to a basis rotation"
213 0 : RETURN
214 : END IF
215 2540 : kind_rot => qs_kpoint%kind_rotmat(rot_slot, :)
216 2540 : natom = SIZE(qs_kpoint%atype)
217 10160 : ALLOCATE (ao_start(natom), ao_size(natom))
218 16572 : irow = 1
219 16572 : DO iatom = 1, natom
220 14032 : ikind = qs_kpoint%atype(iatom)
221 14032 : IF (.NOT. ASSOCIATED(kind_rot(ikind)%rmat)) THEN
222 0 : reason = "a required basis rotation matrix is not available"
223 0 : RETURN
224 : END IF
225 14032 : ao_start(iatom) = irow
226 14032 : ao_size(iatom) = SIZE(kind_rot(ikind)%rmat, 2)
227 16572 : irow = irow + ao_size(iatom)
228 : END DO
229 2540 : IF (irow - 1 /= nao) THEN
230 0 : reason = "atom-resolved AO dimensions do not match the MO coefficient matrix"
231 0 : RETURN
232 : END IF
233 :
234 2540 : time_reversal = kpsym%rotp(isym) < 0
235 2540 : reverse_phase = qs_kpoint%gamma_centered .AND. ANY(MOD(qs_kpoint%nkp_grid, 2) == 0)
236 2540 : IF (ASSOCIATED(kpsym%phase_mode)) THEN
237 2540 : IF (kpsym%phase_mode(isym) > 0) reverse_phase = kpsym%phase_mode(isym) == 2
238 : END IF
239 10160 : xkp_phase(1:3) = kpsym%xkp(1:3, isym)
240 16572 : DO iatom = 1, natom
241 14032 : source_atom = iatom
242 14032 : target_atom = kpsym%f0(iatom, isym)
243 14032 : ikind = qs_kpoint%atype(source_atom)
244 14032 : rotmat => kind_rot(ikind)%rmat
245 14032 : source_dim = ao_size(source_atom)
246 14032 : target_dim = ao_size(target_atom)
247 14032 : IF (SIZE(rotmat, 1) /= target_dim .OR. SIZE(rotmat, 2) /= source_dim) THEN
248 0 : reason = "basis rotation dimensions do not match the atom/AO symmetry transform"
249 0 : RETURN
250 : END IF
251 : arg = REAL(kpsym%fcell_gauge(1, source_atom, isym), KIND=dp)*xkp_phase(1) + &
252 : REAL(kpsym%fcell_gauge(2, source_atom, isym), KIND=dp)*xkp_phase(2) + &
253 14032 : REAL(kpsym%fcell_gauge(3, source_atom, isym), KIND=dp)*xkp_phase(3)
254 : ! Match symtrans_phase: SCF AOs use the lattice-periodic gauge, so reciprocal
255 : ! folding adds no atom-position phase to the coefficients.
256 14032 : IF (reverse_phase) arg = -arg
257 14032 : coskl = COS(twopi*arg)
258 14032 : sinkl = SIN(twopi*arg)
259 72780 : DO imo = 1, nmo
260 383056 : DO irow_work = 1, source_dim
261 312816 : irow_source = ao_start(source_atom) + irow_work - 1
262 312816 : coeff = src_coeff(irow_source, imo)
263 312816 : IF (time_reversal) coeff = CONJG(coeff)
264 312816 : phase = CMPLX(coskl, sinkl, KIND=dp)*coeff
265 2338800 : DO iao = 1, target_dim
266 1969776 : irow_target = ao_start(target_atom) + iao - 1
267 : dst_coeff(irow_target, imo) = dst_coeff(irow_target, imo) + &
268 2282592 : rotmat(iao, irow_work)*phase
269 : END DO
270 : END DO
271 : END DO
272 : END DO
273 :
274 2540 : CALL cp_cfm_set_submatrix(dst, dst_coeff)
275 2540 : success = .TRUE.
276 :
277 8058 : END SUBROUTINE kpoint_transform_scf_mo
278 :
279 : END MODULE kpoint_mo_symmetry_methods
|
