index | modules | gitlab

# Source code for ase.build.tube

    from math import sqrt

    import numpy as np

    from ase.atoms import Atoms
    from ase.utils import gcd

    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:
            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 * l
        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

Documentation v: 3.19.0

Versions
  3.20.0b1 (development)
  3.19.0 (latest stable)