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):
"""Create an atomic structure.
Creates a single-walled nanotube whose structure is specified using the
standardized (n, m) notation.
Parameters
----------
n : int
n in the (n, m) notation.
m : int
m in the (n, m) notation.
length : int, optional
Length (axial repetitions) of the nanotube.
bond : float, optional
Bond length between neighboring atoms.
symbol : str, optional
Chemical element to construct the nanotube from.
verbose : bool, optional
If True, will display key geometric parameters.
Returns
-------
ase.atoms.Atoms
An ASE Atoms object corresponding to the specified molecule.
Examples
--------
>>> from ase.build import nanotube
>>> atoms1 = nanotube(6, 0, length=4)
>>> atoms2 = nanotube(3, 3, length=6, bond=1.4, symbol='Si')
"""
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
l1 = 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:
raise RuntimeError('not found p, q strange!!')
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 * l1
t = sq3 * c / ndr
if 2 * nn > nk:
raise RuntimeError('parameter nk is too small!')
rs = c / (2.0 * np.pi)
if verbose:
print('radius=', rs, t)
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