Spurious Phosphorus Pyramidalization Induced by Some DFT Functionals

The molecular geometries of boraphosphabenzene (BP-benzene) and boron phosphorus coronene (BP-coronene) can be misinterpreted when they are obtained from density functional theory (DFT) calculations. In this study, we found that some exchange-correlation (XC) functionals yielded a distorted geometry of the above molecules when P atoms are present in their resonance structures. This phosphorus pyramidalization may be due to spurious errors caused by using these functionals. To verify this behavior, the electronic structures of BP-benzene and BP-coronene were studied using sixteen functionals (B3LYP, B97D, BHLYP, BP96, PBE, PBE0, PWLDA, Slater-Dirac-exchange, TPSS, M05, M06, M062X, M08HX, M11, wB97 and wB97X-D) with the SVP or TZVPP basis sets. The calculations were carried out using the TURBOMOLE and GAMESS programs. The geometry optimization calculations were carried out for each functional using both of the basis sets. Two different initial geometries, plane (D3h symmetry) and distorted (C1 symmetry) were considered. The optimized geometries of the BP-systems obtained at the MP2/TZVPP and CC2/TZVPP levels of theory exhibited D3h symmetry. These calculations were used as a reference and compared with those obtained from DFT. The optimized geometries obtained from DFT were found to exhibit C1 symmetry for the majority of the XC functionals.


Introduction
We report here a study of benzene and coronene (which comprises six peri-fused benzene rings) with full replacement of the carbon atoms by boron and phosphorus (BP-benzene and BP-coronene).The resultant molecules have resonance structures that are similar to those of their carbon analogues and are isovalent to them. 1 The geometry of BP-benzene (or boraphosphabenzene or s-triphosphatriborin) has been experimentally shown to be planar. 2,3Phosphorus pyramidalization has been described in the literature in several contexts, such as the case of a phosphorus atom inducing curvature in planar polycyclic hydrocarbons. 4Other examples include an observed decrease in the pyramidalization at the phosphorus atom by an increase in the metal-phosphorus bond length 5 and phosphorus pyramidalization in the diphosphine radical cation. 6he coronene molecule is considered to be a simple model for describing graphene.Graphenes have received much attention from the scientific community 7 because of their optoelectronic properties and potential applications in devices such as sensors.The particular use of graphene systems in electronic devices is determined by the size of their band-gap.Because graphene itself has no band-gap, it is necessary to modify its structure with vacancies, impurities or defects to create one.An alternative approach would be to use systems similar to graphene such as BP-graphene. 8hemical functionalization of coronene and corannulene also shows promise for electronics applications because this process allows tunable optoelectronic properties and improves charge-carrier mobility. 9In general, calculations have indicated that a class of 2D group-III/V (boron nitride (BN), boron phosphide (BP), aluminum nitride (AlN), and aluminum phosphide (AIP)) systems can also exhibit magnetic properties (ferromagnetism and half metallicity) induced via charge injection. 10he use of density functional theory (DFT) calculations in studying coronene derivatives (with BN, BP, AlN, etc.) has been frequently reported in the literature.Some examples include studies of the carrier-transport properties of coronene derivatives, 11 of their interactions with beryllium, magnesium and calcium, 12 as well as interactions between coronene derivative dimers. 13Additional studies Vol. 26, No. 8, 2015   have dealt with their bond dissociation energies, 14 calculation of their thermal and magnetic properties, 15 simulation of their vibrational spectra 16 and calculations of atomic hydrogen chemisorption. 17he first step toward understanding the properties of these systems is to determine their molecular geometries.In this work, we observed that the molecular geometry has a very strong dependence on the exchange-correlation (XC) functional 18 used in the calculation.The optimized geometries of the BP-systems obtained from calculations with correlation effects exhibit D 3h symmetry.However, molecular structures with C 1 symmetry can be obtained when any type of DFT exchange-correlation functionals are employed. 1,8here are several DFT-based theoretical studies in the literature on molecules that contain P atoms in their resonance structures.The flat structure leads to a higher electron delocalization by the system.And this resonance effect further stabilizes the molecular system.5][6] To investigate this behavior, the BP-benzene and BP-coronene systems were studied using sixteen different exchangecorrelation functionals.

