Coverage for gpaw/test/response/test_heisenberg.py: 100%
92 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-09 00:21 +0000
1"""Test the Heisenberg model based methodology of the response code."""
3import numpy as np
4import pytest
6from gpaw.response.heisenberg import (calculate_fm_magnon_energies,
7 calculate_single_site_magnon_energies)
10@pytest.mark.ci
11def test_single_site_magnons():
12 """Check the single site magnon dispersion functionality."""
13 rng = get_rng()
14 # ---------- Inputs ---------- #
16 # Magnetic moment
17 mm = 1.
18 # q-point grid
19 nq = 11
20 q_qc = get_randomized_qpoints(nq, rng)
22 # Random J_q, with J=0 at q=0
23 J_q = rng.rand(q_qc.shape[0])
24 J_q[list(q_qc[:, 2]).index(0.)] = 0.
26 # Cosine J_qD with different spin wave stiffnesses D
27 D_D = np.linspace(400., 800., 5)
28 J_qD = np.outer(np.cos(q_qc[:, 2]), D_D)
30 # ---------- Script ---------- #
32 # Calculate magnon energies
33 E_q = calculate_single_site_magnon_energies(J_q, q_qc, mm)
34 E_qD = calculate_single_site_magnon_energies(J_qD, q_qc, mm)
36 # Check dimensions of arrays
37 assert E_q.shape == (q_qc.shape[0],)
38 assert E_qD.shape == J_qD.shape
40 # Check versus formulas
41 assert np.allclose(E_q, -2. / mm * J_q) # Remember: J(0) = 0
42 assert np.allclose(E_qD, 2. / mm * D_D[np.newaxis, :]
43 * (1. - np.cos(q_qc[:, 2]))[:, np.newaxis])
46@pytest.mark.ci
47def test_single_site_magnons_consistency():
48 """Check that the generalized magnon dispersion calculation is consistent
49 for a single site system with the simple analytical formula valid in that
50 case."""
51 rng = get_rng()
52 # ---------- Inputs ---------- #
54 # Magnetic moment
55 mm = 1.
56 # q-point grid
57 nq = 11
58 q_qc = get_randomized_qpoints(nq, rng)
60 # Random isotropic exchange constants
61 nJsamples = 6 # sample some different random combinations
62 J_qx = rng.rand(q_qc.shape[0], nJsamples)
64 # ---------- Script ---------- #
66 # Calculate assuming a single site
67 E_qx = calculate_single_site_magnon_energies(J_qx, q_qc, mm)
69 # Calcualte using generalized functionality
70 E_qnx = calculate_fm_magnon_energies(J_qx[:, np.newaxis, np.newaxis, :],
71 q_qc, mm * np.ones((1, nJsamples)))
73 # Test self-consistency
74 assert E_qnx.shape[0] == E_qx.shape[0]
75 assert E_qnx.shape[1] == 1
76 assert E_qnx.shape[-1] == E_qx.shape[-1]
77 assert np.allclose(E_qnx[:, 0, :], E_qx, atol=1e-8)
80@pytest.mark.ci
81def test_fm_random_magnons():
82 """Check that the functionality to calculate the magnon dispersion of a
83 ferromagnetic system with multiple sites works for a randomized system with
84 three sites."""
85 rng = get_rng()
86 # ---------- Inputs ---------- #
88 # Magnetic moments
89 nsites = 3
90 mm_a = 5. * rng.rand(nsites)
91 # q-point grid
92 nq = 11
93 q_qc = get_randomized_qpoints(nq, rng)
95 # Random isotropic exchange constants
96 J_qab = 1.j * rng.rand(q_qc.shape[0], nsites, nsites)
97 J_qab += rng.rand(q_qc.shape[0], nsites, nsites)
98 # Take the Hermitian part of random tensor
99 J_qab = (J_qab + np.transpose(np.conjugate(J_qab), (0, 2, 1))) / 2.
100 # The q=0 component should furthermore be real
101 J_qab[list(q_qc[:, 2]).index(0.)].imag = 0.
103 # ---------- Script ---------- #
105 # Calculate magnon energies
106 E_qn = calculate_fm_magnon_energies(J_qab, q_qc, mm_a)
107 E_qn = np.sort(E_qn, axis=1) # Make sure the eigenvalues are sorted
109 # Calculate the magnon energies manually
110 mm_inv_ab = 2. / np.sqrt(np.outer(mm_a, mm_a))
111 J0_ab = np.diag(np.sum(J_qab[list(q_qc[:, 2]).index(0.)], axis=1))
112 H_qab = mm_inv_ab[np.newaxis] * (J0_ab[np.newaxis] - J_qab)
113 test_E_qn, _ = np.linalg.eig(H_qab)
115 assert E_qn.shape == (q_qc.shape[0], nsites)
116 assert np.allclose(test_E_qn.imag, 0.)
117 assert np.allclose(E_qn, np.sort(test_E_qn.real, axis=1))
120@pytest.mark.ci
121def test_fm_vectorized_magnons():
122 """Check that the functionality to calculate the magnon dispersion of a
123 ferromagnetic system with multiple sites works when supplying multiple
124 sets of parameters for the same two-site systems."""
125 rng = get_rng()
126 # ---------- Inputs ---------- #
128 # Magnetic moments
129 nsites = 2
130 nmagmoms = 4 # Test the same J_qab, but with different site magnetizations
131 mm_ax = 5. * rng.rand(nsites, nmagmoms)
132 # q-point grid
133 nq = 11
134 q_qc = get_randomized_qpoints(nq, rng)
136 # Use a fixed structure for J_qab with known eigenvalues
137 cos_q = np.cos(q_qc[:, 2])
138 sin_q = np.sin(q_qc[:, 2])
139 J_qab = np.empty((nq, nsites, nsites), dtype=complex)
140 J_qab[:, 0, 0] = cos_q
141 J_qab[:, 0, 1] = 1. + 1.j * sin_q
142 J_qab[:, 1, 0] = 1. - 1.j * sin_q
143 J_qab[:, 1, 1] = 2. * cos_q
145 # Test different energy scales for the exchange interactions
146 nJscales = 6
147 Jscale_y = 800. * rng.rand(nJscales)
149 # Combine different magnetic moments and scale for the exchange
150 J_qabxy = np.empty(J_qab.shape + (nmagmoms, nJscales,), dtype=complex)
151 J_qabxy[:] = np.tensordot(J_qab, Jscale_y,
152 axes=((), ()))[..., np.newaxis, :]
153 mm_axy = np.moveaxis(np.tile(mm_ax, (nJscales, 1, 1)), 0, -1)
155 # ---------- Script ---------- #
157 # Calculate magnon energies
158 E_qnxy = calculate_fm_magnon_energies(J_qabxy, q_qc, mm_axy)
159 E_qnxy = np.sort(E_qnxy, axis=1) # Make sure the eigenvalues are sorted
161 # Calculate magnon energies analytically
162 H_diag1_qxy = np.sqrt(mm_axy[1][np.newaxis] / mm_axy[0][np.newaxis])\
163 * (2. - cos_q[:, np.newaxis, np.newaxis])
164 H_diag2_qxy = np.sqrt(mm_axy[0][np.newaxis] / mm_axy[1][np.newaxis])\
165 * (3. - 2. * cos_q[:, np.newaxis, np.newaxis])
166 H_diag_avg_qxy = (H_diag1_qxy + H_diag2_qxy) / 2.
167 H_diag_diff_qxy = (H_diag1_qxy - H_diag2_qxy) / 2.
168 pm_n = np.array([-1., 1.])
169 E_test_qnxy = H_diag_avg_qxy[:, np.newaxis]\
170 + pm_n[np.newaxis, :, np.newaxis, np.newaxis]\
171 * np.sqrt(H_diag_diff_qxy[:, np.newaxis]**2.
172 + (1 + sin_q[:, np.newaxis, np.newaxis, np.newaxis]**2.))
173 E_test_qnxy *= 2. / np.sqrt(np.prod(mm_axy, axis=0))
174 E_test_qnxy *= Jscale_y
176 assert np.allclose(E_qnxy, E_test_qnxy)
179# ---------- Test functionality ---------- #
182def get_randomized_qpoints(nq, rng):
183 """Make a simple, but shuffled, q-point array."""
184 q_qc = np.zeros((nq, 3), dtype=float)
185 q_qc[:, 2] = np.linspace(0., np.pi, nq)
186 rng.shuffle(q_qc[:, 2])
188 return q_qc
191def get_rng():
192 """Choose a specific random seed to make the tests reproducible."""
193 rng = np.random.RandomState(23)
195 return rng