• Overview and goals
• Building a nanotube
• The cylinder and its flat equivalent
• The helical symmetry
• Input and output
• Carbon nanotubes: a simple case
• Stability and relaxation energy
• Electronic structure
• Chrysotile: a more complex case
• Structure
• Input and output
• Small collection of CRYSTAL inputs and outputs
• References

For a correct visualisation of figures and animations present in this tutorial we advise the use of Mozilla Firefox or Safari browsers. Javascript (Java plug-in with Java Runtime Environment) must be enabled on your machine. The following box informs you if Javascript is enabled.

Java is NOT enabled on this machine!

If the test fails, click here to solve this problem. For further information click here.

Overview and goals

This tutorial describes how nanotube systems (1D periodicity) can be simulated with the CRYSTAL code. First, we focus on building a nanotube by wrapping the corresponding 2D flat structure. In this part we introduce the main features of nanotubes for the simple case of graphite → graphene → SWCNT (Single Walled Carbon Nanotube). Then, a few more general examples will be given.
Nanotubes can contain hundreds or thousands of atoms in their unit cell and are characterised by helical symmetry (cyclic groups, only rotation and roto-translation operators, high number of symmetry operators). In CRYSTAL, the full exploitation of point symmetry in direct and reciprocal space (integral calculation, Fock matrix diagonalisation) is implemented, so that these calculations are affordable with a relatively small amount of CPU time and resources.

Building a nanotube

The cylinder and its flat equivalent

The vectors and the helical symmetry

Nanotubes are cylindrical structures, which can be periodic along a single direction (single lattice parameter). They can be easily modelled by wrapping the corresponding 2D flat layer (see FIG. 1) according to the rolling vector, R, defined as R= n₁a₁ + n₂a₂ (where a₁ and a₂ are the slab cell vectors). The (n₁,n₂) indices fully define the nanotube, providing diameter (|R| is the circumference) and chirality, and are used in the literature to characterise the systems (for example the (4,3) nanotube is a nanotube with a rolling vector R=4a₁+3a₂).

If n₁ and n₂ have a common divisor, N, as is the case shown in FIGs. 1 and 2, R passes N times through the 2D lattice nodes, so that a N-order rotational axis perpendicular to R exists in the cylinder.

In order to establish a correspondence between the translational symmetry of the flat slab and the translational and roto-translational symmetry of the curved surface, two more vectors are defined on the basis of R: the longitudinal vector (translational symmetry, 1D periodicity) and the helical vector (helical symmetry, from translation in the slab to roto-translation in the cylinder). In FIGs. from 1 to 8, the rolling and longitudinal vectors are indicated in red and blue colour, respectively. In FIGs. 3 and 4 the helical vector is indicated in green.

The longitudinal vector, L = l₁a₁ + l₂a₂, is defined as the shortest lattice vector normal to R. In the nanotube, it becomes the 1D lattice parameter and gives the 1D periodicity along the tube axis (FIG. 2).

FIG. 1 : From the flat nanotube to the cylinder: once the nanotube is wrapped, the rolling vector becomes a circle normal to the cylinder axis, its norm corresponding to the cylinder perimeter. The (6,3) carbon nanotube example. Click here or on the image to see rolling up animation.	FIG. 2 : R and L vectors for the (6,3) carbon nanotube. The green area represents the graphene unit cell, the green vectors its cell parameters. Click here or on the image to see the animation.

The helical vector, H=h₁a₁ + h₂a₂, is a lattice vector defining with R an area which is N times the area of the unit cell of the flat slab (see FIG. 3, remember that N is the order of the rotational axis). It satisfies, then, the following relationship

where S(v_i,v_j) is the surface defined by the v_i and v_j vectors. In particular, in Cartesian coordinates, we can write:

With the H vector we want to establish a correspondence between a translation in the flat slab and a roto-translation in the cylinder. In the following, we are going to show how the definition of the H vector given in the previous equations is obtained, why the choice of H is not unique and which is the convention adopted in the CRYSTAL code.

The definition of the H vector comes out from two main considerations: 1) in order to maintain the 1D periodicity along L, the S(R,L) area must be multiple of the S(R,H) area; 2) as a consequence, since S(R,L) defines a supercell of the unit cell of the flat slab, it must be a multiple of the unit cell surface. If m, n and k are integers, we can write If R passes N times through the lattice nodes, R is a multiple of a translation vector in the flat slab (that generates a N-order rotational axis in the tube), and crosses N lattice nodes. As a consequence, the surface defined by any vector of the lattice (then also H) and R is always larger or equal to N times S(a₁,a₂):

with a being an integer. According to the previous equations, as k and N are fixed for a given nanotube, a and n can assume all the values such that

Let us consider, for example, the (4,2) SWCNT, where N=2 and k=28 (FIG. 3). Four sets of (a,n) values are possible:

• a=1, n=14: maximum symmetry order, minimum S(R,H)
• a=2, n=7
• a=7, n=2
• a=14, n=1: H=L

The choice of H is not unique, because the same area can be obtained by choosing different lattice vectors forming with R different angles. For example, consider the flat layer mapped as a (4,2) SWCNT represented in FIG. 3; two choices of H are shown: H(a)=(1,1) and H(b)=(3,2). In both cases |n₁h₂-n₂h₁|=2. It is convenient to choose H as the shortest vector in the unit cell satisfying S(R,H)=N S(a₁,a₂) (i.e. H(a) in this example).