Methodology
The methodology employed in this work is similar to that presented in previously published reports on this subject. 1,85][36][37][38][39] The chosen exchange-correlation functionals are classified as local, pure or hybrid DFT GGA, and pure or hybrid DFT meta-GGA (these XC functionals are constructed considering empirical fitting of their parameters, otherwise non-empirical (ab initio) approaches where all parameters are fixed via exact conditions).The Grimme's functional, including dispersion effects, was also considered in this study. 40The calculations performed at the DFT level used two different function basis sets, SVP and TZVPP. 41he molecular geometry optimizations were carried out without symmetry restrictions starting from the two initial geometries (D 3h and C 1 ) of BP-benzene and BP-coronene.All of the calculations were carried out using the TURBOMOLE (B3LYP, B97D, BHLYP, BP96, PBE, PBE0, PWLDA, Slater-Dirac-exchange and TPSS) 42,43 and GAMESS (M05, M06, M062X, M08HX, M11, wB97 and wB97X-D) 44 software packages.
The equilibrium geometries were checked by analyzing the harmonic vibrational frequencies, which are crucial for verifying that the obtained molecular geometries are local minima on the potential energy surface (PES).The optimization procedure stops when the convergence criteria (a change in the total energy of 10 -6 a.u. and a maximum norm of the gradient of 10 -3 a.u.) are obtained.The numerical integration for DFT calculation was carried out considering default precision of the TURBOMOLE, i.e., the grid step was in average 10 -7 a.u.This default grid provides sufficient numerical accuracy for DFT energies, geometries and frequencies.
The CC2 ground-state energy (E CC ) 45,46 is obtained from equation 1, where the cluster operator T is expanded to only single and double excitation operators, i.e., T ≈ T 1 + T 2 .The single and double excitation are given by equation 2 and equation 3, where a c and a k are creation and annihilation operators, respectively, and summation Einstein notation has been used.The indices (i,j) are associated with occupied orbitals and indices (a,b) are associated with unoccupied orbitals.The approximate solution is obtained by finding the unknown coefficients (cluster amplitudes) t i a and t ij ab .The MP2 (second order Møller-Plesset perturbation) 47 approximation is implemented by setting t i a = 0 in the CC2 equations.

Results and Discussion
The initial geometries, with C 1 and D 3h symmetries, were fully optimized.The optimization procedures were performed with the sixteen different exchange-correlation functionals using the SVP and TZVPP basis sets.At the MP2 and CC2 levels of theory, only the TZVPP basis set was employed.The molecular geometries obtained at the MP2 and CC2 levels were employed as reference values.The most important aspect of the initial geometry is its symmetry group.In general, the DFT calculation converges to the same symmetry of the initial geometry considering the procedure previously described.In this section, we first propose a warning that the most stable structure obtained by DFT calculations does not always coincide with the correct one.We then discuss why the DFT calculations converge to two energy minima, which can depend on the initial geometry.
In general, a flat molecular geometry can be obtained by DFT calculations for molecules with resonance structures.
The DFT calculations for these molecules eventually lead to stable structures with a distorted molecular geometry, i.e., a local minimum on the PES; however, as discussed below, this geometry is not correct.
The spurious distortion is not present in the pure carbon benzene or in borazine (BN-benzene) molecules.For comparison, several calculations beginning from a distorted geometry were performed for benzene and BN-benzene, resulting in dihedral angles of less than 0.03º and 0.3º, respectively, considering the same convergence criteria mentioned above.When the geometry optimization criteria were changed from 10 -7 to 10 -4 , change in the dihedral angles for benzene was not observed, but for BN-benzene the dihedral angles were decreased by 0.09º.

