Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2017; 7: 44645.
Published online 2017 March 16. doi:  10.1038/srep44645
PMCID: PMC5353613

Quantum mechanical calculation of nanomaterial-ligand interaction energies by molecular fractionation with conjugated caps method


Molecular fractionation with conjugate caps (MFCC) method is introduced for the efficient estimation of quantum mechanical (QM) interaction energies between nanomaterial (carbon nanotube, fullerene, and graphene surface) and ligand (charged and neutral). In the calculations, nanomaterials are partitioned into small fragments and conjugated caps that are properly capped, and the interaction energies can be obtained through the summation of QM calculations of the fragments from which the contribution of the conjugated caps is removed. All the calculations were performed by density functional theory (DFT) and dispersion contributions for the attractive interactions were investigated by dispersion corrected DFT method. The predicted interaction energies by MFCC at each computational level are found to give excellent agreement with full system (FS) calculations with the mean energy deviation just a fractional kcal/mol. The accurate determination of nanomaterial-ligand interaction energies by MFCC suggests that it is an effective method for performing QM calculations on nanomaterial-ligand systems.

The growing interest to study the electronic structures of large and complex systems such as protein, DNA, RNA, nanomaterials and other polymers in an efficient manner leads to the development of many computational schemes, which can cope with the limitation of computational resources. Linear scaling methods have been of particular interest because they allow large systems to be treated by quantum mechanics (QM). In particular, the use of quantum fragmentation is a very attractive tool to deal with large macromolecular systems1,2,3,4,5.

In quantum fragment-based approach, a macromolecule is divided into small fragments and simple computations of individual fragments are combined to achieve properties of a whole molecule. The earlier approach of quantum fragmentation for large systems is divide-and-conquer (DAC) method6,7,8, which was introduced by Yang et al. in 1991. In this method, an entire system is divided into subsystems and a buffer region is introduced into each subsystem to account for the interaction between the adjacent subsystems. Later, fragment molecular orbital (FMO) approach9,10,11,12,13 was proposed for ab initio calculations of large biomolecules. In this approach, a biomolecule is divided into a collection of small fragments and MO calculations are performed for each fragment and fragment pairs. One advantage of FMO method is that it can be performed with a number of QM techniques, thus having flexibility in terms of choice of methods for each system. However, FMO method may meet difficulties in the treatment of transition metal bearing systems and elaborate manual efforts are required.

Theoretical chemists have also come up with many other approaches to achieve QM calculations for large molecules. One famous example for calculating large molecules is fragment-based hybrid quantum mechanics/molecular mechanics (QM/MM) approach14,15,16,17 which uses force fields for the larger part of the system and performs more expensive QM calculations for the smaller remaining part. Later, the QM/MM technique within ONIOM framework was developed to carry out large scale modeling in a cost-effective manner18,19,20,21. However, the subsystem that can be treated by QM is limited in most currently available applications.

In recent years, there has been a surge of interests in developing other categories of linear scaling methods to treat complex systems quantum mechanically. Several non-fragmentation-based linear-scaling DFT methods have been developed to overcome the computational scaling of conventional DFT methods and at the same time provide the accuracy of conventional first-principle methods at a fraction of the computational cost for large systems. Examples of widely-used codes include CP2K22, ONETEP23, SIESTA24 and OpenMX25. These approaches have the added advantage of being able to compute forces and perform geometry optimizations and ab initio dynamics. As an alternative, MFCC based approach for the QM calculations of biomaterials was developed by Zhang et al. in 200326. The basic idea of MFCC method is dividing a whole system into fragments and conjugated caps. The properties of the whole molecule can be obtained from the summation of standard QM calculations of the fragments from which the contribution of the conjugated caps is removed. The computational effort in the MFCC approach scales linearly with the molecule size, making it practical to deal with realistic large macromolecular systems. This method has been used for the calculation of interaction energies of protein with water27, drugs28, and ligands29. As well as, this method has been successfully applied to compute energies of charged molecules and large sized molecules such as long oligomers of trans-polyacetylene with up to 400 atoms at the MP2 and HF levels30,31. In a recent study, ab-initio calculation of large protein-ligand systems with 3680, 1798 and 1060 atoms which are beyond the reach of traditional computational methods has been made possible by using the MFCC approach32. Initially, MFCC approach was employed only to study biomolecules, but later on MFCC based attempt to study nanomaterials was reported by Li et al. in 2005 in which BN nanotubes were optimized using both MP2/6-31G* and HF/3-21G levels and the predicted total energy of unit BN and infinite BN nanotubes by MFCC were found to be almost equal to the extrapolated values from conventional MP2 and HF calculations on smaller tubes30.

Studying nanomaterials by QM fragmentation is highly beneficial from the point of view of accurate analyzing and understanding the chemistry of these large systems. Several other QM fragmentation methods have been introduced in the past for the calculation of properties of these large molecules. In 2006, molecular tailoring approach (MTA) was developed by Gadre et al.33,34 to investigate the structures, energetics and reactivity of boric acid nanotubes. Later on, MTA approach for large conjugated systems was reported for one-dimensional and two-dimensional π conjugated molecules by introducing a systematic algorithm for tailoring such systems35. The geometry optimization of these systems within MTA framework was attempted and the generated geometries were found to be in good agreement with their actual counterparts. In 2007, another new fragmentation approach based on DAC strategy was developed by Nakai research group36. They performed DFT calculations of polyene chains CnHn+2 (n = 30–240) using a new DAC method at HF, B3LYP and BLYP levels and found out absolute energy differences between conventional and the new DAC method are comparatively small. In 2009, FMO method was applied to study geometry optimization of a silicon nanowire Si224H162 at several levels of theory37. During the study, the nanowire was divided into 6 fragments of similar size and optimized with a good agreement to experimental results. When reducing the fragment size by dividing the nanowire into 12 fragments, the energy was not well reproduced. However, most of these QM fragmentation methods were applied only for linear molecules, indicating the necessity of the development of an appropriate method which can be employed to study various forms of nanomaterials.

