Coverage for /builds/ase/ase/ase/build/root.py: 100.00%
71 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3from math import atan2, cos, log10, sin
5import numpy as np
8def hcp0001_root(symbol, root, size, a=None, c=None,
9 vacuum=None, orthogonal=False):
10 """HCP(0001) surface maniupulated to have a x unit side length
11 of *root* before repeating. This also results in *root* number
12 of repetitions of the cell.
15 The first 20 valid roots for nonorthogonal are...
16 1, 3, 4, 7, 9, 12, 13, 16, 19, 21, 25,
17 27, 28, 31, 36, 37, 39, 43, 48, 49"""
18 from ase.build import hcp0001
19 atoms = hcp0001(symbol=symbol, size=(1, 1, size[2]),
20 a=a, c=c, vacuum=vacuum, orthogonal=orthogonal)
21 atoms = root_surface(atoms, root)
22 atoms *= (size[0], size[1], 1)
23 return atoms
26def fcc111_root(symbol, root, size, a=None,
27 vacuum=None, orthogonal=False):
28 """FCC(111) surface maniupulated to have a x unit side length
29 of *root* before repeating. This also results in *root* number
30 of repetitions of the cell.
32 The first 20 valid roots for nonorthogonal are...
33 1, 3, 4, 7, 9, 12, 13, 16, 19, 21, 25, 27,
34 28, 31, 36, 37, 39, 43, 48, 49"""
35 from ase.build import fcc111
36 atoms = fcc111(symbol=symbol, size=(1, 1, size[2]),
37 a=a, vacuum=vacuum, orthogonal=orthogonal)
38 atoms = root_surface(atoms, root)
39 atoms *= (size[0], size[1], 1)
40 return atoms
43def bcc111_root(symbol, root, size, a=None,
44 vacuum=None, orthogonal=False):
45 """BCC(111) surface maniupulated to have a x unit side length
46 of *root* before repeating. This also results in *root* number
47 of repetitions of the cell.
50 The first 20 valid roots for nonorthogonal are...
51 1, 3, 4, 7, 9, 12, 13, 16, 19, 21, 25,
52 27, 28, 31, 36, 37, 39, 43, 48, 49"""
53 from ase.build import bcc111
54 atoms = bcc111(symbol=symbol, size=(1, 1, size[2]),
55 a=a, vacuum=vacuum, orthogonal=orthogonal)
56 atoms = root_surface(atoms, root)
57 atoms *= (size[0], size[1], 1)
58 return atoms
61def point_in_cell_2d(point, cell, eps=1e-8):
62 """This function takes a 2D slice of the cell in the XY plane and calculates
63 if a point should lie in it. This is used as a more accurate method of
64 ensuring we find all of the correct cell repetitions in the root surface
65 code. The Z axis is totally ignored but for most uses this should be fine.
66 """
67 # Define area of a triangle
68 def tri_area(t1, t2, t3):
69 t1x, t1y = t1[0:2]
70 t2x, t2y = t2[0:2]
71 t3x, t3y = t3[0:2]
72 return abs(t1x * (t2y - t3y) + t2x *
73 (t3y - t1y) + t3x * (t1y - t2y)) / 2
75 # c0, c1, c2, c3 define a parallelogram
76 c0 = (0, 0)
77 c1 = cell[0, 0:2]
78 c2 = cell[1, 0:2]
79 c3 = c1 + c2
81 # Get area of parallelogram
82 cA = tri_area(c0, c1, c2) + tri_area(c1, c2, c3)
84 # Get area of triangles formed from adjacent vertices of parallelogram and
85 # point in question.
86 pA = tri_area(point, c0, c1) + tri_area(point, c1, c2) + \
87 tri_area(point, c2, c3) + tri_area(point, c3, c0)
89 # If combined area of triangles from point is larger than area of
90 # parallelogram, point is not inside parallelogram.
91 return pA <= cA + eps
94def _root_cell_normalization(primitive_slab):
95 """Returns the scaling factor for x axis and cell normalized by that
96 factor"""
98 xscale = np.linalg.norm(primitive_slab.cell[0, 0:2])
99 cell_vectors = primitive_slab.cell[0:2, 0:2] / xscale
100 return xscale, cell_vectors
103def _root_surface_analysis(primitive_slab, root, eps=1e-8):
104 """A tool to analyze a slab and look for valid roots that exist, up to
105 the given root. This is useful for generating all possible cells
106 without prior knowledge.
108 *primitive slab* is the primitive cell to analyze.
110 *root* is the desired root to find, and all below.
112 This is the internal function which gives extra data to root_surface.
113 """
115 # Setup parameters for cell searching
116 logeps = int(-log10(eps))
117 _xscale, cell_vectors = _root_cell_normalization(primitive_slab)
119 # Allocate grid for cell search search
120 points = np.indices((root + 1, root + 1)).T.reshape(-1, 2)
122 # Find points corresponding to full cells
123 cell_points = [cell_vectors[0] * x + cell_vectors[1] * y for x, y in points]
125 # Find point close to the desired cell (floating point error possible)
126 roots = np.around(np.linalg.norm(cell_points, axis=1)**2, logeps)
128 valid_roots = np.nonzero(roots == root)[0]
129 if len(valid_roots) == 0:
130 raise ValueError(
131 "Invalid root {} for cell {}".format(
132 root, cell_vectors))
133 int_roots = np.array([int(this_root) for this_root in roots
134 if this_root.is_integer() and this_root <= root])
135 return cell_points, cell_points[np.nonzero(
136 roots == root)[0][0]], set(int_roots[1:])
139def root_surface_analysis(primitive_slab, root, eps=1e-8):
140 """A tool to analyze a slab and look for valid roots that exist, up to
141 the given root. This is useful for generating all possible cells
142 without prior knowledge.
144 *primitive slab* is the primitive cell to analyze.
146 *root* is the desired root to find, and all below."""
147 return _root_surface_analysis(
148 primitive_slab=primitive_slab, root=root, eps=eps)[2]
151def root_surface(primitive_slab, root, eps=1e-8):
152 """Creates a cell from a primitive cell that repeats along the x and y
153 axis in a way consisent with the primitive cell, that has been cut
154 to have a side length of *root*.
156 *primitive cell* should be a primitive 2d cell of your slab, repeated
157 as needed in the z direction.
159 *root* should be determined using an analysis tool such as the
160 root_surface_analysis function, or prior knowledge. It should always
161 be a whole number as it represents the number of repetitions."""
163 atoms = primitive_slab.copy()
165 xscale, cell_vectors = _root_cell_normalization(primitive_slab)
167 # Do root surface analysis
168 cell_points, root_point, _roots = _root_surface_analysis(
169 primitive_slab, root, eps=eps)
171 # Find new cell
172 root_angle = -atan2(root_point[1], root_point[0])
173 root_rotation = [[cos(root_angle), -sin(root_angle)],
174 [sin(root_angle), cos(root_angle)]]
175 root_scale = np.linalg.norm(root_point)
177 cell = np.array([np.dot(x, root_rotation) *
178 root_scale for x in cell_vectors])
180 # Find all cell centers within the cell
181 shift = cell_vectors.sum(axis=0) / 2
182 cell_points = [
183 point for point in cell_points if point_in_cell_2d(
184 point + shift, cell, eps=eps)]
186 # Setup new cell
187 atoms.rotate(root_angle, v="z")
188 atoms *= (root, root, 1)
189 atoms.cell[0:2, 0:2] = cell * xscale
190 atoms.center()
192 # Remove all extra atoms
193 del atoms[[atom.index for atom in atoms if not point_in_cell_2d(
194 atom.position, atoms.cell, eps=eps)]]
196 # Rotate cell back to original orientation
197 standard_rotation = [[cos(-root_angle), -sin(-root_angle), 0],
198 [sin(-root_angle), cos(-root_angle), 0],
199 [0, 0, 1]]
201 new_cell = np.array([np.dot(x, standard_rotation) for x in atoms.cell])
202 new_positions = np.array([np.dot(x, standard_rotation)
203 for x in atoms.positions])
205 atoms.cell = new_cell
206 atoms.positions = new_positions
207 return atoms