Reliability of the density functional approximation to describe the charge transfer and electrostatic complexes involved in the modeling of organic conducting polymers

Carlos Alemán,* David Curcó, and Jordi Casanovas Departament d’Enginyeria Química, E.T.S. d’Enginyers Industrials de Barcelona, Universitat Politècnica de Catalunya, Diagonal 647, Barcelona E-08028, Spain Departament d’Enginyeria Química, Facultat de Química, Universitat de Barcelona, Marí i Franques 1, Barcelona E-08028, Spain Departament de Química, Escola Politècnica Superior, Universitat de Lleida, c/ Jaume II No. 69, Lleida E-25001, Spain Received 30 March 2005; revised manuscript received 13 May 2005; published 17 August 2005


I. INTRODUCTION
The interactions of -conjugated molecules, such as pyrrole and thiophene, with metallic atoms either in neutral or positively charged states ͑M and M n+ , respectively͒ are essential in the fields of polymer and material sciences.Thus M ¯ and M n+ ¯ interactions are involved in the doping of organic conducting polymers.Within this context, model complexes containing such interactions have been satisfactorily used to model the doping of organic conducting polymers through quantum mechanical calculations ͓1-7͔.
Density functional theory ͑DFT͒ offers a very interesting alternative to conventional ab initio methods, such as MP2 and MP4, for describing the interaction of metallic species with -conjugated molecules.The most reliable advantage of DFT with respect to these sophisticated ab initio methods is that the former includes correlation effects at a lower computational cost, i.e., the cost is N 3 , N 5 , and N 7 for DFT, MP2, and MP4, respectively, where N is the number of basis functions.However, the ability of DFT calculations to describe noncovalent interactions has been shown to depend on the form of the functional.In the last decade, the suitability of the more popular functionals to predict intermolecular interaction energies and geometries for hydrogen bonded ͓8-11͔, van der Waals ͓12-14͔, and aromatic − ͓15-17͔ complexes has been checked.
In this work we evaluate the performance of the B3LYP, B3PW91, PBE, and MPW91PW91 functionals to describe both M ¯ and M n+ ¯ interactions.Calculations have been performed on systems formed by thiophene ͑Thp͒ and pyrrole ͑Py͒ interacting with M = Li, Na, K, Ca, and Mg; and M n+ =Li + , Na + , K + , Mg 2+ , and Ca 2+ .Results provided by DFT calculations have been compared to those obtained at the HF, MP2, and MP4 levels.

II. METHODS
All ab initio and DFT calculations were performed using the GAUSSIAN 98 ͓18͔ program.Ab initio calculations were carried out at the HF, MP2, and MP4 ͓19͔ levels of theory.DFT calculations were performed using four combinations: the Becke's three-parameter hybrid functional ͑B3͒ ͓20͔ with the Lee, Yang, and Parr ͑LYP͒ ͓21͔ expression for the nonlocal correlation ͑B3LYP͒; the same functional with the nonlocal correlation provided by Perdew and Wang ͑B3PW91͒ ͓22͔; the approximation of Perdew, Burke, and Ernzerhof ͑PBE͒ ͓23,24͔, which is a simpler version of PW91; and the modified Perdew-Wang exchange ͓25͔ with the PW91 gradient-corrected correlation ͓22͔ ͑MPW1PW91͒.Calculations on closed-shell and open-shell systems were done considering the restricted and unrestricted quantum-chemical formalisms, respectively.Geometry optimizations were performed at the HF, MP2, B3LYP, B3PW91, PBE, and MPW1PW91 levels of theory, while MP4 energies ͑including all single, double, triple, and quadruple excitations͒ were derived from single point calculations on geometries optimized at the MP2 level.The 6-31+ G͑d , p͒ basis set ͓26͔ was used in all the ab initio and DFT calculations.Frequency analyses were carried out to verify the nature of the minimum state of all the stationary points located during geometry optimizations.
The interaction energies ͑⌬E int ͒ were calculated as the difference between the energy of the optimized complex and the sum of the energies of the isolated fragments ͑-conjugated ring and metal atom͒.The counterpoise method ͓27͔ ͑CP͒ was applied to correct the basis set superposition error ͑BSSE͒ in the interaction energies.
In order to compare the molecular geometries predicted by the different computational methods, the same starting arrangement was considered for all the complexes.This con-

