Coverage for /builds/kinetik161/ase/ase/build/tube.py: 75.63%

1from math import sqrt, gcd 

2 

3import numpy as np 

4 

5from ase.atoms import Atoms 

6 

7 

8def nanotube(n, m, length=1, bond=1.42, symbol='C', verbose=False, 

9 vacuum=None): 

10 """Create an atomic structure. 

11 

12 Creates a single-walled nanotube whose structure is specified using the 

13 standardized (n, m) notation. 

14 

15 Parameters 

16 ---------- 

17 n : int 

18 n in the (n, m) notation. 

19 m : int 

20 m in the (n, m) notation. 

21 length : int, optional 

22 Length (axial repetitions) of the nanotube. 

23 bond : float, optional 

24 Bond length between neighboring atoms. 

25 symbol : str, optional 

26 Chemical element to construct the nanotube from. 

27 verbose : bool, optional 

28 If True, will display key geometric parameters. 

29 

30 Returns 

31 ------- 

32 ase.atoms.Atoms 

33 An ASE Atoms object corresponding to the specified molecule. 

34 

35 Examples 

36 -------- 

37 >>> from ase.build import nanotube 

38 >>> atoms1 = nanotube(6, 0, length=4) 

39 >>> atoms2 = nanotube(3, 3, length=6, bond=1.4, symbol='Si') 

40 """ 

41 if n < m: 

42 m, n = n, m 

43 sign = -1 

44 else: 

45 sign = 1 

46 

47 nk = 6000 

48 sq3 = sqrt(3.0) 

49 a = sq3 * bond 

50 l2 = n * n + m * m + n * m 

51 l1 = sqrt(l2) 

52 

53 nd = gcd(n, m) 

54 if (n - m) % (3 * nd) == 0: 

55 ndr = 3 * nd 

56 else: 

57 ndr = nd 

58 

59 nr = (2 * m + n) // ndr 

60 ns = -(2 * n + m) // ndr 

61 nn = 2 * l2 // ndr 

62 

63 ichk = 0 

64 if nr == 0: 

65 n60 = 1 

66 else: 

67 n60 = nr * 4 

68 

69 absn = abs(n60) 

70 nnp = [] 

71 nnq = [] 

72 for i in range(-absn, absn + 1): 

73 for j in range(-absn, absn + 1): 

74 j2 = nr * j - ns * i 

75 if j2 == 1: 

76 j1 = m * i - n * j 

77 if j1 > 0 and j1 < nn: 

78 ichk += 1 

79 nnp.append(i) 

80 nnq.append(j) 

81 

82 if ichk == 0: 

83 raise RuntimeError('not found p, q strange!!') 

84 if ichk >= 2: 

85 raise RuntimeError('more than 1 pair p, q strange!!') 

86 

87 nnnp = nnp[0] 

88 nnnq = nnq[0] 

89 

90 if verbose: 

91 print('the symmetry vector is', nnnp, nnnq) 

92 

93 lp = nnnp * nnnp + nnnq * nnnq + nnnp * nnnq 

94 r = a * sqrt(lp) 

95 c = a * l1 

96 t = sq3 * c / ndr 

97 

98 if 2 * nn > nk: 

99 raise RuntimeError('parameter nk is too small!') 

100 

101 rs = c / (2.0 * np.pi) 

102 

103 if verbose: 

104 print('radius=', rs, t) 

105 

106 q1 = np.arctan((sq3 * m) / (2 * n + m)) 

107 q2 = np.arctan((sq3 * nnnq) / (2 * nnnp + nnnq)) 

108 q3 = q1 - q2 

109 

110 q4 = 2.0 * np.pi / nn 

111 q5 = bond * np.cos((np.pi / 6.0) - q1) / c * 2.0 * np.pi 

112 

113 h1 = abs(t) / abs(np.sin(q3)) 

114 h2 = bond * np.sin((np.pi / 6.0) - q1) 

115 

116 ii = 0 

117 x, y, z = [], [], [] 

118 for i in range(nn): 

119 x1, y1, z1 = 0, 0, 0 

120 

121 k = np.floor(i * abs(r) / h1) 

122 x1 = rs * np.cos(i * q4) 

123 y1 = rs * np.sin(i * q4) 

124 z1 = (i * abs(r) - k * h1) * np.sin(q3) 

125 kk2 = abs(np.floor((z1 + 0.0001) / t)) 

126 if z1 >= t - 0.0001: 

127 z1 -= t * kk2 

128 elif z1 < 0: 

129 z1 += t * kk2 

130 ii += 1 

131 

132 x.append(x1) 

133 y.append(y1) 

134 z.append(z1) 

135 z3 = (i * abs(r) - k * h1) * np.sin(q3) - h2 

136 ii += 1 

137 

138 if z3 >= 0 and z3 < t: 

139 x2 = rs * np.cos(i * q4 + q5) 

140 y2 = rs * np.sin(i * q4 + q5) 

141 z2 = (i * abs(r) - k * h1) * np.sin(q3) - h2 

142 x.append(x2) 

143 y.append(y2) 

144 z.append(z2) 

145 else: 

146 x2 = rs * np.cos(i * q4 + q5) 

147 y2 = rs * np.sin(i * q4 + q5) 

148 z2 = (i * abs(r) - (k + 1) * h1) * np.sin(q3) - h2 

149 kk = abs(np.floor(z2 / t)) 

150 if z2 >= t - 0.0001: 

151 z2 -= t * kk 

152 elif z2 < 0: 

153 z2 += t * kk 

154 x.append(x2) 

155 y.append(y2) 

156 z.append(z2) 

157 

158 ntotal = 2 * nn 

159 X = [] 

160 for i in range(ntotal): 

161 X.append([x[i], y[i], sign * z[i]]) 

162 

163 if length > 1: 

164 xx = X[:] 

165 for mnp in range(2, length + 1): 

166 for i in range(len(xx)): 

167 X.append(xx[i][:2] + [xx[i][2] + (mnp - 1) * t]) 

168 

169 transvec = t 

170 numatom = ntotal * length 

171 diameter = rs * 2 

172 chiralangle = np.arctan((sq3 * n) / (2 * m + n)) / np.pi * 180 

173 

174 cell = [[0, 0, 0], [0, 0, 0], [0, 0, length * t]] 

175 atoms = Atoms(symbol + str(numatom), 

176 positions=X, 

177 cell=cell, 

178 pbc=[False, False, True]) 

179 if vacuum: 

180 atoms.center(vacuum, axis=(0, 1)) 

181 if verbose: 

182 print('translation vector =', transvec) 

183 print('diameter = ', diameter) 

184 print('chiral angle = ', chiralangle) 

185 return atoms