Coverage for gpaw/response/site_kernels.py: 99%
152 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-09 00:21 +0000
1"""Compute site-kernels. Used for computing Heisenberg exchange constants.
2Specifically, one maps a DFT calculations onto a Heisenberg lattice model,
3where the site kernels define the lattice sites and magnetic moments."""
5import numpy as np
6from scipy.special import jv
7from gpaw.response.pair_functions import get_pw_coordinates
8from ase.units import Bohr
11class SiteKernels:
12 """Factory for calculating sublattice site kernels
14 1 /
15 K_aGG'(q) = ‾‾ | dr e^(-i[G-G'+q].r) θ(r∊V_a)
16 V0 /
18 where V_a denotes the integration volume of site a, centered at the site
19 position τ_a, and V0 denotes the unit cell volume."""
21 def __init__(self, positions, partitions):
22 """Construct the site kernel factory.
24 Parameters
25 ----------
26 positions : np.ndarray
27 Site positions in a.u. (Bohr). Shape: (nsites, 3).
28 partitions : list
29 List (len=npartitions) of lists (len=nsites) with site geometries
30 where the geometry arguments are given in atomic units.
31 For more information on site geometries, see calculate_site_kernels
32 """
33 assert isinstance(partitions, list)
34 assert all([isinstance(geometries, list)
35 and len(geometries) == positions.shape[0]
36 for geometries in partitions])
37 assert positions.shape[1] == 3
39 self.positions = positions
40 self.partitions = partitions
42 @property
43 def npartitions(self):
44 return len(self.partitions)
46 @property
47 def nsites(self):
48 return self.positions.shape[0]
50 @property
51 def shape(self):
52 return self.npartitions, self.nsites
54 @property
55 def geometry_shapes(self):
56 """If all sites of a given partition has the same geometry, the
57 partition is said to have a geometry shape. Otherwise, the
58 geometry shape is None."""
60 geometry_shapes = []
61 for geometries in self.partitions:
62 # Record the geometry shape, if they are all the same
63 if all([shape == geometries[0][0] for shape, _ in geometries]):
64 geometry_shapes.append(geometries[0][0])
65 else:
66 geometry_shapes.append(None)
68 return geometry_shapes
70 def calculate(self, qpd):
71 """Generate the site kernels of each partition.
73 Returns
74 -------
75 K_aGG : np.ndarray (dtype=complex)
76 Site kernels of sites a and plane wave components G and G'.
77 """
78 for geometries in self.partitions:
79 # We yield one set of site kernels at a time, because they can be
80 # memory intensive
81 yield calculate_site_kernels(qpd, self.positions, geometries)
83 def __add__(self, sitekernels):
84 """Add the sites from two SiteKernels instances to a new joint
85 SiteKernels instance with nsites = nsites1 + nsites2."""
86 assert isinstance(sitekernels, SiteKernels)
87 assert self.npartitions == sitekernels.npartitions
89 # Join positions
90 nsites = self.nsites + sitekernels.nsites
91 positions = np.append(self.positions,
92 sitekernels.positions).reshape(nsites, 3)
94 # Join partitions
95 partitions = [geometries1 + geometries2
96 for geometries1, geometries2
97 in zip(self.partitions, sitekernels.partitions)]
99 return SiteKernels(positions, partitions)
101 def append(self, sitekernels):
102 """Append the partitions of another array with identical site
103 positions such that npartitions = npartitions1 + npartitions2."""
104 assert isinstance(sitekernels, SiteKernels)
105 assert self.nsites == sitekernels.nsites
106 assert np.allclose(self.positions, sitekernels.positions)
108 self.partitions += sitekernels.partitions
110 def copy(self):
111 return SiteKernels(self.positions.copy(), self.partitions.copy())
114class SphericalSiteKernels(SiteKernels):
116 def __init__(self, positions, radii):
117 """Construct a site kernel factory with spherical kernels.
119 Parameters
120 ----------
121 positions : np.ndarray
122 Site positions in Angstrom (Å). Shape: (nsites, 3).
123 radii : list or np.ndarray
124 Spherical radii of the sites in Angstrom (Å).
125 Shape: (npartitions, nsites), where the individual spherical radii
126 can be varried for each spatial partitioning.
127 """
128 positions = np.asarray(positions)
130 # Parse the input spherical radii
131 rc_pa = np.asarray(radii)
132 assert len(rc_pa.shape) == 2
133 assert rc_pa.shape[1] == positions.shape[0]
135 # Convert radii to internal units (Å to Bohr)
136 positions = positions / Bohr
137 rc_pa = rc_pa / Bohr
139 # Generate partitions as list of lists of geometries
140 partitions = [[('sphere', (rc,)) for rc in rc_a]
141 for rc_a in rc_pa]
143 SiteKernels.__init__(self, positions, partitions)
146class CylindricalSiteKernels(SiteKernels):
148 def __init__(self, positions, directions, radii, heights):
149 """Construct a site kernel factory with cylindrical kernels.
151 Parameters
152 ----------
153 positions : np.ndarray
154 Site positions in Angstrom (Å). Shape: (nsites, 3).
155 directions : np.ndarray
156 Normalized directions of the cylindrical axes.
157 Shape: (npartitions, nsites, 3), where the direction of each
158 individual cylinder can be varried (along with the radius and
159 height) for each spatial partitioning.
160 radii : np.ndarray
161 Cylindrical radii of the sites in Angstrom (Å).
162 Shape: (npartitions, nsites).
163 heights : list or np.ndarray
164 Cylinder heights in Angstrom (Å). Shape: (npartitions, nsites).
165 """
166 positions = np.asarray(positions)
168 # Parse the cylinder geometry arguments
169 ez_pav = np.asarray(directions)
170 rc_pa = np.asarray(radii)
171 hc_pa = np.asarray(heights)
172 nsites = positions.shape[0]
173 npartitions = ez_pav.shape[0]
174 assert ez_pav.shape == (npartitions, nsites, 3)
175 assert np.allclose(np.linalg.norm(ez_pav, axis=-1), 1., atol=1.e-8)
176 assert rc_pa.shape == (npartitions, nsites)
177 assert hc_pa.shape == (npartitions, nsites)
179 # Convert to internal units (Å to Bohr)
180 positions = positions / Bohr
181 rc_pa = rc_pa / Bohr
182 hc_pa = hc_pa / Bohr
184 # Generate partitions as list of lists of geometries
185 partitions = [[('cylinder', (ez_v, rc, hc))
186 for ez_v, rc, hc in zip(ez_av, rc_a, hc_a)]
187 for ez_av, rc_a, hc_a in zip(ez_pav, rc_pa, hc_pa)]
189 SiteKernels.__init__(self, positions, partitions)
192class ParallelepipedicSiteKernels(SiteKernels):
194 def __init__(self, positions, cells):
195 """Construct a site kernel factory with parallelepipedic kernels.
197 Parameters
198 ----------
199 positions : np.ndarray
200 Site positions in Angstrom (Å). Shape: (nsites, 3).
201 cells : np.ndarray
202 Cell vectors of the parallelepiped in Angstrom (Å).
203 Shape: (npartitions, nsites, 3, 3), where the parallelepipedic
204 cell of each site can be varried independently for each spatial
205 partitioning. The second to last entry is the vector index and
206 the last entry indexes the cartesian components.
207 """
208 positions = np.asarray(positions)
210 # Parse the parallelepipeds' cells
211 cell_pacv = np.asarray(cells)
212 assert len(cell_pacv.shape) == 4
213 assert cell_pacv.shape[1:] == (positions.shape[0], 3, 3)
215 # Convert to internal units (Å to Bohr)
216 positions = positions / Bohr
217 cell_pacv = cell_pacv / Bohr
219 # Generate partitions as list of lists of geometries
220 partitions = [[('parallelepiped', (cell_cv,)) for cell_cv in cell_acv]
221 for cell_acv in cell_pacv]
223 SiteKernels.__init__(self, positions, partitions)
226def calculate_site_kernels(qpd, positions, geometries):
227 """Calculate the sublattice site kernel:
229 1 /
230 K_aGG'(q) = ‾‾ | dr e^(-i[G-G'+q].r) θ(r∊V_a)
231 V0 /
233 where V_a denotes the integration volume of site a, centered at the site
234 position τ_a, and V0 denotes the unit cell volume.
236 In the calculation, the kernel is split in two contributions:
238 1) The Fourier component of the site position:
240 τ_a(Q) = e^(-iQ.τ_a)
242 2) The site centered geometry factor
244 /
245 Θ(Q) = | dr e^(-iQ.r) θ(r+τ_a∊V_a)
246 /
248 where Θ(Q) only depends on the geometry of the integration volume.
249 With this:
251 1
252 K_aGG'(q) = ‾‾ τ_a(G-G'+q) Θ(G-G'+q)
253 V0
255 Parameters
256 ----------
257 qpd : SingleQPWDescriptor
258 Plane wave descriptor corresponding to the q wave vector of interest.
259 positions : np.ndarray
260 Site positions. Array of shape (nsites, 3).
261 geometries : list
262 List of site geometries. A site geometry is a tuple of the integration
263 volume shape (str) and arguments (tuple): (shape, args). Valid shapes
264 are 'sphere', 'cylinder' and 'parallelepiped'. The integration volume
265 arguments specify the size and orientation of the integration region.
267 Returns
268 -------
269 K_aGG : np.ndarray (dtype=complex)
270 Site kernels of sites a and plane wave components G and G'.
271 """
272 assert positions.shape[0] == len(geometries)
273 assert positions.shape[1] == 3
275 # Extract unit cell volume
276 V0 = qpd.gd.volume
278 # Construct Q=G-G'+q
279 Q_GGv = construct_wave_vectors(qpd)
281 # Allocate site kernel array
282 nsites = len(geometries)
283 K_aGG = np.zeros((nsites,) + Q_GGv.shape[:2], dtype=complex)
285 # Calculate the site kernel for each site individually
286 for a, (tau_v, (shape, args)) in enumerate(zip(positions, geometries)):
288 # Compute the site centered geometry factor
289 _geometry_factor = create_geometry_factor(shape) # factory pattern
290 Theta_GG = _geometry_factor(Q_GGv, *args)
292 # Compute the Fourier component of the site position
293 tau_GG = np.exp(-1.j * Q_GGv @ tau_v)
295 # Update data
296 K_aGG[a, :, :] = 1 / V0 * tau_GG * Theta_GG
298 return K_aGG
301def construct_wave_vectors(qpd):
302 """Construct wave vectors Q=G1-G2+q corresponding to the q-vector of
303 interest."""
304 G_Gv, q_v = get_plane_waves_and_reduced_wave_vector(qpd)
306 # Allocate arrays for G, G' and q respectively
307 nG = len(G_Gv)
308 G1_GGv = np.tile(G_Gv[:, np.newaxis, :], [1, nG, 1])
309 G2_GGv = np.tile(G_Gv[np.newaxis, :, :], [nG, 1, 1])
310 q_GGv = np.tile(q_v[np.newaxis, np.newaxis, :], [nG, nG, 1])
312 # Contruct the wave vector G1 - G2 + q
313 Q_GGv = G1_GGv - G2_GGv + q_GGv
315 return Q_GGv
318def get_plane_waves_and_reduced_wave_vector(qpd):
319 """Get the reciprocal lattice vectors and reduced wave vector of the plane
320 wave representation corresponding to the q-vector of interest."""
321 # Get the reduced wave vector
322 q_c = qpd.q_c
324 # Get the reciprocal lattice vectors in relative coordinates
325 G_Gc = get_pw_coordinates(qpd)
327 # Convert to cartesian coordinates
328 B_cv = 2.0 * np.pi * qpd.gd.icell_cv # Coordinate transform matrix
329 q_v = np.dot(q_c, B_cv) # Unit = Bohr^(-1)
330 G_Gv = np.dot(G_Gc, B_cv)
332 return G_Gv, q_v
335def create_geometry_factor(shape):
336 """Creator component of the geometry factor factory pattern."""
337 if shape == 'sphere':
338 return spherical_geometry_factor
339 elif shape == 'cylinder':
340 return cylindrical_geometry_factor
341 elif shape == 'parallelepiped':
342 return parallelepipedic_geometry_factor
344 raise ValueError('Invalid site kernel shape:', shape)
347def spherical_geometry_factor(Q_Qv, rc):
348 """Calculate the site centered geometry factor for a spherical site kernel:
350 /
351 Θ(Q) = | dr e^(-iQ.r) θ(|r|<r_c)
352 /
354 4πr_c
355 = ‾‾‾‾‾ [sinc(|Q|r_c) - cos(|Q|r_c)]
356 |Q|^2
358 3 [sinc(|Q|r_c) - cos(|Q|r_c)]
359 = V_sphere ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾
360 (|Q|r_c)^2
362 where the dimensionless geometry factor satisfies:
364 Θ(Q)/V_sphere --> 1 for |Q|r_c --> 0.
366 Parameters
367 ----------
368 Q_Qv : np.ndarray
369 Wave vectors to evaluate the site centered geometry factor at. The
370 cartesian coordinates needs to be the last dimension of the array (v),
371 but the preceeding index/indices Q can have any tensor structure, such
372 that Q_Qv.shape = (..., 3).
373 rc : float
374 Radius of the sphere.
375 """
376 assert Q_Qv.shape[-1] == 3
377 assert isinstance(rc, float) and rc > 0.
379 # Calculate the sphere volume
380 Vsphere = 4 * np.pi * rc**3. / 3
382 # Calculate |Q|r_c
383 Qrc_Q = np.linalg.norm(Q_Qv, axis=-1) * rc
385 # Allocate array with ones to provide the correct dimensionless geometry
386 # factor in the |Q|r_c --> 0 limit.
387 # This is done to avoid division by zero.
388 Theta_Q = np.ones(Q_Qv.shape[:-1], dtype=float)
390 # Calculate the dimensionless geometry factor
391 Qrcs = Qrc_Q[Qrc_Q > 1.e-8]
392 Theta_Q[Qrc_Q > 1.e-8] = 3. * (sinc(Qrcs) - np.cos(Qrcs)) / Qrcs**2.
394 Theta_Q *= Vsphere
396 return Theta_Q
399def cylindrical_geometry_factor(Q_Qv, ez_v, rc, hc):
400 """Calculate site centered geometry factor for a cylindrical site kernel:
402 /
403 Θ(Q) = | dr e^(-iQ.r) θ(ρ<r_c) θ(|z|/2<h_c)
404 /
406 4πr_c
407 = ‾‾‾‾‾‾‾ J_1(Q_ρ r_c) sin(Q_z h_c / 2)
408 Q_ρ Q_z
410 2 J_1(Q_ρ r_c)
411 = V_cylinder ‾‾‾‾‾‾‾‾‾‾‾‾‾‾ sinc(Q_z h_c / 2)
412 Q_ρ r_c
414 where z denotes the cylindrical axis, ρ the radial axis and the
415 dimensionless geometry factor satisfy:
417 Θ(Q)/V_cylinder --> 1 for Q --> 0.
419 Parameters
420 ----------
421 Q_Qv : np.ndarray
422 Wave vectors to evaluate the site centered geometry factor at. The
423 cartesian coordinates needs to be the last dimension of the array (v),
424 but the preceeding index/indices Q can have any tensor structure, such
425 that Q_Qv.shape = (..., 3).
426 ez_v : np.ndarray
427 Normalized direction of the cylindrical axis.
428 rc : float
429 Radius of the cylinder.
430 hc : float
431 Height of the cylinder.
432 """
433 assert Q_Qv.shape[-1] == 3
434 assert ez_v.shape == (3,)
435 assert abs(np.linalg.norm(ez_v) - 1.) < 1.e-8
436 assert isinstance(rc, float) and rc > 0.
437 assert isinstance(hc, float) and hc > 0.
439 # Calculate cylinder volume
440 Vcylinder = np.pi * rc**2. * hc
442 # Calculate Q_z h_c and Q_ρ r_c
443 Qzhchalf_Q = np.abs(Q_Qv @ ez_v) * hc / 2.
444 Qrhorc_Q = np.linalg.norm(np.cross(Q_Qv, ez_v), axis=-1) * rc
446 # Allocate array with ones to provide the correct dimensionless geometry
447 # factor in the Q_ρ r_c --> 0 limit.
448 # This is done to avoid division by zero.
449 Theta_Q = np.ones(Q_Qv.shape[:-1], dtype=float)
451 # Calculate the dimensionless geometry factor
452 Qrhorcs = Qrhorc_Q[Qrhorc_Q > 1.e-8]
453 Theta_Q[Qrhorc_Q > 1.e-8] = 2. * jv(1, Qrhorcs) / Qrhorcs
454 Theta_Q *= Vcylinder * sinc(Qzhchalf_Q)
456 return Theta_Q
459def parallelepipedic_geometry_factor(Q_Qv, cell_cv):
460 """Calculate the site centered geometry factor for a parallelepipedic site
461 kernel:
463 /
464 Θ(Q) = | dr e^(-iQ.r) θ(r∊V_parallelepiped)
465 /
467 = |det[a1, a2, a3]| sinc(Q.a1 / 2) sinc(Q.a2 / 2) sinc(Q.a3 / 2)
469 = V_parallelepiped sinc(Q.a1 / 2) sinc(Q.a2 / 2) sinc(Q.a3 / 2)
471 where a1, a2 and a3 denotes the parallelepipedic cell vectors.
473 Parameters
474 ----------
475 Q_Qv : np.ndarray
476 Wave vectors to evaluate the site centered geometry factor at. The
477 cartesian coordinates needs to be the last dimension of the array (v),
478 but the preceeding index/indices Q can have any tensor structure, such
479 that Q_Qv.shape = (..., 3).
480 cell_cv : np.ndarray, shape=(3, 3)
481 Cell vectors of the parallelepiped, where v denotes the cartesian
482 coordinates.
483 """
484 assert Q_Qv.shape[-1] == 3
485 assert cell_cv.shape == (3, 3)
487 # Calculate the parallelepiped volume
488 Vparlp = abs(np.linalg.det(cell_cv))
489 assert Vparlp > 1.e-8 # Not a valid parallelepiped if volume vanishes
491 # Calculate the site-kernel
492 a1, a2, a3 = cell_cv
493 Theta_Q = Vparlp * sinc(Q_Qv @ a1 / 2) * sinc(Q_Qv @ a2 / 2) * \
494 sinc(Q_Qv @ a3 / 2)
496 return Theta_Q
499def sinc(x):
500 """np.sinc(x) = sin(pi*x) / (pi*x), hence the division by pi"""
501 return np.sinc(x / np.pi)