Comment pubs.acs.org/JPCC
Reply to “Comment on ‘How the Number and Location of Lithium Atoms Affect the First Hyperpolarizability of Graphene’” Yang-Yang Hu, Hong-Liang Xu,* and Zhong-Min Su* Institute of Functional Material Chemistry, Faculty of Chemistry, Northeast Normal University, Changchun 130024, Jilin, People’s Republic of China
J. Phys. Chem. C 2010, 114 (46), 19792−19798. DOI: 10.1021/jp105045j J. Phys. Chem. C 2013, 117. DOI: 10.1021/jp3057256 S Supporting Information *
revision B.03.4 The corresponding input file is presented in the Supporting Information (SI). This file was calculated again by Gaussian 09 revision B.01,5 and the unexpectedly large first hyperpolarizability (4.50 × 106 au) was reproduced. However, according to Dr. P. Karamanis’s comment, the corresponding hyperpolarizability should be on the order of 104. The other reasons would be unearthed for the origin of the difference. The details are discussed as follows.
1. INTRODUCTION We published “How the Number and Location of Lithium Atoms Affect the First Hyperpolarizability of Graphene”1 in J. Phys. Chem. C 2010, 114, 19792. In the paper, a series of graphene nanoribbon multilithium salts Lin@pentacene (n = 1−6) were investigated, in which pentacene was selected as a model of graphene nanoribbon. Results show that first hyperpolarizability (βtot) obviously increases with the increase in the lithium atoms. However, the unexpected first hyperpolarizability (4.50 × 106 au) of Li6@pentacence is 302 times larger than that of Li5@pentacence, which caught Dr. P. Karamanis’s attention. Dr. P. Karamanis and coworker used similar models to repeat our investigations on
[email protected] New configurations (Symmetries in C2v and C2 group) and corresponding first hyperpolarizabilities were obtained. They commented that the βtot of Li6@pentacence should be 1.71 × 104 au, which is lower on the order of 102. We greatly acknowledge the interest and preciseness of Dr. P. Karamanis and coworker. To verify the results, we checked all input and output files. In our original work, all geometrical structures were optimized by the B3LYP method.3 Hereby, the planar configuration of Li6@pentacence (See Figure 1) was obtained, which is very close to the C2v configuration in Dr. P. Karamanis’s work. The βtot was calculated by the BHandHLYP method based on the optimized structure by Gaussian 03
2. CHECK FOR THE STABILITY OF CONFIGURATION Dr. P. Karamanis points out that the configuration of Li6@ pentacence in the original paper is not a genuine minimum due to the existence of imaginary frequencies.2 Their investigation validated that the planar C2v configuration of Li6@pentacence is not stable, which suggests that the distorted structures of the lithiated polycyclic acenes are a relatively stable configuration. After reading the comment, we calculated the frequencies of the planar Li6@pentacence in the original work. Three imaginary frequencies (−126, −84, and −28 cm−1) were found indeed. Therefore, we continued to optimize the geometrical structure and obtained a stable configuration with all real frequencies (the structure shown in Figure 1 and the coordinates shown in the SI). The evaluated βtot at the same theoretical level of the stable configuration is 2.20 × 104 au, which is consistent with the result of Dr. P. Karamanis. Therefore, the failure could be avoided if a careful check on the configuration has been done. However, the βtot of the unstable C2v configuration (in their work) with three imaginary frequencies is still on the order of 104. It means that the error configuration in the original work seems not to be the main reason. More work is needed to dig out other reasons for the error. 3. CHECK FOR THE STABILITY OF NUMERICAL DIFFERENTIATION We also noted that Gaussian 03 revision B.034 was used in our original work to calculate the βtot, but Dr. P. Karamanis’s work was performed by Gaussian 09 Revision B.015 with analytical density functional theory (DFT) method. The analytical DFT method was impossible to calculate in Gaussian 03 revision B.03.4 Therefore, we used the option “polar=enonly” to execute
Figure 1. Structure of Li6@pentacence with three imaginary frequencies (a) and the optimized structure with all real frequencies (b). © XXXX American Chemical Society
Received: September 29, 2012 Revised: November 27, 2012
A
dx.doi.org/10.1021/jp309686h | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Comment
the calculation with finite field method, which introduces a double-numerical differentiation of energies to obtain hyperpolarizability in the original work. The analytical DFT method in Dr. P. Karamanis’s work might be more reasonable to calculate the βtot. (The βtot of the unstable C2v configuration is still on the order of 104.) Mathematically,6 “the main difficulty of numerical computations of the derivatives comes from the ill-posedness, that is, a small perturbation in the measurement data may cause huge error in its derivative by direct computation”. It is suggested that the demerit of the doublenumerical differentiation technique of finite field method might be the core origin of the failure. To validate the hypothesis, we calculate the βtot values of all 14 structures in our original paper with the analytical BHandHLYP method. The results are shown in Figure 2 (data are listed in Table S1 in the SI). The βtot
numerical differentiation technique to break down for the calculation of Li6@pentacence in the original work. The details of the further check are shown as follows. 3.1. Check for the Convergence. We referred to every word in the input file (shown in the SI) of the first hyperpolarizability and then found an option “IOP (5/13=1)” in the command line. However, the calculation without the option for Li6@pentacence terminated abnormally: Illegal MolTyp in RwMol1-- Error termination via Lnk1e in **/l111.exe.8 Thus, the option “IOP (5/13=1)” attracts our attention. According to Gaussian IOP Reference, “IOP (5/ 13=1)” determines the action of the program when it faces a convergence failure.9 The difference of the two output files (with/without the option) is shown in Tables 2 and 3. The calculation normally terminated though three of the finite field solutions did not meet the convergence criteria with the option “IOP (5/13=1)”. The calculation abnormally terminated with 12 converged energies without the option “IOP (5/13=1)”. The energy along the z -axis without the option is not converged under applied fields (−0.0010 au). It is the first time we faced this problem. We learned here that the option “IOP (5/13=1)” should be used cautiously, even though the caused error is generally so small that it could be ignored. 3.2. Check for the Spin-Contaminated Solution. Furthermore, the spin-contaminated solution , has been checked. Results show that the values under static and some applied fields are 0.80 (the value is 0.75 ideally), which means that the spin-contaminant exists in the converged wave function. However, four of the values are close to 0.75 under applied fields, where the wave functions are qualified. As shown in the Table 2, the energies of the less spincontaminated solutions are around 0.03 hartree lower than those of the higher spin-contaminated ones. It shows that the internal instability of wave functions in the original work might the intrinsic reason that causes finite field method failure in the simulation for Li6@pentacence. Normally, the finite field method breaks down when the total energies with respect to the applied field did not vary smoothly. Although the effect from the applied field (0.0010 au) falls in the in the linear range, unfortunately, the unexpected hyperpolarizability of the unstable configuration of Li6@ pentacene also occurred. However, this error was mainly attributed to the inappropriate calculation scheme as well as the insufficient understanding and the stability check of the solutions, which should be avoided. We regret reporting the inaccurately enhanced properties, which were published without further investigation on the enhancement mechanism. The lack in mechanism might increase the breeding place of this type of failure.
Figure 2. Recalculated first hyperpolarizabilities of Lin@pentacene (n = 1−6) by the analytical BHandHLYP method.
values of Lin@pentacence (n = 1−5) evaluated by the analytical method are almost equal to those obtained by doublenumerical differentiation technique. However, a breakdown of double numerical differentiation technique is found for the calculation on the unstable configuration of Li6@pentacence. As shown in Table 1, the result of the analytical7 BHandHLYP Table 1. Static First Hyperpolarizability (βTot) of the Unstable Configuration Calculated by Coupled Perturbed BHandHLYP Method (with the keywords “cphf=rdfreq polar”), the Analytical BHandHLYP Method (with the keywords “polar=analytical”), and BHandHLYP Method with Numerical Differentiation Scheme (with the keywords “polar=enonly”) with Gaussian 09, Revision B.01 keywords
cphf=rdfreq polar
polar=analyticala
βx (au) βy (au) βz (au) βtot (au)
−2.73 × 103 0 1.68 × 104 1.70 × 104
−2.85 × 103 0 1.76 × 104 1.78 × 104
polar=enonly 3.16 3.17 −5.21 4.50
× × × ×
106 106 105 106
4. CONCLUSIONS First, the Li6@pentacene could be explained by a classical donor−acceptor multilithium salt (similar to Lin@pentacene (n = 1−5)), and the evolution of the first hyperpolarizability of Lin@pentacene (n = 1−6) varies close to a linear function of the number of Li atoms. (For n = 1−5, it is almost similar to the original work except for the effect of the sixth Li atom, seen in Figure 2.) Second, the instability of the double numerical differentiation technique of the finite field is the critical factor of the error in our original work; however, this error was mainly attributed to the insufficient check for the employed calculation scheme. We regret inappropriately reporting the unreasonable
a Procedure terminated normally with the option “Int=NoXCTest” in the command line.
method is 1.78 × 104 au for the unstable configuration, that is, 1.70 × 104 au by coupled perturbed BHandHLYP method, and is 4.50 × 106 au by BHandHLYP method with numerical differentiation scheme. The unstable configuration might lead to an unstable wave function or inaccurate quadrature under finite field method. Therefore, the numerical accuracy of XC quadrature might be an important factor that causes the doubleB
dx.doi.org/10.1021/jp309686h | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Comment
Table 2. Total Energies Calculated at UBHandHLYP Level with 6-31+G(d) Basis Set on C and H and 6-311G(3df) Basis Set on Li with the Option “IOP (5/13=1)” Using Gaussian 09, Revision B.01
a
filed+xa
filed+y
filed+z
energy
convergence
s**2
0.0000 0.0010 −0.0010 0.0000 0.0000 0.0000 0.0000 0.0010 −0.0010 0.0010 −0.0010 0.0000 0.0000
0.0000 0.0000 0.0000 0.0010 −0.0010 0.0000 0.0000 0.0010 −0.0010 0.0000 0.0000 0.0010 −0.0010
0.0000 0.0000 0.0000 0.0000 0.0000 0.0010 −0.0010 0.0000 0.0000 0.0010 −0.0010 0.0010 −0.0010
−888.3227091 −888.3245376 −888.3220399 −888.3228017 −888.3228017 −888.3155744 −888.3303042 −888.3246303 −888.3221331 −888.3533515 −888.3576648 −888.3491869 −888.3544168
yes yes yes no yes yes yes no yes no yes yes yes
0.808 0.8086 0.8084 0.808 0.808 0.8075 0.8089 0.8086 0.8083 0.7508 0.7508 0.751 0.7511
Filed+x are the applied field strength along x axes, which are also considered identically in Table 3.
Table 3. Total Energies Calculated at UBHandHLYP Level with 6-31+G(d) Basis Set on C and H and 6-311G(3df) Basis Set on Li without the Option “IOP (5/13=1)” Using G09 rev. B01
a
filed+x
filed+y
filed+z
energy
convergence
s**2
0.0000 0.0010 −0.0010 0.0000 0.0000 0.0000 0.0000 0.0010 −0.0010 0.0010 −0.0010 0.0000 0.0000
0.0000 0.0000 0.0000 0.0010 −0.0010 0.0000 0.0000 0.0010 −0.001 0.0000 0.0000 0.0010 −0.0010
0.0000 0.0000 0.0000 0.0000 0.0000 0.0010 −0.0010a 0.0000 0.0000 0.0010 −0.0010 0.0010 −0.0010
−888.3533515 −888.3233111 −888.3233093 −888.3154702 −888.3304097 −888.3228127
yes yes yes yes yes yes
0.808 0.8083 0.8085 0.8074 0.8089 0.8080
−888.316063 −888.3309969 −888.3234039 −888.3234021 −888.3155636 −888.3305012
yes yes yes yes yes yes
0.8079 0.8093 0.8085 0.8085 0.8075 0.8089
Energy is not calculated due to the convergence failure.
■
calculations and a lack of further investigation on the mechanism.
■
ASSOCIATED CONTENT
S Supporting Information *
The input file for the first hyperpolarizability in the original paper, the reoptimized configuration with all real frequencies, and the static first hyperpolarizability (βtot) calculated by analytical method. This material is available free of charge via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
REFERENCES
(1) Hu, Y.-Y.; Sun, S.-L.; Muhammad, S.; Xu, H.-L.; Su, Z.-M. J. Phys. Chem. C 2010, 114, 19792−19798. (2) Karamanis, P.; Pouchan, C. J. Phys. Chem C 2013, DOI: 10.1021/ jp3057256. (3) Becke, A. D. J. Chem. Phys. 1993, 98, 1372−1377. (4) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A.; Vreven, J., T. ; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision B.03;Gaussian, Inc.: Wallingford, CT, 2004. (5) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; ; Mennucci, B.; ; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H.
ACKNOWLEDGMENTS
We acknowledge the helpful suggestions from a peer reviewer. We also acknowledge the software support (Gaussian 09, revision B.01) from the State Key Lab of Theoretical and Computational Chemistry in Jilin University. C
dx.doi.org/10.1021/jp309686h | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Comment
P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision B.01,Gaussian, Inc.: Wallingford, CT, 2010. To ensure the same Gaussian 09, revision B.01 as Dr. P. Karamanis, we have recourse to State Key Lab of Theoretical and Computational Chemistry in Jilin University. (6) Xu, H.; Liu, J. Adv. Comput. Math. 2010, 33, 431−447. (7) We used the keywords “polar=Analytic BHandHLYP/gen nosymm scf(maxcyc=500)” to execute calculation, but the procedure terminated abnormally: “Inaccurate quadrature in CalDSu.,Error termination via Lnk1e in **/l1002.exe.” Therefore, we used the keywords “polar=Analytic BHandHLYP/gen nosymm scf(maxcyc=500) Int=NoXCTest” to execute calculation. The procedure terminated normally, and the result is 1.78 × 104 au. It is a similar case when the coupled perturbed BHandHLYP method is used for calculations. According to Gaussian reference, the option “Int=NoXCTest” would skip tests of numerical accuracy of XC quadrature, but the stable structure can be calculated normally without the option. It indicates that the unstable configuration might lead to the inaccurate quadrature, which causes the double numerical differentiation technique to lose. (8) ** is the path of our Gaussian revision. (9) Gaussian09 IOps Reference, http://www.gaussian.com/g_tech/ g_iops/ov5.htm, 2009.
D
dx.doi.org/10.1021/jp309686h | J. Phys. Chem. C XXXX, XXX, XXX−XXX