A. Interaction energies
The interaction energies with CP corrections calculated for M ¯Thp and M ¯Py charge transfer complexes are listed in Table I.As expected MP2 and MP4 interaction energies are remarkably similar in all cases.The largest difference ͑2.8 kcal/ mol͒ corresponds to the Li¯Thp complex, where the inclusion of electron correlation effects up to fourth-order results in a counterbalancing of the excess of stability achieved at the MP2 level.For the rest of the complexes MP2 interaction energies are within 0.3 kcal/ mol of MP4.Comparison between the MP4 and HF energies points out the limitations of the latter method indicating that inclusion of electron correlation effects is essential.Thus the HF method provides a notable underestimation of the interaction energies for M ¯Py complexes with M = Li, Na, and K; while no stabilization was achieved for M ¯Py with M = Mg, Ca and M ¯Thp with M Li, i.e., geometry optimization led to a separation of the fragments.Interestingly, the Li¯Thp was the only complex satisfactorily reproduced at the HF level.
On the other hand, comparison between DFT and MP4 interaction energies reveals notable discrepancies, which in turn depend on the functional used in the former calculations.The B3LYP and B3PW91 methods, especially the former one, underestimate the stability of M ¯Thp with M = Na, K and M ¯Py with M Li, and neglect the existence of M ¯Thp with M = Mg, Ca.Furthermore, these functionals reproduce satisfactorily the interaction energy of the Li¯Py complex, while that of Li¯Thp is drastically overestimated.The PBE and MPW1PW91 methods provide the best DFT estimates.Thus the interaction energies predicted these approximations for M ¯Thp and M ¯Py complexes with M Li are in good agreement with the MP4 ones.However, the values predicted for Li¯Thp and Li¯Py are overestimated by 65-75% and 41-48%, respectively.The overall of these results indicate that the PBE and MPW1PW91 are the better DFT methods adequate to describe the energetics involved in M ¯ complexes.
The interaction energies computed at the ab initio and DFT levels for M n+ ¯Thp and M n+ ¯Py complexes are shown in Table II.Comparison of results listed in Tables I  and II shows that the interaction energies are much more negative when the metal atom is ionized.Thus M n+ ¯Thp and M n+ ¯Py are dominated by a strong electrostatic interaction between the cloud of the heterocycle and the positively charged atom.Accordingly, at the highest level of theory the interaction energies range from −30.8 to −13.6 kcal/ mol and from −37.2 to −18.0 kcal/ mol when a singly charged atom M + is complexed with Thp and Py, respectively.These values are one order of magnitude more favorable than those shown in Table I for M ¯Thp and M ¯Py.On the other hand, when the metallic specie presents two positive charges, as in M 2+ ¯Thp͑Py͒ complexes, the interaction energy amounts to −104.2͑−113.1͒and −59.6͑−68.8͒kcal/ mol for M 2+ =Mg +2 and Ca 2+ , respectively.
Analysis of electron correlation effects in Table II indicates that the change from the MP2 to MP4 level leads to very small changes in the interaction energy ͑typically around 0.4 kcal/ mol, i.e., ϳ1%͒ suggesting a good convergence.On the other hand, the HF method tends to overestimate the strength of the interaction in these electrostatically stabilized complexes by only ϳ2.5%, which is a very satisfactory result.Finally, the four functionals investigated provide a reasonable description of the complexes containing a metal cation, although in all cases the strength of the interaction is overestimated with respect to the ab initio values.This overestimation is about 9% for both B3LYP and B3PW91 functionals, while the MPW1PW91 and, particularly, the PBE methods provide the more negative values of the interaction energies, i.e., the overestimation is of about 11%.The similarity between B3LYP and B3PW91 results combined with the enhancement produced by the PBE and MPW1PW91 methods suggest that the exchange functionals are responsible for the DFT overestimation.
The BSSE was corrected using the CP method.
b All calculations were carried out using the 6 -31+ G͑d , p͒ basis set.The BSSE was corrected using the CP method.
b All calculations were carried out using the 6-31+ G͑d , p͒ basis set.