BP-benzene
A preliminary study was carried out with BP-benzene to understand the phosphorus pyramidalization.In general, the DFT calculations for BP-benzene revealed the presence of two minima, corresponding to planar and distorted structures.The two minima were discovered when considering two different types of initial geometries.Tables S1 and S2 in the Supplementary Information (SI) show the planar (obtained from CC2/TZVPP) and distorted (obtained from DFT/B3LYP/TZVPP) initial geometries, respectively, obtained from the DFT optimization procedure.
Table 1 contains the energy differences (the absolute energies are provided in the Table S3) and torsion angles (given as average values) of the optimized geometries, obtained using the two different initial geometries (D 3h and C 1 ).Nearly all of the converged geometries (specifically, those obtained by all of the XC functionals except for M05 and M06) represent electronic states that are local minima: i.e., all of the calculated vibrational frequencies are real positive values.The M05 and M06 functionals each exhibit three imaginary frequencies for the planar structures in the intervals [114.97:36.03]and [117.23:64.53],respectively.All of the BP-benzene results were obtained using the TZVPP basis set.The torsion angles are displayed in Figure 1 as average values.
The histograms displayed in Figures 2 and 3 summarize the results of Table 1.The results in Figure 2 show the final energy differences between the flat and distorted initial geometries.If the difference is positive, the geometry obtained from an initial flat structure is less stable.Note that a distorted final geometry was obtained from the M08HX, wB97 and wB97X-D functionals, though the initial geometry was flat.Figure 3 shows the final torsion angle for a distorted initial geometry.The final distorted geometries are characterized by the presence of all hydrogen atoms on the same side of the molecule and a zig-zag conformation of the BP-ring.
Only the B97D, BHLYP, TPSS, M05 and M06 functionals diverge in the calculation of the final energy    1).
differences.When the initial geometry is flat, all of the XC functionals converge to a flat symmetry (the angle of 0.03° can be assumed to be a numerical error) except for M08HX, wB97, and wB97X-D, which yielded torsion angles of 37, 34 and 9°, respectively.The calculation using the M11 functional also does not rigorously converge to a flat geometry because the torsion angle is 0.15°.Although no energy difference was observed between the flat and distorted geometries obtained using the B3LYP functional, the initial and final torsion angles are very different.In other words, the two structures have different symmetries but are nearly degenerate in energy.The B97D, BHLYP, M05 and M06 functionals yielded significant energy differences between the flat and distorted geometries of 0.003, 0.02, 0.07 and 0.05 eV, respectively, as well as structures exhibiting high torsion angles of 18, 32, 38 and 35°, respectively.Because the molecule must converge to the flat conformation, these results indicate one possible shortcoming of the XC functionals.Moreover, the positive energy differences indicate an unexpected stability of the molecules with distorted geometries.
For the BP96, PBE, PBE0, PWLDA, Slater, M062X and M11 functionals, the calculated energy differences are very small, which is consistent with the very small observed torsion angles (less than 2°).Therefore, these seven XC functionals perfectly represent the BP-benzene system: i.e., a flat molecule.In particular, the two geometries obtained with the TPSS functional are stable, meaning that there are no observed imaginary vibrational frequencies.Additionally, the negative energy difference indicates that the flat geometry has the lower energy and is therefore more stable.
In summary, the B3LYP, B97D, BHLYP M05, M06, M08HX, wB97 and wB97X-D functionals do not adequately describe the geometry of the BP-benzene system.When considering both of the initial distorted geometries, these XC functionals yielded high torsion angles compared with those obtained from the other XC functionals.
Because similar calculations on benzene and BN-benzene converge to a flat geometry, we can infer that the presence of phosphorus in the molecule causes a spurious distortion for some of the XC functionals.This result is independent of whether the initial geometry is flat or distorted.In addition, if the geometry optimization criteria are considered more rigorous by one order of magnitude, the obtained dihedral angles for BP-benzene remain almost unchanged.

