Coverage for gpaw/test/test_ibz2bz.py: 97%
78 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-19 00:19 +0000
1import numpy as np
2import pytest
3from gpaw import GPAW
4from gpaw.ibz2bz import IBZ2BZMaps, get_overlap, get_overlap_coefficients
5import gpaw.mpi as mpi
6from gpaw.test.gpwfile import response_band_cutoff
9@pytest.mark.serial
10@pytest.mark.parametrize(
11 'gs',
12 ['bi2i6_pw',
13 'fancy_si_pw',
14 'al_pw',
15 'fe_pw',
16 'co_pw',
17 'gaas_pw',
18 pytest.param(
19 'v2br4_pw',
20 marks=[pytest.mark.old_gpaw_only]), # interpolation=3 not implemented
21 'srvo3_pw'])
22def test_ibz2bz(in_tmp_dir, gpw_files, gs):
23 """ Tests gpaw.ibz2bz.py
24 Tests functionalities to take wavefunction and projections from
25 ibz to full bz by comparing calculations with and without symmetry.
26 """
28 atol = 5e-03 # Tolerance when comparing wfs and projections
29 atol_eig = 2e-04 # Tolerance when comparing eigenvalues
30 atol_deg = 5e-3 # Tolerance for checking degenerate states
32 # Loading calc with symmetry
33 calc = GPAW(gpw_files[gs],
34 communicator=mpi.serial_comm)
35 wfs = calc.wfs
36 nconv = response_band_cutoff[gs]
38 # setting basic stuff
39 nbands = wfs.bd.nbands if nconv == -1 else nconv
40 nbzk = wfs.kd.nbzkpts
41 ibz2bz = IBZ2BZMaps.from_calculator(calc)
43 # Loading calc without symmetry
44 calc_nosym = GPAW(gpw_files[gs + '_nosym'],
45 communicator=mpi.serial_comm)
46 wfs_nosym = calc_nosym.wfs
48 # Check some basic stuff
49 assert wfs_nosym.kd.nbzkpts == wfs_nosym.kd.nibzkpts
51 # Loop over spins and k-points
52 for s in range(wfs.nspins):
53 for K in range(nbzk):
54 ik = wfs.kd.bz2ibz_k[K] # IBZ k-point
56 # Check so that BZ kpoints are the same
57 assert np.allclose(wfs.kd.bzk_kc[K], wfs_nosym.kd.bzk_kc[K])
58 assert np.allclose(wfs_nosym.kd.ibzk_kc[K], wfs_nosym.kd.bzk_kc[K])
60 # Get data for calc without symmetry at BZ kpt K
61 eps_n_nosym, ut_nR_nosym, proj_nosym, dO_aii_nosym = \
62 get_ibz_data_from_wfs(wfs_nosym, nbands, K, s)
64 # Get data for calc with symmetry at ibz kpt ik
65 eps_n, ut_nR, proj, dO_aii = get_ibz_data_from_wfs(wfs,
66 nbands,
67 ik, s)
69 # Map projections and u:s from ik to K
70 proj_sym = ibz2bz[K].map_projections(proj)
71 ut_nR_sym = np.array([ibz2bz[K].map_pseudo_wave_to_BZ(
72 ut_nR[n]) for n in range(nbands)])
74 # Check so that eigenvalues are the same
75 assert np.allclose(eps_n[:nbands],
76 eps_n_nosym[:nbands],
77 atol=atol_eig)
79 # Check so that overlaps are the same for both calculations
80 assert equal_dicts(dO_aii,
81 dO_aii_nosym,
82 atol)
84 # Here starts the actual test
85 # Loop over all bands
86 n = 0
87 while n < nbands:
88 dim = find_degenerate_subspace(eps_n, n, nbands, atol_deg)
89 if dim == 1:
91 # First check untransformed quantities for ibz k-points
92 if np.allclose(wfs.kd.bzk_kc[K],
93 wfs.kd.ibzk_kc[ik]):
94 # Compare untransformed projections
95 compare_projections(proj, proj_nosym, n, atol)
96 # Compare untransformed wf:s
97 assert np.allclose(abs(ut_nR[n]),
98 abs(ut_nR_nosym[n]), atol=atol)
100 # Then check so that absolute value of transformed
101 # projections are the same
102 compare_projections(proj_sym, proj_nosym, n, atol)
104 # Check so that periodic part of pseudo is same,
105 # up to a phase
106 assert np.allclose(abs(ut_nR_sym[n]),
107 abs(ut_nR_nosym[n]), atol=atol)
109 # For degenerate states check transformation
110 # matrix is unitary,
111 # For non-degenerate states check so that all-electron wf:s
112 # are the same up to phase
113 bands = range(n, n + dim)
115 check_all_electron_wfs(bands, wfs.gd,
116 ut_nR_sym, ut_nR_nosym,
117 proj_sym, proj_nosym, dO_aii,
118 atol)
119 n += dim
122def equal_dicts(dict_1, dict_2, atol):
123 """ Checks so that two dicts with np.arrays are
124 equal"""
125 assert len(dict_1.keys()) == len(dict_2.keys())
126 for key in dict_1:
127 # Make sure the dictionaries contain the same set of keys
128 if key not in dict_2:
129 return False
130 # Make sure that the arrays are identical
131 if not np.allclose(dict_1[key], dict_2[key], atol=atol):
132 return False
133 return True
136def find_degenerate_subspace(eps_n, n_start, nbands, atol_eig):
137 # Find degenerate eigenvalues
138 n = n_start
139 dim = 1
140 while n < nbands - 1 and abs(eps_n[n] - eps_n[n + 1]) < atol_eig:
141 dim += 1
142 n += 1
143 return dim
146def compare_projections(proj_sym, proj_nosym, n, atol):
147 # compares so that projections at given k and band index n
148 # differ by at most a phase
149 for a, P_ni in proj_sym.items():
150 for j in range(P_ni.shape[1]):
151 # Check so that absolute values of projections are the same
152 assert np.isclose(abs(P_ni[n, j]),
153 abs(proj_nosym[a][n, j]),
154 atol=atol)
157def check_all_electron_wfs(bands, gd, ut1_nR, ut2_nR,
158 proj_sym, proj_nosym, dO_aii,
159 atol):
160 """sets up transformation matrix between symmetry
161 transformed u:s and normal u:s in degenerate subspace
162 and asserts that it is unitary. It also checks that
163 the pseudo wf:s transform according to the same
164 transformation.
166 Let |ψ^1_i> denote the all electron wavefunctions
167 from the calculation with symmetry and |ψ^2_i>
168 the corresponding wavefunctions from the calculation
169 without symmetry.
170 If the set {|ψ^1_i>} span the same subspace as the set
171 {|ψ^2_i>} they fulfill the following where summation
172 over repeated indexes is assumed:
174 |ψ^2_i> = |ψ^1_k> <ψ^1_k |ψ^2_i> == M_ki |ψ^1_k>
175 and M_ki = <ψ^1_k |ψ^2_i> is a unitary transformation.
176 M_ki is only unitary if the two sets of wfs span the
177 same subspace.
179 Parameters
180 ---------
181 bands: list of ints
182 band indexes in degenerate subspace
183 ut1_nR: np.array
184 ut2_nR: np.array
185 Periodic part of pseudo wave function for two calculations
186 proj_sym: Projections object
187 proj_nosym: Projections object
188 Projections for two calculations
189 dO_aii: dict with np.arrays
190 see ibz2bz.get_overlap_coefficients
191 dv: float
192 calc.wfs.gd.dv
193 atol: float
194 absolute tolerance when comparing arrays
195 """
196 M_nn = get_overlap(bands,
197 gd,
198 ut1_nR,
199 ut2_nR,
200 proj_sym,
201 proj_nosym,
202 dO_aii)
204 # Check so that transformation matrix is unitary
205 MMdag_nn = M_nn @ M_nn.T.conj()
206 assert np.allclose(np.eye(len(MMdag_nn)), MMdag_nn, atol=atol)
208 # Check so that M_nn transforms pseudo wf:s, see docs
209 ut2_from_transform_nR = np.einsum('ji,jklm->iklm', M_nn, ut1_nR[bands])
210 assert np.allclose(ut2_from_transform_nR, ut2_nR[bands], atol=atol)
213def get_ibz_data_from_wfs(wfs, nbands, ik, s):
214 """ gets data at ibz k-point ik
215 """
216 # get energies and wfs
217 kpt = wfs.kpt_qs[ik][s]
218 psit_nG = kpt.psit_nG
219 eps_n = kpt.eps_n
221 # Get periodic part of pseudo wfs
222 ut_nR = np.array([wfs.pd.ifft(
223 psit_nG[n], ik) for n in range(nbands)])
225 # Get projections
226 proj = kpt.projections.new(nbands=nbands, bcomm=None)
227 proj.array[:] = kpt.projections.array[:nbands]
229 # Get overlap coefficients from PAW setups
230 dO_aii = get_overlap_coefficients(wfs)
231 return eps_n, ut_nR, proj, dO_aii