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