1 Introduction
Nowadays, the atomic theory developed by Greek philosophers Leucippus and Democritus has become general knowledge and the most fruitful scientific idea. Thanks to the rapidly developing science and technology, even the manipulation of a single atom has become possible. The accumulation of knowledge and technology has finally led to the nanoscale world. Nanotechnology that fabricates and characterizes materials with a size several hundred nanometers or less has emerged as one of the frontiers in science and technology. Different kinds of nanomaterials, such as carbon nanotubes (CNT) and quantum dots, can now be routinely made. A better and deeper understanding – especially from the first principles – of the properties of nanomaterials is highly desired for the future development of nanotechnology.
The conventional quantum chemical methods, such as Hartree–Fock (HF) method and density functional theory (DFT), have achieved great success in predicting molecular properties. However, when we directly apply these methods to large-scale systems like nanomaterials, the computational time required becomes unacceptable. This is mainly due to the fact that these methods have O
or worse scaling behavior. Computational methods with linear-scaling behavior have attracted great attention in recent years since the first proposal of the so-called “divide-and-conquer” scheme [
1,
2]. Recently, we have proposed a “central insertion scheme” (CIS) [
3] which also follows the spirit of “divide-and-conquer” scheme. This method takes advantage of the fact that the interaction between different units in the middle of a finite periodic system should be converged if the system is large enough. These units in the middle thus become identical, so that one can easily elongate the initial system by adding the identical units in the middle of the system continuously. Since the elongated system is obtained by inserting new units in the middle of an already very large initial system, the so-called self-consistent field (SCF) procedure for the elongated system can generally be eliminated. CIS can save significant computational time without loss of accuracy. It was successfully applied to study electronic structures of sub-120 nm long conjugated polymers, sub-60 nm long single-walled carbon nanotubes (SWCNTs), and sub-30 base pairs DNA segments at the DFT level [
3].
By employing advanced numerical methods, such as the implicitly restarted Arnoldi method (IRAM) [
4–
7] and spectral transformation to solve the large-scale eigenvalue problems, we have further improved the capabilities of CIS [
8] and made an efficient parallel program package named “BioNano-LEGO” [
9]. Such improvements in CIS enable us to treat very large-scale systems at the first principles level with good accuracy and very low computational cost. By using CIS and the tool package, we have also studied electronic structures of different nanomaterials, such as sub-100 nm long semiconductor SWCNTs [
10], hydrogen-terminated nanodiamond clusters up to 7.3 nm (20 959 carbons) in diameter [
8]. These nanomaterials are probably the largest systems that have ever been calculated at the first principles level. Furthermore, by the combination of CIS and properties calculations, we have successfully investigated the electron transport properties of polymers, SWCNTs and DNA [
3], the X-ray spectroscopies of SWCNTs with different diameters and lengths [
10]. Such results have made it possible to interpret the existing experimental observations, and well demonstrated the applicability of the CIS.
In this contribution, we will first present the theoretical background and basic idea of CIS (Section 2), followed by the approximations within CIS, and the corresponding error control (Section 3). As mentioned before, what we finally encounter is a large-scale eigenvalue problem. Such a problem requires special numerical techniques like IRAM [
4–
6] and spectral transformation, which will be discussed in Section 4. The applications of CIS are given in Section 5. Finally, we draw our conclusions in Section 6.
2 CIS procedure
In general, the most time consuming part in conventional computational methods is the SCF procedure to obtain the accurate Fock/Kohn–Sham matrix for the investigated system. Due to the enormous number of electrons involved in the large-scale systems, the number of integrals needed to be computed in each SCF iteration is extremely huge, and become the bottleneck for the investigation of large-scale systems. It could be expected that the computational time can be significantly reduced if the number of SCF iterations decreases. One feasible way to reduce the number of SCF iterations is to obtain a good initial guess which is close to the final result. This is what CIS aims to achieve: providing a better initial guess.
The basic idea of the CIS procedure is quite straightforward. Take one-dimensional system as an example, we first do the traditional electronic structure calculations for an initial system as shown in Fig. 1(a), which should be at least partly periodic and large enough to have well-converged middle parts. We then make another system which is identical with the first one but with exactly one unit mismatched, as described in Fig. 1(b). This could usually be done by performing the so-called screw operator on the first system:
where and are the position vectors of atoms in the first and second systems, respectively. and are the rotation and translation operators along a symmetrical axis of the system.
Notice that units and are exactly the same, then a new and large system containing units can be obtained by combining units L, U1, (or ), and as displayed in Fig. 1(c). This process can also be viewed as inserting a new unit at the center of the first system.
The way of constructing the Fock/Kohn–Sham matrix of the new system follows the same process as for the geometrical structure. As shown in Fig. 1(d), we divide the Fock/Kohn–Sham matrices of the first and second systems into sub-matrices, denoted as and , respectively. Then we enlarge these two matrices by adding zero sub-matrices to the end parts of each matrix. Thus, by putting these two systems together with one unit mismatched, a new matrix with the structure of Fig. 1(e) is obtained for the new system. The new sub-matrix (i and j run over units , and ) can be computed through the following equation:
where
and
are the sub-matrices of the first and second systems. Moreover, it can be proved that [
8] the Fock/Kohn–Sham matrix
of the second system can be obtained by several matrix multiplications,
, without redoing electronic structure calculations. Here,
F and
S are the Fock/Kohn–Sham and overlap matrices of the first system,
is the overlap matrix of the second system, and
RC is the rotation matrix for coefficient matrices of molecular orbitals.
The parameter α is a weighting factor depending on the system under investigation. At present, we determine it according to the distance di between unit i and the center of mass of the system,
where d0 is a threshold depending on the investigated system, and ε is a small fraction number.
We should mention that the CIS procedure does not restrict in one-dimensional case, the central insertion processes of geometrical structure and the Fock/Kohn–Sham matrix can be applied along any different direction [
8]. By repeating this central insertion process, even larger systems can be constructed. With these new Fock/Kohn–Sham matrices one can obtain molecular orbital energies and wave functions of the large-scale systems with high accuracy.
3 Approximations and error control
During the CIS procedure, there are some matrix elements that cannot be generated by the combination of two sub-systems as shown in Fig. 1(e). Notice that these elements represent the long-range (e.g.,
Å) interaction between the two ends [
and
as described in Fig. 1(e)] of the system, we can thus simply set them as zero. This approximation will result in the error
during CIS calculations. Another important approximation made in CIS procedure is that the central parts of a large enough initial system are assumed to have converged electronic structures. For instance, in the middle of the system as described in Fig. 1(a), the wave function, charge density as well as the Hamiltonian of the unit
is approximately the same as that of unit
[
3,
8]. We can thus approximate the interaction between two units with short range in the new large system [Fig. 1(c)] as that of initial systems [Fig. 1(a) and (b)]. We name this kind of error as
.
Let F be the true Fock/Kohn–Sham matrix of the new large system and be the Fock/Kohn–Sham matrix obtained from the CIS procedure, then the error matrix is
where is the total error.
Let
and
be the eigenvalue and eigenvector of
, and the Euclidean norm of vector
satisfies
. Let
be the residual vector. Then there exists an eigenvalue
λ of
F such that [
11]
where is the spectral radius of .
An approximated estimation of
can be obtained using Gerschgorin’s theorem [
11]:
where is the infinite norm of matrix .
From Eq. (6), it can be concluded that if a large enough system whose central parts have converged interaction energies among them has been chosen as the initial system, the insertion of a new unit at the center of the initial structure would hardly affect the convergence and consequently reduce remarkably the errors caused by the two aforementioned approximations. This will eventually ensure that the error will not be cumulative with increasing the size of the system, which is the whole purpose of our error control. Therefore, with the error control, one can always find the best initial system for CIS procedure to ensure the accuracy of large-scale systems.
Furthermore, we have also proposed several numerical procedures to verify these two approximations and to choose a suitable initial system. These procedures can be found in Refs. [
3,
8].
4 Numerical techniques and codes
Notice that the new large system is obtained by inserting new units at the center of an already very large initial system, the SCF steps for the large system can thus be eliminated whilst keeping high accuracy. What one only needs to consider is the so-called generalized eigenvalue problem
where λ and c are the eigenvalue and eigenvector, respectively.
Obtaining all eigenvalues and eigenvectors of Eq. (7) is time-consuming and even not feasible for large-scale systems. Moreover, we are usually just concerned with a few eigenvalues around a specific number
σ. For that, we can use the spectral transformation [
6]:
Thus,
which has transferred Eq. (7) into the standard eigenvalue problem.
For large-scale eigenvalue problems, there usually exist three main techniques: Krylov subspace methods, moment-matching approaches and Padé approximations [
7]. At present, we have used a category of Krylov subspace methods – the implicitly restarted Arnoldi method [
4–
6] which reduces the storage and computational requirements a lot, and is suitable for large-scale eigenvalue problems.
The CIS has been implemented in an efficient tool package named “BioNano-LEGO” [
9] using the above numerical methods and efficient parallelization techniques. This tool package is divided into four steps:
(1) Obtain the overlap and Fock/Kohn–Sham matrices S0 and F0 of the initial system, and calculate .
(2) Compute the Fock/Kohn–Sham matrix Fi of the ith translated and rotated system: , where . Si is the overlap matrix of the ith system, and is the rotation matrix for coefficient matrices of molecular orbitals of the i-th system.
(3) Generate the geometrical structure of the large system by using “central insertion scheme”; combine the Fock/Kohn–Sham matrices Fi () according to Eq. (2) to obtain the Fock/Kohn-Sham matrix of this large system; compute the overlap matrix SN of the large system.
(4) Compute a few interested eigenvalues around a specific number σ of the general eigenvalue problem by using the IRAM and spectral transformation.
5 Applications
In this section, we demonstrate the developed CIS at work by investigating electronic structures of several different nanomaterials – carbon nanotubes, semiconductor nanoclusters and metal doped nanocluster. These materials play an important role not only in the scientific field but also in the industrial world. The understanding of their properties at a high theoretical level thus becomes essential for their further applications. As will be demonstrated below, the developed CIS provides us an efficient way to understand the behaviors of large-scale nanomaterials at the first principles level.
5.1 SWCNT(10, 0)
Carbon nanotubes which were discovered in 1991 [
12] have attracted much attention because of their promising mechanical and electrical properties. From the theoretical side, many studies have contributed to the understanding of structures and properties of single-walled carbon nanotubes with infinite length [
13–
19]. However, from the experimental side, the lengths of synthesized SWCNTs are usually finite with an average length of several hundred nanometers [
20]. Recently, the 10 to 50 nm long SWCNTs devices [
21,
22], sub-20 nm short channel carbon nanotube transistors [
23] and a 40-nm-gate carbon nanotube field-effect transistor (CNFET) [
24] have been fabricated experimentally. These finite SWCNTs possess unique physical and chemical properties that cannot be thoroughly explained by the infinite model.
It is therefore desirable to have a good theoretical understanding of the properties of finite SWCNTs. Moreover, as revealed by recent experimental [
25–
27] and theoretical [
17] studies, highly correlated first principles methods are essential to provide accurate band-gap energies as well as exciton binding energies. However, these finite but long SWCNTs present a great challenge for conventional computational methods, due to the involvement of a huge number of electrons and the breakdown of periodic boundary conditions. Our developed computational approach, CIS, however could overcome such difficulities. By using CIS, we have successfully investigated the electronic and optical properties of SWCNTs (5, 5) [
3,
8,
10], (6, 5) and (8, 3) [
10] up to 200 nm at B3LYP/6-31G level. In this study, we will extend the applications of CIS to a zigzag SWCNT, (10, 0).
The structure of SWCNT(10, 0) is obtained from a non-orthogonal tight-binding model [
28,
29]. The dangling bonds have been saturated by hydrogen atoms and optimized at the PM3 level [
30]. In order to obtain a large initial system and keep high accuracy, we have adopted the following steps during calculations: we first chose a model C
220H
20 as the initial system and performed quantum chemical calculation for it at B3LYP/3–21G level [
30]. By performing the CIS calculation, we obtained the electronic structure of a larger SWCNT(10, 0), C
620H
20. Then, we used the wave functions of C
620H
20 from CIS calculation as an initial guess, and fed them into the Gaussian 03 program to do further SCF calculations. We obtained the converged wave functions immediately after three SCF iterations. The total energy changes from -23 510.146 to -23 510.149 Hartree, which confirms the results from CIS calculation, are indeed a good initial guess and very close to the final result.
Next, we performed CIS calculations by using C
620H
20 as an initial system to keep higher accuracy. The CIS calculated length-dependent band-gap energies, molecular orbitals around the band gap, and length-dependent density of states (DOS) of SWCNT(10, 0) are displayed in Fig. 2 (a), (b), and (c), respectively. Similarly with the case of SWCNT(8, 3) [
10], we can find gap states in SWCNT(10, 0) as denoted as red dot in Fig. 2(b). These states are mainly arisen from the hydrogen-terminated edges, and their contributions to the DOS get smaller and smaller when the length of the tube increases. At the length of 100 nm, we can hardly see any contributions from these gap states. Hence, the true band gap of (10, 0) as indicated in the DOS distributions is controlled by the energy difference between HOMO-3 (highest occupied molecular orbital, HOMO) and LUMO+3 (lowest unoccupied molecular orbital, LUMO).
Moreover, we can find that with the increase of the length, the band-gap energy of SWCNT(10, 0) decreases and reaches a finite gap energy limit. The band-gap energy of SWCNT(10, 0) at the length of 100 nm is 1.098 eV. A recent experimental work [
31] has found that the optical gap of SWCNT(10, 0) is 1.077 eV. By further considering the empirical relationship between the exciton binding energy and nanotube diameter (in nanometers),
[
26], the estimated experimental band-gap energy could be 1.505 eV. The difference between our calculated and experiment results is only around 0.4 eV. This deviation might come from the exciton model, density functional or basis sets that we have adopted.
5.2 Semiconductor nanoclusters
The main challenge nowadays in the semiconductor industry is to demonstrate the control and utilization of devices at the nanoscale with smaller dimensions. Knowledge of electronic structures of semiconductor nanoclusters is thus of crucial importance for the development of the semiconductor industry [
32–
42]. Other than bulk materials, the electronic energy band structure in a nanocluster is characterized by a series of discrete energy states, which are governed by the quantum confinement effect, which states that the band gap should increase with decreasing particle size. This information has been very useful for understanding the electronic structures of nanostructures and for the design of nanoscale devices. Theoretical studies of semiconductor nanoclusters are even more difficult than that of carbon-based nanomaterials, considering the enormous number of electrons involved. We will show here that with the CIS method, it becomes possible to obtain accurate first principle descriptions of many core materials in the semiconductor industry, such as Si and Ge nanoparticles, GaAs and InAs quantum dots, and GaN nanowires.
Si nanoparticles under investigation have the spherical diamond structure with hydrogen termination, as shown schematically in Fig. 3. Structures of these spherical particles consist of several inner shells of silicons and an outer shell of hydrogens, which were labeled as
n-shell Si nanoclusters (Si
n). The size of Si
n is represented by its diameter (
R), which follows the relationship
nm. Small clusters, such as Si4, Si6, Si8, and Si10, are calculated using the Gaussian 03 program at the B3LYP/LanL2DZ level. The computed wave functions for HOMO and LUMO of Si10 (3.4 nm) are plotted in Fig. 3(a). It is found that both orbitals are mainly localized inside the core of the cluster, which is consistent with previous findings [
43]. Based on the electronic structure of Si10, the CIS method has been applied to study larger clusters, Si12, Si14, Si16, and Si18. The largest cluster, Si18, is 6.5 nm in diameter and consists of 4 313 silicon atoms. The computed wave functions for HOMO and LUMO of Si12 and Si18 are shown in Fig. 3 (b) and (c), respectively, which are consistent with those of Si10, indicating that these gap orbitals are always confined within the core of the cluster.
Size dependent energies of HOMO and LUMO for Si clusters are given in Fig. 4(a). It seems that both HOMO and LUMO orbitals are converged when the system size goes beyond 5 nm. The energy gap between HOMO and LUMO (
Eg) has also been plotted as a function of system size in Fig. 4(b). It decreases from 5.33 eV at 1.1 nm to 2.59 eV at 6.5 nm. Moreover, our calculated energy gaps are found to be consistent with previous reported experimental [
34–
37] and theoretical [
32,
38–
42] results. We should mention that, the energy gap has not converged completely at the size of 6.5 nm, although the change of energy gap gets smaller for bigger clusters. The final conclusion about the convergence of gap energy will be obtained from calculations of even larger Si nanoparticles which are still in progress.
The same approach has been applied to Ge nanoparticles, GaAs quantum dots, and GaN nanowires. Similar quantum confinement effect has been found. for instance, we have calculated large size GaAs quantum dots up to a 38-shell (16.7 nm in diameter) at the B3LYP/LanL2DZ level. Size-dependent core-HOMOs and core-LUMOs (HOMO and LUMO with wave function localized inside the core of the particles), and energy gaps of GaAs nanoclusters are shown in Fig. 5. Experimental results [
44,
45] are also plotted for comparison, which demonstrate the accuracy of our CIS calculations.
5.3 Metal doped nanoclusters
It is well known that the doping process is the most important step in building a semiconductor based microelectronic device. Based on the solid state physics theory, we understand that P-type (N-type) doping into semiconductor bulk materials induces doping energy levels right above (below) the valence (conduction) band edge of the system. While for semiconductor nanoclusters, we found that the CIS method makes it possible to investigate the doping effect from a quantum chemical point of view. Based on the CIS principle, we can “dope” the Fock/Kohn–Sham matrix of a small doped cluster into the big Fock/Kohn–Sham matrix of a very large pure nanocluster, hence achieve an accurate description of the electronic structure of doped semiconductor nanoclusters.
Here we have performed simulations for a 12-shell Si nanoparticle (4 nm in diameter) with one aluminum atom doped in the core. The model system was defined as Al@Si12. With the CIS method, theoretical calculations were carried out at B3LYP level, using 3–21G and 3–21G* basis set to describe the silicon atoms and aluminum atom, respectively. We have obtained the electronic structure, as shown in Fig. 6 (a). One can still observe a big energy gap in the energy structure, which is however not the gap between the HOMO and LUMO of the system. It is found that the LUMO level is only 0.29 eV above the HOMO level, which means electrons in the occupied molecular orbitals can easily be excited to that of LUMO. The wave function of LUMO has been plotted in Fig. 6 (b), showing that LUMO is mainly associated with the doped aluminum atom in the core. Therefore, it is confirmed that doping of an aluminum atom into the Si nanoparticle results in an electronic-unoccupied doping level closed to the upper boundary of occupied molecular orbitals, and will finally induce P-type doping effect in larger and larger systems.
6 Conclusions
We have demonstrated that our recently proposed “central insertion scheme” is indeed an effective method for studying electronic structures of nanomaterials at first principles levels. The accurate results obtained for SWCNT, semiconductor nanoclusters and metal doped nanocluster allow us to not only interpret the corresponding experimental results, but also reveal details that are not accessible to the experiments. It is our hope that the CIS principle can be adopted by other computational programs and be applied to a wide range of nanomaterials.
©Higher Education Press and Springer-Verlag 2009