Current methods for computation of nanomaterial ligand interaction energies may be limited and complicated. In a recent study, Car-Parrinello molecular dynamics (CPMD) method was used to calculate interaction energies of isolated and periodic boron nitride nanotubes (BNNTs) with and without water respectively38. When using a force field with the existing Lennard-Jones parameters, the QM interaction energies of water with BNNT can not be reproduced. It was clear that both boron and nitrogen parameters need to be adjusted to obtain a good agreement with the QM energies. Modification and fitting of parameters is complicated and a time consuming process, hence development of a simple and efficient way to calculate nanomaterial-ligand interaction energies is becoming more and more important. Here we present MFCC as an alternative approach for the efficient computation of the interaction energies of tube (carbon nanotube), sphere (fullerene) and surface (graphene) type nanomaterials with both neutral and charged ligands. The main purpose of this study is to prove the applicability and efficiency of MFCC approach to study the intermolecular interactions of different types of nanomaterials by partioning them into smaller fragments. This paper is organized as follows. Section II gives details of MFCC approach and fragmentation procedures of carbon nanotube (CNT), fullerene and graphene surface. In section III, we test the performance of MFCC and make comparison with full system (FS) calculations. In this section, we will test the accuracy and applicability of the MFCC approach within different computational levels such as B3LYP/6-31G* and B97D/6-311G*. The purpose of these calculations is to verify the accuracy of MFCC approach to provide information on different computational levels to account for the nanomaterial-ligand interactions. Finally, a summary is given in section IV.

Computational Methods

Below we briefly outline the overall procedure for constructing subsystems for a given nanomaterial-ligand system. The theoretical background of MFCC approach for the interaction energy calculations of nanomaterials is explained using the CNT-ligand, fullerene-ligand and graphene-ligand systems respectively.

CNT-ligand system

CNT has an infinite long structure comparable to large macromolecules, which gives an ideal example to test the accuracy of MFCC method. During the fragmentation process, we selected a part of long (6, 6) armchair CNT with diameter 8.0 Å and length 9.8 Å without losing the generality. Figure 1 illustrates the MFCC scheme in which the CNT is cut from the middle to produce A/B components, a pair of caps c/c* (colored green) is then inserted to cap the cutoff components to form Ac/c*B fragments, and the caps are fused to form a molecular species c*c (concap). The caps not only perform the task to close the open valency of the components from the cut, they also represent the chemical environment of the part being cut away. Hydrogen atoms are also added to the terminal atoms in the molecular caps to avoid dangling bonds. Now the original CNT system of AB is replaced by three new subsystems, Ac, c*B and the concap c*c.

Figure 1
The MFCC scheme of CNT in which (a) a CNT is cut along the red dashed line, (b) a pair of concap (color green) is inserted to cap the fragments and (c) the concap is fused to form molecular concap species.

By using the MFCC ansatz23, the interaction energy between the CNT and an arbitrary molecule M can be computed by the sum of the interaction energies between the capped fragments and the molecule M from which the contribution of the concap is removed, as expressed by

An external file that holds a picture, illustration, etc.
Object name is srep44645-m1.jpg

where An external file that holds a picture, illustration, etc.
Object name is srep44645-m2.jpg and An external file that holds a picture, illustration, etc.
Object name is srep44645-m3.jpg are the interaction energy of M with Ac and c*B molecular fragments and An external file that holds a picture, illustration, etc.
Object name is srep44645-m4.jpgis the interaction energy of M with the concap c*c.

The above procedure can be easily extended into larger CNT systems, which can be decomposed into multi fragments. The fragmentation procedure of long CNT is shown in Fig. 2. By partitioning the long CNT into N fragments and N-1 concaps, the interaction energy of the CNT system with an arbitrary molecule M can be calculated by the formula

Figure 2
The MFCC scheme of long CNT in which a long CNT is cut into N fragments and N-1 concap species.
An external file that holds a picture, illustration, etc.
Object name is srep44645-m5.jpg

where An external file that holds a picture, illustration, etc.
Object name is srep44645-m6.jpg and An external file that holds a picture, illustration, etc.
Object name is srep44645-m7.jpg denote the interaction energy between the ith fragment and M, and the ith concap and M respectively. Using Eq. (2), the interaction energy can be obtained by simple summation over individual interaction energies between the molecule M and the capped CNT fragments that can be determined by ab-initio calculations such as HF, DFT or even higher level QM method. Obviously, the above method scales linearly with the size of the CNT and therefore applicable to large CNT interaction systems.

Then we extend our MFCC method to study hybrid nanomaterials. Encapsulated C60@CNT is a novel hybrid carbon material, which traps fullerene inside the CNT with interesting molecular properties. In order to study encapsulated C60@CNT system, we applied decomposition procedure as shown in Fig. 3. The CNT is decomposed as shown by red vertical dotted lines, and it would produce Ac, c*B fragments and c*c concap. Using the MFCC calculations on these formed fragments, the interaction energy can be determined by the summation of the contributions of all capped fragments with the contributions of all the conjugated caps being removed, which is given by Eq. (1).

Figure 3
The MFCC scheme of CNT with a fullerene inside.

Fullerene-ligand system

