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