Coverage for gpaw/test/response/test_site_kernels.py: 100%
270 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-20 00:19 +0000
1"""Test the site kernel calculation functionality of the response code"""
3# General modules
4import pytest
5import numpy as np
6import scipy.special as sc
8# Script modules
9from ase.build import bulk
11from gpaw import GPAW, PW
12from gpaw.response.site_kernels import (SphericalSiteKernels,
13 CylindricalSiteKernels,
14 ParallelepipedicSiteKernels,
15 sinc,
16 spherical_geometry_factor,
17 cylindrical_geometry_factor,
18 parallelepipedic_geometry_factor)
19from gpaw.response.pair_functions import get_pw_coordinates
22# ---------- Actual tests ---------- #
25pytestmark = pytest.mark.kspair
28@pytest.mark.ci
29def test_spherical_kernel(rng):
30 """Check the numerics of the spherical kernel"""
31 # ---------- Inputs ---------- #
33 # Relative wave vector lengths to check (relative to 1/rc)
34 Qrel_Q = np.array([0.,
35 np.pi / 2,
36 np.pi])
38 # Expected results (assuming |Q| * rc = 1.)
39 Vsphere = 4. * np.pi / 3.
40 test_K_Q = Vsphere * np.array([1.,
41 3 / (np.pi / 2)**3.,
42 3 / (np.pi)**2.])
44 # Spherical radii to check
45 nr = 5
46 rc_r = rng.random(nr)
48 # Wave vector directions to check
49 nd = 41
50 Q_dv = 2. * rng.random((nd, 3)) - 1.
51 Q_dv /= np.linalg.norm(Q_dv, axis=1)[:, np.newaxis] # normalize
53 # ---------- Script ---------- #
55 # Set up wave vectors
56 Q_Qdv = Qrel_Q[:, np.newaxis, np.newaxis] * Q_dv[np.newaxis, ...]
58 for rc in rc_r:
59 # Calculate site centered geometry factor with rescaled wave vector
60 K_Qd = spherical_geometry_factor(Q_Qdv / rc, rc)
61 # Check against expected result
62 assert np.allclose(K_Qd, rc**3. * test_K_Q[:, np.newaxis])
65@pytest.mark.ci
66def test_cylindrical_kernel(rng):
67 """Check the numerics of the spherical kernel"""
68 # ---------- Inputs ---------- #
70 # Relative wave vectors (relative to 1/rc) in radial direction
71 Qrhorel_Q1 = np.array([0.] + list(sc.jn_zeros(1, 4))) # Roots of J1(x)
73 # Relative wave vectors (relative to 2/hc) in cylindrical direction
74 Qzrel_Q2 = list(np.pi * np.arange(5)) # Roots of sin(x)
75 Qzrel_Q2 += list(np.pi * np.arange(4) + np.pi / 2) # Extrema of sin(x)
76 Qzrel_Q2 = np.array(Qzrel_Q2)
78 # Expected results for roots of J1 (assuming rc=1. and hc=2.)
79 Vcylinder = 2. * np.pi
80 nQ2 = 13 # Choose random Q_z r_z
81 test_Krho_Q1 = np.array([1., 0., 0., 0., 0.])
82 Qzrand_Q2 = 10. * rng.random(nQ2)
83 sinc_zrand_Q2 = np.sin(Qzrand_Q2) / Qzrand_Q2
84 test_Krho_Q1Q2 = Vcylinder * test_Krho_Q1[:, np.newaxis]\
85 * sinc_zrand_Q2[np.newaxis, :]
87 # Expected results for roots and extrema of sin (assuming rc=1. and hc=2.)
88 nQ1 = 15 # Choose random Q_ρ h_c
89 test_Kz_Q2 = [1., 0., 0., 0., 0.] # Nodes in sinc(Q_z h_c)
90 test_Kz_Q2 += list(np.array([1., -1., 1., -1.]) / Qzrel_Q2[5:]) # Extrema
91 test_Kz_Q2 = np.array(test_Kz_Q2)
92 Qrhorand_Q1 = 10. * rng.random(nQ1)
93 J1term_rhorand_Q1 = 2. * sc.jv(1, Qrhorand_Q1) / Qrhorand_Q1
94 test_Kz_Q1Q2 = Vcylinder * J1term_rhorand_Q1[:, np.newaxis]\
95 * test_Kz_Q2[np.newaxis, :]
97 # Cylinder radii to check
98 nr = 5
99 rc_r = 3. * rng.random(nr)
101 # Cylinder height to check
102 nh = 3
103 hc_h = 4. * rng.random(nh)
105 # Cylindrical axes to check
106 nc = 7
107 ez_cv = 2. * rng.random((nc, 3)) - 1.
108 ez_cv /= np.linalg.norm(ez_cv, axis=1)[:, np.newaxis]
110 # Wave vector directions in-plane to check. Generated through the cross
111 # product of a random direction with the cylindrical axis
112 nd = 11
113 Qrho_dv = 2. * rng.random((nd, 3)) - 1.
114 Qrho_cdv = np.cross(Qrho_dv[np.newaxis, ...], ez_cv[:, np.newaxis, :])
115 Qrho_cdv /= np.linalg.norm(Qrho_cdv, axis=-1)[..., np.newaxis] # normalize
117 # ---------- Script ---------- #
119 for rc in rc_r:
120 for hc in hc_h:
121 # Set up wave vectors for radial tests
122 Qrho_cdQ1v = Qrhorel_Q1[np.newaxis, np.newaxis, :, np.newaxis]\
123 * Qrho_cdv[..., np.newaxis, :] / rc
124 Qrho_cQ2v = Qzrand_Q2[np.newaxis, :, np.newaxis]\
125 * ez_cv[:, np.newaxis, :] / (hc / 2.)
126 Qrho_cdQ1Q2v = Qrho_cdQ1v[..., np.newaxis, :]\
127 + Qrho_cQ2v[:, np.newaxis, np.newaxis, ...]
129 # Set up wave vectors for cylindrical tests
130 Qz_cdQ1v = Qrhorand_Q1[np.newaxis, np.newaxis, :, np.newaxis]\
131 * Qrho_cdv[..., np.newaxis, :] / rc
132 Qz_cQ2v = Qzrel_Q2[np.newaxis, :, np.newaxis]\
133 * ez_cv[:, np.newaxis, :] / (hc / 2.)
134 Qz_cdQ1Q2v = Qz_cdQ1v[..., np.newaxis, :]\
135 + Qz_cQ2v[:, np.newaxis, np.newaxis, ...]
137 # Test one cylindrical direction at a time
138 for ez_v, Qrho_dQ1Q2v, Qz_dQ1Q2v in zip(ez_cv,
139 Qrho_cdQ1Q2v, Qz_cdQ1Q2v):
140 # Calculate geometry factors
141 Krho_dQ1Q2 = cylindrical_geometry_factor(Qrho_dQ1Q2v,
142 ez_v, rc, hc)
143 Kz_dQ1Q2 = cylindrical_geometry_factor(Qz_dQ1Q2v,
144 ez_v, rc, hc)
146 # Check against expected result
147 assert np.allclose(Krho_dQ1Q2, rc**2. * hc / 2.
148 * test_Krho_Q1Q2[np.newaxis, ...],
149 atol=1.e-8)
150 assert np.allclose(Kz_dQ1Q2, rc**2. * hc / 2.
151 * test_Kz_Q1Q2[np.newaxis, ...],
152 atol=1.e-8)
155@pytest.mark.ci
156def test_parallelepipedic_kernel(rng):
157 """Check the numerics of the parallelepipedic site kernel."""
158 # ---------- Inputs ---------- #
160 # Relative wave vectors to check and corresponding sinc(x/2)
161 Qrel_Q = np.pi * np.arange(5)
162 sinchalf_Q = np.array([1., 2. / np.pi, 0., - 2. / (3. * np.pi), 0.])
164 # Random parallelepipedic cell vectors to check
165 nC = 9
166 cell_Ccv = 2. * rng.random((nC, 3, 3)) - 1.
167 volume_C = np.abs(np.linalg.det(cell_Ccv))
168 # Normalize the cell volume
169 cell_Ccv /= (volume_C**(1 / 3))[:, np.newaxis, np.newaxis]
171 # Transverse wave vector components to check. Generated through the cross
172 # product of a random direction with the first cell axis.
173 v0_Cv = cell_Ccv[:, 0, :].copy()
174 v0_C = np.linalg.norm(v0_Cv, axis=-1) # Length of primary vector
175 v0n_Cv = v0_Cv / v0_C[:, np.newaxis] # Normalize
176 nd = 11
177 Q_dv = 2. * rng.random((nd, 3)) - 1.
178 Q_dv[0, :] = np.array([0., 0., 0.]) # Check also parallel Q-vector
179 Q_Cdv = np.cross(Q_dv[np.newaxis, ...], v0_Cv[:, np.newaxis, :])
181 # Volumes to test
182 nV = 7
183 Vparlp_V = 10. * rng.random(nV)
185 # ---------- Script ---------- #
187 # Rescale cell
188 cell_CVcv = cell_Ccv[:, np.newaxis, ...]\
189 * (Vparlp_V**(1 / 3.))[np.newaxis, :, np.newaxis, np.newaxis]
191 # Rescale primary vector to let Q.a follow Qrel
192 Qrel_CQ = Qrel_Q[np.newaxis, :] / v0_C[:, np.newaxis]
193 Qrel_CVQ = Qrel_CQ[:, np.newaxis, :]\
194 / (Vparlp_V**(1 / 3.))[np.newaxis, :, np.newaxis]
195 # Generate Q-vectors
196 Q_CVdQv = Qrel_CVQ[..., np.newaxis, :, np.newaxis]\
197 * v0n_Cv[:, np.newaxis, np.newaxis, np.newaxis, :]\
198 + Q_Cdv[:, np.newaxis, :, np.newaxis, :]
200 # Generate test values
201 sinchalf_CVdQ = sinc(np.sum(cell_CVcv[..., np.newaxis, np.newaxis, 1, :]
202 * Q_CVdQv, axis=-1) / 2)\
203 * sinc(np.sum(cell_CVcv[..., np.newaxis, np.newaxis, 2, :]
204 * Q_CVdQv, axis=-1) / 2)
205 test_Theta_CVdQ = Vparlp_V[np.newaxis, :, np.newaxis, np.newaxis]\
206 * sinchalf_Q[np.newaxis, np.newaxis, np.newaxis, :]\
207 * sinchalf_CVdQ
209 for Q_VdQv, test_Theta_VdQ, cell_Vcv in zip(Q_CVdQv, test_Theta_CVdQ,
210 cell_CVcv):
211 for Q_dQv, test_Theta_dQ, cell_cv in zip(Q_VdQv, test_Theta_VdQ,
212 cell_Vcv):
213 for _ in range(3): # Check that primary axis can be anywhere
214 # Slide the cell axis indices
215 cell_cv[:, :] = cell_cv[[2, 0, 1], :]
216 # Calculate geometry factors
217 Theta_dQ = parallelepipedic_geometry_factor(Q_dQv, cell_cv)
219 # Check against expected results
220 assert np.allclose(Theta_dQ, test_Theta_dQ, atol=1.e-8)
223@pytest.mark.ci
224def test_Co_hcp_site_kernels():
225 """Check that the site kernel interface works on run-time inputs."""
226 # ---------- Inputs ---------- #
228 # Part 1: Generate plane wave representation (PWDescriptor)
229 # Atomic configuration
230 a = 2.5071
231 c = 4.0695
232 mm = 1.6
233 # Ground state settings
234 xc = 'LDA'
235 kpts = 4
236 pw = 200
237 # Response settings
238 ecut = 50.
239 gammacentered = False
240 q_c = [0., 0., 0.]
241 qpm_qc = [[0., 0., 1 / 4.],
242 [0., 0., -1. / 4.]]
244 # Part 2: Calculate site kernels
245 # Define partitions to try
246 rc_pa = np.array([[1., 2.], [2., 3.]]) # radii in Å
247 hc_pa = np.array([[2., 3.], [3., 4.]]) # heights in Å
248 ez_pav = np.array([[[0., 0., 1.], [1., 0., 0.]],
249 [[0., 0., 1.], [0., 0., 1.]]])
250 cell_acv = np.array([[[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]],
251 [[1., 0., 0.], [0., 1., 0.], [0., 0., 2.]]])
252 cell_pacv = np.append(cell_acv, cell_acv * 2).reshape(2, 2, 3, 3)
254 # Part 3: Check the calculated kernels
256 # ---------- Script ---------- #
258 # Part 1: Generate plane wave representation (PWDescriptor)
259 atoms = bulk('Co', 'hcp', a=a, c=c)
260 atoms.set_initial_magnetic_moments([mm, mm])
262 calc = GPAW(xc=xc,
263 spinpol=True,
264 mode=PW(pw),
265 kpts={'size': (kpts, kpts, kpts),
266 'gamma': True}
267 )
269 # Perform inexpensive calculator initialization
270 calc.initialize(atoms)
272 qpd0 = get_pw_descriptor(atoms, calc, q_c,
273 ecut=ecut,
274 gammacentered=gammacentered)
276 # Part 2: Calculate site kernels
277 positions = atoms.get_positions()
279 # Generate spherical site kernels instances
280 # Normally
281 sph_sitekernels = SphericalSiteKernels(positions, rc_pa)
282 # Separately as sum of site kernels
283 sph_sitekernels0 = SphericalSiteKernels(positions[:1], rc_pa[:, :1])
284 sph_sitekernels1 = SphericalSiteKernels(positions[1:], rc_pa[:, 1:])
285 sph_sitekernels_sum = sph_sitekernels0 + sph_sitekernels1
286 # Seperately as appended partitions
287 sph_sitekernelsp0 = SphericalSiteKernels(positions, rc_pa[:1, :])
288 sph_sitekernelsp1 = SphericalSiteKernels(positions, rc_pa[1:, :])
289 sph_sitekernels_app = sph_sitekernelsp0.copy()
290 sph_sitekernels_app.append(sph_sitekernelsp1)
292 # Generate cylindrical site kernels instances
293 # Normally
294 cyl_sitekernels = CylindricalSiteKernels(positions, ez_pav, rc_pa, hc_pa)
295 # Separately as a sum of site kernels
296 cyl_sitekernels0 = CylindricalSiteKernels(positions[:1], ez_pav[:, :1, :],
297 rc_pa[:, :1], hc_pa[:, :1])
298 cyl_sitekernels1 = CylindricalSiteKernels(positions[1:], ez_pav[:, 1:, :],
299 rc_pa[:, 1:], hc_pa[:, 1:])
300 cyl_sitekernels_sum = cyl_sitekernels0 + cyl_sitekernels1
301 # Seperately as appended partitions
302 cyl_sitekernelsp0 = CylindricalSiteKernels(positions, ez_pav[:1, :, :],
303 rc_pa[:1, :], hc_pa[:1, :])
304 cyl_sitekernelsp1 = CylindricalSiteKernels(positions, ez_pav[1:, :, :],
305 rc_pa[1:, :], hc_pa[1:, :])
306 cyl_sitekernels_app = cyl_sitekernelsp0.copy()
307 cyl_sitekernels_app.append(cyl_sitekernelsp1)
309 # Generate parallelepipedic site kernels instances
310 # Normally
311 parlp_sitekernels = ParallelepipedicSiteKernels(positions, cell_pacv)
312 # Separately as sum of site kernels
313 parlp_sitekernels0 = ParallelepipedicSiteKernels(positions[:1],
314 cell_pacv[:, :1, ...])
315 parlp_sitekernels1 = ParallelepipedicSiteKernels(positions[1:],
316 cell_pacv[:, 1:, ...])
317 parlp_sitekernels_sum = parlp_sitekernels0 + parlp_sitekernels1
318 # Seperately as appended partitions
319 parlp_sitekernelsp0 = ParallelepipedicSiteKernels(positions,
320 cell_pacv[:1, :, ...])
321 parlp_sitekernelsp1 = ParallelepipedicSiteKernels(positions,
322 cell_pacv[1:, :, ...])
323 parlp_sitekernels_app = parlp_sitekernelsp0.copy()
324 parlp_sitekernels_app.append(parlp_sitekernelsp1)
326 # Collect all unique site kernels as both a sum and as different partitions
327 # to show, that we can compute them in parallel
328 all_sitekernels_sum = sph_sitekernels + cyl_sitekernels + parlp_sitekernels
329 all_sitekernels_app = sph_sitekernels.copy()
330 all_sitekernels_app.append(cyl_sitekernels)
331 all_sitekernels_app.append(parlp_sitekernels)
333 # Calculate spherical site kernels
334 Ksph_paGG = np.array([K_aGG for K_aGG in
335 sph_sitekernels.calculate(qpd0)])
336 Ksph0_paGG = np.array([K_aGG for K_aGG in
337 sph_sitekernels0.calculate(qpd0)])
338 Ksph1_paGG = np.array([K_aGG for K_aGG in
339 sph_sitekernels1.calculate(qpd0)])
340 Ksph_sum_paGG = np.array([K_aGG for K_aGG in
341 sph_sitekernels_sum.calculate(qpd0)])
342 Ksphp0_paGG = np.array([K_aGG for K_aGG in
343 sph_sitekernelsp0.calculate(qpd0)])
344 Ksphp1_paGG = np.array([K_aGG for K_aGG in
345 sph_sitekernelsp1.calculate(qpd0)])
346 Ksph_app_paGG = np.array([K_aGG for K_aGG in
347 sph_sitekernels_app.calculate(qpd0)])
349 # Calculate cylindrical site kernels
350 Kcyl_paGG = np.array([K_aGG for K_aGG in
351 cyl_sitekernels.calculate(qpd0)])
352 Kcyl0_paGG = np.array([K_aGG for K_aGG in
353 cyl_sitekernels0.calculate(qpd0)])
354 Kcyl1_paGG = np.array([K_aGG for K_aGG in
355 cyl_sitekernels1.calculate(qpd0)])
356 Kcyl_sum_paGG = np.array([K_aGG for K_aGG in
357 cyl_sitekernels_sum.calculate(qpd0)])
358 Kcylp0_paGG = np.array([K_aGG for K_aGG in
359 cyl_sitekernelsp0.calculate(qpd0)])
360 Kcylp1_paGG = np.array([K_aGG for K_aGG in
361 cyl_sitekernelsp1.calculate(qpd0)])
362 Kcyl_app_paGG = np.array([K_aGG for K_aGG in
363 cyl_sitekernels_app.calculate(qpd0)])
365 # Calculate parallelepipedic site kernels
366 Kparlp_paGG = np.array([K_aGG for K_aGG in
367 parlp_sitekernels.calculate(qpd0)])
368 Kparlp0_paGG = np.array([K_aGG for K_aGG in
369 parlp_sitekernels0.calculate(qpd0)])
370 Kparlp1_paGG = np.array([K_aGG for K_aGG in
371 parlp_sitekernels1.calculate(qpd0)])
372 Kparlp_sum_paGG = np.array([K_aGG for K_aGG in
373 parlp_sitekernels_sum.calculate(qpd0)])
374 Kparlpp0_paGG = np.array([K_aGG for K_aGG in
375 parlp_sitekernelsp0.calculate(qpd0)])
376 Kparlpp1_paGG = np.array([K_aGG for K_aGG in
377 parlp_sitekernelsp1.calculate(qpd0)])
378 Kparlp_app_paGG = np.array([K_aGG for K_aGG in
379 parlp_sitekernels_app.calculate(qpd0)])
381 # Calculate all site kernels together
382 Kall_sum_paGG = np.array([K_aGG for K_aGG in
383 all_sitekernels_sum.calculate(qpd0)])
384 Kall_app_paGG = np.array([K_aGG for K_aGG in
385 all_sitekernels_app.calculate(qpd0)])
387 # Calculate all site kernels at opposite qs
388 qpd_q = [get_pw_descriptor(atoms, calc, qpm_c,
389 ecut=ecut,
390 gammacentered=gammacentered)
391 for qpm_c in qpm_qc]
392 Kall_pm_qpaGG = [np.array([K_aGG for K_aGG in
393 all_sitekernels_app.calculate(qpd)])
394 for qpd in qpd_q]
396 # Part 4: Check the calculated kernels
398 # Check geometry shapes of basic arrays
399 assert all([gs == 'sphere' for gs in sph_sitekernels.geometry_shapes])
400 assert all([gs == 'cylinder' for gs in cyl_sitekernels.geometry_shapes])
401 assert all([gs == 'parallelepiped'
402 for gs in parlp_sitekernels.geometry_shapes])
404 # Check geometry shapes of summed arrays
405 assert all([gs == 'sphere' for gs in sph_sitekernels_sum.geometry_shapes])
406 assert all([gs == 'cylinder'
407 for gs in cyl_sitekernels_sum.geometry_shapes])
408 assert all([gs == 'parallelepiped'
409 for gs in parlp_sitekernels_sum.geometry_shapes])
410 assert all([gs is None for gs in all_sitekernels_sum.geometry_shapes])
412 # Check geometry shapes of appended arrays
413 assert all([gs == 'sphere' for gs in sph_sitekernels_app.geometry_shapes])
414 assert all([gs == 'cylinder'
415 for gs in cyl_sitekernels_app.geometry_shapes])
416 assert all([gs == 'parallelepiped'
417 for gs in parlp_sitekernels_app.geometry_shapes])
418 gs_refs = 2 * ['sphere'] + 2 * ['cylinder'] + 2 * ['parallelepiped']
419 assert all([gs == ref for gs, ref in
420 zip(all_sitekernels_app.geometry_shapes, gs_refs)])
422 # Check shape of spherical kernel arrays
423 nG = len(get_pw_coordinates(qpd0))
424 assert sph_sitekernels.shape == Ksph_paGG.shape[:2]
425 assert Ksph_paGG.shape == rc_pa.shape + (nG, nG)
426 assert Ksph0_paGG.shape == (rc_pa.shape[0], 1) + (nG, nG)
427 assert Ksph0_paGG.shape == Ksph1_paGG.shape
428 assert Ksph_sum_paGG.shape == rc_pa.shape + (nG, nG)
429 assert Ksphp0_paGG.shape == (1, rc_pa.shape[1]) + (nG, nG)
430 assert Ksphp0_paGG.shape == Ksphp1_paGG.shape
431 assert Ksph_app_paGG.shape == rc_pa.shape + (nG, nG)
433 # Check shape of cylindrical kernel arrays
434 assert cyl_sitekernels.shape == Kcyl_paGG.shape[:2]
435 assert Kcyl_paGG.shape == rc_pa.shape + (nG, nG)
436 assert Kcyl0_paGG.shape == (rc_pa.shape[0], 1) + (nG, nG)
437 assert Kcyl0_paGG.shape == Kcyl1_paGG.shape
438 assert Kcyl_sum_paGG.shape == rc_pa.shape + (nG, nG)
439 assert Kcylp0_paGG.shape == (1, rc_pa.shape[1]) + (nG, nG)
440 assert Kcylp0_paGG.shape == Kcylp1_paGG.shape
441 assert Kcyl_app_paGG.shape == rc_pa.shape + (nG, nG)
443 # Check shape of parallelepipedic kernel arrays
444 assert parlp_sitekernels.shape == Kparlp_paGG.shape[:2]
445 assert Kparlp_paGG.shape == cell_pacv.shape[:2] + (nG, nG)
446 assert Kparlp0_paGG.shape == (cell_pacv.shape[0], 1) + (nG, nG)
447 assert Kparlp0_paGG.shape == Kparlp1_paGG.shape
448 assert Kparlp_sum_paGG.shape == cell_pacv.shape[:2] + (nG, nG)
449 assert Kparlpp0_paGG.shape == (1, cell_pacv.shape[1]) + (nG, nG)
450 assert Kparlpp0_paGG.shape == Kparlpp1_paGG.shape
451 assert Kparlp_app_paGG.shape == cell_pacv.shape[:2] + (nG, nG)
453 # Check shape of array calculated in parallel
454 assert Kall_sum_paGG.shape == (2, 6, nG, nG)
455 assert Kall_app_paGG.shape == (6, 2, nG, nG)
457 # Check self-consitency of spherical arrays
458 assert np.allclose(Ksph0_paGG[:, 0, ...], Ksph_sum_paGG[:, 0, ...])
459 assert np.allclose(Ksph1_paGG[:, 0, ...], Ksph_sum_paGG[:, 1, ...])
460 assert np.allclose(Ksph_paGG, Ksph_sum_paGG)
461 assert np.allclose(Ksphp0_paGG[0], Ksph_app_paGG[0])
462 assert np.allclose(Ksphp1_paGG[0], Ksph_app_paGG[1])
463 assert np.allclose(Ksph_paGG, Ksph_app_paGG)
465 # Check self-consistency of cylindrical arrays
466 assert np.allclose(Kcyl0_paGG[:, 0, ...], Kcyl_sum_paGG[:, 0, ...])
467 assert np.allclose(Kcyl1_paGG[:, 0, ...], Kcyl_sum_paGG[:, 1, ...])
468 assert np.allclose(Kcyl_paGG, Kcyl_sum_paGG)
469 assert np.allclose(Kcylp0_paGG[0], Kcyl_app_paGG[0])
470 assert np.allclose(Kcylp1_paGG[0], Kcyl_app_paGG[1])
471 assert np.allclose(Kcyl_paGG, Kcyl_app_paGG)
473 # Check self-consistency of parallelepipedic arrays
474 assert np.allclose(Kparlp0_paGG[:, 0, ...], Kparlp_sum_paGG[:, 0, ...])
475 assert np.allclose(Kparlp1_paGG[:, 0, ...], Kparlp_sum_paGG[:, 1, ...])
476 assert np.allclose(Kparlp_paGG, Kparlp_sum_paGG)
477 assert np.allclose(Kparlpp0_paGG[0], Kparlp_app_paGG[0])
478 assert np.allclose(Kparlpp1_paGG[0], Kparlp_app_paGG[1])
479 assert np.allclose(Kparlp_paGG, Kparlp_app_paGG)
481 # Check self-consistency of kernels calculated in parallel
482 assert np.allclose(Ksph_paGG, Kall_sum_paGG[:, :2])
483 assert np.allclose(Kcyl_paGG, Kall_sum_paGG[:, 2:4])
484 assert np.allclose(Kparlp_paGG, Kall_sum_paGG[:, -2:])
485 assert np.allclose(Ksph_paGG, Kall_app_paGG[:2])
486 assert np.allclose(Kcyl_paGG, Kall_app_paGG[2:4])
487 assert np.allclose(Kparlp_paGG, Kall_app_paGG[-2:])
489 # Check that K_G=G'(q=0) gives Vint / V0 (fractional integration volume)
490 # Volume of unit cell in Å^3
491 V0 = atoms.get_volume()
492 # Calculate integration volumes in Å^3
493 Vsphere_pa = 4 / 3 * np.pi * rc_pa**3
494 Vcylinder_pa = np.pi * rc_pa**2 * hc_pa
495 Vparlp_pa = np.abs(np.linalg.det(cell_pacv))
496 assert np.allclose(np.diagonal(Ksph_paGG, axis1=2, axis2=3),
497 Vsphere_pa[..., np.newaxis] / V0)
498 assert np.allclose(np.diagonal(Kcyl_paGG, axis1=2, axis2=3),
499 Vcylinder_pa[..., np.newaxis] / V0)
500 assert np.allclose(np.diagonal(Kparlp_paGG, axis1=2, axis2=3),
501 Vparlp_pa[..., np.newaxis] / V0)
503 # Check that K_G=G'(q) gives e^(-iq.τ_a) Θ(q) / V0
504 B_cv = 2.0 * np.pi * atoms.cell.reciprocal() # Coordinate transform
505 for qpm_c, Kall_pm_paGG in zip(qpm_qc, Kall_pm_qpaGG):
506 # Calculate q-vector in Å^(-1)
507 qpm_v = qpm_c @ B_cv
508 exp_a = np.exp(-1.j * positions @ qpm_v)
510 # Calculate site centered geometry factors
511 Theta_pa = []
512 # Spherical geometry factors
513 for rc_a in rc_pa:
514 Theta_a = []
515 for rc in rc_a:
516 Theta_a.append(spherical_geometry_factor(qpm_v, rc))
517 Theta_pa.append(Theta_a)
518 # Cylindrical geometry factors
519 for ez_av, rc_a, hc_a in zip(ez_pav, rc_pa, hc_pa):
520 Theta_a = []
521 for ez_v, rc, hc in zip(ez_av, rc_a, hc_a):
522 Theta_a.append(cylindrical_geometry_factor(qpm_v, ez_v,
523 rc, hc))
524 Theta_pa.append(Theta_a)
525 # Parallelepipedic geometry factors
526 for cell_acv in cell_pacv:
527 Theta_a = []
528 for cell_cv in cell_acv:
529 Theta_a.append(parallelepipedic_geometry_factor(qpm_v,
530 cell_cv))
531 Theta_pa.append(Theta_a)
532 Theta_pa = np.array(Theta_pa)
534 assert np.allclose(np.diagonal(Kall_pm_paGG, axis1=2, axis2=3),
535 Theta_pa[..., np.newaxis]
536 * exp_a[np.newaxis, :, np.newaxis] / V0)
538 # Check that K_GG'(q) is the hermitian conjugate of K_GG'(-q)
539 Kall_mH_paGG = np.conjugate(np.transpose(Kall_pm_qpaGG[1], (0, 1, 3, 2)))
540 assert np.allclose(Kall_pm_qpaGG[0], Kall_mH_paGG)
543# ---------- Test functionality ---------- #
546def get_pw_descriptor(atoms, calc, q_c, ecut=50., gammacentered=False):
547 """Mock-up of ChiKSCalculator.get_pw_descriptor.
549 Works on a bare calculator instance without any actual data in it."""
550 from ase.units import Ha
551 from gpaw.response.qpd import SingleQPWDescriptor
553 # Create the plane wave descriptor
554 q_c = np.asarray(q_c, dtype=float)
555 qpd = SingleQPWDescriptor.from_q(q_c, ecut / Ha, calc.wfs.gd,
556 gammacentered=gammacentered)
557 return qpd