B. Geometries
The intermolecular parameters used to compare the geometries predicted by the different methods for the ten complexes investigated consist of three distances d 1 , d 2,5 , and d 3,4 , which are defined in Fig. 1.It should be noted that we are interested in how the strength predicted by the different methods for the M ¯ and M n+ ¯ interactions correlates with the geometry of the complexes, this information being satisfactorily described by the three parameters considered.Thus other intermolecular parameters as angles or dihedrals, which can be supplied upon request to authors, are not essential in this case because the relative arrangement between the two fragments of the complex was the same in all cases, i.e., the M or M n+ specie above the -system of the heterocycle.
Table III summarizes the optimized geometric parameters for the ten charge transfer complexes.The geometries obtained using the MP2 method are fully consistent with the interaction energies predicted at the same computational level.Thus the three distances increase with the size of M for both M ¯Thp and M ¯Py complexes, and for a given M the distances are larger for M ¯Thp complexes than for M ¯Py complexes.Furthermore, as expected for a weak interaction, the parameters are relatively large ranging from 2.151 to 4.063 Å for M ¯Thp and from 2.521 to 3.483 Å for M ¯Py.
On the other hand, the HF parameters reflect the limitations of this method to describe charge transfer complexes.Thus the resulting HF intermolecular distances show that, in many cases, geometry optimizations led to the disruption of the complex.The four DFT methods tend to overestimate the intermolecular distances with respect to the MP2 values.The largest overestimation was obtained at the B3LYP level.Substitution of the LYP functional by the PW91 one to express correlation induces an improvement in the geometries.As expected, the PBE and MPW1MPW91 provide the best DFT results, the intermolecular distances obtained at these levels being ϳ0.1 Å greater than the MP2 ones.The results are in excellent agreement with the interaction energies displayed in Table I.
Table IV  ¯ at the MP2 level indicates a shortening of ϳ0.3 and ϳ1.0 Å for M + and M 2+ , respectively.On the other hand, the geometries calculated at the different levels for M n+ ¯Thp and M n+ ¯Py complexes show a remarkable agreement.Thus HF and DFT methods are able to describe satisfactorily the geometry of the electrostatically stabilized complexes.Accordingly, although the shortcomings of the exchange functionals pro- duced an overestimation of the corresponding interaction energies ͑Table II͒, they do not affect the geometries.

C. Analysis of the basis set superposition error
The importance of the BSSE in the computed interaction energies has been calibrated by comparing the results displayed in Tables I and II with those listed in Tables V and VI, which show the interaction energies calculated without applying the CP correction.We are aware that the CP method provides only approximated results since it overestimates the size of the BSSE.However, early studies on van der Waals and hydrogen bonded complexes have shown that the CP correction provides interaction energies within the experimental uncertainty ͓28-30͔.
Inspection to the results obtained for M ¯Thp and M ¯Py complexes indicates that the BSSE plays a similar but important role in MP2 and MP4 calculations.Thus the BSSE ranges from 1.2 to 4.0 kcal/ mol for M = Li, Na, and K and from 0.9 to 2.9 kcal/ mol for M = Mg and Ca, these values being about 32-67% of the interaction energy.DFT methods seem to be less susceptible to the basis set truncation than conventional MPn methods for M ¯ complexes.However, this conclusion must be taken with caution since, as was shown above, these complexes are poorly described by DFT methods, especially by B3LYP and B3PW91.The BSSE of the PBE and MPW1PW91 functionals, which provide the more reliable DFT results for M ¯ complexes, range from 0.1 to 0.8 kcal/ mol that corresponds to less than 22% of the interaction energy.Similarly, the BSSE predicted by the four functionals investigated for Li¯Thp and Li¯Py complexes was lower than 12% of the interaction energy.The overall of these results suggest that, in spite of their limitations, the DFT methods are less sensitive to BSSE than conventional MP2 and MP4 methods.However, the magnitude of the BSSE oscillates with the different functionals.
The uncorrected interaction energies calculated for M n+ ¯Thp and M n+ ¯Py at the ab initio MP2, MP4, and HF levels ͑Table VI͒ reveal that the BSSE is considerably smaller than for M ¯ complexes.Moreover, the size of the BSSE decreases when the electrostatic component of the interaction increases.Thus the BSSE is smaller for M n+ ¯Thp than for M n+ ¯Py and, within each group, the influence of the basis set truncation is smaller for M 2+ -than for M + -containing complexes.Examination of the DFT results indicates that, independently of the functional, the BSSE is almost negligible for the electrostatically stabilized complexes.Thus the BSSE is around 2-4% for the complexes TABLE V. BSSE-uncorrected interaction energies a ͑in kcal/ mol͒ calculated for charge transfer complexes M¯thiophene and M¯pyrrole.In every case the level of theory used to calculate the interaction energy is identical to that used to obtain the optimized geometry with the exception of the MP4 level, where the MP2 optimized geometry was used ͑see Sec.formed with M + and smaller than 1% for those containing M 2+ .