Next, we design series of calculations to investigate the interactions of sphere type nanomaterial using C60 fullerene. C60 is similar in chemical structure to CNT, except they have closed ends to form a hollow spherical structure, which is called as buckyball. Unlike the cylindrical CNT, C60 has a limited and confined chemical environment inside its wall. They are able to hold host atoms such as water, metal ions or single molecules in their interior cages to form endohedral fullerenes. Many research groups have extensively studied the host-guest interactions of endohedral fullerenes in the past39,40,41,42,43,44,45,46,47,48,49,50. In this study, an arbitrary molecule such as water or metal cation is placed inside C60 fullerene and MFCC approach is applied to study fullerene-ligand interactions. The decomposition procedure of C60 and the formed Ac, c*B fragments and the concap c*c are shown in Fig. 4. Although C60 is a relatively small system for fragmentation, this procedure can be easily extended into large fullerenes such as C70, C84, C240 or even larger fullerenes. The interaction energy between C60 and the molecule M inside can be computed by separate calculations of individual fragments interacting with M, as given by the Eq. (1).

Figure 4
The MFCC scheme of C60 fullerene in which (a) a fullerene splits in two along the red dashed line, (b) a pair of concap (color green) is inserted to cap the hemisphere fragments and (c) the concap is fused to form molecular concap species.

Graphene-ligand system

Finally, we extend our MFCC study to investigate interactions of surface type carbon nanomaterials. Graphene is a planar form of sp2-bonded carbon material, which emerged as a promising nanoelectronic device along with CNTs and fullerenes. However, QM calculations of graphene sheets with large number of atoms are difficult to perform owing to the huge computational cost. MFCC method can be used as an alternative to carry out the QM studies on graphene-ligand interactions by decomposing graphene into smaller fragments. In present study, three different fragmentation methods are applied to arrange graphene fragments as shown in Fig. 5. (a) with 4 fragments, (b) with 6 fragments and 1 layer conjugated cap, and (c) with 6 fragments and 3 layer conjugated cap.

Figure 5
The MFCC scheme of graphene in which (a) a graphene is partitioned along the red dashed line to form (a) 4 fragments with 1-layer concap, (b) 6 fragments with 1-layer concap, and (c) 6 fragments with 3-layer concap.

First, the graphene surface is decomposed into 4 fragments (A + B + C + D) as shown in Fig. 5a in which the interaction energy can be given by following formula.

An external file that holds a picture, illustration, etc.
Object name is srep44645-m8.jpg

where An external file that holds a picture, illustration, etc.
Object name is srep44645-m9.jpg, An external file that holds a picture, illustration, etc.
Object name is srep44645-m10.jpg,An external file that holds a picture, illustration, etc.
Object name is srep44645-m11.jpg, An external file that holds a picture, illustration, etc.
Object name is srep44645-m12.jpg are the interaction energies of M with A, B, C, D molecular fragments and An external file that holds a picture, illustration, etc.
Object name is srep44645-m13.jpg, An external file that holds a picture, illustration, etc.
Object name is srep44645-m14.jpg, An external file that holds a picture, illustration, etc.
Object name is srep44645-m15.jpg, An external file that holds a picture, illustration, etc.
Object name is srep44645-m16.jpg, are the interaction energies of M with the ab, ac, cd, bd conjugated caps and An external file that holds a picture, illustration, etc.
Object name is srep44645-m17.jpg is the interaction energy of overlapped area of surface e (colored yellow).

Subsequently, the graphene surface is decomposed into 6 fragments (A + B + C + D + E + F) as shown in Fig. 5b, which the interaction energy can be given by following formula.

An external file that holds a picture, illustration, etc.
Object name is srep44645-m18.jpg

where ΔEAM, ΔEBM, ΔECM, ΔEDM, ΔEEM, ΔEFM are the interaction energies of M with A, B, C, D, E, F fragments and ΔEabM, ΔEbcM, ΔEadM, ΔEdeM, ΔEbeM, ΔEefM, ΔEcfM are interaction energies of M with ab, bc, ad, de, be, ef, cf conjugated caps and ΔEgM, ΔEhM are the interaction energies of M with the overlapped area of surface g and h (colored yellow).

In order to introduce different type of conjugated caps, we design series of calculations by increasing conjugated cap size from 1-layer benzene cap to 3-layer benzene cap as shown in Fig. 5c. The interaction energy between graphene and molecule M can be computed by separate calculations of individual fragments interacting with M, as same as given by the Eq. (4).

After the cut, hydrogen atoms are added in each fragment and concap using GaussView program. It should be pointed out that the geometries of the cap atoms in fragments are kept exactly the same with the geometries of the atoms in corresponding concap to ensure that the artifitial interactions between the molecule M and the concaps are cancelled. Both FS and MFCC single point energy calculations were performed with SCF convergence criterion of 10−8 by using the method contained ground state, DFT, default spin, B3LYP hybrid functional, with basis set 6-31G* within GAUSSIAN 09 program package51. During the calculation, ‘No symmetry’ keyword is used to prevent molecule reorientation so that all computations are performed in the input orientation. Since B3LYP functional lacks of the long range dispersion energy term, CNT-ligand calculations were recalculated by dispersion corrected DFT method of Grimme at B97D/6-311G* level52,53. MFCC interaction energies of each system were obtained by above given formulas and traditional single point FS computations at the same level were used as references to assess the accuracy of MFCC approach.

Results and Discussion

MFCC study on CNT

