# Source code for ase.build.tube

from math import sqrt, gcd

import numpy as np

from ase.atoms import Atoms

[docs]def nanotube(n, m, length=1, bond=1.42, symbol='C', verbose=False,
vacuum=None):
if n < m:
m, n = n, m
sign = -1
else:
sign = 1

nk = 6000
sq3 = sqrt(3.0)
a = sq3 * bond
l2 = n * n + m * m + n * m
l = sqrt(l2)

nd = gcd(n, m)
if (n - m) % (3 * nd) == 0:
ndr = 3 * nd
else:
ndr = nd

nr = (2 * m + n) // ndr
ns = -(2 * n + m) // ndr
nn = 2 * l2 // ndr

ichk = 0
if nr == 0:
n60 = 1
else:
n60 = nr * 4

absn = abs(n60)
nnp = []
nnq = []
for i in range(-absn, absn + 1):
for j in range(-absn, absn + 1):
j2 = nr * j - ns * i
if j2 == 1:
j1 = m * i - n * j
if j1 > 0 and j1 < nn:
ichk += 1
nnp.append(i)
nnq.append(j)

if ichk == 0:
if ichk >= 2:
raise RuntimeError('more than 1 pair p, q strange!!')

nnnp = nnp[0]
nnnq = nnq[0]

if verbose:
print('the symmetry vector is', nnnp, nnnq)

lp = nnnp * nnnp + nnnq * nnnq + nnnp * nnnq
r = a * sqrt(lp)
c = a * l
t = sq3 * c / ndr

if 2 * nn > nk:
raise RuntimeError('parameter nk is too small!')

rs = c / (2.0 * np.pi)

if verbose:

q1 = np.arctan((sq3 * m) / (2 * n + m))
q2 = np.arctan((sq3 * nnnq) / (2 * nnnp + nnnq))
q3 = q1 - q2

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

h1 = abs(t) / abs(np.sin(q3))
h2 = bond * np.sin((np.pi / 6.0) - q1)

ii = 0
x, y, z = [], [], []
for i in range(nn):
x1, y1, z1 = 0, 0, 0

k = np.floor(i * abs(r) / h1)
x1 = rs * np.cos(i * q4)
y1 = rs * np.sin(i * q4)
z1 = (i * abs(r) - k * h1) * np.sin(q3)
kk2 = abs(np.floor((z1 + 0.0001) / t))
if z1 >= t - 0.0001:
z1 -= t * kk2
elif z1 < 0:
z1 += t * kk2
ii += 1

x.append(x1)
y.append(y1)
z.append(z1)
z3 = (i * abs(r) - k * h1) * np.sin(q3) - h2
ii += 1

if z3 >= 0 and z3 < t:
x2 = rs * np.cos(i * q4 + q5)
y2 = rs * np.sin(i * q4 + q5)
z2 = (i * abs(r) - k * h1) * np.sin(q3) - h2
x.append(x2)
y.append(y2)
z.append(z2)
else:
x2 = rs * np.cos(i * q4 + q5)
y2 = rs * np.sin(i * q4 + q5)
z2 = (i * abs(r) - (k + 1) * h1) * np.sin(q3) - h2
kk = abs(np.floor(z2 / t))
if z2 >= t - 0.0001:
z2 -= t * kk
elif z2 < 0:
z2 += t * kk
x.append(x2)
y.append(y2)
z.append(z2)

ntotal = 2 * nn
X = []
for i in range(ntotal):
X.append([x[i], y[i], sign * z[i]])

if length > 1:
xx = X[:]
for mnp in range(2, length + 1):
for i in range(len(xx)):
X.append(xx[i][:2] + [xx[i][2] + (mnp - 1) * t])

transvec = t
numatom = ntotal * length
diameter = rs * 2
chiralangle = np.arctan((sq3 * n) / (2 * m + n)) / np.pi * 180

cell = [[0, 0, 0], [0, 0, 0], [0, 0, length * t]]
atoms = Atoms(symbol + str(numatom),
positions=X,
cell=cell,
pbc=[False, False, True])
if vacuum:
atoms.center(vacuum, axis=(0, 1))
if verbose:
print('translation vector =', transvec)
print('diameter = ', diameter)
print('chiral angle = ', chiralangle)
return atoms