'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''                                                 ''''''''''
''''''''''                  TubeVBS                        ''''''''''
''''''''''                                                 '''''''''' 
'''''''''' Coded by Roberto Veiga (raveiga@yahoo.com)      ''''''''''
'''''''''' Corrected by David Tomanek (tomanek@pa.msu.edu) ''''''''''
''''''''''                                                 ''''''''''
''''''''''           Version: 20 April 2007                ''''''''''          
''''''''''                                                 ''''''''''
'''''''''' Free non-commercial use.                        ''''''''''
'''''''''' No warranties whatsoever.                       ''''''''''
'''''''''' No source modification allowed without the      ''''''''''
'''''''''' consent of the authors.                         ''''''''''
'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'''''''''''''''''''''''''''''''''''''''''''''''''''''''['''''''''''''

option explicit

'Declare a constant
const pi=3.141592654

'Declare variables
dim n,m,b,ncells
dim level
dim hcd
dim radius
dim a1,a2
dim p1,p2
dim Rm
dim phi,t
dim alpha,h
dim L
dim Nc,dR
dim i
dim file
dim xyz
dim fso,f

'Declare an array of coordinates
dim coord

'Info message
MsgBox "TubeVBS is a Windows Script Host (WSH) script that generates the atomic coordinates in an (n,m) single-walled carbon nanotube (SWCNT). It is based on the work of C. T. White et al., Phys. Rev. B 47, 5485 (1993). Written by Roberto Veiga (raveiga@yahoo.com at Instituto de Física - Universidade Federal de Uberlândia -- http://www.ufu.br). Corrected by David Tomanek at MSU. NO WARRANTY IS PROVIDED FOR THE CORRECTNESS OF THE RESULTS.",0,"TubeVBS"

'Input
n=inputbox("First chiral index 'n' (default=10): ")
m=inputbox("Second chiral index 'm' (default=10): ")
b=inputbox("C-C bond length in graphene sheet (default=1.42 Angstroms): ")
level=inputbox("Level of radial deformation (default: no deformation = 0): ")
ncells=inputbox("Number of translational unit cells (default=1): ")
file=inputbox("File name to store the coordinates in XYZ format (default=nanotube.xyz): ")

if n="" then
    n=10
else
    n=cint(n)
end if

if m="" then
    m=10
else
    m=cint(m)
end if

if b="" then
    b=1.42
else
    b=cdbl(b)
end if

if level="" then
    level=0
else
    level=cint(level)
end if

if ncells="" then
    ncells=1
else
    ncells=cint(ncells)
end if

if file="" then
    file="nanotube.xyz"
end if

'Get the highest common divisor of both n and m
if m=0 then
    hcd=n
else
    for hcd=m to 1 step -1
        if (n mod hcd)=0 then
            if (m mod hcd)=0 then
                exit for
            end if
        end if
    next
end if

'Get L
L=cdbl(m^2+n*m+n^2)

'Get |R|
Rm=b*sqr(3.0*L)

'Get the tubule radius
radius=Rm/(2.0*pi)

'Get a1 and a2
a1=radius+((radius*level)/100.0)
a2=radius-((radius*level)/100.0)

'Get dR
if ((n-m) mod 3*hcd)=0 then
    dR=3.0*hcd
else
    dR=cdbl(hcd)
end if

'Get the number of two-atom unit cells in the translational cell
Nc=2*L/dR

'Redimension the array of coordinates
redim coord(4,2*Nc*ncells)

'Get phi
phi=(pi*(m+n))/L

'Get t
t=b*(n-m)/(2.0*sqr(cdbl(L)))

'Get p1 and p2
for p1=0 to n
    p2=(hcd+p1*m)/n

    if p2=int(p2) then
        exit for
    end if
next

'Get alpha
alpha=pi*(m*(2*p2+p1)+n*(2*p1+p2))/L

'Get h
h=(3.0*hcd*b)/(2.0*sqr(L))

'Set the first atom position in the two-atom unit cell
coord(1,1)=a1
coord(2,1)=0.000000
coord(3,1)=0.000000
coord(4,1)=0.000000

'Set the second atom position in the two-atom unit cell
coord(1,2)=a1*cos(phi)
coord(2,2)=a2*sin(phi)
coord(3,2)=t
coord(4,2)=coord(4,1)+phi

'Complete the tubule helical motif
for i=3 to 2*hcd
    coord(1,i)=a1*cos(coord(4,i-2)+(2.0*pi)/hcd)
    coord(2,i)=a2*sin(coord(4,i-2)+(2.0*pi)/hcd)
    coord(3,i)=coord(3,i-2)
    coord(4,i)=coord(4,i-2)+(2.0*pi)/hcd
next

'Complete the translational unit cell
for i=2*hcd+1 to 2*Nc*ncells
    coord(1,i)=a1*cos((coord(4,i-(2*hcd)))+alpha)
    coord(2,i)=a2*sin((coord(4,i-(2*hcd)))+alpha)
    coord(3,i)=coord(3,i-(2*hcd))+h
    coord(4,i)=coord(4,i-(2*hcd))+alpha
next

'Set XYZ coordinates
xyz=cstr(2*Nc*ncells) & vbCRLF
xyz = xyz & "(" & n & "," & m & ") carbon nanotube structure generated by TubeVBS" & vbCRLF

for i=1 to 2*Nc*ncells
    xyz=xyz & "C          " & cstr(formatnumber(coord(1,i),7,-1)) & "          " & cstr(formatnumber(coord(2,i),7,-1)) & "          " & cstr(formatnumber(coord(3,i),7,-1)) & vbCRLF
next

'Save the coordinates
set fso=CreateObject("Scripting.FileSystemObject")
set f=fso.OpenTextFile(file, 2, true)
f.Write xyz
f.Close

MsgBox "Coordinates were succesfully generated and saved to the specified file!",64,"TubeVBS"
MsgBox "Please acknowledge TubeVBS when using structures generated by this software in your work:" & vbCRLF & vbCRLF & "TubeVBS, R. G. A. Veiga, Universidade Federal de Uberlândia, Brazil, and David Tomanek, MSU, 2007.",64,"TubeVBS"