As the first task in our study, MFCC method was applied to study CNT-water system. The penetration of water into the hydrophobic pore of CNT has been attracted considerable research attentions over the years, ultimately leading to the development of promising technologies for drug delivery, selective ion transport, and nano-filtration methods54,55,56. The behavior of water inside CNT can provide important insight into the properties of nano-confined interfacial water on hydrophobic surface. In this study, we used MFCC approach to investigate the interaction between water and (6, 6) CNT channel wall, one of the smallest diameter CNTs known to be water permeable. The water molecule was placed inside the CNT and moved inside along a random pathway as shown in Fig. 6. In the beginning, we calculated MFCC and FS interaction energies by B3LYP/6-31G* computational level. Since the dispersion energy is an important part of noncovalent interactions for CNT-water system, calculations were repeated by dispersion corrected DFT using B97D/6-311G*. As shown in Fig. 6, the MFCC and FS calculated interaction energy curves by both computational levels show good agreement with each other. The mean variation of MFCC and FS energies by B3LYP/6-31G* was calculated as 1.38 kcal/mol while the mean variation by B97D/6-311G* as 0.08 kcal/mol, showing that MFCC approach gives a good performance on both computational levels. As listed in Table 1, the average interaction energy calculated by dispersion corrected DFT is −8.0 kcal/mol, which is closer to the previously calculated energy of water molecule adsorbed in (6, 6) CNT (−8.89 kcal/mol) using two phase thermodynamics (2PT) method57. Evidently, B3LYP functional underestimates the CNT-water interaction, while interaction energies can be improved significantly with the dispersion corrected DFT due to the inclusion of van der Waals interaction between CNT wall and water.

Figure 6
A molecule (water or Na+) travels inside CNT along a random pathway to generate the interaction energy curves obtained by MFCC and FS calculations using B3LYP/6-31G* and B97D/6-311G* level respectively.
Table 1
Interaction energies of CNT-water system (Energies are in kcal/mol).

The following task in our study is to examine the reliability of MFCC method in the calculation of interaction energies of CNT with charged ligand. Single charged ions marching through the CNT channel have been observed by Lee et al. in 201058. They found that charged molecules such as sodium or chloride ions can flow rapidly through CNT and such channels can be useful as sensitive detectors or water desalination systems. In our study, MFCC approach was used to investigate the interaction of Na+ with CNT channel by moving Na+ cation inside the CNT channel. Table 2 and Fig. 6 reports the interaction energies and energy curves calculated by B3LYP/6-31G* and B97D/6-311G* computational levels. The mean deviation by B3LYP/6-31G* was calculated as 0.96 kcal/mol, while dispersion corrected DFT gave the mean deviation of 0.35 kcal/mol. Although CNT-water system shows larger interaction energy difference between two computational levels, CNT-Na+ system does not show such a trend, showing that motion of cationic Na+ inside CNT is obviously dominated by electrostatic interactions between CNT wall and charged metal cation, rather than van der Waals interactions.

Table 2
Interaction energies of CNT-Na+ system (Energies are in kcal/mol).

For further extension of our study, MFCC calculations were performed on hybrid nanomaterial system of C60@CNT. C60 trapped inside CNT has been experimentally detected by Smith et al.59,60, showing that stable and closed carbon shells exist inside single-walled CNTs. Measurements of the diameter of these carbon shells suggested that many of them are C60 fullerenes. After that, various research groups have studied the C60@CNT systems61,62,63,64. The energetics of encapsulation of C60 inside CNT and electronic structures has been studied by Okada et al.65 and found out encapsulating process is exothermic for (10, 10) CNT, whereas the processes are endothermic for both (8, 8) and (9, 9) nanotubes. In our study, we selected a part of (15, 15) armchair CNT with radius of 10.1 Å and length of 9.8 Å with C60 molecule inside CNT and moving along a random pathway shown in Fig. 7 to generate interaction energy curves. The MFCC fragmentation procedure for CNT is exactly the same as previously mentioned as shown in Fig. 3. At first, calculations were performed at B3LYP/6-31G* computational level and then calculations were repeated at B97D/6-311G* level to include the dispersion contribution. As shown in Table 3, the B3LYP/6-31G* calculation shows the CNT-C60 interactions are almost zero. Although a change of sign between the energies calculated via FS and MFCC is observed at several positions, absolute energy value differences are still in good control to be less than 0.05 kcal/mol. When employing dispersion corrected DFT, interaction energy minimum for all the considered geometries is clearly visible. The energetics of encapsulation of C60 inside CNT has been previously studied by A. Rochefort using MM3 molecular mechanics force field and found out the reaction is exothermic with interaction energy of 12.00 kcal/mol66, which is similar to B97D/6-311G* calculated interaction energy of our most stable C60@CNT structure (−12.10 kcal/mol). Obviously, the B3LYP/6-31 G* calculated energies are once again underestimated significantly, while dispersion corrected DFT gives much reliable energies. The calculated interaction energies by MFCC and FS calculations using two computational levels show less than 0.01 kcal/mol deviations, showing that MFCC method is a good choice to perform QM studies on tube type nanomaterials.

Figure 7
A C60 molecule travels inside CNT along a random pathway to generate the interaction energy curves obtained by MFCC and FS calculations using B3LYP/6-31G* and B97D/6-311G* level respectively.
Table 3
Interaction energies of C60@CNT system (Energies are in kcal/mol).

MFCC study on C60 Fullerene