BP-coronene
Both flat (obtained with CC2/TZVPP) and distorted (obtained with B3 LYP/TZVPP) initial geometries were considered in the calculations on BP-coronene.Tables S4 and S5 in the SI display the flat and distorted initial geometries, respectively.The relative energies (the absolute energies are shown in the SI, Table S6) of the fully optimized geometries obtained with the sixteen XC functionals using the SVP and TZVPP basis sets are shown in Tables 2 and 3, respectively.Additionally, some results obtained using the QZVPP basis set are presented in Table S7 of the SI.The torsion angles between the external BPBP atoms can be found in these tables, as well as in Figure 4.
Tables 2 and 3 show the effect of changing the basis set functionals on the molecular geometries.For BP-coronene, the PBE functional yields either a flat or a distorted structure when using the SVP or TZVPP basis sets, respectively.The same behavior is observed for the M062X functional.The primary difference observed between the SVP and TZVPP basis sets is due to their number of d-functions or f-functions.SVP basis sets do not have any f-functions and contain a smaller number of d-functions.In this study, although these two XC functionals clearly provided good results when using a small basis set, they may be fortuitous.By contrast, when using the M06 and M11 functionals, the   The histogram in Figure 5 displays the final energy differences between the two converged geometries obtained using the sixteen XC functionals and two basis sets.The torsion angle (Figure 4) obtained from calculations with the distorted initial geometry can be observed in the histogram shown in Figure 6.Histograms that include the relative energies and torsion angles obtained from the QZVPP basis sets are shown in the SI as Figures S1 and S2, respectively.
When the calculations were carried out from a flat initial geometry, all of the XC functionals retained the flat geometry except for M05, M06, M062X, M08HX, M11, wB97 and wB97X-D.However, for most of the XC functionals considered in this work, this geometry does not correspond to the energy minimum.When the calculations were initiated from a distorted geometry, most of the XC functionals also yielded a distorted molecular geometry.When this occurs, the resulting distorted geometry corresponds to one absolute minimum of energy.
contrast to BP-benzene, the converged flat coronene structure obtained using DFT functionals implemented   2 and 3).  2 and 3).For the SVP basis set, the energy difference between the distorted and flat geometries is in the range of 0.05 to 0.17 eV for the B3LYP, B97D, BHLYP and M05 functionals.For the BP96, PBE0, TPSS, M08HX, M11 and wB97X-D functionals, the energy difference is approximately 10 -3 eV.The other XC functionals (PBE, PWLDA, Slater, M062X and wB97) correctly describe the minimum energy as corresponding to a flat conformation; therefore, there is no significant energy difference for these functionals.
When using a better atomic basis, TZVPP, the energy differences between the flat and distorted geometries tend to increase.The QZVPP basis sets provide similar results to those obtained from TZVPP (see SI, Figure S1).For example, the energy difference between the distorted and flat geometries is in the range of 0.06 to 0.16 eV for the B3LYP, B97D, BHLYP, M05, M06, M08HX, wB97 and wB97X-D functionals, while for the BP96, PBE0, TPSS, M062X and M11 functionals, it is approximately 0.01 eV.For the PBE functional, the energy difference is 0.002 eV.The rest of the XC functionals (PWLDA and Slater) correctly describe the minimum energy as corresponding to a flat geometry.Therefore, as was the case when using SVP, there is no significant energy difference between the distorted and flat geometries calculated from these two XC functionals.
The increment in the energy differences observed for the TZVPP (and QZVPP, see SI) basis sets indicates that the problem of convergence for the distorted geometries is related to the XC functionals rather than the basis set because the differences in energies and torsion angles increase when improving the latter.Although the PBE and M062X functionals yield a flat final geometry for the SVP basis, the distortion problem appears upon improving the basis set.Therefore, the fact that the calculations converge to a distorted geometry for some XC functionals indicates a possible error related to the functionals themselves.
We also note that for the BP96, PBE, PBE0, TPSS, M062X and M11 functionals, the different geometries have very close energies, which may also indicate that the problem lies in the XC functionals.
The molecular geometry converges to a flat structure for the isovalent systems coronene and BN-coronene.Therefore, phosphorus, which is a third-row element, cannot be correctly described for some XC functionals.Alternatively, despite its similarity to the other two coronenes, the majority of XC functionals do not adequately describe the energies and the geometries of this system due to the presence of the P atom in the resonance structures.
Therefore, there is a likely spurious pyramidalization of the phosphorus because the optimized geometry of BP-coronene obtained through more sophisticated calculations such as MP2/TZVPP and CC2/TZVPP converges to D 3h symmetry.The calculated energy and final geometry obtained for BP-benzene and BP-coronene using CC2 and MP2 are provided in the SI in Tables S1 and S4 and Tables S8 and S9, respectively.
In summary, the main DFT functionals provide a distorted structure (C 1 ) rather than a planar geometry (D 3h ) as the global minimum.Most importantly, the majority of the functionals do not yield the correct, planar, structure when P atoms are present (e.g., in BP-coronene).The efficiency of devices based on these chemical systems is determined by the mobility of the charge carriers, i.e., how quickly a particle carrying a charge moves across the active layer upon application of an external electric field. 48,49Thus, it is important that the theoretical model correctly describes the molecular geometry.For example, a distorted geometry predicted by a theoretical calculation can lead to estimates of electronic properties that markedly differ from those obtained experimentally.

