Coverage for /builds/ase/ase/ase/build/general_surface.py: 83.61%
61 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3from math import gcd
5import numpy as np
6from numpy.linalg import norm, solve
8from ase.build import bulk
9from ase.build.surface import create_tags
12def surface(lattice, indices, layers, vacuum=None, tol=1e-10, periodic=False):
13 """Create surface from a given lattice and Miller indices.
15 lattice: Atoms object or str
16 Bulk lattice structure of alloy or pure metal. Note that the
17 unit-cell must be the conventional cell - not the primitive cell.
18 One can also give the chemical symbol as a string, in which case the
19 correct bulk lattice will be generated automatically.
20 indices: sequence of three int
21 Surface normal in Miller indices (h,k,l).
22 layers: int
23 Number of equivalent layers of the slab.
24 vacuum: float
25 Amount of vacuum added on both sides of the slab.
26 periodic: bool
27 Whether the surface is periodic in the normal to the surface
28 """
30 indices = np.asarray(indices)
32 if indices.shape != (3,) or not indices.any() or indices.dtype != int:
33 raise ValueError(f'{indices} is an invalid surface type')
35 if isinstance(lattice, str):
36 lattice = bulk(lattice, cubic=True)
38 h, k, l = indices # noqa (E741, the variable l)
39 h0, k0, l0 = (indices == 0)
41 if h0 and k0 or h0 and l0 or k0 and l0: # if two indices are zero
42 if not h0:
43 c1, c2, c3 = [(0, 1, 0), (0, 0, 1), (1, 0, 0)]
44 if not k0:
45 c1, c2, c3 = [(0, 0, 1), (1, 0, 0), (0, 1, 0)]
46 if not l0:
47 c1, c2, c3 = [(1, 0, 0), (0, 1, 0), (0, 0, 1)]
48 else:
49 p, q = ext_gcd(k, l)
50 a1, a2, a3 = lattice.cell
52 # constants describing the dot product of basis c1 and c2:
53 # dot(c1,c2) = k1+i*k2, i in Z
54 k1 = np.dot(p * (k * a1 - h * a2) + q * (l * a1 - h * a3),
55 l * a2 - k * a3)
56 k2 = np.dot(l * (k * a1 - h * a2) - k * (l * a1 - h * a3),
57 l * a2 - k * a3)
59 if abs(k2) > tol:
60 i = -int(round(k1 / k2)) # i corresponding to the optimal basis
61 p, q = p + i * l, q - i * k
63 a, b = ext_gcd(p * k + q * l, h)
65 c1 = (p * k + q * l, -p * h, -q * h)
66 c2 = np.array((0, l, -k)) // abs(gcd(l, k))
67 c3 = (b, a * p, a * q)
69 surf = build(lattice, np.array([c1, c2, c3]), layers, tol, periodic)
70 if vacuum is not None:
71 surf.center(vacuum=vacuum, axis=2)
72 return surf
75def build(lattice, basis, layers, tol, periodic):
76 surf = lattice.copy()
77 scaled = solve(basis.T, surf.get_scaled_positions().T).T
78 scaled -= np.floor(scaled + tol)
79 surf.set_scaled_positions(scaled)
80 surf.set_cell(np.dot(basis, surf.cell), scale_atoms=True)
81 surf *= (1, 1, layers)
82 surf.set_tags(create_tags((1, len(lattice), layers)))
84 a1, a2, a3 = surf.cell
85 surf.set_cell([a1, a2,
86 np.cross(a1, a2) * np.dot(a3, np.cross(a1, a2)) /
87 norm(np.cross(a1, a2))**2])
89 # Change unit cell to have the x-axis parallel with a surface vector
90 # and z perpendicular to the surface:
91 a1, a2, a3 = surf.cell
92 surf.set_cell([(norm(a1), 0, 0),
93 (np.dot(a1, a2) / norm(a1),
94 np.sqrt(norm(a2)**2 - (np.dot(a1, a2) / norm(a1))**2), 0),
95 (0, 0, norm(a3))],
96 scale_atoms=True)
98 surf.pbc = (True, True, periodic)
100 # Move atoms into the unit cell:
101 scaled = surf.get_scaled_positions()
102 scaled[:, :2] %= 1
103 surf.set_scaled_positions(scaled)
105 if not periodic:
106 surf.cell[2] = 0.0
108 return surf
111def ext_gcd(a, b):
112 if b == 0:
113 return 1, 0
114 elif a % b == 0:
115 return 0, 1
116 else:
117 x, y = ext_gcd(b, a % b)
118 return y, x - y * (a // b)