Since our major aim is to explore the applicability of MFCC approach for different types of carbon nanomaterials, next we extend our study to investigate interactions of sphere type nanomaterials. The finding of endohedral fullerenes which encapsulate single molecule or metal ion inside fullerene cage is one of the most exciting finding in nanomaterial chemistry67,68. Doping various elements into C60 fullerene gave rise to the materials with novel properties and developed an entirely new branch of chemistry. Recently, Kurotabi et al. found that single water molecule can be confined in fullerene and intrinsic properties of water molecule, as well as physical properties of fullerene cage can be controlled by this method69. Here we applied MFCC approach to study the interactions of water, Li+ and K+ with C60 fullerene. The ligands were kept at the center of C60 and moved inside along one of the coordinate axes. The interaction energy variation as a function of the distance for water, Li+, K+ ligands by B3LYP/6-31G* computational level is shown in Fig. 8, while the calculated energies are summarized in Table 4. The zero position in the figure refers to the central position of the buckyball. The curves in Fig. 8 represent the one-dimensional potential energy profile when the small molecule is approaching the central position of a 6-member ring on the wall of fullerene from the central position of fullerene. As can be seen, both water and K+ are most stable in the center of C60 and interaction energy gradually increases with closer to the wall of C60, while Li+ is most stable off the center, approximately 1.2 Å away from the center of C60. It is obvious that the relative interaction energy curves computed by MFCC are in good agreement with the standard FS results with a mean deviation of 0.21, 1.00, 1.27 kcal/mol for C60-water, C60-Li+, and C60-K+ respectively. Large energy differences (above 2 kcal/mol) often come from when small molecule is very close to the inside wall of fullerene, leading to high repulsion energy. The MFCC calculation can be improved by introducing a larger cap to incorporate many-body effects more or less. The reasonable consistency of MFCC calculated energies with FS calculated values proves that MFCC method can be successfully applied to study QM interactions of sphere type nanomaterial-ligand systems.

Figure 8
From the central position of C60, a molecule (water, Li+ or K+) moves towards the center of a 6-member ring on the C60 to generate the interaction energy curves obtained by MFCC and FS calculations using B3LYP/6-31G* and B97D/6-311G* level respectively. ...
Table 4
Interaction energies of C60-water/Li+/K+ systems (Energies are in kcal/mol).

MFCC study on Graphene

As an effort to extend the MFCC study into larger nanomaterial systems, calculations were performed to calculate graphene-ligand interaction energies. Graphene sheets can be considered as an unrolled CNT with a monolayer of carbon atoms and similar stable physical properties to CNT. Due to its extremely large surface area and high porosity, graphene is ideal for adsorption of gases such as H2, CH4 and CO. The interaction between graphene and CO has been studied extensively due to its application in graphene-based gas sensors. In this study, we constructed a large graphene surface with CO as the ligand and MFCC calculations were performed by decomposing the surface into smaller fragments. Since the most stable adsorption site for CO on pristine graphene has been found to be at hollow site70, we started our study by keeping CO 1.00 Å away from the hollow site of one benzene ring. In order to generate interaction energy curves, CO molecule was gradually moved away from the surface. All the calculated energies are summarized in Table 5, which shows the CO-graphene interaction is highly repulsive when CO is closer to the surface and gradually stabilizes when moving away from the surface. As shown in the Fig. 9a, the calculated graphene-CO interaction energies by both FS and MFCC methods show a good agreement with only 0.13 kcal/mol mean deviation for energy values below 10 kcal/mol, confirming the reliability of MFCC approach.

Figure 9
A CO molecule moves on graphene surface as shown to generate the interaction energy curves obtained by MFCC and FS calculations using B3LYP/6-31G* level.
Table 5
Interaction energies of graphene-CO/4 fragment system (Energies are in kcal/mol).

Next, we continued the MFCC study by decomposing the large graphene sheet into smaller fragments to check the effect of fragment size on calculated interaction energies. To this aim, we decomposed the graphene surface into 6 fragments as shown in Fig. 5b and extended our study by moving CO ligand vertically with slight oscillation on the graphene surface. The interaction energies calculated by MFCC and FS are plotted in Fig. 9b and listed in Table 6. As can be seen, interaction energy values depend on the graphene-CO distance and ranges from −2.8 to −3.8 kcal/mol. It is obvious that our MFCC method well reproduced the FS interaction energies with mean deviation only a fractional kcal/mol, which proves the reliability of MFCC approach for studying large graphene surfaces using smaller fragments. Although the MFCC fragmentation of graphene surface into a number of smaller fragments is a little complicated, it decreases the computational cost of each fragment and allows faster evaluation of MFCC interaction energies.

Table 6
Interaction energies of graphene-CO/6 fragment system (Energies are in kcal/mol).

In order to verify the effect of the size of conjugated caps on the calculated interaction energies, additional series of calculations were designed by increasing the conjugated cap size into 3 benzene layers as shown in Fig. 5c. The graphene surface was decomposed into 6 fragments and CO ligand moved diagonally on the surface with a higher oscillation movement. As can be seen in Table 7 and Fig. 9c, the interaction energies are fluctuating between highly positive and negative values, depending on the graphene-CO distance. The largest deviation of graphene-CO interaction energies between FS and MFCC calculations is less than 0.5 kcal/mol, which suggests that accuracy of MFCC is not affected by the size of conjugated caps. Since larger cap makes MFCC computation more expensive because their capping enlarges the size of the fragments and increases the computational cost, employing a smaller cap as possible is enough to get accurate MFCC interaction energies. Obviously, all the three MFCC fragmentation procedures used on graphene surface give satisfactory results, indicating our fragmentation and capping procedure provides a new alternative approach for QM interaction energy calculations of surface type nanomaterials.

Table 7
Interaction energies of graphene-CO/6 fragment system with 3-benzene layer conjugated cap (Energies are in kcal/mol).


We have performed a series of MFCC benchmark calculations to study nanomaterial-ligand interactions by decomposing three types of nanomaterials, namely carbon nanotube, fullerene, and graphene into small fragments respectively. Unlike the application of MFCC approach on large biomolecules which is more complicated and difficult to assess the percentage error with FS results, its application on carbon nanomaterials can be an ideal benchmark system to evaluate the accuracy of MFCC approach due to simplicity and orderness of nanomaterial structure. Obviously, our benchmark tests demonstrate that MFCC on various forms of nanomaterials gives consistently accurate interaction energies with very small deviations compared to the corresponding FS calculated energies. Even so, the MFCC scheme still brings some errors. Considering the interaction energy curves (Figs 6, ,7,7, ,8,8, ,9),9), we can see roughly that the MFCC scheme gives smaller or more negative interaction energies compared with the FS calculations. The part of the missing values is probably due to not fully considering the many-body effects. When the FS calculation produces a very small positive energy value, the MFCC scheme is likely to give a negative energy value, resulting in a change of sign between the energies calculated via FS and MFCC.