Conclusions
In contrast to the optimized geometry of BP-coronene (D 3h ) obtained through more sophisticated calculations (e.g., CC2/TZVPP), the majority of XC functionals converge to a distorted geometry corresponding to a minor energy minimum (i.e., more stable).
This finding indicates that for BP-coronene, owing to the presence of phosphorus the majority of the XC functionals studied do not appropriately describe the molecular structure, which appears to exhibit spurious pyramidalization.
Calculations performed on BP-benzene exhibited a similar problem, although to a lesser extent.In this particular case, the problem of distortion appeared for eight (B3LYP, B97D, BHLYP, M05, M06, M08HX, wB97 and wB97X-D) of the sixteen XC functionals employed.
Many theoretical studies of molecules containing a phosphorus atom (when considering the XC functionals studied here) can be found in the literature and based on the described results may be incorrect.We have shown in this study that geometrical distortions are obtained for some XC functionals, regardless of the basis set used.Additionally, the pyramidalization obtained for the phosphorus atom in the resonance systems is likely to be spurious, and it may be an intrinsic error of some of the XC functionals.
We have also verified that for BP-benzene, the final two converged geometries resemble two stable minima, which are degenerate in the case of B3LYP.For BP-coronene, only one stable minimum is obtained for the B3LYP functional, corresponding to the distorted geometry. 1or BP-benzene and BP-coronene, the majority of the XC functionals employed does not produce either the experimental geometry or the theoretical geometry obtained through more sophisticated methods.The effective presence of d orbitals on the atoms makes the barrier involved in the breaking of symmetry much smaller. 50,51his probably occurs because of the known shortcomings of the XC functionals used.We have observed the sensitive dependence of the calculated molecular geometry on the XC functional for BP-systems.Based on our findings, we recommend that careful studies be performed prior to choosing the type of density functional to be used in subsequent calculations.This issue is an intrinsic problem of the DFT methodology.

Figure 1 .
Figure 1.The labeled atoms delimit the torsion angle.

Figure 3 .
Figure 3. Final torsion angle for the structures obtained from an initial distorted geometry (see Table1).
Figure 3. Final torsion angle for the structures obtained from an initial distorted geometry (see Table1).

Figure 4 .
Figure 4. Torsion angles obtained based on the marked atoms.

Figure 5 .
Figure 5. Energy differences (in 0.10 eV) between the two converged geometries obtained from the various XC functionals (see Tables2 and 3).

Figure 6 .
Figure 6.Torsion angles for the final distorted geometries obtained from the various XC functionals (see Tables2 and 3).

Table 1 .
Torsion angles and final energy differences (ΔE) between the calculated BP-benzene structures based on a flat initial geometry and an initial distorted geometry a Flat initial geometry; b distorted initial geometry.

Table 2 .
Torsion angles and energy differences between the flat and distorted structures obtained from several XC functionals using the SVP basis set for BP-coronene a Flat initial geometry; b distorted initial geometry.

Table 3 .
Torsion angles and energy differences between the flat and distorted structures obtained from several XC functionals using the TZVPP basis set for BP-coronene converges to a distorted structure for the SVP basis set and to a nearly flat structure for the TZVPP basis set when the calculation begins from a flat geometry.in the TURBOMOLE software contains imaginary harmonic vibrational frequencies (if).Specifically, B3LYP contains six ifs in the interval [74.16:42.57],B97D contains six ifs in [90.07:59.47],BHLYP contains six ifs in [74.06:43.20],BP96 contains five ifs in [49.97:17.66],PBE contains one if = 20.97,PBE0, PWLDA and Slater-Dirac exchange contain zero ifs and TPSS contains two ifs in [29.26:8.69].These results are in disagreement with the more sophisticated calculations performed at the MP2 and CC2 levels of theory.
a Flat initial geometry; b distorted initial geometry; c irregular geometry.calculation