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

1# fmt: off 

2 

3from math import gcd, sqrt 

4 

5import numpy as np 

6 

7from ase.atoms import Atoms 

8 

9 

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

11 vacuum=None): 

12 """Create an atomic structure. 

13 

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

15 standardized (n, m) notation. 

16 

17 Parameters 

18 ---------- 

19 n : int 

20 n in the (n, m) notation. 

21 m : int 

22 m in the (n, m) notation. 

23 length : int, optional 

24 Length (axial repetitions) of the nanotube. 

25 bond : float, optional 

26 Bond length between neighboring atoms. 

27 symbol : str, optional 

28 Chemical element to construct the nanotube from. 

29 verbose : bool, optional 

30 If True, will display key geometric parameters. 

31 

32 Returns 

33 ------- 

34 ase.atoms.Atoms 

35 An ASE Atoms object corresponding to the specified molecule. 

36 

37 Examples 

38 -------- 

39 >>> from ase.build import nanotube 

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

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

42 """ 

43 if n < m: 

44 m, n = n, m 

45 sign = -1 

46 else: 

47 sign = 1 

48 

49 nk = 6000 

50 sq3 = sqrt(3.0) 

51 a = sq3 * bond 

52 l2 = n * n + m * m + n * m 

53 l1 = sqrt(l2) 

54 

55 nd = gcd(n, m) 

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

57 ndr = 3 * nd 

58 else: 

59 ndr = nd 

60 

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

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

63 nn = 2 * l2 // ndr 

64 

65 ichk = 0 

66 if nr == 0: 

67 n60 = 1 

68 else: 

69 n60 = nr * 4 

70 

71 absn = abs(n60) 

72 nnp = [] 

73 nnq = [] 

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

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

76 j2 = nr * j - ns * i 

77 if j2 == 1: 

78 j1 = m * i - n * j 

79 if j1 > 0 and j1 < nn: 

80 ichk += 1 

81 nnp.append(i) 

82 nnq.append(j) 

83 

84 if ichk == 0: 

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

86 if ichk >= 2: 

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

88 

89 nnnp = nnp[0] 

90 nnnq = nnq[0] 

91 

92 if verbose: 

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

94 

95 lp = nnnp * nnnp + nnnq * nnnq + nnnp * nnnq 

96 r = a * sqrt(lp) 

97 c = a * l1 

98 t = sq3 * c / ndr 

99 

100 if 2 * nn > nk: 

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

102 

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

104 

105 if verbose: 

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

107 

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

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

110 q3 = q1 - q2 

111 

112 q4 = 2.0 * np.pi / nn 

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

114 

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

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

117 

118 ii = 0 

119 x, y, z = [], [], [] 

120 for i in range(nn): 

121 x1, y1, z1 = 0, 0, 0 

122 

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

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

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

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

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

128 if z1 >= t - 0.0001: 

129 z1 -= t * kk2 

130 elif z1 < 0: 

131 z1 += t * kk2 

132 ii += 1 

133 

134 x.append(x1) 

135 y.append(y1) 

136 z.append(z1) 

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

138 ii += 1 

139 

140 if z3 >= 0 and z3 < t: 

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

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

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

144 x.append(x2) 

145 y.append(y2) 

146 z.append(z2) 

147 else: 

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

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

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

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

152 if z2 >= t - 0.0001: 

153 z2 -= t * kk 

154 elif z2 < 0: 

155 z2 += t * kk 

156 x.append(x2) 

157 y.append(y2) 

158 z.append(z2) 

159 

160 ntotal = 2 * nn 

161 X = [] 

162 for i in range(ntotal): 

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

164 

165 if length > 1: 

166 xx = X[:] 

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

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

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

170 

171 transvec = t 

172 numatom = ntotal * length 

173 diameter = rs * 2 

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

175 

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

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

178 positions=X, 

179 cell=cell, 

180 pbc=[False, False, True]) 

181 if vacuum: 

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

183 if verbose: 

184 print('translation vector =', transvec) 

185 print('diameter = ', diameter) 

186 print('chiral angle = ', chiralangle) 

187 return atoms