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.

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   26 . 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 water 27 , drugs 28 , and ligands 29 . 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 levels 30,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 approach 32 . 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 tubes 30 .
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 systems 35 . 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 group 36 . They performed DFT calculations of polyene chains C n H n+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 Si 224 H 162 at several levels of theory 37 . 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 respectively 38 . 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.
By using the MFCC ansatz 23 , 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 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). 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 past [39][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). 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  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.
are the interaction energies of M with the ab, ac, cd, bd conjugated caps and ∆ − E e M 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.  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 package 51 . 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* level [52][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 methods [54][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) method 57 . 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.
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 2010 58 . 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.
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 systems [61][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/mol 66 , 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.
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 Scientific RepoRts | 7:44645 | DOI: 10.1038/srep44645 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 chemistry 67,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 method 69 . 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 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 2. Interaction energies of CNT-Na + system (Energies are in kcal/mol). The maximum absolute difference between the two approaches is represented with bold text. † Absolute difference = |MFCC − FS|. ‡ Percentage absolute difference (%) = |Absolute difference/FS|. extremely large surface area and high porosity, graphene is ideal for adsorption of gases such as H 2 , CH 4 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  Table 3. Interaction energies of C60@CNT system (Energies are in kcal/mol). The maximum absolute difference between the two approaches is represented with bold text. The percentage absolute differences for the calculations based on B3LYP/6-31G* are neglected since a change of sign between the energies calculated via FS and MFCC is observed at several positions, leading to abnormally large percentage error. † Absolute difference = |MFCC − FS|. ‡ Percentage absolute difference (%) = |Absolute difference/FS|. 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 site 70 , 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. 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.
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.

Conclusion
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    Table 6. Interaction energies of graphene-CO/6 fragment system (Energies are in kcal/mol). The maximum absolute difference between the two approaches is represented with bold text. † Absolute difference = |MFCC − FS|. ‡ Percentage absolute difference (%) = |Absolute difference/FS|. 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-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 N 3 for HF/DFT method and N 5 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.
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.  Table 7. Interaction energies of graphene-CO/6 fragment system with 3-benzene layer conjugated cap (Energies are in kcal/mol). † Absolute difference = |MFCC − FS|. ‡ Percentage absolute difference (%) = |Absolute difference/FS|. The maximum absolute difference between the two approaches is represented with bold text.