 r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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),

17 a=a, c=c, vacuum=vacuum, orthogonal=orthogonal)

18 atoms = root_surface(atoms, root)

19 atoms *= (size, size, 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),

33 a=a, vacuum=vacuum, orthogonal=orthogonal)

34 atoms = root_surface(atoms, root)

35 atoms *= (size, size, 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),

50 a=a, vacuum=vacuum, orthogonal=orthogonal)

51 atoms = root_surface(atoms, root)

52 atoms *= (size, size, 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 * x + cell_vectors * 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)

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)], 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)

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, root_point)

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