H can be decomposed into two vectors, one parallel to R and one to L, the two components being H_R and H_L. The translation along H in the flat plane corresponds to a rotation along R and a translation along L in the tube; H defines, then, the roto-translation operator in the tube. The order M (corresponding to a in previous equations) of this roto-translation is obtained by considering how many times the S(R,H) area must be translated along the L vector in order to fill the nanotube unit cell.

In the example of FIG. 3, R=(4,2), H(a)=(1,1), L=(-4,5), so that M=14.

The two sets of operators (N-order rotational axis and M-order roto-translation), lead to the symmetry group of the tube through the direct product of the related groups. The total number of symmetry operators, S is then given by

FIG. 3: Graphene mapped to build a (4,2) SWCNT. Lattice nodes are represented in solid circles. R, L and two different choices of the H vector are reported. The green area represents the irreducible unit cell. As R passes 2 times through the nodes, an axis of second order is present, whereas the translation along the tube axis is 1/14 (order of the roto-translation operator: M=14). The direct product of the rotational times the roto-translational groups gives a helical group with order 28. Note that the unit cell area is 28 times larger than the green area. (a) H=(1,1), the shortest: HR=9/28|R| and HL=1/14|L|. The progressive generation of the atoms by the roto-translation is shown; the direction of the helical vector is shown for the first atoms with green dot-lines. (b) H=(3,2), an alternative choice: HR=23/28|R| and HL=1/14|L|. The rest as in (a).

In summary, given a 2D unit cell with a1 and a2 basis vectors, three vectors define: the nanotube: R= n1a1 + n2a2 (it defines the nanotube radius and chirality, i.e. how the flat slab is wrapped into a cylinder; |R| is the nanotube circumference). The definition of L and H depends on R. the 1D periodicity: L= l1a1 + l2a2 (it is the nanotube lattice parameter, chosen as the shortest lattice vector perpendicular to R) the helical symmetry: H= h1a1 + h2a2 (it defines a correspondence between a translation in the flat slab and a roto-translation in the curved surface; H has a rotational component along the circumference vector, R, and a translational component along the lattice parameter, L) The described vectors and quantities and their relationship with the helical symmetry are represented for the (4,2) SWCNT in FIG. 3 and a 3D interactive animation for the (6,3) SWCNT is presented in FIG. 4.

FIG. 4: The nanotube symmetry. Click here or on the image to load the animation (if you are using Mozilla Firefox on Mac OS X you should have problems to visualise this animation).

Chirality

In case of hexagonal lattices, due to the observed C-C bond motives on the carbon nanotube walls, when n1=n2 the nanotubes belong to the so-called armchair family, whereas when n2=0 or n1=0 they are called zigzag tubes. For all the other combinations of n1 and n1 , they are said chiral. The chiral angle is defined as the angle between R and the a1 unit cell vector of the flat slab, and is indicated with θ (see FIG. 5). This angle ranges from 0° (zigzag) to 30° (armchair) and can be experimentally measured. It depends on n1 and n2 only, and is obtained with the following equation: A few 3D interactive examples are present in FIG. 6.

FIG. 5: The graphene lattice and the rolling directions to generate nanotubes. In this figure, the unit cell vectors are indicated in green, the (n,0) (θ=0°) and (n,n) (θ=30°) rolling directions in red. The shaded area contains all possible symmetry nonequivalent rolling directions.

(6,0) θ=0.0° Zigzag structure	(6,2) θ=13.9°	(6,4) θ=23.4°	(6,6) θ=30.0° Armchair structure
FIG. 6: Chiral angle (θ) for 4 different rolling directions (indicated into parentheses). θ is the angle between the violet line and the R rolling direction (in red). Drag with the mouse to rotate.

The 1D periodicity

The periodicity along the tube axis, i.e. the existence of the longitudinal vector, is not satisfied for all possible 2D lattices. The orthogonality condition between R and L provides the following equation:

from which one obtains:

• Hexagonal lattice: a₁=a₂, cos(γ)=±1/2. The equation becomes
  
  Any roll-up vector is possible.

• Square lattice: a₁=a₂, cos(γ)=0.
  
  Any roll-up vector is possible.

• Rectangular lattice: a₁≠a₂, cos(γ)=0
  
  In this case the right term corresponds to a rational number either if |a₁|=n|a₂|, with n being a rational number, or if n₂=0. More generally, for rectangular lattices, the periodicity along the tube axis is always satisfied for (n,0) and (0,n) rolling vectors.

• Rhombohedral (or centred rectangular) lattice: a₁=a₂, any cos(γ)=0.
  
  The right term provides a rational number only when n₁=n₂ or n₁=-n₂, so that:
  
• Oblique lattice: a₁≠a₂, any cos(γ)=0. The orthogonality equation remains as such, and the right term is, in general, an irrational number.

Input and output

The CRYSTAL keywords associated to nanotubes are listed in this section; they must be inserted in the geometry input block, following the specifications of the starting slab structure.

This keyword builds a nanotube from a slab. The components of the R vector are expressed in the slab basis vectors.
1 - WARNING: in CRYSTAL if the 2D lattice is hexagonal, γ, the angle between a₁ and a₂, is 120°. To use a γ=60° lattice as reference (as done, for example, in the SWCNT literature), use SWCNT (see below).
2 - WARNING: the nanotube features require the condition N1≥N2. In order to explore rolling directions with N1<N2, it is sufficient to exchange the x/a₁ and y/a₂ columns in the slab input, and substitute N1 with N2 and vice versa.
3 - WARNING: at present, a maximum of 48 symmetry operators is possible. In case of larger number of symmetry operators, it is sufficient to enlarge the asymmetric unit in order to fulfill this condition.
4 - WARNING: at present, this strategy allows to build single-walled nanotubes only. Multi-walled nanotubes can be simulated with the CRYSTAL code without the automatic structure generation and the exploitation of the helical symmetry.