IV. CONCLUSIONS
In this work we investigated M ¯ and M n+ ¯ complexes using DFT and ab initio methods.The performance of four DFT functionals ͑B3LYP, B3PW91, PBE, and MPW1PW91͒ was tested considering 20 different complexes.Comparison with MP2 and MP4 results indicated that the B3LYP and B3PW91 methods are completely inadequate in prediction of M ¯ interaction energy and geometry.However, the MPW1PW91 functional provided a good description of all the charge transfer complexes investigated.Furthermore, this DFT method predicts reliable interaction energies even when the molecular geometries were optimized using a different computational procedure.The BSSE estimated for M ¯Thp and M ¯Py was found to be consid-erably smaller for DFT methods than for conventional ab initio ones.Interestingly, the failure of the B3LYP and B3PW91 was also reported for electron donor acceptor systems formed from ethylene interacting with a halogen molecule ͓31͔, where the interaction is expected to be weaker than in the M ¯ complexes presented in this work.
On the other hand, electrostatically stabilized M n+ ¯ complexes are rightly described by all the methods considered in this work, although DFT methods overestimate the strength of the interaction with respect to the ab initio ones.In general the size of the BSSE was found to be small in M n+ ¯Thp and M n+ ¯Py complexes, being almost negligible for the four DFT methods.

FIG. 1 .
FIG. 1. Structure of the complexes investigated in this work.The geometric parameters ͑d 1 , d 2,5 , and d 3,4 ͒ used to compare the geometries provided by the different methods ͑see text͒ are defined.

TABLE I .
Interaction energies a ͑in kcal/ mol͒ calculated for charge transfer complexes M¯thiophene and M¯pyrrole.In every case the level of theory used to calculate the interaction energy is identical to that used to obtain the optimized geometry with exception of the MP4 level, where the MP2 optimized geometry was used ͑see Sec.II͒ b

TABLE II .
Interaction energies a ͑in kcal/ mol͒ calculated for complexes M n+ ¯thiophene and M n+ ¯pyrrole.In every case the level of theory used to calculate the interaction energy is identical to that used to obtain the optimized geometry with exception of the MP4 level, where the MP2 optimized geometry was used ͑see Sec.II͒.b

TABLE III .
Optimized intermolecular parameters ͑in Å͒ for charge transfer complexes M¯thiophene and M¯pyrrole.The parameters are defined in Fig. 1.

TABLE IV .
Optimized intermolecular parameters ͑in Å͒ for charge transfer complexes M n+ ¯thiophene and M n+ ¯pyrrole.The parameters are defined in Fig. 1.

TABLE VI .
BSSE-uncorrected interaction energies a ͑in kcal/ mol͒ calculated for charge transfer complexes M n+ ¯thiophene and M n+ ¯pyrrole.In every case the level of theory used to calculate the interaction energy is identical to that used to obtain the optimized geometry with exception of the MP4 level, where the MP2 optimized geometry was used ͑see Sec.II͒. a a All calculations were carried out using the 6 -31+ G͑d , p͒ basis set.