To illustrate the computational expense of our MFCC approach, information about average calculation time of FS and MFCC calculations is given in Table 8. As an example, the computational time for a single point FS calculation on graphene surface is about two days, while less than 6 hours for MFCC fragment calculations. And the advantage on computing time becomes more prominent with higher correlation methods and for larger nanomaterial systems due to the linear scaling feature of the MFCC method versus standard N3 for HF/DFT method and N5 for MP2 method. An additional attractive feature of the MFCC method is that its ab initio calculation can be easily parallelized to run on multi-node computer clusters that could dramatically speed up the computation.

Table 8
Average computational time of FS and MFCC fragments (B3LYP/6-31G*).

The excellent agreement of MFCC calculated energies with FS values indicates that MFCC method is a promising theoretical tool for fast and reliable calculation of ab-initio interaction energies of nanomaterials with both neutral and charged ligands. This method is almost linear scaling with system size and thus makes it possible to obtain interaction energies and molecular properties even for large nanomaterial-ligand systems.

Additional Information

How to cite this article: Zhang, D. Quantum mechanical calculation of nanomaterial-ligand interaction energies by molecular fractionation with conjugated caps method. Sci. Rep. 7, 44645; doi: 10.1038/srep44645 (2017).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


D.W. Z. is supported by Henan University of Science and Technology start-up grant. D.W. Z. also thanks Prof. Ye Mei at East China Normal University to provide the computational programs.


The authors declare no competing financial interests.