This keyword build a nanotube from a hexagonal slab. In the literature, the rolling direction are commonly specified in a &gamma=60° lattice. When the SWCNT keyword is used instead of NANOTUBE (see previous keyword), the R vector components are specified in this basis. CRYSTAL will automatically converts the R direction in the working basis (&gamma=120°). The relationship between the indices defined in the &gamma=120° and the &gamma=60° basis is the following:
N1(120)=N1(60)+N2(60)
N2(120)=N2(60)
and can be obtained by considering the flat slab unit cells as represented in FIG. 7.

FIG. 7: Representation of the &gamma=60° (black) and &gamma=120° (red) unit cells. The equivalence of the two choices is shown in FIG. 8 for the (6,3) SWCNT; examples of possible input structure to generate the (6,3) SWCNT are supplied in the following table.

FIG. 8: The (6,3) carbon nanotube obtained from &gamma=60° and &gamma=120° unit cell flat slabs. Click here or on the image to see the animation.

Four totally equivalent inputs for the (6,3) SWCNT are reported in the following table.

These input blocks can be tested using the TEST keyword at the end of the geometry section. After printing the information relative to graphene (or graphite and graphene, if we choose to start from the 3D lattice), the vectors described in the previous paragraphs are listed, followed by the coordinate of all the generated atoms and operators.

The following table shows the part of the output listing the nanotube information.

 *******************************************************************************
 *                 CONSTRUCTION OF A NANOTUBE FROM A SLAB                      *
 * Y. Noel, P. D'Arco, R. Demichelis, C.M. Zicovich-Wilson, R. Dovesi          *
 * On the use of symmetry in the ab initio quantum mechanical simulation of    *
 * nanotubes and related materials                                             *
 * J. Comput. Chem., 31, 855-862 (2010)                                        *
 *******************************************************************************

 ********** SWCNT **********
 ROLL UP DIRECTION    6   3 SPECIFIED IN INPUT IS DEFINED FOR
 A GAMMA=60 DEG CELL REF. THESE INDICES ARE TRANSFORMED INTO    9   3 IN
 A GAMMA=120 DEG CELL REF.
 ***************************

 DISTINCTIVE NANOTUBE VECTORS EXPRESSED IN THE SLAB UNIT CELL BASIS A AND B
                                                   A   B
       ROLL UP DIRECTION, R                        9   3
       HELICAL DIRECTION, H                        2   1
       CYLINDER AXIS, L                            1   5

 LENGTHS IN ANGSTROM
       CYLINDER RADIUS                             3.112
       PERIOD ALONG THE CYLINDER AXIS, |L|        11.289
       THICKNESS OF THE SLAB                       0.000

 NUMBER OF ATOMS                                  84
 CHIRAL ANGLE                                     19.107 DEG

 NANOTUBE SYMMETRY
       NUMBER OF SYMMETRY OPERATORS                   42
       NUMBER OF ATOMS IN THE ASYMMETRIC UNIT          2
       THE HELICAL OPERATOR:
         ORDER                                        42
         ANGLE                                 (  23/ 42)x2PI RAD=>197.143 DEG
         TRANSLATION                               1/ 14 OF THE CYLINDER PERIOD
                                                   0.8085 ANGSTROM

In the following part of the output, four sets of coordinates are printed for each atom in the nanotube unit cell: X/A, Y, Z and R. The first one is fractional (X is the only periodic direction), the second and third ones Cartesian. The latter stands for ''radius'' and indicates the distance of the atom from the tube axis. The importance of R will be commented later on, in the example of inorganic nanotubes.

 *******************************************************************************
 ATOMS IN THE ASYMMETRIC UNIT    2 - ATOMS IN THE UNIT CELL:   84
     ATOM              X/A             Y(ANGSTROM)         Z(ANGSTROM)   R(ANGS)
 *******************************************************************************
   1 T   6 C     0.00000000000E+00  0.31121000486E+01  0.00000000000E+00   3.112
   2 T   6 C     0.95238095238E-01  0.29738381753E+01 -0.91730759260E+00   3.112
   3 F   6 C     0.71428571429E-01 -0.29738381753E+01 -0.91730759260E+00   3.112
   4 F   6 C     0.16666666667E+00 -0.31121000486E+01  0.11102230246E-15   3.112
   5 F   6 C     0.14285714286E+00  0.25713377297E+01  0.17531083801E+01   3.112
   6 F   6 C     0.23809523810E+00  0.29738381753E+01  0.91730759260E+00   3.112
   7 F   6 C     0.21428571429E+00 -0.19403626426E+01 -0.24331377946E+01   3.112
   .
   .
   .

With the NANOJMOL keyword, a file named NANOJMOL.DAT containing the nanotube structure is generated. This file can be directly used with the Jmol 3D structure viewer (www.jmol.org).
1 - WARNING: NANOJMOL must appear in the geometry block BEFORE the NANOTUBE or SWCNT keywords.
2 - WARNING: note that by default, if polyhedra are present, Jmol will try to show them; this can be turned off in jmol (see Jmol documentation for details).

The NANOJMOL.DAT file can be, for example, directly inserted in web page as shown in the animation of FIG. 9.

