WWW.C-CHEM.ORG
FULL PAPER
Application of Polarizable Ellipsoidal Force Field Model to Pnicogen Bonds Fang Liu,[a] Likai Du,*[b] Jun Gao,*[a] Lili Wang,[a] Bo Song,[c] and Chengbu Liu[a] Noncovalent interactions, such as hydrogen bonds and halogen bonds, are frequently used in drug designing and crystal engineering. Recently, a novel noncovalent pnicogen bonds have been identified as an important driving force in crystal structures with similar bonding mechanisms as hydrogen bond and halogen bond. Although the pnicogen bond is highly anisotropic, the pnicogen bond angles range from 160 to 180 due to the complicated substituent effects. To understand the anisotropic characters of pnicogen bond, a modification of the polarizable ellipsoidal force field (PEff ) model previously used to define halogen bonds was proposed in this work. The potential energy surfaces (PESs) of mono- and poly-
substituted PH3–NH3 complexes were calculated at CCSD(T), MP2, and density functional theory levels and were used to examine the modified PEff model. The results indicate that the modified PEff model can precisely characterize pnicogen bond. The root mean squared error of PES obtained with PEff model is less than 0.5 kcal/mol, compared with MP2 results. In addition, the modified PEff model may be applied to other noncovalent bond interactions, which is important to understand the role of intermolecular interactions in the self-assembly C 2015 Wiley Periodicals, Inc. structures. V
Introduction
180 . In addition, the r-hole, a concept proposed by Politzer and Murray, can also be found in pnicogen bond.[26,27] Some properties of pnicogen bond in a few model complexes, such as bond length and binding energy have been theoretically determined through quantum mechanics calculations.[24,28] The binding energies of pnicogen bonds and halogen bonds are 1–10[13] and 1–45 kcal/mol,[29] respectively. Theoretically, the pnicogen bond is stabilized partially by electrostatics and polarization. In addition, charge-transfer also contributes to the stability of pnicogen bonds.[13] The molecular mechanics (MM) method is a very useful tool to study condensed systems. However, the traditional force fields that assign simple isotropic point charges to atoms cannot be applied to pnicogen bond complexes due to their unique structural and thermodynamic properties. The traditional point charge model-based force fields, such as MM3,[30]
Noncovalent interactions play important roles in the functions and properties of condensed phases, biomolecular crystals and solutions.[1–4] The hydrogen bond (H-bond) is one of the most well known and understood noncovalent interactions.[5–8] The substitution of the H atom of H-bond by a halogen atom leads to another kind of noncovalent interaction, the halogen bond (X-bond).[9–12] Recent work has shown that it is not only hydrogen and halogen atoms can interact with an electron-donor, but also adjacent main-group elements, such as pnicogen and chalcogen atoms.[13,14] The noncovalent pnicogen bond was first found in crystal structures. They are as strong as hydrogen bonds and can act as a molecular linker in different chemical systems.[15] Since then, the pnicogen bond has received considerable attentions. Tuononen and coworkers[16] showed that the pnicogen bonding interactions, especially P P interactions, were strong enough to be considered in the solid state super molecular structures. Solimannejad and coworkers[17,18] reported that the structure formed with P N interaction is a global minimum and more stable than the structure with a PH N hydrogen bond. A search in the Cambridge structural database (CSD) database indicates that pnicogen bond interaction exists widely in solid phase and shows a geometrical preference for the shorter contacts.[19] In addition, the pnicogen–p interaction is found in many biological systems where it likely participates in the inhibition of Sb-based drugs for the leishmaniasis treatment.[20] Other studies also showed its potential application in the self-assembly controlling.[21–23] In a pnicogen bond,[24] the pnicogen atom (P, As, or Sb) acts as Lewis acid and interacts with electron donor molecule,[25] which is similar to the halogen atom in halogen bonds. Unlike the bond angle of halogen bonds that is always 180 , the bond angle of the pnicogen bond usually ranges from 160 to
DOI: 10.1002/jcc.23819
[a] F. Liu, J. Gao, L. Wang, C. Liu Key Lab of Colloid and Interface Chemistry, Ministry of Education, Institute of Theoretical Chemistry, School of Chemistry & Chemical Engineering, Shandong University, Jinan, 250100, People’s Republic of China E-mail: [email protected] (J. Gao) [b] L. Du Laboratory of Bio-based Materials, Qingdao Institute of Bio-energy and Bioprocess Technology, Chinese Academy of Sciences, Qingdao, 266101 Shandong, People’s Republic of China E-mail: [email protected] (L. Du) [c] B. Song Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Laboratory of Physical Biology, Shanghai, 201800 People’s Republic of China Contract grant sponsor: National Natural Science Foundation of China; Contract grant number: 91127014 and 21373124; Contract grant sponsor: Doctoral Fund of Ministry of Education of China; Contract grant number: 20110131120010; Contract grant sponsor: National Supercomputer Center in Jinan and Shandong University High Performance Computing Center C 2015 Wiley Periodicals, Inc. V
Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
1
FULL PAPER
WWW.C-CHEM.ORG
Figure 1. Electrostatic potential isosurfaces of PH2F, PHF2, and PF3 molecules (blue area: positive potential; red area: negative potential). The electrostatic potential is mapped on the surface of molecular electron density at 0.001 e/au3. The calculations were performed at MP2/aug-cc-pVDZ level.
OPLS,[31] AMBER,[32] and CHARMM[33] ignore the anisotropies of the electronic densities of the covalently bonded atoms. Recently, the positive crown (r-hole) of halogens was approximated with a massless and charged pseudoatom.[34,35] This approach is usually implemented as positive extra-point or explicit r-hole.[35,36] A more sophisticated electrostatic model, the polarizable multipoles model was proposed recently to accurately describe the electrostatics in these noncovalent bonds.[37] We previously reported a polarizable ellipsoidal force field (PEff ) model to determine the halogen bonding interaction where the anisotropic charge distribution was represented with a negatively charged sphere and a positively charged ellipsoid.[38] The results indicate that the model can correctly reproduce the potential energy surface of halogen bonds at MP2 level. In this work, we mainly focus on the theoretical analysis of the nature of pnicogen bonds and the possibility of applying the PEff model to pnicogen bonds, especially to those involving phosphorus (P) atoms, and the mono-, di-, and trisubstitution effects are carefully examined. Although the pnicogen bonds formed by other pnicogen atoms (As or Sb) are also possible, the P atom is usually taken as a typical element in discussion of pnicogen bonds. In addition, the formula and parameters fitting protocol in PEff model can be directly applied to other pnicogen atoms. The shapes of the electrostatic potential around pnicogen bonds are normally different from those of the electrostatic potential around halogen bonds. This may be attributed to the three covalently bonded atoms with the P atom (Fig. 1) and the one atom noncovalently bonded to the halogen atoms. Thus, the r-hole around P atom cannot extend along any A–P bond axis. Therefore, a modification was made to account for the difference between halogen bond and pnicogen bond. In addition, modified model can be applied to other types of noncovalent interactions other than pnicogen bond.
Theory and Methods Ab initio calculations To validate the PEff model for pnicogen bond, We selected 24 typical complexes from mono-, di-, and tri-substituted PH3– NH3 systems. The geometries of all the complexes used in this work were optimized at the MP2/aug-cc-pVDZ level of theory with a Gaussian 09 program,[39] and the optimized geometries and their Cartesian coordinates were given in Supporting Information. This level of theory is consistent with available 2
Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
experimental quantities of pnicogen bond[14,40] and the gold standard CCSD(T) method with larger basis sets.[41] The potential energy surfaces (PESs) were also generated at same level of theory with basis set superposition error (BSSE) correction[42] by varying the pnicogen bond distance (r) and angle (h). The PESs obtained with MP2 method was used as a benchmark to assess the accuracy of various density functional theory (DFT)-based PESs. The symmetry adapted perturbation theory (SAPT) method was used to determine the relative interaction energies of hydrogen bond,[43,44] cation–p interaction[45] and halogen bond.[10] The DFT-SAPT calculations at PBE0/aug-cc-VTZ level were performed with the MOLPRO2010.1 program.[46] This method promised to yield a reliable energy component that agrees remarkably well with SAPT method and with a much cheaper computational cost.[47–49] The localized molecular orbital energy decomposition analysis (LMO–EDA)[50] implemented in GAMESS(US)[51] was also used to characterize the pnicogen-bonding interactions. The modified PEff model The original PEff model was used to describe the typical noncovalent halogen-bonding interactions in our previous work.[38] In this work, a modification of the model was proposed to explore its suitability to the pnicogen bond. The detailed description of the original PEff model can be found in our previous report. The equations of the bonded interactions are kept to be the same as the general AMBER force field.[52] The nonbonded potential (V) consists of three fundamental interactions, for example, the electrostatic interaction (ESI; Velst), the two-body repulsive and dispersion interaction (Vrd), and the polarization interaction (Vpol) [eq. (1)]. V5V elst 1V rep2disp 1V pol
(1)
In a pnicogen bond, three atoms are bonded to a pnicogen atom directly, which is different from halogen bonds. This structural feature of pnicogen can strongly affect its electron density distribution and the electrostatic potential around the pnicogen atom. Compared with those of the covalently bonded halogens, the electrostatic potential distributions of covalently bonded pnicogen atoms are more complicated (Fig. 1). Therefore, the r-hole position of pnicogen may not be easily defined. In this work, the direction along the pnicogen atom and the r-hole was designated as Z-axis. The electrostatic potential distribution then becomes relatively symmetrical along Z-axis (Fig. 2). An extensive negative charge accumulation occurs in the equatorial area that is perpendicular to the Z-axis. The Z-axis is defined with K parameters and the bond vectors (V1, V2, and V3) of the three atoms covalently bonded to the pnicogen atom [eq. (2)]. Z5k1 V1 1k2 V2 1k3 V3
(2)
It can be seen from eq. (2) that the position of the r-hole is very flexible. Here, the vector (k1, k2, k3) is donated as K, which is decided by the electronegativity of the related atoms. The WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
Further detailed discussions were also collected in the Supporting Information.
Results and Discussion The performance of the improved PEff model was validated with a PH2F–NH3 complex that has been well studied as a model system for pnicogen bonds.[13,40,53,54] The critical characters of its PES at ab initio level and the accuracy of the PES at DFT level were discussed. The quality of the PES produced by the PEff model was examined. Various mono- and polysubstituted pnicogen-bonded complexes were used to illustrate the capability of the PEff model. The root mean squared derivation of the total interaction energy PESs obtained with PEff model from these obtained with MP2 method is less than 0.5 kcal/mol. Ab initio and DFT calculations Figure 2. The electrostatic potential isosurface of PH2F. The direction along the atom P and r-hole is designated as Z-axis. The electrostatic potential distribution is relatively symmetrical along the Z-axis.
PEff model for halogen bonds can be considered as a special status where K 5 (0.0, 0.0, 1.0). In addition, this modified PEff model can be applied to other parallel intermolecular interactions, such as chalcogen bond. Based on the definition of Zaxis, the electrostatic potential (Velst) of the covalent pnicogen atom can be determined with the same equation for halogen atoms [eq. (3)]. V elst ðr1 ; R; rÞ5Q ½expð2ar1 2bRÞ2expð21rÞ=r
(3)
The parameters, a, b, f, and Q in eq. (3) are interpreted from the ab initio electrostatic potential to mimic the real charge distributions. The variables r1 and R are the longitudinal distance (along the Z-axis) and the transverse distance (perpendicular to the Z-axis). The variable r is the pnicogen bond length. These parameters can be optimized with the same nonlinear least-squares fit procedure as those for halogen bonds,[38] and the anisotropic charge distribution is simplified as a negatively charged sphere and a positively charged ellipsoid (Fig. S1, Supporting Information). The two-body repulsive and dispersion interaction (Vrd), and the polarization interaction (Vpol) are the same as those in the original PEff model.[38] The steepness parameter k was set to 1.3, which was very close to that for halogen bond. Due to the complex environment of covalently bonded pnicogen atom, the scaling of vdW radius was considered and included in the calculation. Furthermore, the atomic polarizability is an essential parameter in the polarization part of the force field, and the derivation of the polarizability is not a trivial task. Currently, the parameter fitting protocol has been automated by a series of homemade Python scripts, which is available on request. The transferable issues can also be partially resolved by an automated parameter fitting procedure. In our experience, the obtained parameter sets could be transferable with similar substituted groups as the example of halogen bonds.
The PH2F–NH3 complexes were optimized at the MP2/aug-ccpVDZ level. Its PESs were generated at different theoretical levels by altering the pnicogen bond distance (r) and angle (h). Figures 3a and 3c show the distance-dependent energy curves of PH2F–NH3 dimer at MP2 and CCSD(T) levels. It can be seen that the predicted bond length and bond angle at MP2 level are consistent with those at CCSD(T) level. The distancedependent error is less than 0.5 kcal/mol. Therefore, parameters obtained at MP2/aug-cc-pVDZ were used as a benchmark for the calculations at the following DFT and force field levels. The similar results of polysubstituted pnicogen-bonding complexes are given in Supporting Information, Figure S2. To evaluate the reliability of several different DFT (B3LYP,[55,56] M06,[57] B97-D,[58] BP86,[59] and PBE[60] methods, we have computed not only the binding energies and equilibrium geometrical parameters, but also the PES of the pnicogen bonds (Figs. 3 and S3, Supporting Information). The PESs error between MP2 and these DFT methods is given in Figures 3b and 3d. It can be seen that the M06 and PBE functional overestimate the binding energies by 1.5 kcal/ mol at equilibrium distance (Fig. 3b), while the B3LYP and BP86 functional underestimate the binding energy around the equilibrium distance. For the angular dependent PES (Fig. 3d), the B97D and BP86 functional seem to provide good estimation of the binding energy around the equilibrium bond angle. For both radical and angular PES, the B97D give similar results as the MP2, which may attribute a combined effect from the reoptimized parameter of B97 functional, as well as the associated dispersion correction.[58,61] In full, different calculation models may give different predicted binding energies. Therefore, the functional of DFT method should be carefully chosen and validated to characterize the pnicogen bond interactions. Substitution effect of the pnicogen bonds The PH3 molecule has a very flat electrostatic potential surface. The binding energy between the two hydrides (PH3–NH3) is less than 2.0 kcal/mol and no positive charge region (r-hole) is Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
3
FULL PAPER
WWW.C-CHEM.ORG
Figure 3. a) Distance dependent PESs with equilibrium bond angle (h 5 167.5 ) at CCSD(T) and MP2 levels. b) The bias of distance-dependent PESs ˚ ) at CCSD(T) and MP2 levels. d) The bias of between MP2 and DFT methods. c) Angular-dependent energy curves at equilibrium bond distance (r 5 2.62 A angular-dependent PESs between MP2 and DFT methods. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
around either molecule. However, the binding energy can be magnified as much as five times when the H atom of PH3 is substituted by an electron-withdrawing group, such as F, Cl, and NO2. Electrostatic effects often dominate the weak intermolecular interactions.[62,63] Therefore, the substitution could contribute to intermolecular interaction (P N) by altering their Coulombic forces. As shown in Figure 1a, the substitution of the H atom with an F atom leads to the positive charge surface (r-hole) rotated away from the electronegative group. As a result, the broad binding energy of pnicogen bonds can range from 1 to 10 kcal/mol. Figures 1b and 1c show the electrostatic potential isosurfaces of polysubstituted PHF2 and PF3 complexes. The substitution causes multiple r-hole regions that cannot be easily distinguished. However, it is unclear whether di- or trielectronegative substitutions can increase the stability of the complex. Scheiner[18] argued that the binding energy is insensitive to the number of the substituent when PH3 molecule is substituted by F, Cl, and Br atoms paired with NH3 molecules. The additional halogen atom can slightly reduce the XAP N energy and elongate the intermolecular distance. Similar phenomenon was observed in our calculations (Figs. 4, S4, and S5, Supporting Information). The distance between P atom and N atom was elongated with the increasing number of the same substituent atom covalently bonded to P atom.
the atom directly contacts H atom of PH2F) could affect the bond angles and bond energies when the trivalent P atom interacts with an electron donor, such as NH3 molecule. The rhole of halogen bonds is usually at the 180 along AAX bond direction. However, the r-hole of pnicogen bonds is not at
The modified PEff model for pnicogen bonds
Figure 4. Optimized geometries of the monosubstituted complexes (PH2A) where A 5 a) F, b) Cl, c) Br, d) H, e) OH, f ) CN, g) CH3, and h) NO2. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
Compared with halogens, phosphorus atom is usually trivalent that can bond with three atoms. Thus, any substituent (i.e., 4
Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
˚ of PH2F–Ne and PH2F–NH3 dimers at MP2/aug-cc-pVDZ level. Figure 5. a) The angular dependencies of potential energy curves at r 5 2.6, 2.8, and 3.2 A b) The angular dependencies of the energy decomposition analysis of PH2F–NH3 with DFT–SAPT methods at PBE0/aug-cc-VTZ level, they are electrostatic interaction energy (Es), exchange energy (Eexch), induction energy (Eind), dispersion energy (Edisp), and the total interaction energy (Eint). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
any AAP bond direction. For example, the r-hole of PH2F is smaller than 170 along FAP direction. Thus, a Z-axis is introduced to locate the direction of the r-hole and to determine the contribution of the three atoms bonded to the P atom to the pnicogen bonds (Fig. 2). The Z-axis is the only modification of the original PEff model for halogen bond [eq. (2)]. In this work, the electrostatic parameters of PEff model were determined with a nonlinear least-squares optimization. Compared with the point charge model the modified PEff model can reduce the standard error of estimated ESP tens of times (Fig. S6, Supporting Information). For example, the fitted exponential parameters a, b, f, and Q of phosphorus hydrogen fluoride (PH2F) are 0.038 a.u., 0.010 a.u., 8.50 a.u., and 0.09 e. The relationship of a 6¼ b reveals that the anisotropic charge density of the pnicogen atom is higher and the electrostatic potential along the Z-axis is more positive than those in the equatorial area, consistent with theoretical electrostatic potential distribution (Fig. 2). Therefore, this model could precisely reflect the existence of the r-hole. The exponential parameters a, b, and f can be used as an ESI index that is related to the substituent effect and can be used to qualitatively determine the electrostatic properties of the covalent pnicogen atom. The effective shape of the covalently bonded phosphorus (P) atom is usually defined by the attractive dispersive London forces acting at long distance and the steric repulsive forces of the van der Waals interaction at short distance. The effective shape of the P atom in the PH2F molecule with a simple PH2F Ne interaction pair was mapped by setting the Ne ˚ away from P atom with the angles atom at 2.6, 2.8, and 3.2 A increased from 80 to 260 by 5 steps (Fig. 5a). Due to the nonpolarization of Ne atom, high level ab initio calculations in this model complex allow us to focus on the distance and angle distortion of the effect shape.[64] As shown in Figure 5a, the total interaction energy of P Ne is strongly dependent on the angle, indicating that the shape of covalent phosphorus is aspherical. The shape of the angular-dependent potential curves is asymmetric at 80 and 260 , which is different from that of halogen bond. This may be caused by the more complicated covalent conditions of the P atom. As shown in
Figure 5b, the repulsive force from DFT–SAPT calculation of PH2F–NH3 is much stronger at 80 , which can attribute to the short distance between the N atom of NH3 and the H atoms of PH2F. That is to say the electrostatic energy (Es) at 80 is lower due to the positive charge distribution around the H atom. The contribution from the exchange energy (Eexch) varies much at short distances, which can be understood by the effects of the short range contact of P atom and NH3 complex, such as overlap of electron density and charge transfer effects in pnicogen-bonded dimers. In addition, the precise determination of the polarization energy can also improve the accuracy of the force field calculation, and Figure 6 is the energy decomposition curve of PH2F–NH3 dimer derived from LMO–EDA calculations. It is clear that the polarization energy is highly dependent on the bond length and slightly dependent on the bond angle. Therefore, the polarization interaction may also contribute to the stability of pnicogen bonds. Performance of the modified PEff model in characterization of pnicogen bonds We selected these monosubstituted and polysubstituted PH3– NH3 dimers (Figs. 4, S4, and S5, Supporting Information) to examine the capability of the modified PEff model. The pnicogen bond distance (r) was varied by steps of 0.2 A˚ and the bond angle (h) was increased by steps of 5 to calculate PESs at both MP2/aug-cc-pVDZ and force field levels. The BSSE correction was included in the MP2 calculations. The fitted electrostatic parameters (a, b, f, and Q), K parameters (k1, k2, and k3), polar flatting value (Rdiff ) and isotropic polarizability (a) of these complexes are summarized in Supporting Information Tables S1 and S2. The units of these parameters are given in atomic unit unless stated otherwise. Figure 7 is the potential energy curves and dissociation energies of the monosubstituted PH2F pnicogen bond at MP2 and PEff levels. It is clear that the critical characters of the pnicogen bond obtained with PEff model (Figs. 7b and 7d), including unsymmetrical character of angular-dependent curves, are consistent to those obtained with MP2 level (Figs. Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
5
FULL PAPER
WWW.C-CHEM.ORG
Figure 6. The energy curves of the binding energy (E1), polarization energy (Epol), and binding energy without polarization (E0) of PH2F–NH3 dimer derived from LMO–EDA calculations at MP2/aug-cc-pVDZ level. a) Distance-dependent energy curves with h 5 167.5 ; b) Angular-dependent energy curves at r 5 2.62 A˚. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
7a and 7c). Compared with the point charge model (Fig. S6b, Supporting Information), the PEff model can provide correct binding energy as well as angular-dependent PESs. The difference of the potential energy curves at shorter distance is mainly caused by the overlap of electron density and HOMO/ LUMO charge transfer effects averaged in the PEff model without explicit expression. These effects have been carefully discussed in the recent work of Scheiner.[13] However, it does not affect the performance of the PEff model since extremely high PES energy region is very rare in classical molecule dynamics. In full, the modified PEff model can be used to determine the PES of the pnicogen bond interaction at MP2 level with a relative simple expression. For the monosubstituted pnicogen bonded complexes (Fig. 4), the K parameters of the two hydrogen atoms are set to 20.2, which contribute equally to the direction of the r-hole. The stable P N distances (r0) of PH3–NH3 and PH2CH3–NH3
˚ . The PESs of pnicogencomplexes are 3.3 A˚ and 2.8 A ˚ bonding interaction (r0 6 1.5 A) were then statistically analyzed. The root mean squared error (RMSE), mean unsigned/ absolute error, and mean signed error for the PES of different dimers are shown in Figure 8a. The predicted errors of the PESs of pnicogen bond obtained with the modified PEff model are less than 0.5 kcal/mol. The RMSE is 0.3 kcal/mol in the ˚ ), typical range of pnicogen-bonding interactions (r0 6 1.0 A which is acceptable since the biases of predicted binding energies of most DFT functional are higher than 0.5 kcal/mol. The suitability of the modified PEff model to di- and trisubstituted complexes was further examined. Their optimized structures are shown in Supporting Information, Figures S4 and S5. Generally speaking, it is also possible to describe the polysubstituted pnicogen-bonding complex with different substitutions with PEff model. For the case of polysubstituted systems, only one ellipsoid seems to have some limitation to
Figure 7. The angular-dependent PESs of the PH2F–NH3 dimer at a) MP2/aug-cc-pVDZ and b) PEff level; the distance-dependent PES of the PH2F–NH3 dimer at c) MP2/aug-cc-pVDZ and d) PEff level. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
6
Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
of the coexistence of various noncovalent interactions. A careful benchmark study is an important issue to provide a physical meaningful analysis of the obtained results.
Conclusions Pnicogen bonds are highly anisotropic. Their distortions away from the preferred geometries cause a more rapid loss of stability than those of H-bonds,[54] and we have compared the angle distortion between halogen bonds and H-bonds in our recent work.[65] Further comparison among hydrogen bonds, halogen bonds and pnicogen bonds can be found in the Supporting Information (Fig. S7). The PEff model that was proposed for halogen bonds is modified with a defined Z-axis to account for the substitution effects of pnicogen bonds. The original PEff model for halogen bond could be regarded as a special case of the modified model where the anisotropic charge distribution of the covalent halogen atom is simplified as a negatively charged sphere and a positively charged ellipsoid. Mono-, di-, and trisubstituted complexes were used to assess the modified PEff model. The results indicate that the model can accurately characterize ab initio PES of pnicogen bonds with the relative error of PESs less than 0.5 kcal/mol. In addition, the modified PEff model may be applied to other noncovalent interactions, such as chalcogen bond. Further work is required to understand how the noncovalent interactions cooperatively affect the characters of the compound complexes. Keywords: pnicogen bond noncovalent interaction polarizable force field electrostatic potential anisotropic Figure 8. The error analysis (in kcal/mol) of the calculated PESs between MP2/aug-cc-pVDZ method and those determined with PEff model. The column of a), b), and c) are for mono-, di-, and trisubstituted, respectively. * RMSE, root mean squared error; MUE, mean unsigned/absolute error; MSE, mean signed error.
mimic the complicated electrostatic distribution around P atom. In addition, the trends of electrostatic potential parameters and the K parameters of different complexes are harder to be obtained with the original PEff model. However, the ESI and the binding energy can be successfully determined with the properly defined Z-axis (Figs. 8b and 8c). All the error analyses are summarized in the Supporting Information (Tables S3 and S4). These results indicate that the PEff model can precisely describe the interaction between polysubstituted P atom and NH3 molecules. ˚ ) is less than 0.5 kcal/mol in The RMSE of PES (r0 6 1.0 A most pnicogen complexes (Fig. 8). Higher RMSEs of PESs (r0 6 ˚ ) were found in PHF(CN), PH(NO2)Br, and PF(CN)2 com1.0 A plexes, which might attribute to the larger substituent in these complexes. In full, the modified PEff model can be applied to various pnicogen-bonded complexes. The averages of RMSEs for various mono-, di-, and trisubstituted complexes are 0.22, 0.51, and 0.41 kcal/mol, respectively. Other weak interactions may exist and affect the result. The cooperative and competition effects are very important in the explanation of the nature
How to cite this article: F. Liu, L. Du, J. Gao, L. Wang, B. Song, C. Liu. J. Comput. Chem. 2015, DOI: 10.1002/jcc.23819
]
Additional Supporting Information may be found in the online version of this article.
[1] H. L. Nguyen, P. N. Horton, M. B. Hursthouse, A. C. Legon, D. W. Bruce, J. Am. Chem. Soc. 2003, 126, 16. [2] A. C. Carlsson, J. Gr€afenstein, A. Budnjo, J. L. Laurila, J. Bergquist, A. Karim, R. Kleinmaier, U. Brath, M. Erd elyi, J. Am. Chem. Soc. 2012, 134, 5706. [3] M. G. Chudzinski, M. S. Taylor, J. Org. Chem. 2012, 77, 3483. [4] J. S. Murray, P. Politzer, J. Mol. Struct. 1998, 425, 107. [5] A. D. Buckingham, J. E. Del Bene, S. A. C. McDowell, Chem. Phys. Lett. 2008, 463, 1. [6] C. M. Huggins, G. C. Pimentel, J. N. Shoolery, J. Chem. Phys. 1955, 23, 1244. [7] P. A. Kollman, L. C. Allen, Chem. Rev. 1972, 72, 283. [8] R. Taylor, O. Kennard, W. Versichel, J. Am. Chem. Soc. 1983, 105, 5761. [9] A. C. Legon, Phys. Chem. Chem. Phys. 2010, 12, 7736. [10] K. E. Riley, P. Hobza, J. Chem. Theory Comput. 2008, 4, 232. [11] I. Alkorta, F. Blanco, M. Solimannejad, J. Elguero, J. Phys. Chem. A. 2008, 112, 10856. [12] E. Parisini, P. Metrangolo, T. Pilati, G. Resnati, G. Terraneo, Chem. Soc. Rev. 2011, 40, 2267. [13] S. Scheiner, Acc. Chem. Res. 2012, 46, 280. [14] S. Scheiner, Int. J. Quantum Chem. 2013, 113, 1609. [15] S. Zahn, R. Frank, E. Hey-Hawkins, B. Kirchner, Chem. Eur. J. 2011, 17, 6034.
Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
7
FULL PAPER
WWW.C-CHEM.ORG
[16] J. Moilanen, C. Ganesamoorthy, M. S. Balakrishna, H. M. Tuononen, Inorg. Chem. 2009, 48, 6740. [17] M. Solimannejad, M. Gharabaghi, S. Scheiner, J. Chem. Phys. 2011, 134, 024312. [18] S. Scheiner, Chem. Phys. 2011, 387, 79. [19] G. Sanchez-Sanz, C. Trujillo, M. Solimannejad, I. Alkorta, J. Elguero, Phys. Chem. Chem. Phys. 2013, 15, 14310. [20] P. Baiocco, G. Colotti, S. Franceschini, A. Ilari, J. Med. Chem. 2009, 52, 2603. [21] W. J. Vickaryous, R. Herges, D. W. Johnson, Angew. Chem. Int. Ed. 2004, 43, 5831. [22] W. J. Vickaryous, E. R. Healey, O. B. Berryman, D. W. Johnson, Inorg. Chem. 2005, 44, 9247. [23] V. M. Cangelosi, L. N. Zakharov, D. W. Johnson, Angew. Chem. Int. Ed. 2010, 49, 1248. [24] S. Scheiner, U. Adhikari, J. Phys. Chem. A 2011, 115, 11101. [25] J. E. Del Bene, I. Alkorta, G. Sanchez-Sanz,, J. Elguero, J. Phys. Chem. A 2013, 117, 3133. [26] P. Politzer, J. S. Murray, P. Lane, Int. J. Quantum Chem. 2007, 107, 3046. [27] T. Clark, M. Hennemann, J. Murray, P. Politzer, J. Mol. Model. 2007, 13, 291. [28] U. Adhikari, S. Scheiner, J. Phys. Chem. A 2012, 116, 3487. [29] P. Metrangolo, H. Neukirch, T. Pilati, G. Resnati, Acc. Chem. Res. 2005, 38, 386. [30] N. L. Allinger, Y. H. Yuh, J. H. Lii, J. Am. Chem. Soc. 1989, 111, 8551. [31] W. L. Jorgensen, J. Tirado-Rives, J. Am. Chem. Soc. 1988, 110, 1657. [32] D. A. Pearlman, D. A. Case, J. W. Caldwell, W. S. Ross, T. E. Cheatham Iii, S. DeBolt, D. Ferguson, G. Seibel, P. Kollman, Comput. Phys. Commun. 1995, 91, 1. [33] A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wi orkiewicz-Kuczera, D. Yin, M. Karplus, J. Phys. Chem. B 1998, 102, 3586. [34] M. A. Ibrahim, J. Mol. Model. 2012, 18, 4625. [35] M. Kolar, P. Hobza, A. K. Bronowska, Chem. Commun. (Camb) 2013, 49, 981. [36] M. A. A. Ibrahim, J. Comput. Chem. 2011, 32, 2564. [37] X. Mu, Q. Wang, L.-P. Wang, S. D. Fried, J.-P. Piquemal, K. N. Dalby, P. Ren, J. Phys. Chem. B 2014, 118, 6456. [38] L. Du, J. Gao, F. Bi, L. Wang, C. Liu, J. Comput. Chem. 2013, 34, 2032. [39] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N.
8
Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
[40] [41] [42] [43] [44] [45] [46]
[47] [48] [49] [50] [51]
[52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65]
Rega, N. J. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. € Farkas, J. B. Foresman, Dannenberg, S. Dapprich, A. D. Daniels, O. J. V. Ortiz, J. Cioslowski, D. J. Fox, Gaussian: Wallingford, CT, 2009. S. Scheiner, J. Phys. Chem. A 2011, 115, 11202. J. A. Pople, M. Head-Gordon, K. Raghavachari, J. Chem. Phys. 1987, 87, 5968. S. F. Boys, F. Bernardi, Mol. Phys. 1970, 19, 553. A. Fiethen, G. Jansen,, A. Hesselmann,, M. Sch€ utz, J. Am. Chem. Soc. 2008, 130, 1802. H. Szatyłowicz, T. M. Krygowski, J. J. Panek, A. Jezierska, J. Phys. Chem. A 2008, 112, 9895. I. Soteras, M. Orozco, F. J. Luque, Phys. Chem. Chem. Phys. 2008, 10, 2616. MOLPRO, version 2010.1, a package of ab initio programs, H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Sch€ utz, and others, see http://www.molpro.net. A. Heßelmann, G. Jansen, Phys. Chem. Chem. Phys. 2003, 5, 5010. A. Heßelmann, G. Jansen, Chem. Phys. Lett. 2003, 367, 778. H. L. Williams, C. F. Chabalowski, J. Phys. Chem. A 2000, 105, 646. P. Su, H. Li, J. Chem. Phys. 2009, 131, 014102. M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, J. A. Montgomery, J. Comput. Chem. 1993, 14, 1347. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, D. A. Case, J. Comput. Chem. 2004, 25, 1157. U. Adhikari, S. Scheiner, Chem. Phys. Lett. 2012, 536, 30. U. Adhikari, S. Scheiner, Chem. Phys. Lett. 2012, 532, 31. C. Lee, W. Yang, R. G. Parr, Phys. Rev. B 1988, 37, 785. A. D. Becke, Phys. Rev. A 1988, 38, 3098. Y. Zhao, D. Truhlar, Theor. Chem. Acc. 2008, 120, 215. S. Grimme, J. Comput. Chem. 2006, 27, 1787. J. P. Perdew, Phys. Rev. B 1986, 33, 8822. J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 1996, 77, 3865. A. Bauza, I. Alkorta, A. Frontera, J. Elguero, J. Chem. Theory. Comput. 2013, 9, 5201. J. E. Del Bene, I. Alkorta, G. Sanchez-Sanz, J. Elguero, J. Phys. Chem. A 2011, 115, 13724. J. E. Del Bene, I. Alkorta, G. Sanchez-Sanz, J. Elguero, J. Phys. Chem. A 2012, 116, 3056. S. A. Peebles, P. W. Fowler, A. C. Legon, Chem. Phys. Lett. 1995, 240, 130. L. Wang, J. Gao, F. Bi, B. Song, C. Liu, J. Phys. Chem. A 2014, 118, 9140.
Received: 15 October 2014 Revised: 27 November 2014 Accepted: 6 December 2014 Published online on 00 Month 2014
WWW.CHEMISTRYVIEWS.COM
FULL PAPER
Application of Polarizable Ellipsoidal Force Field Model to Pnicogen Bonds Fang Liu,[a] Likai Du,*[b] Jun Gao,*[a] Lili Wang,[a] Bo Song,[c] and Chengbu Liu[a] Noncovalent interactions, such as hydrogen bonds and halogen bonds, are frequently used in drug designing and crystal engineering. Recently, a novel noncovalent pnicogen bonds have been identified as an important driving force in crystal structures with similar bonding mechanisms as hydrogen bond and halogen bond. Although the pnicogen bond is highly anisotropic, the pnicogen bond angles range from 160 to 180 due to the complicated substituent effects. To understand the anisotropic characters of pnicogen bond, a modification of the polarizable ellipsoidal force field (PEff ) model previously used to define halogen bonds was proposed in this work. The potential energy surfaces (PESs) of mono- and poly-
substituted PH3–NH3 complexes were calculated at CCSD(T), MP2, and density functional theory levels and were used to examine the modified PEff model. The results indicate that the modified PEff model can precisely characterize pnicogen bond. The root mean squared error of PES obtained with PEff model is less than 0.5 kcal/mol, compared with MP2 results. In addition, the modified PEff model may be applied to other noncovalent bond interactions, which is important to understand the role of intermolecular interactions in the self-assembly C 2015 Wiley Periodicals, Inc. structures. V
Introduction
180 . In addition, the r-hole, a concept proposed by Politzer and Murray, can also be found in pnicogen bond.[26,27] Some properties of pnicogen bond in a few model complexes, such as bond length and binding energy have been theoretically determined through quantum mechanics calculations.[24,28] The binding energies of pnicogen bonds and halogen bonds are 1–10[13] and 1–45 kcal/mol,[29] respectively. Theoretically, the pnicogen bond is stabilized partially by electrostatics and polarization. In addition, charge-transfer also contributes to the stability of pnicogen bonds.[13] The molecular mechanics (MM) method is a very useful tool to study condensed systems. However, the traditional force fields that assign simple isotropic point charges to atoms cannot be applied to pnicogen bond complexes due to their unique structural and thermodynamic properties. The traditional point charge model-based force fields, such as MM3,[30]
Noncovalent interactions play important roles in the functions and properties of condensed phases, biomolecular crystals and solutions.[1–4] The hydrogen bond (H-bond) is one of the most well known and understood noncovalent interactions.[5–8] The substitution of the H atom of H-bond by a halogen atom leads to another kind of noncovalent interaction, the halogen bond (X-bond).[9–12] Recent work has shown that it is not only hydrogen and halogen atoms can interact with an electron-donor, but also adjacent main-group elements, such as pnicogen and chalcogen atoms.[13,14] The noncovalent pnicogen bond was first found in crystal structures. They are as strong as hydrogen bonds and can act as a molecular linker in different chemical systems.[15] Since then, the pnicogen bond has received considerable attentions. Tuononen and coworkers[16] showed that the pnicogen bonding interactions, especially P P interactions, were strong enough to be considered in the solid state super molecular structures. Solimannejad and coworkers[17,18] reported that the structure formed with P N interaction is a global minimum and more stable than the structure with a PH N hydrogen bond. A search in the Cambridge structural database (CSD) database indicates that pnicogen bond interaction exists widely in solid phase and shows a geometrical preference for the shorter contacts.[19] In addition, the pnicogen–p interaction is found in many biological systems where it likely participates in the inhibition of Sb-based drugs for the leishmaniasis treatment.[20] Other studies also showed its potential application in the self-assembly controlling.[21–23] In a pnicogen bond,[24] the pnicogen atom (P, As, or Sb) acts as Lewis acid and interacts with electron donor molecule,[25] which is similar to the halogen atom in halogen bonds. Unlike the bond angle of halogen bonds that is always 180 , the bond angle of the pnicogen bond usually ranges from 160 to
DOI: 10.1002/jcc.23819
[a] F. Liu, J. Gao, L. Wang, C. Liu Key Lab of Colloid and Interface Chemistry, Ministry of Education, Institute of Theoretical Chemistry, School of Chemistry & Chemical Engineering, Shandong University, Jinan, 250100, People’s Republic of China E-mail: [email protected] (J. Gao) [b] L. Du Laboratory of Bio-based Materials, Qingdao Institute of Bio-energy and Bioprocess Technology, Chinese Academy of Sciences, Qingdao, 266101 Shandong, People’s Republic of China E-mail: [email protected] (L. Du) [c] B. Song Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Laboratory of Physical Biology, Shanghai, 201800 People’s Republic of China Contract grant sponsor: National Natural Science Foundation of China; Contract grant number: 91127014 and 21373124; Contract grant sponsor: Doctoral Fund of Ministry of Education of China; Contract grant number: 20110131120010; Contract grant sponsor: National Supercomputer Center in Jinan and Shandong University High Performance Computing Center C 2015 Wiley Periodicals, Inc. V
Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
1
FULL PAPER
WWW.C-CHEM.ORG
Figure 1. Electrostatic potential isosurfaces of PH2F, PHF2, and PF3 molecules (blue area: positive potential; red area: negative potential). The electrostatic potential is mapped on the surface of molecular electron density at 0.001 e/au3. The calculations were performed at MP2/aug-cc-pVDZ level.
OPLS,[31] AMBER,[32] and CHARMM[33] ignore the anisotropies of the electronic densities of the covalently bonded atoms. Recently, the positive crown (r-hole) of halogens was approximated with a massless and charged pseudoatom.[34,35] This approach is usually implemented as positive extra-point or explicit r-hole.[35,36] A more sophisticated electrostatic model, the polarizable multipoles model was proposed recently to accurately describe the electrostatics in these noncovalent bonds.[37] We previously reported a polarizable ellipsoidal force field (PEff ) model to determine the halogen bonding interaction where the anisotropic charge distribution was represented with a negatively charged sphere and a positively charged ellipsoid.[38] The results indicate that the model can correctly reproduce the potential energy surface of halogen bonds at MP2 level. In this work, we mainly focus on the theoretical analysis of the nature of pnicogen bonds and the possibility of applying the PEff model to pnicogen bonds, especially to those involving phosphorus (P) atoms, and the mono-, di-, and trisubstitution effects are carefully examined. Although the pnicogen bonds formed by other pnicogen atoms (As or Sb) are also possible, the P atom is usually taken as a typical element in discussion of pnicogen bonds. In addition, the formula and parameters fitting protocol in PEff model can be directly applied to other pnicogen atoms. The shapes of the electrostatic potential around pnicogen bonds are normally different from those of the electrostatic potential around halogen bonds. This may be attributed to the three covalently bonded atoms with the P atom (Fig. 1) and the one atom noncovalently bonded to the halogen atoms. Thus, the r-hole around P atom cannot extend along any A–P bond axis. Therefore, a modification was made to account for the difference between halogen bond and pnicogen bond. In addition, modified model can be applied to other types of noncovalent interactions other than pnicogen bond.
Theory and Methods Ab initio calculations To validate the PEff model for pnicogen bond, We selected 24 typical complexes from mono-, di-, and tri-substituted PH3– NH3 systems. The geometries of all the complexes used in this work were optimized at the MP2/aug-cc-pVDZ level of theory with a Gaussian 09 program,[39] and the optimized geometries and their Cartesian coordinates were given in Supporting Information. This level of theory is consistent with available 2
Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
experimental quantities of pnicogen bond[14,40] and the gold standard CCSD(T) method with larger basis sets.[41] The potential energy surfaces (PESs) were also generated at same level of theory with basis set superposition error (BSSE) correction[42] by varying the pnicogen bond distance (r) and angle (h). The PESs obtained with MP2 method was used as a benchmark to assess the accuracy of various density functional theory (DFT)-based PESs. The symmetry adapted perturbation theory (SAPT) method was used to determine the relative interaction energies of hydrogen bond,[43,44] cation–p interaction[45] and halogen bond.[10] The DFT-SAPT calculations at PBE0/aug-cc-VTZ level were performed with the MOLPRO2010.1 program.[46] This method promised to yield a reliable energy component that agrees remarkably well with SAPT method and with a much cheaper computational cost.[47–49] The localized molecular orbital energy decomposition analysis (LMO–EDA)[50] implemented in GAMESS(US)[51] was also used to characterize the pnicogen-bonding interactions. The modified PEff model The original PEff model was used to describe the typical noncovalent halogen-bonding interactions in our previous work.[38] In this work, a modification of the model was proposed to explore its suitability to the pnicogen bond. The detailed description of the original PEff model can be found in our previous report. The equations of the bonded interactions are kept to be the same as the general AMBER force field.[52] The nonbonded potential (V) consists of three fundamental interactions, for example, the electrostatic interaction (ESI; Velst), the two-body repulsive and dispersion interaction (Vrd), and the polarization interaction (Vpol) [eq. (1)]. V5V elst 1V rep2disp 1V pol
(1)
In a pnicogen bond, three atoms are bonded to a pnicogen atom directly, which is different from halogen bonds. This structural feature of pnicogen can strongly affect its electron density distribution and the electrostatic potential around the pnicogen atom. Compared with those of the covalently bonded halogens, the electrostatic potential distributions of covalently bonded pnicogen atoms are more complicated (Fig. 1). Therefore, the r-hole position of pnicogen may not be easily defined. In this work, the direction along the pnicogen atom and the r-hole was designated as Z-axis. The electrostatic potential distribution then becomes relatively symmetrical along Z-axis (Fig. 2). An extensive negative charge accumulation occurs in the equatorial area that is perpendicular to the Z-axis. The Z-axis is defined with K parameters and the bond vectors (V1, V2, and V3) of the three atoms covalently bonded to the pnicogen atom [eq. (2)]. Z5k1 V1 1k2 V2 1k3 V3
(2)
It can be seen from eq. (2) that the position of the r-hole is very flexible. Here, the vector (k1, k2, k3) is donated as K, which is decided by the electronegativity of the related atoms. The WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
Further detailed discussions were also collected in the Supporting Information.
Results and Discussion The performance of the improved PEff model was validated with a PH2F–NH3 complex that has been well studied as a model system for pnicogen bonds.[13,40,53,54] The critical characters of its PES at ab initio level and the accuracy of the PES at DFT level were discussed. The quality of the PES produced by the PEff model was examined. Various mono- and polysubstituted pnicogen-bonded complexes were used to illustrate the capability of the PEff model. The root mean squared derivation of the total interaction energy PESs obtained with PEff model from these obtained with MP2 method is less than 0.5 kcal/mol. Ab initio and DFT calculations Figure 2. The electrostatic potential isosurface of PH2F. The direction along the atom P and r-hole is designated as Z-axis. The electrostatic potential distribution is relatively symmetrical along the Z-axis.
PEff model for halogen bonds can be considered as a special status where K 5 (0.0, 0.0, 1.0). In addition, this modified PEff model can be applied to other parallel intermolecular interactions, such as chalcogen bond. Based on the definition of Zaxis, the electrostatic potential (Velst) of the covalent pnicogen atom can be determined with the same equation for halogen atoms [eq. (3)]. V elst ðr1 ; R; rÞ5Q ½expð2ar1 2bRÞ2expð21rÞ=r
(3)
The parameters, a, b, f, and Q in eq. (3) are interpreted from the ab initio electrostatic potential to mimic the real charge distributions. The variables r1 and R are the longitudinal distance (along the Z-axis) and the transverse distance (perpendicular to the Z-axis). The variable r is the pnicogen bond length. These parameters can be optimized with the same nonlinear least-squares fit procedure as those for halogen bonds,[38] and the anisotropic charge distribution is simplified as a negatively charged sphere and a positively charged ellipsoid (Fig. S1, Supporting Information). The two-body repulsive and dispersion interaction (Vrd), and the polarization interaction (Vpol) are the same as those in the original PEff model.[38] The steepness parameter k was set to 1.3, which was very close to that for halogen bond. Due to the complex environment of covalently bonded pnicogen atom, the scaling of vdW radius was considered and included in the calculation. Furthermore, the atomic polarizability is an essential parameter in the polarization part of the force field, and the derivation of the polarizability is not a trivial task. Currently, the parameter fitting protocol has been automated by a series of homemade Python scripts, which is available on request. The transferable issues can also be partially resolved by an automated parameter fitting procedure. In our experience, the obtained parameter sets could be transferable with similar substituted groups as the example of halogen bonds.
The PH2F–NH3 complexes were optimized at the MP2/aug-ccpVDZ level. Its PESs were generated at different theoretical levels by altering the pnicogen bond distance (r) and angle (h). Figures 3a and 3c show the distance-dependent energy curves of PH2F–NH3 dimer at MP2 and CCSD(T) levels. It can be seen that the predicted bond length and bond angle at MP2 level are consistent with those at CCSD(T) level. The distancedependent error is less than 0.5 kcal/mol. Therefore, parameters obtained at MP2/aug-cc-pVDZ were used as a benchmark for the calculations at the following DFT and force field levels. The similar results of polysubstituted pnicogen-bonding complexes are given in Supporting Information, Figure S2. To evaluate the reliability of several different DFT (B3LYP,[55,56] M06,[57] B97-D,[58] BP86,[59] and PBE[60] methods, we have computed not only the binding energies and equilibrium geometrical parameters, but also the PES of the pnicogen bonds (Figs. 3 and S3, Supporting Information). The PESs error between MP2 and these DFT methods is given in Figures 3b and 3d. It can be seen that the M06 and PBE functional overestimate the binding energies by 1.5 kcal/ mol at equilibrium distance (Fig. 3b), while the B3LYP and BP86 functional underestimate the binding energy around the equilibrium distance. For the angular dependent PES (Fig. 3d), the B97D and BP86 functional seem to provide good estimation of the binding energy around the equilibrium bond angle. For both radical and angular PES, the B97D give similar results as the MP2, which may attribute a combined effect from the reoptimized parameter of B97 functional, as well as the associated dispersion correction.[58,61] In full, different calculation models may give different predicted binding energies. Therefore, the functional of DFT method should be carefully chosen and validated to characterize the pnicogen bond interactions. Substitution effect of the pnicogen bonds The PH3 molecule has a very flat electrostatic potential surface. The binding energy between the two hydrides (PH3–NH3) is less than 2.0 kcal/mol and no positive charge region (r-hole) is Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
3
FULL PAPER
WWW.C-CHEM.ORG
Figure 3. a) Distance dependent PESs with equilibrium bond angle (h 5 167.5 ) at CCSD(T) and MP2 levels. b) The bias of distance-dependent PESs ˚ ) at CCSD(T) and MP2 levels. d) The bias of between MP2 and DFT methods. c) Angular-dependent energy curves at equilibrium bond distance (r 5 2.62 A angular-dependent PESs between MP2 and DFT methods. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
around either molecule. However, the binding energy can be magnified as much as five times when the H atom of PH3 is substituted by an electron-withdrawing group, such as F, Cl, and NO2. Electrostatic effects often dominate the weak intermolecular interactions.[62,63] Therefore, the substitution could contribute to intermolecular interaction (P N) by altering their Coulombic forces. As shown in Figure 1a, the substitution of the H atom with an F atom leads to the positive charge surface (r-hole) rotated away from the electronegative group. As a result, the broad binding energy of pnicogen bonds can range from 1 to 10 kcal/mol. Figures 1b and 1c show the electrostatic potential isosurfaces of polysubstituted PHF2 and PF3 complexes. The substitution causes multiple r-hole regions that cannot be easily distinguished. However, it is unclear whether di- or trielectronegative substitutions can increase the stability of the complex. Scheiner[18] argued that the binding energy is insensitive to the number of the substituent when PH3 molecule is substituted by F, Cl, and Br atoms paired with NH3 molecules. The additional halogen atom can slightly reduce the XAP N energy and elongate the intermolecular distance. Similar phenomenon was observed in our calculations (Figs. 4, S4, and S5, Supporting Information). The distance between P atom and N atom was elongated with the increasing number of the same substituent atom covalently bonded to P atom.
the atom directly contacts H atom of PH2F) could affect the bond angles and bond energies when the trivalent P atom interacts with an electron donor, such as NH3 molecule. The rhole of halogen bonds is usually at the 180 along AAX bond direction. However, the r-hole of pnicogen bonds is not at
The modified PEff model for pnicogen bonds
Figure 4. Optimized geometries of the monosubstituted complexes (PH2A) where A 5 a) F, b) Cl, c) Br, d) H, e) OH, f ) CN, g) CH3, and h) NO2. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
Compared with halogens, phosphorus atom is usually trivalent that can bond with three atoms. Thus, any substituent (i.e., 4
Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
˚ of PH2F–Ne and PH2F–NH3 dimers at MP2/aug-cc-pVDZ level. Figure 5. a) The angular dependencies of potential energy curves at r 5 2.6, 2.8, and 3.2 A b) The angular dependencies of the energy decomposition analysis of PH2F–NH3 with DFT–SAPT methods at PBE0/aug-cc-VTZ level, they are electrostatic interaction energy (Es), exchange energy (Eexch), induction energy (Eind), dispersion energy (Edisp), and the total interaction energy (Eint). [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
any AAP bond direction. For example, the r-hole of PH2F is smaller than 170 along FAP direction. Thus, a Z-axis is introduced to locate the direction of the r-hole and to determine the contribution of the three atoms bonded to the P atom to the pnicogen bonds (Fig. 2). The Z-axis is the only modification of the original PEff model for halogen bond [eq. (2)]. In this work, the electrostatic parameters of PEff model were determined with a nonlinear least-squares optimization. Compared with the point charge model the modified PEff model can reduce the standard error of estimated ESP tens of times (Fig. S6, Supporting Information). For example, the fitted exponential parameters a, b, f, and Q of phosphorus hydrogen fluoride (PH2F) are 0.038 a.u., 0.010 a.u., 8.50 a.u., and 0.09 e. The relationship of a 6¼ b reveals that the anisotropic charge density of the pnicogen atom is higher and the electrostatic potential along the Z-axis is more positive than those in the equatorial area, consistent with theoretical electrostatic potential distribution (Fig. 2). Therefore, this model could precisely reflect the existence of the r-hole. The exponential parameters a, b, and f can be used as an ESI index that is related to the substituent effect and can be used to qualitatively determine the electrostatic properties of the covalent pnicogen atom. The effective shape of the covalently bonded phosphorus (P) atom is usually defined by the attractive dispersive London forces acting at long distance and the steric repulsive forces of the van der Waals interaction at short distance. The effective shape of the P atom in the PH2F molecule with a simple PH2F Ne interaction pair was mapped by setting the Ne ˚ away from P atom with the angles atom at 2.6, 2.8, and 3.2 A increased from 80 to 260 by 5 steps (Fig. 5a). Due to the nonpolarization of Ne atom, high level ab initio calculations in this model complex allow us to focus on the distance and angle distortion of the effect shape.[64] As shown in Figure 5a, the total interaction energy of P Ne is strongly dependent on the angle, indicating that the shape of covalent phosphorus is aspherical. The shape of the angular-dependent potential curves is asymmetric at 80 and 260 , which is different from that of halogen bond. This may be caused by the more complicated covalent conditions of the P atom. As shown in
Figure 5b, the repulsive force from DFT–SAPT calculation of PH2F–NH3 is much stronger at 80 , which can attribute to the short distance between the N atom of NH3 and the H atoms of PH2F. That is to say the electrostatic energy (Es) at 80 is lower due to the positive charge distribution around the H atom. The contribution from the exchange energy (Eexch) varies much at short distances, which can be understood by the effects of the short range contact of P atom and NH3 complex, such as overlap of electron density and charge transfer effects in pnicogen-bonded dimers. In addition, the precise determination of the polarization energy can also improve the accuracy of the force field calculation, and Figure 6 is the energy decomposition curve of PH2F–NH3 dimer derived from LMO–EDA calculations. It is clear that the polarization energy is highly dependent on the bond length and slightly dependent on the bond angle. Therefore, the polarization interaction may also contribute to the stability of pnicogen bonds. Performance of the modified PEff model in characterization of pnicogen bonds We selected these monosubstituted and polysubstituted PH3– NH3 dimers (Figs. 4, S4, and S5, Supporting Information) to examine the capability of the modified PEff model. The pnicogen bond distance (r) was varied by steps of 0.2 A˚ and the bond angle (h) was increased by steps of 5 to calculate PESs at both MP2/aug-cc-pVDZ and force field levels. The BSSE correction was included in the MP2 calculations. The fitted electrostatic parameters (a, b, f, and Q), K parameters (k1, k2, and k3), polar flatting value (Rdiff ) and isotropic polarizability (a) of these complexes are summarized in Supporting Information Tables S1 and S2. The units of these parameters are given in atomic unit unless stated otherwise. Figure 7 is the potential energy curves and dissociation energies of the monosubstituted PH2F pnicogen bond at MP2 and PEff levels. It is clear that the critical characters of the pnicogen bond obtained with PEff model (Figs. 7b and 7d), including unsymmetrical character of angular-dependent curves, are consistent to those obtained with MP2 level (Figs. Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
5
FULL PAPER
WWW.C-CHEM.ORG
Figure 6. The energy curves of the binding energy (E1), polarization energy (Epol), and binding energy without polarization (E0) of PH2F–NH3 dimer derived from LMO–EDA calculations at MP2/aug-cc-pVDZ level. a) Distance-dependent energy curves with h 5 167.5 ; b) Angular-dependent energy curves at r 5 2.62 A˚. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
7a and 7c). Compared with the point charge model (Fig. S6b, Supporting Information), the PEff model can provide correct binding energy as well as angular-dependent PESs. The difference of the potential energy curves at shorter distance is mainly caused by the overlap of electron density and HOMO/ LUMO charge transfer effects averaged in the PEff model without explicit expression. These effects have been carefully discussed in the recent work of Scheiner.[13] However, it does not affect the performance of the PEff model since extremely high PES energy region is very rare in classical molecule dynamics. In full, the modified PEff model can be used to determine the PES of the pnicogen bond interaction at MP2 level with a relative simple expression. For the monosubstituted pnicogen bonded complexes (Fig. 4), the K parameters of the two hydrogen atoms are set to 20.2, which contribute equally to the direction of the r-hole. The stable P N distances (r0) of PH3–NH3 and PH2CH3–NH3
˚ . The PESs of pnicogencomplexes are 3.3 A˚ and 2.8 A ˚ bonding interaction (r0 6 1.5 A) were then statistically analyzed. The root mean squared error (RMSE), mean unsigned/ absolute error, and mean signed error for the PES of different dimers are shown in Figure 8a. The predicted errors of the PESs of pnicogen bond obtained with the modified PEff model are less than 0.5 kcal/mol. The RMSE is 0.3 kcal/mol in the ˚ ), typical range of pnicogen-bonding interactions (r0 6 1.0 A which is acceptable since the biases of predicted binding energies of most DFT functional are higher than 0.5 kcal/mol. The suitability of the modified PEff model to di- and trisubstituted complexes was further examined. Their optimized structures are shown in Supporting Information, Figures S4 and S5. Generally speaking, it is also possible to describe the polysubstituted pnicogen-bonding complex with different substitutions with PEff model. For the case of polysubstituted systems, only one ellipsoid seems to have some limitation to
Figure 7. The angular-dependent PESs of the PH2F–NH3 dimer at a) MP2/aug-cc-pVDZ and b) PEff level; the distance-dependent PES of the PH2F–NH3 dimer at c) MP2/aug-cc-pVDZ and d) PEff level. [Color figure can be viewed in the online issue, which is available at wileyonlinelibrary.com.]
6
Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
WWW.CHEMISTRYVIEWS.COM
WWW.C-CHEM.ORG
FULL PAPER
of the coexistence of various noncovalent interactions. A careful benchmark study is an important issue to provide a physical meaningful analysis of the obtained results.
Conclusions Pnicogen bonds are highly anisotropic. Their distortions away from the preferred geometries cause a more rapid loss of stability than those of H-bonds,[54] and we have compared the angle distortion between halogen bonds and H-bonds in our recent work.[65] Further comparison among hydrogen bonds, halogen bonds and pnicogen bonds can be found in the Supporting Information (Fig. S7). The PEff model that was proposed for halogen bonds is modified with a defined Z-axis to account for the substitution effects of pnicogen bonds. The original PEff model for halogen bond could be regarded as a special case of the modified model where the anisotropic charge distribution of the covalent halogen atom is simplified as a negatively charged sphere and a positively charged ellipsoid. Mono-, di-, and trisubstituted complexes were used to assess the modified PEff model. The results indicate that the model can accurately characterize ab initio PES of pnicogen bonds with the relative error of PESs less than 0.5 kcal/mol. In addition, the modified PEff model may be applied to other noncovalent interactions, such as chalcogen bond. Further work is required to understand how the noncovalent interactions cooperatively affect the characters of the compound complexes. Keywords: pnicogen bond noncovalent interaction polarizable force field electrostatic potential anisotropic Figure 8. The error analysis (in kcal/mol) of the calculated PESs between MP2/aug-cc-pVDZ method and those determined with PEff model. The column of a), b), and c) are for mono-, di-, and trisubstituted, respectively. * RMSE, root mean squared error; MUE, mean unsigned/absolute error; MSE, mean signed error.
mimic the complicated electrostatic distribution around P atom. In addition, the trends of electrostatic potential parameters and the K parameters of different complexes are harder to be obtained with the original PEff model. However, the ESI and the binding energy can be successfully determined with the properly defined Z-axis (Figs. 8b and 8c). All the error analyses are summarized in the Supporting Information (Tables S3 and S4). These results indicate that the PEff model can precisely describe the interaction between polysubstituted P atom and NH3 molecules. ˚ ) is less than 0.5 kcal/mol in The RMSE of PES (r0 6 1.0 A most pnicogen complexes (Fig. 8). Higher RMSEs of PESs (r0 6 ˚ ) were found in PHF(CN), PH(NO2)Br, and PF(CN)2 com1.0 A plexes, which might attribute to the larger substituent in these complexes. In full, the modified PEff model can be applied to various pnicogen-bonded complexes. The averages of RMSEs for various mono-, di-, and trisubstituted complexes are 0.22, 0.51, and 0.41 kcal/mol, respectively. Other weak interactions may exist and affect the result. The cooperative and competition effects are very important in the explanation of the nature
How to cite this article: F. Liu, L. Du, J. Gao, L. Wang, B. Song, C. Liu. J. Comput. Chem. 2015, DOI: 10.1002/jcc.23819
]
Additional Supporting Information may be found in the online version of this article.
[1] H. L. Nguyen, P. N. Horton, M. B. Hursthouse, A. C. Legon, D. W. Bruce, J. Am. Chem. Soc. 2003, 126, 16. [2] A. C. Carlsson, J. Gr€afenstein, A. Budnjo, J. L. Laurila, J. Bergquist, A. Karim, R. Kleinmaier, U. Brath, M. Erd elyi, J. Am. Chem. Soc. 2012, 134, 5706. [3] M. G. Chudzinski, M. S. Taylor, J. Org. Chem. 2012, 77, 3483. [4] J. S. Murray, P. Politzer, J. Mol. Struct. 1998, 425, 107. [5] A. D. Buckingham, J. E. Del Bene, S. A. C. McDowell, Chem. Phys. Lett. 2008, 463, 1. [6] C. M. Huggins, G. C. Pimentel, J. N. Shoolery, J. Chem. Phys. 1955, 23, 1244. [7] P. A. Kollman, L. C. Allen, Chem. Rev. 1972, 72, 283. [8] R. Taylor, O. Kennard, W. Versichel, J. Am. Chem. Soc. 1983, 105, 5761. [9] A. C. Legon, Phys. Chem. Chem. Phys. 2010, 12, 7736. [10] K. E. Riley, P. Hobza, J. Chem. Theory Comput. 2008, 4, 232. [11] I. Alkorta, F. Blanco, M. Solimannejad, J. Elguero, J. Phys. Chem. A. 2008, 112, 10856. [12] E. Parisini, P. Metrangolo, T. Pilati, G. Resnati, G. Terraneo, Chem. Soc. Rev. 2011, 40, 2267. [13] S. Scheiner, Acc. Chem. Res. 2012, 46, 280. [14] S. Scheiner, Int. J. Quantum Chem. 2013, 113, 1609. [15] S. Zahn, R. Frank, E. Hey-Hawkins, B. Kirchner, Chem. Eur. J. 2011, 17, 6034.
Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
7
FULL PAPER
WWW.C-CHEM.ORG
[16] J. Moilanen, C. Ganesamoorthy, M. S. Balakrishna, H. M. Tuononen, Inorg. Chem. 2009, 48, 6740. [17] M. Solimannejad, M. Gharabaghi, S. Scheiner, J. Chem. Phys. 2011, 134, 024312. [18] S. Scheiner, Chem. Phys. 2011, 387, 79. [19] G. Sanchez-Sanz, C. Trujillo, M. Solimannejad, I. Alkorta, J. Elguero, Phys. Chem. Chem. Phys. 2013, 15, 14310. [20] P. Baiocco, G. Colotti, S. Franceschini, A. Ilari, J. Med. Chem. 2009, 52, 2603. [21] W. J. Vickaryous, R. Herges, D. W. Johnson, Angew. Chem. Int. Ed. 2004, 43, 5831. [22] W. J. Vickaryous, E. R. Healey, O. B. Berryman, D. W. Johnson, Inorg. Chem. 2005, 44, 9247. [23] V. M. Cangelosi, L. N. Zakharov, D. W. Johnson, Angew. Chem. Int. Ed. 2010, 49, 1248. [24] S. Scheiner, U. Adhikari, J. Phys. Chem. A 2011, 115, 11101. [25] J. E. Del Bene, I. Alkorta, G. Sanchez-Sanz,, J. Elguero, J. Phys. Chem. A 2013, 117, 3133. [26] P. Politzer, J. S. Murray, P. Lane, Int. J. Quantum Chem. 2007, 107, 3046. [27] T. Clark, M. Hennemann, J. Murray, P. Politzer, J. Mol. Model. 2007, 13, 291. [28] U. Adhikari, S. Scheiner, J. Phys. Chem. A 2012, 116, 3487. [29] P. Metrangolo, H. Neukirch, T. Pilati, G. Resnati, Acc. Chem. Res. 2005, 38, 386. [30] N. L. Allinger, Y. H. Yuh, J. H. Lii, J. Am. Chem. Soc. 1989, 111, 8551. [31] W. L. Jorgensen, J. Tirado-Rives, J. Am. Chem. Soc. 1988, 110, 1657. [32] D. A. Pearlman, D. A. Case, J. W. Caldwell, W. S. Ross, T. E. Cheatham Iii, S. DeBolt, D. Ferguson, G. Seibel, P. Kollman, Comput. Phys. Commun. 1995, 91, 1. [33] A. D. MacKerell, D. Bashford, M. Bellott, R. L. Dunbrack, J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wi orkiewicz-Kuczera, D. Yin, M. Karplus, J. Phys. Chem. B 1998, 102, 3586. [34] M. A. Ibrahim, J. Mol. Model. 2012, 18, 4625. [35] M. Kolar, P. Hobza, A. K. Bronowska, Chem. Commun. (Camb) 2013, 49, 981. [36] M. A. A. Ibrahim, J. Comput. Chem. 2011, 32, 2564. [37] X. Mu, Q. Wang, L.-P. Wang, S. D. Fried, J.-P. Piquemal, K. N. Dalby, P. Ren, J. Phys. Chem. B 2014, 118, 6456. [38] L. Du, J. Gao, F. Bi, L. Wang, C. Liu, J. Comput. Chem. 2013, 34, 2032. [39] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N.
8
Journal of Computational Chemistry 2015, DOI: 10.1002/jcc.23819
[40] [41] [42] [43] [44] [45] [46]
[47] [48] [49] [50] [51]
[52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64] [65]
Rega, N. J. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. € Farkas, J. B. Foresman, Dannenberg, S. Dapprich, A. D. Daniels, O. J. V. Ortiz, J. Cioslowski, D. J. Fox, Gaussian: Wallingford, CT, 2009. S. Scheiner, J. Phys. Chem. A 2011, 115, 11202. J. A. Pople, M. Head-Gordon, K. Raghavachari, J. Chem. Phys. 1987, 87, 5968. S. F. Boys, F. Bernardi, Mol. Phys. 1970, 19, 553. A. Fiethen, G. Jansen,, A. Hesselmann,, M. Sch€ utz, J. Am. Chem. Soc. 2008, 130, 1802. H. Szatyłowicz, T. M. Krygowski, J. J. Panek, A. Jezierska, J. Phys. Chem. A 2008, 112, 9895. I. Soteras, M. Orozco, F. J. Luque, Phys. Chem. Chem. Phys. 2008, 10, 2616. MOLPRO, version 2010.1, a package of ab initio programs, H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Sch€ utz, and others, see http://www.molpro.net. A. Heßelmann, G. Jansen, Phys. Chem. Chem. Phys. 2003, 5, 5010. A. Heßelmann, G. Jansen, Chem. Phys. Lett. 2003, 367, 778. H. L. Williams, C. F. Chabalowski, J. Phys. Chem. A 2000, 105, 646. P. Su, H. Li, J. Chem. Phys. 2009, 131, 014102. M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. Su, T. L. Windus, M. Dupuis, J. A. Montgomery, J. Comput. Chem. 1993, 14, 1347. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, D. A. Case, J. Comput. Chem. 2004, 25, 1157. U. Adhikari, S. Scheiner, Chem. Phys. Lett. 2012, 536, 30. U. Adhikari, S. Scheiner, Chem. Phys. Lett. 2012, 532, 31. C. Lee, W. Yang, R. G. Parr, Phys. Rev. B 1988, 37, 785. A. D. Becke, Phys. Rev. A 1988, 38, 3098. Y. Zhao, D. Truhlar, Theor. Chem. Acc. 2008, 120, 215. S. Grimme, J. Comput. Chem. 2006, 27, 1787. J. P. Perdew, Phys. Rev. B 1986, 33, 8822. J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 1996, 77, 3865. A. Bauza, I. Alkorta, A. Frontera, J. Elguero, J. Chem. Theory. Comput. 2013, 9, 5201. J. E. Del Bene, I. Alkorta, G. Sanchez-Sanz, J. Elguero, J. Phys. Chem. A 2011, 115, 13724. J. E. Del Bene, I. Alkorta, G. Sanchez-Sanz, J. Elguero, J. Phys. Chem. A 2012, 116, 3056. S. A. Peebles, P. W. Fowler, A. C. Legon, Chem. Phys. Lett. 1995, 240, 130. L. Wang, J. Gao, F. Bi, B. Song, C. Liu, J. Phys. Chem. A 2014, 118, 9140.
Received: 15 October 2014 Revised: 27 November 2014 Accepted: 6 December 2014 Published online on 00 Month 2014
WWW.CHEMISTRYVIEWS.COM