Author Contributions D.Z. designed, performed the calculations and wrote the manuscript.


  • Gordon M. S., Fedorov D. G., Pruitt S. R. & Slipchenko L. V. Fragmentation Methods: A Route to Accurate Calculations on Large Systems. Chem. Rev. 112, 632–672 (2012). [PubMed]
  • Yam C. Y., Zhang Q., Wang F. & Chen G. H. Linear-scaling Quantum Mechanical methods for Exited States. Chem. Soc. Rev. 41, 3821–3838 (2012). [PubMed]
  • Li S. & Li W. Fragment Energy Approach to Hartree-Fock Calculations of Macromolecules. Annual Reports Section 104, 256–271 (2008).
  • Jacobson L. D. & Herbert J. M. An Efficient Fragment-based Electronic Structure Method for Molecular Systems: Self-consistent Polarization with Perturbative Two-body Exchange and Dispersion. J. Chem. Phys. 134, 094118 (2011). [PubMed]
  • Gordon M. S., Mullin J. M., Pruitt S. R., Roskop L. B., Slipchenko L. V. et al. . Accurate Methods for Large Molecular Systems. J. Phys. Chem. B 113, 9646–9663 (2009). [PubMed]
  • Yang W. T. Direct Calculation of Electron Density in Density Functional Theory. Phys. Rev. Lett. 66, 1438–1441 (1991). [PubMed]
  • He X. & Merz K. M. Divide and Conquer Hartree-Fock Calculations on Proteins. J. Chem. Theory Comput. 6, 405–411 (2010). [PMC free article] [PubMed]
  • Yang W. T. & Lee T. S. A Density Matrix Divide and Conquer Approach for Electronic Structure Calculations of Large Molecules. J. Chem. Phys. 103, 5674–5681 (1995).
  • Kitaura K., Ikeo E., Asada T., Nakano T. M. & Uebayasi M. Fragment Molecular Orbital Method: An Approximate Computational Method for Large Molecules. Chem. Phys. Lett. 313, 701–706 (1999).
  • Nakano T., Kaminuma T., Sato T., Akiyama Y., Uebayasi M. et al. . Fragment Molecular Orbital Method: Application to Polypeptides, Chem. Phys. Lett. 318, 614–618 (2000).
  • Nakano T., Kaminuma T., Sato T., Fukuzawa K., Akiyama Y. et al. . Fragment Molecular Orbital Method: Use of Approximate Electrostatic Potential. Chem. Phys. Lett. 351, 475–480 (2002).
  • Kitaura K. & Fedorov D. The Fragment Molecular Orbital Method Practical Applications to Large Molecular Systems. CRC, (Boca Rotan, FL) Chap. 2, 5–36 (2009).
  • Fedorov D. G., Nagata T. & Kitaura K. Exploring Chemistry with the Fragment Molecular Orbital Method. Phys. Chem. Chem. Phys. 14, 7562–7577 (2012). [PubMed]
  • Gao J. Towards a Molecular Orbital Derived Empirical Potential for Liquid Simulations. J. Phys. Chem. B 101, 657–663 (1997).
  • Xie W. & Gao J. Design of a Next Generation Force Field: The X-POL Potential. J. Chem. Theory Comput. 3, 1890–1900 (2007). [PMC free article] [PubMed]
  • Trajbl M., Hong G. Y. & Warshel A. Ab Initio QM/MM Simulation with Proper Sampling: First Principle Calculations of the Free Energy of the Autodissociation of Water in Aqueous Solution. J. Phys. Chem. B 106, 13333–13343 (2002).
  • Shurki A. & Warshel A. Structure/Function Correlations of Proteins using MM, QM/MM and Related Approaches: Methods, Concepts, Pitfalls and Current Progress. Adv. Protein Chem. 66, 249–313 (2003). [PubMed]
  • Murphy R. B., Philipp D. M. & Friesner R. A. A Mixed Quantum Mechanics/Molecular Mechanics (QM/MM) Method for Large-scale Modeling of Chemistry in Protein Environments. J. Comput. Chem. 21, 1442–1457 (2000).
  • Vreven T., Morokuma K., Farkas O., Schlegel H. B. & Frisch M. J. Geometry Optimization with QM/MM, ONIOM and other Combined Methods. I. Microiterations and Constraints. J. Comput. Chem. 24, 760–769 (2003). [PubMed]
  • Canfield P., Dahlbom M. G., Hush N. S. & Reimers J. R. Density Functional Geometry Optimization of the 150000 atom Photosystem-I Trimer. J. Chem. Phys. 124, 024301 (2006). [PubMed]
  • Hratchian H. P., Parandekar P. V., Raghavachari K., Frisch M. J. & Vreven T. QM:QM Electronic Embedding using Mulliken Atomic Charges: Energies and Analytic Gradients in an ONIOM Framework. J. Chem. Phys. 128, 034107 (2008). [PubMed]
  • Hutter J., Iannuzzi M., Schiffmann F. & VandeVondele J. cp2k: atomistic simulations of condensed matter systems. WIREs Comput. Mol. Sci. 4, 15–25 (2014).
  • Haynes P. D., Skylaris C., Mostofi A. A. & Payne M. C. ONETEP: linear-scaling density-functional theory with local orbitals and plane waves. Physica Status Solidi B. 243, 2489–2499 (2006).
  • Soler J. M., Artacho E., Gale J. D., García A., Junquera J. et al. . The SIESTA method for ab initio order-N materials simulation. J. Phys. Condens. Matter 14, 2745–2779 (2002).
  • Boker S., Neale M., Maes H., Wilde M., Spiegel M. et al. . OpenMx: An Open Source Extended Structural Equation Modeling Framework. Psychometrika. 76(2), 306–317 (2011). [PMC free article] [PubMed]
  • Zhang D. W. & Zhang J. Z. H. Molecular Fractionation with Conjugate Caps for Full Quantum Mechanical Calculation of Protein-molecule Interaction Energy. J. Chem. Phys. 119, 3599–3605 (2003).
  • Zhang D. W., Chen X. H. & Zhang J. Z. H. Molecular Caps for Full Quantum Mechanical Computation of Peptide-Water Interaction Energy. J. Comput. Chem. 24, 1846–1852 (2003). [PubMed]
  • Zhang D. W. & Zhang J. Z. H. Full Quantum Mechanical Study of Binding of HIV-1 Protease Drugs. Int. J. Quant. Chem. 103, 246–257 (2005).
  • Zhang D. W., Xiang Y., Gao A. M. & Zhang J. Z. H. Quantum Mechanical Map for Protein-ligand Binding with Application to β-Trypsin/benzamidine Complex. J. Chem. Phys. 120, 1145–1148 (2004). [PubMed]
  • Li S., Li W. & Fang T. An Efficient Fragment-Based Approach for Predicting the Ground-State Energies and Structures of Large Molecules. J. Am. Chem. Soc. 127, 7215–7226 (2005). [PubMed]
  • Li W., Li S. & Jiang Y. Generalized Energy-Based Fragmentation Approach for Computing the Ground-State Energies and Properties of Large Molecules. J. Phys. Chem. A 111, 2193–2199 (2007). [PubMed]
  • Antony J. & Grimme S. Fully Ab Initio Protein-Ligand Interaction Energies with Dispersion Corrected Density Functional Theory. J. Comput. Chem. 33, 1730–1739 (2012). [PubMed]
  • Gadre S. R. & Ganesh V. Molecular Tailoring Approach: Towards PC-based Ab Initio Treatment of Large Molecules. J. Theor. Comput. Chem. 5, 835–856 (2006).
  • Elango M., Subramanian V., Rahalkar A. P., Gadre S. R. & Sathyamurthy N. Structure, Energetics, and reactivity of Boric Acid Nanotubes: A Molecular Tailoring approach. J. Phys. Chem. A 112, 7699–7704 (2008). [PubMed]
  • Yeole S. D. & Gadre S. R. On the Applicability of Fragmentation Methods to Conjugated π systems within Density Functional Framework. J. Chem. Phys. 132, 094102 (2010). [PubMed]
  • Akama T., Kobayashi M. & Nakai H. Implementation of Divide-and-Conquer Method including Hartree-Fock Exchange Interaction. J. Comput. Chem. 28, 2003–2012 (2007). [PubMed]
  • Fedorov D. G., Avramov P. V., Jensen J. H. & Kitaura K. Analytic Gradient for the Adaptive Frozen Orbital Bond Detachment in the Fragment Molecular Orbital Method. Chem. Phys. Lett. 477, 169–175 (2009).
  • Hilder T. A., Yang R., Ganesh V., Gordon D., Bliznyuk A. et al. . Validity of Current Force Fields for Simulations of Boron Nitride Nanotubes. Micro Nano Lett. 5, 150–156 (2010).
  • Maheshwari S., Chakraborty D. & Sathyamurthy N. Possibility of proton motion through buckminsterfullerene. Chem. Phys. Lett. 315, 181–186 (1999).
  • Patchkovskii S. & Thiel W. Equilibrium yield for helium incorporation into buckminsterfullerene: Quantum-chemical evaluation. J. Chem. Phys. 106, 1796–1799 (1997).
  • Buckingham A. D. & Read J. P. Degeneracy loss contributions to the stabilisation of the eccentric position of lithium in Li@C60. Chem. Phys. Lett. 253, 414–419 (1996).
  • Cioslowski J. & Fleischmann E. D. Endohedral complexes: Atoms and ions inside the C60 cage. J. Chem. Phys. 94, 3730–3734 (1991).
  • Cioslowski J. & Nanayakkara A. Endohedral effect in inclusion complexes of the C60 cluster. J. Chem. Phys. 96, 8354–8362 (1992).
  • Williams C. I., Whitehead M. A. & Pang L. Interaction and dynamics of endohedral gas molecules in fullerene C60 isomers and C70. J. Phys. Chem. 97, 11652–11656 (1993).
  • Dunlap B. I., Ballester J. L. & Schmidt P. P. Interactions between fullerene (C60) and endohedral alkali atoms. J. Phys. Chem. 96, 9781–9787 (1992).
  • Park J. M., Tarakeshwar P., Kim K. S. & Clark T. Nature of the interaction of paramagnetic atoms (A = 4N, 4P, 3O, 3S) with π systems and C60: A theoretical investigation of A•••C6H6 and endohedral fullerenes A@C60. J. Chem. Phys. 116, 10684–10691 (2002).
  • Pyykkö P., Wang C., Straka M. & Vaara J. A London-type formula for the dispersion interactions of endohedral A@B systems. Phys. Chem. Chem. Phys. 9, 2954–2958 (2007). [PubMed]
  • Cross R. J. Does H2 Rotate Freely Inside Fullerenes? J. Phys. Chem. A 105, 6943–6944 (2001).
  • Dodziuk H. Modeling complexes of H2 molecules in fullerenes. Chem. Phys. Lett. 410, 39–41 (2005).
  • Shameema O., Ramachandran C. N. & Sathyamurthy N. Blue Shift in X-H Stretching Frequency of Molecules Due to Confinement. J. Phys. Chem. A. 110, 2–4 (2006). [PubMed]
  • Frisch M. J., Trucks G. W., Schlegel H. B., Scuseria G. E., Robb M. A. et al. . Gaussian 09, Revision A.02, Gaussian, Inc., Wallingford, CT (2009).
  • Grimme S. Accurate Description of Van der Waals Complexes by Density Functional Theory including Empirical Corrections. J. Comput. Chem. 25, 1463–1473 (2004). [PubMed]
  • Grimme S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 27, 1787–1799 (2006). [PubMed]
  • Zahab A., Spina L., Poncharal P. & Marlière C. Water-vapor Effect on the Electrical Conductivity of a Single-Walled Carbon Nanotube Mat. Phys. Rev. B 62, 10000–10003 (2000).
  • Hummer G., Rasaiah J. C. & Noworyta J. P. Water Conduction through the Hydrophobic Channel of a Carbon Nanotube. Nature 414, 188–190 (2001). [PubMed]
  • Sansom S. M. P. & Biggin P. C. Water at the Nanoscale. Nature 414, 156–159 (2001). [PubMed]
  • Kumar H., Mukherjee B., Lin S. T., Dasgupta C., Sood A. K. et al. . Thermodynamics of Water Entry in Hydrophobic Channels of Carbon Nanotubes. J. Chem. Phys. 134, 124105 (2011). [PubMed]
  • Lee C. Y., Choi W., Han J. H. & Strano M. S. Coherence Resonance in a Single-Walled Carbon Nanotube Ion Channel. Science 329, 1320–1324 (2010). [PubMed]
  • Smith B. W., Monthioux M. & Luzzi D. E. Encapsulated C60 in Carbon Nanotubes. Nature 396, 323–324 (1998).
  • Smith B. W. & Luzzi D. E. Formation Mechanism of Fullerene Peapods and Coaxial Tubes: A Path to large Scale Synthesis. Chem. Phys. Lett. 321, 169–174 (2000).
  • Wang Q., Kitaura R., Yamamoto Y., Arai S. & Shinohara H. Synthesis and TEM structural characterization of C60-flattened carbon nanotube nanopeapods. Nano Res. 7, 1843–1848 (2014).
  • Chadli H., Rahmani A. & Sauvajol J.-L. Raman spectra of C60 dimer and C60 polymer confined inside a (10, 10) single-walled carbon nanotube. Journal of Physics Condensed Matter 22(14), 145303 (2010). [PubMed]
  • Hodak M. & Girifalco L. A. Systems of C60 molecules inside (10, 10) and (15, 15) nanotube: A Monte Carlo study. Physical Review B 68(8), 85405 (2003).
  • Yu H. Y., Lee D. S., Lee S. H., Kim S. S., Lee S. W. et al. . Single-electron transistor mediated by C60 insertion inside a carbon nanotube. Applied Physics Letters 87(16), 163118–163120 (2005).
  • Okada S., Saito S. & Oshiyama A. Energetics and Electronic Structures of Encapsulated C60 in a Carbon Nanotube. Phys. Rev. Lett. 86, 3835–3838 (2001). [PubMed]
  • Rochefort A. Electronic and Transport Properties of Carbon Nanotube Peapods. Phys. Rev. B 67, 115401 (2003).
  • Heath J. R., O´Brien S. C., Zhang Q., Liu Y., Curl R. F. et al. . Lanthanum Complexes of Spheroidal Carbon Shells. J. Am. Chem. Soc. 107, 7779–7780 (1985).
  • Chai Y., Guo T., Jin C., Haufler R. E., Chibante L. P. F. et al. . Fullerenes with Metals Inside. J. Phys. Chem. 95, 7564–7568 (1991).
  • Kurotabi K. & Murata Y. A Single Molecule of Water Encapsulated in Fullerene C60. Science 333, 613–616 (2011). [PubMed]
  • Schedin F., Geim A. K., Morozov S. V., Hill E. W., Blake P. et al. . Detection of Individual Gas Molecules Adsorbed on Graphene. Nat. Mater. 6, 652–655 (2007). [PubMed]

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group