Coverage for /builds/kinetik161/ase/ase/build/surfaces_with_termination.py: 98.72%

1import numpy as np 

2from ase.build.general_surface import surface 

3from ase.geometry import get_layers 

4from ase.symbols import string2symbols 

5 

6 

7def surfaces_with_termination(lattice, indices, layers, vacuum=None, tol=1e-10, 

8 termination=None, return_all=False, 

9 verbose=False): 

10 """Create surface from a given lattice and Miller indices with a given 

11 termination 

12 

13 Parameters 

14 ========== 

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. (not the same as the layers 

24 you choose from for terminations) 

25 vacuum: float 

26 Amount of vacuum added on both sides of the slab. 

27 termination: str 

28 the atoms you wish to be in the top layer. There may be many such 

29 terminations, this function returns all terminations with the same 

30 atomic composition. 

31 e.g. 'O' will return oxygen terminated surfaces. 

32 e.g.'TiO' returns surfaces terminated with layers containing both 

33 O and Ti 

34 Returns: 

35 return_surfs: List 

36 a list of surfaces that match the specifications given 

37 

38 """ 

39 lats = translate_lattice(lattice, indices) 

40 return_surfs = [] 

41 check = [] 

42 check2 = [] 

43 for item in lats: 

44 too_similar = False 

45 surf = surface(item, indices, layers, vacuum=vacuum, tol=tol) 

46 surf.wrap(pbc=[True] * 3) # standardize slabs 

47 

48 positions = surf.get_scaled_positions().flatten() 

49 for i, value in enumerate(positions): 

50 if value >= 1 - tol: # move things closer to zero within tol 

51 positions[i] -= 1 

52 surf.set_scaled_positions(np.reshape(positions, (len(surf), 3))) 

53 # rep = find_z_layers(surf) 

54 z_layers, hs = get_layers(surf, (0, 0, 1)) # just z layers matter 

55 # get the indicies of the atoms in the highest layer 

56 top_layer = [ 

57 i for i, val in enumerate( 

58 z_layers == max(z_layers)) if val] 

59 

60 if termination is not None: 

61 comp = [surf.get_chemical_symbols()[a] for a in top_layer] 

62 term = string2symbols(termination) 

63 # list atoms in top layer and not in requested termination 

64 check = [a for a in comp if a not in term] 

65 # list of atoms in requested termination and not in top layer 

66 check2 = [a for a in term if a not in comp] 

67 if len(return_surfs) > 0: 

68 pos_diff = [a.get_positions() - surf.get_positions() 

69 for a in return_surfs] 

70 for i, su in enumerate(pos_diff): 

71 similarity_test = su.flatten() < tol * 1000 

72 if similarity_test.all(): 

73 # checks if surface is too similar to another surface 

74 too_similar = True 

75 if too_similar: 

76 continue 

77 if return_all is True: 

78 pass 

79 elif check != [] or check2 != []: 

80 continue 

81 return_surfs.append(surf) 

82 return return_surfs 

83 

84 

85def translate_lattice(lattice, indices, tol=10**-3): 

86 """translates a bulk unit cell along a normal vector given by the a set of 

87 miller indices to the next symetric position. This is used to control the 

88 termination of the surface in the smart_surface command 

89 Parameters: 

90 ========== 

91 lattice: Atoms object 

92 atoms object of the bulk unit cell 

93 indices: 1x3 list,tuple, or numpy array 

94 the miller indices you wish to cut along. 

95 returns: 

96 lattice_list: list of Atoms objects 

97 a list of all the different translations of the unit cell that will 

98 yield different terminations of a surface cut along the miller 

99 indices provided. 

100 """ 

101 lattice_list = [] 

102 cell = lattice.get_cell() 

103 pt = [0, 0, 0] 

104 h, k, l = indices # noqa (E741 ambiguous name 'l') 

105 millers = list(indices) 

106 for index, item in enumerate(millers): 

107 if item == 0: 

108 millers[index] = 10**9 # make zeros large numbers 

109 elif pt == [0, 0, 0]: # for numerical stability 

110 pt = list(cell[index] / float(item) / np.linalg.norm(cell[index])) 

111 h1, k1, l1 = millers 

112 N = np.array(cell[0] / h1 + cell[1] / k1 + cell[2] / l1) 

113 n = N / np.linalg.norm(N) # making a unit vector normal to cut plane 

114 # finding distance from cut plan vector 

115 d = [np.round(np.dot(n, (a - pt)) * n, 5) for 

116 a in lattice.get_scaled_positions()] 

117 duplicates = [] 

118 for i, item in enumerate(d): 

119 g = [True for a in d[i + 1:] if np.linalg.norm(a - item) < tol] 

120 if g != []: 

121 duplicates.append(i) 

122 duplicates.reverse() 

123 for i in duplicates: 

124 del d[i] 

125 # put distance to the plane at the end of the array 

126 for i, item in enumerate(d): 

127 d[i] = np.append(item, 

128 np.dot(n, (lattice.get_scaled_positions()[i] - pt))) 

129 d = np.array(d) 

130 d = d[d[:, 3].argsort()] # sort by distance to the plane 

131 d = [a[:3] for a in d] # remove distance 

132 d = list(d) # make it a list again 

133 for i in d: 

134 """ 

135 The above method gives you the boundries of between terminations that 

136 will allow you to build a complete set of terminations. However, it 

137 does not return all the boundries. Thus you must check both above and 

138 below the boundary, and not stray too far from the boundary. If you move 

139 too far away, you risk hitting another boundary you did not find. 

140 """ 

141 lattice1 = lattice.copy() 

142 displacement = (h * cell[0] + k * cell[1] + l * cell[2]) \ 

143 * (i + 10 ** -8) 

144 lattice1.positions -= displacement 

145 lattice_list.append(lattice1) 

146 lattice1 = lattice.copy() 

147 displacement = (h * cell[0] + k * cell[1] + l * cell[2]) \ 

148 * (i - 10 ** -8) 

149 lattice1.positions -= displacement 

150 lattice_list.append(lattice1) 

151 return lattice_list