Two restart keywords, NANORE and SWCNTRE, allow to build a (N1,N2) nanotube by starting from the structure of another one (a previously (N1_old,N2_old) optimised one, read from unit fort.35). The rolling direction of the two tubes must be the same. With NANORE and SWCNTRE, the ''old'' nanotube is unrolled and re-rolled according to a new R vector, with minor modifications to the structure.
1 - WARNING: the ''old'' and ''new'' R vectors must have the same direction.
2 - WARNING: NANORE and SWCNTRE are exclusive with the NANOTUBE and SWCNT keywords, respectively.

The syntax is the same for both keywords, and the difference between them is related to the choice of the hexagonal reference system of the flat slab.

Let us consider, for example, the (8,0) and the (10,0) carbon nanotubes. We have optimised the structure of the former, and we want to build the latter starting from its geometry. With SWCNTRE the (8,0) SWCNT is unrolled and re-rolled as (10,0). In order to do this, the information on geometry of both the starting slab (graphene) and the (8,0) SWCNT is required. The first one is given in input (or read with EXTERNAL from unit fort.34), the second one is read with an EXTERNAL strategy from unit fort.35. The input syntax is shown in the following table.

SLAB
77
2.47
1
6 0.333333  0.666667  0.000000
SWCNTRE
8 0      (old nanotube)
10 0      (new nanotube)

The output contains the features of both the nanotubes. The following information is printed, together with the ones previously shown.

 ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
 ooo                NANOTUBE REBUILD: BUILDING A NANOTUBE                    ooo
 ooo    STARTING FROM THE STRUCTURE OF A PREVIOUSLY ROLLED NANOTUBE          ooo
 ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo
 o
 o REFERENCE NANOTUBE
 o ==================
 o REF NANOTUBE CHARACTERISTICS:
 o ROLL UP NANOTUBE VECTOR (R) EXPRESSED IN THE SLAB UNIT CELL BASIS A AND B
 o                             A   B
 o     R                       8   0
 o     
 o LENGTHS IN ANGSTROMS
 o     CYLINDER RADIUS                             3.145
 o     PERIOD ALONG THE CYLINDER AXIS, |L|         4.269
 o     
 o  ATOMS IN THE REF ASYMMETRIC UNIT:      2
 o ATOM            X/A              Y(ANGSTROM)        Z(ANGSTROM)       R(ANGS)
 o  ****************************************************************************
 o   1    6    0.66606574844E+00  0.31968074414E+01  0.15055500000E-06   3.197
 o   2    6    0.83393425156E+00  0.29534620698E+01 -0.12233642082E+01   3.197
 o *****************************************************************************
 o
 o  ACTUAL NANOTUBE
 o ===============
 o RATIO NEW/REF   1.250
 o
 o FROM NOW, ALL CHARACTERISTICS REFER TO THE NEW NANOTUBE (  10,   0)
 o
 o  ATOMS IN THE NEW ASYMMETRIC UNIT:      2
 o ATOM            X/A              Y(ANGSTROM)        Z(ANGSTROM)       R(ANGS)
 o  ****************************************************************************
 o   1    6    0.66606574844E+00  0.37880893100E+01  0.12308249857E+01   3.983
 o   2    6    0.83393425156E+00  0.39830297894E+01 -0.15006580322E-06   3.983
 o
 ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo

Note that, in the previous example, NANORE and SWCNTRE are equivalent, as N2=0 (see FIGs. 7 and 8).

This option is particularly helpful for the geometry optimisation of inorganic nanotubes (thick slabs, large systems, the geometry of the tubes is very different from the one of the slab), as the number of optimisation steps is reduced. As an example, in the following table the effect of the NANORE strategy on the optimisation of imogolite is shown. Imogolite, Al₂(OH)₃SiO₃OH, is a tubular hydrated aluminosilicate. It was built by wrapping a hexagonal layer of gibbsite where three out of six H atoms are replaced by one Si-OH group (Si results tetracoordinated). The SiO₄ isolated tetrahedra (and, thus, the silanolic groups) are facing towards the inside of the tube, whereas AlO₆ edge-sharing octahedra form the external wall.
N_A and N_S are the number of atoms in the unit cell and the number of symmetry operators, respectively; t_SCF and t_G [s] are the cost of a single step of the SCF cycle and of the first gradient calculation, respectively (16 processors - Dual Core AMD Opteron(tm) Processor 875, 2.21GHz, RAM 16Gb). Beginning indicates the starting structure for the optimisation (the optimised slab or a previously optimised tube), OPT steps the number of SCF+gradient calculations in the optimisation process.
Note that the number of optimisation steps when NANORE is used (all but the (10,0) nanotube) is strongly reduced (by about a half).

Exercise

Build the (5,2) carbon nanotube from graphene by using SWCNT and NANOTUBE and run the SCF calculation for both. Compare the results (test_nano01 and test_nano02). Build the (5,2) carbon nanotube by starting from the bulk, save the input for Jmol and visualise the structure (test_nano03 and test_nano04).

Optimise the structure of the (8,0) SWCNT (test_nano05). Build the (10,0) SWCNT by starting from the (8,0) (test_nano06).

FIG. 10 shows that both ΔE and δE regularly decrease as n increases, as expected and in line with previous findings. Actually, a finer analysis shows that δE data belongs to three different curves, describing the (3n,0), (3n+1,0) and (3n+2,0) cases; however, differences are quite small and tend to disappear with the increasing of the nanotube diameter.

FIG. 10: ΔE and δE for (n,0) carbon nanotubes. The curves are built with n from 8 to 24, a 6-1111G* basis set and the hybrid B3LYP Hamiltonian; energies are normalised to to two carbon atoms.

Exercise

