Y. Noel and R. Demichelis
For a correct visualisation of
figures and animations present in this tutorial
we advise the use of Mozilla Firefox or Safari
must be enabled on your machine.
|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|
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= n1a1 + n2a2 (where a1 and a2 are the slab cell vectors). The (n1,n2) 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=4a1+3a2).
If n1 and n2 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 = l1a1 + l2a2, 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=h1a1 + h2a2,
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(vi,vj) is the surface defined by the vi and vj 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
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:
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 |n1h2-n2h1|=2. It is convenient to choose H as the shortest vector in the unit cell satisfying S(R,H)=N S(a1,a2) (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 HR and HL.
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
Graphene mapped to build a (4,2) SWCNT.
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
basis vectors, three vectors define:
|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).|
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.
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 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:
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.
|•||N1,N2||Components of the R roll-up vector. N1 and N2 are integers (indicated as n1 and n2 in the previous equations). See previous section for details.|
This keyword builds a nanotube from a slab. The components of the
vector are expressed in the slab basis vectors.
1 - WARNING: in CRYSTAL if the 2D lattice is hexagonal, γ, the angle between a1 and a2, 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/a1 and y/a2 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.
|•||N1,N2||ONLY FOR HEXAGONAL SYSTEMS. Components of the R roll-up vector expressed in a γ=60° lattice (instead of γ=120°). N1 and N2 are integers. See previous section for details.|
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°
&gamma=60° basis is the following:
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.
|E.g.: Nanotube from 3D
0 0 0
6 0.00000 0.00000 0.00000
6 0.33333 0.66667 0.00000
0 0 1
E.g.: Nanotube from slab I
6 0.33333 0.66667 0.00000
E.g.: Nanotube from slab II
6 0.33333 0.66667 0.00000
E.g.: Nanotube from slab III
2.47 2.47 60.00
6 0.33333 0.33333 0.00000
6 0.66667 0.66667 0.00000
OPTIONAL GEOMETRY KEYWORDS
OPTIONAL GEOMETRY KEYWORDS
OPTIONAL GEOMETRY KEYWORDS
OPTIONAL GEOMETRY KEYWORDS
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
containing the nanotube structure is generated.
This file can be directly used with the Jmol 3D structure viewer
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.
|FIG. 9: (6,3) SWCNT structure generated by the NANOJMOL keyword. Drag with the mouse to rotate.|
Two restart keywords, NANORE and SWCNTRE, allow to
build a (N1,N2) nanotube by starting from the structure of another one (a previously
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.
|NANORE (or SWCNTRE)|
|•||N1_old,N2_old||Indices of the starting nanotube.|
|•||N1,N2||New indices of the rolling vector.|
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, Al2(OH)3SiO3OH, 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 SiO4 isolated tetrahedra (and, thus,
the silanolic groups) are facing towards the inside of the tube, whereas
AlO6 edge-sharing octahedra form the external wall.
NA and NS are the number of atoms in the unit cell and the number of symmetry operators, respectively; tSCF and tG [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).
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).NANORE (or SWCNTRE) strategy is very small.
|Carbon nanotubes: a simple case|
The stability of nanotubes with respect to
the corresponding flat surface (ΔE) and the relaxation
energy after the rigid rolling (δE)
allow to deepen the knowledge of the systems and
the relationship between structure and properties of
the matter, as well as
bring forward hypotheses on possible applications.
In the case of carbon nanotubes, it is clear that a graphene layer will always be more stable than the most stable carbon nanotube, whereas in other cases, for example in some naturally observed inorganic nanotubes, the wrapped system should result more stable than the corresponding slab.
The example of the (n,0) carbon nanotube family - namely zigzag - is here presented. We want to analyse ΔE and δE as a function of the nanotube radius (and so, of the n index). Let us consider n from 8 to 24. In order to calculate ΔE, we need the equilibrium energy of graphene, E(slab), and of the nanotubes, E(nano) (the latter normalised to two C atoms).
We run, then, geometry optimisations for graphene and the selected nanotubes. δE is obtained as the difference between the energy corresponding to a rigid rolling of graphene, E(0) (first SCF of the optimisation), and the energy corresponding to the optimized structure, E(nano):
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.|
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).
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 n1and n2). 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:
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.|
Calculate Eg from test_nano07_2,
where IS=10. Then, run properties with IS=60 and BWIDTH and calculate Eg in these conditions
Then, compute the band structure
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|
Chrysotile is a mineral of the serpentine group (Mg3Si2O5(OH)4, 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:
|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.
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.
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)6 sheet (brucite) is larger than the (theoretical) SiO4H 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:
|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.|
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 (Cx,Cy). 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 Cx +y Cy on R and L. In particular:
where Rx, Ry, Lx and Ly 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 RS from the cylinder axis is such that 2πRS 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 (PCR,PCL, δ)
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.
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.|
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.
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.
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 C1 operator only) are considered for comparison. In going from the slab to the largest nanotube, NA increases by a factor 48 (NS of the (24,-24) nanotube), and so is for the number of NAO. The NR ratio (the irreducible Fock matrix contains about 59'000 elements), is extremely close to NS.
Times in table refer to a single step of the SCF cycle. tB 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. tF+P and tSCF (time for F diagonalisation plus density matrix construction and cost of a single SCF step) increase for the same reasons discussed for carbon nanotubes. tSCF 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. tG (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, tDFT, remains nearly constant.
|Small collection of CRYSTAL inputs and outputs|
|Single point calculation of (5,2) nanotube from graphene with SWCNT keyword|
|Single point calculation of (5,2) nanotube from graphene with NANOTUBE keyword|
|(5,2) nanotube from bulk with SWCNT keyword, generation of a Jmol file|
|(5,2) nanotube from bulk with NANOTUBE keyword, generation of a Jmol file|
|Geometry optimisation of the (8,0) zigzag SWCNT|
|Mulliken population analysis of the optimised (8,0) SWCNT, from test_nano05|
|Band analysis of the optimised (8,0) SWCNT, from test_nano05|
|Geometry optimisation of the (10,0) SWCNT starting from the (8,0) nanotube (SWCNTRE: rebuild from test_nano05)|
|Same as test_nano06, with external input for slab|
|Single point calculation of the (5,5) armchair SWCNT|
|Geometry optimisation of the (5,5) SWCNT|
|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|
|Band structure of the (5,5) SWCNT from test_nano07_2|
|Single point calculation of the (6,6) armchair SWCNT|
|Single point calculation of the (7,7) armchair SWCNT|
|Geometry optimisation of graphene for energy comparison with test_nano07_1, test_nano07_2 test_nano07_3, test_nano07_4|
|Geometry optimisation of the (8,0) nanotube with PBE functional|
|Frequency calculation of the (8,0) SWCNT optimised in test_nano05|
|First order polarizability of the (8,0) SWCNT optimised in test_nano05|
|Geometry of the (14,-14) inversely wrapped chrysotile|
|Geometry of the (14,-14) directly wrapped chrysotile|
For higher quality basis sets click here.
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)
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).