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_fm_types, ONLY: cp_fm_copy_general,&
13 : cp_fm_get_info,&
14 : cp_fm_get_submatrix,&
15 : cp_fm_set_submatrix,&
16 : cp_fm_type
17 : USE kinds, ONLY: dp
18 : USE kpoint_types, ONLY: kind_rotmat_type,&
19 : kpoint_sym_type,&
20 : kpoint_type
21 : USE mathconstants, ONLY: twopi
22 : USE message_passing, ONLY: mp_para_env_type
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
32 :
33 : CONTAINS
34 :
35 : ! **************************************************************************************************
36 : !> \brief Compare two fractional k-points modulo reciprocal lattice vectors.
37 : !> \param xkp_a first k-point
38 : !> \param xkp_b second k-point
39 : !> \return true if the k-points are equivalent
40 : ! **************************************************************************************************
41 139520 : LOGICAL FUNCTION kpoint_same_periodic(xkp_a, xkp_b) RESULT(same)
42 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: xkp_a, xkp_b
43 :
44 : REAL(KIND=dp), DIMENSION(3) :: diff
45 :
46 558080 : diff(1:3) = xkp_a(1:3) - xkp_b(1:3)
47 558080 : diff(1:3) = diff(1:3) - REAL(NINT(diff(1:3)), KIND=dp)
48 558080 : same = SUM(ABS(diff(1:3))) < 1.0e-8_dp
49 :
50 139520 : END FUNCTION kpoint_same_periodic
51 :
52 : ! **************************************************************************************************
53 : !> \brief Transform one SCF MO coefficient matrix to an equivalent full-mesh k-point.
54 : !> \param src_real real part of source MO coefficients
55 : !> \param src_imag imaginary part of source MO coefficients
56 : !> \param dst_real real part of transformed MO coefficients
57 : !> \param dst_imag imaginary part of transformed MO coefficients
58 : !> \param qs_kpoint SCF k-point object containing symmetry operations
59 : !> \param ikred representative k-point index
60 : !> \param isym symmetry operation index; 0 direct copy, -1 time reversal
61 : !> \param para_env global parallel environment
62 : !> \param success true if the transformation was performed
63 : !> \param reason diagnostic message
64 : ! **************************************************************************************************
65 1912 : SUBROUTINE kpoint_transform_scf_mo(src_real, src_imag, dst_real, dst_imag, qs_kpoint, ikred, &
66 : isym, para_env, success, reason)
67 : TYPE(cp_fm_type), INTENT(IN) :: src_real, src_imag, dst_real, dst_imag
68 : TYPE(kpoint_type), POINTER :: qs_kpoint
69 : INTEGER, INTENT(IN) :: ikred, isym
70 : TYPE(mp_para_env_type), POINTER :: para_env
71 : LOGICAL, INTENT(OUT) :: success
72 : CHARACTER(LEN=*), INTENT(OUT) :: reason
73 :
74 : INTEGER :: iao, iatom, ikind, imo, irow, irow_source, irow_target, irow_work, nao, natom, &
75 : nmo, rot_slot, source_atom, source_dim, target_atom, target_dim
76 1912 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ao_size, ao_start
77 : LOGICAL :: reverse_phase, time_reversal
78 : REAL(KIND=dp) :: arg, coeff_imag, coeff_real, coskl, &
79 : phase_imag, phase_real, sinkl
80 1912 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dst_i, dst_r, src_i, src_r
81 : REAL(KIND=dp), DIMENSION(3) :: xkp_phase
82 1912 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rotmat
83 1912 : TYPE(kind_rotmat_type), DIMENSION(:), POINTER :: kind_rot
84 : TYPE(kpoint_sym_type), POINTER :: kpsym
85 :
86 1912 : success = .FALSE.
87 1912 : reason = ""
88 1912 : IF (isym == 0) THEN
89 60 : CALL cp_fm_copy_general(src_real, dst_real, para_env)
90 60 : CALL cp_fm_copy_general(src_imag, dst_imag, para_env)
91 60 : success = .TRUE.
92 60 : RETURN
93 : END IF
94 :
95 1852 : CALL cp_fm_get_info(src_real, nrow_global=nao, ncol_global=nmo)
96 18520 : ALLOCATE (src_r(nao, nmo), src_i(nao, nmo), dst_r(nao, nmo), dst_i(nao, nmo))
97 1852 : CALL cp_fm_get_submatrix(src_real, src_r)
98 1852 : CALL cp_fm_get_submatrix(src_imag, src_i)
99 1852 : dst_r(:, :) = 0.0_dp
100 1852 : dst_i(:, :) = 0.0_dp
101 :
102 1852 : IF (isym == -1) THEN
103 10812 : dst_r(1:nao, 1:nmo) = src_r(1:nao, 1:nmo)
104 10812 : dst_i(1:nao, 1:nmo) = -src_i(1:nao, 1:nmo)
105 60 : CALL cp_fm_set_submatrix(dst_real, dst_r)
106 60 : CALL cp_fm_set_submatrix(dst_imag, dst_i)
107 60 : DEALLOCATE (src_r, src_i, dst_r, dst_i)
108 60 : success = .TRUE.
109 60 : RETURN
110 : END IF
111 :
112 1792 : IF (.NOT. ASSOCIATED(qs_kpoint%kp_sym)) THEN
113 0 : reason = "SCF symmetry operation data are not available"
114 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i)
115 0 : RETURN
116 : END IF
117 1792 : kpsym => qs_kpoint%kp_sym(ikred)%kpoint_sym
118 1792 : IF (.NOT. ASSOCIATED(kpsym)) THEN
119 0 : reason = "SCF k-point symmetry operation is not available"
120 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i)
121 0 : RETURN
122 : END IF
123 1792 : IF (isym > kpsym%nwred) THEN
124 0 : reason = "SCF k-point symmetry operation index is out of range"
125 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i)
126 0 : RETURN
127 : END IF
128 1792 : IF (.NOT. ASSOCIATED(qs_kpoint%atype) .OR. .NOT. ASSOCIATED(qs_kpoint%kind_rotmat)) THEN
129 0 : reason = "SCF atom mappings or basis rotation matrices are not available"
130 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i)
131 0 : RETURN
132 : END IF
133 :
134 1792 : rot_slot = find_kpoint_rotation_slot(qs_kpoint, kpsym%rotp(isym))
135 1792 : IF (rot_slot == 0) THEN
136 0 : reason = "could not match the SCF symmetry operation to a basis rotation"
137 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i)
138 0 : RETURN
139 : END IF
140 1792 : kind_rot => qs_kpoint%kind_rotmat(rot_slot, :)
141 1792 : natom = SIZE(qs_kpoint%atype)
142 7168 : ALLOCATE (ao_start(natom), ao_size(natom))
143 14724 : irow = 1
144 14724 : DO iatom = 1, natom
145 12932 : ikind = qs_kpoint%atype(iatom)
146 12932 : IF (.NOT. ASSOCIATED(kind_rot(ikind)%rmat)) THEN
147 0 : reason = "a required basis rotation matrix is not available"
148 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i, ao_start, ao_size)
149 0 : RETURN
150 : END IF
151 12932 : ao_start(iatom) = irow
152 12932 : ao_size(iatom) = SIZE(kind_rot(ikind)%rmat, 2)
153 14724 : irow = irow + ao_size(iatom)
154 : END DO
155 1792 : IF (irow - 1 /= nao) THEN
156 0 : reason = "atom-resolved AO dimensions do not match the MO coefficient matrix"
157 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i, ao_start, ao_size)
158 0 : RETURN
159 : END IF
160 :
161 1792 : time_reversal = kpsym%rotp(isym) < 0
162 1792 : reverse_phase = qs_kpoint%gamma_centered .AND. ANY(MOD(qs_kpoint%nkp_grid, 2) == 0)
163 1792 : IF (ASSOCIATED(kpsym%phase_mode)) THEN
164 1792 : IF (kpsym%phase_mode(isym) > 0) reverse_phase = kpsym%phase_mode(isym) == 2
165 : END IF
166 7168 : xkp_phase(1:3) = kpsym%xkp(1:3, isym)
167 14724 : DO iatom = 1, natom
168 12932 : source_atom = iatom
169 12932 : target_atom = kpsym%f0(iatom, isym)
170 12932 : ikind = qs_kpoint%atype(source_atom)
171 12932 : rotmat => kind_rot(ikind)%rmat
172 12932 : source_dim = ao_size(source_atom)
173 12932 : target_dim = ao_size(target_atom)
174 12932 : IF (SIZE(rotmat, 1) /= target_dim .OR. SIZE(rotmat, 2) /= source_dim) THEN
175 0 : reason = "basis rotation dimensions do not match the atom/AO symmetry transform"
176 0 : DEALLOCATE (src_r, src_i, dst_r, dst_i, ao_start, ao_size)
177 0 : RETURN
178 : END IF
179 : arg = REAL(kpsym%fcell_gauge(1, source_atom, isym), KIND=dp)*xkp_phase(1) + &
180 : REAL(kpsym%fcell_gauge(2, source_atom, isym), KIND=dp)*xkp_phase(2) + &
181 12932 : REAL(kpsym%fcell_gauge(3, source_atom, isym), KIND=dp)*xkp_phase(3)
182 12932 : IF (ASSOCIATED(kpsym%kgphase)) THEN
183 12932 : arg = arg + kpsym%kgphase(source_atom, isym)
184 : END IF
185 12932 : IF (reverse_phase) arg = -arg
186 12932 : coskl = COS(twopi*arg)
187 12932 : sinkl = SIN(twopi*arg)
188 66968 : DO imo = 1, nmo
189 338684 : DO irow_work = 1, source_dim
190 273508 : irow_source = ao_start(source_atom) + irow_work - 1
191 273508 : coeff_real = src_r(irow_source, imo)
192 273508 : coeff_imag = src_i(irow_source, imo)
193 273508 : IF (time_reversal) coeff_imag = -coeff_imag
194 273508 : phase_real = coskl*coeff_real - sinkl*coeff_imag
195 273508 : phase_imag = sinkl*coeff_real + coskl*coeff_imag
196 1853036 : DO iao = 1, target_dim
197 1527284 : irow_target = ao_start(target_atom) + iao - 1
198 : dst_r(irow_target, imo) = dst_r(irow_target, imo) + &
199 1527284 : rotmat(iao, irow_work)*phase_real
200 : dst_i(irow_target, imo) = dst_i(irow_target, imo) + &
201 1800792 : rotmat(iao, irow_work)*phase_imag
202 : END DO
203 : END DO
204 : END DO
205 : END DO
206 :
207 1792 : CALL cp_fm_set_submatrix(dst_real, dst_r)
208 1792 : CALL cp_fm_set_submatrix(dst_imag, dst_i)
209 1792 : DEALLOCATE (src_r, src_i, dst_r, dst_i, ao_start, ao_size)
210 1792 : success = .TRUE.
211 :
212 5676 : END SUBROUTINE kpoint_transform_scf_mo
213 :
214 : ! **************************************************************************************************
215 : !> \brief Locate the basis-rotation slot corresponding to a signed k-point symmetry operation.
216 : !> \param qs_kpoint SCF k-point object
217 : !> \param rotp signed operation identifier
218 : !> \return rotation slot, or zero when no slot matches
219 : ! **************************************************************************************************
220 1792 : INTEGER FUNCTION find_kpoint_rotation_slot(qs_kpoint, rotp) RESULT(rot_slot)
221 : TYPE(kpoint_type), POINTER :: qs_kpoint
222 : INTEGER, INTENT(IN) :: rotp
223 :
224 : INTEGER :: irot, rot_abs
225 :
226 1792 : rot_slot = 0
227 1792 : rot_abs = ABS(rotp)
228 1792 : IF (.NOT. ASSOCIATED(qs_kpoint%ibrot)) RETURN
229 38536 : DO irot = 1, SIZE(qs_kpoint%ibrot)
230 38536 : IF (rot_abs == qs_kpoint%ibrot(irot)) THEN
231 1792 : rot_slot = irot
232 : RETURN
233 : END IF
234 : END DO
235 :
236 : END FUNCTION find_kpoint_rotation_slot
237 :
238 : END MODULE kpoint_mo_symmetry_methods
|