Optimise the geometry of graphene (test_nano07_5). Run the SCF of the (5,5), (6,6) and (7,7) SWCNT and build the curve total energy (with respect to graphene)vs radius. Optimise the structure of the (5,5) SWCNT and calculate the energy gain due to the optimisation (test_nano07_1, test_nano07_2 test_nano07_3 test_nano07_4).

Electronic structure

One of the most interesting properties of carbon nanotubes is the band gap: depending on the radius and rolling direction very different electronic properties can be obtained (e.g. iin principle, band gap can be tuned by choosing different combinations of n₁and n₂). The energies of the top of the valence band and the bottom of the conducting band are automatically printed in the CRYSTAL output at each SCF step. The following example refers to the (19,0) carbon nanotube (6-1111G* basis set, B3PW Hamiltonian); the values printed at the end of the SCF cycle must be considered:

 INSULATING STATE - LEVEL SHIFTER    0.60au
 TOP OF VALENCE BANDS -    BAND    228; K    1; EIG -7.6454992E-01 AU
 BOTTOM OF VIRTUAL BANDS - BAND    229; K    1; EIG -1.4087489E-01 AU
 TTTTTTTTTTTTTTTTTTTTTTTTTTTTTT 

The band gap, Eg, is obtained as follows:

Eg=E(BOTTOM)-E(TOP)-E(LEVSH)=0.0237 A.U.

where E(BOTTOM) and E(TOP) are the energies of the bottom of the virtual and the top of the valence bands, respectively, and E(LEVSH) (LEVEL SHIFTER in the output) is an optional computational parameter added to help the SCF convergence (see LEVSHIFT keyword).

In FIG. 11 the band gap of the (n,0) nanotube family is represented as a function of the radius. Information about the Fermi level, the band gap and the energy of the valence and virtual bands also can be obtained as one-electron properties (BWIDTH keyword).

WARNING: the value of the band gap is strongly dependent on the choice of the Hamiltonian.

FIG. 11: Band gap (Eg) of (n,0) nanotubes as a function of the radius (n ranges from 8 to 24). The values are obtained with a 6-1111G* basis set and the B3PW hybrid Hamiltonian. For (3n,0) nanotubes the results are only qualitative: Eg is in these cases much smaller than 0.1 eV (see band structure in FIG. 12).

FIG. 12 shows the band structure around the valence region (6-1111G* basis set, B3PW hybrid Hamiltonian) of a (3n,0) (narrow gap, zigzag), a (3n+1,0) (semi-conductors, zigzag, similar band structure as (3n+2,0)) and a (n,n) (metallic, armchair) carbon nanotubes with similar radius (about 6 Angstrom). Further details on the band structure of SWCNTs in references Samsonidze et al. (2003) and Hamada et al. (1992) - see end of this page.

FIG. 12: Band structure of a two zigzag and one armchair SWCNT. The (15,0) nanotube exhibit a very small (but non-null) band gap. The first irreducible Brillouin zone is represented, ranging from 0 to π/a, where a is the lattice parameter. Energy in eV on the ordinate axis. The degeneracy point in the (9,9) SWCNT is at k=2π/3a.

Exercise

Optimize the geometry of the (8,0) SWCNT with the PBE functional (test_nano08) and compare the resulting Eg with the one obtained in test_nano05 with B3LYP. The experimental value is around 1.5 eV.

Calculate Eg from test_nano07_2, where IS=10. Then, run properties with IS=60 and BWIDTH and calculate Eg in these conditions (test_nano07_2_bwidth). Then, compute the band structure (test_nano07_2_band).
WARNING: the basis set used in these tests is limited, results are only qualitative. Details on the band structure of metallic and semiconducting nanotubes in References Samsonidze et al. (2003) and Hamada et al. (1992).

Chrysotile: a more complex case

Structure

Chrysotile is a mineral of the serpentine group (Mg₃Si₂O₅(OH)₄, phyllosilicate). A single fibre can be considered as a lizardite single layer wrapped along a given direction. Let us start, then, from the ortho-hexagonal lizardite slab. It consists of two sheets:

• a brucite-type octahedral sheet (MgO₆ octahedra sharing edges)
• a tetrahedral sheet built of vertices-sharing SiO₄ tetrahedra forming hexagonal motif.

The two sheets share one oxygen: the apical oxygen of tetrahedra (see FIG. 13). In our convention the lizardite layer lies on the xy plane.

FIG. 13: Lizardite structure. Mg06 octahedra (in green), share edges forming a brucite-like sheet; SiO4 tetrahedra (in beige) form a hexagonal sheet. Some of the SiO4 tetrahedra are translucent to show Si-O bonding. The apical oxygens of tetrahedra link both sheets (hidden by silicon on top view). The hexagonal and the rectangular cells (the second twice larger) cells are in blue and yellow, respectively. Drag with the mouse to rotate and use buttons to change views.

Once the lizardite structure is known, two main aspects must be considered for rolling it up.

• Rolling direction and radius, depending on R
• Position of the atoms with respect to the radius

The approach to the first point is similar to the one described for carbon nanotubes, and is directly related to the choice of R. The last point is related to the non null thickness of the starting slab, and will be introduced with the chrysotile example.

Inside and outside of the cylinder

As mentioned before, chrysotile is made of two sheets (octahedra and tetrahedra). After rolling one will be the inner side and the other one the outer side of the tube. The first one will be shortened, whereas the second one will be elongated with respect to the slab, as shown on FIG. 14. When free sheets are individually studied, it appears that Mg(OH)₆ sheet (brucite) is larger than the (theoretical) SiO₄H one. This misfit, together with other factors not discussed here (see D'Arco et al. (2009) in References), induces a preferential position of each one in the tube:

• Octahedra are outside
• Tetrahedra are inside

This wrapping is called direct and the opposite one inverse.

inverse rolling	slab	direct rolling
FIG. 14: Inverse and direct chrysotile wrappings compared to the lizardite slab. Click here or on the image to open a 3D view in a new page.

Atomic coordinates from flat to cylindrical surfaces

The position of the atoms in the curved structure results obviously different from that in the flat slab. In order to show how they transform, let us start from FIG. 15, where the projection of a generic unit cell on the xy Cartesian plane (black reference vectors on the left) is represented. In case of non-null thickness of the slab, i.e. more than one layer of atoms are present, a ''third non-periodic dimension'' (along the z axis) must be considered.

The green point in FIG. 15 represents an atom in the slab unit cell with coordinates (C_x,C_y). The orientation of the nanotube cell in the flat slab is represented with red and blue dashed-lines (R and L, respectively). It is convenient to rotate the reference Cartesian axis, so that L and R are oriented along x and y, respectively (red reference vectors, on the right). The coordinates of the green atom in the new reference system are then given by the projection of C=x C_x +y C_y on R and L. In particular:

where R_x, Ry, L_x and L_y are the Cartesian components of R and L in the ''black'' reference system of FIG. 15. Note that the same quantity can be calculated as:

FIG. 15: Generic 2D unit cell in the xy Cartesian plane. When the nanotube is built, it is convenient to rotate the reference system (from the black axis on the left to the red ones on the right) in order to orientate L and R along x and y, respectively. The coordinates of the C atom are shown in the slab (Cx,Cy) and in the nanotube (PCR,PCL) reference systems. φ1 and φ2 are the angles of the C position with R and L, respectively.

FIG. 16: From the flat slab (a) to the nanotube (b). Representation of the quantities discussed in the text and reference axis. See text for details.

When the nanotube is build, the tube axis coincides with the x Cartesian axis, and the tube circumference lies in the yz plane, as illustrated in FIG. 16. Upon rolling, atoms located on a given z surface of the 2D slab are mapped on the same cylindrical surface: for example, the blue and green atoms in FIG. 16(a) are on the same surface in FIG. 16(b). For z=0, distances measured on the ''flat'' plane and the ones on the corresponding cylindrical surface (than on the arc) are equal. This cylinder defines the so called neutral surface: its distance R_S from the cylinder axis is such that 2πR_S equals the norm of the rolling vector, |R|.

In FIG. 16 the angle between the y axis and the position of the green atom (with (P_CR,P_CL, δ) coordinates in the nanotube flat unit cell) in the curved structure, φ₃, is proportional to the ratio between the projection of C on R and |R|.

The coordinates of the atom in the curved structure are then easily obtained:

As a result, distances among atoms on cylindrical planes at z=+|δ| increase, and distances at z=-|δ| decrease, the expansion and shrink being proportional to δ; radial distances (orthogonal then to the rolling vector) remain constant.

During optimisation atoms can move in any direction, and then also radially with respect to the tube axis (see FIG. 17). At the end of the optimisation the distance of each atom from the neutral surface can be different from the one just after rolling (see the R(ANGS) coordinate printed for eaxh atom with the X/A, Y, Z set of coordinates shown in a previous table). Information on the radial relaxation can be obtained by comparing for the various tubes the distance of a given atom from the tube axis with the radius of the neutral surface.

FIG. 17: Position of the wrapped slab with respect to the neutral surface for the (14,-14) chrysotile. Click here or on the image to open a 3D animated view in a new page.

WARNING: the choice of the position of the nautral surface can influence the SCF convergence, the number of optimisation steps and the relaxation mechanism of the structure. For example, in FIG. 17, if we place the neutral surface on the external wall of the tube, the structure results very compressed, bond distances are strongly reduced and a larger number of optimisation steps would be required. By default, the nautral surface is placed in the middle of the flat slab.

Input and Output

The following table contains the input for the construction of a (14,-14) chrysotile fibre starting from the lizardite slab.

chrysotile 14,-14 SLAB 1 5.31004945 5.31004945 120.000000 18 1 -3.351279436683E-01 0.000000000000E+00 2.637304577944E+00 1 -1.541822225448E-13 -3.351279436686E-01 2.637304577944E+00 1 3.351279436684E-01 3.351279436686E-01 2.637304577944E+00 8 -3.334648646031E-01 0.000000000000E+00 1.710965746561E+00 8 -1.534050664276E-13 -3.334648646035E-01 1.710965746561E+00 8 3.334648646033E-01 3.334648646035E-01 1.710965746561E+00 12 3.314444744801E-01 0.000000000000E+00 6.708390562340E-01 12 1.524891324323E-13 3.314444744804E-01 6.708390562340E-01 12 -3.314444744802E-01 -3.314444744804E-01 6.708390562340E-01 8 0.000000000000E+00 0.000000000000E+00 -4.076681137970E-01 8 3.333333333330E-01 -3.333333333337E-01 -4.010114366710E-01 8 -3.333333333330E-01 3.333333333334E-01 -4.010114366710E-01 1 0.000000000000E+00 0.000000000000E+00 -1.357208074319E+00 14 3.333333333330E-01 -3.333333333337E-01 -2.033436938287E+00 14 -3.333333333330E-01 3.333333333334E-01 -2.033436938287E+00 8 4.814202807694E-01 0.000000000000E+00 -2.606101339892E+00 8 -2.385869279919E-13 4.814202807689E-01 -2.606101339892E+00 8 -4.814202807692E-01 -4.814202807689E-01 -2.606101339892E+00 NANOTUBE 14 -14 TEST END

In order to inverse the rolling (swap inside and outside the cylinder), the ROTCRY keyword must be used (see manual for keyword description).

ROTCRY MATROT 1 0 0 0 1 0 0 0 -1 NANOTUBE 14 -14

To change the position of the atoms with respect to the neutral surface, ATOMROT keyword can be used (see manual for details). A shift of +0.37 Angstroms along the z axis will be applied to all atoms in the following example.

ATOMROT 0 -1 999 0. 0. 0.37 NANOTUBE 14 -14

The general strategy of a geometry optimisation, in CRYSTAL, is the following. The classification of integrals to be evaluated or approximated is performed on the starting geometry (at the first cycle) and kept fix during the optimisation. Once the convergence on the structure is reached, the integral classification is performed according to the new geometry: if the structure is not in the minimum, the optimisation restarts (see FINALRUN keyword).

When a large number of atoms is present in the unit cell, as in the case of chrysotile, and the geometry varies a lot during the optimisation, this strategy leads to a very large number of optimisation cycles (more than 100 in the case of chrysotile). The FIXDELTE keyword defines a threshold on the total energy change for redefining the geometry to which the integral classification is referred. This tool, together with an accurate technique to control the step size (Trust Radius, default choice), allows to faster and more efficiently reach the minimum on the energy surface.

OPTGEOM FULLOPTG BFGS FIXDELTE 3 FINALRUN 4 END

Block of instruction for geometry optimisation of chrysotile.

Exercise

Build the (14,-14) chrysotile inversely and directly wrapped. Visualise the structure with Jmol (test_nano11 and test_nano12).

Time scaling

In the CRYSTAL code, the helical symmetry of nanotubes is exploited at three levels. We have already shown the automatic generation of the nanotube starting from a two-dimensional structure in the input section. The exploitation of point symmetry for the calculation of the mono- and bi-electronic integrals and for the diagonalisation of the Fock matrix allows to drastically reduce the computational cost of nanotube simulations.

• Point symmetry is exploited to reduce to an irreducible set the number of Fock matrix elements to be calculated, so that only the corresponding irreducible subset of mono- and bi-electronic integrals is evaluated. Once calculated, the irreducible elements of the Fock matrix are rotated by the symmetry operators of the point group to generate the full Fock matrix.
• Point symmetry is exploited in the diagonalisation of the Fock matrix, F(k), in the irreducible part of the first Brillouin zone. By transforming the Bloch basis into a symmetry adapted basis (SACO), F(k) becomes block-diagonal, each block corresponding to an irreducible representation (IR) of the point group. In the present case, where the system is periodic only along one direction (x in our convention) and the rotation part of the symmetry operators acts on the Cartesian coordinates (y and z) orthogonal to the translation direction, the point group of each k point is the same, and contains the full set of point operators. The Fock matrix factorizes then exactly in the same way at each k point. Being the group Abelian, the number of IRs (and then of blocks to be diagonalised) is equal to the number of roto-translation operators.

The effect of such a scheme on the calculation cost is documented in the following figure and table, for the case of (n,0) carbon nanotubes (B3LYP, 6-1111G*) and chrysotile (B3LYP, 6-31G*), respectively (single processor - Intel Xeon, 1.86GHz, RAM 8Gb).

In the case of carbon nanotubes, in going from (8,0) to (24,0), the number of atoms in the unit cell, NA, basis functions in the unit cell, NAO, and symmetry operators, NS, increases exactly by a factor three (from 32 to 96, from 704 to 2112 and from 16 to 48, respectively). The complete Fock matrix for the largest case is then about three times larger than for the (8,0) case, whereas the irreducible Fock matrix remains roughly constant (about 11'000 elements); the ratio between the size of these two matrices, NR, is close to the number of symmetry operators.

The time (t, seconds) required for the calculation of bi-electronic integrals, the most expensive part of the SCF calculation, decreases for a while (due to interactions across the tube, not negligible for small radii) and then remains nearly constant. The gradient, not shown in the figure, exhibit same behaviour: the cost from (12,0) to (24,0) varies by less than 1%. The overall cost of the SCF cycle increases (slowly) for three reasons.

• The diagonalisation step scales linearly with the number of symmetry operators.
• The AO→Bloch→SACO basis transformations and the back transformations to the AO basis for building the density matrix are matrix multiplications of sparse matrices and time scaling is close to N_AO².
• The overhead for the symmetry analysis increases with the number of symmetry operators.

Overall, the total cost increases by less than 50% for the SCF and remains constant for the gradient, whereas in a perfect linear scaling code, not exploiting the symmetry, the cost would increase by a factor three.

Regarding chrysotile, it is important to notice that the first ab initio simulation was performed with the CRYSTAL code. The (14,-14), (19,-19), (24,-24) nanotubes and the lizardite slab (in principle, a nanotube with an infinite radius and the C₁ operator only) are considered for comparison. In going from the slab to the largest nanotube, N_A increases by a factor 48 (N_S of the (24,-24) nanotube), and so is for the number of N_AO. The N_R ratio (the irreducible Fock matrix contains about 59'000 elements), is extremely close to N_S.

	R	NA	NAO	NR	NS	tB	tF+P	tDFT	tSCF	tG
[001] slab	infinite	18	236	1.0	1	949.0	0.8	80.5	1033.4	9801.4
(14,-14)	20.493	504	6'608	27.7	28	1210.1	98.8	111.3	1453.2	11435.7
(19,-19)	27.812	684	8'968	37.6	38	1272.5	205.4	115.1	1637.7	11690.8
(24,-24)	35.131	864	11'328	47.5	48	1330.0	385.4	129.2	1900.8	12024.6

(R: Angstrom, t: seconds)

Times in table refer to a single step of the SCF cycle. t_B is the cost of the evaluation of the bi-electronic integrals. It increases by a factor 1.5 when passing from the slab to the (24,-24) nanotube. t_F+P and t_SCF (time for F diagonalisation plus density matrix construction and cost of a single SCF step) increase for the same reasons discussed for carbon nanotubes. t_SCF increases by about a factor 2 in going from the slab to the (24,-24) nanotube; in a perfectly linear scaling code, not exploiting the symmetry, the factor would be 48. t_G (cost of a gradient calculation) is nearly constant for the nanotubes, whereas it increases by a factor 1.2 when passing from the slab to the (14,-14) nanotube. The time required by the calculation of the DFT contribution, t_DFT, remains nearly constant.

Small collection of CRYSTAL inputs and outputs

Single-Walled Carbon Nanotubes

test_nano01.d12	test_nano01.out	Single point calculation of (5,2) nanotube from graphene with SWCNT keyword
test_nano02.d12	test_nano02.out	Single point calculation of (5,2) nanotube from graphene with NANOTUBE keyword
test_nano03.d12	test_nano03.out	(5,2) nanotube from bulk with SWCNT keyword, generation of a Jmol file
test_nano04.d12	test_nano04.out	(5,2) nanotube from bulk with NANOTUBE keyword, generation of a Jmol file
test_nano05.d12	test_nano05.out	Geometry optimisation of the (8,0) zigzag SWCNT
test_nano05_ppan.d3	test_nano05_ppan.outp	Mulliken population analysis of the optimised (8,0) SWCNT, from test_nano05
test_nano05_band.d3	test_nano05_band.outp	Band analysis of the optimised (8,0) SWCNT, from test_nano05
test_nano06.d12 test_nano06.f35	test_nano06.out	Geometry optimisation of the (10,0) SWCNT starting from the (8,0) nanotube (SWCNTRE: rebuild from test_nano05)
test_nano06_external.d12 test_nano06_external.f35 test_nano06_external.f34	test_nano06_external.out	Same as test_nano06, with external input for slab
test_nano07_1.d12	test_nano07_1.out	Single point calculation of the (5,5) armchair SWCNT
test_nano07_2.d12	test_nano07_2.out	Geometry optimisation of the (5,5) SWCNT
test_nano07_2_bwidth.d3	test_nano07_2_bwidth.outp	Band gap of the (5,5) SWCNT from test_nano07_2: the special point of the first Brillouin Zone where the valence and the virtual bands are overlapped is included in the calculation
test_nano07_2_band.d3	test_nano07_2_band.outp	Band structure of the (5,5) SWCNT from test_nano07_2
test_nano07_3.d12	test_nano07_3.out	Single point calculation of the (6,6) armchair SWCNT
test_nano07_4.d12	test_nano07_4.out	Single point calculation of the (7,7) armchair SWCNT
test_nano07_5.d12	test_nano07_5.out	Geometry optimisation of graphene for energy comparison with test_nano07_1, test_nano07_2 test_nano07_3, test_nano07_4
test_nano08.d12	test_nano08.out	Geometry optimisation of the (8,0) nanotube with PBE functional
test_nano09.d12 test_nano09.f34	test_nano09.out	Frequency calculation of the (8,0) SWCNT optimised in test_nano05
test_nano10.d12 test_nano10.f34	test_nano10.out	First order polarizability of the (8,0) SWCNT optimised in test_nano05

Chrysotile

test_nano11.d12	test_nano11.out	Geometry of the (14,-14) inversely wrapped chrysotile
test_nano12.d12	test_nano12.out	Geometry of the (14,-14) directly wrapped chrysotile

For higher quality basis sets click here.

References

Y. Noel, P. D'Arco, R. Demichelis, C. M. Zicovich-Wilson and R. Dovesi
On the use of symmetry in the ab initio quantum mechanical simulation of nanotubes and related materials
J. Comput. Chem., 31, 855-862 (2010)

i

P. D'Arco, Y. Noel, R. Demichelis and R. Dovesi
Single-layered chrysotile nanotubes: a quantum mechanical ab initio simulation
J. Chem. Phys., 131, 204701 (2009)

C. T. White, D. H. Robertson and J. W. Mintmire
Helical and rotational symmetries of nano-scale graphitic tubules
Phys. Rev. B, 47, 5485-5488 (1993)

E. B. Barrosa, A. Joriob, G. G. Samsonidzef, R. B. Capazc, A. G. Souza Filhoa, J. Mendes Filhoa, G. Dresselhaus and M. S. Dresselhaus
Review on the symmetry-related properties of carbon nanotubes
Phys. Rep., 431, 261-302 (2006), and references therein.

Ge. G. Samsonidze, R. Saito, A. Jorio, M. A. Pimenta, A. G. Souza Filho, A. Grüneis, G. Dresselhaus and M. S. Dresselhaus
The concept of cutting lines in carbon nanotube science
J. Nanoscience and Nanotechnology, 3, 431-458 (2003).

N. Hamada, S. Sawada, and A. Oshiyama
New one-dimensional conductors: graphitic microtubules
Phys. Rev. Lett., 68, 1579-1581 (1992).

---